PostGIS  2.5.0beta2dev-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 1748 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().

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