pointcloud/lib/pc_patch.c

/***********************************************************************
* pc_patch.c
*
*  Pointclound patch handling. Create, get and set values from the
*  basic PCPATCH structure.
*
*  PgSQL Pointcloud is free and open source software provided
*  by the Government of Canada
*  Copyright (c) 2013 Natural Resources Canada
*
***********************************************************************/

#include <math.h>
#include <assert.h>
#include "pc_api_internal.h"

int
pc_patch_compute_extent(PCPATCH *pa)
{
	if ( ! pa ) return PC_FAILURE;
	switch ( pa->type )
	{
	case PC_NONE:
		return pc_patch_uncompressed_compute_extent((PCPATCH_UNCOMPRESSED*)pa);
	case PC_DIMENSIONAL:
		return pc_patch_dimensional_compute_extent((PCPATCH_DIMENSIONAL*)pa);
	case PC_LAZPERF:
		return pc_patch_lazperf_compute_extent((PCPATCH_LAZPERF*)pa);
	}
	return PC_FAILURE;
}

/**
* Calculate or re-calculate statistics for a patch.
*/
int
pc_patch_compute_stats(PCPATCH *pa)
{
	if ( ! pa ) return PC_FAILURE;

	switch ( pa->type )
	{
	case PC_NONE:
		return pc_patch_uncompressed_compute_stats((PCPATCH_UNCOMPRESSED*)pa);

	case PC_DIMENSIONAL:
	{
		PCPATCH_UNCOMPRESSED *pu = pc_patch_uncompressed_from_dimensional((PCPATCH_DIMENSIONAL*)pa);
		pc_patch_uncompressed_compute_stats(pu);
		pa->stats = pu->stats; pu->stats = NULL;
		pc_patch_uncompressed_free(pu);
		return PC_SUCCESS;
	}
	case PC_LAZPERF:
	{
		PCPATCH_UNCOMPRESSED *pu = pc_patch_uncompressed_from_lazperf((PCPATCH_LAZPERF*)pa);
		pc_patch_uncompressed_compute_stats(pu);
		pa->stats = pc_stats_clone(pu->stats);
		pc_patch_uncompressed_free(pu);
		return PC_SUCCESS;
	}
	default:
	{
		pcerror("%s: unknown compression type", __func__, pa->type);
		return PC_FAILURE;
	}
	}

	return PC_FAILURE;
}

void
pc_patch_free_stats(PCPATCH *patch)
{
	if ( patch->stats )
	{
		pc_stats_free( patch->stats );
		patch->stats = NULL;
	}
}

void
pc_patch_free(PCPATCH *patch)
{
	switch( patch->type )
	{
	case PC_NONE:
	{
		pc_patch_uncompressed_free((PCPATCH_UNCOMPRESSED*)patch);
		break;
	}
	case PC_DIMENSIONAL:
	{
		pc_patch_dimensional_free((PCPATCH_DIMENSIONAL*)patch);
		break;
	}
	case PC_LAZPERF:
	{
		pc_patch_lazperf_free((PCPATCH_LAZPERF*)patch);
		break;
	}
	default:
	{
		pcerror("%s: unknown compression type %d", __func__, patch->type);
		break;
	}
	}
}


PCPATCH *
pc_patch_from_pointlist(const PCPOINTLIST *ptl)
{
	return (PCPATCH*)pc_patch_uncompressed_from_pointlist(ptl);
}


PCPATCH *
pc_patch_compress(const PCPATCH *patch, void *userdata)
{
	uint32_t schema_compression = patch->schema->compression;
	uint32_t patch_compression = patch->type;

	switch ( schema_compression )
	{
	case PC_DIMENSIONAL:
	{
		if ( patch_compression == PC_NONE )
		{
			/* Dimensionalize, dimensionally compress, return */
			PCPATCH_DIMENSIONAL *pcdu = pc_patch_dimensional_from_uncompressed((PCPATCH_UNCOMPRESSED*)patch);
			PCPATCH_DIMENSIONAL *pcdd = pc_patch_dimensional_compress(pcdu, (PCDIMSTATS*)userdata);
			pc_patch_dimensional_free(pcdu);
			return (PCPATCH*)pcdd;
		}
		else if ( patch_compression == PC_DIMENSIONAL )
		{
			/* Make sure it's compressed, return */
			return (PCPATCH*)pc_patch_dimensional_compress((PCPATCH_DIMENSIONAL*)patch, (PCDIMSTATS*)userdata);
		}
		else if ( patch_compression == PC_LAZPERF )
		{
			PCPATCH_UNCOMPRESSED *pcu = pc_patch_uncompressed_from_lazperf( (PCPATCH_LAZPERF*) patch );
			PCPATCH_DIMENSIONAL *pal = pc_patch_dimensional_from_uncompressed( pcu );
			PCPATCH_DIMENSIONAL *palc = pc_patch_dimensional_compress( pal, NULL );
			pc_patch_dimensional_free(pal);
			return (PCPATCH*) palc;
		}
		else
		{
			pcerror("%s: unknown patch compression type %d", __func__, patch_compression);
		}
	}
	case PC_NONE:
	{
		if ( patch_compression == PC_NONE )
		{
			return (PCPATCH*)patch;
		}
		else if ( patch_compression == PC_DIMENSIONAL )
		{
			PCPATCH_UNCOMPRESSED *pcu = pc_patch_uncompressed_from_dimensional((PCPATCH_DIMENSIONAL*)patch);
			return (PCPATCH*)pcu;

		}
		else if ( patch_compression == PC_LAZPERF )
		{
			PCPATCH_UNCOMPRESSED *pcu = pc_patch_uncompressed_from_lazperf( (PCPATCH_LAZPERF*)patch );
			return (PCPATCH*)pcu;
		}
		else
		{
			pcerror("%s: unknown patch compression type %d", __func__, patch_compression);
		}
	}
	case PC_LAZPERF:
	{
		if ( patch_compression == PC_NONE )
		{
			PCPATCH_LAZPERF *pgc = pc_patch_lazperf_from_uncompressed((PCPATCH_UNCOMPRESSED*)patch);
			if ( ! pgc ) pcerror("%s: lazperf compression failed", __func__);
			return (PCPATCH*)pgc;
		}
		else if ( patch_compression == PC_DIMENSIONAL )
		{
			PCPATCH_UNCOMPRESSED *pad = pc_patch_uncompressed_from_dimensional((PCPATCH_DIMENSIONAL*)patch);
			PCPATCH_LAZPERF *pal = pc_patch_lazperf_from_uncompressed( (PCPATCH_UNCOMPRESSED*) pad );
			pc_patch_uncompressed_free( pad );
			return (PCPATCH*)pal;
		}
		else if ( patch_compression == PC_LAZPERF )
		{
			return (PCPATCH*)patch;
		}
		else
		{
			pcerror("%s: unknown patch compression type %d", __func__, patch_compression);
		}
	}
	default:
	{
		pcerror("%s: unknown schema compression type %d", __func__, schema_compression);
	}
	}

	pcerror("%s: fatal error", __func__);
	return NULL;
}


PCPATCH *
pc_patch_uncompress(const PCPATCH *patch)
{
	uint32_t patch_compression = patch->type;

	if ( patch_compression == PC_DIMENSIONAL )
	{
		PCPATCH_UNCOMPRESSED *pu = pc_patch_uncompressed_from_dimensional((PCPATCH_DIMENSIONAL*)patch);
		return (PCPATCH*)pu;
	}

	if ( patch_compression == PC_NONE )
	{
		return (PCPATCH*)patch;
	}

	if ( patch_compression == PC_LAZPERF )
	{
		PCPATCH_UNCOMPRESSED *pu = pc_patch_uncompressed_from_lazperf( (PCPATCH_LAZPERF*)patch );
		return (PCPATCH*) pu;
	}

	return NULL;
}



PCPATCH *
pc_patch_from_wkb(const PCSCHEMA *s, uint8_t *wkb, size_t wkbsize)
{
	/*
	byte:	  endianness (1 = NDR, 0 = XDR)
	uint32:   pcid (key to POINTCLOUD_SCHEMAS)
	uint32:   compression (0 = no compression, 1 = dimensional, 2 = lazperf)
	uchar[]:  data (interpret relative to pcid and compression)
	*/
	uint32_t compression, pcid;
	PCPATCH *patch;

	if ( ! wkbsize )
	{
		pcerror("%s: zero length wkb", __func__);
	}

	/*
	* It is possible for the WKB compression to be different from the
	* schema compression at this point. The schema compression is only
	* forced at serialization time.
	*/
	pcid = pc_wkb_get_pcid(wkb);
	compression = wkb_get_compression(wkb);

	if ( pcid != s->pcid )
	{
		pcerror("%s: wkb pcid (%d) not consistent with schema pcid (%d)", __func__, pcid, s->pcid);
	}

	switch ( compression )
	{
	case PC_NONE:
	{
		patch = pc_patch_uncompressed_from_wkb(s, wkb, wkbsize);
		break;
	}
	case PC_DIMENSIONAL:
	{
		patch = pc_patch_dimensional_from_wkb(s, wkb, wkbsize);
		break;
	}
	case PC_LAZPERF:
	{
		patch = pc_patch_lazperf_from_wkb(s, wkb, wkbsize);
		break;
	}
	default:
	{
		/* Don't get here */
		pcerror("%s: unknown compression '%d' requested", __func__, compression);
		return NULL;
	}
	}

	if ( PC_FAILURE == pc_patch_compute_extent(patch) )
		pcerror("%s: pc_patch_compute_extent failed", __func__);

	if ( PC_FAILURE == pc_patch_compute_stats(patch) )
		pcerror("%s: pc_patch_compute_stats failed", __func__);

	return patch;

}



uint8_t *
pc_patch_to_wkb(const PCPATCH *patch, size_t *wkbsize)
{
	/*
	byte:	  endianness (1 = NDR, 0 = XDR)
	uint32:   pcid (key to POINTCLOUD_SCHEMAS)
	uint32:   compression (0 = no compression, 1 = dimensional, 2 = lazperf)
	uchar[]:  data (interpret relative to pcid and compression)
	*/
	switch ( patch->type )
	{
	case PC_NONE:
	{
		return pc_patch_uncompressed_to_wkb((PCPATCH_UNCOMPRESSED*)patch, wkbsize);
	}
	case PC_DIMENSIONAL:
	{
		return pc_patch_dimensional_to_wkb((PCPATCH_DIMENSIONAL*)patch, wkbsize);
	}
	case PC_LAZPERF:
	{
		return pc_patch_lazperf_to_wkb((PCPATCH_LAZPERF*)patch, wkbsize);
	}
	}
	pcerror("%s: unknown compression requested '%d'", __func__, patch->schema->compression);
	return NULL;
}

char *
pc_patch_to_string(const PCPATCH *patch)
{
	switch( patch->type )
	{
	case PC_NONE:
		return pc_patch_uncompressed_to_string((PCPATCH_UNCOMPRESSED*)patch);
	case PC_DIMENSIONAL:
		return pc_patch_dimensional_to_string((PCPATCH_DIMENSIONAL*)patch);
	case PC_LAZPERF:
		return pc_patch_lazperf_to_string( (PCPATCH_LAZPERF*)patch );
	}
	pcerror("%s: unsupported compression %d requested", __func__, patch->type);
	return NULL;
}





PCPATCH *
pc_patch_from_patchlist(PCPATCH **palist, int numpatches)
{
	int i;
	uint32_t totalpoints = 0;
	PCPATCH_UNCOMPRESSED *paout;
	const PCSCHEMA *schema = NULL;
	uint8_t *buf;

	assert(palist);
	assert(numpatches);

	/* All schemas better be the same... */
	schema = palist[0]->schema;

	/* How many points will this output have? */
	for ( i = 0; i < numpatches; i++ )
	{
		if ( schema->pcid != palist[i]->schema->pcid )
		{
			pcerror("%s: inconsistent schemas in input", __func__);
			return NULL;
		}
		totalpoints += palist[i]->npoints;
	}

	/* Blank output */
	paout = pc_patch_uncompressed_make(schema, totalpoints);
	buf = paout->data;

	/* Uncompress dimensionals, copy uncompressed */
	for ( i = 0; i < numpatches; i++ )
	{
		const PCPATCH *pa = palist[i];

		/* Update bounds */
		pc_bounds_merge(&(paout->bounds), &(pa->bounds));

		switch ( pa->type )
		{
		case PC_DIMENSIONAL:
		{
			PCPATCH_UNCOMPRESSED *pu = pc_patch_uncompressed_from_dimensional((const PCPATCH_DIMENSIONAL*)pa);
			size_t sz = pu->schema->size * pu->npoints;
			memcpy(buf, pu->data, sz);
			buf += sz;
			pc_patch_free((PCPATCH*)pu);
			break;
		}
		case PC_NONE:
		{
			PCPATCH_UNCOMPRESSED *pu = (PCPATCH_UNCOMPRESSED*)pa;
			size_t sz = pu->schema->size * pu->npoints;
			memcpy(buf, pu->data, sz);
			buf += sz;
			break;
		}
		case PC_LAZPERF:
		{
			PCPATCH_UNCOMPRESSED *pu = pc_patch_uncompressed_from_lazperf((const PCPATCH_LAZPERF*)pa);
			size_t sz = pu->schema->size * pu->npoints;
			memcpy(buf, pu->data, sz);
			buf += sz;
			pc_patch_uncompressed_free(pu);
			break;
		}
		default:
		{
			pcerror("%s: unknown compression type (%d)", __func__, pa->type);
			break;
		}
		}
	}

	paout->npoints = totalpoints;

	if ( PC_FAILURE == pc_patch_uncompressed_compute_stats(paout) )
	{
		pcerror("%s: stats computation failed", __func__);
		return NULL;
	}

	return (PCPATCH*)paout;
}

// first: the first element to select (1-based indexing)
// count: the number of points to select
PCPATCH *
pc_patch_range(const PCPATCH *pa, int first, int count)
{
	PCPATCH_UNCOMPRESSED *paout, *pu;
	int countmax;
	uint8_t *buf;
	size_t size;
	size_t start;

	assert(pa);

	first--;
	countmax = pa->npoints - first;

	if ( count > countmax )
		count = countmax;

	if ( first < 0 || count <= 0 )
		return NULL;

	if ( count == pa->npoints )
		return (PCPATCH *) pa;

	paout = pc_patch_uncompressed_make(pa->schema, count);
	if ( !paout )
		return NULL;
	paout->npoints = count;

	pu = (PCPATCH_UNCOMPRESSED *) pc_patch_uncompress(pa);
	if ( !pu )
	{
		pc_patch_free((PCPATCH *) paout);
		return NULL;
	}

	buf = paout->data;
	start = pa->schema->size * first;
	size = pa->schema->size * count;

	memcpy(buf, pu->data + start, size);

	if ( ((PCPATCH *) pu) != pa )
		pc_patch_free((PCPATCH *) pu);

	if ( PC_FAILURE == pc_patch_uncompressed_compute_extent(paout) )
	{
		pcerror("%s: extent computation failed", __func__);
		pc_patch_free((PCPATCH *) paout);
		return NULL;
	}
	if ( PC_FAILURE == pc_patch_uncompressed_compute_stats(paout) )
	{
		pcerror("%s: stats computation failed", __func__);
		pc_patch_free((PCPATCH *) paout);
		return NULL;
	}

	return (PCPATCH *) paout;
}

/** get point n from patch */
/** positive 1-based:  1=first point,  npoints=last  point */
/** negative 1-based: -1=last  point, -npoints=first point */
PCPOINT *pc_patch_pointn(const PCPATCH *patch, int n)
{
	if(!patch) return NULL;
	if(n<0) n = patch->npoints+n; // negative indices count a backward
	else --n; // 1-based => 0-based indexing
	if(n<0 || n>= patch->npoints) return NULL;

	switch( patch->type )
	{
	case PC_NONE:
		return pc_patch_uncompressed_pointn((PCPATCH_UNCOMPRESSED*)patch,n);
	case PC_DIMENSIONAL:
		return pc_patch_dimensional_pointn((PCPATCH_DIMENSIONAL*)patch,n);
	case PC_LAZPERF:
		return pc_patch_lazperf_pointn((PCPATCH_LAZPERF*)patch, n);
	}
	pcerror("%s: unsupported compression %d requested", __func__, patch->type);
	return NULL;
}


static void
pc_patch_point_set(
		PCPOINT *p, const uint8_t *data, PCDIMENSION **dims, const uint8_t *def)
{
	size_t i;
	for ( i = 0; i < p->schema->ndims; i++ )
	{
		const PCDIMENSION *ddim = dims[i];
		const PCDIMENSION *pdim = p->schema->dims[i];
		uint8_t *pdata = p->data + pdim->byteoffset;
		const uint8_t *ddata = ddim ?
			data + ddim->byteoffset : def + pdim->byteoffset;
		memcpy(pdata, ddata, pdim->size);
	}
}


/** set schema for patch */
PCPATCH*
pc_patch_set_schema(PCPATCH *patch, const PCSCHEMA *new_schema, double def)
{
	PCDIMENSION** new_dimensions = new_schema->dims;
	PCDIMENSION* old_dimensions[new_schema->ndims];
	const PCSCHEMA *old_schema = patch->schema;
	PCPATCH_UNCOMPRESSED *paout;
	PCPOINT opt, npt;
	PCPATCH *pain;
	PCPOINT *dpt;
	size_t i, j;

	// create a point for storing the default values
	dpt = pc_point_make(new_schema);

	for ( j = 0; j < new_schema->ndims; j++ )
	{
		PCDIMENSION *ndim = new_dimensions[j];
		PCDIMENSION *odim = pc_schema_get_dimension_by_name(
				old_schema, ndim->name);
		old_dimensions[j] = odim;
		if ( odim )
		{
			if ( ndim->interpretation != odim->interpretation )
			{
				pcerror("dimension interpretations are not matching");
				pc_point_free(dpt);
				return NULL;
			}
		}
		else
		{
			pc_point_set_double(dpt, ndim, def);
		}
	}

	pain = pc_patch_uncompress(patch);
	paout = pc_patch_uncompressed_make(new_schema, patch->npoints);
	paout->npoints = pain->npoints;

	opt.schema = old_schema;
	npt.schema = new_schema;
	opt.readonly = PC_TRUE;
	npt.readonly = PC_TRUE;

	opt.data = ((PCPATCH_UNCOMPRESSED *) pain)->data;
	npt.data = paout->data;

	for ( i = 0; i < patch->npoints; i++ )
	{
		pc_patch_point_set(&npt, opt.data, old_dimensions, dpt->data);
		opt.data += old_schema->size;
		npt.data += new_schema->size;
	}

	if ( patch->stats )
	{
		paout->stats = pc_stats_new(new_schema);

		opt.data = patch->stats->min.data;
		npt.data = paout->stats->min.data;
		pc_patch_point_set(&npt, opt.data, old_dimensions, dpt->data);

		opt.data = patch->stats->max.data;
		npt.data = paout->stats->max.data;
		pc_patch_point_set(&npt, opt.data, old_dimensions, dpt->data);

		opt.data = patch->stats->avg.data;
		npt.data = paout->stats->avg.data;
		pc_patch_point_set(&npt, opt.data, old_dimensions, dpt->data);

		pc_point_get_x(&paout->stats->min, &paout->bounds.xmin);
		pc_point_get_y(&paout->stats->min, &paout->bounds.ymin);
		pc_point_get_x(&paout->stats->max, &paout->bounds.xmax);
		pc_point_get_y(&paout->stats->max, &paout->bounds.ymax);
	}
	else
	{
		double xscale = npt.schema->xdim->scale / opt.schema->xdim->scale;
		double yscale = npt.schema->ydim->scale / opt.schema->ydim->scale;
		double xoffset = npt.schema->xdim->offset - opt.schema->xdim->offset;
		double yoffset = npt.schema->ydim->offset - opt.schema->ydim->offset;

		paout->bounds.xmin = patch->bounds.xmin * xscale + xoffset;
		paout->bounds.xmax = patch->bounds.xmax * xscale + xoffset;
		paout->bounds.ymin = patch->bounds.ymin * yscale + yoffset;
		paout->bounds.xmax = patch->bounds.ymax * yscale + yoffset;
	}

	pc_point_free(dpt);

	if ( pain != patch )
		pc_patch_free(pain);

	return (PCPATCH*) paout;
}


