PostGIS  2.5.0dev-r@@SVN_REVISION@@
static double ptarray_distance_spheroid ( const POINTARRAY pa1,
const POINTARRAY pa2,
const SPHEROID s,
double  tolerance,
int  check_intersection 
)
static

Definition at line 1744 of file lwgeodetic.c.

References SPHEROID::a, SPHEROID::b, distance(), edge_distance_to_edge(), edge_distance_to_point(), edge_intersects(), GEOGRAPHIC_EDGE::end, geog2cart(), geographic_point_init(), getPoint2d_cp(), GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, LW_ON_INTERRUPT, LWDEBUG, LWDEBUGF, POINTARRAY::npoints, SPHEROID::radius, sphere_distance(), spheroid_distance(), GEOGRAPHIC_EDGE::start, POINT2D::x, and POINT2D::y.

Referenced by lwgeom_distance_spheroid().

1745 {
1746  GEOGRAPHIC_EDGE e1, e2;
1747  GEOGRAPHIC_POINT g1, g2;
1748  GEOGRAPHIC_POINT nearest1, nearest2;
1749  POINT3D A1, A2, B1, B2;
1750  const POINT2D *p;
1751  double distance;
1752  uint32_t i, j;
1753  int use_sphere = (s->a == s->b ? 1 : 0);
1754 
1755  /* Make result really big, so that everything will be smaller than it */
1756  distance = FLT_MAX;
1757 
1758  /* Empty point arrays? Return negative */
1759  if ( pa1->npoints == 0 || pa2->npoints == 0 )
1760  return -1.0;
1761 
1762  /* Handle point/point case here */
1763  if ( pa1->npoints == 1 && pa2->npoints == 1 )
1764  {
1765  p = getPoint2d_cp(pa1, 0);
1766  geographic_point_init(p->x, p->y, &g1);
1767  p = getPoint2d_cp(pa2, 0);
1768  geographic_point_init(p->x, p->y, &g2);
1769  /* Sphere special case, axes equal */
1770  distance = s->radius * sphere_distance(&g1, &g2);
1771  if ( use_sphere )
1772  return distance;
1773  /* Below tolerance, actual distance isn't of interest */
1774  else if ( distance < 0.95 * tolerance )
1775  return distance;
1776  /* Close or greater than tolerance, get the real answer to be sure */
1777  else
1778  return spheroid_distance(&g1, &g2, s);
1779  }
1780 
1781  /* Handle point/line case here */
1782  if ( pa1->npoints == 1 || pa2->npoints == 1 )
1783  {
1784  /* Handle one/many case here */
1785  uint32_t i;
1786  const POINTARRAY *pa_one;
1787  const POINTARRAY *pa_many;
1788 
1789  if ( pa1->npoints == 1 )
1790  {
1791  pa_one = pa1;
1792  pa_many = pa2;
1793  }
1794  else
1795  {
1796  pa_one = pa2;
1797  pa_many = pa1;
1798  }
1799 
1800  /* Initialize our point */
1801  p = getPoint2d_cp(pa_one, 0);
1802  geographic_point_init(p->x, p->y, &g1);
1803 
1804  /* Initialize start of line */
1805  p = getPoint2d_cp(pa_many, 0);
1806  geographic_point_init(p->x, p->y, &(e1.start));
1807 
1808  /* Iterate through the edges in our line */
1809  for ( i = 1; i < pa_many->npoints; i++ )
1810  {
1811  double d;
1812  p = getPoint2d_cp(pa_many, i);
1813  geographic_point_init(p->x, p->y, &(e1.end));
1814  /* Get the spherical distance between point and edge */
1815  d = s->radius * edge_distance_to_point(&e1, &g1, &g2);
1816  /* New shortest distance! Record this distance / location */
1817  if ( d < distance )
1818  {
1819  distance = d;
1820  nearest2 = g2;
1821  }
1822  /* We've gotten closer than the tolerance... */
1823  if ( d < tolerance )
1824  {
1825  /* Working on a sphere? The answer is correct, return */
1826  if ( use_sphere )
1827  {
1828  return d;
1829  }
1830  /* Far enough past the tolerance that the spheroid calculation won't change things */
1831  else if ( d < tolerance * 0.95 )
1832  {
1833  return d;
1834  }
1835  /* On a spheroid and near the tolerance? Confirm that we are *actually* closer than tolerance */
1836  else
1837  {
1838  d = spheroid_distance(&g1, &nearest2, s);
1839  /* Yes, closer than tolerance, return! */
1840  if ( d < tolerance )
1841  return d;
1842  }
1843  }
1844  e1.start = e1.end;
1845  }
1846 
1847  /* On sphere, return answer */
1848  if ( use_sphere )
1849  return distance;
1850  /* On spheroid, calculate final answer based on closest approach */
1851  else
1852  return spheroid_distance(&g1, &nearest2, s);
1853 
1854  }
1855 
1856  /* Initialize start of line 1 */
1857  p = getPoint2d_cp(pa1, 0);
1858  geographic_point_init(p->x, p->y, &(e1.start));
1859  geog2cart(&(e1.start), &A1);
1860 
1861 
1862  /* Handle line/line case */
1863  for ( i = 1; i < pa1->npoints; i++ )
1864  {
1865  p = getPoint2d_cp(pa1, i);
1866  geographic_point_init(p->x, p->y, &(e1.end));
1867  geog2cart(&(e1.end), &A2);
1868 
1869  /* Initialize start of line 2 */
1870  p = getPoint2d_cp(pa2, 0);
1871  geographic_point_init(p->x, p->y, &(e2.start));
1872  geog2cart(&(e2.start), &B1);
1873 
1874  for ( j = 1; j < pa2->npoints; j++ )
1875  {
1876  double d;
1877 
1878  p = getPoint2d_cp(pa2, j);
1879  geographic_point_init(p->x, p->y, &(e2.end));
1880  geog2cart(&(e2.end), &B2);
1881 
1882  LWDEBUGF(4, "e1.start == GPOINT(%.6g %.6g) ", e1.start.lat, e1.start.lon);
1883  LWDEBUGF(4, "e1.end == GPOINT(%.6g %.6g) ", e1.end.lat, e1.end.lon);
1884  LWDEBUGF(4, "e2.start == GPOINT(%.6g %.6g) ", e2.start.lat, e2.start.lon);
1885  LWDEBUGF(4, "e2.end == GPOINT(%.6g %.6g) ", e2.end.lat, e2.end.lon);
1886 
1887  if ( check_intersection && edge_intersects(&A1, &A2, &B1, &B2) )
1888  {
1889  LWDEBUG(4,"edge intersection! returning 0.0");
1890  return 0.0;
1891  }
1892  d = s->radius * edge_distance_to_edge(&e1, &e2, &g1, &g2);
1893  LWDEBUGF(4,"got edge_distance_to_edge %.8g", d);
1894 
1895  if ( d < distance )
1896  {
1897  distance = d;
1898  nearest1 = g1;
1899  nearest2 = g2;
1900  }
1901  if ( d < tolerance )
1902  {
1903  if ( use_sphere )
1904  {
1905  return d;
1906  }
1907  else
1908  {
1909  d = spheroid_distance(&nearest1, &nearest2, s);
1910  if ( d < tolerance )
1911  return d;
1912  }
1913  }
1914 
1915  /* Copy end to start to allow a new end value in next iteration */
1916  e2.start = e2.end;
1917  B1 = B2;
1918  }
1919 
1920  /* Copy end to start to allow a new end value in next iteration */
1921  e1.start = e1.end;
1922  A1 = A2;
1923  LW_ON_INTERRUPT(return -1.0);
1924  }
1925  LWDEBUGF(4,"finished all loops, returning %.8g", distance);
1926 
1927  if ( use_sphere )
1928  return distance;
1929  else
1930  return spheroid_distance(&nearest1, &nearest2, s);
1931 }
double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Given two points on a unit sphere, calculate their distance apart in radians.
Definition: lwgeodetic.c:913
double edge_distance_to_edge(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2)
Calculate the distance between two edges.
Definition: lwgeodetic.c:1234
Two-point great circle segment from a to b.
Definition: lwgeodetic.h:56
double b
Definition: liblwgeom.h:313
double radius
Definition: liblwgeom.h:317
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
#define LW_ON_INTERRUPT(x)
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:47
unsigned int uint32_t
Definition: uthash.h:78
double x
Definition: liblwgeom.h:327
GEOGRAPHIC_POINT start
Definition: lwgeodetic.h:58
double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest)
Definition: lwgeodetic.c:1183
uint32_t edge_intersects(const POINT3D *A1, const POINT3D *A2, const POINT3D *B1, const POINT3D *B2)
Returns non-zero if edges A and B interact.
Definition: lwgeodetic.c:3411
GEOGRAPHIC_POINT end
Definition: lwgeodetic.h:59
double y
Definition: liblwgeom.h:327
Datum distance(PG_FUNCTION_ARGS)
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesion coordinates on unit sphere.
Definition: lwgeodetic.c:369
double a
Definition: liblwgeom.h:312
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:171
double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
Computes the shortest distance along the surface of the spheroid between two points.
Definition: lwspheroid.c:186
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:364
uint32_t npoints
Definition: liblwgeom.h:370

Here is the call graph for this function:

Here is the caller graph for this function: