/src/postgis/liblwgeom/measures.c

Source
/**********************************************************************
 *
 * PostGIS - Spatial Types for PostgreSQL
 * http://postgis.net
 *
 * PostGIS is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 2 of the License, or
 * (at your option) any later version.
 *
 * PostGIS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
 *
 **********************************************************************
 *
 * Copyright 2001-2006 Refractions Research Inc.
 * Copyright 2010 Nicklas Avén
 * Copyright 2012 Paul Ramsey
 * Copyright 2019 Darafei Praliaskouski
 *
 **********************************************************************/

#include <string.h>
#include <stdlib.h>

#include "../postgis_config.h"
//#define POSTGIS_DEBUG_LEVEL 5

#include "measures.h"
#include "lwgeom_log.h"

/*------------------------------------------------------------------------------------------------------------
Initializing functions
The functions starting the distance-calculation processes
--------------------------------------------------------------------------------------------------------------*/

LWGEOM *
lwgeom_closest_line(const LWGEOM *lw1, const LWGEOM *lw2)
{
  return lw_dist2d_distanceline(lw1, lw2, lw1->srid, DIST_MIN);
}

LWGEOM *
lwgeom_furthest_line(const LWGEOM *lw1, const LWGEOM *lw2)
{
  return lw_dist2d_distanceline(lw1, lw2, lw1->srid, DIST_MAX);
}

LWGEOM *
lwgeom_closest_point(const LWGEOM *lw1, const LWGEOM *lw2)
{
  return lw_dist2d_distancepoint(lw1, lw2, lw1->srid, DIST_MIN);
}

LWGEOM *
lwgeom_furthest_point(const LWGEOM *lw1, const LWGEOM *lw2)
{
  return lw_dist2d_distancepoint(lw1, lw2, lw1->srid, DIST_MAX);
}

void
lw_dist2d_distpts_init(DISTPTS *dl, int mode)
{
  dl->twisted = -1;
  dl->p1.x = dl->p1.y = 0.0;
  dl->p2.x = dl->p2.y = 0.0;
  dl->mode = mode;
  dl->tolerance = 0.0;
  if (mode == DIST_MIN)
    dl->distance = FLT_MAX;
  else
    dl->distance = -1 * FLT_MAX;
}

static void
lw_dist2d_distpts_set(DISTPTS *dl, double distance, const POINT2D *p1, const POINT2D *p2)
{
  int update = (dl->mode == DIST_MIN) ? (distance < dl->distance) : (distance > dl->distance);
  if (update)
  {
    dl->distance = distance;
    dl->p1 = *p1;
    dl->p2 = *p2;
  }
}

/**
Function initializing shortestline and longestline calculations.
*/
LWGEOM *
lw_dist2d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int32_t srid, int mode)
{
  double x1, x2, y1, y2;

  double initdistance = (mode == DIST_MIN ? FLT_MAX : -1.0);
  DISTPTS thedl;
  LWPOINT *lwpoints[2];
  LWGEOM *result;

  thedl.mode = mode;
  thedl.distance = initdistance;
  thedl.tolerance = 0.0;

  LWDEBUG(2, "lw_dist2d_distanceline is called");

  if (!lw_dist2d_comp(lw1, lw2, &thedl))
  {
    /*should never get here. all cases ought to be error handled earlier*/
    lwerror("Some unspecified error.");
    result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
  }

  /*if thedl.distance is unchanged there where only empty geometries input*/
  if (thedl.distance == initdistance)
  {
    LWDEBUG(3, "didn't find geometries to measure between, returning null");
    result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
  }
  else
  {
    x1 = thedl.p1.x;
    y1 = thedl.p1.y;
    x2 = thedl.p2.x;
    y2 = thedl.p2.y;

    lwpoints[0] = lwpoint_make2d(srid, x1, y1);
    lwpoints[1] = lwpoint_make2d(srid, x2, y2);

    result = (LWGEOM *)lwline_from_ptarray(srid, 2, lwpoints);
  }
  return result;
}

/**
Function initializing closestpoint calculations.
*/
LWGEOM *
lw_dist2d_distancepoint(const LWGEOM *lw1, const LWGEOM *lw2, int32_t srid, int mode)
{
  double x, y;
  DISTPTS thedl;
  double initdistance = FLT_MAX;
  LWGEOM *result;

  thedl.mode = mode;
  thedl.distance = initdistance;
  thedl.tolerance = 0;

  LWDEBUG(2, "lw_dist2d_distancepoint is called");

  if (!lw_dist2d_comp(lw1, lw2, &thedl))
  {
    /*should never get here. all cases ought to be error handled earlier*/
    lwerror("Some unspecified error.");
    result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
  }
  if (thedl.distance == initdistance)
  {
    LWDEBUG(3, "didn't find geometries to measure between, returning null");
    result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
  }
  else
  {
    x = thedl.p1.x;
    y = thedl.p1.y;
    result = (LWGEOM *)lwpoint_make2d(srid, x, y);
  }
  return result;
}

/**
Function initializing max distance calculation
*/
double
lwgeom_maxdistance2d(const LWGEOM *lw1, const LWGEOM *lw2)
{
  LWDEBUG(2, "lwgeom_maxdistance2d is called");

  return lwgeom_maxdistance2d_tolerance(lw1, lw2, 0.0);
}

/**
Function handling max distance calculations and dfullywithin calculations.
The difference is just the tolerance.
*/
double
lwgeom_maxdistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
{
  /*double thedist;*/
  DISTPTS thedl;
  LWDEBUG(2, "lwgeom_maxdistance2d_tolerance is called");
  thedl.mode = DIST_MAX;
  thedl.distance = -1;
  thedl.tolerance = tolerance;
  if (lw_dist2d_comp(lw1, lw2, &thedl))
    return thedl.distance;

  /*should never get here. all cases ought to be error handled earlier*/
  lwerror("Some unspecified error.");
  return -1;
}

/**
  Function initializing min distance calculation
*/
double
lwgeom_mindistance2d(const LWGEOM *lw1, const LWGEOM *lw2)
{
  return lwgeom_mindistance2d_tolerance(lw1, lw2, 0.0);
}

/**
  Function handling min distance calculations and dwithin calculations.
  The difference is just the tolerance.
*/
double
lwgeom_mindistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
{
  DISTPTS thedl;
  LWDEBUG(2, "lwgeom_mindistance2d_tolerance is called");
  thedl.mode = DIST_MIN;
  thedl.distance = FLT_MAX;
  thedl.tolerance = tolerance;
  if (lw_dist2d_comp(lw1, lw2, &thedl))
    return thedl.distance;
  /*should never get here. all cases ought to be error handled earlier*/
  lwerror("Some unspecified error.");
  return FLT_MAX;
}

/*------------------------------------------------------------------------------------------------------------
End of Initializing functions
--------------------------------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------------------------------------------
Preprocessing functions
Functions preparing geometries for distance-calculations
--------------------------------------------------------------------------------------------------------------*/

/**
  This function just deserializes geometries
  Bboxes is not checked here since it is the subgeometries
  bboxes we will use anyway.
*/
int
lw_dist2d_comp(const LWGEOM *lw1, const LWGEOM *lw2, DISTPTS *dl)
{
  return lw_dist2d_recursive(lw1, lw2, dl);
}

static int
lw_dist2d_is_collection(const LWGEOM *g)
{
  /* Differs from lwgeom_is_collection by not treating CURVEPOLYGON as collection */
  switch (g->type)
  {
  case TINTYPE:
  case MULTIPOINTTYPE:
  case MULTILINETYPE:
  case MULTIPOLYGONTYPE:
  case COLLECTIONTYPE:
  case MULTICURVETYPE:
  case MULTISURFACETYPE:
  case COMPOUNDTYPE:
  case POLYHEDRALSURFACETYPE:
    return LW_TRUE;
    break;

  default:
    return LW_FALSE;
  }
}


/* Check for overlapping bboxes */
static int
lw_dist2d_check_overlap(const LWGEOM *lwg1, const LWGEOM *lwg2)
{
  assert(lwg1 && lwg2 && lwg1->bbox && lwg2->bbox);

  /* Check if the geometries intersect. */
  if ((lwg1->bbox->xmax < lwg2->bbox->xmin || lwg1->bbox->xmin > lwg2->bbox->xmax ||
       lwg1->bbox->ymax < lwg2->bbox->ymin || lwg1->bbox->ymin > lwg2->bbox->ymax))
  {
    return LW_FALSE;
  }
  return LW_TRUE;
}

/**
This is a recursive function delivering every possible combination of subgeometries
*/
int
lw_dist2d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl)
{
  int i, j;
  int n1 = 1;
  int n2 = 1;
  LWGEOM *g1 = NULL;
  LWGEOM *g2 = NULL;
  LWCOLLECTION *c1 = NULL;
  LWCOLLECTION *c2 = NULL;

  LWDEBUGF(2, "lw_dist2d_comp is called with type1=%d, type2=%d", lwg1->type, lwg2->type);

  if (lw_dist2d_is_collection(lwg1))
  {
    LWDEBUG(3, "First geometry is collection");
    c1 = lwgeom_as_lwcollection(lwg1);
    n1 = c1->ngeoms;
  }
  if (lw_dist2d_is_collection(lwg2))
  {
    LWDEBUG(3, "Second geometry is collection");
    c2 = lwgeom_as_lwcollection(lwg2);
    n2 = c2->ngeoms;
  }

  for (i = 0; i < n1; i++)
  {

    if (lw_dist2d_is_collection(lwg1))
      g1 = c1->geoms[i];
    else
      g1 = (LWGEOM *)lwg1;

    if (!g1) continue;

    if (lwgeom_is_empty(g1))
      continue;

    if (lw_dist2d_is_collection(g1))
    {
      LWDEBUG(3, "Found collection inside first geometry collection, recursing");
      if (!lw_dist2d_recursive(g1, lwg2, dl))
        return LW_FALSE;
      continue;
    }
    for (j = 0; j < n2; j++)
    {
      if (lw_dist2d_is_collection(lwg2))
        g2 = c2->geoms[j];
      else
        g2 = (LWGEOM *)lwg2;

      if (!g2) continue;

      if (lw_dist2d_is_collection(g2))
      {
        LWDEBUG(3, "Found collection inside second geometry collection, recursing");
        if (!lw_dist2d_recursive(g1, g2, dl))
          return LW_FALSE;
        continue;
      }

      if (!g1->bbox)
        lwgeom_add_bbox(g1);

      if (!g2->bbox)
        lwgeom_add_bbox(g2);

      /* If one of geometries is empty, skip */
      if (lwgeom_is_empty(g1) || lwgeom_is_empty(g2))
        continue;

      if ((dl->mode != DIST_MAX) && (!lw_dist2d_check_overlap(g1, g2)) &&
          (g1->type == LINETYPE || g1->type == POLYGONTYPE || g1->type == TRIANGLETYPE) &&
          (g2->type == LINETYPE || g2->type == POLYGONTYPE || g2->type == TRIANGLETYPE))
      {
        if (!lw_dist2d_distribute_fast(g1, g2, dl))
          return LW_FALSE;
      }
      else
      {
        if (!lw_dist2d_distribute_bruteforce(g1, g2, dl))
          return LW_FALSE;
        LWDEBUGF(2, "Distance so far: %.15g (%.15g tolerated)", dl->distance, dl->tolerance);
        if (dl->distance <= dl->tolerance && dl->mode == DIST_MIN)
          return LW_TRUE; /*just a check if the answer is already given*/
        LWDEBUG(2, "Not below tolerance yet");
      }
    }
  }
  return LW_TRUE;
}

int
lw_dist2d_distribute_bruteforce(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS *dl)
{

  int t1 = lwg1->type;
  int t2 = lwg2->type;

  switch (t1)
  {
  case POINTTYPE:
  {
    dl->twisted = 1;
    switch (t2)
    {
    case POINTTYPE:
      return lw_dist2d_point_point((LWPOINT *)lwg1, (LWPOINT *)lwg2, dl);
    case LINETYPE:
      return lw_dist2d_point_line((LWPOINT *)lwg1, (LWLINE *)lwg2, dl);
    case TRIANGLETYPE:
      return lw_dist2d_point_tri((LWPOINT *)lwg1, (LWTRIANGLE *)lwg2, dl);
    case POLYGONTYPE:
      return lw_dist2d_point_poly((LWPOINT *)lwg1, (LWPOLY *)lwg2, dl);
    case CIRCSTRINGTYPE:
      return lw_dist2d_point_circstring((LWPOINT *)lwg1, (LWCIRCSTRING *)lwg2, dl);
    case CURVEPOLYTYPE:
      return lw_dist2d_point_curvepoly((LWPOINT *)lwg1, (LWCURVEPOLY *)lwg2, dl);
    default:
      lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(t2));
      return LW_FALSE;
    }
  }
  case LINETYPE:
  {
    dl->twisted = 1;
    switch (t2)
    {
    case POINTTYPE:
      dl->twisted = -1;
      return lw_dist2d_point_line((LWPOINT *)lwg2, (LWLINE *)lwg1, dl);
    case LINETYPE:
      return lw_dist2d_line_line((LWLINE *)lwg1, (LWLINE *)lwg2, dl);
    case TRIANGLETYPE:
      return lw_dist2d_line_tri((LWLINE *)lwg1, (LWTRIANGLE *)lwg2, dl);
    case POLYGONTYPE:
      return lw_dist2d_line_poly((LWLINE *)lwg1, (LWPOLY *)lwg2, dl);
    case CIRCSTRINGTYPE:
      return lw_dist2d_line_circstring((LWLINE *)lwg1, (LWCIRCSTRING *)lwg2, dl);
    case CURVEPOLYTYPE:
      return lw_dist2d_line_curvepoly((LWLINE *)lwg1, (LWCURVEPOLY *)lwg2, dl);
    default:
      lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(t2));
      return LW_FALSE;
    }
  }
  case TRIANGLETYPE:
  {
    dl->twisted = 1;
    switch (t2)
    {
    case POINTTYPE:
      dl->twisted = -1;
      return lw_dist2d_point_tri((LWPOINT *)lwg2, (LWTRIANGLE *)lwg1, dl);
    case LINETYPE:
      dl->twisted = -1;
      return lw_dist2d_line_tri((LWLINE *)lwg2, (LWTRIANGLE *)lwg1, dl);
    case TRIANGLETYPE:
      return lw_dist2d_tri_tri((LWTRIANGLE *)lwg1, (LWTRIANGLE *)lwg2, dl);
    case POLYGONTYPE:
      return lw_dist2d_tri_poly((LWTRIANGLE *)lwg1, (LWPOLY *)lwg2, dl);
    case CIRCSTRINGTYPE:
      return lw_dist2d_tri_circstring((LWTRIANGLE *)lwg1, (LWCIRCSTRING *)lwg2, dl);
    case CURVEPOLYTYPE:
      return lw_dist2d_tri_curvepoly((LWTRIANGLE *)lwg1, (LWCURVEPOLY *)lwg2, dl);
    default:
      lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(t2));
      return LW_FALSE;
    }
  }
  case CIRCSTRINGTYPE:
  {
    dl->twisted = 1;
    switch (t2)
    {
    case POINTTYPE:
      dl->twisted = -1;
      return lw_dist2d_point_circstring((LWPOINT *)lwg2, (LWCIRCSTRING *)lwg1, dl);
    case LINETYPE:
      dl->twisted = -1;
      return lw_dist2d_line_circstring((LWLINE *)lwg2, (LWCIRCSTRING *)lwg1, dl);
    case TRIANGLETYPE:
      dl->twisted = -1;
      return lw_dist2d_tri_circstring((LWTRIANGLE *)lwg2, (LWCIRCSTRING *)lwg1, dl);
    case POLYGONTYPE:
      return lw_dist2d_circstring_poly((LWCIRCSTRING *)lwg1, (LWPOLY *)lwg2, dl);
    case CIRCSTRINGTYPE:
      return lw_dist2d_circstring_circstring((LWCIRCSTRING *)lwg1, (LWCIRCSTRING *)lwg2, dl);
    case CURVEPOLYTYPE:
      return lw_dist2d_circstring_curvepoly((LWCIRCSTRING *)lwg1, (LWCURVEPOLY *)lwg2, dl);
    default:
      lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(t2));
      return LW_FALSE;
    }
  }
  case POLYGONTYPE:
  {
    dl->twisted = -1;
    switch (t2)
    {
    case POINTTYPE:
      return lw_dist2d_point_poly((LWPOINT *)lwg2, (LWPOLY *)lwg1, dl);
    case LINETYPE:
      return lw_dist2d_line_poly((LWLINE *)lwg2, (LWPOLY *)lwg1, dl);
    case TRIANGLETYPE:
      return lw_dist2d_tri_poly((LWTRIANGLE *)lwg2, (LWPOLY *)lwg1, dl);
    case CIRCSTRINGTYPE:
      return lw_dist2d_circstring_poly((LWCIRCSTRING *)lwg2, (LWPOLY *)lwg1, dl);
    case POLYGONTYPE:
      dl->twisted = 1;
      return lw_dist2d_poly_poly((LWPOLY *)lwg1, (LWPOLY *)lwg2, dl);
    case CURVEPOLYTYPE:
      dl->twisted = 1;
      return lw_dist2d_poly_curvepoly((LWPOLY *)lwg1, (LWCURVEPOLY *)lwg2, dl);
    default:
      lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(t2));
      return LW_FALSE;
    }
  }
  case CURVEPOLYTYPE:
  {
    dl->twisted = -1;
    switch (t2)
    {
    case POINTTYPE:
      return lw_dist2d_point_curvepoly((LWPOINT *)lwg2, (LWCURVEPOLY *)lwg1, dl);
    case LINETYPE:
      return lw_dist2d_line_curvepoly((LWLINE *)lwg2, (LWCURVEPOLY *)lwg1, dl);
    case TRIANGLETYPE:
      return lw_dist2d_tri_curvepoly((LWTRIANGLE *)lwg2, (LWCURVEPOLY *)lwg1, dl);
    case POLYGONTYPE:
      return lw_dist2d_poly_curvepoly((LWPOLY *)lwg2, (LWCURVEPOLY *)lwg1, dl);
    case CIRCSTRINGTYPE:
      return lw_dist2d_circstring_curvepoly((LWCIRCSTRING *)lwg2, (LWCURVEPOLY *)lwg1, dl);
    case CURVEPOLYTYPE:
      dl->twisted = 1;
      return lw_dist2d_curvepoly_curvepoly((LWCURVEPOLY *)lwg1, (LWCURVEPOLY *)lwg2, dl);
    default:
      lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(t2));
      return LW_FALSE;
    }
  }
  default:
  {
    lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(t1));
    return LW_FALSE;
  }
  }

  return LW_FALSE;
}


/** Geometries are distributed for the new faster distance-calculations */
int
lw_dist2d_distribute_fast(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS *dl)
{
  POINTARRAY *pa1, *pa2;
  int type1 = lwg1->type;
  int type2 = lwg2->type;

  LWDEBUGF(2, "lw_dist2d_distribute_fast is called with typ1=%d, type2=%d", lwg1->type, lwg2->type);

  switch (type1)
  {
  case LINETYPE:
    pa1 = ((LWLINE *)lwg1)->points;
    break;
  case POLYGONTYPE:
    pa1 = ((LWPOLY *)lwg1)->rings[0];
    break;
  case TRIANGLETYPE:
    pa1 = ((LWTRIANGLE *)lwg1)->points;
    break;
  default:
    lwerror("Unsupported geometry1 type: %s", lwtype_name(type1));
    return LW_FALSE;
  }
  switch (type2)
  {
  case LINETYPE:
    pa2 = ((LWLINE *)lwg2)->points;
    break;
  case POLYGONTYPE:
    pa2 = ((LWPOLY *)lwg2)->rings[0];
    break;
  case TRIANGLETYPE:
    pa2 = ((LWTRIANGLE *)lwg2)->points;
    break;
  default:
    lwerror("Unsupported geometry2 type: %s", lwtype_name(type1));
    return LW_FALSE;
  }
  dl->twisted = 1;
  return lw_dist2d_fast_ptarray_ptarray(pa1, pa2, dl, lwg1->bbox, lwg2->bbox);
}

/*------------------------------------------------------------------------------------------------------------
End of Preprocessing functions
--------------------------------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------------------------------------------
Brute force functions
The old way of calculating distances, now used for:
1)  distances to points (because there shouldn't be anything to gain by the new way of doing it)
2)  distances when subgeometries geometries bboxes overlaps
--------------------------------------------------------------------------------------------------------------*/

/**
point to point calculation
*/
int
lw_dist2d_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS *dl)
{
  const POINT2D *p1 = getPoint2d_cp(point1->point, 0);
  const POINT2D *p2 = getPoint2d_cp(point2->point, 0);
  return lw_dist2d_pt_pt(p1, p2, dl);
}
/**
point to line calculation
*/
int
lw_dist2d_point_line(LWPOINT *point, LWLINE *line, DISTPTS *dl)
{
  const POINT2D *p = getPoint2d_cp(point->point, 0);
  return lw_dist2d_pt_ptarray(p, line->points, dl);
}

int
lw_dist2d_point_tri(LWPOINT *point, LWTRIANGLE *tri, DISTPTS *dl)
{
  const POINT2D *pt = getPoint2d_cp(point->point, 0);
  /* Is point inside triangle? */
  if (dl->mode == DIST_MIN && ptarray_contains_point(tri->points, pt) != LW_OUTSIDE)
  {
    lw_dist2d_distpts_set(dl, 0.0, pt, pt);
    return LW_TRUE;
  }

  return lw_dist2d_pt_ptarray(pt, tri->points, dl);
}

int
lw_dist2d_point_circstring(LWPOINT *point, LWCIRCSTRING *circ, DISTPTS *dl)
{
  const POINT2D *p;
  p = getPoint2d_cp(point->point, 0);
  return lw_dist2d_pt_ptarrayarc(p, circ->points, dl);
}

/**
 * 1. see if pt in outer boundary. if no, then treat the outer ring like a line
 * 2. if in the boundary, test to see if its in a hole.
 *    if so, then return dist to hole, else return 0 (point in polygon)
 */
int
lw_dist2d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl)
{
  const POINT2D *p = getPoint2d_cp(point->point, 0);

  /* Max distance? Check only outer ring.*/
  if (dl->mode == DIST_MAX)
    return lw_dist2d_pt_ptarray(p, poly->rings[0], dl);

  /* Return distance to outer ring if not inside it */
  if (ptarray_contains_point(poly->rings[0], p) == LW_OUTSIDE)
    return lw_dist2d_pt_ptarray(p, poly->rings[0], dl);

  /*
   * Inside the outer ring.
   * Scan though each of the inner rings looking to
   * see if its inside.  If not, distance==0.
   * Otherwise, distance = pt to ring distance
   */
  for (uint32_t i = 1; i < poly->nrings; i++)
    if (ptarray_contains_point(poly->rings[i], p) != LW_OUTSIDE)
      return lw_dist2d_pt_ptarray(p, poly->rings[i], dl);

  /* Is inside the polygon */
  lw_dist2d_distpts_set(dl, 0.0, p, p);
  return LW_TRUE;
}

int
lw_dist2d_point_curvepoly(LWPOINT *point, LWCURVEPOLY *poly, DISTPTS *dl)
{
  const POINT2D *p = getPoint2d_cp(point->point, 0);

  if (dl->mode == DIST_MAX)
    lwerror("lw_dist2d_point_curvepoly cannot calculate max distance");

  /* Return distance to outer ring if not inside it */
  if (lwgeom_contains_point(poly->rings[0], p) == LW_OUTSIDE)
    return lw_dist2d_recursive((LWGEOM *)point, poly->rings[0], dl);

  /* Inside the outer ring.
   * Scan though each of the inner rings looking to see if its inside.  If not, distance==0.
   * Otherwise, distance = pt to ring distance.
   */
  for (uint32_t i = 1; i < poly->nrings; i++)
    if (lwgeom_contains_point(poly->rings[i], p) == LW_INSIDE)
      return lw_dist2d_recursive((LWGEOM *)point, poly->rings[i], dl);

  /* Is inside the polygon */
  lw_dist2d_distpts_set(dl, 0.0, p, p);
  return LW_TRUE;
}

int
lw_dist2d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl)
{
  POINTARRAY *pa1 = line1->points;
  POINTARRAY *pa2 = line2->points;
  return lw_dist2d_ptarray_ptarray(pa1, pa2, dl);
}

int
lw_dist2d_line_tri(LWLINE *line, LWTRIANGLE *tri, DISTPTS *dl)
{
  const POINT2D *pt = getPoint2d_cp(line->points, 0);
  /* Is there a point inside triangle? */
  if (dl->mode == DIST_MIN && ptarray_contains_point(tri->points, pt) != LW_OUTSIDE)
  {
    lw_dist2d_distpts_set(dl, 0.0, pt, pt);
    return LW_TRUE;
  }

  return lw_dist2d_ptarray_ptarray(line->points, tri->points, dl);
}

int
lw_dist2d_line_circstring(LWLINE *line1, LWCIRCSTRING *line2, DISTPTS *dl)
{
  return lw_dist2d_ptarray_ptarrayarc(line1->points, line2->points, dl);
}

/**
 * line to polygon calculation
 * Brute force.
 * Test line-ring distance against each ring.
 * If there's an intersection (distance==0) then return 0 (crosses boundary).
 * Otherwise, test to see if any point is inside outer rings of polygon,
 * but not in inner rings.
 * If so, return 0  (line inside polygon),
 * otherwise return min distance to a ring (could be outside
 * polygon or inside a hole)
 */
int
lw_dist2d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl)
{
  POINTARRAY *pa = line->points;
  const POINT2D *pt = getPoint2d_cp(pa, 0);

  /* Line has a point outside poly. Check distance to outer ring only. */
  if (ptarray_contains_point(poly->rings[0], pt) == LW_OUTSIDE || dl->mode == DIST_MAX)
    return lw_dist2d_ptarray_ptarray(pa, poly->rings[0], dl);

  for (uint32_t i = 1; i < poly->nrings; i++)
  {
    if (!lw_dist2d_ptarray_ptarray(pa, poly->rings[i], dl))
      return LW_FALSE;

    /* just a check if the answer is already given */
    if (dl->distance <= dl->tolerance && dl->mode == DIST_MIN)
      return LW_TRUE;
  }

  /* It's inside a hole, then the actual distance is the min ring distance */
  for (uint32_t i = 1; i < poly->nrings; i++)
    if (ptarray_contains_point(poly->rings[i], pt) != LW_OUTSIDE)
      return LW_TRUE;

  /* Not in hole, so inside polygon */
  if (dl->mode == DIST_MIN)
  {
    lw_dist2d_distpts_set(dl, 0.0, pt, pt);
  }
  return LW_TRUE;
}

int
lw_dist2d_line_curvepoly(LWLINE *line, LWCURVEPOLY *poly, DISTPTS *dl)
{
  const POINT2D *pt = getPoint2d_cp(line->points, 0);

  /* Line has a point outside curvepoly. Check distance to outer ring only. */
  if (lwgeom_contains_point(poly->rings[0], pt) == LW_OUTSIDE)
    return lw_dist2d_recursive((LWGEOM *)line, poly->rings[0], dl);

  for (uint32_t i = 1; i < poly->nrings; i++)
  {
    if (!lw_dist2d_recursive((LWGEOM *)line, poly->rings[i], dl))
      return LW_FALSE;

    /* just a check if the answer is already given */
    if (dl->distance <= dl->tolerance && dl->mode == DIST_MIN)
      return LW_TRUE;
  }

  /* It's inside a hole, then distance is actual min distance */
  for (uint32_t i = 1; i < poly->nrings; i++)
    if (lwgeom_contains_point(poly->rings[i], pt) != LW_OUTSIDE)
      return LW_TRUE;

  /* Not in hole, so inside polygon */
  if (dl->mode == DIST_MIN)
  {
    lw_dist2d_distpts_set(dl, 0.0, pt, pt);
  }
  return LW_TRUE;
}

int
lw_dist2d_tri_tri(LWTRIANGLE *tri1, LWTRIANGLE *tri2, DISTPTS *dl)
{
  POINTARRAY *pa1 = tri1->points;
  POINTARRAY *pa2 = tri2->points;
  const POINT2D *pt = getPoint2d_cp(pa2, 0);
  if (dl->mode == DIST_MIN && ptarray_contains_point(pa1, pt) != LW_OUTSIDE)
  {
    lw_dist2d_distpts_set(dl, 0.0, pt, pt);
    return LW_TRUE;
  }

  pt = getPoint2d_cp(pa1, 0);
  if (dl->mode == DIST_MIN && ptarray_contains_point(pa2, pt) != LW_OUTSIDE)
  {
    lw_dist2d_distpts_set(dl, 0.0, pt, pt);
    return LW_TRUE;
  }

  return lw_dist2d_ptarray_ptarray(pa1, pa2, dl);
}

int
lw_dist2d_tri_poly(LWTRIANGLE *tri, LWPOLY *poly, DISTPTS *dl)
{
  POINTARRAY *pa = tri->points;
  const POINT2D *pt = getPoint2d_cp(pa, 0);

  /* If we are looking for maxdistance, just check the outer rings.*/
  if (dl->mode == DIST_MAX)
    return lw_dist2d_ptarray_ptarray(pa, poly->rings[0], dl);

  /* Triangle has a point outside poly. Check distance to outer ring only. */
  if (ptarray_contains_point(poly->rings[0], pt) == LW_OUTSIDE)
  {
    if (!lw_dist2d_ptarray_ptarray(pa, poly->rings[0], dl))
      return LW_FALSE;

    /* just a check if the answer is already given */
    if (dl->distance <= dl->tolerance)
      return LW_TRUE;

    /* Maybe poly is inside triangle? */
    const POINT2D *pt2 = getPoint2d_cp(poly->rings[0], 0);
    if (ptarray_contains_point(pa, pt2) != LW_OUTSIDE)
    {
      lw_dist2d_distpts_set(dl, 0.0, pt2, pt2);
      return LW_TRUE;
    }
  }

  for (uint32_t i = 1; i < poly->nrings; i++)
  {
    if (!lw_dist2d_ptarray_ptarray(pa, poly->rings[i], dl))
      return LW_FALSE;

    /* just a check if the answer is already given */
    if (dl->distance <= dl->tolerance && dl->mode == DIST_MIN)
      return LW_TRUE;
  }

  /* It's inside a hole, then the actual distance is the min ring distance */
  for (uint32_t i = 1; i < poly->nrings; i++)
    if (ptarray_contains_point(poly->rings[i], pt) != LW_OUTSIDE)
      return LW_TRUE;

  /* Not in hole, so inside polygon */
  lw_dist2d_distpts_set(dl, 0.0, pt, pt);
  return LW_TRUE;
}
static const POINT2D *
lw_curvering_getfirstpoint2d_cp(LWGEOM *geom, POINT2D *scratch)
{
  switch (geom->type)
  {
  case LINETYPE:
    return getPoint2d_cp(((LWLINE *)geom)->points, 0);
  case CIRCSTRINGTYPE:
    return getPoint2d_cp(((LWCIRCSTRING *)geom)->points, 0);
  case COMPOUNDTYPE:
  {
    LWCOMPOUND *comp = (LWCOMPOUND *)geom;
    return lw_curvering_getfirstpoint2d_cp(comp->geoms[0], scratch);
  }
  case NURBSCURVETYPE:
  {
    LWPOINT *pt = lwnurbscurve_evaluate((LWNURBSCURVE *)geom, 0.0);
    if (!scratch || !pt || !lwpoint_getPoint2d_p(pt, scratch))
    {
      lwpoint_free(pt);
      lwerror("lw_curvering_getfirstpoint2d_cp: invalid NURBSCURVE ring");
    }
    lwpoint_free(pt);
    return scratch;
  }
  default:
    lwerror("lw_curvering_getfirstpoint2d_cp: unknown type");
  }
  return NULL;
}

int
lw_dist2d_tri_curvepoly(LWTRIANGLE *tri, LWCURVEPOLY *poly, DISTPTS *dl)
{
  const POINT2D *pt = getPoint2d_cp(tri->points, 0);
  POINT2D ringpt;

  /* If we are looking for maxdistance, just check the outer rings.*/
  if (dl->mode == DIST_MAX)
    return lw_dist2d_recursive((LWGEOM *)tri, poly->rings[0], dl);

  /* Line has a point outside curvepoly. Check distance to outer ring only. */
  if (lwgeom_contains_point(poly->rings[0], pt) == LW_OUTSIDE)
  {
    if (lw_dist2d_recursive((LWGEOM *)tri, poly->rings[0], dl))
      return LW_TRUE;
    /* Maybe poly is inside triangle? */
    if (lwgeom_contains_point((LWGEOM *)tri, lw_curvering_getfirstpoint2d_cp(poly->rings[0], &ringpt)) != LW_OUTSIDE)
    {
      lw_dist2d_distpts_set(dl, 0.0, pt, pt);
      return LW_TRUE;
    }
  }

  for (uint32_t i = 1; i < poly->nrings; i++)
  {
    if (!lw_dist2d_recursive((LWGEOM *)tri, poly->rings[i], dl))
      return LW_FALSE;

    /* just a check if the answer is already given */
    if (dl->distance <= dl->tolerance && dl->mode == DIST_MIN)
      return LW_TRUE;
  }

  /* It's inside a hole, then distance is actual min distance */
  for (uint32_t i = 1; i < poly->nrings; i++)
    if (lwgeom_contains_point(poly->rings[i], pt) != LW_OUTSIDE)
      return LW_TRUE;

  /* Not in hole, so inside polygon */
  lw_dist2d_distpts_set(dl, 0.0, pt, pt);
  return LW_TRUE;
}

