static LWGEOM*
parse_geojson_linestring(json_object *geojson, int *hasz, int root_srid)
{
	LWGEOM *geom;
	POINTARRAY *pa;
	json_object* points = NULL;
	int i = 0;

	LWDEBUG(2, "parse_geojson_linestring called.");

	points = findMemberByName( geojson, "coordinates" );
	if ( ! points ) 
	{
		geojson_lwerror("Unable to find 'coordinates' in GeoJSON string", 4);
	return NULL;
	}

	pa = ptarray_construct_empty(1, 0, 1);

	if( json_type_array == json_object_get_type( points ) )
	{
		const int nPoints = json_object_array_length( points );
		for(i = 0; i < nPoints; ++i)
		{
			json_object* coords = NULL;
			coords = json_object_array_get_idx( points, i );
			parse_geojson_coord(coords, hasz, pa);
		}
	}

	geom = (LWGEOM *) lwline_construct(root_srid, NULL, pa);

	LWDEBUG(2, "parse_geojson_linestring finished.");
	return geom;
}

int lwcompound_add_lwgeom(LWCOMPOUND *comp, LWGEOM *geom)
{
	LWCOLLECTION *col = (LWCOLLECTION*)comp;
	
	/* Empty things can't continuously join up with other things */
	if ( lwgeom_is_empty(geom) )
	{
		LWDEBUG(4, "Got an empty component for a compound curve!");
		return LW_FAILURE;
	}
	
	if( col->ngeoms > 0 )
	{
		POINT4D last, first;
		/* First point of the component we are adding */
		LWLINE *newline = (LWLINE*)geom;
		/* Last point of the previous component */
		LWLINE *prevline = (LWLINE*)(col->geoms[col->ngeoms-1]);

		getPoint4d_p(newline->points, 0, &first);
		getPoint4d_p(prevline->points, prevline->points->npoints-1, &last);
		
		if ( !(FP_EQUALS(first.x,last.x) && FP_EQUALS(first.y,last.y)) )
		{
			LWDEBUG(4, "Components don't join up end-to-end!");
			LWDEBUGF(4, "first pt (%g %g %g %g) last pt (%g %g %g %g)", first.x, first.y, first.z, first.m, last.x, last.y, last.z, last.m);			
			return LW_FAILURE;
		}
	}
	
	col = lwcollection_add_lwgeom(col, geom);
	return LW_SUCCESS;
}

char *
lwcollection_summary(LWCOLLECTION *col, int offset)
{
	size_t size = 128;
	char *result;
	char *tmp;
	int i;
	char *pad="";

	LWDEBUG(2, "lwcollection_summary called");

	result = (char *)lwalloc(size);

	sprintf(result, "%*.s%s[%s] with %d elements\n",
	        offset, pad, lwgeom_typename(TYPE_GETTYPE(col->type)),
	        lwgeom_typeflags(col->type),
	        col->ngeoms);

	for (i=0; i<col->ngeoms; i++)
	{
		tmp = lwgeom_summary(col->geoms[i], offset+2);
		size += strlen(tmp)+1;
		result = lwrealloc(result, size);

		LWDEBUGF(4, "Reallocated %d bytes for result", size);

		strcat(result, tmp);
		lwfree(tmp);
	}

	LWDEBUG(3, "lwcollection_summary returning");

	return result;
}

char *
lwpoly_summary(LWPOLY *poly, int offset)
{
	char tmp[256];
	size_t size = 64*(poly->nrings+1)+128;
	char *result;
	int i;
	char *pad="";

	LWDEBUG(2, "lwpoly_summary called");

	result = lwalloc(size);

	sprintf(result, "%*.s%s[%s] with %i rings\n",
	        offset, pad, lwgeom_typename(TYPE_GETTYPE(poly->type)),
	        lwgeom_typeflags(poly->type),
	        poly->nrings);

	for (i=0; i<poly->nrings; i++)
	{
		sprintf(tmp,"%s   ring %i has %i points\n",
		        pad, i, poly->rings[i]->npoints);
		strcat(result,tmp);
	}

	LWDEBUG(3, "lwpoly_summary returning");

	return result;
}

/**
Function initializing closestpoint calculations.
*/
LWGEOM *
lw_dist2d_distancepoint(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode)
{
	double x,y;
	DISTPTS thedl;
	double initdistance = MAXFLOAT;
	LWGEOM *result;

	thedl.mode = mode;
	thedl.distance= initdistance;
	thedl.tolerance = 0;

	LWDEBUG(2, "lw_dist2d_distancepoint is called");

	if (!lw_dist2d_comp( lw1,lw2,&thedl))
	{
		/*should never get here. all cases ought to be error handled earlier*/
		lwerror("Some unspecified error.");
		result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
	}
	if (thedl.distance == initdistance)
	{
		LWDEBUG(3, "didn't find geometries to measure between, returning null");
		result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
	}
	else
	{
		x=thedl.p1.x;
		y=thedl.p1.y;
		result = (LWGEOM *)lwpoint_make2d(srid, x, y);
	}
	return result;
}

/**

Checking if the point projected on the plane of the polygon actually is inside that polygon. 
If so the mindistance is between that projected point and our original point.
If not we check from original point to the bounadary.
If the projected point is inside a hole of the polygon we check the distance to the boudary of that hole.
*/
int
lw_dist3d_pt_poly(POINT3DZ *p, LWPOLY *poly, PLANE3D *plane,POINT3DZ *projp, DISTPTS3D *dl)
{	
	int i;
	
	LWDEBUG(2, "lw_dist3d_point_poly called");

	
	if(pt_in_ring_3d(projp, poly->rings[0], plane))
	{
		for (i=1; i<poly->nrings; i++)
		{
			/* Inside a hole. Distance = pt -> ring */
			if ( pt_in_ring_3d(projp, poly->rings[i], plane ))
			{
				LWDEBUG(3, " inside an hole");
				return lw_dist3d_pt_ptarray(p, poly->rings[i], dl);
			}
		}		
		
		return lw_dist3d_pt_pt(p,projp,dl);/* If the projected point is inside the polygon the shortest distance is between that point and the inputed point*/
	}
	else
	{
		return lw_dist3d_pt_ptarray(p, poly->rings[0], dl); /*If the projected point is outside the polygon we search for the closest distance against the boundarry instead*/
	}	
	
	return LW_TRUE;
	
}

/**
* POLYGON
*/
static LWPOLY* lwpoly_from_twkb_state(twkb_parse_state *s)
{
	uint32_t nrings;
	int i;
	LWPOLY *poly;

	LWDEBUG(2,"Entering lwpoly_from_twkb_state");

	if ( s->is_empty )
		return lwpoly_construct_empty(SRID_UNKNOWN, s->has_z, s->has_m);

	/* Read number of rings */
	nrings = twkb_parse_state_uvarint(s);

	/* Start w/ empty polygon */
	poly = lwpoly_construct_empty(SRID_UNKNOWN, s->has_z, s->has_m);

	LWDEBUGF(4,"Polygon has %d rings", nrings);

	/* Empty polygon? */
	if( nrings == 0 )
		return poly;

	for( i = 0; i < nrings; i++ )
	{
		/* Ret number of points */
		uint32_t npoints = twkb_parse_state_uvarint(s);
		POINTARRAY *pa = ptarray_from_twkb_state(s, npoints);

		/* Skip empty rings */
		if( pa == NULL )
			continue;

		/* Force first and last points to be the same. */
		if( ! ptarray_is_closed_2d(pa) )
		{
			POINT4D pt;
			getPoint4d_p(pa, 0, &pt);
			ptarray_append_point(pa, &pt, LW_FALSE);
		}

		/* Check for at least four points. */
		if( s->check & LW_PARSER_CHECK_MINPOINTS && pa->npoints < 4 )
		{
			LWDEBUGF(2, "%s must have at least four points in each ring", lwtype_name(s->lwtype));
			lwerror("%s must have at least four points in each ring", lwtype_name(s->lwtype));
			return NULL;
		}

		/* Add ring to polygon */
		if ( lwpoly_add_ring(poly, pa) == LW_FAILURE )
		{
			LWDEBUG(2, "Unable to add ring to polygon");
			lwerror("Unable to add ring to polygon");
		}

	}
	return poly;
}

/*
 * Given the LWPOINT serialized form (or a pointer into a muli* one)
 * construct a proper LWPOINT.
 * serialized_form should point to the 8bit type format (with type = 1)
 * See serialized form doc
 */
LWPOINT *
lwpoint_deserialize(uchar *serialized_form)
{
	uchar type;
	int geom_type;
	LWPOINT *result;
	uchar *loc = NULL;
	POINTARRAY *pa;

	LWDEBUG(2, "lwpoint_deserialize called");

	result = (LWPOINT*) lwalloc(sizeof(LWPOINT)) ;

	type = serialized_form[0];
	geom_type = lwgeom_getType(type);

	if ( geom_type != POINTTYPE)
	{
		lwerror("lwpoint_deserialize: attempt to deserialize a point which is really a %s", lwgeom_typename(geom_type));
		return NULL;
	}
	result->type = type;

	loc = serialized_form+1;

	if (lwgeom_hasBBOX(type))
	{
		LWDEBUG(3, "lwpoint_deserialize: input has bbox");

		result->bbox = lwalloc(sizeof(BOX2DFLOAT4));
		memcpy(result->bbox, loc, sizeof(BOX2DFLOAT4));
		loc += sizeof(BOX2DFLOAT4);
	}
	else
	{
		result->bbox = NULL;
	}

	if ( lwgeom_hasSRID(type))
	{
		LWDEBUG(3, "lwpoint_deserialize: input has SRID");

		result->SRID = lw_get_int32(loc);
		loc += 4; /* type + SRID */
	}
	else
	{
		result->SRID = -1;
	}

	/* we've read the type (1 byte) and SRID (4 bytes, if present) */

	pa = pointArray_construct(loc, TYPE_HASZ(type), TYPE_HASM(type), 1);

	result->point = pa;

	return result;
}

int lwcurvepoly_add_ring(LWCURVEPOLY *poly, LWGEOM *ring)
{
	int i;
	
	/* Can't do anything with NULLs */
	if( ! poly || ! ring ) 
	{
		LWDEBUG(4,"NULL inputs!!! quitting");
		return LW_FAILURE;
	}

	/* Check that we're not working with garbage */
	if ( poly->rings == NULL && (poly->nrings || poly->maxrings) )
	{
		LWDEBUG(4,"mismatched nrings/maxrings");
		lwerror("Curvepolygon is in inconsistent state. Null memory but non-zero collection counts.");
	}

	/* Check that we're adding an allowed ring type */
	if ( ! ( ring->type == LINETYPE || ring->type == CIRCSTRINGTYPE || ring->type == COMPOUNDTYPE ) )
	{
		LWDEBUGF(4,"got incorrect ring type: %s",lwtype_name(ring->type));
		return LW_FAILURE;
	}

		
	/* In case this is a truly empty, make some initial space  */
	if ( poly->rings == NULL )
	{
		poly->maxrings = 2;
		poly->nrings = 0;
		poly->rings = (LWGEOM **)lwalloc(poly->maxrings * sizeof(LWGEOM*));
	}

	/* Allocate more space if we need it */
	if ( poly->nrings == poly->maxrings )
	{
		poly->maxrings *= 2;
		poly->rings = (LWGEOM **)lwrealloc(poly->rings, sizeof(LWGEOM*) * poly->maxrings);
	}

	/* Make sure we don't already have a reference to this geom */
	for ( i = 0; i < poly->nrings; i++ )
	{
		if ( poly->rings[i] == ring )
		{
			LWDEBUGF(4, "Found duplicate geometry in collection %p == %p", poly->rings[i], ring);
			return LW_SUCCESS;
		}
	}

	/* Add the ring and increment the ring count */
	poly->rings[poly->nrings] = (LWGEOM*)ring;
	poly->nrings++;
	return LW_SUCCESS;	
}

int gbox_from_gserialized(const GSERIALIZED *g, GBOX *gbox)
{

	/* Null input! */
	if ( ! g ) return G_FAILURE;

	/* Initialize the flags on the box */
	gbox->flags = g->flags;

	if ( FLAGS_GET_BBOX(g->flags) )
	{
		int i = 0;
		float *fbox = (float*)(g->data);
		gbox->xmin = fbox[i];
		i++;
		gbox->xmax = fbox[i];
		i++;
		gbox->ymin = fbox[i];
		i++;
		gbox->ymax = fbox[i];
		i++;
		if ( FLAGS_GET_GEODETIC(g->flags) )
		{
			gbox->zmin = fbox[i];
			i++;
			gbox->zmax = fbox[i];
			i++;
			return G_SUCCESS;
		}
		if ( FLAGS_GET_Z(g->flags) )
		{
			gbox->zmin = fbox[i];
			i++;
			gbox->zmax = fbox[i];
			i++;
		}
		if ( FLAGS_GET_M(g->flags) )
		{
			gbox->mmin = fbox[i];
			i++;
			gbox->mmax = fbox[i];
			i++;
		}
		return G_SUCCESS;
	}

	LWDEBUG(4, "calculating new box from scratch");
	if ( gserialized_calculate_gbox_geocentric_p(g, gbox) == G_FAILURE )
	{
		LWDEBUG(4, "calculated null bbox, returning failure");
		return G_FAILURE;
	}
	return G_SUCCESS;
}

/**
 * @brief Add a point in a pointarray.
 *
 * @param pa the source POINTARRAY
 * @param p the point to add
 * @param pdims number of ordinates in p (2..4)
 * @param where to insert the point. 0 prepends, pa->npoints appends
 *
 * @returns a newly constructed POINTARRAY using a newly allocated buffer
 *          for the actual points, or NULL on error.
 */
POINTARRAY *
ptarray_addPoint(const POINTARRAY *pa, uint8_t *p, size_t pdims, uint32_t where)
{
	POINTARRAY *ret;
	POINT4D pbuf;
	size_t ptsize = ptarray_point_size(pa);

	LWDEBUGF(3, "pa %x p %x size %d where %d",
	         pa, p, pdims, where);

	if ( pdims < 2 || pdims > 4 )
	{
		lwerror("ptarray_addPoint: point dimension out of range (%d)",
		        pdims);
		return NULL;
	}

	if ( where > pa->npoints )
	{
		lwerror("ptarray_addPoint: offset out of range (%d)",
		        where);
		return NULL;
	}

	LWDEBUG(3, "called with a %dD point");

	pbuf.x = pbuf.y = pbuf.z = pbuf.m = 0.0;
	memcpy((uint8_t *)&pbuf, p, pdims*sizeof(double));

	LWDEBUG(3, "initialized point buffer");

	ret = ptarray_construct(FLAGS_GET_Z(pa->flags),
	                        FLAGS_GET_M(pa->flags), pa->npoints+1);

	if ( where == -1 ) where = pa->npoints;

	if ( where )
	{
		memcpy(getPoint_internal(ret, 0), getPoint_internal(pa, 0), ptsize*where);
	}

	memcpy(getPoint_internal(ret, where), (uint8_t *)&pbuf, ptsize);

	if ( where+1 != ret->npoints )
	{
		memcpy(getPoint_internal(ret, where+1),
		       getPoint_internal(pa, where),
		       ptsize*(pa->npoints-where));
	}

	return ret;
}

/**
* Byte
* Read a byte and advance the parse state forward.
*/
static char byte_from_wkb_state(wkb_parse_state *s)
{
	char char_value = 0;
	LWDEBUG(4, "Entered function");

	wkb_parse_state_check(s, WKB_BYTE_SIZE);
	LWDEBUG(4, "Passed state check");

	char_value = s->pos[0];
	LWDEBUGF(4, "Read byte value: %x", char_value);
	s->pos += WKB_BYTE_SIZE;

	return char_value;
}

/**
 * 1. see if pt in outer boundary. if no, then treat the outer ring like a line
 * 2. if in the boundary, test to see if its in a hole.
 *    if so, then return dist to hole, else return 0 (point in polygon)
 */
int
lw_dist2d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl)
{
	POINT2D p;
	int i;

	LWDEBUG(2, "lw_dist2d_point_poly called");

	getPoint2d_p(point->point, 0, &p);

	if (dl->mode == DIST_MAX)
	{
		LWDEBUG(3, "looking for maxdistance");
		return lw_dist2d_pt_ptarray(&p, poly->rings[0], dl);
	}
	/* Return distance to outer ring if not inside it */
	if ( ! pt_in_ring_2d(&p, poly->rings[0]) )
	{
		LWDEBUG(3, "first point not inside outer-ring");
		return lw_dist2d_pt_ptarray(&p, poly->rings[0], dl);
	}

	/*
	 * Inside the outer ring.
	 * Scan though each of the inner rings looking to
	 * see if its inside.  If not, distance==0.
	 * Otherwise, distance = pt to ring distance
	 */
	for (i=1; i<poly->nrings; i++)
	{
		/* Inside a hole. Distance = pt -> ring */
		if ( pt_in_ring_2d(&p, poly->rings[i]) )
		{
			LWDEBUG(3, " inside an hole");
			return lw_dist2d_pt_ptarray(&p, poly->rings[i], dl);
		}
	}

	LWDEBUG(3, " inside the polygon");
	if (dl->mode == DIST_MIN)
	{
		dl->distance=0.0;
		dl->p1.x=p.x;
		dl->p1.y=p.y;
		dl->p2.x=p.x;
		dl->p2.y=p.y;
	}
	return LW_TRUE; /* Is inside the polygon */
}

static void
ptarray_dp_findsplit(POINTARRAY *pts, int p1, int p2, int *split, double *dist)
{
	int k;
	const POINT2D *pk, *pa, *pb;
	double tmp, d;

	LWDEBUG(4, "function called");

	*split = p1;
	d = -1;

	if (p1 + 1 < p2)
	{

		pa = getPoint2d_cp(pts, p1);
		pb = getPoint2d_cp(pts, p2);

		LWDEBUGF(4, "P%d(%f,%f) to P%d(%f,%f)",
		         p1, pa->x, pa->y, p2, pb->x, pb->y);

		for (k=p1+1; k<p2; k++)
		{
			pk = getPoint2d_cp(pts, k);

			LWDEBUGF(4, "P%d(%f,%f)", k, pk->x, pk->y);

			/* distance computation */
			tmp = distance2d_sqr_pt_seg(pk, pa, pb);

			if (tmp > d)
			{
				d = tmp;	/* record the maximum */
				*split = k;

				LWDEBUGF(4, "P%d is farthest (%g)", k, d);
			}
		}
		*dist = d;

	} /* length---should be redone if can == 0 */
	else
	{
		LWDEBUG(3, "segment too short, no split/no dist");
		*dist = -1;
	}

}

char compound_is_closed(LWCOMPOUND *compound)
{
	POINT3DZ sp, ep;
	LWGEOM *tmp;

	LWDEBUG(2, "compound_is_closed called.");

	tmp = compound->geoms[0];
	if (lwgeom_getType(tmp->type) == LINETYPE)
	{
		getPoint3dz_p(((LWLINE *)tmp)->points, 0, &sp);
	}
	else
	{
		getPoint3dz_p(((LWCIRCSTRING *)tmp)->points, 0, &sp);
	}

	tmp = compound->geoms[compound->ngeoms - 1];
	if (lwgeom_getType(tmp->type) == LINETYPE)
	{
		getPoint3dz_p(((LWLINE *)tmp)->points, ((LWLINE *)tmp)->points->npoints - 1, &ep);
	}
	else
	{
		getPoint3dz_p(((LWCIRCSTRING *)tmp)->points, ((LWCIRCSTRING *)tmp)->points->npoints - 1, &ep);
	}

	if (sp.x != ep.x) return 0;
	if (sp.y != ep.y) return 0;
	if (TYPE_HASZ(compound->type))
	{
		if (sp.z != ep.z) return 0;
	}
	return 1;
}

/*
 * Find bounding box (standard one) 
 *   zmin=zmax=NO_Z_VALUE if 2d 
 */
BOX3D *
lwpoint_compute_box3d(LWPOINT *point)
{
	LWDEBUGF(2, "lwpoint_compute_box3d called with point %p", point);

	if (point == NULL)
	{
		LWDEBUG(3, "lwpoint_compute_box3d returning NULL");

		return NULL;
	}

	LWDEBUG(3, "lwpoint_compute_box3d returning ptarray_compute_box3d return");

	return ptarray_compute_box3d(point->point);
}

/**
Function initializing 3d max distance calculation
*/
double
lwgeom_maxdistance3d(const LWGEOM *lw1, const LWGEOM *lw2)
{
	LWDEBUG(2, "lwgeom_maxdistance3d is called");

	return lwgeom_maxdistance3d_tolerance( lw1, lw2, 0.0 );
}

LWGEOM *
lwmpolygon_desegmentize(LWMPOLY *mpoly)
{
	LWGEOM **geoms;
	int i, hascurve = 0;

	LWDEBUG(2, "lwmpoly_desegmentize called.");

	geoms = lwalloc(sizeof(LWGEOM *)*mpoly->ngeoms);
	for (i=0; i<mpoly->ngeoms; i++)
	{
		geoms[i] = lwpolygon_desegmentize((LWPOLY *)mpoly->geoms[i]);
		if (geoms[i]->type == CURVEPOLYTYPE)
		{
			hascurve = 1;
		}
	}
	if (hascurve == 0)
	{
		for (i=0; i<mpoly->ngeoms; i++)
		{
			lwfree(geoms[i]);
		}
		return lwgeom_clone((LWGEOM *)mpoly);
	}
	return (LWGEOM *)lwcollection_construct(MULTISURFACETYPE, mpoly->srid, NULL, mpoly->ngeoms, geoms);
}

/**

polygon to polygon calculation
*/
int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl)
{		
	PLANE3D plane;		
	LWDEBUG(2, "lw_dist3d_poly_poly is called");
	if (dl->mode == DIST_MAX)
	{
		return lw_dist3d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl);
	}
	
	if(!define_plane(poly2->rings[0], &plane))
		return LW_FALSE;
	
	/*What we do here is to compare the bondary of one polygon with the other polygon 
	and then take the second boudary comparing with the first polygon*/
	dl->twisted=1;
	if(!lw_dist3d_ptarray_poly(poly1->rings[0], poly2,&plane, dl))
		return LW_FALSE;
	if(dl->distance==0.0) /*Just check if the answer already is given*/
		return LW_TRUE;
	
	if(!define_plane(poly1->rings[0], &plane))
		return LW_FALSE;
	dl->twisted=-1; /*because we swithc the order of geometries we swithch "twisted" to -1 which will give the right order of points in shortest line.*/
	return lw_dist3d_ptarray_poly(poly2->rings[0], poly1,&plane, dl);
}

LWGEOM *
lwmline_desegmentize(LWMLINE *mline)
{
	LWGEOM **geoms;
	int i, hascurve = 0;

	LWDEBUG(2, "lwmline_desegmentize called.");

	geoms = lwalloc(sizeof(LWGEOM *)*mline->ngeoms);
	for (i=0; i<mline->ngeoms; i++)
	{
		geoms[i] = lwline_desegmentize((LWLINE *)mline->geoms[i]);
		if (geoms[i]->type == CIRCSTRINGTYPE || geoms[i]->type == COMPOUNDTYPE)
		{
			hascurve = 1;
		}
	}
	if (hascurve == 0)
	{
		for (i=0; i<mline->ngeoms; i++)
		{
			lwfree(geoms[i]);
		}
		return lwgeom_clone((LWGEOM *)mline);
	}
	return (LWGEOM *)lwcollection_construct(MULTICURVETYPE, mline->srid, NULL, mline->ngeoms, geoms);
}