#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>

/*
 * Local function prototypes
 */

static int _calc_wkblen (
    GeomType type,
    int      num_points,
    int      num_parts,
    int      num_subparts
    );

static int _get_wkb_polys (
    char     *cp,
    int       swap_bytes,
    Geometry *geom
    );

static int _get_wkb_header (
    char    **buf_ptr, 
    int      *swap_bytes,
    WkbType  *wkb_type
    );

static int _get_wkb_points (
    char     *cp,
    int       swap_bytes,
    Geometry *geom
    );

static int _get_wkb_lines (
    char     *cp,
    int       swap_bytes,
    Geometry *geom
    );

static int _get_wkb_polys (
    char     *cp,
    int       swap_bytes,
    Geometry *geom
    );


/***************************************************************************
*
*  geom_to_wkb  -  takes a geometry of any type and converts it 
*                  into a Well Known Binary representation.
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Purpose:
*   
*    Converts a geometry of any type into a WKB buffer. Since the structure
*    of a geometry varies by data type, the parameters of this function do
*    as well. The parts are associated with the multilinestring, and 
*    multipolygon datatypes where each part represents a linestring and 
*    a polygon, respectively. Subparts are associated with various rings
*    of a polygon.  
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Parameters:
*    geom          <Input>  == (Geometry *) The geometry data
*    max_alloced   <Output> == (int *)
*    data_len      <Output> == (int *) The number of bytes in the WKB.
*    binary        <Output> == (char **) The WKB buffer.
*
*    RETURN        <Output> == (void)
*
***************************************************************************/

void geom_to_wkb (
    Geometry  *geom,
    int       *max_alloced,
    int       *data_len, 
    char     **binary
    )
{
    int           type;
    char         *cp;
    int           i, j;
    int           tmp_num_subparts;
    int           end_subparts;
    int           tmp_num_points;
    unsigned char byte_order;

    /* Calculate data_len. */

    *data_len = _calc_wkblen (geom->type, geom->num_points,
                              geom->num_parts, geom->num_subparts);

    /* Allocate memory to the binary character array, set max_alloced
       equal to data_len if this is the first call to put_wkb.
       Otherwise, if the data_len is larger then max_alloced,
       reallocate the larger data_len to the binary character array. */

    if (*max_alloced == 0)
    {
        *binary = (char *)malloc (*data_len);
        *max_alloced = *data_len;
    }
    else if (*max_alloced < *data_len)
    {
        *binary = (char *) realloc (*binary, *data_len);
        *max_alloced = *data_len;
    }

    cp = *binary;

    ENDIAN_TEST (i);  /* Sets i = 0 if big-endian, 1 if little-endian */
    byte_order = (unsigned char) i;

    *cp = byte_order;
    cp++;

    switch (geom->type)
    {
        case geomPoint:
            type = wkbPoint;
            put_integer (&cp, 0, 1, (int *)&type);
            put_xy (&cp, 0, 1, geom->pt);
            break;

        case geomMultiPoint:
            type = wkbMultiPoint;
            put_integer (&cp, 0, 1, (int *)&type);
            put_integer (&cp, 0, 1, &geom->num_points);
            type = wkbPoint; /* Each part is a wkb point */
            for (i = 0; i < geom->num_points; i++)
            {
                *cp = byte_order;
                cp++;
                put_integer (&cp, 0, 1, (int *)&type);
                put_xy (&cp, 0, 1, &geom->pt[i]);
            }
            break;

        case geomLineString:
            type = wkbLineString;
            put_integer (&cp, 0, 1, (int *)&type);
            put_integer (&cp, 0, 1, &geom->num_points);
            put_xy (&cp, 0, geom->num_points, geom->pt);
            break;

        case geomMultiLineString:
            type = wkbMultiLineString;
            put_integer (&cp, 0, 1, (int *)&type);
            put_integer (&cp, 0, 1, &geom->num_parts);
            type = wkbLineString; /* Each part is a wkb line */
            for (i = 0; i < geom->num_parts; i++)
            {
                *cp = byte_order;
                cp++;
                put_integer (&cp, 0, 1, (int *)&type);

                if (i == geom->num_parts - 1)
                    tmp_num_points = geom->num_points - geom->offsets[i];
                else
                    tmp_num_points = geom->offsets[i+1] - geom->offsets[i];

                put_integer (&cp, 0, 1, &tmp_num_points);
                put_xy (&cp, 0, tmp_num_points, &geom->pt[geom->offsets[i]]);
            }
            break;

        case geomPolygon:
            type = wkbPolygon;
            put_integer (&cp, 0, 1, (int *)&type);
            put_integer (&cp, 0, 1, &geom->num_subparts);
            for (i = 0; i < geom->num_subparts; i++)
            {
                if (i == geom->num_subparts - 1)
                {
                    if (NULL != geom->suboffsets)
                    {
                        tmp_num_points = geom->num_points - geom->suboffsets[i];
                    }
                    else
                        tmp_num_points = geom->num_points;
                }
                else
                    tmp_num_points = geom->suboffsets[i+1] - geom->suboffsets[i];

                put_integer (&cp, 0, 1, &tmp_num_points);
                put_xy (&cp, 0, tmp_num_points, &geom->pt[geom->suboffsets[i]]);
            }
            break;

        case geomMultiPolygon:
            type = wkbMultiPolygon;
            put_integer (&cp, 0, 1, (int *)&type);
            put_integer (&cp, 0, 1, &geom->num_parts);
            type = wkbPolygon; /* Each part is a wkb polygon */

            for (j = 0; j < geom->num_parts; j++)
            {
                *cp = byte_order;
                cp++;
                put_integer (&cp, 0, 1, (int *)&type);

                if (j == geom->num_parts - 1)
                    end_subparts = geom->num_subparts;
                else
                    end_subparts = geom->offsets[j+1];
                tmp_num_subparts = end_subparts - geom->offsets[j];
                put_integer (&cp, 0, 1, &tmp_num_subparts);

                for (i = geom->offsets[j]; i < end_subparts; i++)
                {
                    if (i == geom->num_subparts - 1)
                        tmp_num_points = geom->num_points - geom->suboffsets[i];
                    else
                        tmp_num_points = geom->suboffsets[i+1] - geom->suboffsets[i];
  
                    put_integer (&cp, 0, 1, &tmp_num_points);
                    put_xy (&cp, 0, tmp_num_points, &geom->pt[geom->suboffsets[i]]);
                }
            }
            break;

    }
}



/***************************************************************************
*
*  _calc_wkblen  -  Calculates the number of bytes in the well-known
*                   binary representation of a geometry.
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Purpose:
*    Calculate the number of bytes in the shape.
*  
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Parameters:
*    type          <Input>  == (geomType) The geometry data type.
*    num_points    <Input>  == (int) The number of points in the geometry.
*    num_parts     <Input>  == (int) The number of parts in the geometry.
*    num_subparts  <Input>  == (int) The number of subparts in the geometry.
*
*    RETURN        <Output> == (int) The required binary size
*
***************************************************************************/

static int _calc_wkblen (
    GeomType type,
    int      num_points,
    int      num_parts,
    int      num_subparts
    )
{

    int data_len;
 
    data_len = 1 + sizeof (int);  /* Add the size of the data type */

    switch (type)      
    {
        case geomPoint:
            data_len += (2 * sizeof (double));
            break;

        case geomMultiPoint:
            /*
             *   num_points + num_points * (byte order + type + (x,y))
             */
            data_len += sizeof(int);
            data_len += num_points * (1 + sizeof(int) + 2 * sizeof(double));
            break;

        case geomLineString:
            /*
             *   num_points + num_points * (x,y)
             */
            data_len += sizeof(int);
            data_len += num_points * (2 * sizeof(double));
            break;

        case geomMultiLineString:
            /*
             *    num_parts + num_parts * (byte order + type + num_points)
             *  + num_points * (x,y)
             */
            data_len += sizeof(int);
            data_len += num_parts * (1 + 2 * sizeof(int));
            data_len += num_points * (2 * sizeof(double));
            break;

        case geomPolygon:
            /*
             *    num_subparts + num_subparts * (num_points)
             *  + num_points * (x,y)
             */
            data_len += sizeof(int);
            data_len += num_subparts * sizeof (int);
            data_len += num_points * (2 * sizeof(double));
            break;

        case geomMultiPolygon:
            /*
             *    num_parts + num_parts * (byte order + type + num_subparts)
             *  + num_subparts * (num_points) + num_points * (x,y)
             */
            data_len += sizeof(int);
            data_len += num_parts * (1 + 2 * sizeof (int));
            data_len += num_subparts * sizeof (int);
            data_len += num_points * (2 * sizeof(double));
            break;
    }

    return data_len;
}


/***************************************************************************
*
*  _get_wkb_header  -  Gets the wkb byte order and shape type
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Purpose:
*    This function extracts the byte order of the OGIS WKB
*    byte stream, as well as the shape type.
*    The pointer into the buffer will be advanced to the position
*    immediately following the shape type.
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Parameters:
*    buf_ptr    <In/Out> == (char **) Address of pointer into byte stream
*    swap_bytes <Output> == (int) Do bytes need to be swapped
*    wkb_type   <Output> == (WkbType *) OGIS shape type
*
*    RETURN     <Output> ==  (int) Error code
*    
***************************************************************************/

static int _get_wkb_header (
    char    **buf_ptr, 
    int      *swap_bytes,
    WkbType  *wkb_type
    )
{
    int       i;
    ByteOrder byte_order;
    
    switch ((*buf_ptr)[0])
    {
        case 0:
            byte_order = BigEndian;
            break;

        case 1:
            byte_order = LittleEndian;
            break;

        default:
            return GEOM_INVALID_WKB_BYTE_ORDER;
    }

    ENDIAN_TEST (i);  /* Sets i = 0 if big-endian, 1 if little-endian */
    *swap_bytes = 1 - i;

    (*buf_ptr)++;
    get_integer (buf_ptr, *swap_bytes, 1, (int *) wkb_type);

    return GEOM_SUCCESS;
}


/***************************************************************************
*
*   wkb_to_geom  -  Constructs a geometry from an OGIS byte stream
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Purpose:
*    This function builds a shape from the OGIS compliant Well Known Binary 
*    representation of a geometry.
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Parameters:
*     binary       <Input> == (char *) The geometry's wkb representation 
*     geom        <Output> == (Geometry *) A geometry object
*
*     RETURN      <Output> == (int) Error code
*
***************************************************************************/

int wkb_to_geom (
    char     *binary,
    Geometry *geom
    )
{
    WkbType   wkb_type;
    int       rc;
    char     *cp;
    int       swap_bytes;

    /* Get the byte order and geometry type */
    cp = binary;
    rc = _get_wkb_header (&cp, &swap_bytes, &wkb_type);
    if (rc != GEOM_SUCCESS)
        return rc;
    
    switch (wkb_type)
    {
        /* copy the points based on shape type */

        case wkbPoint:
            geom->type = geomPoint;
            geom->num_points = 1;

            /* Allocate space for the shape point.  */
            rc = geom_allocate(geom, 1);
            if (rc != GEOM_SUCCESS)
                return rc;

            /* Get an xy point */
            get_xy (&cp, swap_bytes, 1, geom->pt);
            break;

        case wkbMultiPoint:
            geom->type = geomMultiPoint;
            _get_wkb_points (cp, swap_bytes, geom);
            break;

        case wkbLineString:
            geom->type = geomLineString;
            _get_wkb_lines (cp, swap_bytes, geom);
            break;

        case wkbMultiLineString:
            geom->type = geomMultiLineString;
            _get_wkb_lines (cp, swap_bytes, geom);
            break;
  
        case wkbPolygon:
            geom->type = geomPolygon;
            _get_wkb_polys (cp, swap_bytes, geom);
            break;
  
        case wkbMultiPolygon:
            geom->type = geomMultiPolygon;
            _get_wkb_polys (cp, swap_bytes, geom);
            break;

        case wkbCollection:
            return GEOM_UNSUPPORTED_WKB_TYPE;
            break;
            
        default:
            return GEOM_INVALID_WKB_TYPE;
    }

    return GEOM_SUCCESS;
}


