VoronoiDelaunay.jl

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 ([email protected])
# 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