diff --git a/Project.toml b/Project.toml
index 732f4ac..9a32beb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,15 +4,17 @@ authors = ["Steve Kelly <kd2cca@gmail.com>"]
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118"
 Meshing = "e6723b4c-ebff-59f1-b4b7-d97aa5274f73"
+Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-Meshing = "0.5"
 GeometryBasics = "0.2"
+Meshing = "0.5"
diff --git a/src/Descartes.jl b/src/Descartes.jl
index 440a73a..f546169 100644
--- a/src/Descartes.jl
+++ b/src/Descartes.jl
@@ -1,15 +1,18 @@
 module Descartes
 
 import Base: *, union, diff, intersect
-import GeometryBasics: HyperRectangle,
+import GeometryBasics: AbstractMesh,
+                       HyperRectangle,
+                       Vec,
                        Mesh
 
-using GeometryBasics,
-      FileIO,
+using FileIO,
       StaticArrays,
       Meshing,
       MeshIO,
-      LinearAlgebra
+      LinearAlgebra,
+      CoordinateTransformations,
+      Rotations
 
 include("types.jl")
 include("constructors.jl")
@@ -17,22 +20,15 @@ include("operations.jl")
 include("frep.jl")
 include("hyperrectangles.jl")
 include("meshes.jl")
-#include("hashing.jl")
-include("cache.jl")
-include("visualize.jl")
-
-# GeometryTypes
-export Point
 
 # 3d primitives
 export Cuboid, Cylinder, Sphere, Piping
 
 # 2d primitives
-
 export Square, Circle
 
 # transforms
-export Transform, Translation, translate, rotate
+export translate, rotate
 
 # operations
 export Shell, LinearExtrude
diff --git a/src/constructors.jl b/src/constructors.jl
index e913894..b681b65 100644
--- a/src/constructors.jl
+++ b/src/constructors.jl
@@ -3,7 +3,7 @@ function Square(dims;center=false)
     y = dims[2]
     lb = SVector(0.,0.)
     center && (lb = -SVector{2,Float64}(x,y)/2)
-    Square(SVector{2,Float64}(x,y), lb, SMatrix{3,3}(1.0*I), SMatrix{3,3}(1.0*I))
+    Square(SVector{2,Float64}(x,y), lb)
 end
 
 Square(x,y,z;center=false) = Square([x,y],center=center)
@@ -11,8 +11,8 @@ Square(x,y,z;center=false) = Square([x,y],center=center)
 function Cuboid(dims;center=false)
     @assert length(dims) == 3
     lb = SVector(0.,0.,0.)
-    center && (lb = -SVector{3,Float64}(dims)/2)
-    Cuboid(SVector{3,Float64}(dims), lb, SMatrix{4,4}(1.0*I), SMatrix{4,4}(1.0*I))
+    center && (lb = -SVector{3,Float64}(dims...)/2)
+    Cuboid(SVector{3,Float64}(dims...), lb)
 end
 
 Cuboid(x,y,z;center=false) = Cuboid([x,y,z],center=center)
@@ -26,16 +26,7 @@ function Cylinder(r, h; center=false)
     rn, hn = Float64(r), Float64(h)
     b = 0.0
     center && (b = -h/2)
-    Cylinder(rn, hn, b, SMatrix{4,4}(one(rn)*I), SMatrix{4,4}(one(rn)*I))
-end
-
-function Circle(r)
-    rn = Float64(r)
-    Circle(rn, SMatrix{3,3}(one(rn)*I), SMatrix{3,3}(one(rn)*I))
-end
-
-function Sphere(r)
-    Sphere(Float64(r), SMatrix{4,4}(1.0*I), SMatrix{4,4}(1.0*I))
+    Cylinder(rn, hn, b)
 end
 
 #function PrismaticCylinder(sides, height::T, radius::T) where {T}
@@ -44,7 +35,7 @@ end
 
 function Piping(r, pts)
     ptsc = [SVector{3,Float64}(pt...) for pt in pts]
-    Piping(Float64(r), ptsc, SMatrix{4,4}(1.0*I), SMatrix{4,4}(1.0*I))
+    Piping(Float64(r), ptsc)
 end
 
 #
@@ -130,5 +121,6 @@ function Shell(r)
 end
 
 function LinearExtrude(p::AbstractPrimitive{2,T}, d) where {T}
-    LinearExtrude{3, T, typeof(p)}(p, convert(T, d), SMatrix{4,4}(1.0*I), SMatrix{4,4}(1.0*I))
+    #TODO promote type
+    LinearExtrude{3, T, typeof(p)}(p, convert(T, d))
 end
diff --git a/src/frep.jl b/src/frep.jl
index d8efa8b..4a2e5cc 100644
--- a/src/frep.jl
+++ b/src/frep.jl
@@ -8,35 +8,37 @@ function _radius(a,b,r)
     end
 end
 
+function FRep(p::MapContainer{3,T,P}, v) where {T,P}
+    FRep(p.primitive, p.inv(v))
+end
+
+
+function FRep(p::MapContainer{2,T,P}, v) where {T,P}
+    FRep(p.primitive, p.inv(SVector(v[1],v[2])))
+end
+
 function FRep(p::Sphere, v)
-    it = p.inv_transform
-    x = v[1]*it[1,1]+v[2]*it[1,2]+v[3]*it[1,3]+it[1,4]
-    y = v[1]*it[2,1]+v[2]*it[2,2]+v[3]*it[2,3]+it[2,4]
-    z = v[1]*it[3,1]+v[2]*it[3,2]+v[3]*it[3,3]+it[3,4]
+    x = v[1]
+    y = v[2]
+    z = v[3]
     sqrt(x*x + y*y + z*z) - p.radius
 end
 
 function FRep(p::Cylinder, v)
-    it = p.inv_transform
-    x = v[1]*it[1,1]+v[2]*it[1,2]+v[3]*it[1,3]+it[1,4]
-    y = v[1]*it[2,1]+v[2]*it[2,2]+v[3]*it[2,3]+it[2,4]
-    z = v[1]*it[3,1]+v[2]*it[3,2]+v[3]*it[3,3]+it[3,4]
+    x = v[1]
+    y = v[2]
+    z = v[3]
     max(-z+p.bottom, z-p.height-p.bottom, sqrt(x*x + y*y) - p.radius)
 end
 
 function FRep(p::Circle, v)
-    it = p.inv_transform
-    x = v[1]*it[1,1]+v[2]*it[1,2]+it[1,3]
-    y = v[1]*it[2,1]+v[2]*it[2,2]+it[2,3]
-    sqrt(x*x + y*y) - p.radius
+    norm(v) - p.radius
 end
 
-
 function FRep(p::Cuboid, v)
-    it = p.inv_transform
-    x = v[1]*it[1,1]+v[2]*it[1,2]+v[3]*it[1,3]+it[1,4]
-    y = v[1]*it[2,1]+v[2]*it[2,2]+v[3]*it[2,3]+it[2,4]
-    z = v[1]*it[3,1]+v[2]*it[3,2]+v[3]*it[3,3]+it[3,4]
+    x = v[1]
+    y = v[2]
+    z = v[3]
     dx, dy, dz = p.dimensions
     lbx, lby,lbz = p.lowercorner
     max(-x+lbx, x-dx-lbx,
@@ -45,9 +47,8 @@ function FRep(p::Cuboid, v)
 end
 
 function FRep(p::Square, v)
-    it = p.inv_transform
-    x = v[1]*it[1,1]+v[2]*it[1,2]+it[1,3]
-    y = v[1]*it[2,1]+v[2]*it[2,2]+it[2,3]
+    x = v[1]
+    y = v[2]
     dx, dy = p.dimensions
     lbx, lby = p.lowercorner
     max(-x+lbx, x-dx-lbx,
@@ -80,7 +81,6 @@ end
 
 function FRep(p::Piping{T}, v) where {T}
     num_pts = length(p.points)
-    pt = Point(v[1],v[2],v[3])
 
     val = typemax(T)
 
@@ -88,13 +88,13 @@ function FRep(p::Piping{T}, v) where {T}
         e1 = p.points[i]
         e2 = p.points[i+1]
         v = e2 - e1
-        w = pt - e1
+        w = v - e1
         if dot(w,v) <= 0
-            nv = norm(pt - e1)
+            nv = norm(v - e1)
         elseif dot(v,v) <= dot(w,v)
-            nv = norm(pt - e2)
+            nv = norm(v - e2)
         else
-            nv = norm(cross(pt-e1,pt-e2))/norm(e2-e1)
+            nv = norm(cross(v-e1,v-e2))/norm(e2-e1)
         end
         val = min(nv, val)
     end
@@ -102,10 +102,9 @@ function FRep(p::Piping{T}, v) where {T}
 end
 
 function FRep(p::LinearExtrude, v)
-    it = p.inv_transform
-    x = v[1]*it[1,1]+v[2]*it[1,2]+v[3]*it[1,3]+it[1,4]
-    y = v[1]*it[2,1]+v[2]*it[2,2]+v[3]*it[2,3]+it[2,4]
-    z = v[1]*it[3,1]+v[2]*it[3,2]+v[3]*it[3,3]+it[3,4]
+    x = v[1]
+    y = v[2]
+    z = v[3]
     r = FRep(p.primitive, v)
     max(max(-z,z-p.distance), r)
 end
diff --git a/src/hyperrectangles.jl b/src/hyperrectangles.jl
index 9d521d5..eb5f37d 100644
--- a/src/hyperrectangles.jl
+++ b/src/hyperrectangles.jl
@@ -1,34 +1,28 @@
 function HyperRectangle(square::Square{T}) where {T}
-    orig = HyperRectangle{2,T}(square.lowercorner, square.dimensions)
-    transform(square.transform, orig)
+    HyperRectangle{2,T}(square.lowercorner, square.dimensions)
 end
 
 function HyperRectangle(cube::Cuboid{T}) where {T}
-    orig = HyperRectangle{3,T}(cube.lowercorner, cube.dimensions)
-    transform(cube.transform, orig)
+    HyperRectangle{3,T}(cube.lowercorner, cube.dimensions)
 end
 
 function HyperRectangle(sphere::Sphere{T}) where {T}
-    #TODO?
-    orig = HyperRectangle{3,T}(fill(-sphere.radius,3), fill(sphere.radius*2,3))
-    transform(sphere.transform, orig)
+    HyperRectangle{3,T}(fill(-sphere.radius,3), fill(sphere.radius*2,3))
 end
 
 function HyperRectangle(p::Cylinder{T}) where {T}
-    orig = HyperRectangle{3,T}(Vec(-p.radius,-p.radius,0), Vec(p.radius*2,p.radius*2,p.height))
-    transform(p.transform, orig)
+    HyperRectangle{3,T}(Vec(-p.radius,-p.radius,0), Vec(p.radius*2,p.radius*2,p.height))
 end
 
 function HyperRectangle(p::Circle{T}) where {T}
-    orig = HyperRectangle{2,T}(Vec(-p.radius,-p.radius), Vec(p.radius*2,p.radius*2))
-    transform(p.transform, orig)
+    HyperRectangle{2,T}(Vec(-p.radius,-p.radius), Vec(p.radius*2,p.radius*2))
 end
 
 function HyperRectangle(p::Piping{T}) where {T}
     maxx, maxy, maxz = typemin(Float64), typemin(Float64), typemin(Float64)
     minx, miny, minz = typemax(Float64), typemax(Float64), typemax(Float64)
     for pt in p.points
-        newx, newy, newz = transform(p.transform, pt)
+        newx, newy, newz = pt
         maxx = max(newx,maxx)
         maxy = max(newy,maxy)
         maxz = max(newz,maxz)
@@ -41,6 +35,10 @@ function HyperRectangle(p::Piping{T}) where {T}
     HyperRectangle(minv, maxv-minv)
 end
 
+function HyperRectangle(m::MapContainer)
+    transform(m.map, HyperRectangle(m.primitive))
+end
+
 function HyperRectangle(csg::CSGUnion)
     union(HyperRectangle(csg.left), HyperRectangle(csg.right))
 end
@@ -53,7 +51,7 @@ function HyperRectangle(csg::CSGIntersect)
     intersect(HyperRectangle(csg.left), HyperRectangle(csg.right))
 end
 
-function HyperRectangle(csg::CSGDiff{N, T, L, R}) where {N, T, L, R}
+function HyperRectangle(csg::CSGDiff)
     diff(HyperRectangle(csg.left), HyperRectangle(csg.right))
 end
 
diff --git a/src/operations.jl b/src/operations.jl
index 59a27ac..6df871f 100644
--- a/src/operations.jl
+++ b/src/operations.jl
@@ -2,96 +2,32 @@ function *(s::Shell{Nothing}, obj::AbstractPrimitive{N,T}) where {N,T}
     Shell(obj, s.distance)
 end
 
-function *(a::Transform{N,T}, b::Transform{N,T}) where {N,T}
-    Transform{N,T}(a.transform*b.transform)
+function translate(vect::AbstractVector)
+    Translation(vect...)
 end
 
-"""
-rotate around the given axis, by angle
-"""
-function rotate(ang, axis::Vector)
-    n = length(axis)
-    N = n+1 # homogenous coords
-
-    # handle undefined axis e.g. zero array
-    if norm(axis) == 0
-        return Transform{N,Float64}(SMatrix{N,N}(1.0*I))
-    end
-
-    u = normalize(axis) # normalize the specified axis
-
-    # TODO handle 2D case
-    if n == 3
-        R = [cos(ang)+u[1]^2*(1-cos(ang)) u[1]*u[2]*(1-cos(ang))-u[3]*sin(ang) u[1]*u[3]*(1-cos(ang))+u[2]*sin(ang) 0
-            u[2]*u[1]*(1-cos(ang))+u[3]*sin(ang) cos(ang)+u[2]^2*(1-cos(ang)) u[2]*u[3]*(1-cos(ang))-u[1]*sin(ang) 0
-            u[3]*u[1]*(1-cos(ang))-u[2]*sin(ang) u[3]*u[2]*(1-cos(ang))+u[1]*sin(ang) cos(ang)+u[3]^2*(1-cos(ang)) 0
-            0 0 0 1]
-            return Transform{4, Float64}(R)
-    else
-        error("Axis rotation or order $n, is not implemented yet")
-    end
-end
-
-function translate(vect::Vector)
-    n = length(vect)
-    N = n + 1
-    transform = MMatrix{N,N}(1.0*I)
-    transform[1:n, N] = vect
-    Transform{N, Float64}(SMatrix{N,N}(transform))
-end
-
-function translate(vect::SVector)
-    n = length(vect)
-    N = n + 1
-    transform = MMatrix{N,N}(1.0*I)
-    for i = 1:n
-        transform[i, N] = vect[i]
-    end
-    Transform{N, Float64}(SMatrix{N,N}(transform))
+function *(transform::AbstractAffineMap, obj::PT) where {PT<:AbstractPrimitive}
+    MapContainer(transform, inv(transform), obj)
 end
 
-function *(transform::Transform, obj::PT) where {PT<:AbstractPrimitive}
-    nt= transform.transform*obj.transform
-    nit = inv(nt)
-    PT((getfield(obj,i) for i in fieldnames(PT)[1:end-2])..., nt, nit)
+function *(t::AbstractAffineMap, obj::MapContainer)
+    c = compose(obj.map, t)
+    MapContainer(c, inv(c), obj.primitive)
 end
 
-# commute the transform over each leaf in the CSG Tree
-function *(transform::Transform, obj::CSGTy) where {CSGTy <: AbstractCSGTree}
-    l = transform*obj.left
-    r = transform*obj.right
-    CSGTy(l, r)
+function *(a::AbstractAffineMap, b::AbstractAffineMap)
+    compose(a,b)
 end
 
-function *(transform::Transform, obj::RadiusedCSGUnion)
-    l = transform*obj.left
-    r = transform*obj.right
-    RadiusedCSGUnion(obj.r, l, r)
-end
-
-"""
-Apply a homogenous tranform matrix to the points x, y, z
-"""
-function transform(it::SMatrix, _x, _y, _z) where {T}
-    @inbounds x = _x*it[1,1]+_y*it[1,2]+_z*it[1,3]+it[1,4]
-    @inbounds y = _x*it[2,1]+_y*it[2,2]+_z*it[2,3]+it[2,4]
-    @inbounds z = _x*it[3,1]+_y*it[3,2]+_z*it[3,3]+it[3,4]
-    x, y, z
-end
-
-function transform(it::SMatrix, x::AbstractVector) where {T}
-    transform(it, x...)
-end
-
-function transform(t::SMatrix, h::HyperRectangle{3,T}) where {T}
-    p_1 = t*SVector(h.origin... , 1)
-    p_2 = t*SVector(h.widths... , 1)
-    p_3 = t*SVector(h.origin[1]+h.widths[1],h.origin[2],h.origin[3], 1)
-    p_4 = t*SVector(h.origin[1],h.origin[2]+h.widths[2],h.origin[3], 1)
-    p_5 = t*SVector(h.origin[1],h.origin[2],h.origin[3]+h.widths[3], 1)
-    p_6 = t*SVector(h.origin[1]+h.widths[1],h.origin[2],h.origin[3]+h.widths[3], 1)
-    p_7 = t*SVector(h.origin[1],h.origin[2]+h.widths[2],h.origin[3]+h.widths[3], 1)
-    p_8 = t*SVector(h.origin[1]+h.widths[1],h.origin[2]+h.widths[2],h.origin[3], 1)
+function transform(t::AbstractAffineMap, h::HyperRectangle{3,T}) where {T}
+    p_1 = t(h.origin)
+    p_2 = t(h.widths)
+    p_3 = t(SVector(h.origin[1]+h.widths[1],h.origin[2],h.origin[3]))
+    p_4 = t(SVector(h.origin[1],h.origin[2]+h.widths[2],h.origin[3]))
+    p_5 = t(SVector(h.origin[1],h.origin[2],h.origin[3]+h.widths[3]))
+    p_6 = t(SVector(h.origin[1]+h.widths[1],h.origin[2],h.origin[3]+h.widths[3]))
+    p_7 = t(SVector(h.origin[1],h.origin[2]+h.widths[2],h.origin[3]+h.widths[3]))
+    p_8 = t(SVector(h.origin[1]+h.widths[1],h.origin[2]+h.widths[2],h.origin[3]))
     x_o = min(p_1[1],p_2[1],p_3[1],p_4[1],p_5[1],p_6[1],p_7[1],p_8[1])
     y_o = min(p_1[2],p_2[2],p_3[2],p_4[2],p_5[2],p_6[2],p_7[2],p_8[2])
     z_o = min(p_1[3],p_2[3],p_3[3],p_4[3],p_5[3],p_6[3],p_7[3],p_8[3])
@@ -101,11 +37,23 @@ function transform(t::SMatrix, h::HyperRectangle{3,T}) where {T}
     HyperRectangle(x_o, y_o, z_o, x_w-x_o, y_w-y_o, z_w-z_o)
 end
 
-function transform(t::SMatrix, h::HyperRectangle{2,T}) where {T}
-    p_1 = t*SVector(h.origin... , 1)
-    p_2 = t*SVector(h.widths... , 1)
-    p_3 = t*SVector(h.origin[1]+h.widths[1],h.origin[2], 1)
-    p_4 = t*SVector(h.origin[1],h.origin[2]+h.widths[2], 1)
+function rotate(ang, axis::Vector)
+    if !iszero(axis[1])
+        return LinearMap(RotX(ang))
+    elseif !iszero(axis[2])
+        return LinearMap(RotY(ang))
+    elseif !iszero(axis[3])
+        return LinearMap(RotZ(ang))
+    else
+        return LinearMap(RotX(ang)) # default to X? TODO: What does OpenSCAD use?
+    end
+end
+
+function transform(t::AbstractAffineMap, h::HyperRectangle{2,T}) where {T}
+    p_1 = t(h.origin)
+    p_2 = t(h.widths)
+    p_3 = t(SVector(h.origin[1]+h.widths[1],h.origin[2]))
+    p_4 = t(SVector(h.origin[1],h.origin[2]+h.widths[2]))
     x_o = min(p_1[1],p_2[1],p_3[1],p_4[1])
     y_o = min(p_1[2],p_2[2],p_3[2],p_4[2])
     x_w = max(p_1[1],p_2[1],p_3[1],p_4[1])
diff --git a/src/types.jl b/src/types.jl
index 0d74c50..14464b7 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -1,21 +1,16 @@
 abstract type AbstractPrimitive{N, T} end # dimensionality, numeric type
 abstract type AbstractCSGTree{N,T} <: AbstractPrimitive{N, T} end
 abstract type AbstractOperation{N,T} <: AbstractPrimitive{N,T} end
-abstract type AbstractTransform{N,T} end
 
 # 2D Geometric Primitives
 
 struct Square{T} <: AbstractPrimitive{2, T}
     dimensions::SVector{2,T}
     lowercorner::SVector{2,T}
-    transform::SMatrix{3,3,T,9}
-    inv_transform::SMatrix{3,3,T,9}
 end
 
 struct Circle{T} <: AbstractPrimitive{2, T}
     radius::T
-    transform::SMatrix{3,3,T,9}
-    inv_transform::SMatrix{3,3,T,9}
 end
 
 # 3D Geometric Primitives
@@ -23,35 +18,29 @@ end
 struct Cuboid{T} <: AbstractPrimitive{3, T}
     dimensions::SVector{3,T}
     lowercorner::SVector{3,T}
-    transform::SMatrix{4,4,T,16}
-    inv_transform::SMatrix{4,4,T,16}
 end
 
 struct Cylinder{T} <: AbstractPrimitive{3, T}
     radius::T
     height::T
     bottom::T
-    transform::SMatrix{4,4,T,16}
-    inv_transform::SMatrix{4,4,T,16}
 end
 
 struct Sphere{T} <: AbstractPrimitive{3, T}
     radius::T
-    transform::SMatrix{4,4,T,16}
-    inv_transform::SMatrix{4,4,T,16}
 end
 
 struct Piping{T} <: AbstractPrimitive{3, T}
     radius::T
     points::Vector{SVector{3,T}}
-    transform::SMatrix{4,4,T,16}
-    inv_transform::SMatrix{4,4,T,16}
 end
 
 # transforms
 
-struct Transform{N,T} <: AbstractTransform{N,T}
-    transform::SMatrix{N,N,T} #TODO: where L?
+struct MapContainer{N, T, M, P<:AbstractPrimitive{N,T}} <: AbstractPrimitive{N,T}
+    map::M
+    inv::M
+    primitive::P
 end
 
 
@@ -89,6 +78,4 @@ end
 struct LinearExtrude{N, T, P} <: AbstractPrimitive{3,T}
     primitive::P
     distance::T
-    transform::SMatrix{4,4,T,16}
-    inv_transform::SMatrix{4,4,T,16}
 end
diff --git a/test/tranforms.jl b/test/tranforms.jl
index 65d845f..b68fd54 100644
--- a/test/tranforms.jl
+++ b/test/tranforms.jl
@@ -2,15 +2,7 @@
 
 @testset "Rotations" begin
 
-r1 = rotate(pi/6,[1.2,0,0])
-r2 = rotate(pi/6,[0,3,0])
-r3 = rotate(pi/6,[0,0,5.6])
 
-sp = sin(pi/6)
-cp = cos(pi/6)
-@test all(r1.transform .≈ [1.0 0.0 0.0 0.0; 0.0 cp -sp 0.0; 0.0 sp cp 0.0; 0.0 0.0 0.0 1.0])
-@test all(r2.transform .≈ [cp 0.0 sp 0.0; 0.0 1.0 0.0 0.0; -sp 0.0 cp 0.0; 0.0 0.0 0.0 1.0])
-@test all(r3.transform .≈ [cp -sp 0.0 0.0; 0.5 cp 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0])
 end