diff --git a/src/OpticSim.jl b/src/OpticSim.jl
index 8f84670d8..5b6d2a1e1 100644
--- a/src/OpticSim.jl
+++ b/src/OpticSim.jl
@@ -6,6 +6,8 @@ module OpticSim
 
 import Unitful
 using LinearAlgebra: eigen, svd, I, qr, dot, cross, norm, det, normalize, inv
+import LinearAlgebra
+
 using StaticArrays
 using DataFrames: DataFrame
 using Images
diff --git a/src/Optical/ParaxialAnalysis/HeadEyeModel/Arrangement.jl b/src/Optical/ParaxialAnalysis/HeadEyeModel/Arrangement.jl
index b1f7aeb04..bcfa92c2a 100644
--- a/src/Optical/ParaxialAnalysis/HeadEyeModel/Arrangement.jl
+++ b/src/Optical/ParaxialAnalysis/HeadEyeModel/Arrangement.jl
@@ -101,7 +101,7 @@ function build_paraxial_lens(shape::Shape{3, T}; local_center_point = SVector(0.
     # estimate a best fitting plane (least-squares wise)
     # switching to SMatrix format - should consider using SMatrix in the basic shapes instead of a vector of points
     pts_mat = reshape(SVector(reduce(vcat, pts)...), StaticArrays.Size(length(pts[1]),length(pts)))
-    fitted_center, fitted_normal = HeadEye.plane_from_points2(pts_mat)    
+    fitted_center, fitted_normal = OpticSim.plane_from_points(pts_mat)    
 
     # create a local frame for the fitted plane
     local_frame = Transform(fitted_center, fitted_normal)
diff --git a/src/Optical/ParaxialAnalysis/HeadEyeModel/Misc.jl b/src/Optical/ParaxialAnalysis/HeadEyeModel/Misc.jl
index 77823fe30..a4c80395e 100644
--- a/src/Optical/ParaxialAnalysis/HeadEyeModel/Misc.jl
+++ b/src/Optical/ParaxialAnalysis/HeadEyeModel/Misc.jl
@@ -6,9 +6,11 @@
 """
     plane_from_points(points::Vector{SVector{3, Float64}}) ->  centroid, normal
 
-Estimate the best fitting plane for a set of points in 3D    
+Estimate the best fitting plane for a set of points in 3D.
+    
+See utilities.jl for a simpler and possibly more numerically stable version.  
 """
-function plane_from_points(points::Vector{SVector{3, Float64}}) 
+function plane_from_points2(points::Vector{SVector{3, Float64}}) 
     if length(points) < 3 
         return nothing      # At least three points required
     end
@@ -55,26 +57,6 @@ function plane_from_points(points::Vector{SVector{3, Float64}})
     return centroid, normalize(dir) 
 end
 
-"""
-    plane_from_points2(points::SMatrix{3, N, Float64}}) ->  centroid, normal
-
-Estimate the best fitting plane for a set of points in 3D.
-A more efficient version of plane_from_points.
-"""
-function plane_from_points2(points::SMatrix{3, N, Float64} where {N}) 
-    center = mean(points,dims=2)
-
-    u, _, _ = svd(points .- center)
-    normal = u[:,3]             # singular vectors in decending order
-
-    # make sure the normal is pointing consistently to positive Z direction 
-    if (dot(normal, unitZ3()) < 0.0)
-        normal = normal * -1.0
-    end
-
-    return SVector(center), SVector(normal)     # convert from SMatrix to SVector
-end
-
 
 function csg_sphere(;radius = 10.0)
     sph = Sphere(radius)
diff --git a/src/utilities.jl b/src/utilities.jl
index ede20cb81..95878069f 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -2,7 +2,28 @@
 # Copyright (c) Microsoft Corporation. All rights reserved.
 # See LICENSE in the project root for full license information.
 
-# only returns real roots
+
+"""
+    plane_from_points(points::SMatrix{3, N, Float64}}) ->  centroid, normal
+
+Estimate the best fitting plane for a set of points in 3D.
+A more efficient version of plane_from_points.
+"""
+function plane_from_points(points::SMatrix{3, N, Float64} where {N}) 
+    center = Statistics.mean(points,dims=2)
+
+    u, _, _ = svd(points .- center)
+    normal = u[:,3]             # singular vectors in decending order
+
+    # make sure the normal is pointing consistently to positive Z direction 
+    if (dot(normal, unitZ3()) < 0.0)
+        normal = normal * -1.0
+    end
+
+    return SVector(center), SVector(normal)     # convert from SMatrix to SVector
+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 9777abb12..b6f5adcd6 100644
--- a/test/testsets/Repeat.jl
+++ b/test/testsets/Repeat.jl
@@ -14,5 +14,19 @@
         return Repeat.ClusterWithProperties(lattice,properties)
     end
 
+    """ 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))
+        eltlattice = Repeat.HexBasis1()
+        clusterbasis = Repeat.LatticeBasis(( -1,2),(2,-1))
+        return Repeat.LatticeCluster(clusterbasis,eltlattice,clusterelts)
+    end
+
     @test [-1 2;2 -1] == Repeat.basismatrix(Repeat.clusterbasis(hex3RGB()))
+
+    function basistest(a::Repeat.AbstractLatticeCluster)
+        return Repeat.clusterbasis(a)
+    end
+
+    @test basistest(hex3cluster()) == basistest(hex3RGB()) 
 end
\ No newline at end of file