PostGIS  2.2.7dev-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 1731 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().

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

Here is the call graph for this function:

Here is the caller graph for this function: