self-hosted/postgis
self-hosted/postgis
1
0
Fork 0
mirror of https://git.osgeo.org/gitea/postgis/postgis synced 2024-10-24 09:02:37 +00:00
Code Issues Packages Projects Releases Wiki Activity

3D Distance functions, only point-point and point line. #576

Browse source 
git-svn-id: http://svn.osgeo.org/postgis/trunk@5889 b70326c6-7e19-0410-871a-916f4a2858ee
This commit is contained in:
Nicklas Avén 2010-09-01 19:55:11 +00:00
parent 12b8aedf60
commit d5b2e25d6d
9 changed files with 1000 additions and 22 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
liblwgeom/Makefile.in
View file

@ -21,6 +21,7 @@ LEX=@LEX@
SA_OBJS = \
stringbuffer.o \
measures.o \
measures3d.o \
box2d.o \
ptarray.o \
lwgeom_api.o \

20
liblwgeom/liblwgeom.h
View file

@ -1299,10 +1299,11 @@ extern double next_double_up(float d);
#define LW_ABS(a) ((a) < (0) ? -(a) : (a))
/* for the measure functions*/
#define DIST2D_MAX -1
#define DIST2D_MIN 1
#define DIST_MAX -1
#define DIST_MIN 1
/* general utilities */
/* general utilities
2D*/
extern double distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2);
extern double distance2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B);
extern LWGEOM *lw_dist2d_distancepoint(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode);
@ -1311,6 +1312,19 @@ extern double lwgeom_mindistance2d(LWGEOM *lw1, LWGEOM *lw2);
extern double lwgeom_mindistance2d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance);
extern double lwgeom_maxdistance2d(LWGEOM *lw1, LWGEOM *lw2);
extern double lwgeom_maxdistance2d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance);
/*
3D*/
extern double distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2);
extern double distance3d_pt_seg(const POINT3D *p, const POINT3D *A, const POINT3D *B);
extern LWGEOM *lw_dist3d_distancepoint(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode);
extern LWGEOM *lw_dist3d_distanceline(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode);
extern double lwgeom_mindistance3d(LWGEOM *lw1, LWGEOM *lw2);
extern double lwgeom_mindistance3d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance);
extern double lwgeom_maxdistance3d(LWGEOM *lw1, LWGEOM *lw2);
extern double lwgeom_maxdistance3d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance);
extern double lwgeom_polygon_area(const LWPOLY *poly);
extern double lwgeom_polygon_perimeter(const LWPOLY *poly);
extern double lwgeom_polygon_perimeter2d(const LWPOLY *poly);

32
liblwgeom/measures.c
View file

