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
2 changes: 2 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ 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 All @@ -43,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"
Expand Down
1 change: 1 addition & 0 deletions src/Geometry/Geometry.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")

Expand Down
104 changes: 104 additions & 0 deletions src/Geometry/Primitives/NonCSG/ConvexPolygon.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,104 @@
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

import Statistics
import GeometricalPredicates

"""
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{Vector{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.
"""
struct ConvexPolygon{T} <: Surface{T}
plane::Plane{T,3}
local_frame::Transform{T}
local_points::Vector{Vector{T}}
Copy link
Contributor

Choose a reason for hiding this comment

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

local_points should probably be SVector{SVector}. Vectors will be allocated on the heap which will lead to allocations and will probably also be slower. If we expected large convex polygons then Vector might be appropriate. But this class is going to be used for polygons with less than 20-30 points. If you change to SVector probably should document that should only have small number of points.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I replaced it with Vector{SVector{2, T}}.
I can't find a convenient way to define an static array of static arrays. Please send me an example if you have one.

Copy link
Contributor

Choose a reason for hiding this comment

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

you'd have to do something like this:

struct ConvexPolygon{N,T<:Real} <: Surface{T}
local_points::SVector{N,SVector{2,T}}

or maybe even better

struct ConvexPolygon{Npoints,Dim,T<:Real} <: Surface{T}
local_points::SVector{Npoints,SVector{Dim,T}}

where Dim is the dimension of the points, and Npoints is the number of points.

Copy link
Contributor

Choose a reason for hiding this comment

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

@rgal does my suggestion for defined SVector{SVector} work for this case?

Copy link
Collaborator Author

@galran galran Jul 5, 2021

Choose a reason for hiding this comment

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

I'm getting the following error when trying it:
ERROR: LoadError: TypeError: in ConvexPolygon, in N, expected N<:Int64, got a value of type Int64.
this is my entire test:

using StaticArrays

struct ConvexPolygon{N<:Int64, T<:Real} 
    local_points::SVector{N, SVector{2, T}}

    function ConvexPolygon(points::SVector{N, SVector{2, T}}) where {N<:Int64, T<:Real}
        new{N, T}(points)
    end

end

points = SVector(SVector(1.0, 1.0), SVector(2.0, 2.0), SVector(3.0, 3.0))
p = ConvexPolygon{3, Float64}(points)

Copy link
Contributor

@BrianGun BrianGun Jul 5, 2021

Choose a reason for hiding this comment

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

You don't want this:

struct ConvexPolygon{N<:Int64,T<:Real}

you want this:

struct ConvexPolygon{N,T<:Real}

The first statement is saying you want N to be a subtype of Int64, which is the same as saying a type equal to Int64 because there are no subtypes of Int64. But the argument you will use in the constructor is a number, not a type.

you'll want to get rid of the N<:Int64 in the ConvexPolygonn constructor as well. I made those two changes and the code worked. This should do it:

struct ConvexPolygon{N, T<:Real} 
    local_points::SVector{N, SVector{2, T}}

    function ConvexPolygon(points::SVector{N, SVector{2, T}}) where {N, T<:Real}
        new{N, T}(points)
    end

end


# for efficency
_poly2d::GeometricalPredicates.Polygon2D{GeometricalPredicates.Point2D}

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

@assert length(local_polygon_points) > 3 # need at least 3 points to define apolygon
@assert size(local_polygon_points[1])[1] == 2 # check that first point is a 2D point
Copy link
Contributor

Choose a reason for hiding this comment

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

if you use SVector instead of Vector then you can specify vector size variable declaration and can eliminate this assert.


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...)

plane = Plane(forward(local_frame), world_center, interface = interface)
new{T}(plane, local_frame, local_polygon_points, poly)
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)
Copy link
Contributor

Choose a reason for hiding this comment

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

what is this function for? If you really need two then the reason should probably be in a documentation string for the normal function.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I copied it from the Hexagon example. I think that its just a matter of interface support.
You can have a normal for the surface and a normal for each point on the surface in the case when it's not planer.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

