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

Definition at line 741 of file gserialized_gist_nd.c.

References distance(), gidx_distance_leaf_centroid(), 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.

742 {
743  char b1mem[GIDX_MAX_SIZE];
744  GIDX *b1 = (GIDX*)b1mem;
745  char b2mem[GIDX_MAX_SIZE];
746  GIDX *b2 = (GIDX*)b2mem;
747 
748 #if POSTGIS_PGSQL_VERSION < 95
749 
750  /* Centroid-to-centroid distance */
751  Datum gs1 = PG_GETARG_DATUM(0);
752  Datum gs2 = PG_GETARG_DATUM(1);
753  double box_distance = FLT_MAX;
754 
755  /* Must be able to build box for each argument (ie, not empty geometry). */
756  if ( (gserialized_datum_get_gidx_p(gs1, b1) == LW_SUCCESS) &&
757  (gserialized_datum_get_gidx_p(gs2, b2) == LW_SUCCESS) )
758  {
759  box_distance = gidx_distance_leaf_centroid(b1, b2);
760  POSTGIS_DEBUGF(3, "got boxes %s and %s", gidx_to_string(b1), gidx_to_string(b2));
761  }
762  PG_RETURN_FLOAT8(box_distance);
763 
764 #else /* POSTGIS_PGSQL_VERSION >= 96 */
765 
766  /* Feature-to-feature distance */
767  GSERIALIZED *geom1 = PG_GETARG_GSERIALIZED_P(0);
768  GSERIALIZED *geom2 = PG_GETARG_GSERIALIZED_P(1);
769  LWGEOM *lw1 = lwgeom_from_gserialized(geom1);
770  LWGEOM *lw2 = lwgeom_from_gserialized(geom2);
771  LWGEOM *closest;
772  double distance;
773 
774 
775  /* Find an exact shortest line w/ the dimensions we support */
776  if ( lwgeom_has_z(lw1) && lwgeom_has_z(lw2) )
777  {
778  closest = lwgeom_closest_line_3d(lw1, lw2);
779  distance = lwgeom_length(closest);
780  }
781  else
782  {
783  closest = lwgeom_closest_line(lw1, lw2);
784  distance = lwgeom_length_2d(closest);
785  }
786 
787  /* Un-sqrt the distance so we can add extra terms */
788  distance = distance*distance;
789 
790  /* Can only add the M term if both objects have M */
791  if ( lwgeom_has_m(lw1) && lwgeom_has_m(lw2) )
792  {
793  double m1, m2;
794  int usebox = false;
795 
796  if ( lwgeom_get_type(lw1) == POINTTYPE )
797  {
798  POINT4D p;
799  lwpoint_getPoint4d_p((LWPOINT*)lw1, &p);
800  m1 = p.m;
801  }
802  else if ( lwgeom_get_type(lw1) == LINETYPE )
803  {
804  LWPOINT *lwp1 = lwline_get_lwpoint(lwgeom_as_lwline(closest), 0);
805  m1 = lwgeom_interpolate_point(lw1, lwp1);
806  lwpoint_free(lwp1);
807  }
808  else
809  {
810  usebox = true;
811  }
812 
813  if ( lwgeom_get_type(lw2) == POINTTYPE )
814  {
815  POINT4D p;
816  lwpoint_getPoint4d_p((LWPOINT*)lw2, &p);
817  m2 = p.m;
818  }
819  else if ( lwgeom_get_type(lw2) == LINETYPE )
820  {
821  LWPOINT *lwp2 = lwline_get_lwpoint(lwgeom_as_lwline(closest), 1);
822  m2 = lwgeom_interpolate_point(lw2, lwp2);
823  lwpoint_free(lwp2);
824  }
825  else
826  {
827  usebox = true;
828  }
829 
830  if ( usebox )
831  {
832  double d;
833  gserialized_get_gidx_p(geom1, b1);
834  gserialized_get_gidx_p(geom2, b2);
835  d = gidx_distance_m(b1, b2);
836  distance += d*d;
837  }
838  else
839  {
840  distance += (m2-m1)*(m2-m1);
841  }
842  }
843 
844  lwgeom_free(closest);
845 
846  PG_FREE_IF_COPY(geom1, 0);
847  PG_FREE_IF_COPY(geom2, 1);
848  PG_RETURN_FLOAT8(sqrt(distance));
849 #endif /* POSTGIS_PGSQL_VERSION >= 96 */
850 }
#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:842
void lwpoint_free(LWPOINT *pt)
Definition: lwpoint.c:195
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1063
#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:1681
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:849
double lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt)
Find the measure value at the location on the line closest to the point.
static double gidx_distance_leaf_centroid(const GIDX *a, const GIDX *b)
Calculate the centroid->centroid distance between the boxes.
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:102
double lwgeom_length(const LWGEOM *geom)
Definition: lwgeom.c:1659
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
LWPOINT * lwline_get_lwpoint(const LWLINE *line, int where)
Returns freshly allocated LWPOINT that corresponds to the index where.
Definition: lwline.c:324
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:856

Here is the call graph for this function: