Skip to content
This repository has been archived by the owner on Oct 23, 2022. It is now read-only.

Added a new "Surface" type: ConvexPolygon which support planner simple convex polygons #213

Merged
merged 7 commits into from
Jul 9, 2021
1 change: 0 additions & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@ FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
GeometricalPredicates = "fd0ad045-b25c-564e-8f9c-8ef5c5f21267"
HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
ImageView = "86fae568-95e7-573e-a6b2-d8a6b900c9ef"
Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
Expand Down
72 changes: 52 additions & 20 deletions src/Geometry/Primitives/NonCSG/ConvexPolygon.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,8 @@
# See LICENSE in the project root for full license information.

import Statistics
import GeometricalPredicates

const _abs_err_orientation_2d = 2*eps(Float64)

"""
ConvexPolygon{T} <: Surface{T}
Expand All @@ -13,50 +14,62 @@ The rotation of the polygon around its normal is defined by `rotationvec`.
`rotationvec×surfacenormal` is taken as the vector along the u axis.

```julia
ConvexPolygon(local_frame::Transform{T}, local_polygon_points::Vector{Vector{T}}, interface::NullOrFresnel{T} = nullinterface(T))
ConvexPolygon(local_frame::Transform{T}, local_polygon_points::SVector{N, SVector{2, T}}, interface::NullOrFresnel{T} = nullinterface(T))
```

The local frame defines the plane (spans by the right and up vectors) with the plane normal given by the forward vector.
the local_polygon_points are given with respect to the local frame and are 2D points.
NOTE: This class uses static vectors to hold the points which will lead to more efficient performance, but should not be used with polygons with more than 20-30 points.
"""
struct ConvexPolygon{T} <: Surface{T}
struct ConvexPolygon{N, T<:Real} <: Surface{T}
plane::Plane{T,3}
local_frame::Transform{T}
local_points::Vector{Vector{T}}
local_points::SVector{N, SVector{2, T}}

# for efficency
_local_frame_inv::Transform{T} # cache the inverse matrix to avoid computing it for every intersection test
_poly2d::GeometricalPredicates.Polygon2D{GeometricalPredicates.Point2D}
_local_frame_inv::Transform{T} # cache the inverse matrix to avoid computing it for every intersection test
_local_lines::SVector{N, SVector{3, SVector{2, T}}} # defines the edges + a third point representing the slopes in order to save some calculationsduring ray checking
_length::Int64
# _poly2d::GeometricalPredicates.Polygon2D{GeometricalPredicates.Point2D}

function ConvexPolygon(
local_frame::Transform{T},
local_polygon_points::Vector{SVector{2, T}},
local_polygon_points::SVector{N, SVector{2, T}},
interface::NullOrFresnel{T} = NullInterface(T)
) where {T<:Real}
) where {N, T<:Real}

# need at least 3 points to define apolygon
@assert length(local_polygon_points) > 3

local_center = Statistics.mean(local_polygon_points)
world_center = local2world(local_frame) * Vec3(local_center[1], local_center[2], zero(T))

poly_points = [GeometricalPredicates.Point2D(p[1], p[2]) for p in local_polygon_points]
poly = GeometricalPredicates.Polygon2D(poly_points...)

# poly_points = [GeometricalPredicates.Point2D(p[1], p[2]) for p in local_polygon_points]
# poly = GeometricalPredicates.Polygon2D(poly_points...)
pts = local_polygon_points
local_lines = SVector(
[SVector(
pts[i], # A
pts[mod(i, length(pts))+1], # B
# bx = Bx - Ax, by = By - Ay
SVector(pts[mod(i, length(pts))+1][1] - pts[i][1], pts[mod(i, length(pts))+1][2] - pts[i][2]))
for i in 1:length(pts)]...
)

plane = Plane(forward(local_frame), world_center, interface = interface)
new{T}(plane, local_frame, local_polygon_points, inv(local_frame), poly)
new{N, T}(plane, local_frame, local_polygon_points, inv(local_frame), local_lines, length(local_lines))
end
end
export ConvexPolygon

# Base.show(io::IO, poly::ConvexPolygon{T}) where {T<:Real} = print(io, "ConvexPolygon{$T}($(centroid(hex)), $(normal(hex)), $(hex.side_length), $(interface(hex)))")
centroid(poly::ConvexPolygon{T}) where {T<:Real} = poly.plane.pointonplane
interface(poly::ConvexPolygon{T}) where {T<:Real} = interface(poly.plane)
normal(poly::ConvexPolygon{T}) where {T<:Real} = normal(poly.plane)
normal(poly::ConvexPolygon{T}, ::T, ::T) where {T<:Real} = normal(poly)

centroid(poly::ConvexPolygon{N, T}) where {N, T<:Real} = poly.plane.pointonplane
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since N and T are not being used on the right hand side of the = you shouldn't need them. Should still generate efficient code if you remove them.

interface(poly::ConvexPolygon{N, T}) where {N, T<:Real} = interface(poly.plane)
normal(poly::ConvexPolygon{N, T}) where {N, T<:Real} = normal(poly.plane)
normal(poly::ConvexPolygon{N, T}, ::T, ::T) where {N, T<:Real} = normal(poly)

function surfaceintersection(poly::ConvexPolygon{T}, r::AbstractRay{T,3}) where {T<:Real}
function surfaceintersection(poly::ConvexPolygon{N, T}, r::AbstractRay{T,3}) where {N, T<:Real}
interval = surfaceintersection(poly.plane, r)
if interval isa EmptyInterval{T} || isinfiniteinterval(interval)
return EmptyInterval(T) # no ray plane intersection or inside plane but no hit
Expand All @@ -66,9 +79,28 @@ function surfaceintersection(poly::ConvexPolygon{T}, r::AbstractRay{T,3}) where

local_p = poly._local_frame_inv * p

in_poly = GeometricalPredicates.inpolygon(poly._poly2d, GeometricalPredicates.Point(local_p[1], local_p[2]))
@inline function orientation(l::SVector{3, SVector{2, T}})::Int8 where {T<:Real}
cx = local_p[1] - l[1][1] # getx(p) - getx(geta(l))
cy = local_p[2] - l[1][2] # gety(p) - gety(geta(l))
_pr2 = -(l[3][1]*cy) + l[3][2]*cx # _pr2 = -l._bx*cy + l._by*cx
if _pr2 < -_abs_err_orientation_2d
1
elseif _pr2 > _abs_err_orientation_2d
-1
else
0 # can implement something more accurate in the future to minimize numerical errors
end
end

@inline function inpolygon()
side = orientation(poly._local_lines[1])
for i = 2:poly._length
orientation(poly._local_lines[i]) == side || return false
end
return true
end

if !in_poly
if !inpolygon()
return EmptyInterval(T)
else
if dot(normal(poly), direction(r)) < zero(T)
Expand All @@ -86,7 +118,7 @@ end

Create a triangle mesh that can be rendered by iterating on the polygon's edges and for each edge use the centroid as the third vertex of the triangle.
"""
function makemesh(poly::ConvexPolygon{T}, ::Int = 0) where {T<:Real}
function makemesh(poly::ConvexPolygon{N, T}, ::Int = 0) where {N, T<:Real}
c = centroid(poly)

l2w = local2world(poly.local_frame)
Expand Down
6 changes: 3 additions & 3 deletions src/Optical/Paraxial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,10 +17,10 @@ ParaxialLensConvexPoly(focaldistance, local_frame, local_polygon_points, local_c
```
"""
struct ParaxialLens{T} <: Surface{T}
shape::Union{Rectangle{T},Ellipse{T},Hexagon{T},ConvexPolygon{T}}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this will do what you want. I think this will define a ConvexPolygon whose first type argument is T. This might be what you want:

shape::Union{Rectangle{T},Ellipse{T},Hexagon{T},ConvexPolygon{N,T}} where{N}

Then for sure the T type argument to ConvexPolygon is in the right place.

shape::Union{Rectangle{T},Ellipse{T},Hexagon{T},ConvexPolygon}
interface::ParaxialInterface{T}

function ParaxialLens(shape::Union{Rectangle{T},Ellipse{T},Hexagon{T},ConvexPolygon{T}}, interface::ParaxialInterface{T}) where {T<:Real}
function ParaxialLens(shape::Union{Rectangle{T},Ellipse{T},Hexagon{T},ConvexPolygon}, interface::ParaxialInterface{T}) where {T<:Real}
new{T}(shape, interface)
end
end
Expand All @@ -47,7 +47,7 @@ function ParaxialLensHex(focaldistance::T, side_length::T, surfacenormal::SVecto
return ParaxialLens(h, ParaxialInterface(focaldistance, centrepoint, outsidematerial))
end

function ParaxialLensConvexPoly(focaldistance::T, local_frame::Transform{T}, local_polygon_points::Vector{SVector{2, T}}, local_center_point::SVector{2, T}; outsidematerial::OpticSim.GlassCat.AbstractGlass = OpticSim.GlassCat.Air) where {T<:Real}
function ParaxialLensConvexPoly(focaldistance::T, local_frame::Transform{T}, local_polygon_points::SVector{N, SVector{2, T}}, local_center_point::SVector{2, T}; outsidematerial::OpticSim.GlassCat.AbstractGlass = OpticSim.GlassCat.Air) where {N, T<:Real}
poly = ConvexPolygon(local_frame, local_polygon_points)
centrepoint = SVector{3, T}(local2world(local_frame) * Vec3(local_center_point[1], local_center_point[2], zero(T)))
return ParaxialLens(poly, ParaxialInterface(focaldistance, centrepoint, outsidematerial))
Expand Down
4 changes: 2 additions & 2 deletions test/testsets/Intersection.jl
Original file line number Diff line number Diff line change
Expand Up @@ -397,14 +397,14 @@

t = identitytransform()

local_points_2d = [
local_points_2d = SVector(
SVector(-1.0, -1.0),
SVector(1.0, -1.0),
SVector(2.0, 0.0),
SVector(1.0, 1.0),
SVector(-1.0, 1.0),
SVector(-2.0, 0.0)
]
)

@test_nowarn ConvexPolygon(t, local_points_2d)

Expand Down