/***************************************************************************
*
*  _get_wkb_points  -  Extracts a point or multi shape from the WKB binary
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Purpose:
*    This function builds a shape from the WKB binary point or multipoint
*    shape representation of a shape.
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Parameters:
*    swap_bytes   <Input>  == (int) Do bytes need to be swapped
*    cp           <In/Out> == (char *) a pointer into the binary buffer
*    geom         <Output> == (Geometry *) A geometry object
*
*    RETURN       <Output> == (int) Error code
*
***************************************************************************/

static int _get_wkb_points (
    char     *cp,
    int       swap_bytes,
    Geometry *geom
    )
{
    int       rc;
    int       npoints, i;
    WkbType   wkb_type;

    get_integer (&cp, swap_bytes, 1, &npoints);

    if (npoints <= 0)
        return GEOM_TOO_FEW_POINTS;
  
    /* Allocate space for the geometry points. */

    rc = geom_allocate(geom, npoints);
    if (rc != GEOM_SUCCESS)
        return rc;

    geom->num_points = npoints;
    
    /* Loop through the list of wkbPoints */
    for (i = 0; i < npoints; i++)
    {
        rc = _get_wkb_header (&cp, &swap_bytes, &wkb_type);

        if (rc != GEOM_SUCCESS)
            return rc;
        if (wkb_type != wkbPoint)
            return GEOM_INVALID_WKB_TYPE;

        /* Get an xy point */
        get_xy (&cp, swap_bytes, 1, &geom->pt[i]);
    }

    return GEOM_SUCCESS;
}

/***************************************************************************
*
*  _get_wkb_lines  -  Extracts polyline parts from the OGIS binary
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Purpose:
*    This function builds a polyline shape from the OGIS Well
*    Known binary representation of a shape. 
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Parameters:
*    swap_bytes     <Input> == (int) Do bytes need to be swapped
*    cp            <In/Out> == (char *) a pointer into the binary buffer
*    geom          <In/Out> == (Geometry *) Handle to a geometry object.
*
*    RETURN        <Output> == (int) Error code
*
***************************************************************************/

static int _get_wkb_lines (
    char     *cp,
    int       swap_bytes,
    Geometry *geom
    )
{
    int       rc;
    int       current_pt = 0;
    int       i, npoints, nlines;
    WkbType   wkb_type;

    if (geom->type == geomMultiLineString)
    {
        /* Get the number of line elements */
        get_integer (&cp, swap_bytes, 1, &nlines);
        if (nlines < 1)
            return GEOM_INVALID_NUM_PARTS;

        /* Get the first line's header */
        rc = _get_wkb_header (&cp, &swap_bytes, &wkb_type);
        if (rc != GEOM_SUCCESS)
            return rc;
        
        if (wkb_type != wkbLineString)
            return GEOM_INVALID_WKB_TYPE;
    }
    else
    {
        nlines = 1;
    }

    geom->num_points   = 0;
    geom->num_parts    = nlines;
    geom->num_subparts = 0;
    geom->offsets      = malloc (nlines * sizeof(int));
    
    /* Loop through each part and extract the points. */

    for (i = 0; i < nlines; i++)
    {
        geom->offsets[i] = current_pt;

        /* Get the number of points for this part */
        get_integer (&cp, swap_bytes, 1, &npoints);
        if (npoints < 2)
            return GEOM_TOO_FEW_POINTS;

        /* Allocate space for the geometry points for this line */
        rc = geom_allocate (geom, geom->num_points + npoints); 
        if (rc != GEOM_SUCCESS)
            return rc;
        geom->num_points += npoints;

        /* Now, get the points for the line */
        get_xy (&cp, swap_bytes, npoints, &geom->pt[current_pt]);
        current_pt += npoints;

        if (i < nlines - 1)
        {
            /* Read the next line's header */
            rc = _get_wkb_header (&cp, &swap_bytes, &wkb_type);
            if (rc != GEOM_SUCCESS)
                return rc;
            if (wkb_type != wkbLineString)
                return GEOM_INVALID_WKB_TYPE;
        }
    } /* end for */
    
  return GEOM_SUCCESS;
}

/***************************************************************************
*
*  _get_wkb_polys  -  Extracts polygon parts from the OGIS binary
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Purpose:
*    This function builds a polygon shape from the OGIS Well
*    Known binary representation of a shape. 
*
*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
*
*  Parameters:
*    cp            <In/Out> == (char *) a pointer into the binary buffer
*    swap_bytes    <Input>  == (int) Do bytes need to be swapped
*    geom          <Output> == (Geometry *) A geometry object
*
*    RETURN        <Output> == (int) Error code
*
***************************************************************************/

static int _get_wkb_polys (
    char     *cp,
    int       swap_bytes,
    Geometry *geom
    )
{
    int       current_pt = 0;
    int       rc;
    int       i, j, k;
    int       npoints, npolys, nrings;
    WkbType   wkb_type;

    if (geom->type == geomMultiPolygon)
    {
        /* Get the number of polygon elements */
        get_integer (&cp, swap_bytes, 1, &npolys);
        if (npolys < 1)
            return GEOM_INVALID_NUM_PARTS;

        /* Get the first polygon's header */
        rc = _get_wkb_header (&cp, &swap_bytes, &wkb_type);
        if (rc != GEOM_SUCCESS)
            return rc;
    
        if (wkb_type != wkbPolygon)
            return GEOM_INVALID_WKB_TYPE;
    }
    else
    {
        npolys = 1;
    }
  
    geom->num_points    = 0;
    geom->num_parts     = npolys;
    geom->num_subparts  = 0;
    geom->offsets       = malloc (npolys * sizeof(int));
    geom->suboffsets    = malloc (3 * npolys * sizeof(int));
    geom->sub_allocsize = 3 * npolys;
    
    k = 0;  /* suboffset array index */
    
    /* Loop through each part and extract the points. */
    for (i = 0; i < npolys; i++)
    {
        geom->offsets[i] = current_pt;
        
        /* Get the number of rings for this polygon */
        get_integer (&cp, swap_bytes, 1, &nrings);
        if (nrings < 1)
            return GEOM_INVALID_NUM_PARTS;

        /* Reallocate & grow the suboffsets array if necessary */
        if (geom->num_subparts + nrings >= geom->sub_allocsize)
        {
            int *newArray;
            geom->sub_allocsize = geom->num_subparts + 2 * nrings;
            newArray = malloc (geom->sub_allocsize * sizeof(int));
            memcpy (newArray, geom->suboffsets,
                    geom->num_subparts * sizeof(int));
            free (geom->suboffsets);
            geom->suboffsets = newArray;
        }

        geom->num_subparts += nrings;

        /* Loop through each ring to get the points */
        for (j = 0; j < nrings; j++)
        {
            geom->suboffsets[k] = current_pt;
            k++;
            
            /* Get the number of points for this ring */
            get_integer (&cp, swap_bytes, 1, &npoints);
            if (npoints < 4)
                return GEOM_TOO_FEW_POINTS;
            else
            {
                /* Allocate space for the shape points. */
                rc = geom_allocate (geom, geom->num_points + npoints);
                if (rc != GEOM_SUCCESS)
                    return rc;
                geom->num_points += npoints;


                /* Now, get the points for the ring */
                get_xy (&cp, swap_bytes, npoints, &geom->pt[current_pt]);
                current_pt += npoints;
            }
        } /* end ring loop */

        if (i < npolys - 1)
        {
            /* read the next polygon's header */
            rc = _get_wkb_header (&cp, &swap_bytes, &wkb_type);
            if (rc != GEOM_SUCCESS)
                return rc;

            if (wkb_type != wkbPolygon)
                return GEOM_INVALID_WKB_TYPE;
        }
    } /* end polygon loop */
    
  return GEOM_SUCCESS;
}