int
lw_dist2d_tri_circstring(LWTRIANGLE *tri, LWCIRCSTRING *line, DISTPTS *dl)
{
  POINT2D scratch;
  const POINT2D *pt = lw_curvering_getfirstpoint2d_cp((LWGEOM *)line, &scratch);
  if (ptarray_contains_point(tri->points, pt) != LW_OUTSIDE && dl->mode == DIST_MIN)
  {
    lw_dist2d_distpts_set(dl, 0.0, pt, pt);
    return LW_TRUE;
  }

  return lw_dist2d_ptarray_ptarrayarc(tri->points, line->points, dl);
}

/**
Function handling polygon to polygon calculation
1 if we are looking for maxdistance, just check the outer rings.
2 check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings
3 check if first point of poly2 is in a hole of poly1. If so check outer ring of poly2 against that hole of poly1
4 check if first point of poly1 is in a hole of poly2. If so check outer ring of poly1 against that hole of poly2
5 If we have come all the way here we know that the first point of one of them is inside the other ones outer ring
and not in holes so we check which one is inside.
 */
int
lw_dist2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl)
{
  const POINT2D *pt;

  LWDEBUG(2, "lw_dist2d_poly_poly called");

  /*1 if we are looking for maxdistance, just check the outer rings.*/
  if (dl->mode == DIST_MAX)
    return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl);

  /* 2  check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings
  here it would be possible to handle the information about which one is inside which one and only search for the
  smaller ones in the bigger ones holes.*/
  pt = getPoint2d_cp(poly1->rings[0], 0);
  if (ptarray_contains_point(poly2->rings[0], pt) == LW_OUTSIDE)
  {
    pt = getPoint2d_cp(poly2->rings[0], 0);
    if (ptarray_contains_point(poly1->rings[0], pt) == LW_OUTSIDE)
      return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl);
  }

  /*3 check if first point of poly2 is in a hole of poly1. If so check outer ring of poly2 against that hole
   * of poly1*/
  pt = getPoint2d_cp(poly2->rings[0], 0);
  for (uint32_t i = 1; i < poly1->nrings; i++)
    if (ptarray_contains_point(poly1->rings[i], pt) != LW_OUTSIDE)
      return lw_dist2d_ptarray_ptarray(poly1->rings[i], poly2->rings[0], dl);

  /*4 check if first point of poly1 is in a hole of poly2. If so check outer ring of poly1 against that hole
   * of poly2*/
  pt = getPoint2d_cp(poly1->rings[0], 0);
  for (uint32_t i = 1; i < poly2->nrings; i++)
    if (ptarray_contains_point(poly2->rings[i], pt) != LW_OUTSIDE)
      return lw_dist2d_ptarray_ptarray(poly1->rings[0], poly2->rings[i], dl);

  /*5 If we have come all the way here we know that the first point of one of them is inside the other ones
   * outer ring and not in holes so we check which one is inside.*/
  pt = getPoint2d_cp(poly1->rings[0], 0);
  if (ptarray_contains_point(poly2->rings[0], pt) != LW_OUTSIDE)
  {
    lw_dist2d_distpts_set(dl, 0.0, pt, pt);
    return LW_TRUE;
  }

  pt = getPoint2d_cp(poly2->rings[0], 0);
  if (ptarray_contains_point(poly1->rings[0], pt) != LW_OUTSIDE)
  {
    lw_dist2d_distpts_set(dl, 0.0, pt, pt);
    return LW_TRUE;
  }

  lwerror("Unspecified error in function lw_dist2d_poly_poly");
  return LW_FALSE;
}

int
lw_dist2d_poly_curvepoly(LWPOLY *poly1, LWCURVEPOLY *curvepoly2, DISTPTS *dl)
{
  LWCURVEPOLY *curvepoly1 = lwcurvepoly_construct_from_lwpoly(poly1);
  int rv = lw_dist2d_curvepoly_curvepoly(curvepoly1, curvepoly2, dl);
  lwgeom_free((LWGEOM *)curvepoly1);
  return rv;
}

int
lw_dist2d_circstring_poly(LWCIRCSTRING *circ, LWPOLY *poly, DISTPTS *dl)
{
  LWCURVEPOLY *curvepoly = lwcurvepoly_construct_from_lwpoly(poly);
  int rv = lw_dist2d_line_curvepoly((LWLINE *)circ, curvepoly, dl);
  lwgeom_free((LWGEOM *)curvepoly);
  return rv;
}

int
lw_dist2d_circstring_curvepoly(LWCIRCSTRING *circ, LWCURVEPOLY *poly, DISTPTS *dl)
{
  return lw_dist2d_line_curvepoly((LWLINE *)circ, poly, dl);
}

int
lw_dist2d_circstring_circstring(LWCIRCSTRING *line1, LWCIRCSTRING *line2, DISTPTS *dl)
{
  return lw_dist2d_ptarrayarc_ptarrayarc(line1->points, line2->points, dl);
}

int
lw_dist2d_curvepoly_curvepoly(LWCURVEPOLY *poly1, LWCURVEPOLY *poly2, DISTPTS *dl)
{
  const POINT2D *pt;
  POINT2D scratch1, scratch2;

  /*1 if we are looking for maxdistance, just check the outer rings.*/
  if (dl->mode == DIST_MAX)
    return lw_dist2d_recursive(poly1->rings[0], poly2->rings[0], dl);

  /* 2  check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings
  here it would be possible to handle the information about which one is inside which one and only search for the
  smaller ones in the bigger ones holes.*/
  pt = lw_curvering_getfirstpoint2d_cp(poly1->rings[0], &scratch1);
  if (lwgeom_contains_point(poly2->rings[0], pt) == LW_OUTSIDE)
  {
    pt = lw_curvering_getfirstpoint2d_cp(poly2->rings[0], &scratch2);
    if (lwgeom_contains_point(poly1->rings[0], pt) == LW_OUTSIDE)
      return lw_dist2d_recursive(poly1->rings[0], poly2->rings[0], dl);
  }

  /*3 check if first point of poly2 is in a hole of poly1. If so check outer ring of poly2 against that hole
   * of poly1*/
  pt = lw_curvering_getfirstpoint2d_cp(poly2->rings[0], &scratch2);
  for (uint32_t i = 1; i < poly1->nrings; i++)
    if (lwgeom_contains_point(poly1->rings[i], pt) != LW_OUTSIDE)
      return lw_dist2d_recursive(poly1->rings[i], poly2->rings[0], dl);

  /*4 check if first point of poly1 is in a hole of poly2. If so check outer ring of poly1 against that hole
   * of poly2*/
  pt = lw_curvering_getfirstpoint2d_cp(poly1->rings[0], &scratch1);
  for (uint32_t i = 1; i < poly2->nrings; i++)
    if (lwgeom_contains_point(poly2->rings[i], pt) != LW_OUTSIDE)
      return lw_dist2d_recursive(poly1->rings[0], poly2->rings[i], dl);

  /*5 If we have come all the way here we know that the first point of one of them is inside the other ones
   * outer ring and not in holes so we check which one is inside.*/
  pt = lw_curvering_getfirstpoint2d_cp(poly1->rings[0], &scratch1);
  if (lwgeom_contains_point(poly2->rings[0], pt) != LW_OUTSIDE)
  {
    lw_dist2d_distpts_set(dl, 0.0, pt, pt);
    return LW_TRUE;
  }

  pt = lw_curvering_getfirstpoint2d_cp(poly2->rings[0], &scratch2);
  if (lwgeom_contains_point(poly1->rings[0], pt) != LW_OUTSIDE)
  {
    lw_dist2d_distpts_set(dl, 0.0, pt, pt);
    return LW_TRUE;
  }

  lwerror("Unspecified error in function lw_dist2d_curvepoly_curvepoly");
  return LW_FALSE;
}

/**
 * search all the segments of pointarray to see which one is closest to p1
 * Returns minimum distance between point and pointarray
 */
int
lw_dist2d_pt_ptarray(const POINT2D *p, POINTARRAY *pa, DISTPTS *dl)
{
  const POINT2D *start, *end;
  int twist = dl->twisted;

  start = getPoint2d_cp(pa, 0);

  LWDEBUG(2, "lw_dist2d_pt_ptarray enter");

  if (!lw_dist2d_pt_pt(p, start, dl))
    return LW_FALSE;

  LWDEBUGF(2, "lw_dist2d_pt_ptarray: distance from first point ? : %.15g", dl->distance);

  for (uint32_t t = 1; t < pa->npoints; t++)
  {
    dl->twisted = twist;
    end = getPoint2d_cp(pa, t);
    if (!lw_dist2d_pt_seg(p, start, end, dl))
      return LW_FALSE;

    LWDEBUGF(2, "lw_dist2d_pt_ptarray: distance from seg %u ? : %.15g", t, dl->distance);

    if (dl->distance <= dl->tolerance && dl->mode == DIST_MIN)
      return LW_TRUE; /*just a check if the answer is already given*/
    start = end;
  }

  return LW_TRUE;
}

/**
 * Search all the arcs of pointarray to see which one is closest to p1
 * Returns minimum distance between point and arc pointarray.
 */
int
lw_dist2d_pt_ptarrayarc(const POINT2D *p, const POINTARRAY *pa, DISTPTS *dl)
{
  uint32_t t;
  const POINT2D *A1;
  const POINT2D *A2;
  const POINT2D *A3;
  int twist = dl->twisted;

  LWDEBUG(2, "lw_dist2d_pt_ptarrayarc is called");

  if (pa->npoints % 2 == 0 || pa->npoints < 3)
  {
    lwerror("lw_dist2d_pt_ptarrayarc called with non-arc input");
    return LW_FALSE;
  }

  if (dl->mode == DIST_MAX)
  {
    lwerror("lw_dist2d_pt_ptarrayarc does not currently support DIST_MAX mode");
    return LW_FALSE;
  }

  A1 = getPoint2d_cp(pa, 0);

  if (!lw_dist2d_pt_pt(p, A1, dl))
    return LW_FALSE;

  for (t = 1; t < pa->npoints; t += 2)
  {
    dl->twisted = twist;
    A2 = getPoint2d_cp(pa, t);
    A3 = getPoint2d_cp(pa, t + 1);

    if (lw_dist2d_pt_arc(p, A1, A2, A3, dl) == LW_FALSE)
      return LW_FALSE;

    if (dl->distance <= dl->tolerance && dl->mode == DIST_MIN)
      return LW_TRUE; /*just a check if the answer is already given*/

    A1 = A3;
  }

  return LW_TRUE;
}

/**
 * test each segment of l1 against each segment of l2.
 */
int
lw_dist2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS *dl)
{
  uint32_t t, u;
  const POINT2D *start, *end;
  const POINT2D *start2, *end2;
  int twist = dl->twisted;

  LWDEBUGF(2, "lw_dist2d_ptarray_ptarray called (points: %d-%d)", l1->npoints, l2->npoints);

  /*  If we are searching for maxdistance we go straight to point-point calculation since the maxdistance have
   * to be between two vertices*/
  if (dl->mode == DIST_MAX)
  {
    for (t = 0; t < l1->npoints; t++) /*for each segment in L1 */
    {
      start = getPoint2d_cp(l1, t);
      for (u = 0; u < l2->npoints; u++) /*for each segment in L2 */
      {
        start2 = getPoint2d_cp(l2, u);
        lw_dist2d_pt_pt(start, start2, dl);
      }
    }
  }
  else
  {
    start = getPoint2d_cp(l1, 0);
    for (t = 1; t < l1->npoints; t++) /*for each segment in L1 */
    {
      end = getPoint2d_cp(l1, t);
      start2 = getPoint2d_cp(l2, 0);
      for (u = 1; u < l2->npoints; u++) /*for each segment in L2 */
      {
        end2 = getPoint2d_cp(l2, u);
        dl->twisted = twist;
        lw_dist2d_seg_seg(start, end, start2, end2, dl);
        if (dl->distance <= dl->tolerance && dl->mode == DIST_MIN)
          return LW_TRUE; /*just a check if the answer is already given*/
        start2 = end2;
      }
      start = end;
    }
  }
  return LW_TRUE;
}

/**
 * Test each segment of pa against each arc of pb for distance.
 */
int
lw_dist2d_ptarray_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DISTPTS *dl)
{
  uint32_t t, u;
  const POINT2D *A1;
  const POINT2D *A2;
  const POINT2D *B1;
  const POINT2D *B2;
  const POINT2D *B3;
  int twist = dl->twisted;

  LWDEBUGF(2, "lw_dist2d_ptarray_ptarrayarc called (points: %d-%d)", pa->npoints, pb->npoints);

  if (pb->npoints % 2 == 0 || pb->npoints < 3)
  {
    lwerror("lw_dist2d_ptarray_ptarrayarc called with non-arc input");
    return LW_FALSE;
  }

  if (dl->mode == DIST_MAX)
  {
    lwerror("lw_dist2d_ptarray_ptarrayarc does not currently support DIST_MAX mode");
    return LW_FALSE;
  }
  else
  {
    A1 = getPoint2d_cp(pa, 0);
    for (t = 1; t < pa->npoints; t++) /* For each segment in pa */
    {
      A2 = getPoint2d_cp(pa, t);
      B1 = getPoint2d_cp(pb, 0);
      for (u = 1; u < pb->npoints; u += 2) /* For each arc in pb */
      {
        B2 = getPoint2d_cp(pb, u);
        B3 = getPoint2d_cp(pb, u + 1);
        dl->twisted = twist;

        lw_dist2d_seg_arc(A1, A2, B1, B2, B3, dl);

        /* If we've found a distance within tolerance, we're done */
        if (dl->distance <= dl->tolerance && dl->mode == DIST_MIN)
          return LW_TRUE;

        B1 = B3;
      }
      A1 = A2;
    }
  }
  return LW_TRUE;
}

/**
 * Test each arc of pa against each arc of pb for distance.
 */
int
lw_dist2d_ptarrayarc_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DISTPTS *dl)
{
  uint32_t t, u;
  const POINT2D *A1;
  const POINT2D *A2;
  const POINT2D *A3;
  const POINT2D *B1;
  const POINT2D *B2;
  const POINT2D *B3;
  int twist = dl->twisted;

  LWDEBUGF(2, "lw_dist2d_ptarrayarc_ptarrayarc called (points: %d-%d)", pa->npoints, pb->npoints);

  if (dl->mode == DIST_MAX)
  {
    lwerror("lw_dist2d_ptarrayarc_ptarrayarc does not currently support DIST_MAX mode");
    return LW_FALSE;
  }
  else
  {
    A1 = getPoint2d_cp(pa, 0);
    for (t = 1; t < pa->npoints; t += 2) /* For each segment in pa */
    {
      A2 = getPoint2d_cp(pa, t);
      A3 = getPoint2d_cp(pa, t + 1);
      B1 = getPoint2d_cp(pb, 0);
      for (u = 1; u < pb->npoints; u += 2) /* For each arc in pb */
      {
        B2 = getPoint2d_cp(pb, u);
        B3 = getPoint2d_cp(pb, u + 1);
        dl->twisted = twist;

        lw_dist2d_arc_arc(A1, A2, A3, B1, B2, B3, dl);

        /* If we've found a distance within tolerance, we're done */
        if (dl->distance <= dl->tolerance && dl->mode == DIST_MIN)
          return LW_TRUE;

        B1 = B3;
      }
      A1 = A3;
    }
  }
  return LW_TRUE;
}

/**
 * Calculate the shortest distance between an arc and an edge.
 * Line/circle approach from http://stackoverflow.com/questions/1073336/circle-line-collision-detection
 */
int
lw_dist2d_seg_arc(const POINT2D *A1,
      const POINT2D *A2,
      const POINT2D *B1,
      const POINT2D *B2,
      const POINT2D *B3,
      DISTPTS *dl)
{
  POINT2D C;       /* center of arc circle */
  double radius_C; /* radius of arc circle */
  POINT2D D;       /* point on A closest to C */
  double dist_C_D; /* distance from C to D */
  int pt_in_arc, pt_in_seg;
  DISTPTS dltmp;

  /* Bail out on crazy modes */
  if (dl->mode < 0)
    lwerror("lw_dist2d_seg_arc does not support maxdistance mode");

  /* What if the "arc" is a point? */
  if (lw_arc_is_pt(B1, B2, B3))
    return lw_dist2d_pt_seg(B1, A1, A2, dl);

  /* Calculate center and radius of the circle. */
  radius_C = lw_arc_center(B1, B2, B3, &C);

  /* This "arc" is actually a line (B2 is collinear with B1,B3) */
  if (radius_C < 0.0)
    return lw_dist2d_seg_seg(A1, A2, B1, B3, dl);

  /* Calculate distance between the line and circle center */
  lw_dist2d_distpts_init(&dltmp, DIST_MIN);
  if (lw_dist2d_pt_seg(&C, A1, A2, &dltmp) == LW_FALSE)
    lwerror("lw_dist2d_pt_seg failed in lw_dist2d_seg_arc");

  D = dltmp.p1;
  dist_C_D = dltmp.distance;

  /* Line intersects circle, maybe arc intersects edge? */
  /* If so, that's the closest point. */
  /* If not, the closest point is one of the end points of A */
  if (dist_C_D < radius_C)
  {
    double length_A;  /* length of the segment A */
    POINT2D E, F;     /* points of intersection of edge A and circle(B) */
    double dist_D_EF; /* distance from D_line to E or F (same distance both ways) */
    POINT2D D_line;   /* foot of perpendicular from C to the infinite line through A1-A2 */
    double dist_C_D_line; /* distance from C to the infinite line */
    double t;

    length_A = sqrt((A2->x - A1->x) * (A2->x - A1->x) + (A2->y - A1->y) * (A2->y - A1->y));

    /*
     * D (from lw_dist2d_pt_seg) is clamped to the segment, but we need the
     * foot of the perpendicular on the *infinite* line to correctly compute
     * circle-line intersections E and F. Using the clamped endpoint gives wrong
     * intersection points when the perpendicular falls outside the segment.
     */
    t = ((C.x - A1->x) * (A2->x - A1->x) + (C.y - A1->y) * (A2->y - A1->y))
        / (length_A * length_A);
    D_line.x = A1->x + t * (A2->x - A1->x);
    D_line.y = A1->y + t * (A2->y - A1->y);
    dist_C_D_line = sqrt((C.x - D_line.x) * (C.x - D_line.x) + (C.y - D_line.y) * (C.y - D_line.y));

    dist_D_EF = sqrt(radius_C * radius_C - dist_C_D_line * dist_C_D_line);

    /* Point of intersection E */
    E.x = D_line.x - (A2->x - A1->x) * dist_D_EF / length_A;
    E.y = D_line.y - (A2->y - A1->y) * dist_D_EF / length_A;
    /* Point of intersection F */
    F.x = D_line.x + (A2->x - A1->x) * dist_D_EF / length_A;
    F.y = D_line.y + (A2->y - A1->y) * dist_D_EF / length_A;

    /* If E is within A and within B then it's an intersection point */
    pt_in_arc = lw_pt_in_arc(&E, B1, B2, B3);
    pt_in_seg = lw_pt_in_seg(&E, A1, A2);

    if (pt_in_arc && pt_in_seg)
    {
      lw_dist2d_distpts_set(dl, 0.0, &E, &E);
      return LW_TRUE;
    }

    /* If F is within A and within B then it's an intersection point */
    pt_in_arc = lw_pt_in_arc(&F, B1, B2, B3);
    pt_in_seg = lw_pt_in_seg(&F, A1, A2);

    if (pt_in_arc && pt_in_seg)
    {
      lw_dist2d_distpts_set(dl, 0.0, &F, &F);
      return LW_TRUE;
    }
  }

  /* Line grazes circle, maybe arc intersects edge? */
  /* If so, grazing point is the closest point. */
  /* If not, the closest point is one of the end points of A */
  else if (dist_C_D == radius_C)
  {
    /* Closest point D is also the point of grazing */
    pt_in_arc = lw_pt_in_arc(&D, B1, B2, B3);
    pt_in_seg = lw_pt_in_seg(&D, A1, A2);

    /* Is D contained in both A and B? */
    if (pt_in_arc && pt_in_seg)
    {
      lw_dist2d_distpts_set(dl, 0.0, &D, &D);
      return LW_TRUE;
    }
  }
  /* Line misses circle. */
  /* If closest point to A on circle is within B, then that's the closest */
  /* Otherwise, the closest point will be an end point of A */
  else
  {
    POINT2D G; /* Point on circle closest to A */
    G.x = C.x + (D.x - C.x) * radius_C / dist_C_D;
    G.y = C.y + (D.y - C.y) * radius_C / dist_C_D;

    pt_in_arc = lw_pt_in_arc(&G, B1, B2, B3);
    pt_in_seg = lw_pt_in_seg(&D, A1, A2);

    /* Closest point is on the interior of A and B */
    if (pt_in_arc && pt_in_seg)
      return lw_dist2d_pt_pt(&D, &G, dl);
  }

  /* Now we test the many combinations of end points with either */
  /* arcs or edges. Each previous check determined if the closest */
  /* potential point was within the arc/segment inscribed on the */
  /* line/circle holding the arc/segment. */

  /* Closest point is in the arc, but not in the segment, so */
  /* one of the segment end points must be the closest. */
  if (pt_in_arc && !pt_in_seg)
  {
    lw_dist2d_pt_arc(A1, B1, B2, B3, dl);
    lw_dist2d_pt_arc(A2, B1, B2, B3, dl);
    return LW_TRUE;
  }
  /* or, one of the arc end points is the closest */
  else if (pt_in_seg && !pt_in_arc)
  {
    lw_dist2d_pt_seg(B1, A1, A2, dl);
    lw_dist2d_pt_seg(B3, A1, A2, dl);
    return LW_TRUE;
  }
  /* Finally, one of the end-point to end-point combos is the closest. */
  else
  {
    lw_dist2d_pt_pt(A1, B1, dl);
    lw_dist2d_pt_pt(A1, B3, dl);
    lw_dist2d_pt_pt(A2, B1, dl);
    lw_dist2d_pt_pt(A2, B3, dl);
    return LW_TRUE;
  }

  return LW_FALSE;
}

int
lw_dist2d_pt_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, DISTPTS *dl)
{
  double radius_A, d;
  POINT2D C; /* center of circle defined by arc A */
  POINT2D X; /* point circle(A) where line from C to P crosses */

  if (dl->mode < 0)
    lwerror("lw_dist2d_pt_arc does not support maxdistance mode");

  /* What if the arc is a point? */
  if (lw_arc_is_pt(A1, A2, A3))
    return lw_dist2d_pt_pt(P, A1, dl);

  /* Calculate centers and radii of circles. */
  radius_A = lw_arc_center(A1, A2, A3, &C);

  /* This "arc" is actually a line (A2 is colinear with A1,A3) */
  if (radius_A < 0.0)
    return lw_dist2d_pt_seg(P, A1, A3, dl);

  /* Distance from point to center */
  d = distance2d_pt_pt(&C, P);

  /* P is the center of the circle */
  if (FP_EQUALS(d, 0.0))
  {
    lw_dist2d_distpts_set(dl, radius_A, A1, P);
    return LW_TRUE;
  }

  /* X is the point on the circle where the line from P to C crosses */
  X.x = C.x + (P->x - C.x) * radius_A / d;
  X.y = C.y + (P->y - C.y) * radius_A / d;

  /* Is crossing point inside the arc? Or arc is actually circle? */
  if (p2d_same(A1, A3) || lw_pt_in_arc(&X, A1, A2, A3))
  {
    lw_dist2d_pt_pt(P, &X, dl);
  }
  else
  {
    /* Distance is the minimum of the distances to the arc end points */
    lw_dist2d_pt_pt(A1, P, dl);
    lw_dist2d_pt_pt(A3, P, dl);
  }
  return LW_TRUE;
}



/**
 * @brief Calculates the intersection points of a circle and line.
 *
 * This function assumes the center and test point are distinct.
 * Finds the line between the circle center and test point, and
 * the two intersection points on the circle defined by that line.
 * If those points fall within the provided arc (A1,A2,A3) then
 * the points are added to the provided array (I) and the array
 * length counter is incremented.
 *
 * @param A1   Point of arc
 * @param A2   Point of arc
 * @param A3   Point of arc
 * @param center_A   Center of arc A circle
 * @param radius_A   Radius of arc A circle
 * @param P    Point to use in calculating intersection line
 * @param I    [out] Pointer to an return value array
 * @param ni   [out] Pointer to return array size counter
 * @return int
 */
static int
lw_dist2d_circle_intersections(
  const POINT2D *A1, const POINT2D *A2, const POINT2D *A3,
  const POINT2D *center_A,
  double radius_A,
  const POINT2D *P,
  POINT2D *I,
  uint32_t *ni)
{
  POINT2D R;

  // If the test point is on the center of the other
  // arc, some other point has to be closer, by definition.
  if (p2d_same(center_A, P))
    return 0;

  // Calculate vector from the center to the pt
  double dir_x = center_A->x - P->x;
  double dir_y = center_A->y - P->y;
  double dist = sqrt(dir_x * dir_x + dir_y * dir_y);

  // Normalize the direction vector to get a unit vector (length = 1)
  double unit_x = dir_x / dist;
  double unit_y = dir_y / dist;

  // Calculate the two intersection points on the circle
  // Point 1: Move from center_A along the unit vector by distance radius_A
  R.x = center_A->x + unit_x * radius_A;
  R.y = center_A->y + unit_y * radius_A;
  if (lw_pt_in_arc(&R, A1, A2, A3))
    I[(*ni)++] = R;

  // Point 2: Move from center_A against the unit vector by distance radius_A
  R.x = center_A->x - unit_x * radius_A;
  R.y = center_A->y - unit_y * radius_A;
  if (lw_pt_in_arc(&R, A1, A2, A3))
    I[(*ni)++] = R;

  return 0;
}


/**
 * @brief Calculates the intersection points of two overlapping circles.
 *
 * This function assumes the circles are known to intersect at one or two points.
 * Specifically, the distance 'd' between their centers must satisfy:
 * d < (rA + rB)  and  d > fabs(rA - rB)
 * If these conditions are not met, the results are undefined.
 *
 * @param cA   [in] Pointer to the center point of the first circle (A).
 * @param rA   [in] The radius of the first circle (A).
 * @param cB   [in] Pointer to the center point of the second circle (B).
 * @param rB   [in] The radius of the second circle (B).
 * @param I    [out] Pointer to an array of at least 2 POINT2D structs to store the results.
 * I[0] will contain the first intersection point.
 * If a second exists, it will be in I[1].
 * @return int The number of intersection points found (1 or 2). Returns 0 if
 * the centers are coincident or another error occurs.
 */
static uint32_t
lw_dist2d_circle_circle_intersections(
  const POINT2D *cA, double rA,
  const POINT2D *cB, double rB,
  POINT2D *I)
{
  // Vector from center A to center B
  double dx = cB->x - cA->x;
  double dy = cB->y - cA->y;

  // Distance between the centers
  double d = sqrt(dx * dx + dy * dy);

  // 'a' is the distance from the center of circle A to the point P,
  // which is the base of the right triangle formed by cA, P, and an intersection point.
  double a = (rA * rA - rB * rB + d * d) / (2.0 * d);

  // 'h' is the height of that right triangle.
  double h_squared = rA * rA - a * a;

  // Due to floating point errors, h_squared can be slightly negative.
  // This happens when the circles are perfectly tangent. Clamp to 0.
  if (h_squared < 0.0)
    h_squared = 0.0;

  double h = sqrt(h_squared);

  // Find the coordinates of point P
  double Px = cA->x + a * (dx / d);
  double Py = cA->y + a * (dy / d);

  // The two intersection points are found by moving from P by a distance 'h'
  // in directions perpendicular to the line connecting the centers.
  // The perpendicular vector to (dx, dy) is (-dy, dx).

  // Intersection point 1
  I[0].x = Px - h * (dy / d);
  I[0].y = Py + h * (dx / d);

  // If h is very close to 0, the circles are tangent and there's only one intersection point.
  if (FP_IS_ZERO(h))
    return 1;

  // Intersection point 2
  I[1].x = Px + h * (dy / d);
  I[1].y = Py - h * (dx / d);

  return 2;
}


int
lw_dist2d_arc_arc(
  const POINT2D *A1, const POINT2D *A2, const POINT2D *A3,
  const POINT2D *B1, const POINT2D *B2, const POINT2D *B3,
  DISTPTS *dl)
{
  POINT2D pts_A[16];            /* Points on A that might be the nearest */
  POINT2D pts_B[16];            /* Points on B that might be the nearest */
  uint32_t ai = 0, bi = 0;      /* Number of points in pts_A and pts_B */
  POINT2D center_A, center_B;   /* Center points of arcs A and B */
  double radius_A, radius_B, d; /* Radii of arcs A and B */
  int is_disjoint, is_overlapping, is_contained, is_same_center;
  POINT2D intersectionPts[2];

  if (dl->mode != DIST_MIN)
    lwerror("lw_dist2d_arc_arc only supports mindistance");

  /* What if one or both of our "arcs" is actually a point? */
  if (lw_arc_is_pt(B1, B2, B3) && lw_arc_is_pt(A1, A2, A3))
    return lw_dist2d_pt_pt(B1, A1, dl);
  else if (lw_arc_is_pt(B1, B2, B3))
    return lw_dist2d_pt_arc(B1, A1, A2, A3, dl);
  else if (lw_arc_is_pt(A1, A2, A3))
    return lw_dist2d_pt_arc(A1, B1, B2, B3, dl);

  /* Calculate centers and radii of circles. */
  radius_A = lw_arc_center(A1, A2, A3, &center_A);
  radius_B = lw_arc_center(B1, B2, B3, &center_B);

  /* Two co-linear arcs?!? That's two segments. */
  if (radius_A < 0 && radius_B < 0)
    return lw_dist2d_seg_seg(A1, A3, B1, B3, dl);

  /* A is co-linear, delegate to lw_dist_seg_arc here. */
  if (radius_A < 0)
    return lw_dist2d_seg_arc(A1, A3, B1, B2, B3, dl);

  /* B is co-linear, delegate to lw_dist_seg_arc here. */
  if (radius_B < 0)
    return lw_dist2d_seg_arc(B1, B3, A1, A2, A3, dl);

  /* Circle relationships */
  d = distance2d_pt_pt(&center_A, &center_B);
  is_disjoint = (d > (radius_A + radius_B));
  is_contained = (d < fabs(radius_A - radius_B));
  is_same_center = p2d_same(&center_A, &center_B);
  is_overlapping = ! (is_disjoint || is_contained || is_same_center);

  /*
   * Prime the array of potential closest points with the
   * arc end points, which frequently participate in closest
   * points.
   */
  pts_A[ai++] = *A1;
  pts_A[ai++] = *A3;
  pts_B[bi++] = *B1;
  pts_B[bi++] = *B3;

  /*
   * Overlapping circles might have a zero distance
   * case if the circle intersection points are inside both
   * arcs.
   */
  if (is_overlapping)
  {
    /*
     * Find the two points the circles intersect at.
     */
    uint32_t npoints = lw_dist2d_circle_circle_intersections(
      &center_A, radius_A,
      &center_B, radius_B,
      intersectionPts);
    for (uint32_t i = 0; i < npoints; i++)
    {
      /*
       * If an intersection point is contained in both
       * arcs, that is a location of zero distance, so
       * we are done calculating.
       */
      if (lw_pt_in_arc(&intersectionPts[i], A1, A2, A3) &&
        lw_pt_in_arc(&intersectionPts[i], B1, B2, B3))
      {
        lw_dist2d_distpts_set(dl, 0.0, &intersectionPts[i], &intersectionPts[i]);
        return LW_TRUE;
      }
    }
  }

  /*
   * Join the circle centers and find the places that
   * line intersects the circles. Where those places
   * are in the arcs, they are potential sites of the
   * closest points.
   */
  if (is_disjoint || is_contained || is_overlapping)
  {
    if (!is_same_center)
    {
      /* Add points on A that intersect line from center_A to center_B */
      lw_dist2d_circle_intersections(
          A1, A2, A3,
          &center_A, radius_A, &center_B,
          pts_A, &ai);

      /* Add points on B that intersect line from center_B to center_A */
      lw_dist2d_circle_intersections(
          B1, B2, B3,
          &center_B, radius_B, &center_A,
          pts_B, &bi);
    }

    /* Add points on A that intersect line to B1 */
    lw_dist2d_circle_intersections(
        A1, A2, A3,
        &center_A, radius_A, B1,
        pts_A, &ai);

    /* Add points on A that intersect line to B3 */
    lw_dist2d_circle_intersections(
        A1, A2, A3,
        &center_A, radius_A, B3,
        pts_A, &ai);

    /* Add points on B that intersect line to A1 */
    lw_dist2d_circle_intersections(
        B1, B2, B3,
        &center_B, radius_B, A1,
        pts_B, &bi);

    /* Add points on B that intersect line to A3 */
    lw_dist2d_circle_intersections(
        B1, B2, B3,
        &center_B, radius_B, A3,
        pts_B, &bi);
  }

  /*
   * Now just brute force check all pairs of participating
   * points, to find the pair that is closest together.
   */
  for (uint32_t i = 0; i < ai; i++)
    for (uint32_t j = 0; j < bi; j++)
      lw_dist2d_pt_pt(&pts_A[i], &pts_B[j], dl);

  return LW_TRUE;
}


