PostGIS  2.3.8dev-r@@SVN_REVISION@@

◆ ptarray_distance_spheroid()

static double ptarray_distance_spheroid ( const POINTARRAY pa1,
const POINTARRAY pa2,
const SPHEROID s,
double  tolerance,
int  check_intersection 
)
static

Definition at line 1742 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().

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