effectivearea.c

Go to the documentation of this file.

/**********************************************************************
 *
 * PostGIS - Spatial Types for PostgreSQL
 * http://postgis.net
 *
 * PostGIS is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 2 of the License, or
 * (at your option) any later version.
 *
 * PostGIS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
 *
 **********************************************************************
 *
 * Copyright 2014 Nicklas Avén
 *
 **********************************************************************/
 
 
 #include "effectivearea.h"
 
 
EFFECTIVE_AREAS*
initiate_effectivearea(const POINTARRAY *inpts)
{
        LWDEBUG(2, "Entered  initiate_effectivearea");
        EFFECTIVE_AREAS *ea;
        ea=lwalloc(sizeof(EFFECTIVE_AREAS));
        ea->initial_arealist = lwalloc(inpts->npoints*sizeof(areanode));
        ea->res_arealist = lwalloc(inpts->npoints*sizeof(double));
        ea->inpts=inpts;
        return ea;
}
 
 
void destroy_effectivearea(EFFECTIVE_AREAS *ea)
{
        lwfree(ea->initial_arealist);
        lwfree(ea->res_arealist);
        lwfree(ea);
}
 
 
static MINHEAP
initiate_minheap(int npoints)
{
        MINHEAP tree;
        tree.key_array = lwalloc(npoints*sizeof(void*));
        tree.maxSize=npoints;
        tree.usedSize=0;
        return tree;
}
 
 
static void
destroy_minheap(MINHEAP tree)
{
        lwfree(tree.key_array);
}
 
 
static double triarea2d(const double *P1, const double *P2, const double *P3)
{
        return fabs(0.5*((P1[0]-P2[0])*(P3[1]-P2[1])-(P1[1]-P2[1])*(P3[0]-P2[0])));
}
 
static double triarea3d(const double *P1, const double *P2, const double *P3)
{
        LWDEBUG(2, "Entered  triarea3d");
        double ax,bx,ay,by,az,bz,cx,cy,cz, area;
 
        ax=P1[0]-P2[0];
        bx=P3[0]-P2[0];
        ay=P1[1]-P2[1];
        by=P3[1]-P2[1];
        az=P1[2]-P2[2];
        bz=P3[2]-P2[2];
 
        cx = ay*bz - az*by;
        cy = az*bx - ax*bz;
        cz = ax*by - ay*bx;
 
        area = fabs(0.5*(sqrt(cx*cx+cy*cy+cz*cz)));
        return area;
}
 
static int cmpfunc (const void * a, const void * b)
{
        double v1 =  (*(areanode**)a)->area;
        double v2 = (*(areanode**)b)->area;
        /*qsort gives unpredictable results when comparing identical values.
        If two values is the same we force returning the last point in the point array.
        That way we get the same ordering on different machines and platforms*/
        if (v1==v2)
                return (*(areanode**)a)-(*(areanode**)b);
        else
                return  (v1 > v2) ? 1 : ((v1 < v2) ? -1 : 0);
}
 
 
static void down(MINHEAP *tree,areanode *arealist,int parent)
{
        LWDEBUG(2, "Entered  down");
        areanode **treearray=tree->key_array;
        int left=parent*2+1;
        int right = left +1;
        void *tmp;
        int swap=parent;
        double leftarea=0;
        double rightarea=0;
 
        double parentarea=((areanode*) treearray[parent])->area;
 
        if(left<tree->usedSize)
        {
                leftarea=((areanode*) treearray[left])->area;
                if(parentarea>leftarea)
                        swap=left;
        }
        if(right<tree->usedSize)
        {
                rightarea=((areanode*) treearray[right])->area;
                if(rightarea<parentarea&&rightarea<leftarea)
                        swap=right;
        }
        if(swap>parent)
        {
        /*ok, we have to swap something*/
                tmp=treearray[parent];
                treearray[parent]=treearray[swap];
                /*Update reference*/
                ((areanode*) treearray[parent])->treeindex=parent;
                treearray[swap]=tmp;
                /*Update reference*/
                ((areanode*) treearray[swap])->treeindex=swap;
                if(swap<tree->usedSize)
                        down(tree,arealist,swap);
        }
        return;
}
 
 
static void up(MINHEAP *tree, __attribute__((__unused__)) areanode *e,int c)
{
        LWDEBUG(2, "Entered  up");
        void *tmp;
 
        areanode **treearray=tree->key_array;
 
        int parent=floor((c-1)/2);
 
        while(((areanode*) treearray[c])->area<((areanode*) treearray[parent])->area)
        {
                /*ok, we have to swap*/
                tmp=treearray[parent];
                treearray[parent]=treearray[c];
                /*Update reference*/
                ((areanode*) treearray[parent])->treeindex=parent;
                treearray[c]=tmp;
                /*Update reference*/
                ((areanode*) treearray[c])->treeindex=c;
                c=parent;
                parent=floor((c-1)/2);
        }
        return;
}
 
 
static  areanode* minheap_pop(MINHEAP *tree,areanode *arealist )
{
        LWDEBUG(2, "Entered  minheap_pop");
        areanode *res = tree->key_array[0];
 
        /*put last value first*/
        tree->key_array[0]=tree->key_array[(tree->usedSize)-1];
        ((areanode*) tree->key_array[0])->treeindex=0;
 
        tree->usedSize--;
        down(tree,arealist,0);
        return res;
}
 
 
static void minheap_update(MINHEAP *tree,areanode *arealist , int idx)
{
        areanode **treearray=tree->key_array;
        int     parent=floor((idx-1)/2);
 
        if(((areanode*) treearray[idx])->area<((areanode*) treearray[parent])->area)
                up(tree,arealist,idx);
        else
                down(tree,arealist,idx);
        return;
}
 