@ -31,7 +31,7 @@ lw_dist2d_distanceline(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode)
{
double x1,x2,y1,y2;
double initdistance = ( mode == DIST2D_MIN ? MAXFLOAT : -1.0);
double initdistance = ( mode == DIST_MIN ? MAXFLOAT : -1.0);
DISTPTS thedl;
LWPOINT *lwpoints[2];
LWGEOM *result;
@ -129,7 +129,7 @@ lwgeom_maxdistance2d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
/*double thedist;*/
DISTPTS thedl;
LWDEBUG(2, "lwgeom_maxdistance2d_tolerance is called");
thedl.mode = DIST2D_MAX;
thedl.mode = DIST_MAX;
thedl.distance= -1;
thedl.tolerance = tolerance;
if (lw_dist2d_comp( lw1,lw2,&thedl))
@ -160,7 +160,7 @@ lwgeom_mindistance2d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
{
DISTPTS thedl;
LWDEBUG(2, "lwgeom_mindistance2d_tolerance is called");
thedl.mode = DIST2D_MIN;
thedl.mode = DIST_MIN;
thedl.distance= MAXFLOAT;
thedl.tolerance = tolerance;
if (lw_dist2d_comp( lw1,lw2,&thedl))
@ -266,10 +266,10 @@ int lw_dist2d_recursive(const LWCOLLECTION * lwg1,const LWCOLLECTION * lwg2,DIST
/*If one of geometries is empty, return. True here only means continue searching. False would have stoped the process*/
if (lwgeom_is_empty(g1)||lwgeom_is_empty(g2)) return LW_TRUE;
if ((dl->mode==DIST2D_MAX)||(TYPE_GETTYPE(g1->type)==POINTTYPE) ||(TYPE_GETTYPE(g2->type)==POINTTYPE)||lw_dist2d_check_overlap(g1, g2))
if ((dl->mode==DIST_MAX)||(TYPE_GETTYPE(g1->type)==POINTTYPE) ||(TYPE_GETTYPE(g2->type)==POINTTYPE)||lw_dist2d_check_overlap(g1, g2))
{
if (!lw_dist2d_distribute_bruteforce(g1, g2, dl)) return LW_FALSE;
if (dl->distance<=dl->tolerance && dl->mode == DIST2D_MIN) return LW_TRUE; /*just a check if the answer is already given*/
if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
}
else
{
@ -495,7 +495,7 @@ lw_dist2d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl)
getPoint2d_p(point->point, 0, &p);
if (dl->mode == DIST2D_MAX)
if (dl->mode == DIST_MAX)
{
LWDEBUG(3, "looking for maxdistance");
return lw_dist2d_pt_ptarray(&p, poly->rings[0], dl);
@ -524,7 +524,7 @@ lw_dist2d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl)
}
LWDEBUG(3, " inside the polygon");
if (dl->mode == DIST2D_MIN)
if (dl->mode == DIST_MIN)
{
dl->distance=0.0;
dl->p1.x=p.x;
@ -576,7 +576,7 @@ lw_dist2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl)
LWDEBUG(2, "lw_dist2d_poly_poly called");
/*1 if we are looking for maxdistance, just check the outer rings.*/
if (dl->mode == DIST2D_MAX)
if (dl->mode == DIST_MAX)
{
return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[0],dl);
}
@ -667,7 +667,7 @@ lw_dist2d_pt_ptarray(POINT2D *p, POINTARRAY *pa,DISTPTS *dl)
getPoint2d_p(pa, t, &end);
if (!lw_dist2d_pt_seg(p, &start, &end,dl)) return LW_FALSE;
if (dl->distance<=dl->tolerance && dl->mode == DIST2D_MIN) return LW_TRUE; /*just a check if the answer is already given*/
if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
start = end;
}
@ -705,7 +705,7 @@ lw_dist2d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, DISTPTS *dl)
LWDEBUGF(3, " distance from ring %d: %f, mindist: %f",
i, dl->distance, dl->tolerance);
if (dl->distance<=dl->tolerance && dl->mode == DIST2D_MIN) return LW_TRUE; /*just a check if the answer is already given*/
if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
}
/*
@ -738,7 +738,7 @@ lw_dist2d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, DISTPTS *dl)
return LW_TRUE;
}
}
if (dl->mode == DIST2D_MIN)
if (dl->mode == DIST_MIN)
{
dl->distance=0.0;
dl->p1.x=pt.x;
@ -765,7 +765,7 @@ lw_dist2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2,DISTPTS *dl)
LWDEBUGF(2, "lw_dist2d_ptarray_ptarray called (points: %d-%d)",l1->npoints, l2->npoints);
if (dl->mode == DIST2D_MAX)/*If we are searching for maxdistance we go straight to point-point calculation since the maxdistance have to be between two vertexes*/
if (dl->mode == DIST_MAX)/*If we are searching for maxdistance we go straight to point-point calculation since the maxdistance have to be between two vertexes*/
{
for (t=0; t<l1->npoints; t++) /*for each segment in L1 */
{
@ -795,7 +795,7 @@ lw_dist2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2,DISTPTS *dl)
LWDEBUGF(4, "mindist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
t, u, dl->distance, dl->tolerance);
if (dl->distance<=dl->tolerance && dl->mode == DIST2D_MIN) return LW_TRUE; /*just a check if the answer is already given*/
if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
start2 = end2;
}
start = end;
@ -875,7 +875,7 @@ lw_dist2d_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl)
s = s_top/s_bot;
r= r_top/r_bot;
if (((r<0) || (r>1) || (s<0) || (s>1)) || (dl->mode == DIST2D_MAX))
if (((r<0) || (r>1) || (s<0) || (s>1)) || (dl->mode == DIST_MAX))
{
if ((lw_dist2d_pt_seg(A,C,D,dl)) && (lw_dist2d_pt_seg(B,C,D,dl)))
{
@ -889,7 +889,7 @@ lw_dist2d_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl)
}
else
{
if (dl->mode == DIST2D_MIN) /*If there is intersection we identify the intersection point and return it but only if we are looking for mindistance*/
if (dl->mode == DIST_MIN) /*If there is intersection we identify the intersection point and return it but only if we are looking for mindistance*/
{
POINT2D theP;
@ -1283,7 +1283,7 @@ lw_dist2d_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B, DISTPTS *dl)
the maxdistance have to be between two vertexes,
compared to mindistance which can be between
tvo vertexes vertex.*/
if (dl->mode == DIST2D_MAX)
if (dl->mode == DIST_MAX)
{
if (r>=0.5)
{

528
liblwgeom/measures3d.c Normal file
View file

@ -0,0 +1,528 @@
/**********************************************************************
* $Id: measures.c 5439 2010-03-16 03:13:33Z pramsey $
*
* PostGIS - Spatial Types for PostgreSQL
* http://postgis.refractions.net
* Copyright 2001-2006 Refractions Research Inc.
*
* This is free software; you can redistribute and/or modify it under
* the terms of the GNU General Public Licence. See the COPYING file.
*
**********************************************************************/
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "measures3d.h"
/*------------------------------------------------------------------------------------------------------------
Initializing functions
The functions starting the distance-calculation processses
--------------------------------------------------------------------------------------------------------------*/
/**
Function initializing 3dshortestline and 3dlongestline calculations.
*/
LWGEOM *
lw_dist3d_distanceline(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode)
{
double x1,x2,y1,y2, z1, z2;
double initdistance = ( mode == DIST_MIN ? MAXFLOAT : -1.0);
DISTPTS3D thedl;
LWPOINT *lwpoints[2];
LWGEOM *result;
thedl.mode = mode;
thedl.distance = initdistance;
thedl.tolerance = 0.0;
LWDEBUG(2, "lw_dist3d_distanceline is called");
if (!lw_dist3d_recursive( lw1,lw2,&thedl))
{
/*should never get here. all cases ought to be error handled earlier*/
lwerror("Some unspecified error.");
result =(LWGEOM *)lwcollection_construct_empty(srid, 0, 0);
}
/*if thedl.distance is unchanged there where only empty geometries input*/
if (thedl.distance == initdistance)
{
LWDEBUG(3, "didn't find geometries to measure between, returning null");
result =(LWGEOM *)lwcollection_construct_empty(srid, 0, 0);
}
else
{
x1=thedl.p1.x;
y1=thedl.p1.y;
z1=thedl.p1.z;
x2=thedl.p2.x;
y2=thedl.p2.y;
z2=thedl.p2.z;
lwpoints[0] = make_lwpoint3dz(srid, x1, y1, z1);
lwpoints[1] = make_lwpoint3dz(srid, x2, y2, z2);
result = (LWGEOM *)lwline_from_lwpointarray(srid, 2, lwpoints);
}
return result;
}
/**
Function initializing 3dclosestpoint calculations.
*/
LWGEOM *
lw_dist3d_distancepoint(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode)
{
double x,y,z;
DISTPTS3D thedl;
double initdistance = MAXFLOAT;
LWGEOM *result;
thedl.mode = mode;
thedl.distance= initdistance;
thedl.tolerance = 0;
LWDEBUG(2, "lw_dist3d_distancepoint is called");
if (!lw_dist3d_recursive( lw1,lw2,&thedl))
{
/*should never get here. all cases ought to be error handled earlier*/
lwerror("Some unspecified error.");
result =(LWGEOM *)lwcollection_construct_empty(srid, 0, 0);
}
if (thedl.distance == initdistance)
{
LWDEBUG(3, "didn't find geometries to measure between, returning null");
result =(LWGEOM *)lwcollection_construct_empty(srid, 0, 0);
}
else
{
x=thedl.p1.x;
y=thedl.p1.y;
z=thedl.p1.z;
result = (LWGEOM *)make_lwpoint3dz(srid, x, y, z);
}
return result;
}
/**
Function initialazing 3d max distance calculation
*/
double
lwgeom_maxdistance3d(LWGEOM *lw1, LWGEOM *lw2)
{
LWDEBUG(2, "lwgeom_maxdistance3d is called");
return lwgeom_maxdistance3d_tolerance( lw1, lw2, 0.0 );
}
/**
Function handling 3d max distance calculations and dfyllywithin calculations.
The difference is just the tolerance.
*/
double
lwgeom_maxdistance3d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
{
/*double thedist;*/
DISTPTS3D thedl;
LWDEBUG(2, "lwgeom_maxdistance3d_tolerance is called");
thedl.mode = DIST_MAX;
thedl.distance= -1;
thedl.tolerance = tolerance;
if (lw_dist3d_recursive( lw1,lw2,&thedl))
{
return thedl.distance;
}
/*should never get here. all cases ought to be error handled earlier*/
lwerror("Some unspecified error.");
return -1;
}
/**
Function initialazing 3d min distance calculation
*/
double
lwgeom_mindistance3d(LWGEOM *lw1, LWGEOM *lw2)
{
LWDEBUG(2, "lwgeom_mindistance3d is called");
return lwgeom_mindistance3d_tolerance( lw1, lw2, 0.0 );
}
/**
Function handling 3d min distance calculations and dwithin calculations.
The difference is just the tolerance.
*/
double
lwgeom_mindistance3d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
{
DISTPTS3D thedl;
LWDEBUG(2, "lwgeom_mindistance3d_tolerance is called");
thedl.mode = DIST_MIN;
thedl.distance= MAXFLOAT;
thedl.tolerance = tolerance;
if (lw_dist3d_recursive( lw1,lw2,&thedl))
{
return thedl.distance;
}
/*should never get here. all cases ought to be error handled earlier*/
lwerror("Some unspecified error.");
return MAXFLOAT;
}
/*------------------------------------------------------------------------------------------------------------
End of Initializing functions
--------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------------------------------------------
Preprocessing functions
Functions preparing geometries for distance-calculations
--------------------------------------------------------------------------------------------------------------*/
/**
This is a recursive function delivering every possible combinatin of subgeometries
*/
int lw_dist3d_recursive(const LWCOLLECTION * lwg1,const LWCOLLECTION * lwg2,DISTPTS3D *dl)
{
int i, j;
int n1=1;
int n2=1;
LWGEOM *g1, *g2;
LWDEBUGF(2, "lw_dist3d_recursive is called with type1=%d, type2=%d", TYPE_GETTYPE(lwg1->type), TYPE_GETTYPE(lwg2->type));
if (lwgeom_is_collection(TYPE_GETTYPE(lwg1->type)))
{
LWDEBUG(3, "First geometry is collection");
n1=lwg1->ngeoms;
}
if (lwgeom_is_collection(TYPE_GETTYPE(lwg2->type)))
{
LWDEBUG(3, "Second geometry is collection");
n2=lwg2->ngeoms;
}
for ( i = 0; i < n1; i++ )
{
if (lwgeom_is_collection(TYPE_GETTYPE(lwg1->type)))
{
g1=lwg1->geoms[i];
}
else
{
g1=(LWGEOM*)lwg1;
}
if (lwgeom_is_empty(g1)) return LW_TRUE;
if (lwgeom_is_collection(TYPE_GETTYPE(g1->type)))
{
LWDEBUG(3, "Found collection inside first geometry collection, recursing");
if (!lw_dist3d_recursive((LWCOLLECTION*)g1, (LWCOLLECTION*)lwg2, dl)) return LW_FALSE;
continue;
}
for ( j = 0; j < n2; j++ )
{
if (lwgeom_is_collection(TYPE_GETTYPE(lwg2->type)))
{
g2=lwg2->geoms[j];
}
else
{
g2=(LWGEOM*)lwg2;
}
if (lwgeom_is_collection(TYPE_GETTYPE(g2->type)))
{
LWDEBUG(3, "Found collection inside second geometry collection, recursing");
if (!lw_dist3d_recursive((LWCOLLECTION*) g1, (LWCOLLECTION*)g2, dl)) return LW_FALSE;
continue;
}
/*If one of geometries is empty, return. True here only means continue searching. False would have stoped the process*/
if (lwgeom_is_empty(g1)||lwgeom_is_empty(g2)) return LW_TRUE;
if (!lw_dist3d_distribute_bruteforce(g1, g2, dl)) return LW_FALSE;
if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
}
}
return LW_TRUE;
}
/**
This function distributes the "old-type" brut-force for 3D so far the only type, tasks depending on type
*/
int
lw_dist3d_distribute_bruteforce(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS3D *dl)
{
int t1 = TYPE_GETTYPE(lwg1->type);
int t2 = TYPE_GETTYPE(lwg2->type);
LWDEBUGF(2, "lw_dist3d_distribute_bruteforce is called with typ1=%d, type2=%d", TYPE_GETTYPE(lwg1->type), TYPE_GETTYPE(lwg2->type));
if ( t1 == POINTTYPE )
{
if ( t2 == POINTTYPE )
{
dl->twisted=1;
return lw_dist3d_point_point((LWPOINT *)lwg1, (LWPOINT *)lwg2, dl);
}
else if ( t2 == LINETYPE )
{
dl->twisted=1;
return lw_dist3d_point_line((LWPOINT *)lwg1, (LWLINE *)lwg2, dl);
}
else if ( t2 == POLYGONTYPE )
{
lwerror("Polygons are not yet supported for 3d distance calculations");
}
else
{
lwerror("Unsupported geometry type: %s", lwtype_name(t2));
return LW_FALSE;
}
}
else if ( t1 == LINETYPE )
{
if ( t2 == POINTTYPE )
{
dl->twisted=(-1);
return lw_dist3d_point_line((LWPOINT *)lwg2,(LWLINE *)lwg1,dl);
}
else if ( t2 == LINETYPE )
{
lwerror("Linestring to Linestring distance measures is not yet supported for 3d");
}
else if ( t2 == POLYGONTYPE )
{
lwerror("Polygons are not yet supported for 3d distance calculations");
}
else
{
lwerror("Unsupported geometry type: %s", lwtype_name(t2));
return LW_FALSE;
}
}
else if ( t1 == POLYGONTYPE )
{
if ( t2 == POLYGONTYPE )
{
lwerror("Polygons are not yet supported for 3d distance calculations");
}
else if ( t2 == POINTTYPE )
{
lwerror("Polygons are not yet supported for 3d distance calculations");
}
else if ( t2 == LINETYPE )
{
lwerror("Polygons are not yet supported for 3d distance calculations");
}
else
{
lwerror("Unsupported geometry type: %s", lwtype_name(t2));
return LW_FALSE;
}
}
else
{
lwerror("Unsupported geometry type: %s", lwtype_name(t1));
return LW_FALSE;
}
/*You shouldn't being able to get here*/
lwerror("unspecified error in function lw_dist3d_distribute_bruteforce");
return LW_FALSE;
}
/*------------------------------------------------------------------------------------------------------------
End of Preprocessing functions
--------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------------------------------------------
Brute force functions
So far the only way to do 3D-calculations
--------------------------------------------------------------------------------------------------------------*/
/**
point to point calculation
*/
int
lw_dist3d_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS3D *dl)
{
POINT3DZ p1;
POINT3DZ p2;
getPoint3dz_p(point1->point, 0, &p1);
getPoint3dz_p(point2->point, 0, &p2);
return lw_dist3d_pt_pt(&p1, &p2,dl);
}
/**
point to line calculation
*/
int
lw_dist3d_point_line(LWPOINT *point, LWLINE *line, DISTPTS3D *dl)
{
POINT3DZ p;
POINTARRAY *pa = line->points;
LWDEBUG(2, "lw_dist3d_point_line is called");
getPoint3dz_p(point->point, 0, &p);
return lw_dist3d_pt_ptarray(&p, pa, dl);
}
/**
* search all the segments of pointarray to see which one is closest to p1
* Returns minimum distance between point and pointarray
*/
int
lw_dist3d_pt_ptarray(POINT3DZ *p, POINTARRAY *pa,DISTPTS3D *dl)
{
int t;
POINT3DZ start, end;
int twist = dl->twisted;
LWDEBUG(2, "lw_dist3d_pt_ptarray is called");
getPoint3dz_p(pa, 0, &start);
for (t=1; t<pa->npoints; t++)
{
dl->twisted=twist;
getPoint3dz_p(pa, t, &end);
if (!lw_dist3d_pt_seg(p, &start, &end,dl)) return LW_FALSE;
if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
start = end;
}
return LW_TRUE;
}
/*------------------------------------------------------------------------------------------------------------
End of Brute force functions
--------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------------------------------------------
Functions in common for Brute force and new calculation
--------------------------------------------------------------------------------------------------------------*/
/**
This one is sending every occation to lw_dist3d_pt_pt
*/
int
lw_dist3d_pt_seg(POINT3DZ *p, POINT3DZ *A, POINT3DZ *B, DISTPTS3D *dl)
{
POINT3DZ c;
double r;
/*if start==end, then use pt distance */
if ( ( A->x == B->x) && (A->y == B->y) && (A->z == B->z) )
{
return lw_dist3d_pt_pt(p,A,dl);
}
r = ( (p->x-A->x) * (B->x-A->x) + (p->y-A->y) * (B->y-A->y) + ( p->z-A->z) * (B->z-A->z) )/( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y)+(B->z-A->z)*(B->z-A->z) );
/*This is for finding the 3Dmaxdistance.
the maxdistance have to be between two vertexes,
compared to mindistance which can be between
tvo vertexes vertex.*/
if (dl->mode == DIST_MAX)
{
if (r>=0.5)
{
return lw_dist3d_pt_pt(p,A,dl);
}
if (r<0.5)
{
return lw_dist3d_pt_pt(p,B,dl);
}
}
if (r<0) /*If the first vertex A is closest to the point p*/
{
return lw_dist3d_pt_pt(p,A,dl);
}
if (r>1) /*If the second vertex B is closest to the point p*/
{
return lw_dist3d_pt_pt(p,B,dl);
}
/*else if the point p is closer to some point between a and b
then we find that point and send it to lw_dist3d_pt_pt*/
c.x=A->x + r * (B->x-A->x);
c.y=A->y + r * (B->y-A->y);
c.z=A->z + r * (B->z-A->z);
return lw_dist3d_pt_pt(p,&c,dl);
}
/**
Compares incomming points and
stores the points closest to each other
or most far away from each other
depending on dl->mode (max or min)
*/
int
lw_dist3d_pt_pt(POINT3DZ *thep1, POINT3DZ *thep2,DISTPTS3D *dl)
{
double dx = thep2->x - thep1->x;
double dy = thep2->y - thep1->y;
double dz = thep2->z - thep1->z;
double dist = sqrt ( dx*dx + dy*dy + dz*dz);
if (((dl->distance - dist)*(dl->mode))>0) /*multiplication with mode to handle mindistance (mode=1) and maxdistance (mode = (-1)*/
{
dl->distance = dist;
if (dl->twisted>0) /*To get the points in right order. twisted is updated between 1 and (-1) every time the order is changed earlier in the chain*/
{
dl->p1 = *thep1;
dl->p2 = *thep2;
}
else
{
dl->p1 = *thep2;
dl->p2 = *thep1;
}
}
return LW_TRUE;
}
/*------------------------------------------------------------------------------------------------------------
End of Functions in common for Brute force and new calculation
--------------------------------------------------------------------------------------------------------------*/

52
liblwgeom/measures3d.h Normal file
View file

@ -0,0 +1,52 @@
/**********************************************************************
* $Id: measures.h 4715 2009-11-01 17:58:42Z nicklas $
*
* PostGIS - Spatial Types for PostgreSQL
* http://postgis.refractions.net
* Copyright 2001-2006 Refractions Research Inc.
* Copyright 2007-2008 Mark Cave-Ayland
* Copyright 2008 Paul Ramsey <pramsey@cleverelephant.ca>
*
* This is free software; you can redistribute and/or modify it under
* the terms of the GNU General Public Licence. See the COPYING file.
*
**********************************************************************/
#include "liblwgeom.h"
/**
Structure used in distance-calculations
*/
typedef struct
{
double distance; /*the distance between p1 and p2*/
POINT3DZ p1;
POINT3DZ p2;
int mode; /*the direction of looking, if thedir = -1 then we look for 3dmaxdistance and if it is 1 then we look for 3dmindistance*/
int twisted; /*To preserve the order of incoming points to match the first and second point in 3dshortest and 3dlongest line*/
double tolerance; /*the tolerance for 3ddwithin and 3ddfullywithin*/
} DISTPTS3D;
/*
Preprocessing functions
*/
int lw_dist3d_distribute_bruteforce(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS3D *dl);
int lw_dist3d_recursive(const LWCOLLECTION * lwg1,const LWCOLLECTION * lwg2, DISTPTS3D *dl);
int lw_dist3d_distribute_fast(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS3D *dl);
/*
Brute force functions
*/
int lw_dist3d_pt_ptarray(POINT3DZ *p, POINTARRAY *pa, DISTPTS3D *dl);
int lw_dist3d_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS3D *dl);
int lw_dist3d_point_line(LWPOINT *point, LWLINE *line, DISTPTS3D *dl);
/*
Functions in common for Brute force and new calculation
*/
int lw_dist3d_pt_pt(POINT3DZ *p1, POINT3DZ *p2, DISTPTS3D *dl);
int lw_dist3d_pt_seg(POINT3DZ *p, POINT3DZ *A, POINT3DZ *B, DISTPTS3D *dl);

267
postgis/lwgeom_functions_basic.c
View file

@ -45,11 +45,23 @@ Datum LWGEOM_length2d_linestring(PG_FUNCTION_ARGS);
Datum LWGEOM_length_linestring(PG_FUNCTION_ARGS);
Datum LWGEOM_perimeter2d_poly(PG_FUNCTION_ARGS);
Datum LWGEOM_perimeter_poly(PG_FUNCTION_ARGS);
Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS);
Datum LWGEOM_mindistance2d(PG_FUNCTION_ARGS);
Datum LWGEOM_closestpoint(PG_FUNCTION_ARGS);
Datum LWGEOM_shortestline2d(PG_FUNCTION_ARGS);
Datum LWGEOM_longestline2d(PG_FUNCTION_ARGS);
Datum LWGEOM_maxdistance3d(PG_FUNCTION_ARGS);
Datum LWGEOM_mindistance3d(PG_FUNCTION_ARGS);
Datum LWGEOM_closestpoint3d(PG_FUNCTION_ARGS);
Datum LWGEOM_shortestline3d(PG_FUNCTION_ARGS);
Datum LWGEOM_longestline3d(PG_FUNCTION_ARGS);
Datum LWGEOM_dwithin3d(PG_FUNCTION_ARGS);
Datum LWGEOM_dfullywithin3d(PG_FUNCTION_ARGS);
Datum LWGEOM_inside_circle_point(PG_FUNCTION_ARGS);
Datum LWGEOM_collect(PG_FUNCTION_ARGS);
Datum LWGEOM_accum(PG_FUNCTION_ARGS);
@ -1704,7 +1716,7 @@ Datum LWGEOM_closestpoint(PG_FUNCTION_ARGS)
srid = geom1->SRID;
point = lw_dist2d_distancepoint(geom1, geom2, srid, DIST2D_MIN);
point = lw_dist2d_distancepoint(geom1, geom2, srid, DIST_MIN);
if (lwgeom_is_empty(point))
{
PG_RETURN_NULL();
@ -1734,7 +1746,7 @@ Datum LWGEOM_shortestline2d(PG_FUNCTION_ARGS)
srid = geom1->SRID;
theline = lw_dist2d_distanceline(geom1, geom2, srid, DIST2D_MIN);
theline = lw_dist2d_distanceline(geom1, geom2, srid, DIST_MIN);
if (lwgeom_is_empty(theline))
{
PG_RETURN_NULL();
@ -1764,7 +1776,7 @@ Datum LWGEOM_longestline2d(PG_FUNCTION_ARGS)
srid = geom1->SRID;
theline = lw_dist2d_distanceline(geom1, geom2, srid, DIST2D_MAX);
theline = lw_dist2d_distanceline(geom1, geom2, srid, DIST_MAX);
if (lwgeom_is_empty(theline))
{
PG_RETURN_NULL();
@ -1929,6 +1941,255 @@ Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS)
PG_RETURN_NULL();
}
/**
Returns the point in first input geometry that is closest to the second input geometry in 3D
*/
PG_FUNCTION_INFO_V1(LWGEOM_closestpoint3d);
Datum LWGEOM_closestpoint3d(PG_FUNCTION_ARGS)
{
int srid;
LWGEOM *geom1;
LWGEOM *geom2;
LWGEOM *point;
geom1 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0))));
geom2 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1))));
if (geom1->SRID != geom2->SRID)
{
elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n");
PG_RETURN_NULL();
}
srid = geom1->SRID;
point = lw_dist3d_distancepoint(geom1, geom2, srid, DIST_MIN);
if (lwgeom_is_empty(point))
{
PG_RETURN_NULL();
}
PG_RETURN_POINTER(pglwgeom_serialize(point));
}
/**
Returns the shortest line between two geometries in 3D
*/
PG_FUNCTION_INFO_V1(LWGEOM_shortestline3d);
Datum LWGEOM_shortestline3d(PG_FUNCTION_ARGS)
{
int srid;
LWGEOM *geom1;
LWGEOM *geom2;
LWGEOM *theline;
geom1 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0))));
geom2 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1))));
if (geom1->SRID != geom2->SRID)
{
elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n");
PG_RETURN_NULL();
}
srid = geom1->SRID;
theline = lw_dist3d_distanceline(geom1, geom2, srid, DIST_MIN);
if (lwgeom_is_empty(theline))
{
PG_RETURN_NULL();
}
PG_RETURN_POINTER(pglwgeom_serialize(theline));
}
/**
Returns the longest line between two geometries in 3D
*/
PG_FUNCTION_INFO_V1(LWGEOM_longestline3d);
Datum LWGEOM_longestline3d(PG_FUNCTION_ARGS)
{
int srid;
LWGEOM *geom1;
LWGEOM *geom2;
LWGEOM *theline;
geom1 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0))));
geom2 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1))));
if (geom1->SRID != geom2->SRID)
{
elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n");
PG_RETURN_NULL();
}
srid = geom1->SRID;
theline = lw_dist3d_distanceline(geom1, geom2, srid, DIST_MAX);
if (lwgeom_is_empty(theline))
{
PG_RETURN_NULL();
}
PG_RETURN_POINTER(pglwgeom_serialize(theline));
}
/**
Minimum 2d distance between objects in geom1 and geom2 in 3D
*/
PG_FUNCTION_INFO_V1(LWGEOM_mindistance3d);
Datum LWGEOM_mindistance3d(PG_FUNCTION_ARGS)
{
LWGEOM *geom1;
LWGEOM *geom2;
double mindist;
PROFSTART(PROF_QRUN);
geom1 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0))));
geom2 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1))));
if (geom1->SRID != geom2->SRID)
{
elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n");
PG_RETURN_NULL();
}
mindist = lwgeom_mindistance3d(geom1, geom2);
PROFSTOP(PROF_QRUN);
PROFREPORT("dist",geom1, geom2, NULL);
PG_FREE_IF_COPY(geom1, 0);
PG_FREE_IF_COPY(geom2, 1);
/*if called with empty geometries the ingoing mindistance is untouched, and makes us return NULL*/
if (mindist<MAXFLOAT)
{
PG_RETURN_FLOAT8(mindist);
}
PG_RETURN_NULL();
}
/**
Returns boolean describing if
mininimum 3d distance between objects in
geom1 and geom2 is shorter than tolerance
*/
PG_FUNCTION_INFO_V1(LWGEOM_dwithin3d);
Datum LWGEOM_dwithin3d(PG_FUNCTION_ARGS)
{
LWGEOM *geom1;
LWGEOM *geom2;
double mindist, tolerance;
PROFSTART(PROF_QRUN);
geom1 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0))));
geom2 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1))));
tolerance = PG_GETARG_FLOAT8(2);
if ( tolerance < 0 )
{
elog(ERROR,"Tolerance cannot be less than zero\n");
PG_RETURN_NULL();
}
if (geom1->SRID != geom2->SRID)
{
elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n");
PG_RETURN_NULL();
}
mindist = lwgeom_mindistance3d_tolerance(geom1,geom2,tolerance);
PROFSTOP(PROF_QRUN);
PROFREPORT("dist",geom1, geom2, NULL);
PG_FREE_IF_COPY(geom1, 0);
PG_FREE_IF_COPY(geom2, 1);
/*empty geometries cases should be right handled since return from underlying
functions should be MAXFLOAT which causes false as answer*/
PG_RETURN_BOOL(tolerance >= mindist);
}
/**
Returns boolean describing if
maximum 3d distance between objects in
geom1 and geom2 is shorter than tolerance
*/
PG_FUNCTION_INFO_V1(LWGEOM_dfullywithin3d);
Datum LWGEOM_dfullywithin3d(PG_FUNCTION_ARGS)
{
LWGEOM *geom1;
LWGEOM *geom2;
double maxdist, tolerance;
PROFSTART(PROF_QRUN);
geom1 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0))));
geom2 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1))));
tolerance = PG_GETARG_FLOAT8(2);
if ( tolerance < 0 )
{
elog(ERROR,"Tolerance cannot be less than zero\n");
PG_RETURN_NULL();
}
if (geom1->SRID != geom2->SRID)
{
elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n");
PG_RETURN_NULL();
}
maxdist = lwgeom_maxdistance3d_tolerance(geom1, geom2, tolerance);
PROFSTOP(PROF_QRUN);
PROFREPORT("dist",geom1, geom2, NULL);
PG_FREE_IF_COPY(geom1, 0);
PG_FREE_IF_COPY(geom2, 1);
/*If function is feed with empty geometries we should return false*/
if (maxdist>-1)
{
PG_RETURN_BOOL(tolerance >= maxdist);
}
PG_RETURN_BOOL(LW_FALSE);
}
/**
Maximum 3d distance between objects in geom1 and geom2.
*/
PG_FUNCTION_INFO_V1(LWGEOM_maxdistance3d);
Datum LWGEOM_maxdistance3d(PG_FUNCTION_ARGS)
{
LWGEOM *geom1;
LWGEOM *geom2;
double maxdist;
PROFSTART(PROF_QRUN);
geom1 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0))));
geom2 = lwgeom_deserialize(SERIALIZED_FORM((PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1))));
if (geom1->SRID != geom2->SRID)
{
elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n");
PG_RETURN_NULL();
}
maxdist = lwgeom_maxdistance3d(geom1, geom2);
PROFSTOP(PROF_QRUN);
PROFREPORT("maxdist",geom1, geom2, NULL);
PG_FREE_IF_COPY(geom1, 0);
PG_FREE_IF_COPY(geom2, 1);
/*if called with empty geometries the ingoing mindistance is untouched, and makes us return NULL*/
if (maxdist>-1)
{
PG_RETURN_FLOAT8(maxdist);
}
PG_RETURN_NULL();
}
PG_FUNCTION_INFO_V1(LWGEOM_longitude_shift);
Datum LWGEOM_longitude_shift(PG_FUNCTION_ARGS)

