PostGIS  2.5.0dev-r@@SVN_REVISION@@
Datum gserialized_distance_nd ( PG_FUNCTION_ARGS  )

Definition at line 807 of file gserialized_gist_nd.c.

References distance(), gidx_distance_m(), LINETYPE, LW_SUCCESS, lwgeom_as_lwline(), lwgeom_closest_line(), lwgeom_closest_line_3d(), lwgeom_free(), lwgeom_from_gserialized(), lwgeom_get_type(), lwgeom_has_m(), lwgeom_has_z(), lwgeom_interpolate_point(), lwgeom_length(), lwgeom_length_2d(), lwline_get_lwpoint(), lwpoint_free(), lwpoint_getPoint4d_p(), POINT4D::m, and POINTTYPE.

808 {
809  char b1mem[GIDX_MAX_SIZE];
810  GIDX *b1 = (GIDX*)b1mem;
811  char b2mem[GIDX_MAX_SIZE];
812  GIDX *b2 = (GIDX*)b2mem;
813 
814 #if POSTGIS_PGSQL_VERSION < 95
815 
816  /* Centroid-to-centroid distance */
817  Datum gs1 = PG_GETARG_DATUM(0);
818  Datum gs2 = PG_GETARG_DATUM(1);
819  double box_distance = FLT_MAX;
820 
821  /* Must be able to build box for each argument (ie, not empty geometry). */
822  if ( (gserialized_datum_get_gidx_p(gs1, b1) == LW_SUCCESS) &&
823  (gserialized_datum_get_gidx_p(gs2, b2) == LW_SUCCESS) )
824  {
825  box_distance = gidx_distance_leaf_centroid(b1, b2);
826  POSTGIS_DEBUGF(3, "got boxes %s and %s", gidx_to_string(b1), gidx_to_string(b2));
827  }
828  PG_RETURN_FLOAT8(box_distance);
829 
830 #else /* POSTGIS_PGSQL_VERSION >= 96 */
831 
832  /* Feature-to-feature distance */
833  GSERIALIZED *geom1 = PG_GETARG_GSERIALIZED_P(0);
834  GSERIALIZED *geom2 = PG_GETARG_GSERIALIZED_P(1);
835  LWGEOM *lw1 = lwgeom_from_gserialized(geom1);
836  LWGEOM *lw2 = lwgeom_from_gserialized(geom2);
837  LWGEOM *closest;
838  double distance;
839 
840 
841  /* Find an exact shortest line w/ the dimensions we support */
842  if ( lwgeom_has_z(lw1) && lwgeom_has_z(lw2) )
843  {
844  closest = lwgeom_closest_line_3d(lw1, lw2);
845  distance = lwgeom_length(closest);
846  }
847  else
848  {
849  closest = lwgeom_closest_line(lw1, lw2);
850  distance = lwgeom_length_2d(closest);
851  }
852 
853  /* Un-sqrt the distance so we can add extra terms */
854  distance = distance*distance;
855 
856  /* Can only add the M term if both objects have M */
857  if ( lwgeom_has_m(lw1) && lwgeom_has_m(lw2) )
858  {
859  double m1, m2;
860  int usebox = false;
861 
862  if ( lwgeom_get_type(lw1) == POINTTYPE )
863  {
864  POINT4D p;
865  lwpoint_getPoint4d_p((LWPOINT*)lw1, &p);
866  m1 = p.m;
867  }
868  else if ( lwgeom_get_type(lw1) == LINETYPE )
869  {
870  LWPOINT *lwp1 = lwline_get_lwpoint(lwgeom_as_lwline(closest), 0);
871  m1 = lwgeom_interpolate_point(lw1, lwp1);
872  lwpoint_free(lwp1);
873  }
874  else
875  {
876  usebox = true;
877  }
878 
879  if ( lwgeom_get_type(lw2) == POINTTYPE )
880  {
881  POINT4D p;
882  lwpoint_getPoint4d_p((LWPOINT*)lw2, &p);
883  m2 = p.m;
884  }
885  else if ( lwgeom_get_type(lw2) == LINETYPE )
886  {
887  LWPOINT *lwp2 = lwline_get_lwpoint(lwgeom_as_lwline(closest), 1);
888  m2 = lwgeom_interpolate_point(lw2, lwp2);
889  lwpoint_free(lwp2);
890  }
891  else
892  {
893  usebox = true;
894  }
895 
896  if ( usebox )
897  {
898  double d;
899  gserialized_get_gidx_p(geom1, b1);
900  gserialized_get_gidx_p(geom2, b2);
901  d = gidx_distance_m(b1, b2);
902  distance += d*d;
903  }
904  else
905  {
906  distance += (m2-m1)*(m2-m1);
907  }
908  }
909 
910  lwgeom_free(closest);
911 
912  PG_FREE_IF_COPY(geom1, 0);
913  PG_FREE_IF_COPY(geom2, 1);
914  PG_RETURN_FLOAT8(sqrt(distance));
915 #endif /* POSTGIS_PGSQL_VERSION >= 96 */
916 }
#define LINETYPE
Definition: liblwgeom.h:85
double m
Definition: liblwgeom.h:351
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwgeom.c:916
void lwpoint_free(LWPOINT *pt)
Definition: lwpoint.c:213
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1137
#define LW_SUCCESS
Definition: liblwgeom.h:79
LWGEOM * lwgeom_closest_line_3d(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures3d.c:88
double lwgeom_length_2d(const LWGEOM *geom)
Definition: lwgeom.c:1939
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:923
double lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt)
Find the measure value at the location on the line closest to the point.
int lwpoint_getPoint4d_p(const LWPOINT *point, POINT4D *out)
Definition: lwpoint.c:57
Datum distance(PG_FUNCTION_ARGS)
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:169
double lwgeom_length(const LWGEOM *geom)
Definition: lwgeom.c:1917
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
LWGEOM * lwgeom_closest_line(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures.c:42
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:930
static double gidx_distance_m(const GIDX *a, const GIDX *b)
LWPOINT * lwline_get_lwpoint(const LWLINE *line, uint32_t where)
Returns freshly allocated LWPOINT that corresponds to the index where.
Definition: lwline.c:318

Here is the call graph for this function: