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)