66
postgis/postgis.sql.in.c
View file

@ -6313,6 +6313,72 @@ LANGUAGE 'plpgsql' IMMUTABLE STRICT;
#include "sqlmm.sql.in.c"
#include "geography.sql.in.c"
---------------------------------------------------------------
-- 3D-functions
---------------------------------------------------------------
CREATE OR REPLACE FUNCTION ST_3DDistance(geometry,geometry)
RETURNS float8
AS 'MODULE_PATHNAME', 'LWGEOM_mindistance3d'
LANGUAGE 'C' IMMUTABLE STRICT
COST 100;
CREATE OR REPLACE FUNCTION ST_3DMaxDistance(geometry,geometry)
RETURNS float8
AS 'MODULE_PATHNAME', 'LWGEOM_maxdistance3d'
LANGUAGE 'C' IMMUTABLE STRICT
COST 100;
CREATE OR REPLACE FUNCTION ST_3DClosestPoint(geometry,geometry)
RETURNS geometry
AS 'MODULE_PATHNAME', 'LWGEOM_closestpoint3d'
LANGUAGE 'C' IMMUTABLE STRICT
COST 100;
CREATE OR REPLACE FUNCTION ST_3DShortestLine(geometry,geometry)
RETURNS geometry
AS 'MODULE_PATHNAME', 'LWGEOM_shortestline3d'
LANGUAGE 'C' IMMUTABLE STRICT
COST 100;
CREATE OR REPLACE FUNCTION ST_3DLongestLine(geometry,geometry)
RETURNS geometry
AS 'MODULE_PATHNAME', 'LWGEOM_longestline3d'
LANGUAGE 'C' IMMUTABLE STRICT
COST 100;
CREATE OR REPLACE FUNCTION _ST_3DDWithin(geometry,geometry,float8)
RETURNS boolean
AS 'MODULE_PATHNAME', 'LWGEOM_dwithin3d'
LANGUAGE 'C' IMMUTABLE STRICT
COST 100;
CREATE OR REPLACE FUNCTION ST_3DDWithin(geometry,geometry,float8)
RETURNS boolean
AS 'SELECT $1 && ST_Expand($2,$3) AND $2 && ST_Expand($1,$3) AND _ST_3DDWithin($1, $2, $3)'
LANGUAGE 'SQL' IMMUTABLE
COST 100;
CREATE OR REPLACE FUNCTION _ST_3DDFullyWithin(geometry,geometry,float8)
RETURNS boolean
AS 'MODULE_PATHNAME', 'LWGEOM_dfullywithin3d'
LANGUAGE 'C' IMMUTABLE STRICT
COST 100;
CREATE OR REPLACE FUNCTION ST_3DDFullyWithin(geometry,geometry,float8)
RETURNS boolean
AS 'SELECT $1 && ST_Expand($2,$3) AND $2 && ST_Expand($1,$3) AND _ST_3DDFullyWithin($1, $2, $3)'
LANGUAGE 'SQL' IMMUTABLE
COST 100;
CREATE OR REPLACE FUNCTION ST_3DIntersects(geometry,geometry)
RETURNS boolean
AS 'SELECT $1 && $2 AND _ST_3DDWithin($1, $2, 0.0)'
LANGUAGE 'SQL' IMMUTABLE
COST 100;
---------------------------------------------------------------
-- SQL-MM
---------------------------------------------------------------

51
regress/measures.sql
View file

@ -201,8 +201,59 @@ select 'distancepoly6',
geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b
) as foo;
--3D Distance functions
SELECT '3dDistancetest1',
ST_3DDistance(a,b),
ST_3DMaxDistance(a,b),
ST_3DDWithin(a,b,5),
ST_3DDFullyWithin(a,b,5),
ST_ASEWKT(ST_3DShortestline(a,b)),
ST_ASEWKT(ST_3DClosestpoint(a,b)),
ST_ASEWKT(ST_3DLongestline(a,b)) FROM (
SELECT 'POINT(1 1 1)'::geometry as a, 'POINT(3 2 7)'::geometry as b
) as foo;
SELECT '3dDistancetest2',
ST_3DDistance(a,b),
ST_3DMaxDistance(a,b),
ST_3DDWithin(a,b,5),
ST_3DDFullyWithin(a,b,5),
ST_ASEWKT(ST_3DShortestline(a,b)),
ST_ASEWKT(ST_3DClosestpoint(a,b)),
ST_ASEWKT(ST_3DLongestline(a,b)) FROM (
SELECT 'POINT(1 1 1)'::geometry as a, 'LINESTRING(0 0 0, 2 2 2)'::geometry as b
) as foo;
SELECT '3dDistancetest3',
ST_3DDistance(a,b),
ST_3DMaxDistance(a,b),
ST_3DDWithin(a,b,5),
ST_3DDFullyWithin(a,b,5),
ST_ASEWKT(ST_3DShortestline(a,b)),
ST_ASEWKT(ST_3DClosestpoint(a,b)),
ST_ASEWKT(ST_3DLongestline(a,b)) FROM (
SELECT 'POINT(1 1 1)'::geometry as a, 'LINESTRING(5 2 6, -3 -2 4)'::geometry as b
) as foo;
SELECT '3dDistancetest4',
ST_3DDistance(a,b),
ST_3DMaxDistance(a,b),
ST_3DDWithin(a,b,5),
ST_3DDFullyWithin(a,b,5),
ST_ASEWKT(ST_3DShortestline(a,b)),
ST_ASEWKT(ST_3DClosestpoint(a,b)),
ST_ASEWKT(ST_3DLongestline(a,b)) FROM (
SELECT 'LINESTRING(1 1 3, 5 7 8)'::geometry as a, 'POINT(1 1 1)'::geometry as b
) as foo;
SELECT 'unsupported_test1',
ST_3DDistance(a,b) FROM (
SELECT 'LINESTRING(1 1 1 , 2 2 2)'::geometry as a, 'LINESTRING(0 0 0, 2 2 2)'::geometry as b) as foo;
-- Area of an empty polygon
select 'emptyPolyArea', st_area('POLYGON EMPTY');

5
regress/measures_expected
View file

@ -38,6 +38,11 @@ distancepoly3|0|26.9072480941474|LINESTRING(17 19,17 19)|LINESTRING(17 19,17 19)
distancepoly4|0|28.3196045170126|LINESTRING(16 19,16 19)|LINESTRING(16 19,16 19)|LINESTRING(18 20,-1 -1)|LINESTRING(-1 -1,18 20)
distancepoly5|0|26.1725046566048|LINESTRING(17 12,17 12)|LINESTRING(17 12,17 12)|LINESTRING(17 18,-1 -1)|LINESTRING(-1 -1,17 18)
distancepoly6|0|32.5269119345812|LINESTRING(2 2,2 2)|LINESTRING(2 2,2 2)|LINESTRING(2 2,25 25)|LINESTRING(25 25,2 2)
3dDistancetest1|6.40312423743285|6.40312423743285|f|f|LINESTRING(1 1 1,3 2 7)|POINT(1 1 1)|LINESTRING(1 1 1,3 2 7)
3dDistancetest2|0|1.73205080756888|t|t|LINESTRING(1 1 1,1 1 1)|POINT(1 1 1)|LINESTRING(1 1 1,0 0 0)
3dDistancetest3|4.09994192757944|6.48074069840786|t|f|LINESTRING(1 1 1,0.619047619047619 -0.190476190476191 4.90476190476191)|POINT(1 1 1)|LINESTRING(1 1 1,5 2 6)
3dDistancetest4|2|10.0498756211209|t|f|LINESTRING(1 1 3,1 1 1)|POINT(1 1 3)|LINESTRING(5 7 8,1 1 1)
ERROR: Linestring to Linestring distance measures is not yet supported for 3d
emptyPolyArea|0
emptyLineArea|0
emptyPointArea|0