It seems that all the shapes have the same interface - i think that normal without parameters is the one that should be removed but on the other hand what point would you give to a shape that doesn't have a natural parametrization.



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 = world2local(poly.local_frame) * p
@assert abs(local_p[3]) < 0.00001 # need to find a general epsilon - check that the point lies on the plain
Copy link
Contributor

Choose a reason for hiding this comment

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

Not sure this assert is necessary. If we have a ray/plane intersection then local_p should always be a valid point. An absolute test is subject to problems with scale, if the convex polygons become large for some reason. Could lead to erroneous failures.

in_poly = GeometricalPredicates.inpolygon(poly._poly2d, GeometricalPredicates.Point(local_p[1], local_p[2]))

if !in_poly
return EmptyInterval(T)
else
if dot(normal(poly), direction(r)) < zero(T)
return positivehalfspace(intersect)
else
return rayorigininterval(intersect)
end
end
end
end


function makemesh(poly::ConvexPolygon{T}, ::Int = 0) where {T<:Real}
Copy link
Contributor

Choose a reason for hiding this comment

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

is makemesh used for drawing the convex poly in 3D? Maybe a documentation string as to the purpose of this function.

c = centroid(poly)

l2w = local2world(poly.local_frame)
l = length(poly.local_points)

triangles = []
for i in 1:l
p1 = poly.local_points[i]
p2 = poly.local_points[mod(i,l) + 1]
Copy link
Contributor

Choose a reason for hiding this comment

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

might want to use a variable name other than l here since it looks almost exactly the same as 1 in the code font. Maybe len instead of l?


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
10 changes: 10 additions & 0 deletions src/Geometry/Transform.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
# -----------------------------------------------------------------------------------------------
# GEOMETRY
# -----------------------------------------------------------------------------------------------
using ...OpticSim
using StaticArrays
using LinearAlgebra

Expand Down Expand Up @@ -233,6 +234,15 @@ 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}) = normalize(Vec3(t[1,1], t[2,1], t[3,1]))
Copy link
Contributor

Choose a reason for hiding this comment

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

documentation strings would be helpful here to explain what these functions do and what they are used for. Graphics people will know but optics people might not.

up(t::Transform{<:Real}) = normalize(Vec3(t[1,2], t[2,2], t[3,2]))
forward(t::Transform{<:Real}) = normalize(Vec3(t[1,3], t[2,3], t[3,3]))
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

Expand Down
6 changes: 0 additions & 6 deletions src/Optical/Emitters/Emitters.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

"""
Expand Down
15 changes: 12 additions & 3 deletions src/Optical/Paraxial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,13 +13,14 @@ 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}
function ParaxialLens(shape::Union{Rectangle{T},Ellipse{T},Hexagon{T},ConvexPolygon{T}}, interface::ParaxialInterface{T}) where {T<:Real}
new{T}(shape, interface)
end
end
Expand All @@ -46,6 +47,14 @@ 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{Vector{T}}, local_center_point::Vector{T}; outsidematerial::OpticSim.GlassCat.AbstractGlass = OpticSim.GlassCat.Air) where {T<:Real}
@assert size(local_center_point)[1] == 2 # make sure we got a 2D point

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)
Expand All @@ -57,7 +66,7 @@ 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)
Expand Down
4 changes: 4 additions & 0 deletions src/utilities.jl
Original file line number Diff line number Diff line change
Expand Up @@ -71,3 +71,7 @@ end
return acos(x)
end
end


# some place holders for package level function names
function origin end
Copy link
Contributor

Choose a reason for hiding this comment

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

documentation string would be useful otherwise few people will understand what this function definition is for. Maybe show an example of overloading this function?

71 changes: 71 additions & 0 deletions test/testsets/Intersection.jl
Original file line number Diff line number Diff line change
Expand Up @@ -392,6 +392,77 @@
@test surfaceintersection(hex2, raymiss) isa EmptyInterval
end # testset Hexagon


@testset "Convex Polygon" begin

t = identitytransform()

local_points_2d = [
[-1.0, -1.0],
[1.0, -1.0],
[2.0, 0.0],
[1.0, 1.0],
[-1.0, 1.0],
[-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))
Expand Down