diff --git a/Project.toml b/Project.toml
index f18bb9185..61e2aec9c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -43,7 +43,9 @@ PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StringEncodings = "69024149-9ee7-55f6-a4c4-859efe599b68"
@@ -84,7 +86,6 @@ PlutoUI = "0.7"
 Polynomials = "2.0"
 PyCall = "1.92"
 ReverseDiff = "1.7"
-Revise = "3.1"
 StaticArrays = "1.0"
 StringEncodings = "0.3"
 Unitful = "1.6"
diff --git a/src/Geometry/Primitives/NonCSG/ConvexPolygon.jl b/src/Geometry/Primitives/NonCSG/ConvexPolygon.jl
index cb9640e6a..3b8e6f414 100644
--- a/src/Geometry/Primitives/NonCSG/ConvexPolygon.jl
+++ b/src/Geometry/Primitives/NonCSG/ConvexPolygon.jl
@@ -69,6 +69,8 @@ end
 export ConvexPolygon
 
 centroid(poly::ConvexPolygon) = poly.plane.pointonplane
+normal(poly::ConvexPolygon) = normal(poly.plane)
+localframe(poly::ConvexPolygon) = poly.local_frame
 
 #function barrier to make vertices allocate less and be faster.
 function to3d(pts::SMatrix{2,N,T,L}) where{N,L,T}
diff --git a/src/Geometry/SphericalPolygon.jl b/src/Geometry/SphericalPolygon.jl
index eed2621f5..81a3d6ef0 100644
--- a/src/Geometry/SphericalPolygon.jl
+++ b/src/Geometry/SphericalPolygon.jl
@@ -2,7 +2,7 @@
 # Copyright (c) Microsoft Corporation. All rights reserved.
 # See LICENSE in the project root for full license information.
 
-#SphericalPolygon uses StaticArrays to represent vertices. Expect performance degradation for polygons with large numbers of vertices. Performance appears to be good up to perhaps 100 vertices, perhaps as much as 1000 vertices. By 10,000 vertices performance is terrible.
+"""SphericalPolygon uses StaticArrays to represent vertices. Expect performance degradation for polygons with large numbers of vertices. Performance appears to be good up to perhaps 100 vertices, perhaps as much as 1000 vertices. By 10,000 vertices performance is terrible."""
 struct SphericalPolygon{N,T<:Real}
     ptvectors::SMatrix{3,N,T}
     spherecenter::SVector{3,T}
diff --git a/src/Geometry/Transform.jl b/src/Geometry/Transform.jl
index 214722060..d6ca15010 100644
--- a/src/Geometry/Transform.jl
+++ b/src/Geometry/Transform.jl
@@ -135,6 +135,14 @@ Transform(rotation::AbstractArray{S,2}, translation::AbstractArray{S,1})
 """
 Transform{T} = SMatrix{4,4,T,16}
 export Transform
+#TODO: This was a bad idea. Should make Transform a concrete type of its own containing an SMatrix because the * operator for SMatrix*SVector is specialized for static arrays and vectors of particular sizes. Will lead to subtle and hard to fix bugs. Maybe something like this:
+# struct Transform{T} <: SMatrix{4,4,T,16}
+#     transform::SMatrix{4,4,T,16}
+# end
+# :*(a::Transform{T},b::SVector{3,T}) ...
+# :*(a::Transform{T},b::SVector{4,T}) ...
+
+
 
 # for compatability ith the "old" RigidBodyTransform
 """
@@ -433,6 +441,7 @@ function world2local(t::Transform{T}) where {T<:Real}
 end
 export world2local
 
+#TODO: should make Transform a concrete type of its own containing an SMatrix. This is hijackijng the * operator for SMatrix*SVector for arrays and vectors of particular sizes. Will lead to subtle and hard to fix bugs.
 function Base.:*(t::Transform{T}, v::SVector{3,T}) where {T<:Real}
     res = t * Vec4(v)
     if (t[4,4] == one(T))
@@ -442,7 +451,7 @@ function Base.:*(t::Transform{T}, v::SVector{3,T}) where {T<:Real}
     end
 end
 
-
+#TODO: should make Transform a concrete type of its own containing an SMatrix. This is hijackijng the * operator for SMatrix*SVector for arrays and vectors of particular sizes. Will lead to subtle and hard to fix bugs.
 function Base.:*(t::Transform{T}, m::SMatrix{3,N,T}) where{N,T<:Real}
     res = MMatrix{3,N,T}(undef)
 
diff --git a/src/GlassCat/utilities.jl b/src/GlassCat/utilities.jl
index fbeecf0c8..2560a1e8e 100644
--- a/src/GlassCat/utilities.jl
+++ b/src/GlassCat/utilities.jl
@@ -326,14 +326,20 @@ end
 Draw a scatter plot of index vs dispersion (the derivative of index with respect to wavelength). Both index and
 dispersion are computed at wavelength λ.
 
+Choose glasses to graph using the glassfilterprediate argument. This is a function that receives a Glass object and returns true if the glass should be graphed.
+
 If showprefixglasses is true then glasses with names like `F_BAK7` will be displayed. Otherwise glasses that have a
 leading letter prefix followed by an underscore, such as `F_`, will not be displayed.
 
 The index formulas for some glasses may give incorrect results if λ is outside the valid range for that glass. This can
 give anomalous results, such as indices less than zero or greater than 6. To filter out these glasses set maximumindex
 to a reasonable value such as 3.0.
+
+example: plot only glasses that do not contain the strings "E_" and "J_"
+
+drawglassmap(NIKON,showprefixglasses = true,glassfilterpredicate = (x) -> !occursin("J_",string(x)) && !occursin("E_",string(x)))
 """
-function drawglassmap(glasscatalog::Module; λ::Length = 550nm, glassfontsize::Integer = 3, showprefixglasses::Bool = false, minindex = 1.0, maxindex = 3.0, mindispersion = -.3, maxdispersion = 0.0)
+function drawglassmap(glasscatalog::Module; λ::Length = 550nm, glassfontsize::Integer = 3, showprefixglasses::Bool = false, minindex = 1.0, maxindex = 3.0, mindispersion = -.3, maxdispersion = 0.0, glassfilterpredicate = (x)->true)
     wavelength = Float64(ustrip(uconvert(μm, λ)))
     indices = Vector{Float64}(undef,0)
     dispersions = Vector{Float64}(undef,0)
@@ -351,7 +357,7 @@ function drawglassmap(glasscatalog::Module; λ::Length = 550nm, glassfontsize::I
 
             # don't show glasses that have an _ in the name. This prevents cluttering the map with many glasses of
             # similar (index, dispersion).
-            if (mindispersion <= dispersion <= maxdispersion) && (showprefixglasses || !hasprefix)
+            if glassfilterpredicate(glass) && (mindispersion <= dispersion <= maxdispersion) && (showprefixglasses || !hasprefix)
                 push!(indices, index(glass, wavelength))
                 push!(dispersions, dispersion)
                 push!(glassnames, String(name))
@@ -367,7 +373,8 @@ function drawglassmap(glasscatalog::Module; λ::Length = 550nm, glassfontsize::I
         series_annotations,
         markeralpha = 0.0,
         legends = :none,
-        xaxis = "dispersion",
+        xaxis = "dispersion @$λ",
         yaxis = "index",
-        title = "Glass Catalog: $glasscatalog") #should use markershape = :none to prevent markers from being drawn but this option doesn't work. Used markeralpha = 0 so the markers are invisible. A hack which works.
+        title = "Glass Catalog: $glasscatalog",
+        xflip = true) #should use markershape = :none to prevent markers from being drawn but this option doesn't work. Used markeralpha = 0 so the markers are invisible. A hack which works.
 end
diff --git a/src/Optical/Eye.jl b/src/Optical/HumanEye.jl
similarity index 70%
rename from src/Optical/Eye.jl
rename to src/Optical/HumanEye.jl
index e6d3d44df..e20b1bf65 100644
--- a/src/Optical/Eye.jl
+++ b/src/Optical/HumanEye.jl
@@ -2,12 +2,62 @@
 # Copyright (c) Microsoft Corporation. All rights reserved.
 # See LICENSE in the project root for full license information.
 
-using .GlassCat
+module HumanEye
+
+using OpticSim
+using OpticSim.GlassCat
+using OpticSim.Geometry
+import Unitful
+using Unitful.DefaultSymbols:mm,°
+
+# Approximate average properties of the human eye
+
+
+"""
+# Pupil diameter as a function of scene luminance
+https://jov.arvojournals.org/article.aspx?articleid=2279420
+https://en.wikipedia.org/wiki/Orders_of_magnitude_(luminance)
+
+Pupil diameter is approximately 2.8mm at 100cd/m^2. A typical overcast day is 700cd/m^2 
+"""
+
+"""computes pupil diameter as a function of scene luminance `L`, in cd/m², and the angular area, `a`, over which this luminance is presented to the eye."""
+𝐃sd(L,a) = 7.75 - 5.75 * ((L * a / 846)^.41) / ((L * a / 846)^.41 + 2) # the first letter of this function name is \bfD not D.
+
+eyeradius() = 12mm
+export eyeradius
+
+"""Posterior focal length, i.e., optical distance from entrance pupil to the retina. Focal length will change depending on accomodation. This value is for focus at ∞. When the eye is focused at 25cm focal length will be ≈ 22mm. Because the index of refraction of the vitreous humor is approximately 1.33 the physical distance from the entrance pupil to the retina will be 24mm/1.33 = 18mm."""
+eyefocallength() = 24mm
+export eyefocallength
+
+
+vc_epupil() = 3mm #distance from vertex of cornea to entrance pupil
+
+""" distance from vertex of cornea to center of rotation"""
+cornea_to_eyecenter() = 13.5mm
+export cornea_to_eyecenter
+
+entrancepupil_to_retina() = 22.9mm
+
+"""distance from entrance pupil to center of rotation."""
+entrancepupil_to_eyecenter() = entrancepupil_to_retina() - eyeradius()
+export entrancepupil_to_eyecenter
+
+"""average angle, in degrees, the eye will rotate before users will turn their head"""
+comfortable_eye_rotation_angle() = 20°
+
+"""average translation of the entrance pupil associated with comfortable eye rotation"""
+comfortable_entrance_pupil_translation() = sin(comfortable_eye_rotation_angle())*entrancepupil_to_eyecenter()
+export comfortable_entrance_pupil_translation
+
+
+#Model eyes of varying degrees of verisimilitude
 
 """
     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`.
+Geometrically accurate model of the human eye focused 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.
 
 The eye can be positioned in the scene using the `transform` argument and the resolution of the detector specified with `detpixels`.
@@ -76,4 +126,7 @@ function ArizonaEye(::Type{T} = Float64; accommodation::T = 0.0) where {T<:Real}
     )
 end
 export ArizonaEye
+
+end #module 
+export HumanEye
 #! format: on
diff --git a/src/Optical/Optical.jl b/src/Optical/Optical.jl
index 95212e099..aa70bdd2d 100644
--- a/src/Optical/Optical.jl
+++ b/src/Optical/Optical.jl
@@ -15,5 +15,5 @@ include("LensAssembly.jl")
 include("HierarchicalImage.jl")
 include("OpticalSystem.jl")
 include("Lenses.jl")
-include("Eye.jl")
+include("HumanEye.jl")
 include("Fields.jl")
diff --git a/src/Optical/Paraxial.jl b/src/Optical/Paraxial.jl
index 287c77935..3fbd43a25 100644
--- a/src/Optical/Paraxial.jl
+++ b/src/Optical/Paraxial.jl
@@ -91,6 +91,12 @@ function ParaxialLensHex(focaldistance::T, side_length::T, surfacenormal::SVecto
     return ParaxialLens(h, ParaxialInterface(focaldistance, centrepoint, outsidematerial))
 end
 
+function ParaxialLensConvexPoly(focallength::T,convpoly::ConvexPolygon{N,T},local_center_point::SVector{2,T}; outsidematerial::OpticSim.GlassCat.AbstractGlass = OpticSim.GlassCat.Air) where {N, T<:Real}
+    centrepoint = SVector{3, T}(local2world(localframe(convpoly)) * Vec3(local_center_point[1], local_center_point[2], zero(T)))
+    return ParaxialLens(convpoly, ParaxialInterface(focallength, centrepoint, outsidematerial))
+end
+
+
 function ParaxialLensConvexPoly(focaldistance::T, local_frame::Transform{T}, local_polygon_points::Vector{SVector{2, T}}, local_center_point::SVector{2, T}; outsidematerial::OpticSim.GlassCat.AbstractGlass = OpticSim.GlassCat.Air) where {N, T<:Real}
     poly = ConvexPolygon(local_frame, local_polygon_points)
     centrepoint = SVector{3, T}(local2world(local_frame) * Vec3(local_center_point[1], local_center_point[2], zero(T)))
diff --git a/src/RepeatingStructures/Cluster.jl b/src/RepeatingStructures/Cluster.jl
index dc50eb518..c91fc0ebc 100644
--- a/src/RepeatingStructures/Cluster.jl
+++ b/src/RepeatingStructures/Cluster.jl
@@ -68,6 +68,15 @@ function Base.size(::LatticeCluster{N1,N}) where{N1,N}
     return ntuple((i)->Base.IsInfinite(),N)
 end
 
+"""returns the integer coordinates of the tile at tileindex in the cluster at for a cluster with location cᵢ,cⱼ"""
+function tilecoordinates(cluster::LatticeCluster{N1,N},cᵢ,cⱼ,tileindex) where{N1,N}
+    @assert tileindex <= N1
+    clustercenter = basismatrix(clusterbasis(cluster))*SVector(cᵢ,cⱼ)
+    tileoffset = cluster.clusterelements[tileindex]
+    return clustercenter .+ tileoffset
+end
+export tilecoordinates
+
 Base.setindex!(A::AbstractBasis{N}, v, I::Vararg{Int, N}) where{N} = nothing #can't set lattice points. Might want to throw an exception instead.
 
 """ returns the lattice indices of the elements in the cluster. These are generally not the positions of the elements"""
@@ -77,10 +86,33 @@ function clustercoordinates(a::LatticeCluster{N1,N},indices::Vararg{Int,N}) wher
     for i in 1:N1
         temp[:,i] = SVector{N,Int}(a.clusterelements[i]) + clusteroffset
     end
-    return SMatrix{N,N1,Int}(temp) #positions of the cluster elements, not the lattice coordinates.
+    return SMatrix{N,N1,Int}(temp) #the lattice coordinates in the underlying element basis of the cluster elements
 end
 export clustercoordinates
 
+"""Given the (i,j) coordinates of a tile defined in the the underlying lattice basis of elements of the cluster compute the coordinates (cᵢ,cⱼ) of the cluster containing the tile, and the tile number of the tile in that cluster"""
+function cluster_coordinates_from_tile_coordinates(cluster::LatticeCluster{N1,N},i::Int,j::Int) where{N1,N}
+    found = false
+    bmatrix = Rational.(basismatrix(clusterbasis(cluster)))
+
+    local clustercoords::SVector{}
+    tileindex = 0 #initialize to illegal value
+    for coords in eachcol(clustercoordinates(cluster, 0,0)) #don't know which tile in the cluster the i,j coords come from so try all of them until find one that yields an integer value for the cluster coordinate.
+        tileindex += 1
+        possiblecenter = SVector{N,Rational}((i,j)) .- Rational.(coords)
+        clustercoords = bmatrix \  possiblecenter
+        if all(1 .== denominator.(clustercoords)) #If coordinates are integer this is the correct cluster value. If coordinates are not integer keep trying.
+            found = true    #every tile position i,j corresponds to some cluster (cᵢ,cⱼ) so will always find an answer
+            break
+        end
+    end
+    @assert found == true #should always find a cluster corresponding to any i,j. If not then something is seriously wrong.
+
+    return Int64.(clustercoords),tileindex  #return coordinates of the cluster and the index of the tile within that cluster
+end
+export cluster_coordinates_from_tile_coordinates
+
+
 
 """
 May want to have many properties associated with the elements in a cluster, which is why properties is represented as a DataFrame. The DataFrame in the properties field should have as many rows as there are elements in a cluster. At a minimum it must have a :Color and a :Name column.
diff --git a/src/RepeatingStructures/Lattice.jl b/src/RepeatingStructures/Lattice.jl
index a6e627681..dd7bc32c9 100644
--- a/src/RepeatingStructures/Lattice.jl
+++ b/src/RepeatingStructures/Lattice.jl
@@ -37,20 +37,23 @@ Returns the vertices of the unit polygon that tiles the plane for the basis
 tilevertices(a::S) where{S<:AbstractBasis}
 ```
 """
-abstract type AbstractBasis{N,T<:Real} end
+abstract type AbstractBasis{N,T <: Real} end
 export AbstractBasis
 
-function Base.getindex(A::B1, indices::Vararg{Int, N}) where{N,T,B1<:AbstractBasis{N,T}}
-    return basismatrix(A)*SVector{N,Int}(indices)
+function Base.getindex(A::B1, indices::Vararg{Int,N}) where {N,T,B1 <: AbstractBasis{N,T}}
+    return basismatrix(A) * SVector{N,Int}(indices)
 end
 
 #     temp::SVector{N,T} = (basis(A)*SVector{N,Int}(indices))::SVector{N,T}
 #     return temp
 # end
 
-Base.setindex!(A::AbstractBasis{N,T}, v, I::Vararg{Int, N}) where{T,N} = nothing #can't set lattice points. Might want to throw an exception instead.
+Base.setindex!(A::AbstractBasis{N,T}, v, I::Vararg{Int,N}) where {T,N} = nothing # can't set lattice points. Might want to throw an exception instead.
 
-struct LatticeBasis{N,T<:Real} <: AbstractBasis{N,T}
+"""computes the tile vertices at the offset lattice coordinate."""
+tilevertices(latticecoords::Union{AbstractVector,NTuple},lattice::AbstractBasis) = lattice[latticecoords...] .+ tilevertices(lattice)
+
+struct LatticeBasis{N,T <: Real} <: AbstractBasis{N,T}
     basisvectors::SMatrix{N,N,T}
 
     """ Convenience constructor that lets you use tuple arguments to describe the basis instead of SVector 
@@ -61,10 +64,10 @@ struct LatticeBasis{N,T<:Real} <: AbstractBasis{N,T}
     ```
     This allocates 48 bytes for a 2D basis in Julia 1.6.2. It  shouldn't allocate anything but not performance critical.
     """
-    function LatticeBasis(vectors::Vararg{NTuple{N,T},N}) where{T,N}  
+    function LatticeBasis(vectors::Vararg{NTuple{N,T},N}) where {T,N}  
         temp = MMatrix{N,N,T}(undef)
- 
-        for (j,val) in pairs(vectors)
+        
+        for (j, val) in pairs(vectors)
             for i in 1:N
                 temp[i,j] = val[i]
             end
@@ -73,15 +76,101 @@ struct LatticeBasis{N,T<:Real} <: AbstractBasis{N,T}
         return new{N,T}(SMatrix{N,N,T}(temp))
     end
 
-     LatticeBasis(vectors::SMatrix{N,N,T}) where{N,T<:Real} = new{N,T}(vectors)
+     LatticeBasis(vectors::SMatrix{N,N,T}) where {N,T <: Real} = new{N,T}(vectors)
 end
 export LatticeBasis
 
-function basismatrix(a::LatticeBasis{N,T})::SMatrix{N,N,T,N*N} where{N,T} 
+function basismatrix(a::LatticeBasis{N,T})::SMatrix{N,N,T,N * N} where {N,T} 
     return a.basisvectors
 end
 
 """Can access any point in a lattice so the range of indices is unlimited"""
-function Base.size(a::LatticeBasis{N,T}) where{N,T}
-    return ntuple((i)->Base.IsInfinite(),N)
-end
\ No newline at end of file
+function Base.size(a::LatticeBasis{N,T}) where {N,T}
+    return ntuple((i) -> Base.IsInfinite(), N)
+end
+
+# These functions are used to find the lattice tiles that lie inside a polygonal shape
+
+"""computes the matrix that transforms the basis set into the canonical basic vectors, eᵢ"""
+basisnormalization(a::Repeat.AbstractBasis) = Matrix(inv(Repeat.basismatrix(a))) # unfortunately LazySets doesn't work with StaticArrays. If you return a static Array is messes with their functions.
+
+"""warp the matrix into a coordinate frame where the basis vectors are canonical, eᵢ"""
+function normalizedshape(a::Repeat.AbstractBasis, shape::LazySets.VPolygon) 
+    normalizer = basisnormalization(a)
+    return normalizer * shape
+end
+export normalizedshape
+
+""" compute a conservative range of lattice indices that might intersect containingshape"""
+function latticebox(containingshape::LazySets.VPolygon, lattice::Repeat.AbstractBasis)
+    warpedshape = normalizedshape(lattice, containingshape)
+    box = LazySets.interval_hull(warpedshape)
+    return LazySets.Hyperrectangle(round.(box.center), (1, 1) .+ ceil.(box.radius)) # center the hyperrectangle on an integer point and enlarge the radius to take this into account
+end
+export latticebox
+
+
+"""Returns the lattice coordinates of all tiles that lie at least partly inside containingshape. 
+
+Compute the lattice tiles with non-zero intersection with containingshape. First compute a transformation that maps the lattice basis vectors to canonical unit basis vectors eᵢ (this is the inverse of the lattice basis matrix). Then transform containingshape into this coordinate frame and compute a bounding box. Unit steps along the coordinate axes in this space represent unit *lattice* steps in the original space. This makes it simple to determine coordinate bounds in the original space. Then test for intersection in the original space."""
+function tilesinside(containingshape::LazySets.VPolygon, lattice::Repeat.AbstractBasis{2,T}) where {T}
+    coordtype = Int64
+
+    box = latticebox(containingshape, lattice)
+    
+    coords = coordtype.(box.radius)
+    tilevertices = LazySets.VPolygon(Matrix(Repeat.tilevertices(lattice))) # VPolygon will accept StaticArrays but other LazySets function will barf.
+    result = Matrix{coordtype}(undef, 2, 0)
+    
+    for i in -coords[1]:coords[1]
+        for j in -coords[2]:coords[2] 
+            offsetcoords =   (Int64.(box.center) + [i,j])    
+            center = Vector(lattice[offsetcoords...]) #[] indexer for lattices returns SVector which LazySets doesn't handle. Convert to conventional Vector.
+            offsethex = LazySets.translate(tilevertices, center)
+
+            if !isempty(offsethex ∩ containingshape)
+                 result = hcat(result, offsetcoords)
+            end
+        end
+    end
+    return result
+end
+export tilesinside
+
+"""The vertices of the containing shape are the columns of the matrix containingshape"""
+tilesinside(containingshape::AbstractMatrix,lattice::Repeat.AbstractBasis) = tilesinside(LazySets.VPolygon(containingshape), lattice)
+
+tilesinside(xmin::T,ymin::T,xmax::T,ymax::T,lattice::Repeat.AbstractBasis) where {T <: Real} = tilesinside(SMatrix{2,4,T}(xmin, ymin, xmin, ymax, xmax, ymax, xmax, ymin), lattice)
+
+using Plots
+
+""" to see what the objects look like in the warped coordinate frame use inv(basismatrix(lattice)) as the transform"""
+function plotall(containingshape, lattice, transform=[1.0 0.0;0.0 1.0])
+    temp = Vector{SMatrix}(undef, 0)
+
+    for coordinate in eachcol(tilesinside(containingshape, lattice))
+        println(typeof(coordinate))
+        push!(temp, tilevertices(coordinate, lattice))
+    end
+
+    tiles = LazySets.VPolygon.(temp)
+    for tile in tiles
+        plot!(transform * tile, aspectratio=1)
+    end
+    plot!(containingshape, aspectratio=1)
+end
+export plotall
+
+
+function testtilesinside()
+    # poly = LazySets.VPolygon([-3.0 3.0 3.0; -3.0 3.0 -3.0])
+    hex = HexBasis1()
+    poly = LazySets.VPolygon(3 * tilevertices(hex))
+    tilecoords = tilesinside(poly, hex)
+     insidetiles = LazySets.VPolygon.([tilevertices(x,hex) for x in eachcol(tilecoords)])
+    # plotall(poly, hex)
+    for tile in insidetiles
+        plot!(tile)
+end
+end
+export testtilesinside
diff --git a/src/RepeatingStructures/Multilens/Analysis.jl b/src/RepeatingStructures/Multilens/Analysis.jl
new file mode 100644
index 000000000..769435222
--- /dev/null
+++ b/src/RepeatingStructures/Multilens/Analysis.jl
@@ -0,0 +1,346 @@
+
+# MIT license
+# Copyright (c) Microsoft Corporation. All rights reserved.
+# See LICENSE in the project root for full license information.
+
+
+using Roots
+import DataFrames
+
+
+# luminance (cd/m2)	Multiple	Value	Item
+# 10−6	µcd/m2	1 µcd/m2	Absolute threshold of vision[1]
+# 10−5			
+# 10−4		400 µcd/m2	Darkest sky[2]
+# 10−3	mcd/m2	1 mcd/m2	Night sky[3]
+# 1.4   mcd/m2	Typical photographic scene lit by full moon[4]
+# 5     mcd/m2	Approximate scotopic/mesopic threshold[5]
+# 10−2		40 mcd/m2	Phosphorescent markings on a watch dial after 1 h in the dark[6][7]
+# 10−1			
+# 100	cd/m2	2 cd/m2	Floodlit buildings, monuments, and fountains[8]
+# 5     cd/m2	Approximate mesopic/photopic threshold[5]
+# 101		25 cd/m2	Typical photographic scene at sunrise or sunset[4]
+# 30    cd/m2	Green electroluminescent source[2]
+# 102		250 cd/m2	Peak luminance of a typical LCD monitor[10][11]
+# 700   cd/m2	Typical photographic scene on overcast day[4][8][11]
+# 103	kcd/m2	2 kcd/m2	Average cloudy sky[2]
+# 5     kcd/m2	Typical photographic scene in full sunlight[4][8]
+
+
+
+"""This function returns the radius of the longest basis vector of the lattice cluster. Most lattices defined in this project have symmetric basis vectors so the radii of all basis vectors will be identical. This function is used in determing which clusters "fit" within the eye pupil."""
+function latticediameter(basismatrix::SMatrix)
+    maximum(norm.([basismatrix[:,i] for i in 1:size(basismatrix)[2]]))
+end
+export latticediameter
+
+latticediameter(a::Repeat.AbstractLatticeCluster) =   latticediameter(Repeat.basismatrix(Repeat.clusterbasis(a)))
+latticediameter(a::Repeat.AbstractBasis) =  latticediameter(Repeat.basismatrix(a))
+
+const hex3latticeclusterbasis = [2 // 1 -1 // 1;-1 // 1 2 // 1]
+export hex3latticeclusterbasis
+
+"""Only will work properly if lattice basis matrix contains only integer or rational terms. Returns the integer lattice coords of point in the given basis if the point is in the span of latticebasis. Otherwise returns nothing"""
+function latticepoint(latticebasis::AbstractMatrix, origin, point) 
+    Ainv = inv(Rational.(latticebasis))
+    b = [(point .- origin)...]
+    x = Ainv * b
+    if reduce(&, (1, 1) .== denominator.(x))
+        return Integer.(x)
+    else
+        return nothing
+    end
+end
+export latticepoint
+
+colorbasis(::Repeat.HexBasis1) = SMatrix{2,2}(2, -1, 1, 1)
+colorbasis(::Repeat.HexBasis3) = SMatrix{2,2}(2, -1, 1, 1)
+colororigins(::Repeat.HexBasis1) = ((0, 0), (-1, 0), (-1, 1))
+colororigins(::Repeat.HexBasis3) = ((0, 0), (0, -1), (1, -1))
+
+"""To reduce color channel cross talk it may be useful to arrange the color lenslets in a hexagonal lattice so that each color is surrounded by lenslets of a different color. With appropriate color filtering crosstalk between immediately neighboring lenslets can be reduced. This function computes the color which should be assigned to any lattice point to ensure this properly holds."""
+function pointcolor(point, cluster::Repeat.AbstractLatticeCluster)
+    latticematrix = colorbasis(Repeat.elementbasis(cluster))
+    origins = colororigins(Repeat.elementbasis(cluster))
+    colors = zip(origins, ("red", "green", "blue"))
+    for (origin, color) in colors
+        if nothing !== latticepoint(latticematrix, origin, point)
+            return color
+        end
+    end
+end
+export pointcolor
+
+defaultclusterproperties() = (mtf = .2, minfnumber = 2.0, cyclesperdegree = 11, λ = 530nm, pixelpitch = .9μm)
+export defaultclusterproperties
+
+"""maximum allowable display for visibility reasons"""
+maxlensletdisplaysize() = (250μm, 250μm)
+
+
+const ρ_quartervalue = 2.21509 # value of ρ at which the airy disk function has magnitude .25
+export ρ_quartervalue
+const ρ_zerovalue = 3.832 # value of ρ at which the airy disk function has magnitude 0
+
+"""given pixelpitch and angular subtense (in degrees) of pixel returns focal length"""
+focallength(pixelpitch::Unitful.Length,θ) = uconvert(mm, .5 * pixelpitch / tand(θ / 2))
+export focallength
+
+"""returns the diffraction limit frequency in cycles/deg
+
+focal length = 𝒇𝒍
+diffraction cutoff frequency,fc, in cycles/mm = 1/λF# = diameter/λ*𝒇𝒍
+cutoff wavelength, Wc, = 1/cutoff frequency = λ*𝒇𝒍/diameter
+
+angular wavelength, θc, radians/cycle, corresponding to cutoff wavelength:
+
+θ ≈ tanθ for small θ
+θc corresponding to Wc:  θc ≈ tanθ = Wc/𝒇𝒍
+Wc = θc*𝒇𝒍
+
+from equation for Wc:
+
+λ*𝒇𝒍/diameter = Wc = θc*𝒇𝒍
+θc = λ/diameter
+cycles/rad = 1/θc = diameter/λ
+
+"""
+cyclesperdegree(diameter::Unitful.Length,λ::Unitful.Length) = uconvert(Unitful.NoUnits, diameter / (rad2deg(1) * λ))
+export cyclesperdegree
+
+"""diffraction limited response for a circular aperture, normalized by maximum cutoff frequency"""
+function mtfcircular(freq, freqcutoff) 
+    s =  freq / freqcutoff
+    return 2 / π * (acos(s) - s * ((1 - s^2)^.5))
+end
+export mtfcircular
+
+"""returns the diffraction limit frequency in cycles/degree. At this frequeny the response of the system is zero"""
+diffractionlimit(λ::Unitful.Length,diameter::Unitful.Length) = uconvert(Unitful.NoUnits, diameter / λ) / rad2deg(1)
+export diffractionlimit
+
+"""computes the diameter of a diffraction limited lens that has response mtf at frequency cyclesperdeg"""
+function diameter_for_cycles_deg(mtf, cyclesperdeg, λ::Unitful.Length)
+    cyclesperrad = rad2deg(1) * cyclesperdeg
+    f(x) = mtfcircular(x, 1.0) - mtf
+    normalizedfrequency = find_zero(f, (0.0, 1.0))
+    fc = cyclesperrad / normalizedfrequency
+    return uconvert(mm, λ * fc)
+end
+export diameter_for_cycles_deg
+
+airydisk(ρ) = (2 * SpecialFunctions.besselj1(ρ) / ρ)^2
+
+"""The f# required for the first zero of the airy diffraction disk to be at the next sample point"""
+diffractionfnumber(λ::Unitful.Length,pixelpitch::Unitful.Length,indexofrefraction) = uconvert(Unitful.NoUnits, pixelpitch / (2.44λ / indexofrefraction))
+export diffractionfnumber
+
+"""returns ρ, the normalized distance, at which the airy disk will have the value airyvalue"""
+function ρatairyvalue(airyvalue)
+    ρ = 1e-5 # start at small angle
+    while airydisk(ρ) > airyvalue
+        ρ += 1e-5 # numerical precision not an issue here
+    end
+    return ρ
+end
+export ρatairyvalue
+
+""" Spacing between lenslets which guarantess that for any pixel visible in all lenslets every point on the eyebox plane is covered. This is closest packing of circles."""
+closestpackingdistance(pupildiameter::Unitful.Length) = pupildiameter * cosd(30)
+export closestpackingdistance
+
+"""Tries clusters of various sizes to choose the largest one which fits within the eye pupil. Larger clusters allow for greater reduction of the fov each lenslet must cover so it returns the largest feasible cluster"""
+function choosecluster(pupildiameter::Unitful.Length, lensletdiameter::Unitful.Length)
+    clusters = (hex3RGB(), hex4RGB(), hex7RGB(), hex9RGB(), hex12RGB(), hex19RGB())
+    # cdist = closestpackingdistance(pupildiameter)
+    cdist = pupildiameter
+    maxcluster = clusters[1]
+    ratio = 0.0
+
+    for cluster in clusters
+        temp = cdist / (lensletdiameter * latticediameter(cluster))
+            if temp >= 1.0
+            ratio = temp
+            maxcluster = cluster
+        else
+            break
+        end
+    end
+
+    @assert ratio ≥ 1.0 "ratio $ratio cdist $cdist lensletdiameter $lensletdiameter latticediameter $(latticediameter(maxcluster)) scaled=$(lensletdiameter * latticediameter(maxcluster))"
+
+    return (cluster = maxcluster, lensletdiameter = lensletdiameter * ratio, diameteroflattice = latticediameter(maxcluster) / ratio, packingdistance = cdist * ustrip(mm, lensletdiameter))
+end
+export choosecluster
+
+"""Computes the largest feasible cluster size meeting constraints"""
+function choosecluster(pupildiameter::Unitful.Length, λ::Unitful.Length, mtf, cyclesperdeg::T) where {T <: Real} 
+    diam = diameter_for_cycles_deg(mtf, cyclesperdeg, λ)
+    return choosecluster(pupildiameter, diam) # use minimum diameter for now.
+end
+
+"""Computes display size assuming lenslet normal to eyebox plane passes through lenslet. This is an approximation but for the narrow fov we are considering it is accurate enough to estimate pixel redundandcy, etc."""
+function sizeoflensletdisplay(eyerelief::T, eyebox::AbstractVector{T}, ppd, pixelpitch::S) where {T <: Unitful.Length,S <: Unitful.Length}
+    θ = eyeboxangles(eyebox, eyerelief)
+    displaysize = @.  θ * ppd * pixelpitch
+    return uconvert.(mm, displaysize)
+end
+export sizeoflensletdisplay
+
+"""This computes the approximate size of the entire display, not the individual lenslet displays."""
+sizeofdisplay(fov,eyerelief) = @. 2 * tand(fov / 2) * eyerelief
+export sizeofdisplay
+
+"""Computes the number of lenslets required to create a specified field of view at the given eyerelief"""
+function numberoflenslets(fov, eyerelief::Unitful.Length, lensletdiameter::Unitful.Length)
+    lensletarea = π * (lensletdiameter / 2)^2
+    dispsize = sizeofdisplay(fov, eyerelief)
+     return dispsize[1] * dispsize[2] / lensletarea
+end
+export numberoflenslets
+
+
+"""given the angles each lenslet has to cover compute the corresponding display size"""
+sizeoflensletdisplay(angles,ppd,pixelpitch::Unitful.Length) = @. angles * ppd * pixelpitch
+
+"""angular size of the eyebox when viewed from distance eyerelief"""
+eyeboxangles(eyebox,eyerelief) = @. atand(uconvert(Unitful.NoUnits, eyebox / eyerelief))
+export eyeboxangles
+
+"""computes how the fov can be subdivided among lenslets based on cluster size. Assumes the horizontal size of the eyebox is larger than the vertical so the larger number of subdivisions will always be the first number in the returns Tuple."""
+function anglesubdivisions(pupildiameter::Unitful.Length, λ::Unitful.Length, mtf, cyclesperdegree;RGB=true)
+    cluster, _ = choosecluster(pupildiameter, λ, mtf, cyclesperdegree)
+    numelements = Repeat.clustersize(cluster)
+    if numelements == 19
+        return RGB ? (3, 2) : (5, 3)
+    elseif numelements == 12
+        return RGB ? (2, 2) : (4, 3)
+   elseif numelements == 9
+        return RGB ? (3, 1) : (3, 3)
+    elseif numelements == 7
+        return RGB ? (3, 1) : (3, 2) 
+    elseif numelements == 4 || numelements == 3
+        return RGB ? (1, 1) : (3, 1)
+    else throw(ErrorException("this should never happen"))
+    end
+end
+export anglesubdivisions
+
+
+"""Computes the approximate fov required of each lenslet for the given constraints. This is strictly correct only for a lenslet centered in front of the eyebox, but the approximation is good enough for high level analysis"""
+function lensletangles(eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, pupildiameter::Unitful.Length, ppd; clusterproperties=defaultclusterproperties(), RGB=true)
+    cyclesperdegree = ppd / 2.0
+    return eyeboxangles(eyebox, eyerelief) ./ anglesubdivisions(pupildiameter, clusterproperties.λ, clusterproperties.mtf, clusterproperties.cyclesperdegree, RGB=RGB)
+end
+export lensletangles
+
+testangles() = lensletangles(18mm, (10mm, 6mm), 4mm, 45)
+export testangles
+
+lensletresolution(angles,ppd) = angles .* ppd
+
+"""Multilens displays tradeoff pixel redundancy for a reduction in total track of the display, by using many short focal length lenses to cover the eyebox. This function computes the ratio of pixels in the multilens display vs. a conventional display of the same nominal resolution"""
+function pixelredundancy(fov, eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, pupildiameter::Unitful.Length, ppd; RGB=true)
+    lensprops = defaultclusterproperties()
+    clusterdata = choosecluster(pupildiameter, lensprops.λ, lensprops.mtf, lensprops.cyclesperdegree)
+    nominalresolution = fov .* ppd
+    angles = lensletangles(eyerelief, eyebox, pupildiameter, ppd, RGB=RGB)
+    pixelsperlenslet = angles .* ppd
+    numlenses = numberoflenslets(fov, eyerelief, clusterdata.lensletdiameter)
+    return (numlenses * pixelsperlenslet[1] * pixelsperlenslet[2]) / ( nominalresolution[1] * nominalresolution[2])
+end
+export pixelredundancy
+
+testpixelredundancy() = pixelredundancy((55, 35), 18mm, (10mm, 6mm), 4mm, 45, RGB=false)
+export testpixelredundancy
+
+label(color) = color ? "RGB" : "Monochrome"
+
+"""generates a contour plot showing pixel redundancy as a function of ppd and pupil diameter"""
+function redundancy_ppdvspupildiameter()
+    x = 20:2:45
+    y = 3.0:.05:4
+
+    RGB = true
+    Plots.plot(Plots.contour(x, y, (x, y) -> pixelredundancy((50, 35), 18mm, (10mm, 6mm), y * mm, x, RGB=RGB), fill=true, xlabel="pixels per degree", ylabel="pupil diameter", legendtitle="pixel redundancy", title="$(label(RGB)) lenslets"))
+end
+export redundancy_ppdvspupildiameter
+
+"""computes lenslet display size to match the design constraints"""
+function lensletdisplaysize(fov, eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, pupildiameter::Unitful.Length, ppd; RGB=true)
+    lensprops = defaultclusterproperties()
+    angles = lensletangles(eyerelief, eyebox, pupildiameter, ppd, RGB=RGB)
+    return @. angles * ppd * lensprops.pixelpitch
+end
+export lensletdisplaysize
+
+testlensletdisplaysize() = lensletdisplaysize((55, 35), 18mm, (10mm, 6mm), 4mm, 30, RGB=true)
+export testlensletdisplaysize
+
+"""generates a contour plot of lenslet display size as a function of ppd and pupil diameter"""
+function displaysize_ppdvspupildiameter()
+    x = 20:2:45
+    y = 3.0:.05:4
+    RGB = true
+
+    Plots.plot(Plots.contour(x, y, (x, y) -> maximum(ustrip.(μm, lensletdisplaysize((50, 35), 18mm, (10mm, 6mm), y * mm, x, RGB=RGB))), fill=true, xlabel="pixels per degree", ylabel="pupil diameter", legendtitle="display size μm", title="$(label(RGB)) lenslets"))
+end
+export displaysize_ppdvspupildiameter
+
+""" Compute high level properties of a lenslet display system. Uses basic geometric and optical constraints to compute approximate values which should be within 10% or so of a real physical system.
+
+`eyerelief` is in mm: `18mm` not `18`.
+
+`eyebox` is a 2 tuple of lengths, also in mm, that specifies the size of the rectangular eyebox, assumed to lie on vertex of the cornea when the eye is looking directly forward.
+
+`fov` is the field of view of the display as seen by the eye, in degrees (no unit type necessary here).
+
+`pupildiameter` is the diameter of the eye pupil in mm.
+
+`mtf` is the MTF response of the system at the angular frequency `cyclesperdegree`.
+
+`pixelsperdegree` is the number of pixels per degree on the display, as seen by the eye through the lenslets. This is not a specification of the optical resolution of the system (use `mtf` and `cyclesperdegree` for that). It is a physical characteristic of the display.
+
+Example:
+```
+julia> systemproperties(18mm,(10mm,9mm),(55°,45°),4.0mm,.2,11,30)
+(lenslet_diameter = 0.7999999999999999 mm, diffraction_limit = 26.344592482933276, fnumber = 1.9337325040249893, focal_length = 1.5469860032199914 mm, display_size = (18.740413819862866 mm, 14.911688245431423 mm), lenslet_display_size = (261.4914368916943 μm, 358.6281908905529 μm), total_silicon_area = 52.1360390897504 mm^2, number_lenslets = 555.9505147668694, pixel_redundancy = 28.895838544429424 °^-2, eyebox_angles = (29.054604099077146, 26.56505117707799), lenslet_fov = (9.684868033025715, 13.282525588538995), subdivisions = (3, 2))
+```
+"""
+function systemproperties(eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, fov, pupildiameter::Unitful.Length, mtf, cyclesperdegree,pixelsperdegree; minfnumber=2.0,RGB=true,λ=530nm,pixelpitch=.9μm)
+    diameter = diameter_for_cycles_deg(mtf, cyclesperdegree, λ)
+    clusterdata = choosecluster(pupildiameter, diameter)
+    difflimit = diffractionlimit(λ, clusterdata.lensletdiameter)
+    numlenses = numberoflenslets(fov, eyerelief, clusterdata.lensletdiameter)
+    redundancy = pixelredundancy(fov, eyerelief, eyebox, pupildiameter, difflimit, RGB=RGB)
+    subdivisions = anglesubdivisions(pupildiameter, λ, mtf, cyclesperdegree, RGB=RGB)
+    eyebox_angles = eyeboxangles(eyebox,eyerelief)
+    angles = lensletangles(eyerelief, eyebox, pupildiameter, difflimit, clusterproperties=(mtf = mtf, minfnumber = minfnumber, cyclesperdegree = cyclesperdegree, λ = λ, pixelpitch = pixelpitch))
+    dispsize = lensletdisplaysize(angles, eyerelief, eyebox, pupildiameter, pixelsperdegree, RGB=RGB)
+    focal_length = focallength(pixelpitch,1/pixelsperdegree)
+    fnumber = focal_length/clusterdata.lensletdiameter
+    siliconarea = uconvert(mm^2,numlenses  * dispsize[1]*dispsize[2])
+    fulldisplaysize = sizeofdisplay(fov,eyerelief)
+
+    return (lenslet_diameter = clusterdata.lensletdiameter, diffraction_limit = difflimit, fnumber = fnumber, focal_length = focal_length, display_size = fulldisplaysize, lenslet_display_size = dispsize, total_silicon_area = siliconarea, number_lenslets = numlenses, pixel_redundancy = redundancy, eyebox_angles = eyebox_angles, lenslet_fov = angles, subdivisions = subdivisions)
+end
+export systemproperties
+
+"""prints system properties nicely"""
+function printsystemproperties(eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, fov, pupildiameter::Unitful.Length, mtf, cyclesperdegree,pixelsperdegree; minfnumber=2.0,RGB=true,λ=530nm,pixelpitch=.9μm) 
+    println("eye relief = $eyerelief")
+    println("eye box = $eyebox")
+    println("fov = $(fov)°")
+    println("pupil diameter = $pupildiameter")
+    println("mtf = $mtf @ $cyclesperdegree cycles/degree")
+    for (key,value) in pairs(systemproperties(eyerelief, eyebox, fov, pupildiameter, mtf, cyclesperdegree,pixelsperdegree, minfnumber = minfnumber,RGB=RGB,λ=λ,pixelpitch=pixelpitch))
+        if key == :diffraction_limit
+            println("$key = $value cycles/°")
+        else
+            println("$key = $value")
+        end
+    end
+end
+export printsystemproperties
+
+
diff --git a/src/RepeatingStructures/Multilens/DisplayGeneration.jl b/src/RepeatingStructures/Multilens/DisplayGeneration.jl
new file mode 100644
index 000000000..9e616294e
--- /dev/null
+++ b/src/RepeatingStructures/Multilens/DisplayGeneration.jl
@@ -0,0 +1,203 @@
+# MIT license
+# Copyright (c) Microsoft Corporation. All rights reserved.
+# See LICENSE in the project root for full license information.
+
+
+#display lenslet systems: emitter, paraxial lenslet, paraxial eye lens, retina detector.
+#start with eyebox plane and hexagonal tiling. Project onto display surface, generate lenslets automatically.
+
+using OpticSim.Geometry:Transform,world2local
+using OpticSim:plane_from_points,surfaceintersection,closestintersection,Ray,Plane,ConvexPolygon,Sphere 
+using OpticSim
+using OpticSim.Repeat:tilevertices,HexBasis1,tilesinside
+using Unitful:upreferred
+using Unitful.DefaultSymbols:°
+
+# using OpticSim.Vis:drawcells
+
+"""compute the mean of the columns of `a`. If `a` is an `SMatrix` this is very fast and will not allocate."""
+columncentroid(a::AbstractMatrix) = sum(eachcol(a))/size(a)[2] #works except for the case of zero dimensional matrix.
+export columncentroid
+
+function project(point::AbstractVector{T},projectionvector::AbstractVector{T},surface::OpticSim.Surface{T}) where{T}
+    ray = Ray(point, projectionvector)
+    intsct = surfaceintersection(surface, ray)
+    pointintsct = closestintersection(intsct,false)
+    if pointintsct === nothing #point didn't project to a point on the surface
+        return nothing
+    else
+        return OpticSim.point(pointintsct)
+    end
+end
+export project
+
+"""project the vertices of a polygon represented by `vertices` onto `surface` using the point as the origin and `projectionvector` as the projection direction. Return nothing if any of the projected points do not intersect the surface. The projected vertices are not guaranteed to be coplanar."""
+function project(vertices::AbstractMatrix{T}, projectionvector::AbstractVector{T}, surface::OpticSim.Surface{T}) where {T}
+    result = similar(vertices)
+
+    for i in 1:size(vertices)[2]
+        origin = vertices[:, i]
+        pt = project(origin,projectionvector,surface)
+
+        if pt === nothing #one of the polygon points didn't project to a point on the surface so reject the polygon
+            return nothing
+        else
+            result[:, i] = pt
+        end
+    end
+    return result
+end
+
+
+"""Finds the best fit plane to `vertices` then projects `vertices` onto this plane by transforming from the global to the local coordinate frame. The projected points are represented in the local coordinate frame of the plane."""
+function projectonbestfitplane(vertices::AbstractMatrix{T}) where{T}
+    @assert size(vertices)[1] == 3 "projection only works for 3D points"
+
+    center, _, localrotation  = plane_from_points(vertices) 
+    toworld = Transform(localrotation,center) #compute local to world transformation
+    tolocal = world2local(toworld)
+    numpts = size(vertices)[2]
+    result = tolocal * SMatrix{3,numpts}(vertices) 
+
+    return result,toworld,tolocal #project onto best fit plane by setting z coordinate to zero.
+end
+export projectonbestfitplane
+
+function testproject()
+    normal = SVector(0.0, 0, 1)
+    hex = HexBasis1()
+    verts =vcat(tilevertices((0, 0), hex), [0.0 0 0 0 0 0])
+    surf = Plane(normal, SVector(0.0, 0, 10))
+
+    project(verts, normal, surf)
+end
+export testproject
+
+function testprojectonplane()
+    verts = tilevertices(HexBasis1())
+    verts = vcat(verts,[0 for _ in 1:6]')
+    projectonbestfitplane(SMatrix{size(verts)...}(verts))
+end
+export testprojectonplane
+
+"""projects convex polygon, represented by `vertices`, onto `surface` along vector `normal`. Assumes original polygon is convex and that the projection will be convex. No guarantee that this will be true but for smoothly curved surfaces that are not varying too quickly relative to the size of the polygon it should be true."""
+function planarpoly(projectedpoints::AbstractMatrix{T}) where{T}
+    planarpoints,toworld,_ = projectonbestfitplane(projectedpoints)
+    numpts = size(planarpoints)[2]
+    temp = SMatrix{2,numpts}(planarpoints[1:2,:])
+    vecofpts = collect(reinterpret(reshape,SVector{2,T},temp))
+    #this code is inefficient because of different data structures used by convex polygon and most other code for representing lists of vertices.
+     return ConvexPolygon(toworld,vecofpts)
+end
+
+spherepoint(radius,θ,ϕ) = radius .* SVector(cos(θ)sin(ϕ),sin(θ),cos(θ)cos(ϕ))
+export spherepoint
+
+"""Computes points on the edges of the spherical rectangle defined by the range of θ,ϕ. This is used to determine lattice boundaries on the eyebox surface."""
+function spherepoints(radius, θmin,θmax,ϕmin,ϕmax)
+    θedges =  [spherepoint(radius,ϕ,θ) for θ in θmin:.5:θmax, ϕ in (ϕmin,ϕmax)]
+    ϕedges =  [spherepoint(radius,ϕ,θ) for ϕ in ϕmin:.5:ϕmax, θ in (θmin,θmax)]
+    allpoints = vcat(reshape(θedges,reduce(*,size(θedges))),reshape(ϕedges,reduce(*,size(ϕedges))))
+    
+    reshape(reinterpret(Float64,allpoints),3,length(allpoints)) #return points as 3xn matrix with points as columns
+end
+export spherepoints
+
+"""given a total fov in θ  and ϕ compute sample points on the edges of the spherical rectangle."""
+spherepoints(radius,θ,ϕ) = spherepoints(radius,-θ/2,θ/2,-ϕ/2,ϕ/2)
+
+function testprojection()
+    pts = spherepoints(1.0,-.2,-.2,1.0,1.1)
+    surf = Plane(0.0,0.0,-1.0,0.0,0.0,0.0)
+    dir = [0.0,0.0,-1.0]
+    project(pts,dir,surf)
+end
+export testprojection
+
+function bounds(pts::AbstractMatrix{T}) where{T} 
+    return [extrema(row) for row in eachrow(pts)]
+end
+export bounds
+
+function eyeboxbounds(eyebox::OpticSim.Plane,dir::AbstractVector, radius,fovθ,fovϕ) 
+    pts = spherepoints(radius,fovθ,fovϕ)
+    projectedpts = project(pts,dir,eyebox)
+    return bounds(projectedpts)
+end
+export eyeboxbounds
+
+function boxtiles(bbox,lattice)
+    tiles = tilesinside(bbox[1][1],bbox[2][1],bbox[1][2],bbox[2][2],lattice)
+end
+export boxtiles
+
+eyeboxtiles(eyebox,dir,radius,fovθ,fovϕ,lattice) = boxtiles(eyeboxbounds(eyebox,dir,radius,fovθ,fovϕ),lattice)
+export eyeboxtiles
+
+function spherepolygon(vertices,projectiondirection,sph::Sphere)
+    @assert typeof(vertices) <: SMatrix
+    projectedpts = project(vertices,projectiondirection,sph)
+    return planarpoly(projectedpts) #polygon on best fit plane to projectedpts
+end
+
+#TODO need to figure out what to use as the normal (eventually this will need to take into account the part of the eyebox the lenslet should cover) 
+#TODO write test for planarpoly and generation of Paraxial lens system using planar hexagon as lenslet.
+function spherepolygons(eyebox::Plane{T,N},dir,radius,fovθ,fovϕ,lattice) where{T,N}
+    #if fovθ, fovϕ are in degrees convert to radians. If they are unitless then the assumption is that the represent radians
+    eyeboxz = eyebox.pointonplane[3]
+    θ = upreferred(fovθ) #converts to radians if in degrees
+    ϕ = upreferred(fovϕ) #converts to radians if in degrees
+    tiles = eyeboxtiles(eyebox,dir,radius,θ,ϕ,lattice)
+    shapes = Vector{ConvexPolygon{6,T}}(undef,0)
+    coordinates = Vector{Tuple{Int64,Int64}}(undef,0)
+    for coords in eachcol(tiles)
+        twodverts = tilevertices(coords,lattice)
+        numverts = size(twodverts)[2]
+        temp = MMatrix{1,numverts,T}(undef)
+        fill!(temp,eyeboxz)
+        zerorow = SMatrix{1,numverts,T}(temp)
+        vertices = SMatrix{N,numverts,T}(vcat(twodverts,zerorow))
+        @assert typeof(vertices) <: SMatrix
+        push!(shapes,spherepolygon(vertices,-dir,Sphere(radius)))
+        push!(coordinates,Tuple{T,T}(coords))
+    end
+    shapes,coordinates
+end
+
+function spherelenslets(eyebox::Plane{T,N},dir,radius,fovθ,fovϕ,lattice) where{T,N}
+    lenspolys,tilecoords = spherepolygons(eyebox,dir,radius,fovθ,fovϕ,lattice)
+    result = Vector{ParaxialLens{T}}(undef,length(lenspolys))
+    empty!(result)
+    for poly in lenspolys
+        lenslet = ParaxialLensConvexPoly(2.0,poly,SVector{2,T}(T(0),T(0)))
+        push!(result,lenslet)
+    end
+    return result,tilecoords
+end
+export spherelenslets
+
+function testspherelenslets()
+    spherelenslets(Plane(0.0,0.0,1.0,0.0,0.0,12.0),[0.0,0.0,-1.0],30.0,55°,45°,HexBasis1())
+end
+export testspherelenslets
+
+function testspherepolygons()
+    eyebox = Plane(0.0,0.0,1.0,0.0,0.0,12.0)
+    dir = [0.0,0.0,-1.0]
+    radius = 30.0
+    fovθ = 55°
+    fovϕ = 45°
+    lattice = HexBasis1()
+
+    return spherepolygons(eyebox,dir,radius,fovθ,fovϕ,lattice)
+end
+export testspherepolygons
+
+
+# function testeyeboxtiles()
+#     tiles = eyeboxtiles(Plane(0.0,0.0,1.0,0.0,0.0,12.0),[0.0,0.0,-1.0],30,deg2rad(55),deg2rad(45),HexBasis1())
+#     drawcells(HexBasis1(),10,tiles)
+#     tiles = eyeboxtiles(Plane(0.0,0.0,1.0,0.0,0.0,12.0),[0.0,0.0,-1.0],30,deg2rad(25),deg2rad(45),HexBasis1())
+#     drawcells(HexBasis1(),10,tiles)
+# end
+# export testeyeboxtiles
\ No newline at end of file
diff --git a/src/RepeatingStructures/Multilens/HexClusters.jl b/src/RepeatingStructures/Multilens/HexClusters.jl
new file mode 100644
index 000000000..dcaec7aca
--- /dev/null
+++ b/src/RepeatingStructures/Multilens/HexClusters.jl
@@ -0,0 +1,68 @@
+# MIT license
+# Copyright (c) Microsoft Corporation. All rights reserved.
+# See LICENSE in the project root for full license information.
+
+#Various sizes of hexagonal clusters that are useful in lenslet design. To see what they look like use the functions in HexTilings to make lattices with color properties and then draw them with Vis.draw
+#example: Vis.draw(hex3RGB())
+
+function hex3()
+    clusterelements = SVector((0, 0), (-1, 0), (-1, 1))
+    eltlattice = Repeat.HexBasis1()
+    clusterbasis = Repeat.LatticeBasis((-1, 2), (2, -1))
+    return Repeat.LatticeCluster(clusterbasis, eltlattice, clusterelements)
+
+end
+export hex3
+
+function hex4()
+    clusterelements = SVector((-1, 0), (-1, 1), (0, 0), (0, 1))
+    eltlattice = Repeat.HexBasis1()
+    clusterbasis = Repeat.LatticeBasis((2, 0), (0, 2))
+    return Repeat.LatticeCluster(clusterbasis, eltlattice, clusterelements)
+end
+export hex4
+
+function hex7()
+    clusterelements = SVector((1, -1), (0, -1), (-1, 0), (0, 0), (-1, 1), (0, 1), (1, 0), )
+    eltlattice = Repeat.HexBasis1()
+    clusterbasis = Repeat.LatticeBasis((3, 0), (0, 3))
+    clusterbasis = Repeat.LatticeBasis((2, 2), (-3, 2))
+    return Repeat.LatticeCluster(clusterbasis, eltlattice, clusterelements)
+end
+export hex7
+
+
+function hex9()
+    clusterelements = SVector((-1, -1), (-1, 0), (-2, 1), (0, -1), (0, 0), (-1, 1), (1, -1), (1, 0), (0, 1))
+    eltlattice = Repeat.HexBasis1()
+    clusterbasis = Repeat.LatticeBasis((3, 0), (0, 3))
+    return Repeat.LatticeCluster(clusterbasis, eltlattice, clusterelements)
+end
+export hex9
+
+function hex12()
+    clusterelements = SVector(
+        (-1, 1),(0, 1),
+        (-1, 0),(0, 0),(1, 0),
+        (-1, -1),(0, -1),(1, -1),(2, -1),
+        (0, -2),(1, -2),(2, -2)
+    )
+
+    eltlattice = Repeat.HexBasis3()
+    clusterbasis = Repeat.LatticeBasis((2, 2), (-3, 2))
+    return Repeat.LatticeCluster(clusterbasis, eltlattice, clusterelements)
+end
+export hex12
+
+function hex19()
+    clusterelements = SVector(
+        (0, 0), # ring 0
+        (-1, 0),(0, -1),(1, -1),(1, 0),(0, 1),(-1, 1), # ring1
+        (-2, 0),(-1, -1),(0, -2),(1, -2),(2, -2),(2, -1),(2, 0),(1, 1),(0, 2),(-1, 2),(-2, 2),(-2, 1) # ring2
+        )
+    
+    eltlattice = Repeat.HexBasis3()
+    clusterbasis = Repeat.LatticeBasis((5, 0), (-2, 4))
+    return Repeat.LatticeCluster(clusterbasis, eltlattice, clusterelements)
+end
+export hex19
diff --git a/src/RepeatingStructures/Multilens/HexTilings.jl b/src/RepeatingStructures/Multilens/HexTilings.jl
new file mode 100644
index 000000000..ac3750038
--- /dev/null
+++ b/src/RepeatingStructures/Multilens/HexTilings.jl
@@ -0,0 +1,131 @@
+# MIT license
+# Copyright (c) Microsoft Corporation. All rights reserved.
+# See LICENSE in the project root for full license information.
+
+# All of the objects can be displayed with Vis.draw. Example:
+# Vis.draw(hex4RGB())
+
+function clustercolors(lattice)
+    elements = Repeat.clusterelements(lattice)
+    colors = Vector{String}(undef, length(elements))
+
+    for (i, element) in pairs(elements)
+        colors[i] = pointcolor(element, lattice)
+    end
+    return colors
+end
+export clustercolors
+
+function hex3RGB()
+    lattice = hex3()
+    colors = clustercolors(lattice)
+    names = ["R","G","B"]
+    properties =  DataFrames.DataFrame(Color=colors, Name=names)
+    return Repeat.ClusterWithProperties(hex3(), properties)
+end
+export hex3RGB
+
+function hex4RGB()
+    lattice = hex4()
+    colors = clustercolors(lattice)
+    names = ["R","G","B","W"]
+    properties =  DataFrames.DataFrame(Color=colors, Name=names)
+    return Repeat.ClusterWithProperties(hex4(), properties)
+end
+export hex4RGB
+
+function hex7RGB()
+    lattice = hex7()
+    colors = clustercolors(lattice)
+    names = ["R1","G1","B1","W","R2","G2","B2"]
+    properties =  DataFrames.DataFrame(Color=colors, Name=names)
+    return Repeat.ClusterWithProperties(hex7(), properties)
+end
+export hex7RGB
+
+function hex9RGB()
+    lattice = hex9()
+    colors = clustercolors(lattice)
+    names = ["R-1","G-1","B-1","R0","G0","B0","R1","G1","B1"]
+    properties =  DataFrames.DataFrame(Color=colors, Name=names)
+    return Repeat.ClusterWithProperties(hex9(), properties)
+end
+export hex9RGB
+
+function hex12RGB()
+   
+    lattice = hex12()
+
+    colors = clustercolors(lattice)
+
+    names = [
+        "G3","B0",
+        "B2","R1","G2",
+        "R0","G1","B1","R3",
+        "B3","R2","G0"
+    ]
+    properties =  DataFrames.DataFrame(Color=colors, Name=names)
+    return Repeat.ClusterWithProperties(hex12(), properties)
+end
+export hex12RGB
+
+function hex19RGB()
+    lattice = hex19()
+
+    colors = clustercolors(lattice)
+    
+    names = [
+        "W",
+        "B0","G1","R2","B3","G4","R5",
+        "G0","B1","R1","G2","B2","R3","G3","B4","R4","G5","B5","R0"
+    ]
+    properties =  DataFrames.DataFrame(Color=colors, Name=names)
+    # properties =  DataFrames.DataFrame(Color = colors)
+    return Repeat.ClusterWithProperties(lattice, properties)
+end
+export hex19RGB
+
+function hex19fov1()
+    lattice = hex19()
+
+    offset = 1
+    colornames = ["R","G","B"]
+    colors = ["red","green","blue"]
+    fov1 = [8,19,18,17,7,2]
+    fov2 = [16,15,14,13,5,6] 
+    fov3 = [12,11,10,9,3,4]
+
+    allfov = vcat([0],fov1,fov2,fov3)
+    names = [x -> string("fov",x) for x in 1:6]
+    
+    allnames = Vector{String}(undef,19)
+    allcolors = Vector{String}(undef,19)
+
+    for (index,fov) in pairs((fov1,fov2,fov3))
+        for i in eachindex(fov)
+            println("$index $i")
+            latticeindex = (index-1)*length(fov1) + i + offset
+            allcolors[allfov[latticeindex]] = colors[index]
+            allnames[allfov[latticeindex]] = string("fov",string(index))
+        end
+    end
+    allcolors[offset] = "white"
+    allnames[offset] = "W"
+    properties =  DataFrames.DataFrame(Color=allcolors, Name=allnames)
+    # properties =  DataFrames.DataFrame(Color = colors)
+    return Repeat.ClusterWithProperties(lattice, properties)
+end
+export hex19fov1
+
+
+function hex19fov2()
+    lattice = hex19()
+
+    colors = clustercolors(lattice)
+    names = [string("fov",x÷3) for x in 1:18]
+    names = vcat("W",names)
+    properties =  DataFrames.DataFrame(Color=colors, Name=names)
+    # properties =  DataFrames.DataFrame(Color = colors)
+    return Repeat.ClusterWithProperties(lattice, properties)
+end
+export hex19fov2
\ No newline at end of file
diff --git a/src/RepeatingStructures/Multilens/Lenslets.jl b/src/RepeatingStructures/Multilens/Lenslets.jl
new file mode 100644
index 000000000..0f95baf78
--- /dev/null
+++ b/src/RepeatingStructures/Multilens/Lenslets.jl
@@ -0,0 +1,23 @@
+# MIT license
+# Copyright (c) Microsoft Corporation. All rights reserved.
+# See LICENSE in the project root for full license information.
+module Lenslets
+using LinearAlgebra
+import Unitful
+using Unitful:uconvert,ustrip
+using Unitful.DefaultSymbols:mm,μm,nm #Unitful.DefaultSymbols exports a variable T which can cause major confusion with type declarations since T is a common parametric type symbol. Import only what is needed.
+using Colors
+using StaticArrays
+import DataFrames
+import ..Repeat
+import SpecialFunctions
+import Plots
+import ...OpticSim
+
+include("HexClusters.jl")
+include("HexTilings.jl")
+include("Analysis.jl")
+include("DisplayGeneration.jl")
+
+end # module
+export Lenslets
\ No newline at end of file
diff --git a/src/RepeatingStructures/RectangularLattice.jl b/src/RepeatingStructures/RectangularLattice.jl
index fedec98ed..3a82d35ae 100644
--- a/src/RepeatingStructures/RectangularLattice.jl
+++ b/src/RepeatingStructures/RectangularLattice.jl
@@ -13,7 +13,7 @@ basismatrix(::RectangularBasis{2,T}) where{T} = SMatrix{2,2,T}(T(1),T(0),T(0),T(
 # SVector{2,SVector{2,T}}(hexe₁(T),hexe₂(T))
 
 """Returns the vertices of the unit tile polygon for the basis"""
-function tilevertices(::RectangularBasis{2,T}) where{T}
+function vertices(::RectangularBasis{2,T}) where{T}
 	return SMatrix{2,4,T}(
 			-.5, -.5,
 			-.5, .5,
diff --git a/src/RepeatingStructures/Repeat.jl b/src/RepeatingStructures/Repeat.jl
index cb6040ad7..45c5a0134 100644
--- a/src/RepeatingStructures/Repeat.jl
+++ b/src/RepeatingStructures/Repeat.jl
@@ -10,12 +10,15 @@ using Base: offset_if_vec
 using StaticArrays:SVector,MVector,SMatrix,MMatrix
 using DataFrames:DataFrame
 using Colors
+import LazySets
 
 include("Lattice.jl")
 include("HexagonalLattice.jl")
 include("RectangularLattice.jl")
 include("Array.jl")
 include("Cluster.jl")
+include("Multilens/Lenslets.jl")
+
 
 end #module
 export Repeat
\ No newline at end of file
diff --git a/src/Vis/VisRepeatingStructures.jl b/src/Vis/VisRepeatingStructures.jl
index d5d42f37e..6e5921981 100644
--- a/src/Vis/VisRepeatingStructures.jl
+++ b/src/Vis/VisRepeatingStructures.jl
@@ -10,11 +10,11 @@
 
 # lattice visualizations are drawn with Luxor because it is easier to do 2D drawings with Luxor than with Makie.
 
-function draw(tilebasis::AbstractBasis,tilesize,i,j,color,name)
+function draw(tilebasis::AbstractBasis, tilesize, i, j, color, name)
     vertices = Repeat.tilevertices(tilebasis)
-    tile = tilesize*[Luxor.Point(vertices[:,i]...) for i in 1:size(vertices)[2]]
-    pt = tilesize*tilebasis[i,j]
-    offset = Luxor.Point(pt[1],-pt[2]) #flip y so indices show up correctly
+    tile = tilesize * [Luxor.Point(vertices[:,i]...) for i in 1:size(vertices)[2]]
+    pt = tilesize * tilebasis[i,j]
+    offset = Luxor.Point(pt[1], -pt[2]) # flip y so indices show up correctly
     
     Luxor.translate(offset)
     
@@ -25,31 +25,31 @@ function draw(tilebasis::AbstractBasis,tilesize,i,j,color,name)
     Luxor.poly(tile, :stroke, close=true)
     Luxor.sethue("black")
 
-    #scale and offset text so coordinates are readable
-    Luxor.fontsize(tilesize/3)
-    Luxor.text("$i, $j",Luxor.Point(-tilesize/3,tilesize/2.5))
+    # scale and offset text so coordinates are readable
+    Luxor.fontsize(tilesize / 3)
+    Luxor.text("$i, $j", Luxor.Point(-tilesize / 3, tilesize / 2.5))
     if name !== nothing
-        Luxor.text(name,Luxor.Point(-tilesize/4,-tilesize/4))
+        Luxor.text(name, Luxor.Point(-tilesize / 4, -tilesize / 4))
     end
 
     Luxor.translate(-offset)
 end
 
-"""Draws a list of hexagonal cells, represented by their lattice coordinates"""
-function drawcells(tilebasis::AbstractBasis, tilesize,cells; color::Union{AbstractArray,Nothing} = nothing, name::Union{AbstractArray{String},Nothing} = nothing, format=:png, resolution=(500,500))
+"""Draws a list of hexagonal cells, represented by their lattice coordinates, which are represented as a 2D matrix, with each column being one lattice coordinate."""
+function drawcells(tilebasis::AbstractBasis, tilesize, cells::AbstractMatrix; color::Union{AbstractArray,Nothing}=nothing, name::Union{AbstractArray{String},Nothing}=nothing, format=:png, resolution=(500, 500))
 
     numcells = size(cells)[2]
     Luxor.Drawing(resolution[1], resolution[2], format)
     Luxor.origin()
     Luxor.background(Colors.RGBA(0, 1, 1, 0.0))
     if color === nothing
-        color = Colors.distinguishable_colors(numcells,lchoices = 30:250) #type unstable but not performance critical code
+        color = Colors.distinguishable_colors(numcells, lchoices=30:250) # type unstable but not performance critical code
     end
 
     for i in 1:numcells
         cell = cells[:,i]
         cellname = name === nothing ? nothing : name[i]
-        draw(tilebasis,tilesize,cell[1],cell[2],color[i],cellname)
+        draw(tilebasis, tilesize, cell[1], cell[2], color[i], cellname)
     end
  
     Luxor.finish()
@@ -63,19 +63,19 @@ function drawcells(tilebasis::AbstractBasis, tilesize,cells; color::Union{Abstra
 end
 
 """ draw the LatticeCluster offset to (0,0) """
-draw(clstr::LatticeCluster,tilesize = 50.0, cells = hcat([0,0])) = drawcells(clstr.elementbasis,tilesize,clustercoordinates(clstr,0,0))
+draw(clstr::LatticeCluster,tilesize=50.0, cells=hcat([0,0])) = drawcells(clstr.elementbasis, tilesize, clustercoordinates(clstr, 0, 0))
 
 """ draw the ClusterWithProperties at coordinates specified by cells """
-function draw(clstr::Repeat.ClusterWithProperties, cells = hcat([0,0]), scale = 50.0) #this is how you have to make a 2x1 matrix. One would expect [0;0] to work but it doesn't.
+function draw(clstr::Repeat.ClusterWithProperties, cells=hcat([0,0]), scale=50.0) # this is how you have to make a 2x1 matrix. One would expect [0;0] to work but it doesn't.
     dims = size(cells)
     @assert dims[1] == 2 "dims[1] must be 2 instead was $(dims[1])"
     clstrsize = clustersize(clstr)
-    points = Matrix(undef,dims[1],dims[2]*clstrsize)
+    points = Matrix(undef, dims[1], dims[2] * clstrsize)
     for i in 1:dims[2]
-        points[:,(i-1)*clstrsize+1:i*clstrsize] = clustercoordinates(clstr,cells[1,i],cells[2,i])
+        points[:,(i - 1) * clstrsize + 1:i * clstrsize] = clustercoordinates(clstr, cells[1,i], cells[2,i])
     end
-    props = repeat(Repeat.properties(clstr),dims[2])
-    drawcells(elementbasis(cluster(clstr)),scale,points,color = props[:,:Color],name = props[:,:Name])
+    props = repeat(Repeat.properties(clstr), dims[2])
+    drawcells(elementbasis(cluster(clstr)), scale, points, color=props[:,:Color], name=props[:,:Name])
 end
     
 
diff --git a/src/Vis/Visualization.jl b/src/Vis/Visualization.jl
index a6af04d8f..562f87ee5 100644
--- a/src/Vis/Visualization.jl
+++ b/src/Vis/Visualization.jl
@@ -927,7 +927,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 * Transform(r, t))
+            sys = HumanEye.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/src/utilities.jl b/src/utilities.jl
index 95878069f..381c7ac41 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -8,21 +8,31 @@
 
 Estimate the best fitting plane for a set of points in 3D.
 A more efficient version of plane_from_points.
+`D` is the dimension of the plane.
+`N` is the number of points to fit.
+`T` is the type of the number used to represent points.
 """
-function plane_from_points(points::SMatrix{3, N, Float64} where {N}) 
+function plane_from_points(points::SMatrix{D, N, P}) where {D,N,P<:Real}
     center = Statistics.mean(points,dims=2)
 
     u, _, _ = svd(points .- center)
-    normal = u[:,3]             # singular vectors in decending order
+    normal = u[:,3]             # The two largest singular vectors lie in the plane that is the best fit to the points, i.e., that accounts for the largest fraction of variance in the set of points. The smallest singular vector is perpendicular to this plane.
 
     # make sure the normal is pointing consistently to positive Z direction 
-    if (dot(normal, unitZ3()) < 0.0)
-        normal = normal * -1.0
+    if dot(normal, unitZ3()) < P(0.0)
+        normal = normal * P(-1.0)
     end
 
-    return SVector(center), SVector(normal)     # convert from SMatrix to SVector
+    return SVector(center), SVector(normal), u     # convert from SMatrix to SVector, return u matrix to use as local coordinate frame for the plane.
 end
 
+function plane_from_points(points::AbstractMatrix{P}) where{P<:Real}
+    dims = size(points)
+    temp = MMatrix{dims[1],dims[2],P}(points)
+    return plane_from_points(SMatrix(temp))
+end
+
+
 """only returns real roots"""
 function quadraticroots(a::T, b::T, c::T) where {T<:Real}
     temp = b^2 - 4 * a * c
diff --git a/test/testsets/Repeat.jl b/test/testsets/Repeat.jl
index b6f5adcd6..27f4224d6 100644
--- a/test/testsets/Repeat.jl
+++ b/test/testsets/Repeat.jl
@@ -14,6 +14,13 @@
         return Repeat.ClusterWithProperties(lattice,properties)
     end
 
+    #spherepoint tests
+    @test isapprox(Repeat.Lenslets.spherepoint(1,π/2,0.0), [0.0,1.0,0.0])
+    @test isapprox(Repeat.Lenslets.spherepoint(1,0.0,π/2), [1.0,0.0,0.0])
+    @test isapprox(Repeat.Lenslets.spherepoint(1,0,0.0), [0.0,0.0,1.0])
+    @test isapprox(Repeat.Lenslets.spherepoint(1,0.0,π/4), [sqrt(2)/2,0.0,sqrt(2)/2])
+
+
     """ Create a LatticeCluser with three elements at (0,0),(-1,0),(-1,1) coordinates in the HexBasis1 lattice"""
     function hex3cluster()
         clusterelts = SVector((0,0),(-1,0),(-1,1))
@@ -28,5 +35,22 @@
         return Repeat.clusterbasis(a)
     end
 
-    @test basistest(hex3cluster()) == basistest(hex3RGB()) 
+    @test basistest(hex3cluster()) == basistest(hex3RGB())  
+
+    #LatticeCluster testset
+    cluster = Repeat.Lenslets.hex9()
+
+    for iter in 1:100
+        (i,j) = rand.((1:1000,1:1000))
+        coords,tileindex = Repeat.cluster_coordinates_from_tile_coordinates(cluster,i,j)
+        reconstructed = Repeat.tilecoordinates(cluster,coords...,tileindex)
+        @test all((i,j) .== reconstructed)
+    end
+ 
+    #verify that the 0,0 cluster is correct
+    for (index,element) in pairs(Repeat.clusterelements(cluster))
+        coords, tileindex = Repeat.cluster_coordinates_from_tile_coordinates(cluster, element...)
+        @test all(coords .== 0)
+        @test tileindex == index
+    end
 end
\ No newline at end of file