Brian gun/issue310 lenslet analysis

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"

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}

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}

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)
 

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

src/Optical/Eye.jl → 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

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")

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)))

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.