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,10 +107,10 @@ 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}
@@ -108,6 +137,25 @@ 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{Point2D,1} )
+  scaledPoints, ranges = scaleShiftPoints( points )
+  scaledTess = DelaunayTessellation( length( points ) )
+  push!( scaledTess, scaledPoints )  
+  tess = expand( scaledTess, ranges )
+  return tess, ranges
+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 +181,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::NTuple{4,Float64} ) where T<:AbstractPoint2D
+  scaledTess = deepcopy(tess)
+  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)
+    scaledTess._trigs[i]._a._x = ( tess._trigs[i]._a._x - offset ) * scale + xmin
+    scaledTess._trigs[i]._a._y = ( tess._trigs[i]._a._y - offset ) * scale + ymin
+    scaledTess._trigs[i]._b._x = ( tess._trigs[i]._b._x - offset ) * scale + xmin
+    scaledTess._trigs[i]._b._y = ( tess._trigs[i]._b._y - offset ) * scale + ymin
+    scaledTess._trigs[i]._c._x = ( tess._trigs[i]._c._x - offset ) * scale + xmin
+    scaledTess._trigs[i]._c._y = ( tess._trigs[i]._c._y - offset ) * scale + ymin
+    scaledTess._trigs[i]._bx   = tess._trigs[i]._bx  * scale
+    scaledTess._trigs[i]._by   = tess._trigs[i]._by  * scale
+    scaledTess._trigs[i]._cx   = tess._trigs[i]._cx  * scale
+    scaledTess._trigs[i]._cy   = tess._trigs[i]._cy  * scale
+    scaledTess._trigs[i]._px   = tess._trigs[i]._px  * scale ^ 3
+    scaledTess._trigs[i]._py   = tess._trigs[i]._py  * scale ^ 3
+    scaledTess._trigs[i]._pr2  = tess._trigs[i]._pr2 * scale ^ 2
+  end
+  return scaledTess
+end
+
 struct DelaunayEdge{T<:AbstractPoint2D}
     _a::T
     _b::T
@@ -159,12 +232,12 @@ geta(e::VoronoiEdgeWithoutGenerators) = e._a
 getb(e::VoronoiEdgeWithoutGenerators) = e._b
 
 # TODO: is an iterator faster?
-function delaunayedges(t::DelaunayTessellation2D)
+function delaunayedges( t::DelaunayTessellation2D, ranges::NTuple{4,Float64}=(min_coord,max_coord,min_coord,max_coord) )
     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 +325,9 @@ 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), ranges::NTuple{4,Float64}=(min_coord,max_coord,min_coord,max_coord) ) # default it resembles old start
     # 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

src/VoronoiDelaunayExtensions.jl

@@ -0,0 +1,174 @@
+# 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{Point2D,1} )
+  # 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} )
+  scaledPoints = deepcopy(points)
+  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(scaledPoints)
+    scaledPoints[i]._x = ( points[i]._x - offset ) * scale + xmin
+    scaledPoints[i]._y = ( points[i]._y - offset ) * scale + ymin
+  end
+  return scaledPoints
+end
+
+
+
+#=== auxiliary functions and structs ===#
+
+# convex hull of point set
+
+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{Point2D,1}, points::Array{Point2D,1}, A::Point2D, B::Point2D )
+  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{Point2D,1} )
+  pr = Array{Point2D,1}(undef,0)
+  pl = Array{Point2D,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::Point2D,b::Point2D,points::Array{Point2D,1})
+  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::Point2D,b::Point2D,c::Point2D)
+  abs( (b._y-a._y)*c._x - (b._x-a._x)*c._y + b._x*a._y - b._y*a._x ) / sqrt( (b._y-a._y)^2 + (b._x-a._x)^2 )
+end
+
+
+
+mutable struct Circle
+  c::Point2D
+  r::Float64
+end
+
+function circumcircleUnion( hull::Array{Point2D,1}, points::Array{Point2D,1} )
+  ccU = Array{Circle,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 cc3.r > ccU[i].r
+          ccU[i] = cc3
+        end
+      end
+    end
+  end
+  return ccU
+end
+
+function pointsDistance( p1::Point2D, p2::Point2D )
+  sqrt( ( p1._x - p2._x ) ^ 2 + ( p1._y - p2._y ) ^ 2 )
+end
+
+function circumcircle( a::Point2D, b::Point2D )
+  px = (a._x + b._x) / 2.0
+  py = (a._y + b._y) / 2.0
+  dx = (b._x - a._x) / 2.0
+  dy = (b._y - a._y) / 2.0
+  r = sqrt( dx^2 + dy^2 )
+  return Circle( Point2D(px,py), r )
+end
+
+function circumcircle( a::Point2D, b::Point2D, c::Point2D )
+  la = a._x^2 + a._y^2
+  lb = b._x^2 + b._y^2
+  lc = c._x^2 + c._y^2
+  xyy = a._x * (b._y - c._y)
+  yxx = a._y * (b._x - c._x)
+  xy = b._x * c._y
+  yx = b._y * c._x
+  z = 2 * ( xyy - yxx + xy - yx )
+  px = ( la*(b._y-c._y) + lb*(c._y-a._y) + lc*(a._y-b._y) ) / z
+  py = ( la*(c._x-b._x) + lb*(a._x-c._x) + lc*(b._x-a._x) ) / z
+  r = sqrt( (px-a._x)^2 + (py-a._y)^2 )
+  return Circle( Point2D(px,py), r )
+end
+
+getcen(circ::Circle)  = circ.c
+getcenx(circ::Circle) = getx( getcen( circ ) )
+getceny(circ::Circle) = gety( getcen( circ ) )
+getrad(circ::Circle)  = circ.r
+
+
+
+function frameRanges( ccU::Array{Circle,1} )
+  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
+
+
+
+function shrink( points::Array{Point2D,1}, ranges::NTuple{4,Float64} )
+  scaledPoints = deepcopy(points)
+  h = ranges[4] - ranges[3]
+  b = ranges[2] - ranges[1]
+  offset = 1.01
+  scale = 0.98 / max( h, b )
+  for i in 1:length(points)
+    scaledPoints[i]._x = ( points[i]._x - ranges[1] ) * scale + offset
+    scaledPoints[i]._y = ( points[i]._y - ranges[3] ) * scale + offset
+  end
+  return scaledPoints
+end