VoronoiDelaunay_v0.2.9

October 7, 2014 ยท View on GitHub

module VoronoiDelaunay

VoronoiDelaunay_v0.2.9

Fast, robust 2D Voronoi/Delaunay tesselation

Implementation follows algorithms described in http://arxiv.org/abs/0901.4107

and used (for e.g.) in the Illustris Simulation

http://www.illustris-project.org/

Author: Ariel Keselman (skariel@gmail.com)

License: MIT

Bug reportss welcome!

Usage notes:

-------------

- TODO!

Todo

------------------

- Make a module out of this

- Make a package out of this

- More comments in code

- Documentation

export DTess2D, min_coord, max_coord, delaunayedges, voronoiedges, length, getindex, start, next, done, findindex!, push!, pushunsorted!, Point2D, AbstractPoint2D

using GeometricalPredicates import GeometricalPredicates

import Base:push!, length, start, done, next

const min_coord = GeometricalPredicates.min_coord + eps(Float64) const max_coord = GeometricalPredicates.max_coord - eps(Float64)

type VorTrig{T<:AbstractPoint2D} <: AbstractNegativelyOrientedTriangle _a::T; _b::T; _c::T _bx::Float64; _by::Float64 _cx::Float64; _cy::Float64 _px::Float64; _py::Float64 _pr2::Float64 _sz::Float32 _nb_a::Int64 _nb_b::Int64 _nb_c::Int64 function VorTrig(pa::T, pb::T, pc::T, na::Int64, nb::Int64, nc::Int64) t = new(pa, pb, pc, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, na, nb, nc) clean!(t) t end end VorTrig{T<:Point2D}(pa::T, pb::T, pc::T, na::Int64, nb::Int64, nc::Int64) = VorTrig{T}(pa, pb, pc, na, nb, nc)

geta(t::VorTrig) = t._a getb(t::VorTrig) = t._b getc(t::VorTrig) = t._c

seta(t::VorTrig, p::AbstractPoint) = t._a=p setb(t::VorTrig, p::AbstractPoint) = t._b=p setc(t::VorTrig, p::AbstractPoint) = t._c=p setpoints(t::VorTrig, pa::AbstractPoint, pb::AbstractPoint, pc::AbstractPoint) = (t._a=pa; t._b=pb; t._c=pc)

type DTess2D{T<:AbstractPoint2D} _trigs::Array{VorTrig{T}, 1} _last_trig_index::Int64 _edges_to_check::Array{Int64, 1} function DTess2D(n::Int64 = 2) const a = T(GeometricalPredicates.min_coord, GeometricalPredicates.min_coord) const b = T(GeometricalPredicates.min_coord, GeometricalPredicates.max_coord) const c = T(GeometricalPredicates.max_coord, GeometricalPredicates.min_coord) const d = T(GeometricalPredicates.max_coord, GeometricalPredicates.max_coord) const _trigs = VorTrig{T}[VorTrig{T}(d,c,b, 2,1,1), VorTrig{T}(a,b,c, 3,1,1), VorTrig{T}(d,c,b, 2,1,1)]

	for i in [1:2*n+10]
		push!(_trigs, VorTrig{T}(a,b,c, 1,1,1))
	end

	const _edges_to_check = Int64[]
	sizehint(_edges_to_check, 1000)
	new(_trigs, 3, _edges_to_check)
end

end DTess2D(n::Int64) = DTess2D{Point2D}(n) DTess2D{T<:AbstractPoint2D}(n::Int64, ::T) = DTess2D{T}(n)

immutable Edge{T<:AbstractPoint2D} a::T b::T end

function delaunayedges(t::DTess2D) visited = zeros(Bool, t._last_trig_index) function delaunayiterator() for ix in 2:t._last_trig_index visited[ix] && continue visited[ix] = true const tr = t._trigs[ix] const ix_na = tr._nb_a if !visited[ix_na] produce(Edge(getb(tr), getc(tr))) end const ix_nb = tr._nb_b if !visited[ix_nb] produce(Edge(geta(tr), getc(tr))) end const ix_nc = tr._nb_c if !visited[ix_nc] produce(Edge(geta(tr), getb(tr))) end end end Task(delaunayiterator) end

function voronoiedges(t::DTess2D) visited = zeros(Bool, t._last_trig_index) visited[1] = true function voronoiiterator() for ix in 2:t._last_trig_index visited[ix] && continue visited[ix] = true const tr = t._trigs[ix] const cc = circumcenter(tr)

		const ix_na = tr._nb_a
		if !visited[ix_na] 
			produce(Edge(cc, circumcenter(t._trigs[ix_na])))
		end
		const ix_nb = tr._nb_b
		if !visited[ix_nb] 
			produce(Edge(cc, circumcenter(t._trigs[ix_nb])))
		end
		const ix_nc = tr._nb_c
		if !visited[ix_nc] 
			produce(Edge(cc, circumcenter(t._trigs[ix_nc])))
		end
	end
end
Task(voronoiiterator)

end

length(tess::DTess2D) = tess._last_trig_index-1 getindex(tess::DTess2D, n::Int64) = tess._trigs[n+1]

start(t::DTess2D) = 2 next(t::DTess2D, it::Int64) = (t._trigs[it], it+1) done(t::DTess2D, it::Int64) = it > t._last_trig_index

function findindex!{T<:AbstractPoint2D}(tess::DTess2D{T}, p::T) i::Int64 = tess._last_trig_index while true @inbounds const w = intriangle(tess._trigs[i], p) w > 0 && return i @inbounds const tr = tess._trigs[i] i = w==-1? tr._nb_a : w==-2? tr._nb_b : tr._nb_c @assert i!=1 prinln("p=",p) end end

function _pushunfixed!{T<:AbstractPoint2D}(tess::DTess2D{T}, p::T) i = findindex!(tess, p) const ltrigs1::Int64 = tess._last_trig_index+1 const ltrigs2::Int64 = tess._last_trig_index+2

@inbounds const t1 = tess._trigs[i]
old_t1_a = geta(t1)
seta(t1, p)
old_t1_b = t1._nb_b
old_t1_c = t1._nb_c
t1._nb_b = ltrigs1
t1._nb_c = ltrigs2
clean!(t1)

@inbounds const t2 = tess._trigs[ltrigs1]
setpoints(t2, old_t1_a, p, getc(t1))
t2._nb_a = i
t2._nb_b = old_t1_b
t2._nb_c = ltrigs2
clean!(t2)

@inbounds const t3 = tess._trigs[ltrigs2]
setpoints(t3, old_t1_a, getb(t1), p)
t3._nb_a = i
t3._nb_b = ltrigs1
t3._nb_c = old_t1_c
clean!(t3)

@inbounds const nt2 = tess._trigs[t2._nb_b]
if nt2._nb_a==i
	nt2._nb_a = ltrigs1
elseif nt2._nb_b==i
	nt2._nb_b = ltrigs1
else
	nt2._nb_c = ltrigs1
end

@inbounds const nt3 = tess._trigs[t3._nb_c]
if nt3._nb_a==i
	nt3._nb_a = ltrigs2
elseif nt3._nb_b==i
	nt3._nb_b = ltrigs2
else
	nt3._nb_c = ltrigs2
end

tess._last_trig_index += 2

i

end

function _flipa!(tess::DTess2D, ix1::Int64, ix2::Int64) @inbounds const ot1 = tess._trigs[ix1] @inbounds const ot2 = tess._trigs[ix2] if ot2._nb_a == ix1 _flipaa!(tess, ix1, ix2, ot1, ot2) elseif ot2._nb_b == ix1 _flipab!(tess, ix1, ix2, ot1, ot2) else _flipac!(tess, ix1, ix2, ot1, ot2) end end

function _endflipa!(tess::DTess2D, ix1::Int64, ix2::Int64, ot1::VorTrig, ot2::VorTrig) @inbounds const n1 = tess._trigs[ot1._nb_a] if n1._nb_a==ix2 n1._nb_a = ix1 elseif n1._nb_b==ix2 n1._nb_b = ix1 else n1._nb_c = ix1 end @inbounds const n2 = tess._trigs[ot2._nb_c] if n2._nb_a==ix1 n2._nb_a = ix2 elseif n2._nb_b==ix1 n2._nb_b = ix2 else n2._nb_c = ix2 end end

function _flipaa!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64, ot1::VorTrig{T}, ot2::VorTrig{T}) old_ot1_geom_b = getb(ot1) setb(ot1, geta(ot2)) ot1._nb_a = ot2._nb_c old_ot1_nb_c = ot1._nb_c ot1._nb_c = ix2 clean!(ot1)

setc(ot2, geta(ot2))
seta(ot2, geta(ot1))
setb(ot2, old_ot1_geom_b)
ot2._nb_a = ot2._nb_b
ot2._nb_b = ix1
ot2._nb_c = old_ot1_nb_c	
clean!(ot2)

_endflipa!(tess,ix1,ix2,ot1,ot2)

end

function _flipab!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64, ot1::VorTrig{T}, ot2::VorTrig{T}) old_ot1_geom_b = getb(ot1)

setb(ot1, getb(ot2))
ot1._nb_a = ot2._nb_a
old_ot1_nb_c = ot1._nb_c
ot1._nb_c = ix2
clean!(ot1)

setc(ot2, getb(ot2))
seta(ot2, geta(ot1))
setb(ot2, old_ot1_geom_b)
ot2._nb_a = ot2._nb_c
ot2._nb_b = ix1
ot2._nb_c = old_ot1_nb_c
clean!(ot2)

_endflipa!(tess,ix1,ix2,ot1,ot2)

end

function _flipac!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64, ot1::VorTrig{T}, ot2::VorTrig{T}) old_ot1_geom_b = getb(ot1) setb(ot1, getc(ot2)) ot1._nb_a = ot2._nb_b old_ot1_nb_c = ot1._nb_c ot1._nb_c = ix2 clean!(ot1)

setc(ot2, getc(ot2))
seta(ot2,  geta(ot1))
setb(ot2, old_ot1_geom_b)
ot2._nb_b = ix1
ot2._nb_c = old_ot1_nb_c
clean!(ot2)

_endflipa!(tess,ix1,ix2,ot1,ot2)

end

########################

function _flipb!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64) @inbounds const ot1 = tess._trigs[ix1] @inbounds const ot2 = tess._trigs[ix2] if ot2._nb_a == ix1 _flipba!(tess, ix1, ix2, ot1, ot2) elseif ot2._nb_b == ix1 _flipbb!(tess, ix1, ix2, ot1, ot2) else _flipbc!(tess, ix1, ix2, ot1, ot2) end end

function _endflipb!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64, ot1::VorTrig{T}, ot2::VorTrig{T}) @inbounds const n1 = tess._trigs[ot1._nb_b] if n1._nb_a==ix2 n1._nb_a = ix1 elseif n1._nb_b==ix2 n1._nb_b = ix1 else n1._nb_c = ix1 end @inbounds const n2 = tess._trigs[ot2._nb_a] if n2._nb_a==ix1 n2._nb_a = ix2 elseif n2._nb_b==ix1 n2._nb_b = ix2 else n2._nb_c = ix2 end end

function _flipba!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64, ot1::VorTrig{T}, ot2::VorTrig{T}) old_ot1_geom_c = getc(ot1) setc(ot1, geta(ot2)) old_ot1_nb_a = ot1._nb_a ot1._nb_a = ix2 ot1._nb_b = ot2._nb_c clean!(ot1)

setb(ot2, getb(ot1))
setc(ot2, old_ot1_geom_c)
ot2._nb_a = old_ot1_nb_a
ot2._nb_c = ix1
clean!(ot2)

_endflipb!(tess,ix1,ix2,ot1,ot2)

end

function _flipbb!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64, ot1::VorTrig{T}, ot2::VorTrig{T}) old_ot1_geom_c = getc(ot1) setc(ot1, getb(ot2)) old_ot1_nb_a = ot1._nb_a ot1._nb_a = ix2 ot1._nb_b = ot2._nb_a clean!(ot1)

seta(ot2, getb(ot2))
setb(ot2, getb(ot1))
setc(ot2, old_ot1_geom_c)
ot2._nb_a = old_ot1_nb_a
ot2._nb_b = ot2._nb_c
ot2._nb_c = ix1
clean!(ot2)

_endflipb!(tess,ix1,ix2,ot1,ot2)

end

function _flipbc!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64, ot1::VorTrig{T}, ot2::VorTrig{T}) old_ot1_geom_c = getc(ot1) setc(ot1, getc(ot2)) old_ot1_nb_a = ot1._nb_a ot1._nb_a = ix2 ot1._nb_b = ot2._nb_b clean!(ot1)

seta(ot2, getc(ot2))
setb(ot2, getb(ot1))
setc(ot2, old_ot1_geom_c)
ot2._nb_b = ot2._nb_a
ot2._nb_a = old_ot1_nb_a
ot2._nb_c = ix1
clean!(ot2)

_endflipb!(tess,ix1,ix2,ot1,ot2)

end

########################

function _flipc!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64) @inbounds const ot1 = tess._trigs[ix1] @inbounds const ot2 = tess._trigs[ix2] if ot2._nb_a == ix1 _flipca!(tess, ix1, ix2, ot1, ot2) elseif ot2._nb_b == ix1 _flipcb!(tess, ix1, ix2, ot1, ot2) else _flipcc!(tess, ix1, ix2, ot1, ot2) end end

function _endflipc!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64, ot1::VorTrig{T}, ot2::VorTrig{T}) @inbounds const n1 = tess._trigs[ot1._nb_c] if n1._nb_a==ix2 n1._nb_a = ix1 elseif n1._nb_b==ix2 n1._nb_b = ix1 else n1._nb_c = ix1 end @inbounds const n2 = tess._trigs[ot2._nb_b] if n2._nb_a==ix1 n2._nb_a = ix2 elseif n2._nb_b==ix1 n2._nb_b = ix2 else n2._nb_c = ix2 end end

function _flipca!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64, ot1::VorTrig{T}, ot2::VorTrig{T}) old_ot1_geom_a = geta(ot1) seta(ot1, geta(ot2)) old_ot1_nb_b = ot1._nb_b ot1._nb_b = ix2 ot1._nb_c = ot2._nb_c clean!(ot1)

setb(ot2, geta(ot2))
seta(ot2, old_ot1_geom_a)
setc(ot2, getc(ot1))
ot2._nb_a = ix1
ot2._nb_c = ot2._nb_b
ot2._nb_b = old_ot1_nb_b
clean!(ot2)

_endflipc!(tess,ix1,ix2,ot1,ot2)

end

function _flipcb!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64, ot1::VorTrig{T}, ot2::VorTrig{T}) old_ot1_geom_a = geta(ot1) seta(ot1, getb(ot2)) old_ot1_nb_b = ot1._nb_b ot1._nb_b = ix2 ot1._nb_c = ot2._nb_a clean!(ot1)

seta(ot2, old_ot1_geom_a)
setc(ot2, getc(ot1))
ot2._nb_a = ix1
ot2._nb_b = old_ot1_nb_b
clean!(ot2)

_endflipc!(tess,ix1,ix2,ot1,ot2)

end

function _flipcc!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix1::Int64, ix2::Int64, ot1::VorTrig{T}, ot2::VorTrig{T}) old_ot1_geom_a = geta(ot1) seta(ot1, getc(ot2)) old_ot1_nb_b = ot1._nb_b ot1._nb_b = ix2 ot1._nb_c = ot2._nb_b clean!(ot1)

seta(ot2, old_ot1_geom_a)
setb(ot2, getc(ot2))
setc(ot2, getc(ot1))
ot2._nb_c = ot2._nb_a
ot2._nb_a = ix1
ot2._nb_b = old_ot1_nb_b
clean!(ot2)

_endflipc!(tess,ix1,ix2,ot1,ot2)

end

################################################################################

function _fix!{T<:AbstractPoint2D}(tess::DTess2D{T}, ix_trig::Int64) @inbounds const center_pt = geta(tess._trigs[ix_trig])

