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

Project.toml

@@ -44,6 +44,7 @@ ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StringEncodings = "69024149-9ee7-55f6-a4c4-859efe599b68"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008"

src/Geometry/Geometry.jl

@@ -34,6 +34,7 @@ include("Primitives/SphericalCap.jl")
 
 include("Primitives/NonCSG/Rectangle.jl")
 include("Primitives/NonCSG/Hexagon.jl")
+include("Primitives/NonCSG/ConvexPolygon.jl")
 include("Primitives/NonCSG/Ellipse.jl")
 include("Primitives/NonCSG/Stop.jl")
 

src/Geometry/Primitives/NonCSG/ConvexPolygon.jl

@@ -0,0 +1,139 @@
+# MIT license
+# Copyright (c) Microsoft Corporation. All rights reserved.
+# See LICENSE in the project root for full license information.
+
+import Statistics
+
+const _abs_err_orientation_2d = 2*eps(Float64)
+
+"""
+    ConvexPolygon{T} <: Surface{T}
+
+General Convex Polygon surface, not a valid CSG object.
+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{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<:Real}  <: Surface{T} 
+    plane::Plane{T,3}
+    local_frame::Transform{T}
+    local_points::Vector{SVector{2, T}}
+
+    # for efficency
+    _local_frame_inv::Transform{T}                                  # cache the inverse matrix to avoid computing it for every intersection test
+    _local_lines::Vector{SVector{3, SVector{2, T}}}                 # defines the edge points + a third point representing the slopes in order to save some calculationsduring ray checking
+    _length::Int64                                                  # cache the length of the polygon 
+
+    function ConvexPolygon(
+            local_frame::Transform{T},
+            local_polygon_points::Vector{SVector{2, T}}, 
+            interface::NullOrFresnel{T} = NullInterface(T)
+        ) where {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...)
+        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), 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) = poly.plane.pointonplane
+interface(poly::ConvexPolygon) = interface(poly.plane)
+normal(poly::ConvexPolygon) = normal(poly.plane)
+
+function surfaceintersection(poly::ConvexPolygon{T}, r::AbstractRay{T,3}) where {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
+    else
+        intersect = halfspaceintersection(interval)
+        p = point(intersect)
+
+        local_p = poly._local_frame_inv * p
+
+        @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 !inpolygon()
+            return EmptyInterval(T)
+        else
+            if dot(normal(poly), direction(r)) < zero(T)
+                return positivehalfspace(intersect)
+            else
+                return rayorigininterval(intersect)
+            end
+        end
+    end
+end
+
+
+"""
+    makemesh(poly::ConvexPolygon{T}, ::Int = 0) where {T<:Real} -> TriangleMesh
+
+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}
+    c = centroid(poly)
+
+    l2w = local2world(poly.local_frame)
+    len = length(poly.local_points)
+
+    triangles = []
+    for i in 1:len
+        p1 = poly.local_points[i]
+        p2 = poly.local_points[mod(i,len) + 1]
+
+        tri = Triangle(
+            Vector(l2w * Vec3(p2[1], p2[2], zero(T))), 
+            Vector(c), 
+            Vector(l2w * Vec3(p1[1], p1[2], zero(T))))
+        push!(triangles, tri)
+    end
+
+    triangles = Vector{Triangle{T}}(triangles)
+
+    return TriangleMesh(triangles)
+end

src/Geometry/Primitives/NonCSG/Ellipse.jl

@@ -53,7 +53,6 @@ interface(a::Ellipse{T}) where {T<:Real} = interface(a.plane)
 uvrange(::Type{Ellipse{T}}) where {T<:Real} = ((-T(π), T(π)), (zero(T), one(T))) # θ and ρ
 
 normal(r::Ellipse{T}) where {T<:Real} = normal(r.plane)
-normal(r::Ellipse{T}, ::T, ::T) where {T<:Real} = normal(r)
 
 point(r::Ellipse{T}, θ::T, ρ::T) where {T<:Real} = centroid(r) + ρ * (r.halfsizeu * cos(θ) * r.uvec + r.halfsizev * sin(θ) * r.vvec)
 partials(r::Ellipse{T}, θ::T, ρ::T) where {T<:Real} = (ρ * (r.halfsizeu * -sin(θ) * r.uvec + r.halfsizev * cos(θ) * r.vvec), (r.halfsizeu * cos(θ) * r.uvec + r.halfsizev * sin(θ) * r.vvec))

src/Geometry/Primitives/NonCSG/Hexagon.jl

@@ -48,7 +48,6 @@ Base.show(io::IO, hex::Hexagon{T}) where {T<:Real} = print(io, "Hexagon{$T}($(ce
 centroid(hex::Hexagon{T}) where {T<:Real} = hex.plane.pointonplane
 interface(hex::Hexagon{T}) where {T<:Real} = interface(hex.plane)
 normal(hex::Hexagon{T}) where {T<:Real} = normal(hex.plane)
-normal(hex::Hexagon{T}, ::T, ::T) where {T<:Real} = normal(hex)
 
 function surfaceintersection(hex::Hexagon{T}, r::AbstractRay{T,3}) where {T<:Real}
     interval = surfaceintersection(hex.plane, r)

src/Geometry/Primitives/NonCSG/Rectangle.jl

@@ -59,7 +59,6 @@ interface(a::Rectangle{T}) where {T<:Real} = interface(a.plane)
 uvrange(::Type{Rectangle{T}}) where {T<:Real} = ((-one(T), one(T)), (-one(T), one(T)))
 
 normal(r::Rectangle{T}) where {T<:Real} = normal(r.plane)
-normal(r::Rectangle{T}, ::T, ::T) where {T<:Real} = normal(r)
 
 point(r::Rectangle{T}, u::T, v::T) where {T<:Real} = centroid(r) + (r.halfsizeu * u * r.uvec) + (r.halfsizev * v * r.vvec)
 partials(r::Rectangle{T}, ::T, ::T) where {T<:Real} = r.halfsizeu * r.uvec, r.halfsizev * r.vvec

src/Geometry/Transform.jl

@@ -5,6 +5,7 @@
 # -----------------------------------------------------------------------------------------------
 # GEOMETRY
 # -----------------------------------------------------------------------------------------------
+using ...OpticSim
 using StaticArrays
 using LinearAlgebra
 
@@ -233,6 +234,40 @@ function Transform(rotation::AbstractArray{T,2}, translation::AbstractArray{T,1}
         translation[1], translation[2], translation[3], one(T))
 end
 
+
+# define some utility functions 
+
+"""
+    right(t::Transform{<:Real}) -> Vec3
+
+Assuming t is a 3D rigid transform representing a local left-handed coordinate system, this function will return the first column, representing the "X" axis.
+"""
+right(t::Transform{<:Real}) = normalize(Vec3(t[1,1], t[2,1], t[3,1]))
+
+"""
+    up(t::Transform{<:Real}) -> Vec3
+
+Assuming t is a 3D rigid transform representing a local left-handed coordinate system, this function will return the second column, representing the "Y" axis.
+"""
+up(t::Transform{<:Real}) = normalize(Vec3(t[1,2], t[2,2], t[3,2]))
+
+"""
+    forward(t::Transform{<:Real}) -> Vec3
+
+Assuming t is a 3D rigid transform representing a local left-handed coordinate system, this function will return the third column, representing the "Z" axis.
+"""
+forward(t::Transform{<:Real}) = normalize(Vec3(t[1,3], t[2,3], t[3,3]))
+
+"""
+    origin(t::Transform{<:Real}) -> Vec3
+
+Assuming t is a 3D rigid transform representing a local left-handed coordinate system, this function will return the fourth column, containing the translation part of the transform in 3D.
+"""
+OpticSim.origin(t::Transform{<:Real}) = Vec3(t[1,4], t[2,4], t[3,4])
+
+export right, up, forward, origin
+
+
 """
     rotationX(angle::T) where {T<:Real} -> Transform
 

src/Optical/Emitters/Emitters.jl

@@ -18,12 +18,6 @@ using LinearAlgebra
 import Unitful.Length
 using Unitful.DefaultSymbols
 
-# define some utility functions 
-right(t::OpticSim.Geometry.Transform{<:Real}) = normalize(Vec3(t[1,1], t[2,1], t[3,1]))
-up(t::OpticSim.Geometry.Transform{<:Real}) = normalize(Vec3(t[1,2], t[2,2], t[3,2]))
-forward(t::OpticSim.Geometry.Transform{<:Real}) = normalize(Vec3(t[1,3], t[2,3], t[3,3]))
-OpticSim.origin(t::OpticSim.Geometry.Transform{<:Real}) = Vec3(t[1,4], t[2,4], t[3,4])
-
 # defining name placeholders to override in nested modules
 
 """

src/Optical/Paraxial.jl

@@ -13,14 +13,15 @@ Create with the following functions
 ParaxialLensEllipse(focaldistance, halfsizeu, halfsizev, surfacenormal, centrepoint; rotationvec = [0.0, 1.0, 0.0], outsidematerial = OpticSim.GlassCat.Air, decenteruv = (0.0, 0.0))
 ParaxialLensRect(focaldistance, halfsizeu, halfsizev, surfacenormal, centrepoint; rotationvec = [0.0, 1.0, 0.0], outsidematerial = OpticSim.GlassCat.Air, decenteruv = (0.0, 0.0))
 ParaxialLensHex(focaldistance, side_length, surfacenormal, centrepoint; rotationvec = [0.0, 1.0, 0.0], outsidematerial = OpticSim.GlassCat.Air, decenteruv = (0.0, 0.0))
+ParaxialLensConvexPoly(focaldistance, local_frame, local_polygon_points, local_center_point; outsidematerial = OpticSim.GlassCat.Air)
 ```
 """
 struct ParaxialLens{T} <: Surface{T}
-    shape::Union{Rectangle{T},Ellipse{T},Hexagon{T}}
+    shape::Union{Rectangle{T},Ellipse{T},Hexagon{T},ConvexPolygon{T}}
     interface::ParaxialInterface{T}
 
-    function ParaxialLens(shape::Union{Rectangle{T},Ellipse{T},Hexagon{T}}, interface::ParaxialInterface{T}) where {T<:Real}
-        new{T}(shape, interface)
+    function ParaxialLens(shape::Union{Rectangle{T},Ellipse{T},Hexagon{T},ConvexPolygon{T}}, interface::ParaxialInterface{T}) where {T<:Real}
+            new{T}(shape, interface)
     end
 end
 
@@ -46,6 +47,12 @@ 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 {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))
+end
+
 function ParaxialLensEllipse(focaldistance::T, halfsizeu::T, halfsizev::T, surfacenormal::AbstractArray{T,1}, centrepoint::AbstractArray{T,1}; rotationvec::AbstractArray{T,1} = SVector{3,T}(0.0, 1.0, 0.0), outsidematerial::OpticSim.GlassCat.AbstractGlass = OpticSim.GlassCat.Air, decenteruv::Tuple{T,T} = (zero(T), zero(T))) where {T<:Real}
     @assert length(surfacenormal) == 3 && length(centrepoint) == 3
     return ParaxialLensEllipse(focaldistance, halfsizeu, halfsizev, SVector{3,T}(surfacenormal), SVector{3,T}(centrepoint), rotationvec = SVector{3,T}(rotationvec), outsidematerial = outsidematerial, decenteruv = decenteruv)
@@ -57,12 +64,11 @@ function ParaxialLensEllipse(focaldistance::T, halfsizeu::T, halfsizev::T, surfa
     return ParaxialLens(e, ParaxialInterface(focaldistance, centrepoint, outsidematerial))
 end
 
-export ParaxialLens, ParaxialLensRect, ParaxialLensEllipse, ParaxialLensHex
+export ParaxialLens, ParaxialLensRect, ParaxialLensEllipse, ParaxialLensHex, ParaxialLensConvexPoly
 
 interface(r::ParaxialLens{T}) where {T<:Real} = r.interface
 centroid(r::ParaxialLens{T}) where {T<:Real} = centroid(r.shape)
 normal(r::ParaxialLens{T}) where {T<:Real} = normal(r.shape)
-normal(r::ParaxialLens{T}, ::T, ::T) where {T<:Real} = normal(r)
 point(r::ParaxialLens{T}, u::T, v::T) where {T<:Real} = point(r.shape, u, v)
 uv(r::ParaxialLens{T}, x::T, y::T, z::T) where {T<:Real} = uv(r, SVector{3,T}(x, y, z))
 uv(r::ParaxialLens{T}, p::SVector{3,T}) where {T<:Real} = uv(r.shape, p)

src/utilities.jl

@@ -71,3 +71,8 @@ end
         return acos(x)
     end
 end
+
+
+# some place holders for package level function names.
+# these names need to exist before any internal module can override them.
+function origin end    

test/testsets/Intersection.jl

@@ -392,6 +392,77 @@
         @test surfaceintersection(hex2, raymiss) isa EmptyInterval
     end # testset Hexagon
 
+
+    @testset "Convex Polygon" begin
+
+        t = identitytransform()
+
+        local_points_2d = [
+            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)
+
+        rinplane = Ray([0.0, 0.0, 0.0], [1.0, 1.0, 0.0])
+        routside = Ray([0.0, 0.0, 1.0], [1.0, 1.0, 1.0])
+        rinside = Ray([0.0, 0.0, -1.0], [1.0, 1.0, -1.0])
+        routintersects = Ray([0.0, 0.0, 1.0], [0.0, 0.0, -1.0])
+        rinintersects = Ray([0.0, 0.0, -1.0], [0.0, 0.0, 1.0])
+        routmiss = Ray([0.0, 2.0, 1.0], [0.0, 0.0, -1.0])
+        rinmiss = Ray([0.0, 2.0, -1.0], [0.0, 0.0, 1.0])
+        rinbounds = Ray([sqrt(3) / 2, 0.0, -1.0], [0.0, 0.0, 1.0])
+        routbounds = Ray([sqrt(3) / 2, 0.0, 1.0], [0.0, 0.0, -1.0])
+        rasymhit = Ray([0.0, 0.9, 1.0], [0.0, 0.0, -1.0])
+        rasymmiss = Ray([3.0, 0.0, 1.0], [0.0, 0.0, -1.0])
+
+        poly = ConvexPolygon(t, local_points_2d)
+
+        # parallel to plane and outside
+        @test surfaceintersection(poly, routside) isa EmptyInterval
+
+        # parallel and on plane
+        @test surfaceintersection(poly, rinplane) isa EmptyInterval
+
+        # inside but not hitting
+        @test surfaceintersection(poly, rinside) isa EmptyInterval
+
+        # starts outside and hits
+        res = surfaceintersection(poly, routintersects)
+        @test isapprox(point(lower(res)), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && OpticSim.upper(res) isa Infinity
+
+        # starts inside and hits
+        res = surfaceintersection(poly, rinintersects)
+        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)
+
+        # starts inside and misses bounds
+        @test surfaceintersection(poly, rinmiss) isa EmptyInterval
+
+        # starts outside and misses bounds
+        @test surfaceintersection(poly, routmiss) isa EmptyInterval
+
+        # starts outside and hits bounds
+        res = surfaceintersection(poly, routbounds)
+        @test isapprox(point(lower(res)), [sqrt(3) / 2, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity
+
+        # starts inside and hits bounds
+        res = surfaceintersection(poly, rinbounds)
+        @test lower(res) isa RayOrigin && isapprox(point(upper(res)), [sqrt(3) / 2, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE)
+
+        # asymmetric hit
+        res = surfaceintersection(poly, rasymhit)
+        @test isapprox(point(lower(res)), [0.0, 0.9, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) && upper(res) isa Infinity
+
+        # asymmetric miss
+        @test surfaceintersection(poly, rasymmiss) isa EmptyInterval
+
+    end # testset Convex Polygon
+
+
     @testset "Stops" begin
         infiniterect = RectangularAperture(0.4, 0.8, SVector(0.0, 1.0, 0.0), SVector(0.0, 0.0, 1.0))
         finiterect = RectangularAperture(0.4, 0.8, 1.0, 2.0, SVector(0.0, 1.0, 0.0), SVector(0.0, 0.0, 1.0))