static LWPOLY* lwpoly_from_gserialized_buffer(uint8_t *data_ptr, uint8_t g_flags, size_t *g_size)
{
	uint8_t *start_ptr = data_ptr;
	LWPOLY *poly;
	uint8_t *ordinate_ptr;
	uint32_t nrings = 0;
	int i = 0;

	assert(data_ptr);

	poly = (LWPOLY*)lwalloc(sizeof(LWPOLY));
	poly->srid = SRID_UNKNOWN; /* Default */
	poly->bbox = NULL;
	poly->type = POLYGONTYPE;
	poly->flags = g_flags;

	data_ptr += 4; /* Skip past the polygontype. */
	nrings = lw_get_uint32_t(data_ptr); /* Zero => empty geometry */
	poly->nrings = nrings;
	LWDEBUGF(4, "nrings = %d", nrings);
	data_ptr += 4; /* Skip past the nrings. */

	ordinate_ptr = data_ptr; /* Start the ordinate pointer. */
	if ( nrings > 0)
	{
		poly->rings = (POINTARRAY**)lwalloc( sizeof(POINTARRAY*) * nrings );
		ordinate_ptr += nrings * 4; /* Move past all the npoints values. */
		if ( nrings % 2 ) /* If there is padding, move past that too. */
			ordinate_ptr += 4;
	}
	else /* Empty polygon */
	{
		poly->rings = NULL;
	}

	for ( i = 0; i < nrings; i++ )
	{
		uint32_t npoints = 0;

		/* Read in the number of points. */
		npoints = lw_get_uint32_t(data_ptr);
		data_ptr += 4;

		/* Make a point array for the ring, and move the ordinate pointer past the ring ordinates. */
		poly->rings[i] = ptarray_construct_reference_data(FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), npoints, ordinate_ptr);
		
		ordinate_ptr += sizeof(double) * FLAGS_NDIMS(g_flags) * npoints;
	}

	if ( g_size )
		*g_size = ordinate_ptr - start_ptr;

	return poly;
}

/*
 * Construct a triangle from a LWLINE being
 * the shell
 * Pointarray from intput geom are cloned.
 * Input line must have 4 points, and be closed.
 */
LWTRIANGLE *
lwtriangle_from_lwline(const LWLINE *shell)
{
	LWTRIANGLE *ret;
	POINTARRAY *pa;

	if ( shell->points->npoints != 4 )
		lwerror("lwtriangle_from_lwline: shell must have exactly 4 points");

	if (   (!FLAGS_GET_Z(shell->flags) && !ptarray_is_closed_2d(shell->points)) ||
	        (FLAGS_GET_Z(shell->flags) && !ptarray_is_closed_3d(shell->points)) )
		lwerror("lwtriangle_from_lwline: shell must be closed");

	pa = ptarray_clone_deep(shell->points);
	ret = lwtriangle_construct(shell->srid, NULL, pa);

	if (lwtriangle_is_repeated_points(ret))
		lwerror("lwtriangle_from_lwline: some points are repeated in triangle");

	return ret;
}

static size_t gserialized_from_lwline(const LWLINE *line, uint8_t *buf)
{
	uint8_t *loc;
	int ptsize;
	size_t size;
	int type = LINETYPE;

	assert(line);
	assert(buf);

	LWDEBUGF(2, "lwline_to_gserialized(%p, %p) called", line, buf);

	if ( FLAGS_GET_Z(line->flags) != FLAGS_GET_Z(line->points->flags) )
		lwerror("Dimensions mismatch in lwline");

	ptsize = ptarray_point_size(line->points);

	loc = buf;

	/* Write in the type. */
	memcpy(loc, &type, sizeof(uint32_t));
	loc += sizeof(uint32_t);

	/* Write in the npoints. */
	memcpy(loc, &(line->points->npoints), sizeof(uint32_t));
	loc += sizeof(uint32_t);

	LWDEBUGF(3, "lwline_to_gserialized added npoints (%d)", line->points->npoints);

	/* Copy in the ordinates. */
	if ( line->points->npoints > 0 )
	{
		size = line->points->npoints * ptsize;
		memcpy(loc, getPoint_internal(line->points, 0), size);
		loc += size;
	}
	LWDEBUGF(3, "lwline_to_gserialized copied serialized_pointlist (%d bytes)", ptsize * line->points->npoints);

	return (size_t)(loc - buf);
}

/**
 * @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;
}

static size_t
asgeojson_line_size(const LWLINE *line, char *srs, GBOX *bbox, int precision)
{
	int size;

	size = sizeof("{'type':'LineString',");
	if (srs) size += asgeojson_srs_size(srs);
	if (bbox) size += asgeojson_bbox_size(FLAGS_GET_Z(line->flags), precision);
	size += sizeof("'coordinates':[]}");
	size += pointArray_geojson_size(line->points, precision);

	return size;
}

LWLINE *
lwcircstring_segmentize(const LWCIRCSTRING *icurve, uint32_t perQuad)
{
	LWLINE *oline;
	POINTARRAY *ptarray;
	POINTARRAY *tmp;
	uint32_t i, j;
	POINT4D p1, p2, p3, p4;

	LWDEBUGF(2, "lwcircstring_segmentize called., dim = %d", icurve->points->flags);

	ptarray = ptarray_construct_empty(FLAGS_GET_Z(icurve->points->flags), FLAGS_GET_M(icurve->points->flags), 64);

	for (i = 2; i < icurve->points->npoints; i+=2)
	{
		LWDEBUGF(3, "lwcircstring_segmentize: arc ending at point %d", i);

		getPoint4d_p(icurve->points, i - 2, &p1);
		getPoint4d_p(icurve->points, i - 1, &p2);
		getPoint4d_p(icurve->points, i, &p3);
		tmp = lwcircle_segmentize(&p1, &p2, &p3, perQuad);

		if (tmp)
		{
			LWDEBUGF(3, "lwcircstring_segmentize: generated %d points", tmp->npoints);

			for (j = 0; j < tmp->npoints; j++)
			{
				getPoint4d_p(tmp, j, &p4);
				ptarray_append_point(ptarray, &p4, LW_TRUE);
			}
			ptarray_free(tmp);
		}
		else
		{
			LWDEBUG(3, "lwcircstring_segmentize: points are colinear, returning curve points as line");

			for (j = i - 1 ; j <= i ; j++)
			{
				getPoint4d_p(icurve->points, j, &p4);
				ptarray_append_point(ptarray, &p4, LW_TRUE);
			}
		}

	}
	getPoint4d_p(icurve->points, icurve->points->npoints-1, &p1);
	ptarray_append_point(ptarray, &p1, LW_TRUE);
		
	oline = lwline_construct(icurve->srid, NULL, ptarray);
	return oline;
}

/*
 * @param icurve input curve
 * @param tol tolerance, semantic driven by tolerance_type
 * @param tolerance_type see LW_LINEARIZE_TOLERANCE_TYPE
 * @param flags see flags in lwarc_linearize
 *
 * @return a newly allocated LWLINE
 */
static LWLINE *
lwcircstring_linearize(const LWCIRCSTRING *icurve, double tol,
                        LW_LINEARIZE_TOLERANCE_TYPE tolerance_type,
                        int flags)
{
	LWLINE *oline;
	POINTARRAY *ptarray;
	uint32_t i, j;
	POINT4D p1, p2, p3, p4;
	int ret;

	LWDEBUGF(2, "lwcircstring_linearize called., dim = %d", icurve->points->flags);

	ptarray = ptarray_construct_empty(FLAGS_GET_Z(icurve->points->flags), FLAGS_GET_M(icurve->points->flags), 64);

	for (i = 2; i < icurve->points->npoints; i+=2)
	{
		LWDEBUGF(3, "lwcircstring_linearize: arc ending at point %d", i);

		getPoint4d_p(icurve->points, i - 2, &p1);
		getPoint4d_p(icurve->points, i - 1, &p2);
		getPoint4d_p(icurve->points, i, &p3);

		ret = lwarc_linearize(ptarray, &p1, &p2, &p3, tol, tolerance_type, flags);
		if ( ret > 0 )
		{
			LWDEBUGF(3, "lwcircstring_linearize: generated %d points", ptarray->npoints);
		}
		else if ( ret == 0 )
		{
			LWDEBUG(3, "lwcircstring_linearize: points are colinear, returning curve points as line");

			for (j = i - 2 ; j < i ; j++)
			{
				getPoint4d_p(icurve->points, j, &p4);
				ptarray_append_point(ptarray, &p4, LW_TRUE);
			}
		}
		else
		{
			/* An error occurred, lwerror should have been called by now */
			ptarray_free(ptarray);
			return NULL;
		}
	}
	getPoint4d_p(icurve->points, icurve->points->npoints-1, &p1);
	ptarray_append_point(ptarray, &p1, LW_TRUE);

	oline = lwline_construct(icurve->srid, NULL, ptarray);
	return oline;
}