/**
* Read all the points from "patch", and transform them based on "new_schema".
* Return a new patch with the transformed points.
*/
PCPATCH*
pc_patch_transform(const PCPATCH *patch, const PCSCHEMA *new_schema, double def)
{
	PCDIMENSION** new_dimensions = new_schema->dims;
	PCDIMENSION* old_dimensions[new_schema->ndims];
	const PCSCHEMA *old_schema = patch->schema;
	PCPATCH_UNCOMPRESSED *paout;
	PCPOINT opt, npt;
	PCPATCH *pain;
	size_t i, j;

	if ( old_schema->srid != new_schema->srid )
	{
		pcwarn("old and new schemas have different srids, and data "
			   "reprojection is not yet supported");
		return NULL;
	}

	for ( j = 0; j < new_schema->ndims; j++ )
	{
		PCDIMENSION *ndim = new_dimensions[j];
		PCDIMENSION *odim = pc_schema_get_dimension_by_name(
				old_schema, ndim->name);
		old_dimensions[j] = odim;
	}

	pain = pc_patch_uncompress(patch);

	paout = pc_patch_uncompressed_make(new_schema, patch->npoints);
	paout->npoints = pain->npoints;

	opt.schema = old_schema;
	npt.schema = new_schema;
	opt.readonly = PC_TRUE;
	npt.readonly = PC_TRUE;

	opt.data = ((PCPATCH_UNCOMPRESSED *) pain)->data;
	npt.data = paout->data;

	// reinterpret the data and fill the output patch
	//
	// TODO: for the case where the old and new dimension sets don't intersect (all
	// the values in old_dimensions are NULL) a faster path could probably be used
	for ( i = 0; i <patch->npoints; i++ )
	{
		for ( j = 0; j < new_schema->ndims; j++ )
		{
			// pc_point_get_double returns immediately w/o changing val if the
			// dimension it is passed is NULL
			double val = def;
			pc_point_get_double(&opt, old_dimensions[j], &val);
			pc_point_set_double(&npt, new_dimensions[j], val);
		}

		opt.data += old_schema->size;
		npt.data += new_schema->size;
	}

	if ( pain != patch )
		pc_patch_free(pain);

	if ( PC_FAILURE == pc_patch_uncompressed_compute_extent(paout) )
	{
		pcerror("%s: failed to compute patch extent", __func__);
		pc_patch_free((PCPATCH *)paout);
		return NULL;
	}

	if ( PC_FAILURE == pc_patch_uncompressed_compute_stats(paout) )
	{
		pcerror("%s: failed to compute patch stats", __func__);
		pc_patch_free((PCPATCH *)paout);
		return NULL;
	}

	return (PCPATCH*) paout;
}