static void tune_areas(EFFECTIVE_AREAS *ea, int avoid_collaps, int set_area, double trshld)
{
        LWDEBUG(2, "Entered  tune_areas");
        const double *P1;
        const double *P2;
        const double *P3;
        double area;
        int go_on=1;
        double check_order_min_area = 0;
 
        int npoints=ea->inpts->npoints;
        int i;
        int current, before_current, after_current;
 
        MINHEAP tree = initiate_minheap(npoints);
 
        int is3d = FLAGS_GET_Z(ea->inpts->flags);
 
 
        /*Add all keys (index in initial_arealist) into minheap array*/
        for (i=0;i<npoints;i++)
        {
                tree.key_array[i]=ea->initial_arealist+i;
                LWDEBUGF(2, "add nr %d, with area %lf, and %lf",i,ea->initial_arealist[i].area, tree.key_array[i]->area );
        }
        tree.usedSize=npoints;
 
        /*order the keys by area, small to big*/
        qsort(tree.key_array, npoints, sizeof(void*), cmpfunc);
 
        /*We have to put references to our tree in our point-list*/
        for (i=0;i<npoints;i++)
        {
                 ((areanode*) tree.key_array[i])->treeindex=i;
                LWDEBUGF(4,"Check ordering qsort gives, area=%lf and belong to point %ld",((areanode*) tree.key_array[i])->area, tree.key_array[i]-ea->initial_arealist);
        }
        /*Ok, now we have a minHeap, just need to keep it*/
 
        /*for (i=0;i<npoints-1;i++)*/
        i=0;
        while (go_on)
        {
                /*Get a reference to the point with the currently smallest effective area*/
                current=minheap_pop(&tree, ea->initial_arealist)-ea->initial_arealist;
 
                /*We have found the smallest area. That is the resulting effective area for the "current" point*/
                if (i<npoints-avoid_collaps)
                        ea->res_arealist[current]=ea->initial_arealist[current].area;
                else
                        ea->res_arealist[current]=FLT_MAX;
 
                if(ea->res_arealist[current]<check_order_min_area)
                        lwerror("Oh no, this is a bug. For some reason the minHeap returned our points in the wrong order. Please file a ticket in PostGIS ticket system, or send a mail at the mailing list.Returned area = %lf, and last area = %lf",ea->res_arealist[current],check_order_min_area);
 
                check_order_min_area=ea->res_arealist[current];
 
                /*The found smallest area point is now regarded as eliminated and we have to recalculate the area the adjacent (ignoring earlier eliminated points) points gives*/
 
                /*FInd point before and after*/
                before_current=ea->initial_arealist[current].prev;
                after_current=ea->initial_arealist[current].next;
 
                P2= (double*)getPoint_internal(ea->inpts, before_current);
                P3= (double*)getPoint_internal(ea->inpts, after_current);
 
                /*Check if point before current point is the first in the point array. */
                if(before_current>0)
                {
 
                        P1= (double*)getPoint_internal(ea->inpts, ea->initial_arealist[before_current].prev);
                        if(is3d)
                                area=triarea3d(P1, P2, P3);
                        else
                                area=triarea2d(P1, P2, P3);
 
                        ea->initial_arealist[before_current].area = FP_MAX(area,ea->res_arealist[current]);
                        minheap_update(&tree, ea->initial_arealist, ea->initial_arealist[before_current].treeindex);
                }
                if(after_current<npoints-1)/*Check if point after current point is the last in the point array. */
                {
                        P1=P2;
                        P2=P3;
 
                        P3= (double*)getPoint_internal(ea->inpts, ea->initial_arealist[after_current].next);
 
 
                        if(is3d)
                                area=triarea3d(P1, P2, P3);
                        else
                                area=triarea2d(P1, P2, P3);
 
 
                        ea->initial_arealist[after_current].area = FP_MAX(area,ea->res_arealist[current]);
                        minheap_update(&tree, ea->initial_arealist, ea->initial_arealist[after_current].treeindex);
                }
 
                /*rearrange the nodes so the eliminated point will be ignored on the next run*/
                ea->initial_arealist[before_current].next = ea->initial_arealist[current].next;
                ea->initial_arealist[after_current].prev = ea->initial_arealist[current].prev;
 
                /*Check if we are finished*/
                if((!set_area && ea->res_arealist[current]>=trshld) || (ea->initial_arealist[0].next==(npoints-1)))
                        go_on=0;
 
                i++;
        };
        destroy_minheap(tree);
        return;
}
 
 
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 POINTARRAY * ptarray_set_effective_area(POINTARRAY *inpts,int avoid_collaps,int set_area, double trshld)
{
        LWDEBUG(2, "Entered  ptarray_set_effective_area");
        uint32_t 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 threshold*/
                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 threshold*/
                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;
 
}
 
static LWLINE* lwline_set_effective_area(const LWLINE *iline,int set_area, double trshld)
{
        LWDEBUG(2, "Entered  lwline_set_effective_area");
 
                /* Skip empty case or too small to simplify */
        if( lwline_is_empty(iline) || iline->points->npoints<3)
                return lwline_clone(iline);
 
        int set_m;
        if(set_area)
                set_m=1;
        else
                set_m=FLAGS_GET_M(iline->flags);
 
        LWLINE *oline = lwline_construct_empty(iline->srid, FLAGS_GET_Z(iline->flags), set_m);
 
 
 
        oline = lwline_construct(iline->srid, NULL, ptarray_set_effective_area(iline->points,2,set_area,trshld));
 
        oline->type = iline->type;
        return oline;
 
}
 
 
static LWPOLY* lwpoly_set_effective_area(const LWPOLY *ipoly,int set_area, double trshld)
{
        LWDEBUG(2, "Entered  lwpoly_set_effective_area");
        uint32_t i;
        int set_m;
        int avoid_collapse=4;
        if(set_area)
                set_m=1;
        else
                set_m=FLAGS_GET_M(ipoly->flags);
        LWPOLY *opoly = lwpoly_construct_empty(ipoly->srid, FLAGS_GET_Z(ipoly->flags), set_m);
 
        if( lwpoly_is_empty(ipoly) )
                return opoly; /* should we return NULL instead ? */
 
        for (i = 0; i < ipoly->nrings; i++)
        {
                POINTARRAY *pa = ptarray_set_effective_area(ipoly->rings[i],avoid_collapse,set_area,trshld);
                /* Add ring to simplified polygon */
                if(pa->npoints>=4)
                {
                        if( lwpoly_add_ring(opoly,pa ) == LW_FAILURE )
                                return NULL;
                }
                /*Inner rings we allow to collapse and then we remove them*/
                avoid_collapse=0;
        }
 
 
        opoly->type = ipoly->type;
 
        if( lwpoly_is_empty(opoly))
                return NULL;
 
        return opoly;
 
}
 
 
static LWCOLLECTION* lwcollection_set_effective_area(const LWCOLLECTION *igeom,int set_area, double trshld)
{
        LWDEBUG(2, "Entered  lwcollection_set_effective_area");
        uint32_t i;
        int set_m;
        if(set_area)
                set_m=1;
        else
                set_m=FLAGS_GET_M(igeom->flags);
        LWCOLLECTION *out = lwcollection_construct_empty(igeom->type, igeom->srid, FLAGS_GET_Z(igeom->flags), set_m);
 
        if( lwcollection_is_empty(igeom) )
                return out; /* should we return NULL instead ? */
 
        for( i = 0; i < igeom->ngeoms; i++ )
        {
                LWGEOM *ngeom = lwgeom_set_effective_area(igeom->geoms[i],set_area,trshld);
                if ( ngeom ) out = lwcollection_add_lwgeom(out, ngeom);
        }
 
        return out;
}
 
 
LWGEOM* lwgeom_set_effective_area(const LWGEOM *igeom,int set_area, double trshld)
{
        LWDEBUG(2, "Entered  lwgeom_set_effective_area");
        switch (igeom->type)
        {
        case POINTTYPE:
        case MULTIPOINTTYPE:
                return lwgeom_clone(igeom);
        case LINETYPE:
                return (LWGEOM*)lwline_set_effective_area((LWLINE*)igeom,set_area, trshld);
        case POLYGONTYPE:
                return (LWGEOM*)lwpoly_set_effective_area((LWPOLY*)igeom,set_area, trshld);
        case MULTILINETYPE:
        case MULTIPOLYGONTYPE:
        case COLLECTIONTYPE:
                return (LWGEOM*)lwcollection_set_effective_area((LWCOLLECTION *)igeom,set_area, trshld);
        default:
                lwerror("lwgeom_simplify: unsupported geometry type: %s",lwtype_name(igeom->type));
        }
        return NULL;
}