/*
 * Stick an array of points to the given gridspec.
 * Return "gridded" points in *outpts and their number in *outptsn.
 *
 * Two consecutive points falling on the same grid cell are collapsed
 * into one single point.
 *
 */
POINTARRAY *
ptarray_grid(POINTARRAY *pa, gridspec *grid)
{
	POINT4D pbuf;
	int ipn, opn; /* point numbers (input/output) */
	POINTARRAY *dpa;

	POSTGIS_DEBUGF(2, "ptarray_grid called on %p", pa);

	dpa = ptarray_construct_empty(FLAGS_GET_Z(pa->flags),FLAGS_GET_M(pa->flags), pa->npoints);

	for (ipn=0, opn=0; ipn<pa->npoints; ++ipn)
	{

		getPoint4d_p(pa, ipn, &pbuf);

		if ( grid->xsize )
			pbuf.x = rint((pbuf.x - grid->ipx)/grid->xsize) *
			         grid->xsize + grid->ipx;

		if ( grid->ysize )
			pbuf.y = rint((pbuf.y - grid->ipy)/grid->ysize) *
			         grid->ysize + grid->ipy;

		if ( FLAGS_GET_Z(pa->flags) && grid->zsize )
			pbuf.z = rint((pbuf.z - grid->ipz)/grid->zsize) *
			         grid->zsize + grid->ipz;

		if ( FLAGS_GET_M(pa->flags) && grid->msize )
			pbuf.m = rint((pbuf.m - grid->ipm)/grid->msize) *
			         grid->msize + grid->ipm;

		ptarray_append_point(dpa, &pbuf, LW_FALSE);

	}

	return dpa;
}

/*
 * Stick an array of points to the given gridspec.
 * Return "gridded" points in *outpts and their number in *outptsn.
 *
 * Two consecutive points falling on the same grid cell are collapsed
 * into one single point.
 *
 */
POINTARRAY *
ptarray_grid(const POINTARRAY *pa, const gridspec *grid)
{
	POINT4D pt;
	int ipn; /* input point numbers */
	POINTARRAY *dpa;

	LWDEBUGF(2, "ptarray_grid called on %p", pa);

	dpa = ptarray_construct_empty(FLAGS_GET_Z(pa->flags),FLAGS_GET_M(pa->flags), pa->npoints);

	for (ipn=0; ipn<pa->npoints; ++ipn)
	{

		getPoint4d_p(pa, ipn, &pt);

		if ( grid->xsize )
			pt.x = rint((pt.x - grid->ipx)/grid->xsize) *
			         grid->xsize + grid->ipx;

		if ( grid->ysize )
			pt.y = rint((pt.y - grid->ipy)/grid->ysize) *
			         grid->ysize + grid->ipy;

		if ( FLAGS_GET_Z(pa->flags) && grid->zsize )
			pt.z = rint((pt.z - grid->ipz)/grid->zsize) *
			         grid->zsize + grid->ipz;

		if ( FLAGS_GET_M(pa->flags) && grid->msize )
			pt.m = rint((pt.m - grid->ipm)/grid->msize) *
			         grid->msize + grid->ipm;

		ptarray_append_point(dpa, &pt, LW_FALSE);

	}

	return dpa;
}

