PostGIS  2.1.10dev-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 1730 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_p(), GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, LWDEBUG, LWDEBUGF, MAXFLOAT, POINTARRAY::npoints, SPHEROID::radius, sphere_distance(), spheroid_distance(), GEOGRAPHIC_EDGE::start, POINT2D::x, and POINT2D::y.

Referenced by lwgeom_distance_spheroid().

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