Skip to content

Instantly share code, notes, and snippets.

@skariel

skariel/GeometricalPredicates.jl

Last active August 29, 2015 14:01
Show Gist options
  • Save skariel/da85943803a6f57a52fd to your computer and use it in GitHub Desktop.
Save skariel/da85943803a6f57a52fd to your computer and use it in GitHub Desktop.
Download ZIP
Raw
GeometricalPredicates.jl
module GeometricalPredicates
# GeometricalPredicates_v0.2.9
#
# Fast, robust 2D and 3D geometrical predicates on generic point types.
# 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!
#
# Calculations initially performed on Float64 while bounding max
# absolute errors. If unable to determine result fall back to exact
# calculation using BigInts. This is a form of floating point filtering.
# Most calculations are cached for fast repeated testing of
# incircle/intriangle
#
# Usage notes:
# -------------
# - all point types must implement getx(p), gety(p) and (if 3D point) getz(p).
# These methods should return Float64. The predicated are generic in the sense
# that they accept any such point carrying additional user defined data.
# - all point coordinates must be in the range [1, 2), i.e. including 1.0 and
# ending with inclusion of at 2.0-e where thiss represents the Float64 number
# which is closest to 2.0
# - beware of creating points (2D or 3D) using somthing like Point2D(3, 4) -
# it will return the point belonging to a Peano-Hilbert key `3` in
# an 2^4 X 2^4 grid. If you want to create a point located at 3, 4 use
# Point2D(3.0, 4.0)
#
# `incircle` method: determines if the point
# lies inside (ret->1), exactly on (ret->0),
# outside (ret->-1), or NA (ret->2) of a circle or sphere
# defined by the primitive
#
# `intriangle` method: determines if the point
# lies inside (ret->1) or outside (ret->negative num) the primitive.
# The negative number is the negative index ot the triangle point
# which test point is in front of. If test point is in front of two
# triangle points then one is chosen in an undefined mannar.
# If point is exactly on one of the sides it returns the
# index+1 of the point that is in front of said side
# i.e. for a point infront of Triangle.a return 2,
# for a point infront of Triangle.b return 3, etc.
# If test point is exactly on one of the triangle points then one
# side is chosen in an undefined mannar.
#
# `peanokey` method: this is the scale-dependent Peano-Hilbert interface.
# returns the Peano-Hilbert key for the given
# point in an nXn grid starting with 0 at (1.0, 2.0-e) and
# ending at nXn-1 at (2.0-e,2.0-e) for 2D, and starting at
# (1.0, 1.0, 1.0) with zero and ending with nXnXn-1 at
# (1.0, 1.0, 2.0-e) for 3D. `n` here is 2^bits, `bits` being a parameter
# passed to the function. For inverse Peano-Hilbert keys just create
# points using integers such as Point3D(1, 5) this will create a 3D point
# with peano index of 3 on a grid of size 2^5 X 2^5 X 2^5
#
# Todo
# ------------------
# - Make a package out of this
# - More comments in code
# - Documentation:
# - document scale-free spatial ordering
# - document primitive creation
# - add general examples
#
#
export
min_coord, max_coord,
AbstractPoint,
AbstractPoint2D,
AbstractPoint3D,
AbstractPrimitive,
AbstractOrientedPrimitive,
AbstractPositivelyOrientedPrimitive,
AbstractPositivelyOrientedTriangle,
AbstractPositivelyOrientedTetrahedron,
AbstractNegativelyOrientedPrimitive,
AbstractNegativelyOrientedTriangle,
AbstractNegativelyOrientedTetrahedron,
AbstractUnOrientedPrimitive,
AbstractTriangleUnOriented,
AbstractTetrahedronUnOriented,
TriangleTypes, TetrahedronTypes,
positivelyoriented, negativelyoriented, unoriented, orientation,
Point2D, Point3D, getx, gety, getz, geta, getb, getc, getd,
seta, setb, setc, setd, setabc, setabcd,
setab, setbc, setcd, setac, setad, setbd,
setabd, setacd, setbcd,
Triangle, Tetrahedron,
area, volume, centroid, circumcenter, circumradius2, incircle, intriangle,
peanokey, hilbertsort!, mssort!, clean!
const _float_err = eps(Float64)
const _abs_err_incircle_2d = 12*_float_err
const _abs_err_incircle_3d = 48*_float_err
const _abs_err_orientation_2d = 2*_float_err
const _abs_err_orientation_3d = 6*_float_err
const _abs_err_intriangle = 6*_float_err
const _abs_err_intriangle_zero = 2*_float_err
const _abs_err_intetra = 24*_float_err
const _abs_err_intetra_zero = 6*_float_err
const min_coord = 1.0
const max_coord = 2.0 - eps(Float64)
abstract AbstractPoint
abstract AbstractPoint2D <: AbstractPoint
abstract AbstractPoint3D <: AbstractPoint
abstract AbstractPrimitive
abstract AbstractUnOrientedPrimitive <: AbstractPrimitive
abstract AbstractOrientedPrimitive <: AbstractPrimitive
abstract AbstractPositivelyOrientedPrimitive <: AbstractOrientedPrimitive
abstract AbstractNegativelyOrientedPrimitive <: AbstractOrientedPrimitive
abstract AbstractTriangleUnOriented <: AbstractUnOrientedPrimitive
abstract AbstractTetrahedronUnOriented <: AbstractUnOrientedPrimitive
abstract AbstractPositivelyOrientedTriangle <: AbstractPositivelyOrientedPrimitive
abstract AbstractNegativelyOrientedTriangle <: AbstractNegativelyOrientedPrimitive
abstract AbstractPositivelyOrientedTetrahedron <: AbstractPositivelyOrientedPrimitive
abstract AbstractNegativelyOrientedTetrahedron <: AbstractNegativelyOrientedPrimitive
typealias TriangleTypes Union(AbstractTriangleUnOriented, AbstractPositivelyOrientedTriangle, AbstractNegativelyOrientedTriangle)
typealias TetrahedronTypes Union(AbstractTetrahedronUnOriented, AbstractPositivelyOrientedTetrahedron, AbstractNegativelyOrientedTetrahedron)
# standard 2D point
immutable Point2D <: AbstractPoint2D
_x::Float64
_y::Float64
end
getx(p::Point2D) = p._x
gety(p::Point2D) = p._y
# standard 3D point
immutable Point3D <: AbstractPoint3D
_x::Float64
_y::Float64
_z::Float64
end
getx(p::Point3D) = p._x
gety(p::Point3D) = p._y
getz(p::Point3D) = p._z
macro _define_triangle_type(name, abstracttype)
oriented = !contains(string(name), "UnOriented")
esc(parse("""
type $name{T<:AbstractPoint2D} <: $abstracttype
_a::T; _b::T; _c::T
_bx::Float64; _by::Float64
_cx::Float64; _cy::Float64
_px::Float64; _py::Float64
_pr2::Float64
_sz::Float32
$(oriented ? "" : "_o::Int8")
function $name(a::T, b::T, c::T)
t = new(a, b, c, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0$(oriented? "":", 0"))
clean!(t)
t
end
end
"""))
end
@_define_triangle_type(UnOrientedTriangle, AbstractTriangleUnOriented)
@_define_triangle_type(PositivelyOrientedTriangle, AbstractPositivelyOrientedTriangle)
@_define_triangle_type(NegativelyOrientedTriangle, AbstractNegativelyOrientedTriangle)
abstract AbstractOrientation
type PositivelyOriented <: AbstractOrientation; end
type NegativelyOriented <: AbstractOrientation; end
type UnOriented <: AbstractOrientation; end
const positivelyoriented = PositivelyOriented()
const negativelyoriented = NegativelyOriented()
const unoriented = UnOriented()
Triangle{T<:AbstractPoint2D}(a::T, b::T, c::T, ::PositivelyOriented) = PositivelyOrientedTriangle{T}(a, b, c)
Triangle{T<:AbstractPoint2D}(a::T, b::T, c::T, ::NegativelyOriented) = NegativelyOrientedTriangle{T}(a, b, c)
Triangle{T<:AbstractPoint2D}(a::T, b::T, c::T, ::UnOriented) = UnOrientedTriangle{T}(a, b, c)
Triangle{T<:AbstractPoint2D}(a::T, b::T, c::T) = Triangle(a, b, c, unoriented)
Triangle(ax::Float64, ay::Float64, bx::Float64, by::Float64, cx::Float64, cy::Float64, orientation::AbstractOrientation=unoriented) =
Triangle(Point2D(ax, ay), Point2D(bx, by), Point2D(cx, cy), orientation)
area(tr::TriangleTypes) = abs(tr._pr2)/2
centroid(tr::TriangleTypes) =
Point2D(
(getx(geta(tr)) + getx(getb(tr)) + getx(getc(tr))) / 3.0,
(gety(geta(tr)) + gety(getb(tr)) + gety(getc(tr))) / 3.0)
function circumcenter(tr::TriangleTypes)
const d = -2.0 * tr._pr2
Point2D(
tr._px/d + getx(geta(tr)),
tr._py/d + gety(geta(tr))
)
end
function circumradius2(tr::TriangleTypes)
const c = circumcenter(tr)
const x = getx(c) - getx(geta(tr))
const y = gety(c) - gety(geta(tr))
x*x + y*y
end
macro _define_tetrahedron_type(name, abstracttype)
oriented = !contains(string(name), "UnOriented")
esc(parse("""
type $name{T<:AbstractPoint3D} <: $abstracttype
_a::T; _b::T; _c::T; _d::T
_bx::Float64; _by::Float64; _bz::Float64
_cx::Float64; _cy::Float64; _cz::Float64
_dx::Float64; _dy::Float64; _dz::Float64
_px::Float64; _py::Float64; _pz::Float64
_pr2::Float64
_sz::Float32
$(oriented? "":"_o::Int8")
function $name(a::AbstractPoint3D, b::AbstractPoint3D, c::AbstractPoint3D, d::AbstractPoint3D)
t = new(a, b, c, d, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 $(oriented? "" : ", 0"))
clean!(t)
t
end
end
"""))
end
@_define_tetrahedron_type(UnOrientedTetrahedron, AbstractTetrahedronUnOriented)
@_define_tetrahedron_type(PositivelyOrientedTetrahedron, AbstractPositivelyOrientedTetrahedron)
@_define_tetrahedron_type(NegativelyOrientedTetrahedron, AbstractNegativelyOrientedTetrahedron)
Tetrahedron{T<:AbstractPoint3D}(a::T, b::T, c::T, d::T, ::PositivelyOriented) = PositivelyOrientedTetrahedron{T}(a, b, c, d)
Tetrahedron{T<:AbstractPoint3D}(a::T, b::T, c::T, d::T, ::NegativelyOriented) = NegativelyOrientedTetrahedron{T}(a, b, c, d)
Tetrahedron{T<:AbstractPoint3D}(a::T, b::T, c::T, d::T, ::UnOriented) = UnOrientedTetrahedron{T}(a, b, c, d)
Tetrahedron{T<:AbstractPoint3D}(a::T, b::T, c::T, d::T) = Tetrahedron(a, b, c, d, unoriented)
Tetrahedron(ax::Float64, ay::Float64, az::Float64, bx::Float64, by::Float64, bz::Float64,
cx::Float64, cy::Float64, cz::Float64, dx::Float64, dy::Float64, dz::Float64, orientation::AbstractOrientation=unoriented) =
Tetrahedron(Point3D(ax,ay,az), Point3D(bx,by,bz), Point3D(cx,cy,cz), Point3D(dx,dy,dz), orientation)
volume(tr::TetrahedronTypes) = abs(tr._pr2)/2
centroid(tr::TetrahedronTypes) =
Point3D(
(getx(geta(tr)) + getx(getb(tr)) + getx(getc(tr)) + getx(getd(tr))) / 4.0,
(gety(geta(tr)) + gety(getb(tr)) + gety(getc(tr)) + gety(getd(tr))) / 4.0,
(getz(geta(tr)) + getz(getb(tr)) + getz(getc(tr)) + getz(getd(tr))) / 4.0)
function circumcenter(tr::TetrahedronTypes)
const d = -2.0 * tr._pr2
Point3D(
tr._px/d + getx(geta(tr)),
tr._py/d + gety(geta(tr)),
tr._pz/d + getz(geta(tr))
)
end
function circumradius2(tr::TetrahedronTypes)
const c = circumcenter(tr)
const x = getx(c) - getx(geta(tr))
const y = gety(c) - gety(geta(tr))
const z = getz(c) - getz(geta(tr))
x*x + y*y + z*z
end
# extract exact integer representation of float to be used in exact calculations when needed
_extract_int(n::Float64) = reinterpret(Uint64, n) & 0x000fffffffffffff
_extract_bigint(n::Float64) = BigInt(_extract_int(n))
# functions to re-validate cached pre-calculations in primitives
function _clean!(t::TriangleTypes)
t._bx = getx(getb(t))-getx(geta(t)); t._by = gety(getb(t))-gety(geta(t));
t._cx = getx(getc(t))-getx(geta(t)); t._cy = gety(getc(t))-gety(geta(t));
const br2 = t._bx*t._bx+t._by*t._by
const cr2 = t._cx*t._cx+t._cy*t._cy
t._sz = abs(t._bx)+abs(t._by)+abs(t._cx)+abs(t._cy)
t._px = br2*t._cy - t._by*cr2
t._py = -br2*t._cx + t._bx*cr2
t._pr2 = -t._bx *t._cy + t._by*t._cx
end
function clean!(t::AbstractTriangleUnOriented)
_clean!(t)
if t._pr2 < -_abs_err_orientation_2d*t._sz
t._o = 1
elseif t._pr2 > _abs_err_orientation_2d*t._sz
t._o = -1
else
t._o = _exact_sign_orientation_determinant!(
_extract_bigint(getx(geta(t))), _extract_bigint(gety(geta(t))),
_extract_bigint(getx(getb(t))), _extract_bigint(gety(getb(t))),
_extract_bigint(getx(getc(t))), _extract_bigint(gety(getc(t))))
end
end
function _clean!(t::TetrahedronTypes)
t._bx = getx(getb(t))-getx(geta(t)); t._by = gety(getb(t))-gety(geta(t)); t._bz = getz(getb(t))-getz(geta(t))
t._cx = getx(getc(t))-getx(geta(t)); t._cy = gety(getc(t))-gety(geta(t)); t._cz = getz(getc(t))-getz(geta(t))
t._dx = getx(getd(t))-getx(geta(t)); t._dy = gety(getd(t))-gety(geta(t)); t._dz = getz(getd(t))-getz(geta(t))
const br2 = t._bx*t._bx+t._by*t._by+t._bz*t._bz
const cr2 = t._cx*t._cx+t._cy*t._cy+t._cz*t._cz
const dr2 = t._dx*t._dx+t._dy*t._dy+t._dz*t._dz
t._sz = abs(t._bx)+abs(t._by)+abs(t._cx)+abs(t._cy)+abs(t._dx)+abs(t._dy)
t._px = t._by*t._cz*dr2 - br2*t._cz*t._dy - t._by*cr2*t._dz - t._bz*t._cy*dr2 + br2*t._cy*t._dz + t._bz*cr2*t._dy
t._py = br2*t._cz*t._dx + t._bz*t._cx*dr2 - br2*t._cx*t._dz - t._bx*t._cz*dr2 - t._bz*cr2*t._dx + t._bx*cr2*t._dz
t._pz = br2*t._cx*t._dy + t._bx*t._cy*dr2 + t._by*cr2*t._dx - br2*t._cy*t._dx - t._bx*cr2*t._dy - t._by*t._cx*dr2
t._pr2 = -t._bx*t._cy*t._dz + t._bx*t._cz*t._dy + t._by*t._cx*t._dz - t._by*t._cz*t._dx - t._bz*t._cx*t._dy + t._bz*t._cy*t._dx
end
function clean!(t::AbstractTetrahedronUnOriented)
_clean!(t)
# calculate the orientation
if t._pr2 < -_abs_err_orientation_3d
t._o = 1
elseif t._pr2 > _abs_err_orientation_3d
t._o = -1
else
# exact calculation is required
t._o = _exact_sign_orientation_determinant!(
_extract_bigint(getx(geta(t))), _extract_bigint(gety(geta(t))), _extract_bigint(getz(geta(t))),
_extract_bigint(getx(getb(t))), _extract_bigint(gety(getb(t))), _extract_bigint(getz(getb(t))),
_extract_bigint(getx(getc(t))), _extract_bigint(gety(getc(t))), _extract_bigint(getz(getc(t))),
_extract_bigint(getx(getd(t))), _extract_bigint(gety(getd(t))), _extract_bigint(getz(getd(t))))
end
end
clean!(t::AbstractOrientedPrimitive) = _clean!(t)
# getting points of the primitive
geta(t::AbstractPrimitive) = t._a
getb(t::AbstractPrimitive) = t._b
getc(t::AbstractPrimitive) = t._c
getd(t::TetrahedronTypes) = t._d
# changing points in the primitive
seta(t::AbstractPrimitive, p::AbstractPoint) = (t._a=p; clean!(t))
setb(t::AbstractPrimitive, p::AbstractPoint) = (t._b=p; clean!(t))
setc(t::AbstractPrimitive, p::AbstractPoint) = (t._c=p; clean!(t))
setd(t::TetrahedronTypes, p::AbstractPoint3D) = (t._d=p; clean!(t))
setabc(t::AbstractPrimitive, pa::AbstractPoint, pb::AbstractPoint, pc::AbstractPoint) = (t._a=pa; t._b=pb; t._c=pc; clean!(t))
setab(t::AbstractPrimitive, pa::AbstractPoint, pb::AbstractPoint) = (t._a=pa; t._b=pb; clean!(t))
setbc(t::AbstractPrimitive, pb::AbstractPoint, pc::AbstractPoint) = (t._b=pb; t._c=pc; clean!(t))
setac(t::AbstractPrimitive, pa::AbstractPoint, pc::AbstractPoint) = (t._a=pa; t._c=pc; clean!(t))
setabcd(t::TetrahedronTypes, pa::AbstractPoint3D, pb::AbstractPoint3D, pc::AbstractPoint3D, pd::AbstractPoint3D) = (t._a=pa; t._b=pb; t._c=pc; t._d=pd; clean!(t))
setabc(t::TetrahedronTypes, pa::AbstractPoint3D, pb::AbstractPoint3D, pc::AbstractPoint3D) = (t._a=pa; t._b=pb; t._c=pc; clean!(t))
setabd(t::TetrahedronTypes, pa::AbstractPoint3D, pb::AbstractPoint3D, pd::AbstractPoint3D) = (t._a=pa; t._b=pb; t._d=pd; clean!(t))
setacd(t::TetrahedronTypes, pa::AbstractPoint3D, pc::AbstractPoint3D, pd::AbstractPoint3D) = (t._a=pa; t._c=pc; t._d=pd; clean!(t))
setbcd(t::TetrahedronTypes, pb::AbstractPoint3D, pc::AbstractPoint3D, pd::AbstractPoint3D) = (t._b=pb; t._c=pc; t._d=pd; clean!(t))
setad(t::TetrahedronTypes, pa::AbstractPoint3D, pd::AbstractPoint3D) = (t._a=pa; t._d=pd; clean!(t))
setbd(t::TetrahedronTypes, pb::AbstractPoint3D, pd::AbstractPoint3D) = (t._b=pb; t._d=pd; clean!(t))
setcd(t::TetrahedronTypes, pc::AbstractPoint3D, pd::AbstractPoint3D) = (t._c=pc; t._d=pd; clean!(t))
orientation(p::AbstractUnOrientedPrimitive) = p._o
orientation(p::AbstractPositivelyOrientedPrimitive) = 1
orientation(p::AbstractNegativelyOrientedPrimitive) = -1
orientation(ax::Float64, ay::Float64, bx::Float64, by::Float64, cx::Float64, cy::Float64) =
orientation(Triangle(ax, ay, bx, by, cx, cy))
orientation(ax::Float64, ay::Float64, az::Float64, bx::Float64, by::Float64, bz::Float64, cx::Float64, cy::Float64, cz::Float64, dx::Float64, dy::Float64, dz::Float64) =
orientation(Tetrahedron(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz))
# exact orientation for triangle
function _exact_sign_orientation_determinant!(ax::BigInt, ay::BigInt, bx::BigInt, by::BigInt, cx::BigInt, cy::BigInt)
bx -= ax; by -= ay
cx -= ax; cy -= ay
sign(bx*cy - by*cx)
end
# exact orientation for tetrahedron
function _exact_sign_orientation_determinant!(ax::BigInt, ay::BigInt, az::BigInt, bx::BigInt, by::BigInt, bz::BigInt, cx::BigInt, cy::BigInt, cz::BigInt, dx::BigInt, dy::BigInt, dz::BigInt)
bx -= ax; by -= ay; bz -= az
cx -= ax; cy -= ay; cz -= az
dx -= ax; dy -= ay; dz -= az
sign(+bx*cy*dz - bx*cz*dy - by*cx*dz + by*cz*dx + bz*cx*dy - bz*cy*dx)
end
# exact incircle for triangle
function _exact_sign_incircle_determinant!(ax::BigInt, ay::BigInt, bx::BigInt, by::BigInt, cx::BigInt, cy::BigInt, px::BigInt, py::BigInt)
bx -= ax; by -= ay;
cx -= ax; cy -= ay;
px -= ax; py -= ay;
const br2 = bx*bx+by*by
const cr2 = cx*cx+cy*cy
const pr2 = px*px+py*py
sign(-br2*cx*py + br2*cy*px + bx*cr2*py - bx*cy*pr2 - by*cr2*px + by*cx*pr2)
end
# exact incircle for tetrahedron
function _exact_sign_incircle_determinant!(ax::BigInt, ay::BigInt, az::BigInt, bx::BigInt, by::BigInt, bz::BigInt, cx::BigInt, cy::BigInt, cz::BigInt, dx::BigInt, dy::BigInt, dz::BigInt, px::BigInt, py::BigInt, pz::BigInt)
bx -= ax; by -= ay; bz -= az
cx -= ax; cy -= ay; cz -= az
dx -= ax; dy -= ay; dz -= az
px -= ax; py -= ay; pz -= az
const br2 = bx*bx+by*by+bz*bz
const cr2 = cx*cx+cy*cy+cz*cz
const dr2 = dx*dx+dy*dy+dz*dz
const pr2 = px*px+py*py+pz*pz
sign(
+br2*cx*dy*pz - br2*cx*dz*py - br2*cy*dx*pz + br2*cy*dz*px +
br2*cz*dx*py - br2*cz*dy*px - bx*cr2*dy*pz + bx*cr2*dz*py +
bx*cy*dr2*pz - bx*cy*dz*pr2 - bx*cz*dr2*py + bx*cz*dy*pr2 +
by*cr2*dx*pz - by*cr2*dz*px - by*cx*dr2*pz + by*cx*dz*pr2 +
by*cz*dr2*px - by*cz*dx*pr2 - bz*cr2*dx*py + bz*cr2*dy*px +
bz*cx*dr2*py - bz*cx*dy*pr2 - bz*cy*dr2*px + bz*cy*dx*pr2)
end
# filtered (fast and exact) incircle for triangle. This one is exported
function incircle(t::TriangleTypes, p::AbstractPoint2D)
orientation(t) == 0 && return 2
px = getx(p) - getx(geta(t))
py = gety(p) - gety(geta(t))
sz = abs(px)+abs(py)+t._sz
pr2 = px*px + py*py
d = t._px*px + t._py*py + t._pr2*pr2
if d < -_abs_err_incircle_2d*sz
return -orientation(t)
elseif d > _abs_err_incircle_2d*sz
return orientation(t)
end
orientation(t)*_exact_sign_incircle_determinant!(
_extract_bigint(getx(geta(t))), _extract_bigint(gety(geta(t))),
_extract_bigint(getx(getb(t))), _extract_bigint(gety(getb(t))),
_extract_bigint(getx(getc(t))), _extract_bigint(gety(getc(t))),
_extract_bigint(getx(p)) , _extract_bigint(gety(p)))
end
# filtered (fast and exact) incircle for tetrahedron. This one is exported
function incircle(t::TetrahedronTypes, p::AbstractPoint3D)
orientation(t) == 0 && return 2
px = getx(p) - getx(geta(t))
py = gety(p) - gety(geta(t))
pz = getz(p) - getz(geta(t))
sz = abs(px)+abs(py)+abs(pz)+t._sz
pr2 = px*px + py*py + pz*pz
d = t._px*px + t._py*py + t._pz*pz + t._pr2*pr2
if d < -_abs_err_incircle_3d*sz
return -orientation(t)
elseif d > _abs_err_incircle_3d*sz
return orientation(t)
end
orientation(t)*_exact_sign_incircle_determinant!(
_extract_bigint(getx(geta(t))), _extract_bigint(gety(geta(t))), _extract_bigint(getz(geta(t))),
_extract_bigint(getx(getb(t))), _extract_bigint(gety(getb(t))), _extract_bigint(getz(getb(t))),
_extract_bigint(getx(getc(t))), _extract_bigint(gety(getc(t))), _extract_bigint(getz(getc(t))),
_extract_bigint(getx(getd(t))), _extract_bigint(gety(getd(t))), _extract_bigint(getz(getd(t))),
_extract_bigint(getx(p)) , _extract_bigint(gety(p)) , _extract_bigint(getz(p)))
end
# helper methods to use incircle directly with coordinates
incircle(ax::Float64, ay::Float64, bx::Float64, by::Float64, cx::Float64, cy::Float64, dx::Float64, dy::Float64) =
incircle(Triangle(ax, ay, bx, by, cx, cy), Point2D(dx, dy))
incircle(ax::Float64, ay::Float64, az::Float64, bx::Float64, by::Float64, bz::Float64, cx::Float64, cy::Float64, cz::Float64, dx::Float64, dy::Float64, dz::Float64, ex::Float64, ey::Float64, ez::Float64) =
incircle(Tetrahedron(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz), Point3D(ex, ey, ez))
# exact intriangle (slow!)
function _exact_intriangle!(ax::BigInt, ay::BigInt, bx::BigInt, by::BigInt, cx::BigInt, cy::BigInt, px::BigInt, py::BigInt)
bx -= ax; by -= ay;
cx -= ax; cy -= ay;
px -= ax; py -= ay;
const nb = -cx*py + cy*px
const nc = bx*py - by*px
const denom = bx*cy - by*cx
if sign(nb) * sign(denom) < 0
return -2
end
if sign(nc) * sign(denom) < 0
return -3
end
const l = nb+nc - denom
const sl = sign(l)*sign(denom)
if sl > 0
return -1
end
if sign(nb) == 0
return 3
end
if sign(nc) == 0
return 4
end
if sl == 0
return 2
end
-sl
end
function _exact_intriangle!(ax::BigInt, ay::BigInt, az::BigInt, bx::BigInt, by::BigInt, bz::BigInt, cx::BigInt, cy::BigInt, cz::BigInt, dx::BigInt, dy::BigInt, dz::BigInt, px::BigInt, py::BigInt, pz::BigInt)
bx -= ax; by -= ay; bz -= az
cx -= ax; cy -= ay; cz -= az
dx -= ax; dy -= ay; dz -= az
px -= ax; py -= ay; pz -= az
const denom = bx*cy*dz-bx*cz*dy-by*cx*dz+by*cz*dx+bz*cx*dy-bz*cy*dx
const nb = cx*dy*pz-cx*dz*py-cy*dx*pz+cy*dz*px+cz*dx*py-cz*dy*px
if sign(nb) * sign(denom) < 0
return -2
end
const nc = -bx*dy*pz+bx*dz*py+by*dx*pz-by*dz*px-bz*dx*py+bz*dy*px
if sign(nc) * sign(denom) < 0
return -3
end
const nd = bx*cy*pz-bx*cz*py-by*cx*pz+by*cz*px+bz*cx*py-bz*cy*px
if sign(nd) * sign(denom) < 0
return -4
end
const l = (nb+nc+nd - denom) * sign(denom)
const sl = sign(l)
if sl > 0
return -1
end
if sign(nb) == 0
return 3
end
if sign(nc) == 0
return 4
end
if sign(nd) == 0
return 5
end
if sl == 0
return 2
end
-sl
end
# helper methods to use the exact intriangle directly with coordinates
_exact_intriangle(t::TriangleTypes, p::AbstractPoint2D) =
_exact_intriangle!(
_extract_bigint(getx(geta(t))), _extract_bigint(gety(geta(t))),
_extract_bigint(getx(getb(t))), _extract_bigint(gety(getb(t))),
_extract_bigint(getx(getc(t))), _extract_bigint(gety(getc(t))),
_extract_bigint(getx(p)) , _extract_bigint(gety(p)))
_exact_intriangle(t::TetrahedronTypes, p::AbstractPoint3D) =
_exact_intriangle!(
_extract_bigint(getx(geta(t))), _extract_bigint(gety(geta(t))), _extract_bigint(getz(geta(t))),
_extract_bigint(getx(getb(t))), _extract_bigint(gety(getb(t))), _extract_bigint(getz(getb(t))),
_extract_bigint(getx(getc(t))), _extract_bigint(gety(getc(t))), _extract_bigint(getz(getc(t))),
_extract_bigint(getx(getd(t))), _extract_bigint(gety(getd(t))), _extract_bigint(getz(getd(t))),
_extract_bigint(getx(p)), _extract_bigint(gety(p)), _extract_bigint(getz(p)))
# filtered intriangle for triangles (fast, exact). This one is exported
function intriangle(t::TriangleTypes, p::AbstractPoint2D)
const px = getx(p) - getx(geta(t)); const py = gety(p) - gety(geta(t))
sz = abs(px)+abs(py)+t._sz
const nb = (-t._cx*py + t._cy*px) * sign(-t._pr2)
if nb < -_abs_err_intriangle_zero*sz
return -2
elseif nb < _abs_err_intriangle_zero*sz
# we need an exact calculation
return _exact_intriangle(t, p)
end
const nc = (t._bx*py - t._by*px) * sign(-t._pr2)
if nc < -_abs_err_intriangle_zero*sz
return -3
elseif nc < _abs_err_intriangle_zero*sz
# we need an exact calculation
return _exact_intriangle(t, p)
end
const l = nb+nc + t._pr2*sign(-t._pr2)
if l < -_abs_err_intriangle*sz
return 1
elseif l > _abs_err_intriangle*sz
return -1
end
# we need an exact calculation
_exact_intriangle(t, p)
end
function intriangle(t::TetrahedronTypes, p::AbstractPoint3D)
const px = getx(p) - getx(geta(t)); const py = gety(p) - gety(geta(t)); const pz = getz(p) - getz(geta(t))
sz = abs(px)+abs(py)+abs(pz)+t._sz
const nb = (t._cx*t._dy*pz-t._cx*t._dz*py-t._cy*t._dx*pz+t._cy*t._dz*px+t._cz*t._dx*py-t._cz*t._dy*px) * sign(-t._pr2)
if nb < -_abs_err_intetra_zero*sz
return -2
elseif nb < _abs_err_intetra_zero*sz
# we need an exact calculation
return _exact_intriangle(t, p)
end
const nc = (-t._bx*t._dy*pz+t._bx*t._dz*py+t._by*t._dx*pz-t._by*t._dz*px-t._bz*t._dx*py+t._bz*t._dy*px) * sign(-t._pr2)
if nc < -_abs_err_intetra_zero*sz
return -3
elseif nc < _abs_err_intetra_zero*sz
# we need an exact calculation
return _exact_intriangle(t, p)
end
const nd = (t._bx*t._cy*pz-t._bx*t._cz*py-t._by*t._cx*pz+t._by*t._cz*px+t._bz*t._cx*py-t._bz*t._cy*px) * sign(-t._pr2)
if nd < -_abs_err_intetra_zero*sz
return -4
elseif nd < _abs_err_intetra_zero*sz
# we need an exact calculation
return _exact_intriangle(t, p)
end
const l = nb+nc+nd + t._pr2*sign(-t._pr2)
if l < -_abs_err_intetra*sz
return 1
elseif l > _abs_err_intetra*sz
return -1
end
# we need an exact calculation
_exact_intriangle(t, p)
end
# helper methods to use the filtered (fast, exact) intriangle directly with raw coordinates
intriangle(ax::Float64, ay::Float64, bx::Float64, by::Float64, cx::Float64, cy::Float64, px::Float64, py::Float64) =
intriangle(Triangle(ax, ay, bx, by, cx, cy), Point2D(px, py))
intriangle(ax::Float64, ay::Float64, az::Float64, bx::Float64, by::Float64, bz::Float64, cx::Float64, cy::Float64, cz::Float64, dx::Float64, dy::Float64, dz::Float64, px::Float64, py::Float64, pz::Float64) =
intriangle(Tetrahedron(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz), Point3D(px, py, pz))
###################################################################################
#
# Hilbert stuff
#
#
# number of bits to use per dimension in calculating the peano-key
const peano_2D_bits = 31
const peano_3D_bits = 21
# implementing 2D scale dependednt Peano-Hilbert indexing
_extract_peano_bin_num(nbins::Int64, n::Float64) = itrunc( (n-1)*nbins )
# calculate peano key for given point
function peanokey(p::AbstractPoint2D, bits::Int64=peano_2D_bits)
const n = 1 << bits
s = n >> 1; d = 0
x = _extract_peano_bin_num(n, getx(p))
y = _extract_peano_bin_num(n, gety(p))
while true
rx = (x & s) > 0
ry = (y & s) > 0
d += s * s * ((3 * rx) $ ry)
s = s >> 1
(s == 0) && break
if ry == 0
if rx == 1
x = n - 1 - x;
y = n - 1 - y;
end
x, y = y, x
end
end
d
end
# Inverse calculation. I.e. calculate the point that given given peano key
function Point2D(peanokey::Int64, bits::Int64=peano_2D_bits)
const n = 1 << bits
x = 0; y = 0; s=1
while true
rx = 1 & (peanokey >> 1)
ry = 1 & (peanokey $ rx)
if ry == 0
if rx == 1
x = s - 1 - x;
y = s - 1 - y;
end
x, y = y, x
end
x += s * rx
y += s * ry
s = s << 1
(s >= n) && break
peanokey = peanokey >> 2
end
Point2D(1+x/n, 1+y/n)
end
# implementing 3D scaleful Peano-Hilbert indexing
const quadrants_arr = [
0, 7, 1, 6, 3, 4, 2, 5,
7, 4, 6, 5, 0, 3, 1, 2,
4, 3, 5, 2, 7, 0, 6, 1,
3, 0, 2, 1, 4, 7, 5, 6,
1, 0, 6, 7, 2, 3, 5, 4,
0, 3, 7, 4, 1, 2, 6, 5,
3, 2, 4, 5, 0, 1, 7, 6,
2, 1, 5, 6, 3, 0, 4, 7,
6, 1, 7, 0, 5, 2, 4, 3,
1, 2, 0, 3, 6, 5, 7, 4,
2, 5, 3, 4, 1, 6, 0, 7,
5, 6, 4, 7, 2, 1, 3, 0,
7, 6, 0, 1, 4, 5, 3, 2,
6, 5, 1, 2, 7, 4, 0, 3,
5, 4, 2, 3, 6, 7, 1, 0,
4, 7, 3, 0, 5, 6, 2, 1,
6, 7, 5, 4, 1, 0, 2, 3,
7, 0, 4, 3, 6, 1, 5, 2,
0, 1, 3, 2, 7, 6, 4, 5,
1, 6, 2, 5, 0, 7, 3, 4,
2, 3, 1, 0, 5, 4, 6, 7,
3, 4, 0, 7, 2, 5, 1, 6,
4, 5, 7, 6, 3, 2, 0, 1,
5, 2, 6, 1, 4, 3, 7, 0]
quadrants(a::Int64, b::Int64, c::Int64, d::Int64) = (@inbounds const x = quadrants_arr[1+a<<3+b<<2+c<<1+d]; x)
rotxmap_table = [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 0, 1, 2, 3, 17, 18, 19, 16, 23, 20, 21, 22]
rotymap_table = [1, 2, 3, 0, 16, 17, 18, 19, 11, 8, 9, 10, 22, 23, 20, 21, 14, 15, 12, 13, 4, 5, 6, 7]
rotx_table = [3, 0, 0, 2, 2, 0, 0, 1]
roty_table = [0, 1, 1, 2, 2, 3, 3, 0]
sense_table = [-1, -1, -1, +1, +1, -1, -1, -1]
function peanokey(p::AbstractPoint3D, bits::Int64=peano_3D_bits)
const n = 1 << bits
const x = _extract_peano_bin_num(n, getx(p))
const y = _extract_peano_bin_num(n, gety(p))
const z = _extract_peano_bin_num(n, getz(p))
mask = 1 << (bits - 1)
key = 0
rotation = 0
sense = 1
for i in 1:bits
const bitx = (x & mask > 0) ? 1 : 0
const bity = (y & mask > 0) ? 1 : 0
const bitz = (z & mask > 0) ? 1 : 0
const quad = quadrants(rotation, bitx, bity, bitz)
key <<= 3
key += sense == 1 ? quad : 7-quad
@inbounds rotx = rotx_table[quad+1]
@inbounds roty = roty_table[quad+1]
@inbounds sense *= sense_table[quad+1]
while rotx > 0
@inbounds rotation = rotxmap_table[rotation+1]
rotx -= 1
end
while(roty > 0)
@inbounds rotation = rotymap_table[rotation+1]
roty -= 1
end
mask >>= 1
end
key
end
const quadrants_inverse_x = Array(Int64, (24,8))
const quadrants_inverse_y = Array(Int64, (24,8))
const quadrants_inverse_z = Array(Int64, (24,8))
function _init_inv_peano_3d()
for rotation in 0:23
for bitx in 0:1
for bity in 0:1
for bitz in 0:1
quad = quadrants(rotation, bitx, bity, bitz)
quadrants_inverse_x[rotation+1, quad+1] = bitx;
quadrants_inverse_y[rotation+1, quad+1] = bity;
quadrants_inverse_z[rotation+1, quad+1] = bitz;
end
end
end
end
end
_init_inv_peano_3d()
# inverse transformation, give a key get a point
function Point3D(key::Int64, bits::Int64=peano_3D_bits)
const n = 1 << bits
shift = 3*(bits - 1)
mask = 7 << shift
rotation = 0
sense = 1
x = 0
y = 0
z = 0
for i in 1:bits
keypart = (key & mask) >> shift
quad = sense == 1 ? keypart : 7 - keypart
x = (x << 1) + quadrants_inverse_x[rotation+1, quad+1]
y = (y << 1) + quadrants_inverse_y[rotation+1, quad+1]
z = (z << 1) + quadrants_inverse_z[rotation+1, quad+1]
rotx = rotx_table[quad+1]
roty = roty_table[quad+1]
sense *= sense_table[quad+1]
while rotx > 0
rotation = rotxmap_table[rotation+1]
rotx -= 1
end
while roty > 0
rotation = rotymap_table[rotation+1]
roty -= 1
end
mask >>= 3
shift -= 3
end
Point3D(1+x/n, 1+y/n, 1+z/n)
end
# implementing scale-free Hilbert ordering. Real all about it here:
# http://doc.cgal.org/latest/Spatial_sorting/index.html
abstract AbstractCoordinate
type CoordinateX <: AbstractCoordinate end
type CoordinateY <: AbstractCoordinate end
type CoordinateZ <: AbstractCoordinate end
const coordinatex = CoordinateX()
const coordinatey = CoordinateY()
const coordinatez = CoordinateZ()
next2d(::CoordinateX) = coordinatey
next2d(::CoordinateY) = coordinatex
next3d(::CoordinateX) = coordinatey
next3d(::CoordinateY) = coordinatez
next3d(::CoordinateZ) = coordinatex
nextnext3d(::CoordinateX) = coordinatez
nextnext3d(::CoordinateY) = coordinatex
nextnext3d(::CoordinateZ) = coordinatey
abstract AbstractDirection
type Forward <: AbstractDirection end
type Backward <: AbstractDirection end
const forward = Forward()
const backward = Backward()
!(::Forward) = backward
!(::Backward) = forward
compare(::Forward, ::CoordinateX, p1::AbstractPoint, p2::AbstractPoint) = getx(p1) < getx(p2)
compare(::Backward, ::CoordinateX, p1::AbstractPoint, p2::AbstractPoint) = getx(p1) > getx(p2)
compare(::Forward, ::CoordinateY, p1::AbstractPoint, p2::AbstractPoint) = gety(p1) < gety(p2)
compare(::Backward, ::CoordinateY, p1::AbstractPoint, p2::AbstractPoint) = gety(p1) > gety(p2)
compare(::Forward, ::CoordinateZ, p1::AbstractPoint, p2::AbstractPoint) = getz(p1) < getz(p2)
compare(::Backward, ::CoordinateZ, p1::AbstractPoint, p2::AbstractPoint) = getz(p1) > getz(p2)
function select!{T<:AbstractPoint}(direction::AbstractDirection, coordinate::AbstractCoordinate, v::Array{T,1}, k::Int, lo::Int, hi::Int)
lo <= k <= hi || error("select index $k is out of range $lo:$hi")
@inbounds while lo < hi
if hi-lo == 1
if compare(direction, coordinate, v[hi], v[lo])
v[lo], v[hi] = v[hi], v[lo]
end
return v[k]
end
pivot = v[(lo+hi)>>>1]
i, j = lo, hi
while true
while compare(direction, coordinate, v[i], pivot); i += 1; end
while compare(direction, coordinate, pivot, v[j]); j -= 1; end
i <= j || break
v[i], v[j] = v[j], v[i]
i += 1; j -= 1
end
if k <= j
hi = j
elseif i <= k
lo = i
else
return pivot
end
end
return v[lo]
end
function hilbertsort!{T<:AbstractPoint2D}(directionx::AbstractDirection, directiony::AbstractDirection, coordinate::AbstractCoordinate, a::Array{T,1}, lo::Int64, hi::Int64, lim::Int64=4)
hi-lo <= lim && return a
const i2 = (lo+hi)>>>1
const i1 = (lo+i2)>>>1
const i3 = (i2+hi)>>>1
select!(directionx, coordinate, a, i2, lo, hi)
select!(directiony, next2d(coordinate), a, i1, lo, i2)
select!(!directiony, next2d(coordinate), a, i3, i2, hi)
hilbertsort!(directiony, directionx, next2d(coordinate), a, lo, i1, lim)
hilbertsort!(directionx, directiony, coordinate, a, i1, i2, lim)
hilbertsort!(directionx, directiony, coordinate, a, i2, i3, lim)
hilbertsort!(!directiony, !directionx, next2d(coordinate), a, i3, hi, lim)
return a
end
function hilbertsort!{T<:AbstractPoint3D}(directionx::AbstractDirection, directiony::AbstractDirection, directionz::AbstractDirection, coordinate::AbstractCoordinate, a::Array{T,1}, lo::Int64, hi::Int64, lim::Int64=8)
hi-lo <= lim && return a
const i4 = (lo+hi)>>>1
const i2 = (lo+i4)>>>1
const i1 = (lo+i2)>>>1
const i3 = (i2+i4)>>>1
const i6 = (i4+hi)>>>1
const i5 = (i4+i6)>>>1
const i7 = (i6+hi)>>>1
select!(directionx, coordinate, a, i4, lo, hi)
select!(directiony, next3d(coordinate), a, i2, lo, i4)
select!(directionz, nextnext3d(coordinate), a, i1, lo, i2)
select!(!directionz, nextnext3d(coordinate), a, i3, i2, i4)
select!(!directiony, next3d(coordinate), a, i6, i4, hi)
select!(directionz, nextnext3d(coordinate), a, i5, i4, i6)
select!(!directionz, nextnext3d(coordinate), a, i7, i6, hi)
hilbertsort!( directionz, directionx, directiony, nextnext3d(coordinate), a, lo, i1, lim)
hilbertsort!( directiony, directionz, directionx, next3d(coordinate), a, i1, i2, lim)
hilbertsort!( directiony, directionz, directionx, next3d(coordinate), a, i2, i3, lim)
hilbertsort!( directionx, !directiony, !directionz, coordinate, a, i3, i4, lim)
hilbertsort!( directionx, !directiony, !directionz, coordinate, a, i4, i5, lim)
hilbertsort!(!directiony, directionz, !directionx, next3d(coordinate), a, i5, i6, lim)
hilbertsort!(!directiony, directionz, !directionx, next3d(coordinate), a, i6, i7, lim)
hilbertsort!(!directionz, !directionx, directiony, nextnext3d(coordinate), a, i7, hi, lim)
return a
end
hilbertsort!{T<:AbstractPoint2D}(a::Array{T,1}) = hilbertsort!(backward, backward, coordinatey, a, 1, length(a))
hilbertsort!{T<:AbstractPoint2D}(a::Array{T,1}, lo::Int64, hi::Int64, lim::Int64) = hilbertsort!(backward, backward, coordinatey, a, lo, hi, lim)
hilbertsort!{T<:AbstractPoint3D}(a::Array{T,1}) = hilbertsort!(backward, backward, backward, coordinatez, a, 1, length(a))
hilbertsort!{T<:AbstractPoint3D}(a::Array{T,1}, lo::Int64, hi::Int64, lim::Int64) = hilbertsort!(backward, backward, backward, coordinatey, a, lo, hi, lim)
# multi-scale sort. Read all about it here:
# http://doc.cgal.org/latest/Spatial_sorting/classCGAL_1_1Multiscale__sort.html
function _mssort!{T<:AbstractPoint}(a::Array{T,1}, lim_ms::Int64, lim_hl::Int64, rat::Float64)
hi = length(a)
lo = 1
while true
lo = hi - int((1-rat)*hi)
hi-lo <= lim_ms && return a
hilbertsort!(a, lo, hi, lim_hl)
hi = lo-1
end
return a
end
# Utility methods, setting some different defaults for 2D and 3D. These are exported
mssort!{T<:AbstractPoint2D}(a::Array{T,1}; lim_ms::Int64=16, lim_hl::Int64=4, rat::Float64=0.25) =
_mssort!(a, lim_ms, lim_hl, rat)
mssort!{T<:AbstractPoint3D}(a::Array{T,1}; lim_ms::Int64=64, lim_hl::Int64=8, rat::Float64=0.125) =
_mssort!(a, lim_ms, lim_hl, rat)
###################################################################################
#
# Tests!
#
#
using Base.Test
function test()
# 2D orientation
ax = 1.1; ay = 1.1
bx = 1.2; by = 1.2
cx = 1.3; cy = 1.3
@test orientation(ax,ay,bx,by,cx,cy) == 0
@test orientation(bx,by,ax,ay,cx,cy) == 0
@test orientation(bx,by,cx,cy,ax,ay) == 0
ax = 1.1; ay = 1.1
bx = 1.2; by = 1.199999999999
cx = 1.3; cy = 1.3
@test orientation(ax,ay,bx,by,cx,cy) == 1
@test orientation(bx,by,ax,ay,cx,cy) == -1
@test orientation(bx,by,cx,cy,ax,ay) == 1
ax = 1.1; ay = 1.1
bx = 1.2; by = 1.199999999999999
cx = 1.3; cy = 1.3
@test orientation(ax,ay,bx,by,cx,cy) == 1
@test orientation(bx,by,ax,ay,cx,cy) == -1
@test orientation(bx,by,cx,cy,ax,ay) == 1
ax = 1.1; ay = 1.1
bx = 1.2; by = 1.21
cx = 1.3; cy = 1.3
@test orientation(ax,ay,bx,by,cx,cy) == -1
@test orientation(bx,by,ax,ay,cx,cy) == 1
@test orientation(bx,by,cx,cy,ax,ay) == -1
ax = 1.1; ay = 1.1
bx = 1.2; by = 1.20000000000001
cx = 1.3; cy = 1.3
@test orientation(ax,ay,bx,by,cx,cy) == -1
@test orientation(bx,by,ax,ay,cx,cy) == 1
@test orientation(bx,by,cx,cy,ax,ay) == -1
p1 = Point2D(1.0, 1.0)
p2 = Point2D(1.9, 1.5)
p3 = Point2D(1.45, 1.25)
tr = Triangle(p1, p2, p3)
@test orientation(tr) == 0
p3 = Point2D(getx(p3), 1.3)
tr = Triangle(p1, p2, p3)
@test orientation(tr) == 1
p3 = Point2D(getx(p3), 1.2)
tr = Triangle(p1, p2, p3)
@test orientation(tr) == -1
tr = Triangle(p1, p2, p3, positivelyoriented)
@test orientation(tr) == 1
tr = Triangle(p1, p3, p2, negativelyoriented)
@test orientation(tr) == -1
# 2D incircle
ax = 1.1; ay = 1.1
bx = 1.2; by = 1.2
cx = 1.3; cy = 1.3
dx = 1.4; dy = 1.7
@test incircle(ax,ay,bx,by,cx,cy,dx,dy) == 2
ax = 1.1; ay = 1.1
bx = 1.3; by = 1.1
cx = 1.3; cy = 1.3
dx = 1.1; dy = 1.3
@test incircle(ax,ay,bx,by,cx,cy,dx,dy) == 0
@test incircle(bx,by,ax,ay,cx,cy,dx,dy) == 0
@test incircle(bx,by,cx,cy,ax,ay,dx,dy) == 0
ax = 1.1; ay = 1.1
bx = 1.3; by = 1.1
cx = 1.3; cy = 1.3
dx = 1.1; dy = 1.3000001
@test incircle(ax,ay,bx,by,cx,cy,dx,dy) == -1
@test incircle(bx,by,ax,ay,cx,cy,dx,dy) == -1
@test incircle(bx,by,cx,cy,ax,ay,dx,dy) == -1
ax = 1.1; ay = 1.1
bx = 1.3; by = 1.1
cx = 1.3; cy = 1.3
dx = 1.1; dy = 1.30000000000001
@test incircle(ax,ay,bx,by,cx,cy,dx,dy) == -1
@test incircle(bx,by,ax,ay,cx,cy,dx,dy) == -1
@test incircle(bx,by,cx,cy,ax,ay,dx,dy) == -1
ax = 1.1; ay = 1.1
bx = 1.3; by = 1.1
cx = 1.3; cy = 1.3
dx = 1.1; dy = 1.29999
@test incircle(ax,ay,bx,by,cx,cy,dx,dy) == 1
@test incircle(bx,by,ax,ay,cx,cy,dx,dy) == 1
@test incircle(bx,by,cx,cy,ax,ay,dx,dy) == 1
ax = 1.1; ay = 1.1
bx = 1.3; by = 1.1
cx = 1.3; cy = 1.3
dx = 1.1; dy = 1.29999999999999
@test incircle(bx,by,ax,ay,cx,cy,dx,dy) == 1
@test incircle(bx,by,ax,ay,cx,cy,dx,dy) == 1
@test incircle(bx,by,cx,cy,ax,ay,dx,dy) == 1
p1 = Point2D(1.0, 1.0)
p2 = Point2D(1.0625, 1.0)
p3 = Point2D(1.0625, 1.0625)
tr = Triangle(p1, p2, p3)
p4 = Point2D(1.03, 1.003)
p5 = Point2D(1.99, 1.99)
p6 = Point2D(1.0, 1.0625)
@test incircle(tr, p4) == 1
@test incircle(tr, p5) == -1
@test incircle(tr, p6) == 0
@test incircle(tr, p2) == 0
@test incircle(tr, p3) == 0
@test incircle(tr, p1) == 0
# 3D orientation
ax = 1.1; ay = 1.1; az=1.1
bx = 1.2; by = 1.2; bz=1.1
cx = 1.3; cy = 1.3; cz=1.1
dx = 1.4; dy = 1.4; dz=1.1
@test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == 0
@test orientation(dx,dy,dz,bx,by,bz,cx,cy,cz,ax,ay,az) == 0
@test orientation(dx,dy,dz,cx,cy,cz,bx,by,bz,ax,ay,az) == 0
ax = 1.1; ay = 1.1; az=1.1
bx = 1.2; by = 1.2; bz=1.1
cx = 1.3; cy = 1.3; cz=1.1
dx = 1.4; dy = 1.4; dz=1.100001
@test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == 0
@test orientation(dx,dy,dz,bx,by,bz,cx,cy,cz,ax,ay,az) == 0
@test orientation(dx,dy,dz,cx,cy,cz,bx,by,bz,ax,ay,az) == 0
ax = 1.1; ay = 1.1; az=1.1
bx = 1.2; by = 1.2; bz=1.1
cx = 1.3; cy = 1.3; cz=1.1
dx = 1.4; dy = 1.4; dz=1.100000000001
@test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == 0
@test orientation(dx,dy,dz,bx,by,bz,cx,cy,cz,ax,ay,az) == 0
@test orientation(dx,dy,dz,cx,cy,cz,bx,by,bz,ax,ay,az) == 0
ax = 1.1; ay = 1.1; az=1.1
bx = 1.2; by = 1.2; bz=1.1
cx = 1.3; cy = 1.15; cz=1.1
dx = 1.4; dy = 1.17; dz=1.1
@test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == 0
@test orientation(dx,dy,dz,bx,by,bz,cx,cy,cz,ax,ay,az) == 0
@test orientation(dx,dy,dz,cx,cy,cz,bx,by,bz,ax,ay,az) == 0
ax = 1.1; ay = 1.1; az=1.1
bx = 1.3; by = 1.1; bz=1.1
cx = 1.3; cy = 1.3; cz=1.1
dx = 1.2; dy = 1.2; dz=1.5
@test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == 1
@test orientation(dx,dy,dz,bx,by,bz,cx,cy,cz,ax,ay,az) == -1
@test orientation(dx,dy,dz,cx,cy,cz,bx,by,bz,ax,ay,az) == 1
ax = 1.1; ay = 1.1; az=1.1
bx = 1.3; by = 1.1; bz=1.1
cx = 1.3; cy = 1.3; cz=1.1
dx = 1.2; dy = 1.2; dz=1.10000000000001
@test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == 1
@test orientation(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz) == -1
ax = 1.1; ay = 1.1; az=1.1
bx = 1.3; by = 1.1; bz=1.1
cx = 1.3; cy = 1.3; cz=1.1
dx = 1.2; dy = 1.2; dz=1.09
@test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == -1
@test orientation(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz) == 1
@test orientation(ax,ay,az,bx,by,bz,dx,dy,dz,cx,cy,cz) == 1
ax = 1.1; ay = 1.1; az=1.1
bx = 1.3; by = 1.1; bz=1.1
cx = 1.3; cy = 1.3; cz=1.1
dx = 1.2; dy = 1.2; dz=1.099999999999999
@test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == -1
@test orientation(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz) == 1
@test orientation(ax,ay,az,bx,by,bz,dx,dy,dz,cx,cy,cz) == 1
p1 = Point3D(ax, ay, az)
p2 = Point3D(bx, by, bz)
p3 = Point3D(cx, cy, cz)
p4 = Point3D(dx, dy, dz)
t = Tetrahedron(p1, p2, p3, p4)
@test orientation(t) == -1
t = Tetrahedron(p1, p2, p3, p4, positivelyoriented)
@test orientation(t) == 1
t = Tetrahedron(p1, p2, p4, p3, negativelyoriented)
@test orientation(t) == -1
# 3D incircle
ax = 1.1; ay = 1.1; az=1.1
bx = 1.3; by = 1.1; bz=1.1
cx = 1.3; cy = 1.3; cz=1.1
dx = 1.4; dy = 1.4; dz=1.2
ex = 1.1; ey = 1.1; ez=1.1
@test incircle(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,ex,ey,ez) == 0
@test incircle(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz,ex,ey,ez) == 0
@test incircle(ax,ay,az,bx,by,bz,dx,dy,dz,cx,cy,cz,ex,ey,ez) == 0
ax = 1.1; ay = 1.1; az=1.1
bx = 1.3; by = 1.1; bz=1.1
cx = 1.3; cy = 1.3; cz=1.1
dx = 1.2; dy = 1.2; dz=1.5
ex = 1.2; ey = 1.15; ez=1.13
@test incircle(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,ex,ey,ez) == 1
@test incircle(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz,ex,ey,ez) == 1
@test incircle(ax,ay,az,bx,by,bz,dx,dy,dz,cx,cy,cz,ex,ey,ez) == 1
ax = 1.1; ay = 1.1; az=1.1
bx = 1.3; by = 1.1; bz=1.1
cx = 1.3; cy = 1.3; cz=1.1
dx = 1.2; dy = 1.2; dz=1.5
ex = 1.2; ey = 1.15; ez=1.1000000000001
@test incircle(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,ex,ey,ez) == 1
@test incircle(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz,ex,ey,ez) == 1
@test incircle(ax,ay,az,bx,by,bz,dx,dy,dz,cx,cy,cz,ex,ey,ez) == 1
ax = 1.1; ay = 1.1; az=1.1
bx = 1.3; by = 1.1; bz=1.1
cx = 1.3; cy = 1.3; cz=1.1
dx = 1.2; dy = 1.2; dz=1.5
ex = 1.9; ey = 1.15; ez=1.01
@test incircle(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,ex,ey,ez) == -1
@test incircle(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz,ex,ey,ez) == -1
@test incircle(ax,ay,az,bx,by,bz,dx,dy,dz,cx,cy,cz,ex,ey,ez) == -1
# intriangle (2D)!
ax = 1.1; ay = 1.1
bx = 1.4; by = 1.1
cx = 1.1; cy = 1.4
dx = 1.11; dy = 1.11
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == 1
@test intriangle(bx,by,ax,ay,cx,cy,dx,dy) == 1
@test intriangle(bx,by,cx,cy,ax,ay,dx,dy) == 1
ax = 1.1; ay = 1.1
bx = 1.4; by = 1.1
cx = 1.1; cy = 1.4
dx = 1.01; dy = 1.01
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) < 0
@test intriangle(bx,by,ax,ay,cx,cy,dx,dy) < 0
@test intriangle(bx,by,cx,cy,ax,ay,dx,dy) < 0
ax = 1.1; ay = 1.1
bx = 1.4; by = 1.1
cx = 1.1; cy = 1.4
dx = 1.31; dy = 1.01
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) < 0
@test intriangle(bx,by,ax,ay,cx,cy,dx,dy) < 0
@test intriangle(bx,by,cx,cy,ax,ay,dx,dy) < 0
ax = 1.1; ay = 1.1
bx = 1.4; by = 1.1
cx = 1.1; cy = 1.4
dx = 1.1; dy = 1.2
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == 3
@test intriangle(bx,by,ax,ay,cx,cy,dx,dy) == 2
@test intriangle(bx,by,cx,cy,ax,ay,dx,dy) == 2
ax = 1.1; ay = 1.1
bx = 1.4; by = 1.1
cx = 1.1; cy = 1.4
dx = 1.2; dy = 1.1
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == 4
@test intriangle(bx,by,ax,ay,cx,cy,dx,dy) == 4
@test intriangle(bx,by,cx,cy,ax,ay,dx,dy) == 3
ax = 1.1; ay = 1.1
bx = 1.4; by = 1.1
cx = 1.1; cy = 1.4
dx = 1.25; dy = 1.25
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == 2
@test intriangle(bx,by,ax,ay,cx,cy,dx,dy) == 3
@test intriangle(bx,by,cx,cy,ax,ay,dx,dy) == 4
ax = 1.1; ay = 1.1
bx = 1.4; by = 1.1
cx = 1.1; cy = 1.4
dx = 1.25; dy = 1.250000000000001
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) < 0
@test intriangle(bx,by,ax,ay,cx,cy,dx,dy) < 0
@test intriangle(bx,by,cx,cy,ax,ay,dx,dy) < 0
ax = 1.1; ay = 1.1
bx = 1.4; by = 1.1
cx = 1.1; cy = 1.4
dx = 1.25; dy = 1.249999999999999
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == 1
@test intriangle(bx,by,ax,ay,cx,cy,dx,dy) == 1
@test intriangle(bx,by,cx,cy,ax,ay,dx,dy) == 1
ax = 1.1; ay = 1.1
bx = 1.4; by = 1.1
cx = 1.1; cy = 1.4
dx = 1.7; dy = 1.1
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -1
dx = 1.1; dy = 1.7
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -1
dx = 1.01; dy = 1.1
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -2
dx = 1.1; dy = 1.01
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -3
ax = 1.1; ay = 1.1
bx = 1.4; by = 1.1
cx = 1.1; cy = 1.4
dx = 1.55; dy = 1.55
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -1
dx = 1.09; dy = 1.11
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -2
dx = 1.11; dy = 1.09
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -3
ax = min_coord; ay = min_coord
bx = 1.1; by = 1.1
cx = max_coord; cy = min_coord
dx = 1.2; dy = 1.2
@test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -1
# intriangle (3D) !
ax = 1.1; ay = 1.1; az = 1.1
bx = 1.4; by = 1.1; bz = 1.1
cx = 1.1; cy = 1.4; cz = 1.1
dx = 1.1; dy = 1.1; dz = 1.4
px = 1.11; py = 1.11; pz = 1.11
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == 1
@test intriangle(bx, by, bz, ax, ay, az, cx, cy, cz, dx, dy, dz, px, py, pz) == 1
@test intriangle(ax, ay, az, cx, cy, cz, bx, by, bz, dx, dy, dz, px, py, pz) == 1
ax = 1.1; ay = 1.1; az = 1.1
bx = 1.4; by = 1.1; bz = 1.1
cx = 1.1; cy = 1.4; cz = 1.1
dx = 1.1; dy = 1.1; dz = 1.4
px = 1.02; py = 1.2; pz = 1.2
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) < 0
@test intriangle(bx, by, bz, ax, ay, az, cx, cy, cz, dx, dy, dz, px, py, pz) < 0
@test intriangle(ax, ay, az, cx, cy, cz, bx, by, bz, dx, dy, dz, px, py, pz) < 0
ax = 1.1; ay = 1.1; az = 1.1
bx = 1.4; by = 1.1; bz = 1.1
cx = 1.1; cy = 1.4; cz = 1.1
dx = 1.1; dy = 1.1; dz = 1.4
px = 1.09999999999999; py = 1.2; pz = 1.2
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) < 0
@test intriangle(bx, by, bz, ax, ay, az, cx, cy, cz, dx, dy, dz, px, py, pz) < 0
@test intriangle(ax, ay, az, cx, cy, cz, bx, by, bz, dx, dy, dz, px, py, pz) < 0
ax = 1.1; ay = 1.1; az = 1.1
bx = 1.4; by = 1.1; bz = 1.1
cx = 1.1; cy = 1.4; cz = 1.1
dx = 1.1; dy = 1.1; dz = 1.4
px = 1.10000000000001; py = 1.2; pz = 1.2
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == 1
@test intriangle(bx, by, bz, ax, ay, az, cx, cy, cz, dx, dy, dz, px, py, pz) == 1
@test intriangle(ax, ay, az, cx, cy, cz, bx, by, bz, dx, dy, dz, px, py, pz) == 1
ax = 1.1; ay = 1.1; az = 1.1
bx = 1.4; by = 1.1; bz = 1.1
cx = 1.1; cy = 1.4; cz = 1.1
dx = 1.1; dy = 1.1; dz = 1.4
px = 1.11; py = 1.11; pz = 1.1
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == 5
@test intriangle(bx, by, bz, ax, ay, az, cx, cy, cz, dx, dy, dz, px, py, pz) == 5
@test intriangle(ax, ay, az, cx, cy, cz, bx, by, bz, dx, dy, dz, px, py, pz) == 5
@test intriangle(dx, dy, dz, bx, by, bz, cx, cy, cz, ax, ay, az, px, py, pz) == 2
@test intriangle(ax, ay, az, dx, dy, dz, cx, cy, cz, bx, by, bz, px, py, pz) == 3
@test intriangle(ax, ay, az, bx, by, bz, dx, dy, dz, cx, cy, cz, px, py, pz) == 4
ax = 1.1; ay = 1.1; az = 1.1
bx = 1.4; by = 1.1; bz = 1.1
cx = 1.1; cy = 1.4; cz = 1.1
dx = 1.1; dy = 1.1; dz = 1.4
px = 1.55; py = 1.55; pz = 1.55
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
px = 1.09; py = 1.11; pz = 1.11
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -2
px = 1.11; py = 1.09; pz = 1.11
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -3
px = 1.11; py = 1.11; pz = 1.09
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -4
ax = 1.1; ay = 1.1; az = 1.1
bx = 1.4; by = 1.1; bz = 1.1
cx = 1.1; cy = 1.4; cz = 1.1
dx = 1.1; dy = 1.1; dz = 1.4
px = 1.7; py = 1.1; pz = 1.1
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
px = 1.7; py = 1.1; pz = 1.1
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
px = 1.1; py = 1.7; pz = 1.1
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
px = 1.1; py = 1.1; pz = 1.7
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
px = 1.9; py = 1.9; pz = 1.1
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
px = 1.9; py = 1.1; pz = 1.9
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
px = 1.1; py = 1.9; pz = 1.9
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
px = 1.01; py = 1.1; pz = 1.1
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -2
px = 1.1; py = 1.01; pz = 1.1
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -3
px = 1.1; py = 1.1; pz = 1.01
@test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -4
# Scale dependednt Peano-Hilbert stuff
# 2D
@test peanokey(Point2D(1.01,1.01), 2) == 0
@test peanokey(Point2D(1.9901,1.01), 2) == 4*4-1
@test Point2D(15, 2) == Point2D(1.75, 1.0)
@test Point2D(0, 2) == Point2D(1.0, 1.0)
for x in linspace(1.0,1.999999,100), y in linspace(1.0,1.999999,100)
p = Point2D(x, y)
d = peanokey(p)
pp= Point2D(d)
@test abs(getx(p) - getx(pp)) < 1e-7
@test abs(gety(p) - gety(pp)) < 1e-7
end
# 3D
@test peanokey(Point3D(1.01,1.01,1.01), 2) == 0
@test peanokey(Point3D(1.01,1.01,1.9901), 2) == 4*4*4-1
@test Point3D(15, 2) == Point3D(1.25,1.5,1.0)
@test Point3D(0, 2) == Point3D(1.0, 1.0, 1.0)
for x in linspace(1.0,1.99999,30), y in linspace(1.0,1.99999,30), z in linspace(1.0,1.999999,30)
p = Point3D(x, y, z)
d = peanokey(p)
pp= Point3D(d)
@test abs(getx(p) - getx(pp)) < 1e-5
@test abs(gety(p) - gety(pp)) < 1e-5
@test abs(getz(p) - getz(pp)) < 1e-5
end
# Scale independent Peano-Hilbert stuff
a2 = Point2D[]
N=1000
for i in 1:N
push!(a2, Point2D(1.0+rand(), 1.0+rand()))
end
# TODO: not much of a test, I need to betteer think how to test this
# At least this runs the code...
@test length(mssort!(a2)) == N
a3 = Point3D[]
for i in 1:N
push!(a3, Point3D(1.0+rand(), 1.0+rand(), 1.0+rand()))
end
# TODO: not much of a test, I need to betteer think how to test this
# At least this runs the code...
@test length(mssort!(a3)) == N
# Qttys functions -- 2D
a = Point2D(1.0, 1.0)
b = Point2D(1.0, 1.5)
c = Point2D(1.5, 1.0)
t = Triangle(a, b, c)
@test abs(area(t)-1/8) < 1e-7
c = centroid(t)
@test abs(getx(c)-1-1/6) < 1e-7
@test abs(gety(c)-1-1/6) < 1e-7
c = circumcenter(t)
@test abs(getx(c)-1.25) < 1e-7
@test abs(gety(c)-1.25) < 1e-7
@test abs(circumradius2(t)-1/8) < 1e-7
# Qttys functions -- 3D
a = Point3D(1.0, 1.0, 1.0)
b = Point3D(1.0, 1.0, 1.5)
c = Point3D(1.0, 1.5, 1.0)
d = Point3D(1.5, 1.0, 1.0)
t = Tetrahedron(a, b, c, d)
@test abs(volume(t)-1/16) < 1e-7
c = centroid(t)
@test abs(getx(c)-1.125) < 1e-7
@test abs(gety(c)-1.125) < 1e-7
@test abs(getz(c)-1.125) < 1e-7
c = circumcenter(t)
@test abs(getx(c)-1.25) < 1e-7
@test abs(gety(c)-1.25) < 1e-7
@test abs(getz(c)-1.25) < 1e-7
@test abs(circumradius2(t)-0.25^2*3) < 1e-7
# that's it for today!
println("gpred: All tests passed!")
end
test()
end # module GeometicalPredicates
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment