Index: mappostgis.c =================================================================== --- mappostgis.c (revision 10768) +++ mappostgis.c (working copy) @@ -7,7 +7,7 @@ * Dave Blasby * ****************************************************************************** - * Copyright (c) 2008 Paul Ramsey + * Copyright (c) 2010 Paul Ramsey * Copyright (c) 2002 Refractions Research * * Permission is hereby granted, free of charge, to any person obtaining a @@ -66,10 +66,19 @@ #define FLT_MAX 25000000.0 #endif +#define FP_EPSILON 1e-12 +#define FP_EQ(a, b) (fabs((a)-(b)) < FP_EPSILON) +#define FP_LEFT -1 +#define FP_RIGHT 1 +#define FP_COLINEAR 0 + +#define SEGMENT_ANGLE 15.0 + #ifdef USE_POSTGIS -MS_CVSID("$Id$") +MS_CVSID("$Id$"); + /* ** msPostGISCloseConnection() ** @@ -130,232 +139,689 @@ } } -/* -** TODO, review and clean -*/ -static int force_to_points(unsigned char *wkb, shapeObj *shape) +pointArrayObj* +pointArrayNew(int maxpoints) { - /* we're going to make a 'line' for each entity (point, line or ring) in the geom collection */ + pointArrayObj *d = malloc(sizeof(pointArrayObj)); + d->data = malloc(maxpoints * sizeof(pointObj)); + d->npoints = 0; + d->maxpoints = maxpoints; + return d; +} - int offset = 0; - int pt_offset; - int ngeoms; - int t, u, v; - int type, nrings, npoints; - lineObj line = {0, NULL}; +void +pointArrayFree(pointArrayObj *d) +{ + if ( ! d ) return; + if ( d->data ) free(d->data); + free(d); +} - shape->type = MS_SHAPE_NULL; /* nothing in it */ +static int +pointArrayAddPoint(pointArrayObj *d, const pointObj *p) +{ + if ( !p || !d ) return MS_FAILURE; + /* Avoid overwriting memory buffer */ + if ( d->maxpoints - d->npoints == 0 ) { + d->maxpoints *= 2; + d->data = realloc(d->data, d->maxpoints * sizeof(pointObj)); + } + d->data[d->npoints] = *p; + d->npoints++; + return MS_SUCCESS; +} - memcpy(&ngeoms, &wkb[5], 4); - offset = 9; /* where the first geometry is */ - - for(t=0; t < ngeoms; t++) { - memcpy(&type, &wkb[offset + 1], 4); /* type of this geometry */ - if(type == 1) { - /* Point */ - shape->type = MS_SHAPE_POINT; - line.numpoints = 1; - line.point = (pointObj *) malloc(sizeof(pointObj)); +static int +wkbType(wkbObj *w) +{ + int t; + memcpy(&t, (w->ptr + 1), sizeof(int)); + return t; +} - memcpy(&line.point[0].x, &wkb[offset + 5], 8); - memcpy(&line.point[0].y, &wkb[offset + 5 + 8], 8); - offset += 5 + 16; - msAddLine(shape, &line); - free(line.point); - } else if(type == 2) { - /* Linestring */ - shape->type = MS_SHAPE_POINT; +static int +wkbCollectionSubType(wkbObj *w) +{ + int t; + memcpy(&t, (w->ptr + 1 + 4 + 4 + 1), sizeof(int)); + return t; +} - memcpy(&line.numpoints, &wkb[offset+5], 4); /* num points */ - line.point = (pointObj *) malloc(sizeof(pointObj) * line.numpoints); - for(u = 0; u < line.numpoints; u++) { - memcpy( &line.point[u].x, &wkb[offset+9 + (16 * u)], 8); - memcpy( &line.point[u].y, &wkb[offset+9 + (16 * u)+8], 8); - } - offset += 9 + 16 * line.numpoints; /* length of object */ - msAddLine(shape, &line); - free(line.point); - } else if(type == 3) { - /* Polygon */ - shape->type = MS_SHAPE_POINT; - memcpy(&nrings, &wkb[offset+5],4); /* num rings */ - /* add a line for each polygon ring */ - pt_offset = 0; - offset += 9; /* now points at 1st linear ring */ - for(u = 0; u < nrings; u++) { - /* for each ring, make a line */ - memcpy(&npoints, &wkb[offset], 4); /* num points */ - line.numpoints = npoints; - line.point = (pointObj *) malloc(sizeof(pointObj)* npoints); - /* point struct */ - for(v = 0; v < npoints; v++) - { - memcpy(&line.point[v].x, &wkb[offset + 4 + (16 * v)], 8); - memcpy(&line.point[v].y, &wkb[offset + 4 + (16 * v) + 8], 8); - } - /* make offset point to next linear ring */ - msAddLine(shape, &line); - free(line.point); - offset += 4 + 16 *npoints; - } - } - } - - return MS_SUCCESS; +static char +wkbReadChar(wkbObj *w) +{ + char c; + memcpy(&c, w->ptr, sizeof(char)); + w->ptr += sizeof(char); + return c; } -/* convert the wkb into lines */ -/* points-> remove */ -/* lines -> pass through */ -/* polys -> treat rings as lines */ +static inline int +wkbReadInt(wkbObj *w) +{ + int i; + memcpy(&i, w->ptr, sizeof(int)); + w->ptr += sizeof(int); + //printf(" int %d\n", i); + return i; +} -static int force_to_lines(unsigned char *wkb, shapeObj *shape) +static inline double +wkbReadDouble(wkbObj *w) { - int offset = 0; - int pt_offset; - int ngeoms; - int t, u, v; - int type, nrings, npoints; - lineObj line = {0, NULL}; + double d; + memcpy(&d, w->ptr, sizeof(double)); + w->ptr += sizeof(double); + //printf(" dbl %.5g\n", d); + return d; +} - shape->type = MS_SHAPE_NULL; /* nothing in it */ +static inline void +wkbReadPointP(wkbObj *w, pointObj *p) +{ + memcpy(&(p->x), w->ptr, sizeof(double)); + //printf(" pnt %.5g ", p->x); + w->ptr += sizeof(double); + memcpy(&(p->y), w->ptr, sizeof(double)); + //printf("%.5g\n", p->y); + w->ptr += sizeof(double); +} - memcpy(&ngeoms, &wkb[5], 4); - offset = 9; /* were the first geometry is */ - for(t=0; t < ngeoms; t++) { - memcpy(&type, &wkb[offset + 1], 4); /* type of this geometry */ +static inline pointObj +wkbReadPoint(wkbObj *w) +{ + pointObj p; + wkbReadPointP(w, &p); + return p; +} - /* cannot do anything with a point */ +static lineObj* +wkbReadLine(wkbObj *w) +{ + int i; + pointObj p; + lineObj *line = malloc(sizeof(lineObj)); + int npoints = wkbReadInt(w); - if(type == 2) { - /* Linestring */ - shape->type = MS_SHAPE_LINE; + line->numpoints = npoints; + line->point = malloc(npoints * sizeof(pointObj)); + for ( i = 0; i < npoints; i++ ) { + wkbReadPointP(w, &p); + line->point[i] = p; + } + return line; +} - memcpy(&line.numpoints, &wkb[offset + 5], 4); - line.point = (pointObj*) malloc(sizeof(pointObj) * line.numpoints ); - for(u=0; u < line.numpoints; u++) { - memcpy(&line.point[u].x, &wkb[offset + 9 + (16 * u)], 8); - memcpy(&line.point[u].y, &wkb[offset + 9 + (16 * u)+8], 8); - } - offset += 9 + 16 * line.numpoints; /* length of object */ - msAddLine(shape, &line); - free(line.point); - } else if(type == 3) { - /* polygon */ - shape->type = MS_SHAPE_LINE; +static void +wkbSkipGeometry(wkbObj *w) +{ + char endian; + int type, npoints, nrings, ngeoms, i; + endian = wkbReadChar(w); + type = wkbReadInt(w); + switch(type) { + case WKB_POINT: + w->ptr += 2 * sizeof(double); + break; + case WKB_CIRCULARSTRING: + case WKB_LINESTRING: + npoints = wkbReadInt(w); + w->ptr += npoints * 2 * sizeof(double); + break; + case WKB_POLYGON: + nrings = wkbReadInt(w); + for ( i = 0; i < nrings; i++ ) { + npoints = wkbReadInt(w); + w->ptr += npoints * 2 * sizeof(double); + } + break; + case WKB_MULTIPOINT: + case WKB_MULTILINESTRING: + case WKB_MULTIPOLYGON: + case WKB_GEOMETRYCOLLECTION: + case WKB_COMPOUNDCURVE: + case WKB_CURVEPOLYGON: + case WKB_MULTICURVE: + case WKB_MULTISURFACE: + ngeoms = wkbReadInt(w); + for ( i = 0; i < ngeoms; i++ ) { + wkbSkipGeometry(w); + } + } +} - memcpy(&nrings, &wkb[offset + 5], 4); /* num rings */ - /* add a line for each polygon ring */ - pt_offset = 0; - offset += 9; /* now points at 1st linear ring */ - for(u = 0; u < nrings; u++) { - /* for each ring, make a line */ - memcpy(&npoints, &wkb[offset], 4); - line.numpoints = npoints; - line.point = (pointObj*) malloc(sizeof(pointObj) * npoints); - /* point struct */ - for(v = 0; v < npoints; v++) { - memcpy(&line.point[v].x, &wkb[offset + 4 + (16 * v)], 8); - memcpy(&line.point[v].y, &wkb[offset + 4 + (16 * v) + 8], 8); - } - /* make offset point to next linear ring */ - msAddLine(shape, &line); - free(line.point); - offset += 4 + 16 * npoints; - } - } - } +int +wkbConvPointToShape(wkbObj *w, shapeObj *shape) +{ + char endian; + int type; + lineObj *line; - return MS_SUCCESS; + endian = wkbReadChar(w); + type = wkbReadInt(w); + + if( type != WKB_POINT ) return MS_FAILURE; + + if( ! (shape->type == MS_SHAPE_POINT) ) return MS_FAILURE; + line = malloc(sizeof(lineObj)); + line->numpoints = 1; + line->point = malloc(sizeof(pointObj)); + line->point[0] = wkbReadPoint(w); + msAddLineDirectly(shape, line); + return MS_SUCCESS; } -/* point -> reject */ -/* line -> reject */ -/* polygon -> lines of linear rings */ -static int force_to_polygons(unsigned char *wkb, shapeObj *shape) +int +wkbConvLineStringToShape(wkbObj *w, shapeObj *shape) { - int offset = 0; - int pt_offset; - int ngeoms; - int t, u, v; - int type, nrings, npoints; - lineObj line = {0, NULL}; + char endian; + int type; - shape->type = MS_SHAPE_NULL; /* nothing in it */ + endian = wkbReadChar(w); + type = wkbReadInt(w); - memcpy(&ngeoms, &wkb[5], 4); - offset = 9; /* were the first geometry is */ - for(t = 0; t < ngeoms; t++) { - memcpy(&type, &wkb[offset + 1], 4); /* type of this geometry */ + if( type != WKB_LINESTRING ) return MS_FAILURE; + + msAddLineDirectly(shape, wkbReadLine(w)); + + return MS_SUCCESS; +} - if(type == 3) { - /* polygon */ - shape->type = MS_SHAPE_POLYGON; +int +wkbConvPolygonToShape(wkbObj *w, shapeObj *shape) +{ + char endian; + int type; + int i, nrings; - memcpy(&nrings, &wkb[offset + 5], 4); /* num rings */ - /* add a line for each polygon ring */ - pt_offset = 0; - offset += 9; /* now points at 1st linear ring */ - for(u=0; u < nrings; u++) { - /* for each ring, make a line */ - memcpy(&npoints, &wkb[offset], 4); /* num points */ - line.numpoints = npoints; - line.point = (pointObj*) malloc(sizeof(pointObj) * npoints); - for(v=0; v < npoints; v++) { - memcpy(&line.point[v].x, &wkb[offset + 4 + (16 * v)], 8); - memcpy(&line.point[v].y, &wkb[offset + 4 + (16 * v) + 8], 8); - } - /* make offset point to next linear ring */ - msAddLine(shape, &line); - free(line.point); - offset += 4 + 16 * npoints; - } - } - } + endian = wkbReadChar(w); + type = wkbReadInt(w); + + if( type != WKB_POLYGON ) return MS_FAILURE; - return MS_SUCCESS; + /* How many rings? */ + nrings = wkbReadInt(w); + + /* Add each ring to the shape */ + for( i = 0; i < nrings; i++ ) + msAddLineDirectly(shape, wkbReadLine(w)); + + return MS_SUCCESS; } +int +wkbConvCurvePolygonToShape(wkbObj *w, shapeObj *shape) +{ + char endian; + int type, i, ncomponents; + int failures = 0; + int was_poly = ( shape->type == MS_SHAPE_POLYGON ); + + endian = wkbReadChar(w); + type = wkbReadInt(w); + ncomponents = wkbReadInt(w); + + /* Lower the allowed dimensionality so we can + * catch the linear ring components */ + shape->type = MS_SHAPE_LINE; + + for ( i = 0; i < ncomponents; i++ ) { + if ( wkbConvGeometryToShape(w, shape) == MS_FAILURE ) { + wkbSkipGeometry(w); + failures++; + } + } + + /* Go back to expected dimensionality */ + if ( was_poly) shape->type = MS_SHAPE_POLYGON; + + if ( failures == ncomponents ) + return MS_FAILURE; + else + return MS_SUCCESS; +} + +int +wkbConvCircularStringToShape(wkbObj *w, shapeObj *shape) +{ + char endian; + int type; + lineObj line; + + endian = wkbReadChar(w); + type = wkbReadInt(w); + + if( type != WKB_CIRCULARSTRING ) return MS_FAILURE; + + /* Stroke the string into a point array */ + if ( arcStrokeCircularString(w, SEGMENT_ANGLE, &line) == MS_FAILURE ) + return MS_FAILURE; + + /* Fill in the lineObj */ + if ( line.numpoints > 0 ) + msAddLine(shape, &line); + + return MS_SUCCESS; +} + +int +wkbConvCompoundCurveToShape(wkbObj *w, shapeObj *shape) +{ + char endian; + int npoints = 0; + int type, ncomponents, i, j; + lineObj *line; + shapeObj shapebuf; + + endian = wkbReadChar(w); + type = wkbReadInt(w); + + /* Init our shape buffer */ + msInitShape(&shapebuf); + + if( type != WKB_COMPOUNDCURVE ) return MS_FAILURE; + + //printf("type %d\n", type); + + /* How many components in the compound curve? */ + ncomponents = wkbReadInt(w); + + //printf("ncomponents %d\n", ncomponents); + + /* We'll load each component onto a line in a shape */ + for( i = 0; i < ncomponents; i++ ) { + int type = wkbType(w); + //printf(" subtype %d\n", type); + if ( type == WKB_CIRCULARSTRING ) + wkbConvCircularStringToShape(w, &shapebuf); + else if ( type == WKB_LINESTRING ) + wkbConvLineStringToShape(w, &shapebuf); + else + return MS_FAILURE; + } + + /* Do nothing on empty */ + if ( shapebuf.numlines == 0 ) + return MS_FAILURE; + + /* Count the total number of points */ + for( i = 0; i < shapebuf.numlines; i++ ) { + //printf("shapebuf.line[%d].numpoints %d\n",i,shapebuf.line[i].numpoints); + npoints += shapebuf.line[i].numpoints; + } + //printf("npoints %d\n",npoints); + + /* Do nothing on empty */ + if ( npoints == 0 ) + return MS_FAILURE; + + /* Allocate space for the new line */ + line = malloc(sizeof(lineObj)); + line->numpoints = npoints; + line->point = malloc(sizeof(pointObj) * npoints); + + /* Copy in the points */ + npoints = 0; + //printf("shapebuf.numlines %d\n", shapebuf.numlines); + for ( i = 0; i < shapebuf.numlines; i++ ) { + //printf("shapebuf.line[%d].numpoints %d\n",i,shapebuf.line[i].numpoints); + for ( j = 0; j < shapebuf.line[i].numpoints; j++ ) { + /* Don't add a start point that duplicates an endpoint */ + if( j == 0 && i > 0 && + memcmp(&(line->point[npoints - 1]),&(shapebuf.line[i].point[j]),sizeof(pointObj)) == 0 ) { + continue; + } + line->point[npoints++] = shapebuf.line[i].point[j]; + } + } + line->numpoints = npoints; + + /* Clean up */ + msFreeShape(&shapebuf); + + /* Fill in the lineObj */ + msAddLineDirectly(shape, line); + + return MS_SUCCESS; +} + +int +wkbConvCollectionToShape(wkbObj *w, shapeObj *shape) +{ + char endian; + int type, i, ncomponents; + int failures = 0; + + endian = wkbReadChar(w); + type = wkbReadInt(w); + ncomponents = wkbReadInt(w); + + for ( i = 0; i < ncomponents; i++ ) { + if ( wkbConvGeometryToShape(w, shape) == MS_FAILURE ) { + wkbSkipGeometry(w); + failures++; + } + } + if ( failures == ncomponents ) + return MS_FAILURE; + else + return MS_SUCCESS; +} + +int +wkbConvGeometryToShape(wkbObj *w, shapeObj *shape) +{ + int wkbtype = wkbType(w); /* Peak at the type number */ + //printf ("wkbtype %d\n",wkbtype); + /* Unwind collections into primatives */ + if ( wkbtype == WKB_GEOMETRYCOLLECTION ) return wkbConvCollectionToShape(w, shape); + + /* Handle area types */ + if ( wkbtype == WKB_POLYGON ) return wkbConvPolygonToShape(w, shape); + if ( wkbtype == WKB_MULTIPOLYGON ) return wkbConvCollectionToShape(w, shape); + if ( wkbtype == WKB_CURVEPOLYGON ) return wkbConvCurvePolygonToShape(w, shape); + if ( wkbtype == WKB_MULTISURFACE ) return wkbConvCollectionToShape(w, shape); + + /* We can't convert any of the following types into polygons */ + if ( shape->type == MS_SHAPE_POLYGON ) return MS_FAILURE; + + /* Handle linear types */ + if ( wkbtype == WKB_LINESTRING ) return wkbConvLineStringToShape(w, shape); + if ( wkbtype == WKB_CIRCULARSTRING ) return wkbConvCircularStringToShape(w, shape); + if ( wkbtype == WKB_COMPOUNDCURVE ) return wkbConvCompoundCurveToShape(w, shape); + if ( wkbtype == WKB_MULTILINESTRING ) return wkbConvCollectionToShape(w, shape); + if ( wkbtype == WKB_MULTICURVE ) return wkbConvCollectionToShape(w, shape); + + /* We can't convert any of the following types into lines */ + if ( shape->type == MS_SHAPE_LINE ) return MS_FAILURE; + + /* Handle point types */ + if ( wkbtype == WKB_POINT ) return wkbConvPointToShape(w, shape); + if ( wkbtype == WKB_MULTIPOINT ) return wkbConvCollectionToShape(w, shape); + + /* This is a type we don't know about! */ + return MS_FAILURE; +} + + +/* +** Calculate determinant of a 3x3 matrix. Handy for +** the circle center calculation. +*/ +static inline double +arcDeterminant3x3(double *m) +{ + /* This had better be a 3x3 matrix or we'll fall to bits */ + return m[0] * ( m[4] * m[8] - m[7] * m[5] ) - + m[3] * ( m[1] * m[8] - m[7] * m[2] ) + + m[6] * ( m[1] * m[5] - m[4] * m[2] ); +} + +/* +** What side of p1->p2 is q on? +*/ +static inline int +arcSegmentSide(const pointObj *p1, const pointObj *p2, const pointObj *q) +{ + double side = ( (q->x - p1->x) * (p2->y - p1->y) - (p2->x - p1->x) * (q->y - p1->y) ); + if ( FP_EQ(side,0.0) ) { + return FP_COLINEAR; + } + else { + if ( side < 0.0 ) + return FP_LEFT; + else + return FP_RIGHT; + } +} + +/* +** Calculate the center of the circle defined by three points. +** Using matrix approach from http://mathforum.org/library/drmath/view/55239.html +*/ +int +arcCircleCenter(const pointObj *p1, const pointObj *p2, const pointObj *p3, pointObj *center, double *radius) +{ + pointObj c; + double r; + + /* Components of the matrices. */ + double x1sq = p1->x * p1->x; + double x2sq = p2->x * p2->x; + double x3sq = p3->x * p3->x; + double y1sq = p1->y * p1->y; + double y2sq = p2->y * p2->y; + double y3sq = p3->y * p3->y; + double matrix_num_x[9] = { x1sq+y1sq, p1->y, 1.0, x2sq+y2sq, p2->y, 1.0, x3sq+y3sq, p3->y, 1.0 }; + double matrix_num_y[9] = { p1->x, x1sq+y1sq, 1.0, p2->x, x2sq+y2sq, 1.0, p3->x, x3sq+y3sq, 1.0 }; + double matrix_denom[9] = { p1->x, p1->y, 1.0, p2->x, p2->y, 1.0, p3->x, p3->y, 1.0 }; + + /* Circle is closed, so p1 must be opposite p1 & p3. */ + if ( FP_EQ(p1->x,p3->x) && FP_EQ(p1->y,p3->y) ) { + c.x = (p1->x + p2->x) / 2.0; + c.y = (p1->y + p2->y) / 2.0; + r = sqrt( (p1->x - p2->x) * (p1->x - p2->x) + (p1->y - p2->y) * (p1->y - p2->y) ) / 2.0; + } + /* There is no circle here, the points are actually co-linear */ + else if ( arcSegmentSide(p1, p3, p2) == FP_COLINEAR ) { + return MS_FAILURE; + } + /* Calculate the center and radius. */ + else { + double denom = 2.0 * arcDeterminant3x3(matrix_denom); + /* Center components */ + c.x = arcDeterminant3x3(matrix_num_x) / denom; + c.y = arcDeterminant3x3(matrix_num_y) / denom; + + /* Radius */ + r = sqrt((p1->x-c.x) * (p1->x-c.x) + (p1->y-c.y) * (p1->y-c.y)); + } + + if ( radius ) *radius = r; + if ( center ) *center = c; + + return MS_SUCCESS; +} + +/* +** Write a stroked version of the circle defined by three points into a +** point buffer. The segment_angle (degrees) is the coverage of each stroke segment, +** and depending on whether this is the first arc in a circularstring, +** you might want to include_first +*/ +int +arcStrokeCircle(const pointObj *p1, const pointObj *p2, const pointObj *p3, + double segment_angle, int include_first, pointArrayObj *pa) +{ + pointObj center; /* Center of our circular arc */ + double radius; /* Radius of our circular arc */ + double sweep_angle_r; /* Total angular size of our circular arc in radians */ + double segment_angle_r; /* Segment angle in radians */ + double a1, a2, a3; /* Angles represented by p1, p2, p3 relative to center */ + int side = arcSegmentSide(p1, p3, p2); /* What side of p1,p3 is the middle point? */ + int num_edges; /* How many edges we will be generating */ + double current_angle_r; /* What angle are we generating now (radians)? */ + int i; /* Counter */ + pointObj p; /* Temporary point */ + + //printf("side = %d\n",side); + + if ( side == FP_COLINEAR ) { + /* Co-linear case, we just need to write in the end points */ + if ( include_first ) + pointArrayAddPoint(pa, p1); + pointArrayAddPoint(pa, p3); + return MS_SUCCESS; + } + + if ( arcCircleCenter(p1, p2, p3, ¢er, &radius) == MS_FAILURE ) + return MS_FAILURE; + + a1 = atan2(p1->y - center.y, p1->x - center.x); + a2 = atan2(p2->y - center.y, p2->x - center.x); + a3 = atan2(p3->y - center.y, p3->x - center.x); + segment_angle_r = M_PI * segment_angle / 180.0; + + //printf("a1 = %g\n",a1); + //printf("a2 = %g\n",a2); + //printf("a3 = %g\n",a3); + + /* Closed-circle case */ + if ( FP_EQ(p1->x, p3->x) && FP_EQ(p1->y, p3->y) ) { + sweep_angle_r = 2.0 * M_PI; + } + /* Clockwise arc direction */ + else if ( side == FP_LEFT ) { + //printf ("clockwise!\n"); + if ( a3 > a1 ) /* Wrapping past 180? */ + sweep_angle_r = a1 + (2.0 * M_PI - a3); + else + sweep_angle_r = a1 - a3; + } + /* Counter-clockwise arc direction */ + else if ( side == FP_RIGHT ) { + //printf ("counter clockwise!\n"); + if ( a3 > a1 ) /* Wrapping past 180? */ + sweep_angle_r = a3 - a1; + else + sweep_angle_r = a3 + (2.0 * M_PI - a1); + } + else + sweep_angle_r = 0.0; + + //printf("sweep_angle_r = %g\n", sweep_angle_r); + //printf("segment_angle_r = %g\n", segment_angle_r); + + /* We don't have enough resolution to stroke this arc */ + if ( sweep_angle_r < segment_angle_r ) { + if ( include_first ) + pointArrayAddPoint(pa, p1); + pointArrayAddPoint(pa, p3); + return MS_SUCCESS; + } + + /* How many edges to generate (we add the final edge + * by sticking on the last point */ + num_edges = floor(sweep_angle_r / fabs(segment_angle_r)); + + //printf("num_edges = %d\n", num_edges); + + /* Go backwards if we are stroking clockwise */ + if ( side == FP_LEFT ) + segment_angle_r *= -1; + + //printf("segment_angle_r = %g\n", segment_angle_r); + + /* Start with the first point? */ + if( include_first ) { + current_angle_r = a1; + } + else { + current_angle_r = a1 + segment_angle_r; + num_edges--; + } + + //printf("center = %g %g\n", center.x, center.y); + //printf("radius = %g\n", radius); + + for( i = 0; i < num_edges; i++ ) { + if (segment_angle_r > 0.0 && current_angle_r > M_PI) + current_angle_r -= 2*M_PI; + if (segment_angle_r < 0.0 && current_angle_r < -1*M_PI) + current_angle_r -= 2*M_PI; + p.x = center.x + radius*cos(current_angle_r); + p.y = center.y + radius*sin(current_angle_r); + //printf(" point[%d] %g %g\n", i, p.x, p.y); + pointArrayAddPoint(pa, &p); + //printf(" current_angle_r[%d] %g\n",i,current_angle_r); + current_angle_r += segment_angle_r; + } + + /* Add the last point */ + pointArrayAddPoint(pa, p3); + return MS_SUCCESS; +} + +int +arcStrokeCircularString(wkbObj *w, double segment_angle, lineObj *line) +{ + pointObj p1, p2, p3; + int npoints, nedges; + int edge = 0; + pointArrayObj *pa; + + if ( ! w || ! pa ) return MS_FAILURE; + + npoints = wkbReadInt(w); + nedges = npoints / 2; + + /* Make a large guess at how much space we'll need */ + pa = pointArrayNew(nedges * 180 / segment_angle); + + /* All CircularStrings have an odd number of points */ + if ( npoints < 3 || npoints % 2 != 1 ) + return MS_FAILURE; + + wkbReadPointP(w,&p3); + + /* Fill out the point array with stroked arcs */ + while( edge < nedges ) { + p1 = p3; + wkbReadPointP(w,&p2); + wkbReadPointP(w,&p3); + if ( arcStrokeCircle(&p1, &p2, &p3, segment_angle, edge ? 0 : 1, pa) == MS_FAILURE ) { + pointArrayFree(pa); + return MS_FAILURE; + } + edge++; + } + + /* Copy the point array into the line */ + line->numpoints = pa->npoints; + line->point = malloc(line->numpoints * sizeof(pointObj)); + memcpy(line->point, pa->data, line->numpoints * sizeof(pointObj)); + + /* Clean up */ + pointArrayFree(pa); + + return MS_SUCCESS; +} + + /* if there is any polygon in wkb, return force_polygon */ /* if there is any line in wkb, return force_line */ /* otherwise return force_point */ -static int dont_force(unsigned char *wkb, shapeObj *shape) +static int msPostGISFindBestType(wkbObj *w, shapeObj *shape) { - int offset =0; - int ngeoms; - int type, t; - int best_type; + int wkbtype; - best_type = MS_SHAPE_NULL; /* nothing in it */ + wkbtype = wkbCollectionSubType(w); - memcpy(&ngeoms, &wkb[5], 4); - offset = 9; /* were the first geometry is */ - for(t = 0; t < ngeoms; t++) { - memcpy(&type, &wkb[offset + 1], 4); /* type of this geometry */ + switch ( wkbtype ) { + case WKB_POLYGON: + case WKB_CURVEPOLYGON: + case WKB_MULTIPOLYGON: + shape->type = MS_SHAPE_POLYGON; + break; + case WKB_LINESTRING: + case WKB_CIRCULARSTRING: + case WKB_COMPOUNDCURVE: + case WKB_MULTICURVE: + case WKB_MULTILINESTRING: + shape->type = MS_SHAPE_LINE; + break; + case WKB_POINT: + case WKB_MULTIPOINT: + shape->type = MS_SHAPE_POINT; + break; + default: + return MS_FAILURE; + } - if(type == 3) { - best_type = MS_SHAPE_POLYGON; - } else if(type ==2 && best_type != MS_SHAPE_POLYGON) { - best_type = MS_SHAPE_LINE; - } else if(type == 1 && best_type == MS_SHAPE_NULL) { - best_type = MS_SHAPE_POINT; - } - } - - if(best_type == MS_SHAPE_POINT) { - return force_to_points(wkb, shape); - } - if(best_type == MS_SHAPE_LINE) { - return force_to_lines(wkb, shape); - } - if(best_type == MS_SHAPE_POLYGON) { - return force_to_polygons(wkb, shape); - } - - return MS_FAILURE; /* unknown type */ + return wkbConvGeometryToShape(w, shape); } /* find the bounds of the shape */ @@ -1367,6 +1833,7 @@ char *wkbstr = NULL; unsigned char *wkb = NULL; + wkbObj w; msPostGISLayerInfo *layerinfo = NULL; int result = 0; int wkbstrlen = 0; @@ -1399,25 +1866,32 @@ free(wkb); return MS_FAILURE; } + + /* Initialize our wkbObj */ + w.wkb = (char*)wkb; + w.ptr = w.wkb; switch (layer->type) { case MS_LAYER_POINT: - result = force_to_points(wkb, shape); + shape->type = MS_SHAPE_POINT; + result = wkbConvGeometryToShape(&w, shape); break; case MS_LAYER_LINE: - result = force_to_lines(wkb, shape); + shape->type = MS_SHAPE_LINE; + result = wkbConvGeometryToShape(&w, shape); break; case MS_LAYER_POLYGON: - result = force_to_polygons(wkb, shape); + shape->type = MS_SHAPE_POLYGON; + result = wkbConvGeometryToShape(&w, shape); break; case MS_LAYER_ANNOTATION: case MS_LAYER_QUERY: case MS_LAYER_CHART: - result = dont_force(wkb, shape); + result = msPostGISFindBestType(&w, shape); break; case MS_LAYER_RASTER: @@ -1433,7 +1907,10 @@ break; } - if (shape->type != MS_SHAPE_NULL) { + /* All done with geometry, free this */ + free(wkb); + + if (result != MS_FAILURE) { int t; long uid; char *tmp; @@ -1491,7 +1968,6 @@ free(tmp); } - free(wkb); return MS_SUCCESS; } Index: mappostgis.h =================================================================== --- mappostgis.h (revision 10768) +++ mappostgis.h (working copy) @@ -30,13 +30,58 @@ } msPostGISLayerInfo; + /* +** Utility structure for handling the WKB returned by the database while +** reading. +*/ +typedef struct { + char *wkb; /* Pointer to front of WKT */ + char *ptr; /* Pointer to current write point */ + size_t size; /* Size of allocated space */ +} wkbObj; + +/* +** Utility structure used when building up stroked lines while +** handling curved feature types. +*/ +typedef struct { + pointObj *data; /* Pointer to front of WKT */ + int npoints; + int maxpoints; +} pointArrayObj; + +/* +** All the WKB type numbers from the OGC +*/ +typedef enum { + WKB_POINT=1, + WKB_LINESTRING=2, + WKB_POLYGON=3, + WKB_MULTIPOINT=4, + WKB_MULTILINESTRING=5, + WKB_MULTIPOLYGON=6, + WKB_GEOMETRYCOLLECTION=7, + WKB_CIRCULARSTRING=8, + WKB_COMPOUNDCURVE=9, + WKB_CURVEPOLYGON=13, + WKB_MULTICURVE=14, + WKB_MULTISURFACE=15 +} wkb_typenum; + + +/* ** Prototypes */ void msPostGISFreeLayerInfo(layerObj *layer); msPostGISLayerInfo *msPostGISCreateLayerInfo(void); char *msPostGISBuildSQL(layerObj *layer, rectObj *rect, long *uid); int msPostGISParseData(layerObj *layer); +int arcStrokeCircularString(wkbObj *w, double segment_angle, lineObj *line); +int wkbConvGeometryToShape(wkbObj *w, shapeObj *shape); +pointArrayObj* pointArrayNew(int maxpoints); +void pointArrayFree(pointArrayObj *d); +static void wkbSkipGeometry(wkbObj *w); #endif /* USE_POSTGIS */