pointcloud/lib/pc_filter.c

/***********************************************************************
* pc_filter.c
*
*  Pointclound patch filtering.
*
*  Copyright (c) 2013 OpenGeo
*
***********************************************************************/

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

PCBITMAP *
pc_bitmap_new(uint32_t npoints)
{
    PCBITMAP *map = pcalloc(sizeof(PCBITMAP));
    map->map = pcalloc(sizeof(uint8_t)*npoints);
    map->npoints = npoints;
    map->nset = 0;
    return map;
}

void
pc_bitmap_free(PCBITMAP *map)
{
    if (map->map) pcfree(map->map);
    pcfree(map);
}

inline void
pc_bitmap_set(PCBITMAP *map, int i, int val)
{
    uint8_t curval = map->map[i];
    if ( val && (!curval) )
        map->nset++;
    if ( (!val) && curval )
        map->nset--;
    
    map->map[i] = (val!=0);
}

inline uint8_t
pc_bitmap_get(const PCBITMAP *map, int i)
{
    return map->map[i];
}

void
pc_bitmap_filter(PCBITMAP *map, PC_FILTERTYPE filter, int i, double d, double val1, double val2)
{
    switch ( filter )
    {
        case PC_GT:
            pc_bitmap_set(map, i, (d > val1 ? 1 : 0)); 
            break;
        case PC_LT:
            pc_bitmap_set(map, i, (d < val1 ? 1 : 0)); 
            break;
        case PC_EQUAL:
            pc_bitmap_set(map, i, (d == val1 ? 1 : 0)); 
            break;
        case PC_BETWEEN:
            pc_bitmap_set(map, i, (d > val1 && d < val2 ? 1 : 0)); 
            break;
    }
}

static PCBITMAP *
pc_patch_uncompressed_bitmap(const PCPATCH_UNCOMPRESSED *pa, uint32_t dimnum, PC_FILTERTYPE filter, double val1, double val2)
{
    PCPOINT pt;
    uint32_t i = 0;    
    uint8_t *buf = pa->data;
    double d;
    size_t sz = pa->schema->size;
    PCBITMAP *map = pc_bitmap_new(pa->npoints);

    pt.readonly = PC_TRUE;
    pt.schema = pa->schema;
    
    while ( i < pa->npoints )
    {
        pt.data = buf;
        pc_point_get_double(&pt, pa->schema->dims[dimnum], &d);
        /* Apply the filter to the bitmap */
        pc_bitmap_filter(map, filter, i, d, val1, val2);
        /* Advance data pointer and counter */
        buf += sz;
        i++;
    }
    
    return map;
}


static PCPATCH_UNCOMPRESSED *
pc_patch_uncompressed_filter(const PCPATCH_UNCOMPRESSED *pu, const PCBITMAP *map)
{
    int i = 0, j = 0;
    size_t sz = pu->schema->size;
    PCPATCH_UNCOMPRESSED *fpu = pc_patch_uncompressed_make(pu->schema, map->nset);
    uint8_t *buf = pu->data;
    uint8_t *fbuf = fpu->data;
    
    assert(map->npoints == pu->npoints);
    
    while ( i < pu->npoints )
    {
        if ( pc_bitmap_get(map, i) )
        {
            memcpy(fbuf, buf, sz);
            fbuf += sz;
        }
        buf += sz;
        i++;
    }
    
    fpu->maxpoints = fpu->npoints = map->nset;
    
    if ( PC_FAILURE == pc_patch_uncompressed_compute_extent(fpu) )
	{
		pcerror("%s: failed to compute patch extent", __func__);
		return NULL;
	}

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



static PCBITMAP *
pc_patch_dimensional_bitmap(const PCPATCH_DIMENSIONAL *pdl, uint32_t dimnum, PC_FILTERTYPE filter, double val1, double val2)
{
    assert(dimnum < pdl->schema->ndims);
    double unscaled1 = pc_value_unscale_unoffset(val1, pdl->schema->dims[dimnum]);
    double unscaled2 = pc_value_unscale_unoffset(val2, pdl->schema->dims[dimnum]);
    
    return pc_bytes_bitmap(&(pdl->bytes[dimnum]), filter, unscaled1, unscaled2);
}

static PCPATCH_DIMENSIONAL *
pc_patch_dimensional_filter(const PCPATCH_DIMENSIONAL *pdl, const PCBITMAP *map)
{
    int i = 0;
    PCPATCH_DIMENSIONAL *fpdl = pc_patch_dimensional_clone(pdl);
    
    for ( i = 0; i < pdl->schema->ndims; i++ )
    {
        fpdl->bytes[i] = pc_bytes_filter(&(pdl->bytes[i]), map);
    }

    fpdl->npoints = map->nset;
    return fpdl;
}


PCPATCH *
pc_patch_filter(const PCPATCH *pa, uint32_t dimnum, PC_FILTERTYPE filter, double val1, double val2)
{
    if ( ! pa ) return NULL;
    PCPATCH *paout;
    
	switch ( pa->type )
	{
		case PC_NONE:
		{
            PCBITMAP *map = pc_patch_uncompressed_bitmap((PCPATCH_UNCOMPRESSED*)pa, dimnum, filter, val1, val2);
            PCPATCH_UNCOMPRESSED *pu = pc_patch_uncompressed_filter((PCPATCH_UNCOMPRESSED*)pa, map);
            pc_bitmap_free(map);
            paout = (PCPATCH*)pu;
            break;
		}
		case PC_GHT:
		{
            PCPATCH_UNCOMPRESSED *pu = pc_patch_uncompressed_from_ght((PCPATCH_GHT*)pa);
            PCBITMAP *map = pc_patch_uncompressed_bitmap(pu, dimnum, filter, val1, val2);
            PCPATCH_UNCOMPRESSED *pu2 = pc_patch_uncompressed_filter(pu, map);
            PCPATCH_GHT *pgh = pc_patch_ght_from_uncompressed(pu2);
            pc_patch_free((PCPATCH*)pu);
            pc_patch_free((PCPATCH*)pu2);
            paout = (PCPATCH*)pgh;
            break;
        }
		case PC_DIMENSIONAL:
		{
		    PCBITMAP *map = pc_patch_dimensional_bitmap((PCPATCH_DIMENSIONAL*)pa, dimnum, filter, val1, val2);
            PCPATCH_DIMENSIONAL *pdl = pc_patch_dimensional_filter((PCPATCH_DIMENSIONAL*)pa, map);
            pc_bitmap_free(map);
            paout = (PCPATCH*)pdl;
            break;
        }
		default:
            pcerror("%s: failure", __func__);
	}
	
	if ( PC_FAILURE == pc_patch_compute_extent(paout) )
		pcerror("%s: pc_patch_compute_extent failed", __func__);

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

    return paout;
}