Merge branch 'main' into compathelper/new_version/2021-12-04-00-09-50…

…-225-01666007016

samples/notebooks/BasicCSG.jl

@@ -100,7 +100,7 @@ begin
 
 	# two transformed copies of the canonic bezier surface
 	csg1_surf2 = leaf(csg1_surf1, OpticSim.translation(-0.5, -0.5, 0.0))
-	csg1_surf3 = leaf(csg1_surf1, Transform{Float64}(0.0, Float64(π), 0.0, 0.5, -0.5, 0.0))
+	csg1_surf3 = leaf(csg1_surf1, Transform(0.0, Float64(π), 0.0, 0.5, -0.5, 0.0))
 
 	# transformed cilinder
 	csg1_surf4_1 = leaf(Cylinder(0.3, 1.0), OpticSim.translation(0.0, 0.0, 0.0))

samples/notebooks/Samples.jl

@@ -177,7 +177,7 @@ begin
 	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)), OpticSim.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 = (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0))
+	lens = (barrel ∩ topsurface ∩ botsurface)(Transform(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)	
 end

src/Examples/docs_examples.jl

@@ -130,7 +130,7 @@ function draw_lensconstruction(filename::Union{Nothing,AbstractString} = nothing
     barrel = Cylinder(
         9.0, 20.0, interface = FresnelInterface{Float64}(SCHOTT.N_BK7, Air, reflectance=0.0, transmission=0.0)
     )
-    lens = (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0))
+    lens = (barrel ∩ topsurface ∩ botsurface)(Transform(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0))
     detector = Rectangle(15.0, 15.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())
     sys = CSGOpticalSystem(LensAssembly(lens), detector)
 

src/Geometry/Transform.jl

@@ -133,15 +133,84 @@ Transform(rotation::AbstractArray{S,2}, translation::AbstractArray{S,1})
 ```
 `θ`, `ϕ` and `ψ` in first constructor are in **radians**.
 """
-Transform{T} = SMatrix{4,4,T,16}
+struct Transform{T} 
+    matrix::SMatrix{4,4,T,16}
+
+    """ This is a private internal function. In general don't want to allow users to populate Transform matrices with arbitrary elements. Not to be called by code outside of the Transform module. Don't use Transform{Float64}(...) for example. Instead use Transform(..)"""
+    function Transform{T}(a11::T,a21::T,a31::T,a41::T,a12::T,a22::T,a32::T,a42::T,a13::T,a23::T,a33::T,a43::T,a14::T,a24::T,a34::T,a44::T) where{T<:Real} 
+        return new{T}(SMatrix{4,4,T,16}(a11,a21,a31,a41,a12,a22,a32,a42,a13,a23,a33,a43,a14,a24,a34,a44))
+    end
+
+    Transform{T}(mat::SMatrix{4,4,T,16}) where{T<:Real} = new{T}(mat)
+end
 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}) ...
 
+#functions to make Transform compatible with base matrix API
+matrix(a::Transform) = a.matrix
+
+# Base.length(a::Transform) = length(matrix(a))
+Base.getindex(a::Transform, indices::Vararg{Int,N}) where{N} = getindex(matrix(a),indices...)
+# Base.iterate(a::Transform) = iterate(matrix(a))
+# Base.iterate(a::Transform{Float64}, b::Tuple{StaticArrays.SOneTo{16}, Int64}) = iterate(matrix(a),b)
+
+Base.collect(a::Transform) = collect(matrix(a))
+
+Base.:*(transa::Transform{T},transb::Transform{T}) where{T<:Real} = Transform{T}(matrix(transa)*matrix(transb))
+
+function Base.:*(transform::Transform{T}, v::SVector{3,S}) where {T<:Real,S<:Number}
+    t = matrix(transform)
+
+    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
+
+#Transform is not necessarily constrained to be a rigid body transformation so use general invers.
+Base.inv(a::Transform{T}) where{T<:Real}= Transform{T}(inv(matrix(a)))
+
+""" The t and m matrices are allowed to be of different element type. This allows transforming a Unitful matrix for example:
+```
+id = identitytransform()
+m = fill(1mm,3,4)
+id*m #returns a matrix filled with Unitful quantities. If both matrices had to be the same type this would not work
+```
+"""
+function Base.:*(transform::Transform{T}, m::SMatrix{3,N,S}) where{N,T<:Real,S<:Number}
+    res = MMatrix{3,N,T}(undef)
+    t = matrix(transform)
+
+    for outcol in 1:N
+        for row in 1:3
+            sum = T(0)
+            for incol in 1:3
+                sum += t[row,incol]*m[incol,outcol]
+            end
+            #implicit 1 w coordinate value
+            sum += t[row,4]
+            res[row,outcol] = sum
+        end
+        if t[4,4] != 1
+            res[:,outcol] /= t[4,4]
+        end
+    end
+    return SMatrix{3,N,S}(res)
+end
+
+""" The t and m matrices are allowed to be of different element type. This allows transforming a Unitful matrix for example:
+```
+id = identitytransform()
+m = fill(1mm,3,4)
+id*m #returns a matrix filled with Unitful quantities. If both matrices had to be the same type this would not work
+```
+"""
+Base.:*(transform::Transform{T},m::SMatrix{4,N,S}) where{N,T<:Real,S<:Number} = matrix(transform)*m
+
+Base.transpose(a::Transform{T}) where{T<:Real} = Transform{T}(a')
+
+# END of functions for compatibility with base matrix API
 
 
 # for compatability ith the "old" RigidBodyTransform
@@ -158,7 +227,6 @@ identitytransform(::Type{T} = Float64) where {T<:Real} = Transform{T}(
 )
 export identitytransform
 
-
 """
     Transform([S::Type]) -> Transform{S}
 
@@ -174,7 +242,7 @@ end
 Costruct a transform from the input columns.     
 """
 function Transform(colx::Vec3{T}, coly::Vec3{T}, colz::Vec3{T}, colw::Vec3{T} = zero(Vec3{T})) where {T<:Real}
-    return vcat(hcat(colx,coly,colz,colw),SMatrix{1,4,T}(zero(T),zero(T),zero(T),one(T)) )
+    return Transform{T}(vcat(hcat(colx,coly,colz,colw),SMatrix{1,4,T}(zero(T),zero(T),zero(T),one(T)) ))
 end
 
 
@@ -184,23 +252,22 @@ end
 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 hcat(colx,coly,colz,colw)
+    return Transform{T}(hcat(colx,coly,colz,colw))
 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}
+function Transform(origin::Vec3{T}, forward::Vec3{T} = unitZ3()) 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)
+function Transform(θ::T, ϕ::T, ψ::T, x::T, y::T, z::T) where {T<:Number} 
+    return Transform(rotmat(T, θ, ϕ, ψ), Vec3{T}(x, y, z))
 end
 
 """
@@ -209,7 +276,7 @@ end
 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(
+    return Transform{T}(
         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), 
@@ -223,7 +290,7 @@ Returns the [`Transform`](@ref) of type `S` (default `Float64`) created by a rot
 """
 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(
+    return Transform{T}(
         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), 
@@ -369,6 +436,7 @@ export translation, rotation, rotationd
 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)
@@ -441,36 +509,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))
-        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
 
-#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)
-
-    for outcol in 1:N
-        for row in 1:3
-            sum = T(0)
-            for incol in 1:3
-                sum += t[row,incol]*m[incol,outcol]
-            end
-            #implicit 1 w coordinate value
-            sum += t[row,4]
-            res[row,outcol] = sum
-        end
-        if t[4,4] != 1
-            res[:,outcol] /= t[4,4]
-        end
-    end
-    return SMatrix{3,N,T}(res)
-end
 
 """
     decomposeRTS(tr::Transform{T}) where {T<:Real}

src/Optical/Lenses.jl

@@ -104,7 +104,7 @@ function AsphericLens(insidematerial::T, frontvertex::S, frontradius::S, frontco
             backaspherics = Tuple{Int,S}.(backaspherics)
         end
         surf = AcceleratedParametricSurface(AsphericSurface(semidiameter + backdecenter_l + S(0.01), radius = backradius, conic = backconic, aspherics = backaspherics), nsamples, interface = opticinterface(S, insidematerial, nextmaterial, backsurfacereflectance, interfacemode))
-        lens_rear = leaf(surf, Transform{S}(zero(S), S(π), zero(S), backdecenter[1], backdecenter[2], frontvertex - thickness))
+        lens_rear = leaf(surf, Transform(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)

src/RepeatingStructures/Multilens/LensletAssignment.jl

@@ -0,0 +1,55 @@
+# MIT license
+# Copyright (c) Microsoft Corporation. All rights reserved.
+# See LICENSE in the project root for full license information.
+
+
+""" Computes the location of the optical center of the lens that will project the centroid of the display to the centroid of the eyebox. Normally the display centroid will be aligned with the geometric centroid of the lens, rather than the optical center of the lens."""
+function compute_optical_center(eyeboxcentroid,displaycenter,lensplane)
+    v = displaycenter - eyeboxcentroid
+    r = Ray(eyeboxcentroid,v)
+    return closestintersection(surfaceintersection(lensplane))
+end
+
+function subdivide_fov(eyeboxpolygon)
+end
+
+""" eyebox rectangle represented as a 3x4 SMatrix with x,y,z coordinates in rows 1,2,3 respectively."""
+function eyeboxpolygon(xsize,ysize)
+    xext,yext = (xsize,ysize) ./ 2.0 #divide by 2 to get min and max extensions
+    return SMatrix{3,4}([xext;yext;0.0 ;; xext;-yext;0.0 ;; -xext;-yext;0.0 ;; -xext;yext;0.0])
+end
+
+subpoints(pts::AbstractVector,subdivisions) = reshape(collect(range(extrema(pts)...,subdivisions)),1,subdivisions)
+export subpoints
+
+subdivide(poly::AbstractMatrix,xsubdivisions,ysubdivisions) = subdivide(SMatrix{3,4}(poly),xsubdivisions,ysubdivisions)
+
+"""poly is a two dimensional rectangle represented as a 3x4 SMatrix with x values in row 1 and y values in row 2. Returns a vector of length xsubdivisions*ysubdivisions containing all the subdivided rectangles."""
+function subdivide(poly::T, xsubdivisions, ysubdivisions) where{T<:SMatrix{3,4}}
+    result = Vector{T}(undef,0) #efficiency not an issue here so don't need to preallocate vector
+
+    x = subpoints(poly[1,:],xsubdivisions+1) #range function makes one less subdivision that people will be expecting, e.g., collect(range(1,3,2)) returns [1,3] 
+    y = subpoints(poly[2,:],ysubdivisions+1)
+
+    for i in 1:length(x) -1
+        for j in 1:length(y) -1
+            push!(result,T([x[i];y[j];0.0 ;; x[i+1];y[j];0.0 ;; x[i+1];y[j+1];0.0 ;; x[i];y[j+1];0.0]))
+        end
+    end
+    return result
+end
+export subdivide
+
+function setup_system(subdivisions::NTuple{2,Int})
+    eyeballframe = Transform(I(4)) #This establishes the global coordinate frame for the systems. Positive Z axis is assumed to be the forward direction, the default direction of the eye's optical axis when looking directly ahead.
+    corneavertex = OpticSim.HumanEye.cornea_to_eyecenter()
+
+
+    eyeboxframe = eyeballframe*OpticSim.translation(0.0,0.0,ustrip(mm,corneavertex))  #unfortunately can't use unitful values in transforms because the rotation and translation components would have different types which is not allowed in a Matrix.
+
+    println(typeof(eyeboxframe))
+    eyeboxpoly =  eyeboxpolygon(10mm,9mm) #four corners of the eyebox frame which is assumed centered around the positive Z axis. Transformed to the eyeballframe.
+    println(typeof(subdivide(eyeboxpoly,subdivisions...)[1]))
+    subdivided_eyeboxpolys = map(x-> eyeboxframe .* x, subdivide(eyeboxpoly,subdivisions...))  
+end
+export setup_system

src/RepeatingStructures/Multilens/Lenslets.jl

@@ -13,11 +13,14 @@ import ..Repeat
 import SpecialFunctions
 import Plots
 import ...OpticSim
+import ...OpticSim.Geometry
 
 include("HexClusters.jl")
 include("HexTilings.jl")
 include("Analysis.jl")
 include("DisplayGeneration.jl")
+include("LensletAssignment.jl")
+
 
 end # module
 export Lenslets

src/Vis/VisRepeatingStructures.jl

@@ -63,16 +63,16 @@ function drawcells(tilebasis::AbstractBasis, tilesize, cells::AbstractMatrix; co
 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, cluster_coordinate_offset=[0;0;;], scale=50.0) = drawcells(clstr.elementbasis, scale, clustercoordinates(clstr, cluster_coordinate_offset[1,1],cluster_coordinate_offset[2,1]))
 
-""" 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.
-    dims = size(cells)
+""" draw the ClusterWithProperties at coordinates specified by lattice_coordinate_offset """
+function draw(clstr::Repeat.ClusterWithProperties, cluster_coordinate_offset=[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(cluster_coordinate_offset)
     @assert dims[1] == 2 "dims[1] must be 2 instead was $(dims[1])"
     clstrsize = clustersize(clstr)
     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, cluster_coordinate_offset[1,i], cluster_coordinate_offset[2,i])
     end
     props = repeat(Repeat.properties(clstr), dims[2])
     drawcells(elementbasis(cluster(clstr)), scale, points, color=props[:,:Color], name=props[:,:Name])