From 28bf6948685de4b11dd497e48286eae80a56429b Mon Sep 17 00:00:00 2001 From: Ran Gal <79867742+galran@users.noreply.github.com> Date: Fri, 2 Apr 2021 15:08:10 -0700 Subject: [PATCH] Galran/new geometry module (#70) * added the new Geometry module and split the file into the core geometry and geometry functions override (mainly multiplication by a transform). next step is to remove the old RigidBodyTransform and replace it with the new Transform. * remove the RigidBodyransform switched to the new Tansform. package-wide small changee, including tests, CSG and other places. next, will try to find the test issues of unwanted allocations. * HTTP version limit created a bunch of limitation with JSServe and WGLMakie that prevented all kind of other packages versions. * small changes to Transform to negate compiler inability to infer specific sizes which resulted is un-wanted allocations. after this change all tests passed successfully. * remove the old file RigidTransform.jl and Transform2.jl moved all the Transform application functions (overiding the * operators to apply transform to different types) to their appropriate locations (Interval, Triangle, TriangleMesh, Intersection, Rays) * forgot to remove the reference to Transform2.jl - fixed * doing some cleaning up and removing old commented code. will work on the documentation next. * added some documentation and eliminated all the expansion operator uses. * changed the documentation to describe the new Transform, Vec3 and Vec4 types (replacing the RigidBodyTransform). added a new page (basic_types.md). this is my first attempt at documentation, so please check well. * fixed last places that used/referenced the old RigidBodyTRansform, including in the examples file and documentation. * move the Geometry module declaration to the Geometry.jl file (instead of Transformjl). fixed a typo in the documentation. * fix missed RigidBodyTransform * removed the origin3 and origin4 methods - we still have zero3 or zero4 and we can also use the LinearAlgebra methods, such as zeors(Vec3) ones(Vec4). improved some declarations in Transform.jl according to Brian's recomendation. Changed the Transform type to be identical to an SMatrix{4, 4, T 16} removed identityT and left only identitytransform. added the basic_types.md page to the help generation process. Co-authored-by: Alfred Wong --- Project.toml | 2 +- docs/make.jl | 3 +- docs/src/basic_types.md | 55 +++ docs/src/csg.md | 14 - docs/src/examples.md | 8 +- src/Examples/Examples.jl | 9 +- src/Geometry/BoundingBox.jl | 4 +- src/Geometry/CSG/CSG.jl | 48 +- src/Geometry/CSG/Intersection.jl | 15 + src/Geometry/CSG/Interval.jl | 27 ++ src/Geometry/CSG/RigidBodyTransform.jl | 196 -------- src/Geometry/Geometry.jl | 10 +- src/Geometry/Primitives/Curves/Bezier.jl | 2 +- src/Geometry/Primitives/NonCSG/Triangle.jl | 17 + src/Geometry/Ray.jl | 8 + src/Geometry/Transform.jl | 536 +++++++++++++++++++++ src/Optical/Eye.jl | 4 +- src/Optical/Lenses.jl | 2 +- src/Optical/OpticalRay.jl | 2 +- src/Visualization/Visualization.jl | 13 +- test/TestData/TestData.jl | 15 +- test/runtests.jl | 2 + test/testsets/General.jl | 14 +- test/testsets/Visualization.jl | 2 +- 24 files changed, 734 insertions(+), 274 deletions(-) create mode 100644 docs/src/basic_types.md delete mode 100644 src/Geometry/CSG/RigidBodyTransform.jl create mode 100644 src/Geometry/Transform.jl diff --git a/Project.toml b/Project.toml index cfb1434e5..ba67fd0a0 100644 --- a/Project.toml +++ b/Project.toml @@ -47,7 +47,7 @@ DataFrames = "0.22" Distributions = "0.24" FileIO = "1.6" ForwardDiff = "0.10" -HTTP = "0.9" +HTTP = "0.8" ImageView = "0.10" Images = "0.23" Ipopt = "0.6" diff --git a/docs/make.jl b/docs/make.jl index 7cd776349..3819485eb 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -21,7 +21,7 @@ # SOFTWARE using Documenter -using OpticSim +using OpticSim, OpticSim.Geometry makedocs( sitename = "OpticSim.jl", @@ -35,6 +35,7 @@ makedocs( "Examples" => "examples.md", # "Glasses" => cat_pages, "Geometry" => [ + "Basic Types" => "basic_types.md", "Primitives" => "primitives.md", "CSG" => "csg.md" ], diff --git a/docs/src/basic_types.md b/docs/src/basic_types.md new file mode 100644 index 000000000..b1aa1f466 --- /dev/null +++ b/docs/src/basic_types.md @@ -0,0 +1,55 @@ +# Basic Types + +The following are basic geometric types and utilities, such as vectors and transforms, that are used all acoross the package. Using these types requires, in addition to the `using OpticSim` statement to also use the OpticSim.Geoemtry module like: + +```@example +using OpticSim, OpticSim.Geometry +``` + +## Vec3 + +Representing a 3D vector. + +```@docs +Geometry.Vec3 +Geometry.unitX3 +Geometry.unitY3 +Geometry.unitZ3 +Geometry.zero3 +Geometry.one3 +``` + +## Vec4 + +Representing a 4D vector + +```@docs +Geometry.Vec4 +Geometry.unitX4 +Geometry.unitY4 +Geometry.unitZ4 +Geometry.unitW4 +Geometry.zero4 +Geometry.one4 +``` + +## Transform + +Representing a general 3D transform (4x4 matrix). Currently only used as a rigid-body transform. +Transforms are used to position surfaces within the CSG tree, position emitters in 3D, etc. + +```@docs +Geometry.Transform +Geometry.identitytransform +Geometry.rotationX +Geometry.rotationY +Geometry.rotationZ +Geometry.rotation +Geometry.rotationd +Geometry.rotate +Geometry.translation +Geometry.rotmat +Geometry.rotmatd +Geometry.rotmatbetween +``` + diff --git a/docs/src/csg.md b/docs/src/csg.md index 953ea2eb1..b8170cd7c 100644 --- a/docs/src/csg.md +++ b/docs/src/csg.md @@ -46,20 +46,6 @@ TriangularPrism Spider ``` -## Transforms - -Transforms are used to position does within the CSG tree. - -```@docs -RigidBodyTransform -OpticSim.rotation -OpticSim.rotationd -OpticSim.translation -OpticSim.rotmat -OpticSim.rotmatd -OpticSim.rotmatbetween -``` - ## CSG Types These are the types of the primary CSG elements, i.e. the nodes in the CSG tree. diff --git a/docs/src/examples.md b/docs/src/examples.md index c37d28982..2e35965d6 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -72,14 +72,14 @@ nothing # hide ## Schmidt Cassegrain Telescope ```@example -using OpticSim +using OpticSim, OpticSim.Geometry using StaticArrays # glass entrance lens on telescope topsurf = Plane(SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 0.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air), vishalfsizeu = 12.00075, vishalfsizev = 12.00075) botsurf = AcceleratedParametricSurface(ZernikeSurface(12.00075, radius = -1.14659768e+4, aspherics = [(4, 3.68090959e-7), (6, 2.73643352e-11), (8, 3.20036892e-14)]), 17, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)) -coverlens = csgintersection(leaf(Cylinder(12.00075, 1.4)), csgintersection(leaf(topsurf), leaf(botsurf, RigidBodyTransform(OpticSim.rotmatd(0, 180, 0), SVector(0.0, 0.0, -0.65))))) +coverlens = csgintersection(leaf(Cylinder(12.00075, 1.4)), csgintersection(leaf(topsurf), leaf(botsurf, Transform(OpticSim.rotmatd(0, 180, 0), Vec3(0.0, 0.0, -0.65))))) # big mirror with a hole in it bigmirror = ConicLens(OpticSim.GlassCat.SCHOTT.N_BK7, -72.65, -95.2773500000134, 0.077235, Inf, 0.0, 0.2, 12.18263, frontsurfacereflectance = 1.0) bigmirror = csgdifference(bigmirror, leaf(Cylinder(4.0, 0.3, interface = opaqueinterface()), translation(0.0, 0.0, -72.75))) @@ -102,13 +102,13 @@ nothing # hide ## Lens Construction ```@example -using OpticSim +using OpticSim, OpticSim.Geometry using StaticArrays topsurface = leaf(AcceleratedParametricSurface(QTypeSurface(9.0, radius = -25.0, conic = 0.3, αcoeffs = [(1, 0, 0.3), (1, 1, 1.0)], βcoeffs = [(1, 0, -0.1), (2, 0, 0.4), (3, 0, -0.6)], normradius = 9.5), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) -lens = csgintersection(barrel, csgintersection(topsurface, botsurface))(RigidBodyTransform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) +lens = csgintersection(barrel, csgintersection(topsurface, botsurface))(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) sys = CSGOpticalSystem(LensAssembly(lens), Rectangle(15.0, 15.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) Vis.drawtracerays(sys, test = true, trackallrays = true) Vis.save("assets/qtype.png") # hide diff --git a/src/Examples/Examples.jl b/src/Examples/Examples.jl index 462f42bb5..eef970429 100644 --- a/src/Examples/Examples.jl +++ b/src/Examples/Examples.jl @@ -24,6 +24,7 @@ module Examples using ..OpticSim using ..OpticSim.Vis +using ..OpticSim.Geometry # using ..OpticSim.GlassCat use this if you want to type SCHOTT.N_BK7 rather than OpticSim.GlassCat.SCHOTT.N_BK7 using StaticArrays using DataFrames @@ -186,7 +187,7 @@ function SchmidtCassegrainTelescope() # glass entrance lens on telescope topsurf = Plane(SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 0.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air), vishalfsizeu = 12.00075, vishalfsizev = 12.00075) botsurf = AcceleratedParametricSurface(ZernikeSurface(12.00075, radius = -1.14659768e+4, aspherics = [(4, 3.68090959e-7), (6, 2.73643352e-11), (8, 3.20036892e-14)]), 17, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)) - coverlens = csgintersection(leaf(Cylinder(12.00075, 1.4)), csgintersection(leaf(topsurf), leaf(botsurf, RigidBodyTransform(OpticSim.rotmatd(0, 180, 0), SVector(0.0, 0.0, -0.65))))) + coverlens = csgintersection(leaf(Cylinder(12.00075, 1.4)), csgintersection(leaf(topsurf), leaf(botsurf, Transform(OpticSim.rotmatd(0, 180, 0), Vec3(0.0, 0.0, -0.65))))) # big mirror with a hole in it bigmirror = ConicLens(OpticSim.GlassCat.SCHOTT.N_BK7, -72.65, -95.2773500000134, 0.077235, Inf, 0.0, 0.2, 12.18263, frontsurfacereflectance = 1.0) bigmirror = csgdifference(bigmirror, leaf(Cylinder(4.0, 0.3, interface = opaqueinterface()), translation(0.0, 0.0, -72.75))) @@ -349,7 +350,7 @@ function eyetrackHOE(nrays = 5000, det = false, showhead = true, zeroorder = fal barreltop = Plane(camdir_norm, camloc) barrelbot = Plane(-camdir_norm, camloc - 3 * barrellength * camdir_norm) barrelrot = OpticSim.rotmatbetween(SVector(0.0, 0.0, 1.0), camdir_norm) - cambarrel = csgintersection(barrelbot, csgintersection(barreltop, leaf(Cylinder(camrad, barrellength, interface = opaqueinterface(Float64)), RigidBodyTransform(barrelrot, barrelloc))))() + cambarrel = csgintersection(barrelbot, csgintersection(barreltop, leaf(Cylinder(camrad, barrellength, interface = opaqueinterface(Float64)), Transform(barrelrot, barrelloc))))() camdet = Circle(sensorrad, camdir_norm, camloc - barrellength * camdir_norm, interface = opaqueinterface(Float64)) # sourceleft = hoecenter[1] + hoehalfwidth - sourceloc[1] @@ -382,8 +383,8 @@ function eyetrackHOE(nrays = 5000, det = false, showhead = true, zeroorder = fal Vis.draw!((corneavertex - 50 * d, corneavertex), color = :red) end if showhead - Vis.draw!(joinpath(@__DIR__, "../../OBJ/glasses.obj"), scale = 100.0, transform = RigidBodyTransform(OpticSim.rotmatd(90, 0, 0), [27.0, 45.0, -8.0]), color = :black) - Vis.draw!(joinpath(@__DIR__, "../../OBJ/femalehead.obj"), scale = 13.0, transform = RigidBodyTransform(OpticSim.rotmatd(0, 0, 180), [27.0, 105.0, -148.0]), color = :white) + Vis.draw!(joinpath(@__DIR__, "../../OBJ/glasses.obj"), scale = 100.0, transform = Transform(OpticSim.rotmatd(90, 0, 0), [27.0, 45.0, -8.0]), color = :black) + Vis.draw!(joinpath(@__DIR__, "../../OBJ/femalehead.obj"), scale = 13.0, transform = Transform(OpticSim.rotmatd(0, 0, 180), [27.0, 105.0, -148.0]), color = :white) end Vis.display() end diff --git a/src/Geometry/BoundingBox.jl b/src/Geometry/BoundingBox.jl index bcc469a9b..6b0b8c648 100644 --- a/src/Geometry/BoundingBox.jl +++ b/src/Geometry/BoundingBox.jl @@ -28,7 +28,7 @@ Axis-aligned three-dimensional bounding box. ```julia BoundingBox(xmin::T, xmax::T, ymin::T, ymax::T, zmin::T, zmax::T) BoundingBox(s::Surface{T}) -BoundingBox(s::ParametricSurface{T,3}, transform::RigidBodyTransform{T} = identitytransform(T)) +BoundingBox(s::ParametricSurface{T,3}, transform::Transform{T} = identitytransform(T)) BoundingBox(c::CSGTree{T}) BoundingBox(tri::Triangle{T}) BoundingBox(triangles::AbstractVector{Triangle{T}}) @@ -63,7 +63,7 @@ struct BoundingBox{T<:Real} end export BoundingBox -function BoundingBox(s::ParametricSurface{T,3}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} +function BoundingBox(s::ParametricSurface{T,3}, transform::Transform{T} = identitytransform(T)) where {T<:Real} # get the bounding box of a transformed bounding box bbox = BoundingBox(s) if transform == identitytransform(T) diff --git a/src/Geometry/CSG/CSG.jl b/src/Geometry/CSG/CSG.jl index b4b7c9243..68fa96d0f 100644 --- a/src/Geometry/CSG/CSG.jl +++ b/src/Geometry/CSG/CSG.jl @@ -87,11 +87,11 @@ As such, transforming points from the geometry using this transform puts them in """ struct LeafNode{T,S<:ParametricSurface{T,3}} <: CSGTree{T} geometry::S - transform::RigidBodyTransform{T} - invtransform::RigidBodyTransform{T} + transform::Transform{T} + invtransform::Transform{T} bbox::BoundingBox{T} - function LeafNode(a::S, transform::RigidBodyTransform{T}) where {T<:Real,S<:ParametricSurface{T,3}} + function LeafNode(a::S, transform::Transform{T}) where {T<:Real,S<:ParametricSurface{T,3}} # store the transformed bounding box so nodes higher in the tree have correct global space bounding boxes return new{T,S}(a, transform, inv(transform), BoundingBox(a, transform)) end @@ -113,7 +113,7 @@ a = Cylinder(1.0,1.0) b = Plane([0.0,0.0,1.0], [0.0,0.0,0.0]) generator = csgintersection(a,b) # now make a csg object that can be ray traced -csgobj = generator(RigidBodyTransform(1.0,1.0,2.0)) +csgobj = generator(Transform(1.0,1.0,2.0)) ``` """ struct CSGGenerator{T<:Real} @@ -121,7 +121,7 @@ struct CSGGenerator{T<:Real} end export CSGGenerator -function (a::CSGGenerator{T})(transform::RigidBodyTransform{T})::CSGTree{T} where {T<:Real} +function (a::CSGGenerator{T})(transform::Transform{T})::CSGTree{T} where {T<:Real} a.f(transform) end function (a::CSGGenerator{T})()::CSGTree{T} where {T<:Real} @@ -129,19 +129,19 @@ function (a::CSGGenerator{T})()::CSGTree{T} where {T<:Real} end """ - leaf(surf::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) -> CSGGenerator{T} + leaf(surf::ParametricSurface{T}, transform::Transform{T} = identitytransform(T)) -> CSGGenerator{T} Create a leaf node from a parametric surface with a given transform. """ -function leaf(surf::S, transform::RigidBodyTransform{T} = identitytransform(T))::CSGGenerator{T} where {T<:Real,S<:ParametricSurface{T}} +function leaf(surf::S, transform::Transform{T} = identitytransform(T))::CSGGenerator{T} where {T<:Real,S<:ParametricSurface{T}} return CSGGenerator{T}((parenttransform) -> LeafNode(surf, parenttransform * transform)) end """ - leaf(surf::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) -> CSGGenerator{T} + leaf(surf::CSGGenerator{T}, transform::Transform{T} = identitytransform(T)) -> CSGGenerator{T} Create a (pseudo) leaf node from another CSGGenerator, this is useful if you want multiple copies of a premade CSG structure with different transforms, for example in an MLA. """ -function leaf(n::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T))::CSGGenerator{T} where {T<:Real} +function leaf(n::CSGGenerator{T}, transform::Transform{T} = identitytransform(T))::CSGGenerator{T} where {T<:Real} return CSGGenerator{T}((parenttransform) -> n(parenttransform * transform)) end export leaf @@ -151,47 +151,47 @@ export leaf # Wrap the csg operations in function that delays evaluation until the transform has been computed. """ - csgintersection(a::CSGGenerator{T} b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) -> CSGGenerator{T} + csgintersection(a::CSGGenerator{T} b::CSGGenerator{T}, transform::Transform{T} = identitytransform(T)) -> CSGGenerator{T} Create a binary node in the CSG tree representing an intersection between a and b. A shortcut method for `a` and `b` as [`ParametricSurface`](@ref)s is also available. ![Intersect Image](https://upload.wikimedia.org/wikipedia/commons/0/0b/Boolean_intersect.PNG) """ -csgintersection(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = CSGGenerator{T}((parenttransform) -> IntersectionNode(a(parenttransform * transform), b(parenttransform * transform))) -csgintersection(a::ParametricSurface{T}, b::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgintersection(leaf(a), leaf(b), transform) -csgintersection(a::ParametricSurface{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgintersection(leaf(a), b, transform) -csgintersection(a::CSGGenerator{T}, b::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgintersection(a, leaf(b), transform) +csgintersection(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::Transform{T} = identitytransform(T)) where {T<:Real} = CSGGenerator{T}((parenttransform) -> IntersectionNode(a(parenttransform * transform), b(parenttransform * transform))) +csgintersection(a::ParametricSurface{T}, b::ParametricSurface{T}, transform::Transform{T} = identitytransform(T)) where {T<:Real} = csgintersection(leaf(a), leaf(b), transform) +csgintersection(a::ParametricSurface{T}, b::CSGGenerator{T}, transform::Transform{T} = identitytransform(T)) where {T<:Real} = csgintersection(leaf(a), b, transform) +csgintersection(a::CSGGenerator{T}, b::ParametricSurface{T}, transform::Transform{T} = identitytransform(T)) where {T<:Real} = csgintersection(a, leaf(b), transform) export csgintersection """ - csgunion(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) -> CSGGenerator{T} + csgunion(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::Transform{T} = identitytransform(T)) -> CSGGenerator{T} Create a binary node in the CSG tree representing a union between a and b. A shortcut method for `a` and `b` as [`ParametricSurface`](@ref)s is also available. ![Union Image](https://upload.wikimedia.org/wikipedia/commons/4/4a/Boolean_union.PNG) """ -csgunion(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = CSGGenerator{T}((parenttransform) -> UnionNode(a(parenttransform * transform), b(parenttransform * transform))) -csgunion(a::ParametricSurface{T}, b::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgunion(leaf(a), leaf(b), transform) -csgunion(a::ParametricSurface{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = ccsgunion(leaf(a), b, transform) -csgunion(a::CSGGenerator{T}, b::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgunion(a, leaf(b), transform) +csgunion(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::Transform{T} = identitytransform(T)) where {T<:Real} = CSGGenerator{T}((parenttransform) -> UnionNode(a(parenttransform * transform), b(parenttransform * transform))) +csgunion(a::ParametricSurface{T}, b::ParametricSurface{T}, transform::Transform{T} = identitytransform(T)) where {T<:Real} = csgunion(leaf(a), leaf(b), transform) +csgunion(a::ParametricSurface{T}, b::CSGGenerator{T}, transform::Transform{T} = identitytransform(T)) where {T<:Real} = ccsgunion(leaf(a), b, transform) +csgunion(a::CSGGenerator{T}, b::ParametricSurface{T}, transform::Transform{T} = identitytransform(T)) where {T<:Real} = csgunion(a, leaf(b), transform) export csgunion """ - csgdifference(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) -> CSGGenerator{T} + csgdifference(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::Transform{T} = identitytransform(T)) -> CSGGenerator{T} Create a binary node in the CSG tree representing the difference of a and b, essentially a - b. A shortcut method for `a` and `b` as [`ParametricSurface`](@ref)s is also available. ![Difference Image](https://upload.wikimedia.org/wikipedia/commons/8/86/Boolean_difference.PNG) """ -csgdifference(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = CSGGenerator{T}((parenttransform) -> IntersectionNode(a(parenttransform * transform), ComplementNode(b(parenttransform * transform)))) -csgdifference(a::CSGGenerator{T}, b::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgdifference(a, leaf(b), transform) -csgdifference(a::ParametricSurface{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgdifference(leaf(a), b, transform) -csgdifference(a::ParametricSurface{T}, b::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgdifference(leaf(a), leaf(b), transform) +csgdifference(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::Transform{T} = identitytransform(T)) where {T<:Real} = CSGGenerator{T}((parenttransform) -> IntersectionNode(a(parenttransform * transform), ComplementNode(b(parenttransform * transform)))) +csgdifference(a::CSGGenerator{T}, b::ParametricSurface{T}, transform::Transform{T} = identitytransform(T)) where {T<:Real} = csgdifference(a, leaf(b), transform) +csgdifference(a::ParametricSurface{T}, b::CSGGenerator{T}, transform::Transform{T} = identitytransform(T)) where {T<:Real} = csgdifference(leaf(a), b, transform) +csgdifference(a::ParametricSurface{T}, b::ParametricSurface{T}, transform::Transform{T} = identitytransform(T)) where {T<:Real} = csgdifference(leaf(a), leaf(b), transform) export csgdifference diff --git a/src/Geometry/CSG/Intersection.jl b/src/Geometry/CSG/Intersection.jl index 1b9a82995..83cf3e4c8 100644 --- a/src/Geometry/CSG/Intersection.jl +++ b/src/Geometry/CSG/Intersection.jl @@ -183,3 +183,18 @@ israyorigin(::RayOrigin) = true α(::RayOrigin{T}) where {T<:Real} = zero(T) Base.eltype(::RayOrigin{T}) where {T<:Real} = T + + +""" +Apply a Transform to an Intersection object +""" +function Base.:*(a::Transform{T}, intsct::Intersection{T,3})::Intersection{T,3} where {T<:Real} + u, v = uv(intsct) + i = interface(intsct) + if VERSION < v"1.6.0-DEV" + # TODO REMOVE + return @unionsplit OpticalInterface T i Intersection(α(intsct), a * point(intsct), Geometry.rotate(a, normal(intsct)), u, v, i, flippednormal = flippednormal(intsct)) + else + return Intersection(α(intsct), a * point(intsct), Geometry.rotate(a, normal(intsct)), u, v, interface(intsct), flippednormal = flippednormal(intsct)) + end +end diff --git a/src/Geometry/CSG/Interval.jl b/src/Geometry/CSG/Interval.jl index 614cbda11..0a74bcb9b 100644 --- a/src/Geometry/CSG/Interval.jl +++ b/src/Geometry/CSG/Interval.jl @@ -575,3 +575,30 @@ function intervalcomplement(a::Interval{T}) where {T<:Real} end difference(a::DisjointUnion, b::DisjointUnion) = intervalintersection(a, intervalcomplement(b)) + + +""" +Apply a Transform to an Interval object +""" +function Base.:*(transformation::Transform{T}, a::Interval{T}) where {T<:Real} + # looks ridiculous but necessary to dissambiguate the elements of the interval + u = upper(a) + l = lower(a) + if l isa RayOrigin{T} + if u isa Infinity{T} + return Interval(l, u) + else + u = transformation * u + return Interval(l, u) + end + else + l = transformation * l + if u isa Infinity{T} + return Interval(l, u) + + else + u = transformation * u + return Interval(l, u) + end + end +end diff --git a/src/Geometry/CSG/RigidBodyTransform.jl b/src/Geometry/CSG/RigidBodyTransform.jl deleted file mode 100644 index fe61c42a6..000000000 --- a/src/Geometry/CSG/RigidBodyTransform.jl +++ /dev/null @@ -1,196 +0,0 @@ -# MIT License - -# Copyright (c) Microsoft Corporation. - -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: - -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. - -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE - -""" - RigidBodyTransform{S<:Real} - -Transform encapsulating rotation and translation in 3D space. Translation happens **after** rotation. - -```julia -RigidBodyTransform{S}(θ::T, ϕ::T, ψ::T, x::T, y::T, z::T) -RigidBodyTransform(rotation::SMatrix{3,3,S}, translation::SVector{3,S}) -RigidBodyTransform(rotation::AbstractArray{S,2}, translation::AbstractArray{S,1}) -``` -`θ`, `ϕ` and `ψ` in first constructor are in **radians**. -""" -struct RigidBodyTransform{T<:Real} - rotation::SMatrix{3,3,T,9} - translation::SVector{3,T} - - RigidBodyTransform{S}(θ::T, ϕ::T, ψ::T, x::T, y::T, z::T) where {T<:Number,S<:Real} = new{S}(rotmat(S, θ, ϕ, ψ), SVector{3,S}(x, y, z)) - RigidBodyTransform(rotation::SMatrix{3,3,T}, translation::SVector{3,T}) where {T<:Real} = new{T}(rotation, translation) - function RigidBodyTransform(rotation::AbstractArray{T,2}, translation::AbstractArray{T,1}) where {T<:Real} - @assert size(rotation)[1] == size(rotation)[2] == length(translation) == 3 - return new{T}(SMatrix{3,3,T}(rotation), SVector{3,T}(translation)) - end -end -export RigidBodyTransform - -""" - identitytransform([S::Type]) -> RigidBodyTransform{S} - -Returns the [`RigidBodyTransform`](@ref) of type `S` (default `Float64`) representing the identity transform. -""" -identitytransform(::Type{T} = Float64) where {T<:Real} = RigidBodyTransform(SMatrix{3,3,T,9}(I), zeros(SVector{3,T})) -export identitytransform - -""" - rotmatbetween([S::Type], a::SVector{3,T}, b::SVector{3,T}) -> SMatrix{3,3,S} - -Returns the rotation matrix of type `S` (default `Float64`) representing the rotation between vetors `a` and `b`. -""" -rotmatbetween(a::SVector{3,T}, b::SVector{3,T}) where {T<:Number} = rotmatbetween(Float64, a, b) -function rotmatbetween(::Type{S}, a::SVector{3,T}, b::SVector{3,T}) where {T<:Number,S<:Real} - v = cross(a, b) - c = dot(a, b) - V = SMatrix{3,3,S,9}(0, v[3], -v[2], -v[3], 0, v[1], v[2], -v[1], 0) - R = I + V + V^2 * one(T) / (one(T) + c) - return SMatrix{3,3,S,9}(R) -end - -""" - rotmatd([S::Type], θ::T, ϕ::T, ψ::T) -> SMatrix{3,3,S} - -Returns the rotation matrix of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in degrees**. -""" -rotmatd(θ::T, ϕ::T, ψ::T) where {T<:Number} = rotmat(Float64, θ * π / 180, ϕ * π / 180, ψ * π / 180) -rotmatd(::Type{S}, θ::T, ϕ::T, ψ::T) where {T<:Number,S<:Real} = rotmat(S, θ * π / 180, ϕ * π / 180, ψ * π / 180) -""" - rotmat([S::Type], θ::T, ϕ::T, ψ::T) -> SMatrix{3,3,S} - -Returns the rotation matrix of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in radians**. -""" -rotmat(θ::T, ϕ::T, ψ::T) where {T<:Number} = rotmat(Float64, θ, ϕ, ψ) -function rotmat(::Type{S}, θ::T, ϕ::T, ψ::T) where {T<:Number,S<:Real} - sinψ = sin(ψ) - sinϕ = sin(ϕ) - sinθ = sin(θ) - cosψ = cos(ψ) - cosϕ = cos(ϕ) - cosθ = cos(θ) - return SMatrix{3,3,S,9}(cosψ * cosϕ, sinψ * cosϕ, -sinϕ, cosψ * sinϕ * sinθ - sinψ * cosθ, sinψ * sinϕ * sinθ + cosψ * cosθ, cosϕ * sinθ, cosψ * sinϕ * cosθ + sinψ * sinθ, sinψ * sinϕ * cosθ - cosψ * sinθ, cosϕ * cosθ) -end -""" - translation([S::Type], x::T, y::T, z::T) -> RigidBodyTransform{S} - -Returns the [`RigidBodyTransform`](@ref) of type `S` (default `Float64`) representing the translation by `x`, `y` and `z`. -""" -translation(x::T, y::T, z::T) where {T<:Number} = translation(Float64, x, y, z) -translation(::Type{S}, x::T, y::T, z::T) where {T<:Number,S<:Real} = RigidBodyTransform(SMatrix{3,3,S,9}(I), SVector{3,S}(x, y, z)) -""" - rotation([S::Type], θ::T, ϕ::T, ψ::T) -> RigidBodyTransform{S} - -Returns the [`RigidBodyTransform`](@ref) of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in radians**. -""" -rotation(θ::T, ϕ::T, ψ::T) where {T<:Number} = rotation(Float64, θ, ϕ, ψ) -rotation(::Type{S}, θ::T, ϕ::T, ψ::T) where {T<:Number,S<:Real} = RigidBodyTransform(rotmat(S, θ, ϕ, ψ), zeros(SVector{3,S})) -""" - rotationd([S::Type], θ::T, ϕ::T, ψ::T) -> RigidBodyTransform{S} - -Returns the [`RigidBodyTransform`](@ref) of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in degrees**. -""" -rotationd(θ::T, ϕ::T, ψ::T) where {T<:Number} = rotationd(Float64, θ, ϕ, ψ) -rotationd(::Type{S}, θ::T, ϕ::T, ψ::T) where {T<:Number,S<:Real} = RigidBodyTransform(rotmatd(S, θ, ϕ, ψ), zeros(SVector{3,S})) -export translation, rotation, rotationd - -function Base.:*(a::RigidBodyTransform{T}, r::Ray{T,3})::Ray{T,3} where {T} - return Ray(a * origin(r), rotate(a, direction(r))) -end - -function Base.:*(a::RigidBodyTransform{T}, intsct::Intersection{T,3})::Intersection{T,3} where {T<:Real} - u, v = uv(intsct) - i = interface(intsct) - if VERSION < v"1.6.0-DEV" - # TODO REMOVE - return @unionsplit OpticalInterface T i Intersection(α(intsct), a * point(intsct), rotate(a, normal(intsct)), u, v, i, flippednormal = flippednormal(intsct)) - else - return Intersection(α(intsct), a * point(intsct), rotate(a, normal(intsct)), u, v, interface(intsct), flippednormal = flippednormal(intsct)) - end -end - -function Base.:*(transformation::RigidBodyTransform{T}, a::Interval{T}) where {T<:Real} - # looks ridiculous but necessary to dissambiguate the elements of the interval - u = upper(a) - l = lower(a) - if l isa RayOrigin{T} - if u isa Infinity{T} - return Interval(l, u) - else - u = transformation * u - return Interval(l, u) - end - else - l = transformation * l - if u isa Infinity{T} - return Interval(l, u) - else - u = transformation * u - return Interval(l, u) - end - end -end - -function Base.:*(a::RigidBodyTransform{T}, tmesh::TriangleMesh{T})::TriangleMesh{T} where {T<:Real} - newT = Vector{Triangle{T}}(undef, length(tmesh.triangles)) - @inbounds @simd for i in 1:length(tmesh.triangles) - newT[i] = a * tmesh.triangles[i] - end - return TriangleMesh(newT) -end - -Base.:*(a::RigidBodyTransform{T}, t::Triangle{T}) where {T<:Real} = Triangle(a * vertex(t, 1), a * vertex(t, 2), a * vertex(t, 3)) - -function Base.:*(a::RigidBodyTransform{T}, vector::SVector{3,T}) where {T<:Real} - if a.rotation === SMatrix{3,3,T,9}(I) - return vector + a.translation - else - # in julia 0 * Inf = NaN, so if the vector has any infinite value in it then the rotation breaks - # have to manually test and correct this - this can still fail when summing Infs of different signs (e.g. Inf + -Inf = NaN) - # this happens when the rotation isn't axis aligned so the result becomes less meaningful anyway - # in this case we fall back to infinite bounding boxes anyway which are usually clipped later in the csg - # TODO better way to do this? - if any(isinf.(vector)) - x = zero(T) - y = zero(T) - z = zero(T) - @inbounds begin - x += abs(a.rotation[1, 1]) <= eps(T) ? zero(T) : a.rotation[1, 1] * vector[1] - x += abs(a.rotation[1, 2]) <= eps(T) ? zero(T) : a.rotation[1, 2] * vector[2] - x += abs(a.rotation[1, 3]) <= eps(T) ? zero(T) : a.rotation[1, 3] * vector[3] - y += abs(a.rotation[2, 1]) <= eps(T) ? zero(T) : a.rotation[2, 1] * vector[1] - y += abs(a.rotation[2, 2]) <= eps(T) ? zero(T) : a.rotation[2, 2] * vector[2] - y += abs(a.rotation[2, 3]) <= eps(T) ? zero(T) : a.rotation[2, 3] * vector[3] - z += abs(a.rotation[3, 1]) <= eps(T) ? zero(T) : a.rotation[3, 1] * vector[1] - z += abs(a.rotation[3, 2]) <= eps(T) ? zero(T) : a.rotation[3, 2] * vector[2] - z += abs(a.rotation[3, 3]) <= eps(T) ? zero(T) : a.rotation[3, 3] * vector[3] - end - return SVector(x, y, z) + a.translation - else - return a.rotation * vector + a.translation - end - end -end -Base.:*(a::RigidBodyTransform{T}, b::RigidBodyTransform{T}) where {T<:Real} = RigidBodyTransform(a.rotation * b.rotation, a.translation + a.rotation * b.translation) -rotate(a::RigidBodyTransform{T}, vector::SVector{3,T}) where {T<:Real} = a.rotation * vector -Base.inv(a::RigidBodyTransform{T}) where {T<:Real} = RigidBodyTransform(a.rotation', -a.rotation' * a.translation) -Base.print(io::IO, a::RigidBodyTransform) = print(io, hcat(a.rotation, a.translation)) -Base.collect(a::RigidBodyTransform) = hcat(a.rotation, a.translation) diff --git a/src/Geometry/Geometry.jl b/src/Geometry/Geometry.jl index 3a72ed5b7..1d7ac5f53 100644 --- a/src/Geometry/Geometry.jl +++ b/src/Geometry/Geometry.jl @@ -20,6 +20,14 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE +# declaring the Geometry module. in the future we might add or move existing code to it. currently it contains the basic types: Vec3, Vec4 and Transform +module Geometry + +include("Transform.jl") + +end # module Geometry +using .Geometry + include("Ray.jl") include("Surface.jl") @@ -32,8 +40,6 @@ include("CSG/Interval.jl") include("Primitives/NonCSG/Triangle.jl") -include("CSG/RigidBodyTransform.jl") - include("BoundingBox.jl") include("Primitives/Plane.jl") diff --git a/src/Geometry/Primitives/Curves/Bezier.jl b/src/Geometry/Primitives/Curves/Bezier.jl index b4e34a673..da3e01d70 100644 --- a/src/Geometry/Primitives/Curves/Bezier.jl +++ b/src/Geometry/Primitives/Curves/Bezier.jl @@ -108,7 +108,7 @@ export BezierSurface # """this function transforms the control points of the Bezier patch before ray tracing. Probably a better idea to use transforms in CSG operators, # although this could be a little less computation at run time.""" -# function Base.:*(a::RigidBodyTransform{T}, surf::BezierSurface{P,T,N,M}) where {P,T,N,M} +# function Base.:*(a::Transform{T}, surf::BezierSurface{P,T,N,M}) where {P,T,N,M} # result = deepcopy(surf) # for index in CartesianIndices(surf.controlpolygon) # result.controlpolygon[index] = a * result.controlpolygon[index] diff --git a/src/Geometry/Primitives/NonCSG/Triangle.jl b/src/Geometry/Primitives/NonCSG/Triangle.jl index dc33111b2..764118b2b 100644 --- a/src/Geometry/Primitives/NonCSG/Triangle.jl +++ b/src/Geometry/Primitives/NonCSG/Triangle.jl @@ -160,3 +160,20 @@ function makiemesh(tmesh::TriangleMesh{T}) where {T<:Real} end return (points, indices) end + + +""" +Apply a Transform to a TriangleMesh object +""" +function Base.:*(a::Transform{T}, tmesh::TriangleMesh{T})::TriangleMesh{T} where {T<:Real} + newT = Vector{Triangle{T}}(undef, length(tmesh.triangles)) + @inbounds @simd for i in 1:length(tmesh.triangles) + newT[i] = a * tmesh.triangles[i] + end + return TriangleMesh(newT) +end + +""" +Apply a Transform to a Triangle object +""" +Base.:*(a::Transform{T}, t::Triangle{T}) where {T<:Real} = Triangle(a * vertex(t, 1), a * vertex(t, 2), a * vertex(t, 3)) diff --git a/src/Geometry/Ray.jl b/src/Geometry/Ray.jl index 552c08686..9dc22f8fc 100644 --- a/src/Geometry/Ray.jl +++ b/src/Geometry/Ray.jl @@ -120,3 +120,11 @@ end Computes the alpha corresponding to the closest position on the ray to point """ α(ray::AbstractRay{T,N}, point::SVector{N,T}) where {T<:Real,N} = dot(direction(ray), (point .- origin(ray))) + + +""" +Apply a Transform to a Ray object +""" +function Base.:*(a::Transform{T}, r::Ray{T,3})::Ray{T,3} where {T} + return Ray(a * origin(r), Geometry.rotate(a, direction(r))) +end diff --git a/src/Geometry/Transform.jl b/src/Geometry/Transform.jl new file mode 100644 index 000000000..7f3902559 --- /dev/null +++ b/src/Geometry/Transform.jl @@ -0,0 +1,536 @@ +# ----------------------------------------------------------------------------------------------- +# GEOMETRY +# ----------------------------------------------------------------------------------------------- +using StaticArrays +using LinearAlgebra + +#region Vec3 + +""" +`Vec3{T}` provides an immutable vector of fixed length 3 and type `T`. + +`Vec3` defines a series of convenience constructors, so you can just type e.g. `Vec3(1, 2, 3)` or `Vec3([1.0, 2.0, 3.0])`. +It also supports comprehensions, and the `zeros()`, `ones()`, `fill()`, `rand()` and `randn()` functions, such as `Vec3(rand(3))`. +""" +struct Vec3{T} <: FieldVector{3, T} + _x::T + _y::T + _z::T +end +export Vec3 + +# empty constructor - initialized with zeros +Vec3(::Type{T} = Float64) where {T<:Real} = zeros(Vec3{T}) + +""" +returns the unit vector `[1, 0, 0]` +""" +unitX3(::Type{T} = Float64) where {T<:Real} = Vec3{T}(one(T), zero(T), zero(T)) +""" +returns the unit vector `[0, 1, 0]` +""" +unitY3(::Type{T} = Float64) where {T<:Real} = Vec3{T}(zero(T), one(T), zero(T)) +""" +returns the unit vector `[0, 0, 1]` +""" +unitZ3(::Type{T} = Float64) where {T<:Real} = Vec3{T}(zero(T), zero(T), one(T)) +""" +returns the unit vector `[0, 0, 0]` +""" +zero3(::Type{T} = Float64) where {T<:Real} = zeros(Vec3{T}) +""" +returns the unit vector `[1, 1, 1]` +""" +one3(::Type{T} = Float64) where {T<:Real} = Vec3{T}(one(T), one(T), one(T)) + +export unitX3, unitY3, unitZ3, zero3, one3 + +#endregion Vec3 + +#region Vec4 + +""" +`Vec4{T}` provides an immutable vector of fixed length 4 and type `T`. + +`Vec4` defines a series of convenience constructors, so you can just type e.g. `Vec3(1, 2, 3, 4)` or `Vec3([1.0, 2.0, 3.0, 4.0])`. +It also supports comprehensions, and the `zeros()`, `ones()`, `fill()`, `rand()` and `randn()` functions, such as `Vec4(rand(4))`. +""" +struct Vec4{T} <: FieldVector{4, T} + _x::T + _y::T + _z::T + _w::T +end +export Vec4 + +# empty constructor - initialized with zeros +Vec4(::Type{T} = Float64) where {T<:Real} = zeros(Vec4{T}) + +# convert vec3 to vec4 +""" + Vec4(v::Vec3{T}) -> Vec4 + +Accept `Vec3` and create a `Vec4` type [v[1], v[2], v[3], 1] +""" +function Vec4(v::Vec3{T}) where {T<:Real} + return Vec4{T}(v[1], v[2], v[3], one(T)) +end + +""" + Vec4(v::SVector{3, T}) where {T<:Real} -> Vec4{T} + +Accept `SVector` and create a `Vec4` type [v[1], v[2], v[3], 1] +""" +function Vec4(v::SVector{3, T}) where {T<:Real} + return Vec4{T}(v[1], v[2], v[3], one(T)) +end + +""" +returns the unit vector `[1, 0, 0, 0]` +""" +unitX4(::Type{T} = Float64) where {T<:Real} = Vec4{T}(one(T), zero(T), zero(T), zero(T)) +""" +returns the unit vector `[0, 1, 0, 0]` +""" +unitY4(::Type{T} = Float64) where {T<:Real} = Vec4{T}(zero(T), one(T), zero(T), zero(T)) +""" +returns the unit vector `[0, 0, 1, 0]` +""" +unitZ4(::Type{T} = Float64) where {T<:Real} = Vec4{T}(zero(T), zero(T), one(T), zero(T)) +""" +returns the unit vector `[0, 0, 0, 1]` +""" +unitW4(::Type{T} = Float64) where {T<:Real} = Vec4{T}(zero(T), zero(T), zero(T), one(T)) +""" +returns the unit vector `[0, 0, 0, 0]` +""" +zero4(::Type{T} = Float64) where {T<:Real} = zeros(Vec4{T}) +""" +returns the unit vector `[1, 1, 1, 1]` +""" +one4(::Type{T} = Float64) where {T<:Real} = Vec4{T}(one(T), one(T), one(T), one(T)) + +export unitX4, unitY4, unitZ4, unitW4, zero4, one4 +#endregion Vec4 + + + +#region Transform + +# utility function that given a vector, return 2 orthogonal vectors to that one +function get_orthogonal_vectors(direction::Vec3{T}) where {T<:Real} + axis1 = normalize(direction) + + dp = dot(unitX3(T), axis1) + # check the special case where the input vector is parallel to the X axis + if ( dp >= one(T) - eps(T)) + # axis1 = unitX(T); + axis2 = unitY3(T); + axis3 = unitZ3(T); + return (axis2, axis3) + elseif ( dp <= -one(T) + eps(T)) + # axis1 = -unitX(T); + axis2 = -unitY3(T); + axis3 = -unitZ3(T); + return (axis2, axis3) + end + + axis3 = normalize(cross(axis1, unitX3(T))) + axis2 = normalize(cross(axis3, axis1)) + return (axis2, axis3) +end + +#--------------------------------------- +# 3D Transform / Local Frame +#--------------------------------------- +""" + Transform{S<:Real} + +Transform encapsulating rotation, translation and scale in 3D space. Translation happens **after** rotation. + +```julia +Transform{S}(θ::T, ϕ::T, ψ::T, x::T, y::T, z::T) +Transform(rotation::SMatrix{3,3,S}, translation::SVector{3,S}) +Transform(rotation::AbstractArray{S,2}, translation::AbstractArray{S,1}) +``` +`θ`, `ϕ` and `ψ` in first constructor are in **radians**. +""" +Transform{T} = SMatrix{4,4,T,16} +export Transform + + +# for compatability ith the "old" RigidBodyTransform +""" +identitytransform([S::Type]) -> Transform{S} + +Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the identity transform. +""" +identitytransform(::Type{T} = Float64) where {T<:Real} = Transform{T}( + one(T), zero(T), zero(T), zero(T), + zero(T), one(T), zero(T), zero(T), + zero(T), zero(T), one(T), zero(T), + zero(T), zero(T), zero(T), one(T) +) +export identitytransform + + +""" + Transform([S::Type]) -> Transform{S} + +Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the identity transform. +""" +function Transform(::Type{T} = Float64) where {T<:Real} + return identitytransform(T) +end + +""" + Transform(colx::Vec3{T}, coly::Vec3{T},colz::Vec3{T}, colw::Vec3{T}, ::Type{T} = Float64) where {T<:Real} + +Costruct a transform from the input columns. +""" +function Transform(colx::Vec3{T}, coly::Vec3{T}, colz::Vec3{T}, colw::Vec3{T} = zero3(T)) where {T<:Real} + return Transform{T}( + colx[1], colx[2], colx[3], zero(T), + coly[1], coly[2], coly[3], zero(T), + colz[1], colz[2], colz[3], zero(T), + colw[1], colw[2], colw[3], one(T) + ) +end + +""" + Transform(colx::Vec3{T}, coly::Vec3{T},colz::Vec3{T}, colw::Vec3{T}, ::Type{T} = Float64) where {T<:Real} + +Costruct a transform from the input columns. +""" +function Transform(colx::Vec4{T}, coly::Vec4{T}, colz::Vec4{T}, colw::Vec4{T}) where {T<:Real} + return Transform{T}( + colx[1], colx[2], colx[3], colx[4], + coly[1], coly[2], coly[3], coly[4], + colz[1], colz[2], colz[3], colz[4], + colw[1], colw[2], colw[3], colw[4] + ) +end + +""" + Transform(origin, forward) -> Transform{S} + +Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the local frame with origin and forward direction. the other 2 axes are computed automaticlly. +""" +function Transform(origin::Vec3{T}, forward::Vec3{T} = unitZ3(); type::Type{T} = Float64) where {T<:Real} + forward = normalize(forward) + right, up = get_orthogonal_vectors(forward) + return Transform(right, up, forward, origin) +end + +function Transform{S}(θ::T, ϕ::T, ψ::T, x::T, y::T, z::T; type::Type{S} = Float64) where {T<:Number,S<:Real} + temp_transform = Transform(rotmat(S, θ, ϕ, ψ), Vec3{S}(x, y, z)) + return Transform{S}(temp_transform) +end + +""" + Transform(rotation::SMatrix{3,3,T}, translation::SVector{3,T}) where {T<:Real} -> Transform{S} + +Returns the [`Transform`](@ref) of type `S` (default `Float64`) created by a rotation matrix and translation vector. +""" +function Transform(rotation::SMatrix{3,3,T}, translation::SVector{3,T}) where {T<:Real} + return Transform( + rotation[1,1], rotation[2,1], rotation[3,1], zero(T), + rotation[1,2], rotation[2,2], rotation[3,2], zero(T), + rotation[1,3], rotation[2,3], rotation[3,3], zero(T), + translation[1], translation[2], translation[3], one(T)) +end + +""" + Transform(rotation::AbstractArray{T,2}, translation::AbstractArray{T,1}) where {T<:Real} -> Transform{S} + +Returns the [`Transform`](@ref) of type `S` (default `Float64`) created by a rotation matrix (3x3) and translation vector of length 3. +""" +function Transform(rotation::AbstractArray{T,2}, translation::AbstractArray{T,1}) where {T<:Real} + @assert size(rotation)[1] == size(rotation)[2] == length(translation) == 3 + return Transform( + rotation[1,1], rotation[2,1], rotation[3,1], zero(T), + rotation[1,2], rotation[2,2], rotation[3,2], zero(T), + rotation[1,3], rotation[2,3], rotation[3,3], zero(T), + translation[1], translation[2], translation[3], one(T)) +end + + +""" + rotationX(angle::T) where {T<:Real} -> Transform + +Builds a rotation matrix for a rotation around the x-axis. +Parameters: + The counter-clockwise `angle` in radians. +""" +function rotationX(angle::T) where {T<:Real} + c = cos(angle); + s = sin(angle); + + row1 = unitX4() + row2 = Vec4(zero(T), c, s, zero(T)) + row3 = Vec4(zero(T), -s, c, zero(T)) + row4 = unitW4() + + # transposing because the constructors treat these vectors as columns instead of rows + return transpose(Transform(row1, row2, row3, row4)) +end +export rotationX + +""" + rotationY(angle::T) where {T<:Real} -> Transform + +Builds a rotation matrix for a rotation around the y-axis. +Parameters: + The counter-clockwise `angle` in radians. +""" +function rotationY(angle::T) where {T<:Real} + c = cos(angle); + s = sin(angle); + + row1 = Vec4(c, zero(T), -s, zero(T)) + row2 = unitY4() + row3 = Vec4(s, zero(T), c, zero(T)) + row4 = unitW4() + + # transposing because the constructors treat these vectors as columns instead of rows + return transpose(Transform(row1, row2, row3, row4)) +end +export rotationY + +""" + rotationZ(angle::T) where {T<:Real} -> Transform + +Builds a rotation matrix for a rotation around the z-axis. +Parameters: + The counter-clockwise `angle` in radians. +""" +function rotationZ(angle::T) where {T<:Real} + c = cos(angle); + s = sin(angle); + + row1 = Vec4(c, s, zero(T), zero(T)) + row2 = Vec4(-s, c, zero(T), zero(T)) + row3 = unitZ4() + row4 = unitW4() + + # transposing because the constructors treat these vectors as columns instead of rows + return transpose(Transform(row1, row2, row3, row4)) +end +export rotationZ + +""" + rotation(t::Transform{T}) where {T<:Real} -> SMatrix{3,3,T} + +returns the rotation part of the transform `t` - a 3x3 matrix. +""" +function rotation(t::Transform{T}) where {T<:Real} + rot = SMatrix{3, 3, T}( + t[1, 1], t[2, 1], t[3, 1], + t[1, 2], t[2, 2], t[3, 2], + t[1, 3], t[2, 3], t[3, 3]) + return Transform(rot, zero3(T)) +end + +""" + rotate(a::Transform{T}, vector::Union{Vec3{T}, SVector{3,T}}) where {T<:Real} -> Vec3{T} + +apply the rotation part of the transform `a` to the vector `vector` - this operation is usually used to rotate direction vectors. +""" +rotate(a::Transform{T}, vector::Union{Vec3{T}, SVector{3,T}}) where {T<:Real} = rotation(a) * vector + + +""" + rotation([S::Type], θ::T, ϕ::T, ψ::T) -> Transform{S} + +Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in radians**. +""" +rotation(θ::T, ϕ::T, ψ::T) where {T<:Number} = rotation(Float64, θ, ϕ, ψ) +rotation(::Type{S}, θ::T, ϕ::T, ψ::T) where {T<:Number,S<:Real} = Transform(rotmat(S, θ, ϕ, ψ), zeros(SVector{3,S})) +""" + rotationd([S::Type], θ::T, ϕ::T, ψ::T) -> Transform{S} + +Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in degrees**. +""" +rotationd(θ::T, ϕ::T, ψ::T) where {T<:Number} = rotationd(Float64, θ, ϕ, ψ) +rotationd(::Type{S}, θ::T, ϕ::T, ψ::T) where {T<:Number,S<:Real} = Transform(rotmatd(S, θ, ϕ, ψ), zeros(SVector{3,S})) +export translation, rotation, rotationd +""" + translation(x::T, y::T, z::T) where {T<:Real} + +Creates a translation transform +""" +translation(::Type{S}, x::T, y::T, z::T) where {T<:Number,S<:Real} = convert(Transform{S},translation(x, y, z)) +function translation(x::T, y::T, z::T) where {T<:Real} + col1 = unitX4(T) + col2 = unitY4(T) + col3 = unitZ4(T) + col4 = Vec4(x, y, z, one(T)) + + return Transform(col1, col2, col3, col4) +end + +""" + translation(x::T, y::T, z::T) where {T<:Real} + +Creates a translation transform +""" +function translation(t::Vec3{T}) where {T<:Real} + return translation(t[1], t[2], t[3]) +end +export translation + +""" + scale(x::T, y::T, z::T) where {T<:Real} + +Creates a scaling transform +""" +function scale(x::T, y::T, z::T) where {T<:Real} + col1 = unitX4(T) * x + col2 = unitY4(T) * y + col3 = unitZ4(T) * z + col4 = unitW4(T) + + return Transform(col1, col2, col3, col4) +end + +""" + scale(s::T) where {T<:Real} + +Creates a uniform scaling transform +""" +function scale(s::T) where {T<:Real} + return scale(s, s, s) +end + +""" + scale(t::Vec3{T}) where {T<:Real} + +Creates a scaling transform +""" +function scale(t::Vec3{T}) where {T<:Real} + return scale(t[1], [2], [3]) +end +export scale + +""" + local2world(t::Transform{T}) where {T<:Real} + +return the transform matrix that takes a point in the local coordinate system to the global one +""" +function local2world(t::Transform{T}) where {T<:Real} + return t +end +export local2world + +""" + world2local(t::Transform{T}) where {T<:Real} + +return the transform matrix that takes a point in the global coordinate system to the local one +""" +function world2local(t::Transform{T}) where {T<:Real} + return inv(t) +end +export world2local + +function Base.:*(t::Transform{T}, v::Vec3{T}) where {T<:Real} + res = t * Vec4(v) + if (t[4,4] == one(T)) + return Vec3(res[1], res[2], res[3]) + else + return Vec3(res[1]/res[4], res[2]/res[4], res[3]/res[4]) + end +end + +function Base.:*(t::Transform{T}, v::SVector{3,T}) where {T<:Real} + res = t * Vec4(v) + if (t[4,4] == one(T)) + return SVector(res[1], res[2], res[3]) + else + return SVector(res[1]/res[4], res[2]/res[4], res[3]/res[4]) + end +end + +# function Base.:*(t::Transform{T}, v::SVector{3,T}) where {T<:Real} +# res = t * Vec4(v) +# if (t[4,4] == one(T)) +# return Vec3(res[1], res[2], res[3]) +# else +# return Vec3(res[1]/res[4], res[2]/res[4], res[3]/res[4]) +# end +# end + +# function Base.:*(t::Transform{T}, v::Vec4{T}) where {T<:Real} +# res = SMatrix(t) * v +# end + + +""" + decomposeRTS(tr::Transform{T}) where {T<:Real} + +return a touple containing the rotation matrix, the translation vector and the scale vecto represnting the transform. +""" +function decomposeRTS(tr::Transform{T}) where {T<:Real} + t = Vec3(tr[1,4], tr[2,4], tr[3,4]) + sx = norm(Vec3(tr[1,1], tr[2,1], tr[3,1])) + sy = norm(Vec3(tr[1,2], tr[2,2], tr[3,2])) + sz = norm(Vec3(tr[1,3], tr[2,3], tr[3,3])) + s = Vec3(sx, sy, sz) + rot = SMatrix{4, 4, T}(tr[1,1]/sx, tr[2, 1]/sx, tr[3,1]/sx, 0, tr[1,2]/sy, tr[2, 2]/sy, tr[3,2]/sy, 0, tr[1,3]/sz, tr[2, 3]/sz, tr[3,3]/sz, 0, 0, 0, 0, 1) + + return rot, t, s +end +export decomposeRTS + + +""" + rotmatbetween([S::Type], a::SVector{3,T}, b::SVector{3,T}) -> SMatrix{3,3,S} + +Returns the rotation matrix of type `S` (default `Float64`) representing the rotation between vetors `a` and `b`, i.e. rotation(a,b) * a = b. +""" +rotmatbetween(a::Vec3{T}, b::Vec3{T}) where {T<:Real} = rotmatbetween(Float64, a, b) +function rotmatbetween(::Type{S}, a::Vec3{T}, b::Vec3{T}) where {T<:Real,S<:Real} + # TODO: Brian, is there a hidden assumption that a and b are normalized? + v = cross(a, b) + c = dot(a, b) + V = SMatrix{3,3,S,9}(0, v[3], -v[2], -v[3], 0, v[1], v[2], -v[1], 0) + R = I + V + V^2 * one(T) / (one(T) + c) + return SMatrix{3,3,S,9}(R) +end +rotmatbetween(a::SVector{3,T}, b::SVector{3,T}) where {T<:Real} = rotmatbetween(Float64, Vec3(a), Vec3(b)) +function rotmatbetween(type::Type{S}, a::SVector{3,T}, b::SVector{3,T}) where {T<:Real,S<:Real} + return rotmatbetween(type, Vec3(a), Vec3(b)) +end +export rotmatbetween + + +""" + rotmatd([S::Type], θ::T, ϕ::T, ψ::T) -> SMatrix{3,3,S} + +Returns the rotation matrix of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in degrees**. +""" +rotmatd(θ::T, ϕ::T, ψ::T) where {T<:Number} = rotmat(Float64, deg2rad(θ), deg2rad(ϕ), deg2rad(ψ)) +rotmatd(::Type{S}, θ::T, ϕ::T, ψ::T) where {T<:Number,S<:Real} = rotmat(S, deg2rad(θ), deg2rad(ϕ), deg2rad(ψ)) +export rotmatd +""" + rotmat([S::Type], θ::T, ϕ::T, ψ::T) -> SMatrix{3,3,S} + +Returns the rotation matrix of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in radians**. +""" +rotmat(θ::T, ϕ::T, ψ::T) where {T<:Number} = rotmat(Float64, θ, ϕ, ψ) +function rotmat(::Type{S}, θ::T, ϕ::T, ψ::T) where {T<:Number,S<:Real} + sinψ = sin(ψ) + sinϕ = sin(ϕ) + sinθ = sin(θ) + cosψ = cos(ψ) + cosϕ = cos(ϕ) + cosθ = cos(θ) + return SMatrix{3,3,S,9}(cosψ * cosϕ, sinψ * cosϕ, -sinϕ, cosψ * sinϕ * sinθ - sinψ * cosθ, sinψ * sinϕ * sinθ + cosψ * cosθ, cosϕ * sinθ, cosψ * sinϕ * cosθ + sinψ * sinθ, sinψ * sinϕ * cosθ - cosψ * sinθ, cosϕ * cosθ) +end +export rotmat + + +#endregion Transform + + + + diff --git a/src/Optical/Eye.jl b/src/Optical/Eye.jl index 0006d3ae8..cc138103e 100644 --- a/src/Optical/Eye.jl +++ b/src/Optical/Eye.jl @@ -21,7 +21,7 @@ # SOFTWARE """ - ModelEye(assembly::LensAssembly{T}, nsamples::Int = 17; pupil_radius::T = 3.0, detpixels::Int = 1000, transform::RigidBodyTransform{T} = identitytransform(T)) + ModelEye(assembly::LensAssembly{T}, nsamples::Int = 17; pupil_radius::T = 3.0, detpixels::Int = 1000, transform::Transform{T} = identitytransform(T)) Geometrically accurate model of the human eye focussed at infinity with variable `pupil_radius`. The eye is added to the provided `assembly` to create a [`CSGOpticalSystem`](@ref) with the retina of the eye as the detector. @@ -31,7 +31,7 @@ By default the eye is directed along the positive z-axis with the vertex of the `nsamples` determines the resolution at which accelerated surfaces within the eye are triangulated. """ -function ModelEye(assembly::LensAssembly{T}; nsamples::Int = 17, pupil_radius::T = 3.0, detpixels::Int = 1000, transform::RigidBodyTransform{T} = identitytransform(T), detdatatype::Type{D} = Float32) where {T<:Real,D<:Real} +function ModelEye(assembly::LensAssembly{T}; nsamples::Int = 17, pupil_radius::T = 3.0, detpixels::Int = 1000, transform::Transform{T} = identitytransform(T), detdatatype::Type{D} = Float32) where {T<:Real,D<:Real} cornea_front = AcceleratedParametricSurface(ZernikeSurface(6.9, radius = -7.8, conic = -0.5), nsamples, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.CORNEA, OpticSim.GlassCat.Air)) cornea_rear = Plane(SVector{3,T}(0, 0, -1), SVector{3,T}(0, 0, -3.2)) cornea_csg = csgintersection(cornea_front, cornea_rear) diff --git a/src/Optical/Lenses.jl b/src/Optical/Lenses.jl index 824b805a7..4b3c3f06d 100644 --- a/src/Optical/Lenses.jl +++ b/src/Optical/Lenses.jl @@ -119,7 +119,7 @@ function AsphericLens(insidematerial::T, frontvertex::S, frontradius::S, frontco else #conic or aspheric surf = AcceleratedParametricSurface(ZernikeSurface(semidiameter + backdecenter_l + S(0.01), radius = backradius, conic = backconic, aspherics = backaspherics), nsamples, interface = opticinterface(S, insidematerial, nextmaterial, backsurfacereflectance)) - lens_rear = leaf(surf, RigidBodyTransform{S}(zero(S), S(π), zero(S), backdecenter[1], backdecenter[2], frontvertex - thickness)) + lens_rear = leaf(surf, Transform{S}(zero(S), S(π), zero(S), backdecenter[1], backdecenter[2], frontvertex - thickness)) end extra_front = frontradius >= zero(S) || isinf(frontradius) ? zero(S) : abs(frontradius) - sqrt(frontradius^2 - semidiameter^2) extra_back = backradius >= zero(S) || isinf(backradius) ? zero(S) : abs(backradius) - sqrt(backradius^2 - semidiameter^2) diff --git a/src/Optical/OpticalRay.jl b/src/Optical/OpticalRay.jl index 30f037595..3a52f0808 100644 --- a/src/Optical/OpticalRay.jl +++ b/src/Optical/OpticalRay.jl @@ -97,6 +97,6 @@ function Base.print(io::IO, a::OpticalRay{T,N}) where {T,N} end end -function Base.:*(a::RigidBodyTransform{T}, r::OpticalRay{T,N}) where {T,N} +function Base.:*(a::Transform{T}, r::OpticalRay{T,N}) where {T,N} return OpticalRay(a * ray(r), power(r), wavelength(r), opl = pathlength(r), nhits = nhits(r), sourcenum = sourcenum(r), sourcepower = sourcepower(r)) end diff --git a/src/Visualization/Visualization.jl b/src/Visualization/Visualization.jl index fb787e826..8cb58e349 100644 --- a/src/Visualization/Visualization.jl +++ b/src/Visualization/Visualization.jl @@ -24,6 +24,7 @@ module Vis using ..OpticSim using ..OpticSim: euclideancontrolpoints, evalcsg, vertex, makiemesh, detector, centroid, origingen, lower, upper, intervals, α +using ..Geometry using Unitful using ImageView @@ -302,11 +303,11 @@ indexedcolor2(i::Int) = ColorSchemes.hsv[1.0 - rem(i / (2.1 * π), 1.0)] .* 0.5 ## displaying imported meshes -Base.:*(a::RigidBodyTransform, p::Makie.AbstractPlotting.GeometryBasics.PointMeta) = a * p.main +Base.:*(a::Transform, p::Makie.AbstractPlotting.GeometryBasics.PointMeta) = a * p.main Base.:*(a::Real, p::Makie.AbstractPlotting.GeometryBasics.PointMeta{N,S}) where {S<:Real,N} = Makie.AbstractPlotting.GeometryBasics.Point{N,S}((a * SVector{N,S}(p))...) -Base.:*(a::RigidBodyTransform, p::Makie.AbstractPlotting.GeometryBasics.Point{N,S}) where {S<:Real,N} = Makie.AbstractPlotting.GeometryBasics.Point{N,S}((a.rotation * SVector{N,S}(p) + a.translation)...) +Base.:*(a::Transform, p::Makie.AbstractPlotting.GeometryBasics.Point{N,S}) where {S<:Real,N} = Makie.AbstractPlotting.GeometryBasics.Point{N,S}((a.rotation * SVector{N,S}(p) + a.translation)...) -function draw!(scene::MakieLayout.LScene, ob::AbstractString; color = :gray, linewidth = 3, shaded::Bool = true, wireframe::Bool = false, transform::RigidBodyTransform{Float64} = identitytransform(Float64), scale::Float64 = 1.0, kwargs...) +function draw!(scene::MakieLayout.LScene, ob::AbstractString; color = :gray, linewidth = 3, shaded::Bool = true, wireframe::Bool = false, transform::Transform{Float64} = identitytransform(Float64), scale::Float64 = 1.0, kwargs...) if any(endswith(lowercase(ob), x) for x in [".obj", "ply", ".2dm", ".off", ".stl"]) meshdata = FileIO.load(ob) if transform != identitytransform(Float64) || scale != 1.0 @@ -964,7 +965,7 @@ function surfacesag(object::Union{CSGTree{T},Surface{T}}, resolution::Tuple{Int, end """ - eyebox_eval_eye(assembly::LensAssembly{T}, raygen::OpticalRayGenerator{T}, eye_rotation_x::T, eye_rotation_y::T, sample_points_x::Int, sample_points_y::Int; pupil_radius::T = T(2.0), resolution::Int = 512, eye_transform::RigidBodyTransform{T} = identitytransform(T)) + eyebox_eval_eye(assembly::LensAssembly{T}, raygen::OpticalRayGenerator{T}, eye_rotation_x::T, eye_rotation_y::T, sample_points_x::Int, sample_points_y::Int; pupil_radius::T = T(2.0), resolution::Int = 512, eye_transform::Transform{T} = identitytransform(T)) Visualise the images formed when tracing `assembly` with a human eye for an evenly sampled `sample_points_x × sample_points_y` grid in the angular range of eyeball rotations `-eye_rotation_x:eye_rotation_x` and `-eye_rotation_y:eye_rotation_y` in each dimension respectively. `resolution` is the size of the detector image (necessarily square). @@ -975,7 +976,7 @@ By default the eye is directed along the positive z-axis with the vertex of the The result is displayed as a 4D image - the image seen by the eye is shown in 2D as normal with sliders to vary eye rotation in x and y. The idea being that the whole image should be visible for all rotations in the range. """ -function eyebox_eval_eye(assembly::LensAssembly{T}, raygen::OpticalRayGenerator{T}, eye_rotation_x::T, eye_rotation_y::T, sample_points_x::Int, sample_points_y::Int; pupil_radius::T = T(2.0), resolution::Int = 512, eye_transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} +function eyebox_eval_eye(assembly::LensAssembly{T}, raygen::OpticalRayGenerator{T}, eye_rotation_x::T, eye_rotation_y::T, sample_points_x::Int, sample_points_y::Int; pupil_radius::T = T(2.0), resolution::Int = 512, eye_transform::Transform{T} = identitytransform(T)) where {T<:Real} out_image = zeros(Float32, resolution, resolution, sample_points_y, sample_points_x) xrange = sample_points_x > 1 ? range(-eye_rotation_x, stop = eye_rotation_x, length = sample_points_x) : [0.0] yrange = sample_points_y > 1 ? range(-eye_rotation_y, stop = eye_rotation_y, length = sample_points_y) : [0.0] @@ -986,7 +987,7 @@ function eyebox_eval_eye(assembly::LensAssembly{T}, raygen::OpticalRayGenerator{ r = OpticSim.rotmatd(T, ery, erx, 0.0) ov = OpticSim.rotate(eye_transform, SVector{3,T}(0.0, 0.0, 13.0)) t = r * ov - ov - sys = ModelEye(assembly, pupil_radius = pupil_radius, detpixels = resolution, transform = eye_transform * RigidBodyTransform(r, t)) + sys = ModelEye(assembly, pupil_radius = pupil_radius, detpixels = resolution, transform = eye_transform * Transform(r, t)) # Vis.draw(sys) # Vis.draw!((eye_transform.translation - ov - SVector(0.0, 20.0, 0.0), eye_transform.translation - ov + SVector(0.0, 20.0, 0.0))) # Vis.draw!((eye_transform.translation - ov - SVector(20.0, 0.0, 0.0), eye_transform.translation - ov + SVector(20.0, 0.0, 0.0))) diff --git a/test/TestData/TestData.jl b/test/TestData/TestData.jl index 04464b306..30dd83bb3 100644 --- a/test/TestData/TestData.jl +++ b/test/TestData/TestData.jl @@ -27,6 +27,7 @@ module TestData using OpticSim using OpticSim: tobeziersegments # try to use only exported functions so this list should stay short +using OpticSim.Geometry using StaticArrays using DataFrames using Unitful @@ -178,49 +179,49 @@ function zernikelens() topsurface = leaf(AcceleratedParametricSurface(ZernikeSurface(9.5, zcoeff = [(1, 0.0), (4, -1.0), (7, 0.5)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) - return csgintersection(barrel, csgintersection(topsurface, botsurface))(RigidBodyTransform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) + return csgintersection(barrel, csgintersection(topsurface, botsurface))(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) end function asphericlens() topsurface = leaf(AcceleratedParametricSurface(ZernikeSurface(9.5, aspherics = [(2, 0.0), (4, -0.001)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) - return csgintersection(barrel, csgintersection(topsurface, botsurface))(RigidBodyTransform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) + return csgintersection(barrel, csgintersection(topsurface, botsurface))(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) end function coniclensZ() topsurface = leaf(AcceleratedParametricSurface(ZernikeSurface(9.5, radius = -15.0, conic = -5.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) - return csgintersection(barrel, csgintersection(topsurface, botsurface))(RigidBodyTransform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) + return csgintersection(barrel, csgintersection(topsurface, botsurface))(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) end function coniclensQ() topsurface = leaf(AcceleratedParametricSurface(QTypeSurface(9.5, radius = -15.0, conic = -5.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) - return csgintersection(barrel, csgintersection(topsurface, botsurface))(RigidBodyTransform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) + return csgintersection(barrel, csgintersection(topsurface, botsurface))(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) end function qtypelens() topsurface = leaf(AcceleratedParametricSurface(QTypeSurface(9.0, radius = -25.0, conic = 0.3, αcoeffs = [(1, 0, 0.3), (1, 1, 1.0)], βcoeffs = [(1, 0, -0.1), (2, 0, 0.4), (3, 0, -0.6)], normradius = 9.5), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) - return csgintersection(barrel, csgintersection(topsurface, botsurface))(RigidBodyTransform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) + return csgintersection(barrel, csgintersection(topsurface, botsurface))(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) end function chebyshevlens() topsurface = leaf(AcceleratedParametricSurface(ChebyshevSurface(9.0, 9.0, [(1, 0, -0.081), (2, 0, -0.162), (0, 1, 0.243), (1, 1, 0.486), (2, 1, 0.081), (0, 2, -0.324), (1, 2, -0.405), (2, 2, -0.81)], radius = -25.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) - return csgintersection(barrel, csgintersection(topsurface, botsurface))(RigidBodyTransform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) + return csgintersection(barrel, csgintersection(topsurface, botsurface))(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) end function gridsaglens() topsurface = leaf(AcceleratedParametricSurface(GridSagSurface{Float64}(joinpath(@__DIR__, "../../test/test.GRD"), radius = -25.0, conic = 0.4, interpolation = GridSagBicubic), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) botsurface = leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -5.0, vishalfsizeu = 9.5, vishalfsizev = 9.5, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air))) barrel = leaf(Cylinder(9.0, 20.0, interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air, reflectance = zero(Float64), transmission = zero(Float64)))) - return csgintersection(barrel, csgintersection(topsurface, botsurface))(RigidBodyTransform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) + return csgintersection(barrel, csgintersection(topsurface, botsurface))(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) end ### SYSTEMS diff --git a/test/runtests.jl b/test/runtests.jl index 67b9705a9..4621aa78e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,6 +30,8 @@ using Unitful using Plots using OpticSim +# Geometry Module +using OpticSim.Geometry # curve/surface imports using OpticSim: findspan, makemesh, knotstoinsert, coefficients, inside, quadraticroots, tobeziersegments, evalcsg, makiemesh # interval imports diff --git a/test/testsets/General.jl b/test/testsets/General.jl index d14a8fa4e..958a302fe 100644 --- a/test/testsets/General.jl +++ b/test/testsets/General.jl @@ -19,7 +19,7 @@ @test similarroots(-1, -1, x1, x2) end # testset QuadraticRoots - @testset "RigidBodyTransform" begin + @testset "Transform" begin Random.seed!(SEED) @test isapprox(rotmatd(180, 0, 0), [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0], rtol = RTOLERANCE, atol = ATOLERANCE) @test isapprox(rotmatd(0.0, 180.0, 0.0), [-1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 -1.0], rtol = RTOLERANCE, atol = ATOLERANCE) @@ -29,16 +29,16 @@ x, y, z = rand(3) @test isapprox(rotmatd(x * 180 / π, y * 180 / π, z * 180 / π), rotmat(x, y, z), rtol = RTOLERANCE, atol = ATOLERANCE) - @test isapprox(RigidBodyTransform(rotmatd(0, 90, 0), SVector(0.0, 0.0, 1.0)) * SVector(1.0, 0.0, 0.0), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) + @test isapprox(Transform(rotmatd(0, 90, 0), SVector(0.0, 0.0, 1.0)) * SVector(1.0, 0.0, 0.0), [0.0, 0.0, 0.0], rtol = RTOLERANCE, atol = ATOLERANCE) - ta = RigidBodyTransform(rotmatd(0, 90, 0), SVector(0.0, 0.0, 1.0)) - tb = RigidBodyTransform(rotmatd(90, 0, 0), SVector(1.0, 0.0, 0.0)) - @test isapprox(collect(ta * tb), collect(RigidBodyTransform(rotmatd(90, 90, 0), SVector(0.0, 0.0, 0.0))), rtol = RTOLERANCE, atol = ATOLERANCE) + ta = Transform(rotmatd(0, 90, 0), SVector(0.0, 0.0, 1.0)) + tb = Transform(rotmatd(90, 0, 0), SVector(1.0, 0.0, 0.0)) + @test isapprox(collect(ta * tb), collect(Transform(rotmatd(90, 90, 0), SVector(0.0, 0.0, 0.0))), rtol = RTOLERANCE, atol = ATOLERANCE) @test isapprox(collect(ta * inv(ta)), collect(identitytransform()), rtol = RTOLERANCE, atol = ATOLERANCE) @test isapprox(collect(tb * inv(tb)), collect(identitytransform()), rtol = RTOLERANCE, atol = ATOLERANCE) - @test isapprox(collect(inv(ta)), collect(RigidBodyTransform(rotmatd(0, -90, 0), SVector(1.0, 0.0, 0.0))), rtol = RTOLERANCE, atol = ATOLERANCE) - end # testset RigidBodyTransform + @test isapprox(collect(inv(ta)), collect(Transform(rotmatd(0, -90, 0), SVector(1.0, 0.0, 0.0))), rtol = RTOLERANCE, atol = ATOLERANCE) + end # testset Transform @testset "Interval" begin pt1 = Intersection(0.3, [4.0, 5.0, 6.0], normalize(rand(3)), 0.4, 0.5, NullInterface()) diff --git a/test/testsets/Visualization.jl b/test/testsets/Visualization.jl index f2754f1a3..5bc5dd33c 100644 --- a/test/testsets/Visualization.jl +++ b/test/testsets/Visualization.jl @@ -4,7 +4,7 @@ # test that this all at least runs surf1 = AcceleratedParametricSurface(TestData.beziersurface(), 15) surf2 = AcceleratedParametricSurface(TestData.upsidedownbeziersurface(), 15) - m = csgintersection(leaf(surf1, translation(-0.5, -0.5, 0.0)), csgintersection(leaf(Cylinder(0.3, 5.0)), leaf(surf1, RigidBodyTransform{Float64}(0.0, Float64(π), 0.0, 0.5, -0.5, 0.0))))() + m = csgintersection(leaf(surf1, translation(-0.5, -0.5, 0.0)), csgintersection(leaf(Cylinder(0.3, 5.0)), leaf(surf1, Transform{Float64}(0.0, Float64(π), 0.0, 0.5, -0.5, 0.0))))() @test_nowarn Vis.draw(m) m = csgintersection(leaf(surf1, translation(-0.5, -0.5, 0.0)), csgintersection(leaf(Cylinder(0.3, 5.0)), leaf(surf2, translation(-0.5, -0.5, 0.0))))() @test_nowarn Vis.draw(m)