int gbox_overlaps(const GBOX *g1, const GBOX *g2)
{

	/* Make sure our boxes are consistent */
	if ( FLAGS_GET_GEODETIC(g1->flags) != FLAGS_GET_GEODETIC(g2->flags) )
		lwerror("gbox_overlaps: cannot compare geodetic and non-geodetic boxes");

	/* Check X/Y first */
	if ( g1->xmax < g2->xmin || g1->ymax < g2->ymin ||
	     g1->xmin > g2->xmax || g1->ymin > g2->ymax )
		return LW_FALSE;

	/* Deal with the geodetic case special: we only compare the geodetic boxes (x/y/z) */
	/* Never the M dimension */
	if ( FLAGS_GET_GEODETIC(g1->flags) && FLAGS_GET_GEODETIC(g2->flags) )
	{
		if ( g1->zmax < g2->zmin || g1->zmin > g2->zmax )
			return LW_FALSE;
		else
			return LW_TRUE;		
	}
		
	/* If both geodetic or both have Z, check Z */
	if ( FLAGS_GET_Z(g1->flags) && FLAGS_GET_Z(g2->flags) )
	{
		if ( g1->zmax < g2->zmin || g1->zmin > g2->zmax )
			return LW_FALSE;
	}
	
	/* If both have M, check M */
	if ( FLAGS_GET_M(g1->flags) && FLAGS_GET_M(g2->flags) )
	{
		if ( g1->mmax < g2->mmin || g1->mmin > g2->mmax )
			return LW_FALSE;
	}
	
	return LW_TRUE;
}

/**

We calculate the effective area for the first time
*/
void ptarray_calc_areas(EFFECTIVE_AREAS *ea, int avoid_collaps, int set_area, double trshld)
{
	LWDEBUG(2, "Entered  ptarray_calc_areas");
	int i;
	int npoints=ea->inpts->npoints;
	int is3d = FLAGS_GET_Z(ea->inpts->flags);
	double area;
	
	const double *P1;
	const double *P2;
	const double *P3;	
		
	P1 = (double*)getPoint_internal(ea->inpts, 0);
	P2 = (double*)getPoint_internal(ea->inpts, 1);
	
	/*The first and last point shall always have the maximum effective area. We use float max to not make trouble for bbox*/
	ea->initial_arealist[0].area=ea->initial_arealist[npoints-1].area=FLT_MAX;
	ea->res_arealist[0]=ea->res_arealist[npoints-1]=FLT_MAX;
	
	ea->initial_arealist[0].next=1;
	ea->initial_arealist[0].prev=0;
	
	for (i=1;i<(npoints)-1;i++)
	{
		ea->initial_arealist[i].next=i+1;
		ea->initial_arealist[i].prev=i-1;
		P3 = (double*)getPoint_internal(ea->inpts, i+1);

		if(is3d)
			area=triarea3d(P1, P2, P3);
		else
			area=triarea2d(P1, P2, P3);
		
		LWDEBUGF(4,"Write area %lf to point %d on address %p",area,i,&(ea->initial_arealist[i].area));
		ea->initial_arealist[i].area=area;
		P1=P2;
		P2=P3;
		
	}	
		ea->initial_arealist[npoints-1].next=npoints-1;
		ea->initial_arealist[npoints-1].prev=npoints-2;
	
	for (i=1;i<(npoints)-1;i++)
	{
		ea->res_arealist[i]=FLT_MAX;
	}
	
	tune_areas(ea,avoid_collaps,set_area, trshld);
	return ;
}

static size_t
asgeojson_point_size(const LWPOINT *point, char *srs, GBOX *bbox, int precision)
{
	int size;

	size = pointArray_geojson_size(point->point, precision);
	size += sizeof("{'type':'Point',");
	size += sizeof("'coordinates':}");

	if (srs) size += asgeojson_srs_size(srs);
	if (bbox) size += asgeojson_bbox_size(FLAGS_GET_Z(point->flags), precision);

	return size;
}

int gbox_same(const GBOX *g1, const GBOX *g2)
{
	if (FLAGS_GET_ZM(g1->flags) != FLAGS_GET_ZM(g2->flags))
		return LW_FALSE;

	if (!gbox_same_2d(g1, g2)) return LW_FALSE;

	if (FLAGS_GET_Z(g1->flags) && (g1->zmin != g2->zmin || g1->zmax != g2->zmax))
		return LW_FALSE;
	if (FLAGS_GET_M(g1->flags) && (g1->mmin != g2->mmin || g1->mmax != g2->mmax))
		return LW_FALSE;

	return LW_TRUE;
}

