proc xlocatePoint() may fail due to assert not p.rightOf(e)

[StefanSalewski]

We discovered a few random failures, and so commented out that line before the last one in the proc for now. The issue may occur when we insert new points with insertPoint(), and that new point is located on an existing edge. We have to investigate this issue carefully later, for now it seems to work without the assert. Actually we are using the CDT2 variant of this library, but it should be the same for initial CDT.

# Returns an edge e, s.t. either x is on e, or e is an edge of
# a triangle containing x. The search starts from startingEdge
# and proceeds in the general direction of x. Based on the
# pseudocode in Guibas and Stolfi (1985) p.121.
# http://www.karlchenofhell.org/cppswp/lischinski.pdf
# if point is in face, we return the closest edge -- triArea(e.org.point, e.dest.point, p) >= 0
proc xlocatePoint(p: Vector; startEdge: Edge): Edge =
  var e = startEdge
  assert e != nil
  while true: # TODO: can we get an infinite loop -- and how to best fix it?
    if p == e.org.point or p == e.dest.point:
      return e
    elif p.rightOf(e):
      e = e.sym
    elif not p.rightOf(e.onext):
      e = e.onext
    elif not p.rightOf(e.dprev):
      e = e.dprev
    else:
      break
  assert not p.rightOf(e)
  let a = triArea(e.org.point, e.dest.point, p)
  assert a >= 0
  if a > 0:
    e = minValueByIt([e, e.onext.sym, e.dprev.sym], linePointDistSqr(it.org.point, it.dest.point, p))
  ###assert not p.rightOf(e)
  return e

[StefanSalewski]

The reason for this issue seems to be the triArea() proc and its wrong use:

func triArea*(a, b, c: Vector): float {.inline.} =
  (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x)

which is used by other procs

func ccw*(a, b, c: Vector): bool {.inline.} =
  triArea(a, b, c) > 0

func rightOf(x: Vector; e: Edge): bool {.inline.} =
  ccw(e.dest.point, e.org.point, x)

Due to numerical inaccuracy, for three nearly collinear points x, y, z, triArea(x, y, z) may be zero,
while triArea(y, x, z) is not zero.

To avoid inaccuracy issues and full symmetry, we should change rightOf() to

func rightOf(x: Vector; e: Edge): bool {.inline.} =
  triArea(e.org.point, e.dest.point, x) < 0

We will do some more tests, and than apply this fix for sdt and sdt2 package. You may think about this issue yourself, and consult the related papers.