# `A` - edge
push!(tess._edges_to_check, ix_trig)
while length(tess._edges_to_check) > 0
	@inbounds const trix = tess._edges_to_check[end]
	@inbounds const tr_i = tess._trigs[trix]
	@inbounds const nb_a = tr_i._nb_a
	if nb_a > 1
		@inbounds const tr_f = tess._trigs[nb_a]
		if incircle(tr_f, center_pt) > 0
			_flipa!(tess, trix, nb_a)
			push!(tess._edges_to_check, nb_a)
			continue
		end
	end
	pop!(tess._edges_to_check)
end

# `B` - edge
push!(tess._edges_to_check, tess._last_trig_index-1)
while length(tess._edges_to_check) > 0
	@inbounds const trix = tess._edges_to_check[end]
	@inbounds const tr_i = tess._trigs[trix]
	@inbounds const nb_b = tr_i._nb_b
	if nb_b > 1
		@inbounds const tr_f = tess._trigs[nb_b]
		if incircle(tr_f, center_pt) > 0
			_flipb!(tess, trix, nb_b)
			push!(tess._edges_to_check, nb_b)
			continue
		end
	end
	pop!(tess._edges_to_check)
end

# `C` - edge
push!(tess._edges_to_check, tess._last_trig_index)
while length(tess._edges_to_check) > 0
	@inbounds const trix = tess._edges_to_check[end]
	@inbounds const tr_i = tess._trigs[trix]
	@inbounds const nb_c = tr_i._nb_c
	if nb_c > 1
		@inbounds const tr_f = tess._trigs[nb_c]
		if incircle(tr_f, center_pt) > 0
			_flipc!(tess, trix, nb_c)
			push!(tess._edges_to_check, nb_c)
			continue
		end
	end
	pop!(tess._edges_to_check)
end

end

function pushunfixed!{T<:AbstractPoint2D}(tess::DTess2D{T}, a::Array{T, 1}, n::Int64) order = sortperm([peanokey(p,n) for p in a]) for i in order @inbounds _pushunfixed!(tess, a[i]) end end

function push!{T<:AbstractPoint2D}(tess::DTess2D{T}, p::T) const i = _pushunfixed!(tess, p) _fix!(tess, i) end

function pushunsorted!{T<:AbstractPoint2D}(tess::DTess2D{T}, a::Array{T, 1}) for p in a push!(tess, p) end end

function push!{T<:AbstractPoint2D}(tess::DTess2D{T}, a::Array{T, 1}) shuffle!(a) mssort!(a) pushunsorted!(tess, a) end

using Base.Test

function test() tess = DTess2D(100) @test findindex!(tess, Point2D(1.1, 1.1)) == 2 @test findindex!(tess, Point2D(1.9, 1.9)) == 3 @test findindex!(tess, Point2D(1.9, 1.9)) == 3

tess = DTess2D(100)
@test findindex!(tess, Point2D(1.9, 1.9)) == 3
@test findindex!(tess, Point2D(1.1, 1.1)) == 2
@test findindex!(tess, Point2D(1.1, 1.1)) == 2

@test tess._trigs[2]._nb_a == 3

_pushunfixed!(tess, Point2D(1.1, 1.1))

@test tess._trigs[2]._nb_a == 3
@test tess._trigs[2]._nb_b == 4
@test tess._trigs[2]._nb_c == 5

@test tess._trigs[3]._nb_a == 2
@test tess._trigs[3]._nb_b == 1
@test tess._trigs[3]._nb_c == 1

@test tess._trigs[4]._nb_a == 2
@test tess._trigs[4]._nb_b == 1
@test tess._trigs[4]._nb_c == 5

@test tess._trigs[5]._nb_a == 2
@test tess._trigs[5]._nb_b == 4
@test tess._trigs[5]._nb_c == 1

pa = Point2D(GeometricalPredicates.min_coord, GeometricalPredicates.min_coord)
pb = Point2D(GeometricalPredicates.min_coord, GeometricalPredicates.max_coord)
pc = Point2D(GeometricalPredicates.max_coord, GeometricalPredicates.min_coord)
pd = Point2D(GeometricalPredicates.max_coord, GeometricalPredicates.max_coord)
pp = Point2D(1.1,1.1)

@test getc(tess._trigs[5]) == pp
@test geta(tess._trigs[2]) == pp
@test getb(tess._trigs[4]) == pp

@test getb(tess._trigs[5]) == pb
@test getb(tess._trigs[2]) == pb
@test getc(tess._trigs[3]) == pb

@test getc(tess._trigs[4]) == pc
@test getc(tess._trigs[2]) == pc
@test getb(tess._trigs[3]) == pc

@test geta(tess._trigs[5]) == pa
@test geta(tess._trigs[4]) == pa

@test geta(tess._trigs[3]) == pd

@test findindex!(tess, Point2D(1.01, 1.1)) == 5
@test findindex!(tess, Point2D(1.1, 1.01)) == 4

@test findindex!(tess, Point2D(1.11, 1.11)) == 2
@test findindex!(tess, Point2D(1.6, 1.6)) == 3
@test findindex!(tess, Point2D(1.11, 1.1101)) == 2
@test findindex!(tess, Point2D(1.6, 1.601)) == 3
@test findindex!(tess, Point2D(1.11, 1.11)) == 2
@test findindex!(tess, Point2D(1.6, 1.6)) == 3

p2 = Point2D(1.9,1.9)
_pushunfixed!(tess, p2)

@test geta(tess._trigs[7]) == pd
@test getb(tess._trigs[7]) == pc
@test getc(tess._trigs[7]) == p2

@test tess._trigs[7]._nb_a == 3
@test tess._trigs[7]._nb_b == 6
@test tess._trigs[7]._nb_c == 1

@test geta(tess._trigs[6]) == pd
@test getb(tess._trigs[6]) == p2
@test getc(tess._trigs[6]) == pb

@test tess._trigs[6]._nb_a == 3
@test tess._trigs[6]._nb_b == 1
@test tess._trigs[6]._nb_c == 7

@test geta(tess._trigs[3]) == p2
@test getb(tess._trigs[3]) == pc
@test getc(tess._trigs[3]) == pb

@test tess._trigs[3]._nb_a == 2
@test tess._trigs[3]._nb_b == 6
@test tess._trigs[3]._nb_c == 7

@test geta(tess._trigs[5]) == pa
@test getb(tess._trigs[5]) == pb
@test getc(tess._trigs[5]) == pp

@test tess._trigs[5]._nb_a == 2
@test tess._trigs[5]._nb_b == 4
@test tess._trigs[5]._nb_c == 1

@test geta(tess._trigs[4]) == pa
@test getb(tess._trigs[4]) == pp
@test getc(tess._trigs[4]) == pc

@test tess._trigs[4]._nb_a == 2
@test tess._trigs[4]._nb_b == 1
@test tess._trigs[4]._nb_c == 5

@test geta(tess._trigs[2]) == pp
@test getb(tess._trigs[2]) == pb
@test getc(tess._trigs[2]) == pc

@test tess._trigs[2]._nb_a == 3
@test tess._trigs[2]._nb_b == 4
@test tess._trigs[2]._nb_c == 5

@test tess._last_trig_index == 7


_flipa!(tess, 2, 3)

@test geta(tess._trigs[2]) == pp
@test getb(tess._trigs[2]) == p2
@test getc(tess._trigs[2]) == pc

@test tess._trigs[2]._nb_a == 7
@test tess._trigs[2]._nb_b == 4
@test tess._trigs[2]._nb_c == 3

@test geta(tess._trigs[3]) == pp
@test getb(tess._trigs[3]) == pb
@test getc(tess._trigs[3]) == p2

@test tess._trigs[3]._nb_a == 6
@test tess._trigs[3]._nb_b == 2
@test tess._trigs[3]._nb_c == 5

@test geta(tess._trigs[5]) == pa
@test getb(tess._trigs[5]) == pb
@test getc(tess._trigs[5]) == pp

@test geta(tess._trigs[6]) == pd
@test getb(tess._trigs[6]) == p2
@test getc(tess._trigs[6]) == pb

@test geta(tess._trigs[4]) == pa
@test getb(tess._trigs[4]) == pp
@test getc(tess._trigs[4]) == pc

@test geta(tess._trigs[7]) == pd
@test getb(tess._trigs[7]) == pc
@test getc(tess._trigs[7]) == p2

@test tess._trigs[4]._nb_a == 2
@test tess._trigs[4]._nb_b == 1
@test tess._trigs[4]._nb_c == 5

@test tess._trigs[5]._nb_a == 3
@test tess._trigs[5]._nb_b == 4
@test tess._trigs[5]._nb_c == 1

@test tess._trigs[6]._nb_a == 3
@test tess._trigs[6]._nb_b == 1
@test tess._trigs[6]._nb_c == 7

@test tess._trigs[7]._nb_a == 2
@test tess._trigs[7]._nb_b == 6
@test tess._trigs[7]._nb_c == 1

_flipb!(tess, 3, 2)

@test geta(tess._trigs[3]) == pp
@test getb(tess._trigs[3]) == pb
@test getc(tess._trigs[3]) == pc

@test geta(tess._trigs[2]) == pc
@test getb(tess._trigs[2]) == pb
@test getc(tess._trigs[2]) == p2

@test tess._trigs[3]._nb_a == 2
@test tess._trigs[3]._nb_b == 4
@test tess._trigs[3]._nb_c == 5

@test tess._trigs[2]._nb_a == 6
@test tess._trigs[2]._nb_b == 7
@test tess._trigs[2]._nb_c == 3

@test tess._trigs[6]._nb_a == 2
@test tess._trigs[7]._nb_a == 2
@test tess._trigs[5]._nb_a == 3
@test tess._trigs[4]._nb_a == 3

_flipc!(tess,2,3)

@test geta(tess._trigs[2]) == pp
@test getb(tess._trigs[2]) == pb
@test getc(tess._trigs[2]) == p2

@test geta(tess._trigs[3]) == pc
@test getb(tess._trigs[3]) == pp
@test getc(tess._trigs[3]) == p2

@test tess._trigs[2]._nb_a == 6
@test tess._trigs[2]._nb_b == 3
@test tess._trigs[2]._nb_c == 5

@test tess._trigs[3]._nb_a == 2
@test tess._trigs[3]._nb_b == 7
@test tess._trigs[3]._nb_c == 4

@test tess._trigs[6]._nb_a == 2
@test tess._trigs[7]._nb_a == 3
@test tess._trigs[5]._nb_a == 2
@test tess._trigs[4]._nb_a == 3

tess = DTess2D(100)
pp = Point2D(1.45, 1.49)
_pushunfixed!(tess, pp)
_flipa!(tess,2,3)

@test geta(tess._trigs[2]) == pp
@test getb(tess._trigs[2]) == pd
@test getc(tess._trigs[2]) == pc

@test geta(tess._trigs[3]) == pp
@test getb(tess._trigs[3]) == pb
@test getc(tess._trigs[3]) == pd

point_arr = Point2D[]
n=1000
tess = DTess2D(n*10)
for i in [1:n]
	push!(point_arr, Point2D(rand()+1.0, rand()+1.0))
end
push!(tess, point_arr)
@test tess._last_trig_index == n*2+3

for p in point_arr
	for t in tess._trigs[2:tess._last_trig_index]
		try
			@test incircle(t, p) <= 0 
		catch e
			if (p == geta(t)) || (p == getb(t)) || (p == getc(t))
				continue
			end
		end
	end
end

point_arr = Point2D[]
n=10
tess = DTess2D(n*n*10)
for x in linspace(1.001,1.999,n)
	for y in linspace(1.001,1.999,n)
		push!(point_arr, Point2D(x,y))
	end
end
push!(tess, point_arr)
@test tess._last_trig_index == n*n*2+3

for p in point_arr
	for t in tess._trigs[2:tess._last_trig_index]
		try
			@test incircle(t, p) <= 0 
		catch e
			if (p == geta(t)) || (p == getb(t)) || (p == getc(t))
				continue
			end
		end
	end
end

println("vor: All tests passed!")

end

test()

end # module VoronoiDelaunay