LWCOLLECTION *
lwcollection_construct(uint8_t type, int srid, GBOX *bbox,
                       uint32_t ngeoms, LWGEOM **geoms)
{
	LWCOLLECTION *ret;
	int hasz, hasm;
#ifdef CHECK_LWGEOM_ZM
	char zm;
	uint32_t i;
#endif

	LWDEBUGF(2, "lwcollection_construct called with %d, %d, %p, %d, %p.", type, srid, bbox, ngeoms, geoms);

	if( ! lwtype_is_collection(type) )
		lwerror("Non-collection type specified in collection constructor!");

	hasz = 0;
	hasm = 0;
	if ( ngeoms > 0 )
	{
		hasz = FLAGS_GET_Z(geoms[0]->flags);
		hasm = FLAGS_GET_M(geoms[0]->flags);
#ifdef CHECK_LWGEOM_ZM
		zm = FLAGS_GET_ZM(geoms[0]->flags);

		LWDEBUGF(3, "lwcollection_construct type[0]=%d", geoms[0]->type);

		for (i=1; i<ngeoms; i++)
		{
			LWDEBUGF(3, "lwcollection_construct type=[%d]=%d", i, geoms[i]->type);

			if ( zm != FLAGS_GET_ZM(geoms[i]->flags) )
				lwerror("lwcollection_construct: mixed dimension geometries: %d/%d", zm, FLAGS_GET_ZM(geoms[i]->flags));
		}
#endif
	}


	ret = lwalloc(sizeof(LWCOLLECTION));
	ret->type = type;
	ret->flags = gflags(hasz,hasm,0);
	FLAGS_SET_BBOX(ret->flags, bbox?1:0);
	ret->srid = srid;
	ret->ngeoms = ngeoms;
	ret->maxgeoms = ngeoms;
	ret->geoms = geoms;
	ret->bbox = bbox;

	return ret;
}

static size_t
asgeojson_line_buf(const LWLINE *line, char *srs, char *output, GBOX *bbox, int precision)
{
	char *ptr=output;

	ptr += sprintf(ptr, "{\"type\":\"LineString\",");
	if (srs) ptr += asgeojson_srs_buf(ptr, srs);
	if (bbox) ptr += asgeojson_bbox_buf(ptr, bbox, FLAGS_GET_Z(line->flags), precision);
	ptr += sprintf(ptr, "\"coordinates\":[");
	ptr += pointArray_to_geojson(line->points, ptr, precision);
	ptr += sprintf(ptr, "]}");

	return (ptr-output);
}

static size_t
asgeojson_point_buf(const LWPOINT *point, char *srs, char *output, GBOX *bbox, int precision)
{
	char *ptr = output;

	ptr += sprintf(ptr, "{\"type\":\"Point\",");
	if (srs) ptr += asgeojson_srs_buf(ptr, srs);
	if (bbox) ptr += asgeojson_bbox_buf(ptr, bbox, FLAGS_GET_Z(point->flags), precision);

	ptr += sprintf(ptr, "\"coordinates\":");
	ptr += pointArray_to_geojson(point->point, ptr, precision);
	ptr += sprintf(ptr, "}");

	return (ptr-output);
}

static POINTARRAY * ptarray_set_effective_area(POINTARRAY *inpts,int avoid_collaps,int set_area, double trshld)
{
	LWDEBUG(2, "Entered  ptarray_set_effective_area");
	int p;
	POINT4D pt;
	EFFECTIVE_AREAS *ea;
	POINTARRAY *opts;
	int set_m;
	if(set_area)
		set_m=1;
	else
		set_m=FLAGS_GET_M(inpts->flags);
	ea=initiate_effectivearea(inpts);

	opts = ptarray_construct_empty(FLAGS_GET_Z(inpts->flags), set_m, inpts->npoints);

	ptarray_calc_areas(ea,avoid_collaps,set_area,trshld);	
	
	if(set_area)
	{
		/*Only return points with an effective area above the threashold*/
		for (p=0;p<ea->inpts->npoints;p++)
		{
			if(ea->res_arealist[p]>=trshld)
			{
				pt=getPoint4d(ea->inpts, p);
				pt.m=ea->res_arealist[p];
				ptarray_append_point(opts, &pt, LW_TRUE);
			}
		}
	}
	else
	{	
		/*Only return points with an effective area above the threashold*/
		for (p=0;p<ea->inpts->npoints;p++)
		{
			if(ea->res_arealist[p]>=trshld)
			{
				pt=getPoint4d(ea->inpts, p);
				ptarray_append_point(opts, &pt, LW_TRUE);
			}
		}	
	}
	destroy_effectivearea(ea);
	
	return opts;
	
}

/*
 * @param icompound input compound curve
 * @param tol tolerance, semantic driven by tolerance_type
 * @param tolerance_type see LW_LINEARIZE_TOLERANCE_TYPE
 * @param flags see flags in lwarc_linearize
 *
 * @return a newly allocated LWLINE
 */