/**
Finds the shortest distance between two segments.
This function is changed so it is not doing any comparison of distance
but just sending every possible combination further to lw_dist2d_pt_seg
*/
int
lw_dist2d_seg_seg(const POINT2D *A, const POINT2D *B, const POINT2D *C, const POINT2D *D, DISTPTS *dl)
{
  double s_top, s_bot, s;
  double r_top, r_bot, r;

  /*A and B are the same point */
  if ((A->x == B->x) && (A->y == B->y))
  {
    return lw_dist2d_pt_seg(A, C, D, dl);
  }

  /*U and V are the same point */
  if ((C->x == D->x) && (C->y == D->y))
  {
    dl->twisted = ((dl->twisted) * (-1));
    return lw_dist2d_pt_seg(D, A, B, dl);
  }

  /* AB and CD are line segments */
  /* from comp.graphics.algo

  Solving the above for r and s yields
        (Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy)
       r = ----------------------------- (eqn 1)
        (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)

      (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
    s = ----------------------------- (eqn 2)
      (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
  Let P be the position vector of the intersection point, then
    P=A+r(B-A) or
    Px=Ax+r(Bx-Ax)
    Py=Ay+r(By-Ay)
  By examining the values of r & s, you can also determine some other limiting conditions:
    If 0<=r<=1 & 0<=s<=1, intersection exists
    r<0 or r>1 or s<0 or s>1 line segments do not intersect
    If the denominator in eqn 1 is zero, AB & CD are parallel
    If the numerator in eqn 1 is also zero, AB & CD are collinear.

  */
  r_top = (A->y - C->y) * (D->x - C->x) - (A->x - C->x) * (D->y - C->y);
  r_bot = (B->x - A->x) * (D->y - C->y) - (B->y - A->y) * (D->x - C->x);

  s_top = (A->y - C->y) * (B->x - A->x) - (A->x - C->x) * (B->y - A->y);
  s_bot = (B->x - A->x) * (D->y - C->y) - (B->y - A->y) * (D->x - C->x);

  if ((r_bot == 0) || (s_bot == 0))
  {
    if ((lw_dist2d_pt_seg(A, C, D, dl)) && (lw_dist2d_pt_seg(B, C, D, dl)))
    {
      /* change the order of inputted geometries and that we notice by changing sign on dl->twisted*/
      dl->twisted *= -1;
      return ((lw_dist2d_pt_seg(C, A, B, dl)) && (lw_dist2d_pt_seg(D, A, B, dl)));
    }
    else
      return LW_FALSE; /* if any of the calls to lw_dist2d_pt_seg goes wrong we return false*/
  }

  s = s_top / s_bot;
  r = r_top / r_bot;

  if (((r < 0) || (r > 1) || (s < 0) || (s > 1)) || (dl->mode == DIST_MAX))
  {
    if ((lw_dist2d_pt_seg(A, C, D, dl)) && (lw_dist2d_pt_seg(B, C, D, dl)))
    {
      /* change the order of inputted geometries and that we notice by changing sign on dl->twisted*/
      dl->twisted *= -1;
      return ((lw_dist2d_pt_seg(C, A, B, dl)) && (lw_dist2d_pt_seg(D, A, B, dl)));
    }
    else
      return LW_FALSE; /* if any of the calls to lw_dist2d_pt_seg goes wrong we return false*/
  }
  else
  {
    /* If there is intersection we identify the intersection point and return it but only if we are looking
     * for mindistance */
    if (dl->mode == DIST_MIN)
    {
      POINT2D theP;

      if (((A->x == C->x) && (A->y == C->y)) || ((A->x == D->x) && (A->y == D->y)))
      {
        theP.x = A->x;
        theP.y = A->y;
      }
      else if (((B->x == C->x) && (B->y == C->y)) || ((B->x == D->x) && (B->y == D->y)))
      {
        theP.x = B->x;
        theP.y = B->y;
      }
      else
      {
        theP.x = A->x + r * (B->x - A->x);
        theP.y = A->y + r * (B->y - A->y);
      }
      lw_dist2d_distpts_set(dl, 0, &theP, &theP);
    }
    return LW_TRUE;
  }
}

/*------------------------------------------------------------------------------------------------------------
End of Brute force functions
--------------------------------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------------------------------------------
New faster distance calculations
--------------------------------------------------------------------------------------------------------------*/

/**

The new faster calculation comparing pointarray to another pointarray
the arrays can come from both polygons and linestrings.
The naming is not good but comes from that it compares a
chosen selection of the points not all of them
*/
int
lw_dist2d_fast_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS *dl, GBOX *box1, GBOX *box2)
{
  /*here we define two lists to hold our calculated "z"-values and the order number in the geometry*/

  double k, thevalue;
  float deltaX, deltaY, c1m, c2m;
  POINT2D c1, c2;
  const POINT2D *theP;
  float min1X, max1X, max1Y, min1Y, min2X, max2X, max2Y, min2Y;
  int t;
  int n1 = l1->npoints;
  int n2 = l2->npoints;

  LISTSTRUCT *list1, *list2;
  list1 = (LISTSTRUCT *)lwalloc(sizeof(LISTSTRUCT) * n1);
  list2 = (LISTSTRUCT *)lwalloc(sizeof(LISTSTRUCT) * n2);

  LWDEBUG(2, "lw_dist2d_fast_ptarray_ptarray is called");

  max1X = box1->xmax;
  min1X = box1->xmin;
  max1Y = box1->ymax;
  min1Y = box1->ymin;
  max2X = box2->xmax;
  min2X = box2->xmin;
  max2Y = box2->ymax;
  min2Y = box2->ymin;
  /*we want the center of the bboxes, and calculate the slope between the centerpoints*/
  c1.x = min1X + (max1X - min1X) / 2;
  c1.y = min1Y + (max1Y - min1Y) / 2;
  c2.x = min2X + (max2X - min2X) / 2;
  c2.y = min2Y + (max2Y - min2Y) / 2;

  deltaX = (c2.x - c1.x);
  deltaY = (c2.y - c1.y);

  /*Here we calculate where the line perpendicular to the center-center line crosses the axes for each vertex
  if the center-center line is vertical the perpendicular line will be horizontal and we find it's crossing the
  Y-axes with z = y-kx */
  if ((deltaX * deltaX) < (deltaY * deltaY)) /*North or South*/
  {
    k = -deltaX / deltaY;
    for (t = 0; t < n1; t++) /*for each segment in L1 */
    {
      theP = getPoint2d_cp(l1, t);
      thevalue = theP->y - (k * theP->x);
      list1[t].themeasure = thevalue;
      list1[t].pnr = t;
    }
    for (t = 0; t < n2; t++) /*for each segment in L2*/
    {
      theP = getPoint2d_cp(l2, t);
      thevalue = theP->y - (k * theP->x);
      list2[t].themeasure = thevalue;
      list2[t].pnr = t;
    }
    c1m = c1.y - (k * c1.x);
    c2m = c2.y - (k * c2.x);
  }

  /*if the center-center line is horizontal the perpendicular line will be vertical. To eliminate problems with
   dividing by zero we are here mirroring the coordinate-system and we find it's crossing the X-axes with z =
   x-(1/k)y */
  else /*West or East*/
  {
    k = -deltaY / deltaX;
    for (t = 0; t < n1; t++) /*for each segment in L1 */
    {
      theP = getPoint2d_cp(l1, t);
      thevalue = theP->x - (k * theP->y);
      list1[t].themeasure = thevalue;
      list1[t].pnr = t;
      /* lwnotice("l1 %d, measure=%f",t,thevalue ); */
    }
    for (t = 0; t < n2; t++) /*for each segment in L2*/
    {
      theP = getPoint2d_cp(l2, t);
      thevalue = theP->x - (k * theP->y);
      list2[t].themeasure = thevalue;
      list2[t].pnr = t;
      /* lwnotice("l2 %d, measure=%f",t,thevalue ); */
    }
    c1m = c1.x - (k * c1.y);
    c2m = c2.x - (k * c2.y);
  }

  /*we sort our lists by the calculated values*/
  qsort(list1, n1, sizeof(LISTSTRUCT), struct_cmp_by_measure);
  qsort(list2, n2, sizeof(LISTSTRUCT), struct_cmp_by_measure);

  if (c1m < c2m)
  {
    if (!lw_dist2d_pre_seg_seg(l1, l2, list1, list2, k, dl))
    {
      lwfree(list1);
      lwfree(list2);
      return LW_FALSE;
    }
  }
  else
  {
    dl->twisted = ((dl->twisted) * (-1));
    if (!lw_dist2d_pre_seg_seg(l2, l1, list2, list1, k, dl))
    {
      lwfree(list1);
      lwfree(list2);
      return LW_FALSE;
    }
  }
  lwfree(list1);
  lwfree(list2);
  return LW_TRUE;
}

int
struct_cmp_by_measure(const void *a, const void *b)
{
  LISTSTRUCT *ia = (LISTSTRUCT *)a;
  LISTSTRUCT *ib = (LISTSTRUCT *)b;
  return (ia->themeasure > ib->themeasure) ? 1 : ((ia->themeasure < ib->themeasure) ? -1 : 0);
}

/** preparation before lw_dist2d_seg_seg. */
int
lw_dist2d_pre_seg_seg(POINTARRAY *l1, POINTARRAY *l2, LISTSTRUCT *list1, LISTSTRUCT *list2, double k, DISTPTS *dl)
{
  const POINT2D *p1, *p2, *p3, *p4, *p01, *p02;
  int pnr1, pnr2, pnr3, pnr4, n1, n2, i, u, r, twist;
  double maxmeasure;
  n1 = l1->npoints;
  n2 = l2->npoints;

  LWDEBUG(2, "lw_dist2d_pre_seg_seg is called");

  p1 = getPoint2d_cp(l1, list1[0].pnr);
  p3 = getPoint2d_cp(l2, list2[0].pnr);
  lw_dist2d_pt_pt(p1, p3, dl);
  maxmeasure = sqrt(dl->distance * dl->distance + (dl->distance * dl->distance * k * k));
  twist = dl->twisted; /*to keep the incoming order between iterations*/
  for (i = (n1 - 1); i >= 0; --i)
  {
    /*we break this iteration when we have checked every point closer to our perpendicular "checkline" than
     * our shortest found distance*/
    if (((list2[0].themeasure - list1[i].themeasure)) > maxmeasure)
      break;
    /*because we are not iterating in the original point order we have to check the segment before and after
     * every point*/
    for (r = -1; r <= 1; r += 2)
    {
      pnr1 = list1[i].pnr;
      p1 = getPoint2d_cp(l1, pnr1);
      if (pnr1 + r < 0)
      {
        p01 = getPoint2d_cp(l1, (n1 - 1));
        if ((p1->x == p01->x) && (p1->y == p01->y))
          pnr2 = (n1 - 1);
        else
          pnr2 = pnr1; /* if it is a line and the last and first point is not the same we
              avoid the edge between start and end this way*/
      }

      else if (pnr1 + r > (n1 - 1))
      {
        p01 = getPoint2d_cp(l1, 0);
        if ((p1->x == p01->x) && (p1->y == p01->y))
          pnr2 = 0;
        else
          pnr2 = pnr1; /* if it is a line and the last and first point is not the same we
              avoid the edge between start and end this way*/
      }
      else
        pnr2 = pnr1 + r;

      p2 = getPoint2d_cp(l1, pnr2);
      for (u = 0; u < n2; ++u)
      {
        if (((list2[u].themeasure - list1[i].themeasure)) >= maxmeasure)
          break;
        pnr3 = list2[u].pnr;
        p3 = getPoint2d_cp(l2, pnr3);
        if (pnr3 == 0)
        {
          p02 = getPoint2d_cp(l2, (n2 - 1));
          if ((p3->x == p02->x) && (p3->y == p02->y))
            pnr4 = (n2 - 1);
          else
            pnr4 = pnr3; /* if it is a line and the last and first point is not the
                same we avoid the edge between start and end this way*/
        }
        else
          pnr4 = pnr3 - 1;

        p4 = getPoint2d_cp(l2, pnr4);
        dl->twisted = twist;
        if (!lw_dist2d_selected_seg_seg(p1, p2, p3, p4, dl))
          return LW_FALSE;

        if (pnr3 >= (n2 - 1))
        {
          p02 = getPoint2d_cp(l2, 0);
          if ((p3->x == p02->x) && (p3->y == p02->y))
            pnr4 = 0;
          else
            pnr4 = pnr3; /* if it is a line and the last and first point is not the
                same we avoid the edge between start and end this way*/
        }

        else
          pnr4 = pnr3 + 1;

        p4 = getPoint2d_cp(l2, pnr4);
        dl->twisted = twist; /*we reset the "twist" for each iteration*/
        if (!lw_dist2d_selected_seg_seg(p1, p2, p3, p4, dl))
          return LW_FALSE;
        /*here we "translate" the found mindistance so it can be compared to our "z"-values*/
        maxmeasure = sqrt(dl->distance * dl->distance + (dl->distance * dl->distance * k * k));
      }
    }
  }

  return LW_TRUE;
}

/**
  This is the same function as lw_dist2d_seg_seg but
  without any calculations to determine intersection since we
  already know they do not intersect
*/
int
lw_dist2d_selected_seg_seg(const POINT2D *A, const POINT2D *B, const POINT2D *C, const POINT2D *D, DISTPTS *dl)
{
  /*A and B are the same point */
  if ((A->x == B->x) && (A->y == B->y))
  {
    return lw_dist2d_pt_seg(A, C, D, dl);
  }
  /*U and V are the same point */

  if ((C->x == D->x) && (C->y == D->y))
  {
    dl->twisted *= -1;
    return lw_dist2d_pt_seg(D, A, B, dl);
  }

  if ((lw_dist2d_pt_seg(A, C, D, dl)) && (lw_dist2d_pt_seg(B, C, D, dl)))
  {
    /* change the order of inputted geometries and that we notice by changing sign on dl->twisted */
    dl->twisted *= -1;
    return ((lw_dist2d_pt_seg(C, A, B, dl)) && (lw_dist2d_pt_seg(D, A, B, dl)));
  }
  else
    return LW_FALSE; /* if any of the calls to lw_dist2d_pt_seg goes wrong we return false*/
}

/*------------------------------------------------------------------------------------------------------------
End of New faster distance calculations
--------------------------------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------------------------------------------
Functions in common for Brute force and new calculation
--------------------------------------------------------------------------------------------------------------*/

int
lw_dist2d_pt_seg(const POINT2D *C, const POINT2D *A, const POINT2D *B, DISTPTS *dl)
{
  POINT2D P;
  int is_vertical   = (A->x == B->x);
  int is_horizontal = (A->y == B->y);

  /*
   * If A == B, then use pt-pt distance
   */
  if (is_vertical && is_horizontal)
    return lw_dist2d_pt_pt(C, A, dl);

  /*
   * We use comp.graphics.algorithms Frequently Asked Questions method
   *
   *  Let the point be C (Cx,Cy) and the line be AB (Ax,Ay) to (Bx,By).
   *  Let P be the point of perpendicular projection of C on AB.  The parameter
   *  r, which indicates P's position along AB, is computed by the dot product
   *  of AC and AB divided by the square of the length of AB:
   *
   *  (1)     AC dot AB
   *      r = ---------
   *          ||AB||^2
   *
   *  The length of a line segment
   *
   *      L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 )
   *
   *  So (1) expands to:
   *
   *          (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
   *      r = -------------------------------
   *                        L^2
   *
   *  r has the following meaning:
   *
   *      r=0      P = A
   *      r=1      P = B
   *      r<0      P is on the backward extension of AB
   *      r>1      P is on the forward extension of AB
   *      0<r<1    P is interior to AB
   *
   */

  double dx = B->x - A->x;
  double dy = B->y - A->y;
  double L2 = (dx*dx) + (dy*dy);

  double r = ((C->x - A->x) * dx + (C->y - A->y) * dy) / L2;

  /*
   * This is for finding the maxdistance.
   * The maxdistance have to be between two vertices,
   * compared to mindistance which can be between two vertices.
   */
  if (dl->mode == DIST_MAX)
  {
    if (r >= 0.5)
      return lw_dist2d_pt_pt(C, A, dl);
    else /* (r < 0.5) */
      return lw_dist2d_pt_pt(C, B, dl);
  }

  /*
   * If projected point P is outside point A
   */
  if (r < 0)
    return lw_dist2d_pt_pt(C, A, dl);

  /*
   * If projected point P is is outside
   * point B or on point B
   */
  if (r >= 1)
    return lw_dist2d_pt_pt(C, B, dl);

  /*
   * If point C is on the segment AB then
   * distance is zero and nearest point is C
   */
  if ((A->y - C->y) * dx == (A->x - C->x) * dy)
  {
    lw_dist2d_distpts_set(dl, 0.0, C, C);
  }

  /*
   *  The point P can then be found:
   *
   *      Px = Ax + r(Bx-Ax)
   *      Py = Ay + r(By-Ay)
   *
   *  And the distance from A to P = r*L.
   */

  /*
   * If the projection of point C on the segment is
   * between A and B then we find that "point on segment"
   * and send it to lw_dist2d_pt_pt
   */
  if (is_horizontal || is_vertical)
  {
    P.x = A->x + r * dx;
    P.y = A->y + r * dy;
    return lw_dist2d_pt_pt(C, &P, dl);
  }

  /*
   *
   *  Use another parameter s to indicate the location along PC, with the
   *  following meaning:
   *         s<0      C is left of AB
   *         s>0      C is right of AB
   *         s=0      C is on AB
   *
   *  Compute s as follows:
   *
   *          (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
   *      s = -----------------------------
   *                      L^2
   *
   *
   *  Then the distance from C to P = s*L.
   */

  /*
   * Calculate distance without reference to the
   * projected point P.
   */
  double s = ((A->y - C->y) * dx - (A->x - C->x) * dy) / L2;

  double dist = fabs(s) * sqrt(L2);
  // double dist = fabs(s) * hypot(dx, dy);
  if (dist < dl->distance)
  {
    dl->distance = dist;
    {
      P.x = A->x + r * dx;
      P.y = A->y + r * dy;
      if (dl->twisted > 0)
      {
        dl->p1 = *C;
        dl->p2 = P;
      }
      else
      {
        dl->p1 = P;
        dl->p2 = *C;
      }
    }
  }
  return LW_TRUE;
}

/** Compares incoming points and stores the points closest to each other or most far away from each other depending on
 * dl->mode (max or min) */
int
lw_dist2d_pt_pt(const POINT2D *thep1, const POINT2D *thep2, DISTPTS *dl)
{
  double hside = thep2->x - thep1->x;
  double vside = thep2->y - thep1->y;
  double dist = sqrt(hside * hside + vside * vside);

  /*multiplication with mode to handle mindistance (mode=1) and maxdistance (mode = (-1)*/
  if (((dl->distance - dist) * (dl->mode)) > 0)
  {
    dl->distance = dist;

    /* To get the points in right order. twisted is updated between 1 and (-1) every time the order is
     * changed earlier in the chain*/
    if (dl->twisted > 0)
    {
      dl->p1 = *thep1;
      dl->p2 = *thep2;
    }
    else
    {
      dl->p1 = *thep2;
      dl->p2 = *thep1;
    }
  }
  return LW_TRUE;
}

/*------------------------------------------------------------------------------------------------------------
End of Functions in common for Brute force and new calculation
--------------------------------------------------------------------------------------------------------------*/

inline double
distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
{
  double hside = p2->x - p1->x;
  double vside = p2->y - p1->y;

  return hypot(hside, vside);
}

/* return distance squared, useful to avoid sqrt calculations */
double
distance2d_sqr_pt_seg(const POINT2D *C, const POINT2D *A, const POINT2D *B)
{
  /*if start==end, then use pt distance */
  if ((A->x == B->x) && (A->y == B->y))
    return distance2d_sqr_pt_pt(C, A);

  /*
   * otherwise, we use comp.graphics.algorithms
   * Frequently Asked Questions method
   *
   *  (1)       AC dot AB
   *         r = ---------
   *              ||AB||^2
   *  r has the following meaning:
   *  r=0 P = A
   *  r=1 P = B
   *  r<0 P is on the backward extension of AB
   *  r>1 P is on the forward extension of AB
   *  0<r<1 P is interior to AB
   */

  double ba_x = (B->x - A->x);
  double ba_y = (B->y - A->y);
  double ab_length_sqr = (ba_x * ba_x + ba_y * ba_y);
  double ca_x = (C->x - A->x);
  double ca_y = (C->y - A->y);
  double dot_ac_ab = (ca_x * ba_x + ca_y * ba_y);

  if (dot_ac_ab <= 0)
    return distance2d_sqr_pt_pt(C, A);
  if (dot_ac_ab >= ab_length_sqr)
    return distance2d_sqr_pt_pt(C, B);

  /*
   * (2)
   *       (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
   *  s = -----------------------------
   *                L^2
   *
   *  Then the distance from C to P = |s|*L.
   *
   */

  double s_numerator = ca_x * ba_y - ca_y * ba_x;

  /* Distance = (s_num / ab) * (s_num / ab) * ab == s_num * s_num / ab) */
  return s_numerator * s_numerator / ab_length_sqr;
}

/**
 * Compute the azimuth of segment AB in radians.
 * Return 0 on exception (same point), 1 otherwise.
 */
int
azimuth_pt_pt(const POINT2D *A, const POINT2D *B, double *d)
{
  if (A->x == B->x )
  {
    if ( A->y == B->y)
    {
      return LW_FALSE;
    }
    *d = A->y < B->y ? 0.0 : M_PI;
    return LW_TRUE;
  }
  /* TODO: use M_PI_2 instead of M_PI/2 ? */
  *d = fmod(2 * M_PI + M_PI / 2 - atan2(B->y - A->y, B->x - A->x), 2 * M_PI);

  /* TODO: convert 2*M_PI to 0 ? */
  return LW_TRUE;
}

/**
 * Azimuth is angle in radians from vertical axis
 *
 */
int
project_pt(const POINT2D *P, double distance, double azimuth, POINT2D *R)
{
  const double TWOPI = 2.0 * M_PI;
  double slope;
  /* Deal with azimuth out of (-360,360) range */
  int orbits = floor(azimuth / TWOPI);
  azimuth -= TWOPI * orbits;
  /* Convert from azimuth to conventional slope */
  slope = TWOPI - azimuth + M_PI_2;
  if (slope > 0 && slope >  TWOPI) slope -= TWOPI;
  if (slope < 0 && slope < -TWOPI) slope += TWOPI;

  double dx = cos(slope) * distance;
  double dy = sin(slope) * distance;
  R->x = P->x + dx;
  R->y = P->y + dy;
  return LW_TRUE;
}

/**
 * Azimuth is angle in radians from vertical axis
 *
 */
int
project_pt_pt(const POINT4D *A, const POINT4D *B, double distance, POINT4D *R)
{
  /* Convert from azimuth to conventional slope */
  double len = distance2d_pt_pt((const POINT2D *)A, (const POINT2D *)B);
  double prop = distance / len;
  double dx = (B->x - A->x) * prop;
  double dy = (B->y - A->y) * prop;
  double dz = (B->z - A->z) * prop;
  double dm = (B->m - A->m) * prop;
  R->x = B->x + dx;
  R->y = B->y + dy;
  if (isfinite(dz)) R->z = B->z + dz;
  if (isfinite(dm)) R->m = B->m + dm;
  return LW_TRUE;
}

Coverage Report

Created: 2026-07-16 06:52