Ensure convexity of tessellation and enable handling of points outside [1,2]x[1,2]

examples/convexity.jl

@@ -0,0 +1,19 @@
+using VoronoiDelaunay, Gadfly
+
+points = [
+ Point2D(-1.0563841812533212, -1.4606363138997696) 
+ Point2D(0.06312982975128989, -0.48031801152366027)
+ Point2D(0.1624918689993189, -0.19919450833195906) 
+ Point2D(-1.5293344878962758, -0.7657808444340142) 
+ Point2D(0.5319064220493406, 0.6107808132476504)   
+ Point2D(-0.3670342825169435, 0.8207427582546951)  
+ Point2D(-1.9797019290444364, -0.5066353099040788) 
+ Point2D(-1.5811584920674324, 1.0346409888830976)  
+ Point2D(1.2185165319349451, 1.4177209167374605)   
+ Point2D(-1.5991536318191626, -1.3063986775765466) ];
+
+tess, ranges = DelaunayTessellation( points )
+xy = getplotxy( delaunayedges( tess, ranges ) )
+l1 = layer( x=getx.(points), y=gety.(points), Theme(default_color=colorant"orange"), Geom.point )
+l2 = layer( x=xy[1], y=xy[2], Geom.path )
+plot( l1, l2, Coord.cartesian(fixed=true) )

src/VoronoiDelaunay.jl

@@ -10,12 +10,39 @@ module VoronoiDelaunay
 # Bug reports welcome!
 
 export
-DelaunayTessellation, DelaunayTessellation2D, sizehint!, isexternal,
-min_coord, max_coord, locate, movea, moveb, movec,
-delaunayedges, voronoiedges, voronoiedgeswithoutgenerators,
-iterate, findindex, push!,
-Point, Point2D, AbstractPoint2D, getx, gety, geta, getb, getc,
-getgena, getgenb, getplotxy
+DelaunayTessellation,
+DelaunayTessellation2D,
+sizehint!,
+isexternal,
+min_coord,
+max_coord,
+locate,
+movea,
+moveb,
+movec,
+delaunayedges,
+voronoiedges,
+voronoiedgeswithoutgenerators,
+iterate,
+findindex,
+push!,
+Point,
+Point2D,
+AbstractPoint2D,
+getx,
+gety,
+geta,
+getb,
+getc,
+getgena,
+getgenb,
+getplotxy,
+scaleShiftPoints,
+expand,
+Triple,
+DelaunayTriangle
+
+
 
 using GeometricalPredicates
 import GeometricalPredicates: geta, getb, getc
@@ -24,6 +51,8 @@ import Base: push!, iterate, copy, sizehint!
 import Colors: RGB, RGBA
 using Random: shuffle!
 
+include("VoronoiDelaunayExtensions.jl")
+
 const min_coord = GeometricalPredicates.min_coord + eps(Float64)
 const max_coord = GeometricalPredicates.max_coord - eps(Float64)
 
@@ -78,17 +107,18 @@ function copy(t::DelaunayTriangle{T}) where T<:AbstractPoint2D
                      )
 end
 
-function isexternal(t::DelaunayTriangle{T}) where T<:AbstractPoint2D
-    getx(geta(t)) < min_coord || getx(geta(t)) > max_coord ||
-    getx(getb(t)) < min_coord || getx(getb(t)) > max_coord ||
-    getx(getc(t)) < min_coord || getx(getc(t)) > max_coord
+function isexternal( t::DelaunayTriangle{T}, ranges::NTuple{4,Float64} ) where T<:AbstractPoint2D
+    getx(geta(t)) < ranges[1] || getx(geta(t)) > ranges[2] ||
+    getx(getb(t)) < ranges[1] || getx(getb(t)) > ranges[2] ||
+    getx(getc(t)) < ranges[1] || getx(getc(t)) > ranges[2]
 end
 
 mutable struct DelaunayTessellation2D{T<:AbstractPoint2D}
     _trigs::Vector{DelaunayTriangle{T}}
     _last_trig_index::Int64
     _edges_to_check::Vector{Int64}
     _total_points_added::Int64
+    _ranges::NTuple{4,Float64}
 
     function DelaunayTessellation2D{T}(n::Int64 = 100) where T
         a = T(GeometricalPredicates.min_coord, GeometricalPredicates.min_coord)
@@ -102,12 +132,34 @@ mutable struct DelaunayTessellation2D{T<:AbstractPoint2D}
         t = new(_trigs, 3, Int64[], 0)
         sizehint!(t._edges_to_check, 1000)
         sizehint!(t, n)
+        _ranges=(min_coord,max_coord,min_coord,max_coord)
+        t = new(_trigs, 3, Int64[], 0, _ranges)
+        return t
     end
 end
 DelaunayTessellation2D(n::Int64) = DelaunayTessellation2D{Point2D}(n)
 DelaunayTessellation2D(n::Int64, ::T) where {T<:AbstractPoint2D} = DelaunayTessellation2D{T}(n)
 DelaunayTessellation(n::Int64=100) = DelaunayTessellation2D(n)
 
+# tessellation function that can deal with points outside [1,2]x[1,2]
+# and that always returns a convex tessellation
+function DelaunayTessellation( points::Array{T,1} ) where T<:AbstractPoint2D
+  scaledPoints, ranges = scaleShiftPoints( points )
+  scaledTess = DelaunayTessellation2D{T}(length( points ))
+  push!( scaledTess, scaledPoints )
+  tess = expand( scaledTess, ranges )
+  return tess
+end
+# in order to reduce computation, if you are interested in the edges only, you can also
+# apply the expand function to the points of the edges directly.
+# if this is what you want, you can (see VoronoiDelaunayExtensions.jl)
+# 1. apply scaleShiftPoints to your point set
+# 2. do the tessellation and push the scaled and shifted point set
+# 3. get the edges with the delaunayedges function and
+# 4. expand the end points of the edges with the expand function
+
+
+
 function sizehint!(t::DelaunayTessellation2D{T}, n::Int64) where T<:AbstractPoint2D
     required_total_size = 2n + 10
     required_total_size <= length(t._trigs) && return
@@ -133,6 +185,31 @@ function sizefit_at_least(t::DelaunayTessellation2D{T}, n::Int64) where T<:Abstr
     t
 end
 
+# convert the tessellation back to original scale after the tessellation
+function expand( tess::DelaunayTessellation2D{T},ranges) where T<:AbstractPoint2D
+  xmin = ranges[1]
+  ymin = ranges[3]
+  scale = max( ranges[4] - ranges[3], ranges[2] - ranges[1] ) / 0.98
+  offset = 1.01
+  for i in 1:length(tess._trigs)
+    #tess._trigs[i]._a = Point2D( ( tess._trigs[i]._a._x - offset ) * scale + xmin, ( tess._trigs[i]._a._y - offset ) * scale + ymin )
+    tess._trigs[i]._a = expand([tess._trigs[i]._a],ranges)[1]
+    tess._trigs[i]._b = expand([tess._trigs[i]._b],ranges)[1]
+    tess._trigs[i]._c = expand([tess._trigs[i]._c],ranges)[1]
+    #tess._trigs[i]._b = Point2D( ( tess._trigs[i]._b._x - offset ) * scale + xmin, ( tess._trigs[i]._b._y - offset ) * scale + ymin )
+    #tess._trigs[i]._c = Point2D( ( tess._trigs[i]._c._x - offset ) * scale + xmin, ( tess._trigs[i]._c._y - offset ) * scale + ymin )
+    tess._trigs[i]._bx   = tess._trigs[i]._bx  * scale
+    tess._trigs[i]._by   = tess._trigs[i]._by  * scale
+    tess._trigs[i]._cx   = tess._trigs[i]._cx  * scale
+    tess._trigs[i]._cy   = tess._trigs[i]._cy  * scale
+    tess._trigs[i]._px   = tess._trigs[i]._px  * scale ^ 3
+    tess._trigs[i]._py   = tess._trigs[i]._py  * scale ^ 3
+    tess._trigs[i]._pr2  = tess._trigs[i]._pr2 * scale ^ 2
+  end
+  tess._ranges = ranges
+  return tess
+end
+
 struct DelaunayEdge{T<:AbstractPoint2D}
     _a::T
     _b::T
@@ -159,12 +236,13 @@ geta(e::VoronoiEdgeWithoutGenerators) = e._a
 getb(e::VoronoiEdgeWithoutGenerators) = e._b
 
 # TODO: is an iterator faster?
-function delaunayedges(t::DelaunayTessellation2D)
+function delaunayedges( t::DelaunayTessellation2D)
+    ranges = t._ranges
     visited = zeros(Bool, t._last_trig_index)
     function delaunayiterator(c::Channel)
         @inbounds for ix in 2:t._last_trig_index
             tr = t._trigs[ix]
-            isexternal(tr) && continue
+            isexternal( tr, ranges ) && continue
             visited[ix] && continue
             visited[ix] = true
             ix_na = tr._neighbour_a
@@ -252,9 +330,10 @@ mutable struct TrigIter
     ix::Int64
 end
 
-function iterate(t::DelaunayTessellation2D, it::TrigIter=TrigIter(2)) # default it resembles old start
+function iterate( t::DelaunayTessellation2D, it::TrigIter=TrigIter(2) ) # default it resembles old start
+    ranges = t._ranges
     # resembles old done
-    while isexternal(t._trigs[it.ix]) && it.ix <= t._last_trig_index
+    while isexternal(t._trigs[it.ix], ranges) && it.ix <= t._last_trig_index
         it.ix += 1
     end
     if it.ix > t._last_trig_index
@@ -268,7 +347,10 @@ end
 
 function findindex(tess::DelaunayTessellation2D{T}, p::T) where T<:AbstractPoint2D
     i::Int64 = tess._last_trig_index
+    counter::Int64 = 0
+    ntriangles = length(tess._trigs)
     while true
+        counter += 1
         @inbounds w = intriangle(tess._trigs[i], p)
         w > 0 && return i
         @inbounds tr = tess._trigs[i]
@@ -279,6 +361,9 @@ function findindex(tess::DelaunayTessellation2D{T}, p::T) where T<:AbstractPoint
         else
             i = tr._neighbour_c
         end
+        if counter > ntriangles
+            throw(DomainError("Point was not found in any triangle"))
+        end
     end
 end
 

src/VoronoiDelaunayExtensions.jl

@@ -0,0 +1,176 @@
+# scale the point set down such that the union of all circumcircles lies within the square frame.
+# this does not only allow for arbitrary point distributions, but more importantly
+# ensures convexity (and correctness) of the resulting Delaunay triangulation.
+
+function scaleShiftPoints( points::Array{T,1} ) where T<:AbstractPoint2D
+  # convex hull of point set to find the outer edges
+  hull = quickHull( points )
+  # union of the outer circumcircles
+  ccu = circumcircleUnion( hull, points )
+  # frame around the union of circumcircles
+  ranges = frameRanges( ccu )
+  # scale the point set down to the initial square
+  scaledPoints = shrink( points, ranges )
+  # do the tessellation on the scaled point set,
+  # use the ranges to convert back to original scale afterwards
+  return scaledPoints, ranges
+end
+
+# convert the edge points back to original scale after the tessellation
+
+function expand( points::Array{Point2D,1}, ranges::NTuple{4,Float64} )
+  xmin = ranges[1]
+  ymin = ranges[3]
+  scale = max( ranges[4] - ranges[3], ranges[2] - ranges[1] ) / 0.98
+  offset = 1.01
+  scaledPoints = [ Point2D( ( p._x - offset ) * scale + xmin, ( p._y - offset ) * scale + ymin ) for p in points ]
+end
+
+function shrink( points::Array{Point2D,1}, ranges::NTuple{4,Float64} )
+  h = ranges[4] - ranges[3]
+  b = ranges[2] - ranges[1]
+  offset = 1.01
+  scale = 0.98 / max( h, b )
+  scaledPoints = [ Point2D( ( p._x - ranges[1] ) * scale + offset, ( p._y - ranges[3] ) * scale + offset ) for p in points ]
+end
+
+
+
+#=== auxiliary functions and structs ===#
+
+# convex hull of point set
+
+function quickHull(points::Array{T,1}) where T<:AbstractPoint2D
+    return quickHull(Point2D.(points))
+end
+
+function quickHull( points::Array{Point2D,1} )
+  ConvexHull = Array{Point2D,1}(undef,0)
+  A = points[ argmin( getx.( points ) ) ]   # leftmost point
+  B = points[ argmax( getx.( points ) ) ]   # rightmost point
+  push!( ConvexHull, A )
+  push!( ConvexHull, B )
+  pr = Array{Point2D,1}(undef,0)            # points to the right of AB
+  pl = Array{Point2D,1}(undef,0)            # points to the left of AB
+  ( pl, pr ) = dividePointSet( Line2D( A, B ), setdiff( points, ConvexHull ) )
+  findHull!( ConvexHull, pr, A, B )         # divide-and-conquer approach
+  findHull!( ConvexHull, pl, B, A )
+  return ConvexHull
+end
+
+function findHull!( ConvexHull::Array{T,1}, points::Array{T,1}, A::T, B::T ) where T<:AbstractPoint2D
+  if isempty( points )
+    return
+  end
+  C = findFarthestPoint( A, B, points )
+  pos = findfirst( x -> x == A, ConvexHull ) + 1
+  insert!( ConvexHull, pos, C )
+  pAC = dividePointSet( Line2D(A,C), setdiff(points,[C]) )[2]
+  pCB = dividePointSet( Line2D(C,B), setdiff(points,[C]) )[2]
+  findHull!( ConvexHull, pAC, A, C )
+  findHull!( ConvexHull, pCB, C, B )
+end
+
+function dividePointSet( l::Line2D, points::Array{T,1} ) where T <: AbstractPoint2D
+  pr = Array{T,1}(undef,0)
+  pl = Array{T,1}(undef,0)
+  for i in 1:length(points)
+    if orientation(l,points[i]) == -1
+      push!(pr,points[i])
+    else
+      push!(pl,points[i])
+    end
+  end
+  return (pl, pr)
+end
+
+function findFarthestPoint(a::T,b::T,points::Array{T,1}) where T<:AbstractPoint2D
+  distances = Array{Float64,1}(undef,size(points,1))
+  for i in 1:length(distances)
+    distances[i] = pointLineDistance(a,b,points[i])
+  end
+  return points[argmax(distances)]
+end
+
+function pointLineDistance(a::T,b::T,c::T) where T<:AbstractPoint2D
+  return abs( (gety(b)-gety(a))*getx(c) - (getx(b)-getx(a))*gety(c) + getx(b)*gety(a) - gety(b)*getx(a) ) / sqrt( (gety(b)-gety(a))^2 + (getx(b)-getx(a))^2 )
+end
+
+
+
+mutable struct Circle{T<:AbstractPoint2D}
+  c::T
+  r::Float64
+end
+
+function getcen(circ::Circle{T}) where T<:AbstractPoint2D
+    return circ.c
+end
+
+function getcenx(circ::Circle{T}) where T<:AbstractPoint2D
+    return getx( getcen( circ ) )
+end
+function getceny(circ::Circle{T}) where T<:AbstractPoint2D
+    return gety( getcen( circ ) )
+end
+function getrad(circ::Circle{T}) where T<: AbstractPoint2D
+    return circ.r
+end
+
+function circumcircleUnion( hull::Array{T,1}, points::Array{S,1} ) where {T<:AbstractPoint2D,S<:AbstractPoint2D}
+  ccU = Array{Circle{Point2D},1}(undef,length(hull))
+  hullcirc = copy(hull)
+  push!( hullcirc, hull[1] )
+  for i in 1:length( hull )
+    cc2 = circumcircle( hullcirc[i], hullcirc[i+1] )
+    ccU[i] = cc2
+    for j in 1:length(points)
+      if pointsDistance( points[j], cc2.c ) < cc2.r
+        cc3 = circumcircle( hullcirc[i], hullcirc[i+1], points[j] )
+        if getrad(cc3) > getrad(ccU[i])
+          ccU[i] = cc3
+        end
+      end
+    end
+  end
+  return ccU
+end
+
+function pointsDistance( p1::T, p2::S ) where {T<:AbstractPoint2D, S<:AbstractPoint2D}
+  sqrt( ( getx(p1) - getx(p2) ) ^ 2 + ( gety(p1) - gety(p2) ) ^ 2 )
+end
+
+function circumcircle( a::T, b::S ) where {T<:AbstractPoint2D, S<:AbstractPoint2D}
+  px = (getx(a) + getx(b)) / 2.0
+  py = (gety(a) + gety(b)) / 2.0
+  dx = (getx(b) - getx(a)) / 2.0
+  dy = (gety(b) - gety(a)) / 2.0
+  r = sqrt( dx^2 + dy^2 )
+  return Circle{Point2D}( Point2D(px,py), r )
+end
+
+function circumcircle( a::T, b::S, c::U ) where {T<:AbstractPoint2D,S<:AbstractPoint2D, U<:AbstractPoint2D}
+  la = getx(a)^2 + gety(a)^2
+  lb = getx(b)^2 + gety(b)^2
+  lc = getx(c)^2 + gety(c)^2
+  xyy = getx(a) * (gety(b) - gety(c))
+  yxx = gety(a) * (getx(b) - getx(c))
+  xy = getx(b) * gety(c)
+  yx = gety(b) * getx(c)
+  z = 2 * ( xyy - yxx + xy - yx )
+  px = ( la*(gety(b)-gety(c)) + lb*(gety(c)-gety(a)) + lc*(gety(a)-gety(b)) ) / z
+  py = ( la*(getx(c)-getx(b)) + lb*(getx(a)-getx(c)) + lc*(getx(b)-getx(a)) ) / z
+  r = sqrt( (px-getx(a))^2 + (py-gety(a))^2 )
+  return Circle{Point2D}( Point2D(px,py), r )
+end
+
+
+
+
+function frameRanges( ccU::Array{Circle{T},1} ) where T<:AbstractPoint2D
+  xmin = minimum( getcenx.(ccU) .- getrad.(ccU) )
+  xmax = maximum( getcenx.(ccU) .+ getrad.(ccU) )
+  ymin = minimum( getceny.(ccU) .- getrad.(ccU) )
+  ymax = maximum( getceny.(ccU) .+ getrad.(ccU) )
+  return xmin, xmax, ymin, ymax
+end