static LWLINE *
lwcompound_linearize(const LWCOMPOUND *icompound, double tol,
                      LW_LINEARIZE_TOLERANCE_TYPE tolerance_type,
                      int flags)
{
	LWGEOM *geom;
	POINTARRAY *ptarray = NULL, *ptarray_out = NULL;
	LWLINE *tmp = NULL;
	uint32_t i, j;
	POINT4D p;

	LWDEBUG(2, "lwcompound_stroke called.");

	ptarray = ptarray_construct_empty(FLAGS_GET_Z(icompound->flags), FLAGS_GET_M(icompound->flags), 64);

	for (i = 0; i < icompound->ngeoms; i++)
	{
		geom = icompound->geoms[i];
		if (geom->type == CIRCSTRINGTYPE)
		{
			tmp = lwcircstring_linearize((LWCIRCSTRING *)geom, tol, tolerance_type, flags);
			for (j = 0; j < tmp->points->npoints; j++)
			{
				getPoint4d_p(tmp->points, j, &p);
				ptarray_append_point(ptarray, &p, LW_TRUE);
			}
			lwline_free(tmp);
		}
		else if (geom->type == LINETYPE)
		{
			tmp = (LWLINE *)geom;
			for (j = 0; j < tmp->points->npoints; j++)
			{
				getPoint4d_p(tmp->points, j, &p);
				ptarray_append_point(ptarray, &p, LW_TRUE);
			}
		}
		else
		{
			lwerror("Unsupported geometry type %d found.",
			        geom->type, lwtype_name(geom->type));
			return NULL;
		}
	}
	ptarray_out = ptarray_remove_repeated_points(ptarray, 0.0);
	ptarray_free(ptarray);
	return lwline_construct(icompound->srid, NULL, ptarray_out);
}

LWCOLLECTION* lwcollection_simplify(const LWCOLLECTION *igeom, double dist)
{
 	int i;
	LWCOLLECTION *out = lwcollection_construct_empty(igeom->type, igeom->srid, FLAGS_GET_Z(igeom->flags), FLAGS_GET_M(igeom->flags));

	if( lwcollection_is_empty(igeom) )
		return out; /* should we return NULL instead ? */

	for( i = 0; i < igeom->ngeoms; i++ )
	{
		LWGEOM *ngeom = lwgeom_simplify(igeom->geoms[i], dist);
		if ( ngeom ) out = lwcollection_add_lwgeom(out, ngeom);
	}

	return out;
}

/*
 * Construct a LWCIRCSTRING from an array of LWPOINTs
 * LWCIRCSTRING dimensions are large enough to host all input dimensions.
 */
LWCIRCSTRING *
lwcircstring_from_lwpointarray(int srid, uint32_t npoints, LWPOINT **points)
{
	int zmflag=0;
	uint32_t i;
	POINTARRAY *pa;
	uint8_t *newpoints, *ptr;
	size_t ptsize, size;

	/*
	 * Find output dimensions, check integrity
	 */
	for (i = 0; i < npoints; i++)
	{
		if (points[i]->type != POINTTYPE)
		{
			lwerror("lwcurve_from_lwpointarray: invalid input type: %s",
			        lwtype_name(points[i]->type));
			return NULL;
		}
		if (FLAGS_GET_Z(points[i]->flags)) zmflag |= 2;
		if (FLAGS_GET_M(points[i]->flags)) zmflag |= 1;
		if (zmflag == 3) break;
	}

	if (zmflag == 0) ptsize = 2 * sizeof(double);
	else if (zmflag == 3) ptsize = 4 * sizeof(double);
	else ptsize = 3 * sizeof(double);

	/*
	 * Allocate output points array
	 */
	size = ptsize * npoints;
	newpoints = lwalloc(size);
	memset(newpoints, 0, size);

	ptr = newpoints;
	for (i = 0; i < npoints; i++)
	{
		size = ptarray_point_size(points[i]->point);
		memcpy(ptr, getPoint_internal(points[i]->point, 0), size);
		ptr += ptsize;
	}
	pa = ptarray_construct_reference_data(zmflag&2, zmflag&1, npoints, newpoints);
	
	return lwcircstring_construct(srid, NULL, pa);
}