From d4e4c12de6855b3c6582ba31e6af065c65b1adf9 Mon Sep 17 00:00:00 2001 From: mattderstine <59268355+mattderstine@users.noreply.github.com> Date: Thu, 11 Nov 2021 15:39:42 -0800 Subject: [PATCH 1/2] Add aspheric surface with odd and even terms (#291) (#300) * first step of changes to move Aspheric function to it's own Surface type * Fix Geometry.jl to include AsphericSurface.jl * change to include A1 and up * Working AsphericSurface type, tests etc * fix documentation for AsphericSurface * all working, but allocation problems * Simple change to test * Fix variable names of ZernikeSurface * evaluate based on on aspheric pseudo-types * Update src/Geometry/Primitives/AsphericSurface.jl Co-authored-by: Charlie Hewitt * Update src/Geometry/Primitives/Zernike.jl Co-authored-by: Charlie Hewitt * Update src/Geometry/Primitives/Zernike.jl Co-authored-by: Charlie Hewitt * Update src/Geometry/Primitives/AsphericSurface.jl Co-authored-by: Charlie Hewitt * Update src/Geometry/Primitives/AsphericSurface.jl Co-authored-by: Charlie Hewitt * fix argument lists, other minor mistakes * 1st version to pass all tests. Still need benchmark tests and documentation * Ready for final check * fixed: documentation, last tests, axisym optisur * Update ref.md Co-authored-by: mattderstine Co-authored-by: BrianGun <36491850+BrianGun@users.noreply.github.com> Co-authored-by: Charlie Hewitt --- docs/src/primitives.md | 1 + src/Geometry/Geometry.jl | 1 + src/Geometry/Primitives/AsphericSurface.jl | 341 +++++++++++++++++++++ src/Geometry/Primitives/Zernike.jl | 132 +++----- src/Optical/Lenses.jl | 4 +- test/Benchmarks/Benchmarks.jl | 5 +- test/TestData/lenses.jl | 20 +- test/TestData/surfaces.jl | 6 +- test/TestData/systems.jl | 4 +- test/testsets/Comparison.jl | 46 ++- test/testsets/Intersection.jl | 1 - 11 files changed, 460 insertions(+), 101 deletions(-) create mode 100644 src/Geometry/Primitives/AsphericSurface.jl diff --git a/docs/src/primitives.md b/docs/src/primitives.md index 6b3058608..12546a507 100644 --- a/docs/src/primitives.md +++ b/docs/src/primitives.md @@ -60,6 +60,7 @@ Cylinder Plane Sphere SphericalCap +AsphericSurface ZernikeSurface BezierSurface BSplineSurface diff --git a/src/Geometry/Geometry.jl b/src/Geometry/Geometry.jl index 011158618..9f7dfab68 100644 --- a/src/Geometry/Geometry.jl +++ b/src/Geometry/Geometry.jl @@ -51,6 +51,7 @@ include("Primitives/Curves/PowerBasis.jl") include("CSG/CSG.jl") +include("Primitives/AsphericSurface.jl") include("Primitives/Zernike.jl") include("Primitives/Qtype.jl") include("Primitives/Chebyshev.jl") diff --git a/src/Geometry/Primitives/AsphericSurface.jl b/src/Geometry/Primitives/AsphericSurface.jl new file mode 100644 index 000000000..ec19bd92d --- /dev/null +++ b/src/Geometry/Primitives/AsphericSurface.jl @@ -0,0 +1,341 @@ +# MIT license +# Copyright (c) Microsoft Corporation. All rights reserved. +# See LICENSE in the project root for full license information. + + + +######################################################################################################### + +#"pseudo-types" of aspheres +""" +AphericSurfaces polynomial evaluation is optimized for no terms `CONIC`, odd terms `ODD`, even terms `EVEN` or any `ODDEVEN` +""" +@enum AsphericSurfaceType CONIC ODD EVEN ODDEVEN + + +""" + AsphericSurface{T,N,Q,M} <: ParametricSurface{T,N} + +Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter,. +`T` is the datatype, `N` is the dimensionality, `Q` is the number of aspheric terms, and `M` is the type of aspheric polynomial. + + +```julia +AsphericSurface(semidiameter; radius, conic, aspherics=nothing, normradius = semidiameter) +``` + +The surface is centered at the origin and treated as being the cap of an infinite cylinder, thus creating a true half-space. +Outside of `0 <= ρ <= 1` the height of the surface is not necessarily well defined, so NaN may be returned. + +`aspherics` aspherics should be vectors containing tuples of the form (i, v) where i is the polynomial power of the aspheric term + +The sag is defined by the equation + +```math +z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir^{i} +``` + +where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, and ``k = \\texttt{conic}`` . + +The function checks if the aspheric terms are missing, even, odd or both and uses an efficient polynomial evaluation strategy. + +""" +struct AsphericSurface{T,N,Q,M} <: ParametricSurface{T,N} + semidiameter::T + curvature::T + conic::T + aspherics::SVector{Q,T} + normradius::T + boundingcylinder::Cylinder{T,N} + + function AsphericSurface(M::AsphericSurfaceType, semidiameter::T, curvature::T, conic::T, aspherics::SVector{Q,T}, normradius::T, boundingcylinder) where {T<:Real,Q} + new{T,3,Q,M}(semidiameter, curvature, conic, aspherics, normradius, boundingcylinder) + end + + function AsphericSurface(semidiameter::T; radius::T = typemax(T), conic::T = zero(T), aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter) where {T<:Real} + @assert semidiameter > 0 + @assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic) + @assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)" + + acs = [] + if aspherics===nothing + M = CONIC + else + asphericTerms = [i for (i, ) in aspherics] + minAsphericTerm = minimum(asphericTerms) + maxAsphericTerm = maximum(asphericTerms) + @assert minAsphericTerm > 0 "Aspheric Terms must be Order 1 or higher ($minAsphericTerm)" + acs = zeros(T, maxAsphericTerm ) + for (i, k) in aspherics + acs[i] = k + end + odd = any([isodd(i) && a != zero(T) for (i, a) in enumerate(acs)]) + even = any([iseven(i) && a != zero(T) for (i, a) in enumerate(acs)]) + if odd && even + M = ODDEVEN + elseif even + M = EVEN + acs = [acs[2i] for i in 1:(maxAsphericTerm ÷ 2)] + elseif odd + M = ODD + acs = [acs[2i-1] for i in 1:((maxAsphericTerm+1) ÷ 2)] + else #there are no nonzero aspherics terms in the list + M = CONIC + acs = [] + end + end + Q = length(acs) + surf = new{T,3, Q, M}(semidiameter, 1/radius, conic, acs, normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) + return surf + end + +end + +""" + EvenAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) + +Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter. + +`aspherics` should be an array of the even coefficients of the aspheric polynomial starting with A2 + + +""" +function EvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real + Q=length(aspherics) + AsphericSurface(EVEN, semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) +end + +""" + OddAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) + +Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter. + +`aspherics` should be an array of the odd coefficients of the aspheric polynomial starting with A1 + +""" +function OddAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real + Q=length(aspherics) + AsphericSurface(ODD, semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) +end + +""" + OddEvenAsphericSurface(semidiameter, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) + +Surface incorporating an aspheric polynomial - radius, conic and aspherics are defined relative to absolute semi-diameter. + +`aspherics` should be an array of the both odd and even coefficients of the aspheric polynomial starting with A1 + +""" +function OddEvenAsphericSurface(semidiameter::T, curvature::T, conic::T, aspherics::Vector{T}; normradius::T=semidiameter) where T<:Real + Q=length(aspherics) + AsphericSurface(ODDEVEN, semidiameter, curvature, conic, SVector{Q,T}(aspherics), normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) +end + +#don't have a function for CONIC as it would better to directly solve for the surface intersection instead of the rootfinding of a ParametricSurface + +""" + asphericType(surf::AsphericSurface) + +Query the polynomial type of `asp. Returns CONIC, ODD, EVEN, or ODDEVEN. CONIC corresponds to no aspheric terms, ODD +means it only has odd aspheric terms, EVEN means only even aspheric terms and ODDEVEN means both even and odd terms. + +This function is to enable proper interpretation of `surf.aspherics` by any optimization routines that directly query the aspheric coefficients. +""" +asphericType(z::AsphericSurface{T,3,Q,M}) where {T<:Real,Q,M} = M + +export AsphericSurface, EvenAsphericSurface, OddAsphericSurface, OddEvenAsphericSurface, asphericType, EVEN, CONIC, ODD, ODDEVEN + + +uvrange(::AsphericSurface{T,N}) where {T<:Real,N} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ + +semidiameter(z::AsphericSurface{T}) where {T<:Real} = z.semidiameter +halfsizeu(z::AsphericSurface{T}) where {T<:Real} = semidiameter(z) +halfsizev(z::AsphericSurface{T}) where {T<:Real} = semidiameter(z) + +boundingobj(z::AsphericSurface{T}) where {T<:Real} = z.boundingcylinder + +prod_step(z::AsphericSurface{T,N,Q,ODDEVEN}, r, r2) where {T<:Real,N,Q} = r, r +prod_step(z::AsphericSurface{T,N,Q,ODD}, r, r2) where {T<:Real,N,Q} = r, r2 +prod_step(z::AsphericSurface{T,N,Q,EVEN}, r, r2) where {T<:Real,N,Q} = r2, r2 + +function point(z::AsphericSurface{T,3,Q,M}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,Q,M} + rad = z.semidiameter + r = ρ * rad + r2 = r^2 + t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 + if t < zero(T) + return SVector{3,T}(NaN, NaN, NaN) + end + h = z.curvature * r2 / (one(T) + sqrt(t)) + # sum aspheric + if M != CONIC + prod, step = prod_step(z, r, r2) #multiple dispatch on M + asp,rest = Iterators.peel(z.aspherics) + h += asp * prod + for asp in rest + prod *= step + h += asp * prod + end + end + return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) +end + +partial_prod_step(z::AsphericSurface{T,3,Q,EVEN}, r::T, r2::T) where {T<:Real,Q} = r, r2, 2:2:2Q +partial_prod_step(z::AsphericSurface{T,3,Q,ODD}, r::T, r2::T) where {T<:Real,Q} = one(T), r2, 1:2:(2Q-1) +partial_prod_step(z::AsphericSurface{T,3,Q,ODDEVEN}, r::T, r2::T) where {T<:Real,Q} = one(T), r, 1:1:Q + +function partials(z::AsphericSurface{T,3,Q,M}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,Q,M} + rad = z.semidiameter + r = ρ * rad + r2 = r*r + t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 + if t < zero(T) + return SVector{3,T}(NaN, NaN, NaN), SVector{3,T}(NaN, NaN, NaN) + end + dhdρ = rad * z.curvature * r * sqrt(t) / t + # sum aspherics partial + if M != CONIC + prod, step, mIter = partial_prod_step(z, r, r2) + + ((m, asp), rest) = Iterators.peel(zip(mIter,z.aspherics)) + dhdρ += rad * asp * 2 * prod #first term m=1*2 and prod = r + for (m,asp) in rest + prod *= step + dhdρ += rad * m * asp * prod + end + end + dhdϕ = zero(T) + cosϕ = cos(ϕ) + sinϕ = sin(ϕ) + return SVector{3,T}(rad * cosϕ, rad * sinϕ, dhdρ), SVector{3,T}(r * -sinϕ, r * cosϕ, dhdϕ) +end + +function normal(z::AsphericSurface{T,3}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real} + du, dv = partials(z, ρ, ϕ) + if ρ == zero(T) && norm(dv) == zero(T) + # in cases where there is no δϕ at ρ = 0 (i.e. anything which is rotationally symetric) + # then we get some big problems, hardcoding this case solves the problems + return SVector{3,T}(0, 0, 1) + end + return normalize(cross(du, dv)) +end + +function uv(z::AsphericSurface{T,3}, p::SVector{3,T}) where {T<:Real} + # avoid divide by zero for ForwardDiff + ϕ = NaNsafeatan(p[2], p[1]) + if p[1] == zero(T) && p[2] == zero(T) + ρ = zero(T) + else + ρ = sqrt(p[1]^2 + p[2]^2) / semidiameter(z) + end + return SVector{2,T}(ρ, ϕ) +end + +function onsurface(surf::AsphericSurface{T,3}, p::SVector{3,T}) where {T<:Real} + ρ, ϕ = uv(surf, p) + if ρ > one(T) + return false + else + surfpoint = point(surf, ρ, ϕ) + return samepoint(p[3], surfpoint[3]) + end +end + +function inside(surf::AsphericSurface{T,3}, p::SVector{3,T}) where {T<:Real} + ρ, ϕ = uv(surf, p) + if ρ > one(T) + return false + else + surfpoint = point(surf, ρ, ϕ) + return p[3] < surfpoint[3] + end +end + +######################################################################################################### + +# Assumes the ray has been transformed into the canonical coordinate frame which has the vertical axis passing through (0,0,0) and aligned with the z axis. +function surfaceintersection(surf::AcceleratedParametricSurface{T,3,AsphericSurface{T,3}}, r::AbstractRay{T,3}) where {T<:Real} + cylint = surfaceintersection(surf.surface.boundingcylinder, r) + if cylint isa EmptyInterval{T} + return EmptyInterval(T) + else + if doesintersect(surf.triangles_bbox, r) || inside(surf.triangles_bbox, origin(r)) + surfint = triangulatedintersection(surf, r) + if !(surfint isa EmptyInterval{T}) + return intervalintersection(cylint, surfint) + end + end + # hasn't hit the surface + if lower(cylint) isa RayOrigin{T} && upper(cylint) isa Infinity{T} + if inside(surf.surface, origin(r)) + return Interval(RayOrigin(T), Infinity(T)) + else + return EmptyInterval(T) + end + # otherwise check that the intersection is underneath the surface + else + p = point(closestintersection(cylint, false)) + ρ, ϕ = uv(surf, p) + surfpoint = point(surf.surface, ρ, ϕ) + if p[3] < surfpoint[3] + return cylint # TODO!! UV (and interface) issues? + else + return EmptyInterval(T) + end + end + end +end + +function AcceleratedParametricSurface(surf::T, numsamples::Int = 17; interface::NullOrFresnel{S} = NullInterface(S)) where {S<:Real,N,T<:AsphericSurface{S,N}} + # Zernike uses ρ, ϕ uv space so need to modify extension of triangulation + a = AcceleratedParametricSurface(surf, triangulate(surf, numsamples, true, false, true, false), interface = interface) + emptytrianglepool!(S) + return a +end + +function asphericSag(surf::AsphericSurface{T,3,Q, EVEN}) where {T<:Real,Q} + amin = sum(k < zero(T) ? k * surf.semidiameter^(2m) : zero(T) for (m, k) in enumerate(surf.aspherics)) + amax = sum(k > zero(T) ? k * surf.semidiameter^(2m) : zero(T) for (m, k) in enumerate(surf.aspherics)) + return amin, amax +end + +function asphericSag(surf::AsphericSurface{T,3,Q, ODD}) where {T<:Real,Q} + amin = sum(k < zero(T) ? k * surf.semidiameter^(2m-1) : zero(T) for (m, k) in enumerate(surf.aspherics)) + amax = sum(k > zero(T) ? k * surf.semidiameter^(2m-1) : zero(T) for (m, k) in enumerate(surf.aspherics)) + return amin, amax +end + +function asphericSag(surf::AsphericSurface{T,3,Q, ODDEVEN}) where {T<:Real,Q} + amin = sum(k < zero(T) ? k * surf.semidiameter^(m) : zero(T) for (m, k) in enumerate(surf.aspherics)) + amax = sum(k > zero(T) ? k * surf.semidiameter^(m) : zero(T) for (m, k) in enumerate(surf.aspherics)) + return amin, amax +end + +function asphericSag(surf::AsphericSurface{T,3,Q, CONIC}) where {T<:Real,Q} + amin = zero(T) + amax = zero(T) + return amin, amax +end + +function BoundingBox(surf::AsphericSurface{T,3}) where {T<:Real} + xmin = -semidiameter(surf) + xmax = semidiameter(surf) + ymin = -semidiameter(surf) + ymax = semidiameter(surf) + # aspherics can be more comlpicated, so just take sum of all negative and all positives + amin, amax = asphericSag(surf) + zmin = amin + zmax = amax + # curvature only goes one way + q = one(T) - (one(T) + surf.conic) * surf.curvature^2 * surf.semidiameter^2 + if q < zero(T) + throw(ErrorException("The surface is invalid, no bounding box can be constructed")) + end + hmax = surf.curvature * surf.semidiameter^2 / (one(T) + sqrt(q)) + if hmax > zero(T) + zmax += hmax + else + zmin += hmax + end + return BoundingBox(xmin, xmax, ymin, ymax, zmin, zmax) +end diff --git a/src/Geometry/Primitives/Zernike.jl b/src/Geometry/Primitives/Zernike.jl index 3bb0180ff..b3df4a5d5 100644 --- a/src/Geometry/Primitives/Zernike.jl +++ b/src/Geometry/Primitives/Zernike.jl @@ -116,10 +116,10 @@ Either `ZernikeIndexingOSA` or `ZernikeIndexingNoll`, see [Zernike polynomials w export ZernikeIndexType, ZernikeIndexingOSA, ZernikeIndexingNoll """ - ZernikeSurface{T,N,P,Q} <: ParametricSurface{T,N} + ZernikeSurface{T,N,P,Q,M} <: ParametricSurface{T,N} Surface incorporating the Zernike polynomials - radius, conic and aspherics are defined relative to absolute semi-diameter, Zernike terms are normalized according to the `normradius` parameter. -`T` is the datatype, `N` is the dimensionality, `P` is the number of Zernike terms and `Q` is the number of aspheric terms. Only even aspheric terms are supported. +`T` is the datatype, `N` is the dimensionality, `P` is the number of Zernike terms, `Q` is the number of aspheric terms and `M` is the Aspheric Type. The surface is centered at the origin and treated as being the cap of an infinite cylinder, thus creating a true half-space. Outside of `0 <= ρ <= 1` the height of the surface is not necessarily well defined, so NaN may be returned. @@ -130,7 +130,10 @@ For convenience the input `zcoeff` can be indexed using either OSA or Noll conve ZernikeSurface(semidiameter, radius = Inf, conic = 0, zcoeff = nothing, aspherics = nothing, normradius = semidiameter, indexing = ZernikeIndexingOSA) ``` -`zcoeff` and `aspherics` should be vectors containing tuples of the form `(i, v)` where `i` is either the index of the Zernike term for the corresponding `indexing`, or the polynomial power of the aspheric term (must be even) and `v` is the corresponding coefficient ``A_i`` or ``\\alpha_i`` respectively.. +`zcoeff` and `aspherics` should be vectors containing tuples of the form `(i, v)` where `i` is either the index of the Zernike term +for the corresponding `indexing`, or the polynomial power of the aspheric term (may be even or odd) and +`v` is the corresponding coefficient ``A_i`` or ``\\alpha_i`` respectively.. `M` will be determined from the terms entered to optimize +the evaluation of the aspheric polynomial. The sag is defined by the equation @@ -140,19 +143,15 @@ z(r,\\phi) = \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\sum_{i}^{Q}\\alpha_ir where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, ``k = \\texttt{conic}`` and ``Z_n`` is the nᵗʰ Zernike polynomial. """ -struct ZernikeSurface{T,N,P,Q} <: ParametricSurface{T,N} - semidiameter::T +struct ZernikeSurface{T,N,P,Q,M} <: ParametricSurface{T,N} + asp::AsphericSurface{T,N,Q,M} coeffs::SVector{P,Tuple{Int,Int,T}} - aspherics::SVector{Q,Tuple{Int,T}} - curvature::T - conic::T - normradius::T boundingcylinder::Cylinder{T,N} function ZernikeSurface(semidiameter::T; radius::T = typemax(T), conic::T = zero(T), zcoeff::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, aspherics::Union{Nothing,Vector{Tuple{Int,T}}} = nothing, normradius::T = semidiameter, indexing::ZernikeIndexType = ZernikeIndexingOSA) where {T<:Real} - @assert semidiameter > 0 - @assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic) - @assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)" + asp = AsphericSurface(semidiameter; radius, conic, aspherics, normradius) + Q = length(asp.aspherics) #this is not the same as the aspherics variable passed to the function! + zcs = [] if zcoeff !== nothing for (i, k) in zcoeff @@ -167,84 +166,54 @@ struct ZernikeSurface{T,N,P,Q} <: ParametricSurface{T,N} end end P = length(zcs) - acs = [] - if aspherics !== nothing - for (i, k) in aspherics - if abs(k) > zero(T) - @assert i % 2 == 0 "Only even aspheric terms are supported" - push!(acs, (i ÷ 2, k)) - end - end - end - Q = length(acs) - new{T,3,P,Q}(semidiameter, SVector{P,Tuple{Int,Int,T}}(zcs), SVector{Q,Tuple{Int,T}}(acs), 1 / radius, conic, normradius, Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder + M = asphericType(asp) + new{T,3,P,Q,M}(asp::AsphericSurface{T,3,Q,M}, SVector{P,Tuple{Int,Int,T}}(zcs), Cylinder(semidiameter, interface = opaqueinterface(T))) # TODO!! incorrect interface on cylinder end end export ZernikeSurface -uvrange(::Type{ZernikeSurface{T,N,P,Q}}) where {T<:Real,N,P,Q} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ +uvrange(::Type{ZernikeSurface{T,N,P,Q,M}}) where {T<:Real,N,P,Q,M} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and ϕ -semidiameter(z::ZernikeSurface{T}) where {T<:Real} = z.semidiameter +semidiameter(z::ZernikeSurface{T}) where {T<:Real} = z.asp.semidiameter halfsizeu(z::ZernikeSurface{T}) where {T<:Real} = semidiameter(z) halfsizev(z::ZernikeSurface{T}) where {T<:Real} = semidiameter(z) boundingobj(z::ZernikeSurface{T}) where {T<:Real} = z.boundingcylinder -function point(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,P,Q} - rad = z.semidiameter - r = ρ * rad - r2 = r^2 - t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 - if t < zero(T) - return SVector{3,T}(NaN, NaN, NaN) - end - h = z.curvature * r2 / (one(T) + sqrt(t)) - # sum aspheric - @inbounds @simd for i in 1:Q - (m, k) = z.aspherics[i] - h += k * r2^m - end +function point(z::ZernikeSurface{T,3,P,Q,M}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,P,Q,M} + pnt = point(z.asp, ρ, ϕ) + # sum zernike - u = r / z.normradius + rad = semidiameter(z.asp) + r = ρ * rad + u = r / z.asp.normradius + h = zero(T) @inbounds @simd for m in 1:P (R, S, k) = z.coeffs[m] h += k * Zernike.ζ(R, S, u, ϕ) end - return SVector{3,T}(r * cos(ϕ), r * sin(ϕ), h) + return SVector{3,T}(pnt[1], pnt[2], pnt[3] + h) end -function partials(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,P,Q} - rad = z.semidiameter - r = ρ * rad - t = one(T) - z.curvature^2 * (z.conic + one(T)) * r^2 - if t < zero(T) - return SVector{3,T}(NaN, NaN, NaN), SVector{3,T}(NaN, NaN, NaN) - end - dhdρ = rad * z.curvature * r * sqrt(t) / t - # sum aspherics partial - @inbounds @simd for i in 1:Q - (m, k) = z.aspherics[i] - dhdρ += rad * (2 * m) * k * r^(2 * m - 1) - end +function partials(z::ZernikeSurface{T,3,P,Q,M}, ρ::T, ϕ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,P,Q,M} + pρ,pϕ = partials(z.asp, ρ, ϕ) # sum zernike partials - dhdϕ = zero(T) - n = rad / z.normradius + rad=z.asp.semidiameter + n = rad / z.asp.normradius u = ρ * n + dhdρ = zero(T) + dhdϕ = zero(T) @inbounds @simd for m in 1:P (R, S, k) = z.coeffs[m] du, dϕ = Zernike.δζ(R, S, u, ϕ) dhdρ += k * du * n # want the derivative wrt ρ, not u dhdϕ += k * dϕ end - cosϕ = cos(ϕ) - sinϕ = sin(ϕ) - pu = SVector{3,T}(rad * cosϕ, rad * sinϕ, dhdρ) - pv = SVector{3,T}(r * -sinϕ, r * cosϕ, dhdϕ) - return pu, pv + return SVector{3,T}(pρ[1], pρ[2], pρ[3] + dhdρ), SVector{3,T}(pϕ[1], pϕ[2], pϕ[3] + dhdϕ) end -function normal(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,P,Q} +function normal(z::ZernikeSurface{T,3,P,Q,M}, ρ::T, ϕ::T)::SVector{3,T} where {T<:Real,P,Q,M} du, dv = partials(z, ρ, ϕ) if ρ == zero(T) && norm(dv) == zero(T) # in cases where there is no δϕ at ρ = 0 (i.e. anything which is rotationally symetric) @@ -254,7 +223,7 @@ function normal(z::ZernikeSurface{T,3,P,Q}, ρ::T, ϕ::T)::SVector{3,T} where {T return normalize(cross(du, dv)) end -function uv(z::ZernikeSurface{T,3,P,Q}, p::SVector{3,T}) where {T<:Real,P,Q} +function uv(z::ZernikeSurface{T,3,P,Q,M}, p::SVector{3,T}) where {T<:Real,P,Q,M} # avoid divide by zero for ForwardDiff ϕ = NaNsafeatan(p[2], p[1]) if p[1] == zero(T) && p[2] == zero(T) @@ -265,7 +234,7 @@ function uv(z::ZernikeSurface{T,3,P,Q}, p::SVector{3,T}) where {T<:Real,P,Q} return SVector{2,T}(ρ, ϕ) end -function onsurface(surf::ZernikeSurface{T,3,P,Q}, p::SVector{3,T}) where {T<:Real,P,Q} +function onsurface(surf::ZernikeSurface{T,3,P,Q,M}, p::SVector{3,T}) where {T<:Real,P,Q,M} ρ, ϕ = uv(surf, p) if ρ > one(T) return false @@ -275,7 +244,7 @@ function onsurface(surf::ZernikeSurface{T,3,P,Q}, p::SVector{3,T}) where {T<:Rea end end -function inside(surf::ZernikeSurface{T,3,P,Q}, p::SVector{3,T}) where {T<:Real,P,Q} +function inside(surf::ZernikeSurface{T,3,P,Q,M}, p::SVector{3,T}) where {T<:Real,P,Q,M} ρ, ϕ = uv(surf, p) if ρ > one(T) return false @@ -288,7 +257,7 @@ end ######################################################################################################### # Assumes the ray has been transformed into the canonical coordinate frame which has the vertical axis passing through (0,0,0) and aligned with the z axis. -function surfaceintersection(surf::AcceleratedParametricSurface{T,3,ZernikeSurface{T,3,P,Q}}, r::AbstractRay{T,3}) where {T<:Real,P,Q} +function surfaceintersection(surf::AcceleratedParametricSurface{T,3,ZernikeSurface{T,3,P,Q,M}}, r::AbstractRay{T,3}) where {T<:Real,P,Q,M} cylint = surfaceintersection(surf.surface.boundingcylinder, r) if cylint isa EmptyInterval{T} return EmptyInterval(T) @@ -327,31 +296,12 @@ function AcceleratedParametricSurface(surf::T, numsamples::Int = 17; interface:: return a end -function BoundingBox(surf::ZernikeSurface{T,3,P,Q}) where {T<:Real,P,Q} - xmin = -semidiameter(surf) - xmax = semidiameter(surf) - ymin = -semidiameter(surf) - ymax = semidiameter(surf) +function BoundingBox(surf::ZernikeSurface{T,3,P,Q,M}) where {T<:Real,P,Q,M} + bb = BoundingBox(z.asp) # zernike terms have condition than |Zᵢ| <= 1 # so this gives us a (loose) bounding box - ak = P > 0 ? sum(abs.(Zernike.normalisation(T, N, M) * k for (N, M, k) in surf.coeffs)) : zero(T) - zmin = -ak - zmax = ak - # aspherics can be more comlpicated, so just take sum of all negative and all positives - amin = Q > 0 ? sum(k < zero(T) ? k * surf.semidiameter^(2m) : zero(T) for (m, k) in surf.aspherics) : zero(T) - amax = Q > 0 ? sum(k > zero(T) ? k * surf.semidiameter^(2m) : zero(T) for (m, k) in surf.aspherics) : zero(T) - zmin += amin - zmax += amax - # curvature only goes one way - q = one(T) - (one(T) + surf.conic) * surf.curvature^2 * surf.semidiameter^2 - if q < zero(T) - throw(ErrorException("The surface is invalid, no bounding box can be constructed")) - end - hmax = surf.curvature * surf.semidiameter^2 / (one(T) + sqrt(q)) - if hmax > zero(T) - zmax += hmax - else - zmin += hmax - end - return BoundingBox(xmin, xmax, ymin, ymax, zmin, zmax) + ak = P > 0 ? sum(abs.(Zernike.normalisation(T, n, m) * k for (n, m, k) in surf.coeffs)) : zero(T) + bb.zmin -= ak + bb.zmax += ak + return BoundingBox(bb.xmin, bb.xmax, bb.ymin, bb.ymax, bb.zmin, bb.zmax) #could just return bb, but this way is safer end diff --git a/src/Optical/Lenses.jl b/src/Optical/Lenses.jl index 85396cac4..eb81ae196 100644 --- a/src/Optical/Lenses.jl +++ b/src/Optical/Lenses.jl @@ -69,7 +69,7 @@ function AsphericLens(insidematerial::T, frontvertex::S, frontradius::S, frontco if frontaspherics !== nothing frontaspherics = [(i, -k) for (i, k) in frontaspherics] end - surf = AcceleratedParametricSurface(ZernikeSurface(semidiameter + frontdecenter_l + S(0.01), radius = -frontradius, conic = frontconic, aspherics = frontaspherics), nsamples, interface = opticinterface(S, insidematerial, lastmaterial, frontsurfacereflectance, interfacemode)) + surf = AcceleratedParametricSurface(AsphericSurface(semidiameter + frontdecenter_l + S(0.01), radius = -frontradius, conic = frontconic, aspherics = frontaspherics), nsamples, interface = opticinterface(S, insidematerial, lastmaterial, frontsurfacereflectance, interfacemode)) lens_front = leaf(surf, translation(S, frontdecenter[1], frontdecenter[2], frontvertex)) end if isinf(backradius) && (backaspherics === nothing) @@ -103,7 +103,7 @@ function AsphericLens(insidematerial::T, frontvertex::S, frontradius::S, frontco if backaspherics !== nothing backaspherics = Tuple{Int,S}.(backaspherics) end - surf = AcceleratedParametricSurface(ZernikeSurface(semidiameter + backdecenter_l + S(0.01), radius = backradius, conic = backconic, aspherics = backaspherics), nsamples, interface = opticinterface(S, insidematerial, nextmaterial, backsurfacereflectance, interfacemode)) + 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)) end extra_front = frontradius >= zero(S) || isinf(frontradius) ? zero(S) : abs(frontradius) - sqrt(frontradius^2 - semidiameter^2) diff --git a/test/Benchmarks/Benchmarks.jl b/test/Benchmarks/Benchmarks.jl index 014222958..d9713a451 100644 --- a/test/Benchmarks/Benchmarks.jl +++ b/test/Benchmarks/Benchmarks.jl @@ -56,13 +56,16 @@ simple_surface_benchmarks() = plane, rectangle, ellipse, triangle, sphere_outsid zernike_surface1() = surfaceintersection, (AcceleratedParametricSurface(TestData.zernikesurface1(), 30), perturbrayz()) zernike_surface2() = surfaceintersection, (AcceleratedParametricSurface(TestData.zernikesurface3(), 30), perturbrayz()) +evenaspheric_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.evenasphericsurface(), 30), perturbrayz()) +oddaspheric_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.oddasphericsurface(), 30), perturbrayz()) +oddevenaspheric_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.oddevenasphericsurface(), 30), perturbrayz()) qtype_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.qtypesurface1(), 30), perturbrayz()) bezier_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.beziersurface(), 30), bezierray()) chebyshev_surface1() = surfaceintersection, (AcceleratedParametricSurface(TestData.chebyshevsurface1(), 30), perturbrayz()) chebyshev_surface2() = surfaceintersection, (AcceleratedParametricSurface(TestData.chebyshevsurface2(), 30), perturbrayz()) gridsag_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.gridsagsurfacebicubic(), 30), perturbrayz()) -accel_surface_benchmarks() = zernike_surface1, zernike_surface2, qtype_surface, bezier_surface, chebyshev_surface1, chebyshev_surface2, gridsag_surface +accel_surface_benchmarks() = zernike_surface1, zernike_surface2, evenaspheric_surface, oddaspheric_surface, oddevenaspheric_surface, qtype_surface, bezier_surface, chebyshev_surface1, chebyshev_surface2, gridsag_surface all_benchmarks() = (simple_surface_benchmarks()..., accel_surface_benchmarks()..., system_benchmarks()...) diff --git a/test/TestData/lenses.jl b/test/TestData/lenses.jl index 31d4d2f30..b9341fdf7 100644 --- a/test/TestData/lenses.jl +++ b/test/TestData/lenses.jl @@ -9,15 +9,29 @@ function zernikelens() return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) end +function evenasphericlens() + topsurface = leaf(AcceleratedParametricSurface(AsphericSurface(9.5, radius = -50.0,aspherics = [(2, 0.0), (4, -0.001)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), Geometry.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)))) + return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) +end + +function oddasphericlens() + topsurface = leaf(AcceleratedParametricSurface(AsphericSurface(9.5, radius = -50.0,aspherics = [(3, 0.001), (5, -0.0001)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), Geometry.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)))) + return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) +end + function asphericlens() - topsurface = leaf(AcceleratedParametricSurface(ZernikeSurface(9.5, aspherics = [(2, 0.0), (4, -0.001)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) + topsurface = leaf(AcceleratedParametricSurface(AsphericSurface(9.5, aspherics = [(3, 0.0), (4, -0.001)]), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), 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)))) return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) end -function coniclensZ() - topsurface = leaf(AcceleratedParametricSurface(ZernikeSurface(9.5, radius = -15.0, conic = -5.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), translation(0.0, 0.0, 5.0)) +function coniclensA() + topsurface = leaf(AcceleratedParametricSurface(AsphericSurface(9.5, radius = -15.0, conic = -5.0), interface = FresnelInterface{Float64}(OpticSim.GlassCat.SCHOTT.N_BK7, OpticSim.GlassCat.Air)), 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)))) return (barrel ∩ topsurface ∩ botsurface)(Transform{Float64}(0.0, Float64(π), 0.0, 0.0, 0.0, -5.0)) diff --git a/test/TestData/surfaces.jl b/test/TestData/surfaces.jl index 4b24520a8..801c3dc13 100644 --- a/test/TestData/surfaces.jl +++ b/test/TestData/surfaces.jl @@ -83,7 +83,11 @@ zernikesurface2() = ZernikeSurface(1.5, radius = 5.0, conic = 1.0, zcoeff = [(7, zernikesurface3() = ZernikeSurface(9.5, zcoeff = [(1, 0.0), (4, -1.0), (7, 0.5)]) -conicsurface() = AcceleratedParametricSurface(ZernikeSurface(9.5, radius = -15.0, conic = -5.0)) +conicsurface() = AcceleratedParametricSurface(ZernikeSurface(9.5, radius = -15.0, conic = -5.0)) #this also tests the CONIC parameter of AsphericSurface + +oddasphericsurface() = OddAsphericSurface(9.5, 1.0/-50.0, 1.0, [0.0, 0.001, -0.0001]) +evenasphericsurface() = EvenAsphericSurface(9.5, 1.0/-50.0, 1.0, [0.001, -0.0001]) +oddevenasphericsurface() = OddEvenAsphericSurface(9.5, 1.0/-50.0, 1.0, [0.0, 0.001, -0.0001]) function simplesaggrid() points = map( diff --git a/test/TestData/systems.jl b/test/TestData/systems.jl index b8d37bfeb..d5fe9d767 100644 --- a/test/TestData/systems.jl +++ b/test/TestData/systems.jl @@ -148,7 +148,9 @@ end zernikesystem() = CSGOpticalSystem(LensAssembly(zernikelens()), Rectangle(40.0, 40.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) asphericsystem() = CSGOpticalSystem(LensAssembly(asphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) -conicsystemZ() = CSGOpticalSystem(LensAssembly(coniclensZ()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) +evenasphericsystem()=CSGOpticalSystem(LensAssembly(evenasphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) +oddasphericsystem()=CSGOpticalSystem(LensAssembly(oddasphericlens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) +conicsystemA() = CSGOpticalSystem(LensAssembly(coniclensA()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) conicsystemQ() = CSGOpticalSystem(LensAssembly(coniclensQ()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) qtypesystem() = CSGOpticalSystem(LensAssembly(qtypelens()), Rectangle(25.0, 25.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) chebyshevsystem() = CSGOpticalSystem(LensAssembly(chebyshevlens()), Rectangle(40.0, 40.0, [0.0, 0.0, 1.0], [0.0, 0.0, -67.8], interface = opaqueinterface())) diff --git a/test/testsets/Comparison.jl b/test/testsets/Comparison.jl index 00a858638..5eb392fd8 100644 --- a/test/testsets/Comparison.jl +++ b/test/testsets/Comparison.jl @@ -157,7 +157,7 @@ r4 = OpticalRay([0.0, -5.0, 1.0], [0.0, 0.08715574274765818, -0.9961946980917454], 1.0, λ) r5 = OpticalRay([-5.0, -5.0, 1.0], [0.08715574274765818, -0.01738599476176408, -0.9960429728140486], 1.0, λ) - a = TestData.conicsystemZ() + a = TestData.conicsystemA() res = trace(a, r1, test = true) @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE) @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE) @@ -195,6 +195,50 @@ @test isapprox(point(track[end]), [5.49905367197174, 5.66882664623822, 0.0], rtol = COMP_TOLERANCE) @test isapprox(pathlength(track[end]), 47.6825931025333, rtol = COMP_TOLERANCE) + a = TestData.evenasphericsystem() + res = trace(a, r1, test = true) + @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE) + res = trace(a, r2, test = true) + @test isapprox(point(res), [-1.1537225076569997, -1.153722507656998, -67.8], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(res), 74.08191763367685, rtol = COMP_TOLERANCE) + track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0) + res = trace(a, r3, test = true, trackrays = track) + @test (res === nothing) # TIR + @test isapprox(point(track[end]), [-4.684043128742365, -4.684043128742363, 0.0], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(track[end]), 45.610874251552275, rtol = COMP_TOLERANCE) + track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0) + res = trace(a, r4, test = true, trackrays = track) + @test isapprox(point(res), [0.0, 15.887441944041871, -67.8], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(res), 77.14042689028618, rtol = COMP_TOLERANCE) + track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0) + res = trace(a, r5, test = true, trackrays = track) + @test (res === nothing) # TIR + @test isapprox(point(track[end]), [5.095756154983988, 5.24166010736873, 0.0], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(track[end]), 45.618062468023375, rtol = COMP_TOLERANCE) + + a = TestData.oddasphericsystem() + res = trace(a, r1, test = true) + @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE) + res = trace(a, r2, test = true) + @test isapprox(point(res), [0.6256581890889898, 0.6256581890889894, -67.8], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(res), 73.97868203031574, rtol = COMP_TOLERANCE) + track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0) + res = trace(a, r3, test = true, trackrays = track) + @test (res === nothing) # TIR + @test isapprox(point(track[end]), [-6.097868437446996, -6.097868437446997, 0.0], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(track[end]), 48.34832912579036, rtol = COMP_TOLERANCE) + track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0) + res = trace(a, r4, test = true, trackrays = track) + @test isapprox(point(res), [0.0, 7.485281534548376, -67.8], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(res), 75.07125754375438, rtol = COMP_TOLERANCE) + track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0) + res = trace(a, r5, test = true, trackrays = track) + @test (res === nothing) # TIR + @test isapprox(point(track[end]), [6.197436823172213, 6.440993733499131, 0.0], rtol = COMP_TOLERANCE) + @test isapprox(pathlength(track[end]), 48.292537161700146, rtol = COMP_TOLERANCE) + a = TestData.zernikesystem() res = trace(a, r1, test = true) @test isapprox(point(res), [0.0, -8.9787010034042, -67.8], rtol = COMP_TOLERANCE) diff --git a/test/testsets/Intersection.jl b/test/testsets/Intersection.jl index 96682cb6c..f4dc12d9d 100644 --- a/test/testsets/Intersection.jl +++ b/test/testsets/Intersection.jl @@ -658,7 +658,6 @@ @testset "Zernike" begin @test_nowarn ZernikeSurface(1.5) @test_throws AssertionError ZernikeSurface(0.0) - @test_throws AssertionError ZernikeSurface(1.5, aspherics = [(1, 0.1)]) z1 = TestData.zernikesurface1() az1 = AcceleratedParametricSurface(z1, 20) From 72a6abf3f4a083dd9775c59a00d9f15aa2704262 Mon Sep 17 00:00:00 2001 From: BrianGun <36491850+BrianGun@users.noreply.github.com> Date: Fri, 12 Nov 2021 10:51:50 -0800 Subject: [PATCH 2/2] Return centroid from beamenergy function (#322) Fixes #321 --- .../ParaxialAnalysis/LensletAssembly.jl | 97 ++++++++++--------- test/testsets/ParaxialAnalysis.jl | 8 +- 2 files changed, 53 insertions(+), 52 deletions(-) diff --git a/src/Optical/ParaxialAnalysis/LensletAssembly.jl b/src/Optical/ParaxialAnalysis/LensletAssembly.jl index 1cd4a5aa2..bcec13c76 100644 --- a/src/Optical/ParaxialAnalysis/LensletAssembly.jl +++ b/src/Optical/ParaxialAnalysis/LensletAssembly.jl @@ -7,19 +7,19 @@ struct Display{T} <: Surface{T} surface::Rectangle{T} pixellattice::LatticeBasis{2,T} - pixelpitch::Tuple{T,T} #x,y pitch of pixels + pixelpitch::Tuple{T,T} # x,y pitch of pixels xresolution::T yresolution::T transform::Geometry.Transform{T} - function Display(xres::Int64,yres::Int64,xpitch::Unitful.Length,ypitch::Unitful.Length,transform::Geometry.Transform{T}) where{T<:Real} + function Display(xres::Int64, yres::Int64, xpitch::Unitful.Length, ypitch::Unitful.Length, transform::Geometry.Transform{T}) where {T <: Real} # convert xpitch,ypitch to unitless numbers representing μm - pitchx = ustrip(μm,xpitch) - pitchy = ustrip(μm,ypitch) - size = (xres*pitchx,yres*pitchy) - surface = Rectangle(size[1]/T(2),size[2]/T(2),SVector(T(0),T(0),T(1)),SVector(T(0),T(0),T(0))) #surface normal is +z axis - pixellattice = Repeat.rectangularlattice(size[2],size[1]) #lattice takes ypitch,xpitch in that order - return new{T}(surface,pixellattice,(pitchx,pitchy),xres,yres,transform) + pitchx = ustrip(μm, xpitch) + pitchy = ustrip(μm, ypitch) + size = (xres * pitchx, yres * pitchy) + surface = Rectangle(size[1] / T(2), size[2] / T(2), SVector(T(0), T(0), T(1)), SVector(T(0), T(0), T(0))) # surface normal is +z axis + pixellattice = Repeat.rectangularlattice(size[2], size[1]) # lattice takes ypitch,xpitch in that order + return new{T}(surface, pixellattice, (pitchx, pitchy), xres, yres, transform) end end @@ -27,43 +27,43 @@ surfaceintersection(a::Display) = surfaceintersection(a.surface) normal(a::Display) = normal(a.surface) interface(a::Display) = interface(a.surface) makemesh(a::Display) = makemesh(a.surface) -uv(a::Display,p::SVector{3,T}) where{T<:Real} = uv(a.surface,p) -onsurface(a::Display, point::SVector{3,T}) where {T<:Real} = onsurface(a.surface,point) -pixelposition(xindex::Int,yindex::Int) = transform*pixellattice[yindex,xindex] +uv(a::Display,p::SVector{3,T}) where {T <: Real} = uv(a.surface, p) +onsurface(a::Display, point::SVector{3,T}) where {T <: Real} = onsurface(a.surface, point) +pixelposition(xindex::Int,yindex::Int) = transform * pixellattice[yindex,xindex] struct LensletAssembly{T} lens::ParaxialLens{T} transform::Geometry.Transform{T} display::Display{T} worldtodisplay::Geometry.Transform{T} - LensletAssembly(lens::ParaxialLens{T},transform::Geometry.Transform{T},display::Display{T}) where{T} = new{T}(lens,transform,display,display.transform*transform) + LensletAssembly(lens::ParaxialLens{T}, transform::Geometry.Transform{T}, display::Display{T}) where {T} = new{T}(lens, transform, display, display.transform * transform) end lens(a::LensletAssembly) = a.lens display(a::LensletAssembly) = a.display transform(a::LensletAssembly) = a.transform -worldtolens(a::LensletAssembly, pt::SVector{3})= a.transform*pt -worldtolens(a::LensletAssembly,mat::SMatrix{3}) = a.transform*mat -lenstodisplay(a::LensletAssembly,pt::SVector{3,T}) where{T<:Real} = a.display.transform*pt -worldtodisplay(a::LensletAssembly,pt::SVector{3,T}) where{T<:Real} = a.worldtodisplay*pt +worldtolens(a::LensletAssembly, pt::SVector{3}) = a.transform * pt +worldtolens(a::LensletAssembly,mat::SMatrix{3}) = a.transform * mat +lenstodisplay(a::LensletAssembly,pt::SVector{3,T}) where {T <: Real} = a.display.transform * pt +worldtodisplay(a::LensletAssembly,pt::SVector{3,T}) where {T <: Real} = a.worldtodisplay * pt -vertices3d(vertices::SMatrix{2,N,T}) where{N,T<:Real} = vcat(vertices,zero(SMatrix{1,N,T})) +vertices3d(vertices::SMatrix{2,N,T}) where {N,T <: Real} = vcat(vertices, zero(SMatrix{1,N,T})) export vertices3d """Projects vertexpoints in world space onto the lens plane and converts them to two dimensional points represented in the local x,y coordinates of the lens coordinate frame. Used for projecting eye pupil onto lens plane.""" -function project(lenslet::LensletAssembly{T},displaypoint::SVector{3,T},vertexpoints::SMatrix{3,N,T}) where{T<:Real,N} +function project(lenslet::LensletAssembly{T}, displaypoint::SVector{3,T}, vertexpoints::SMatrix{3,N,T}) where {T <: Real,N} projectedpoints = MMatrix{2,N,T}(undef) - locvertices = worldtolens(lenslet,vertexpoints) #transform pupil vertices into local coordinate frame of lens - virtpoint = point(virtualpoint(lens(lenslet), displaypoint)) #compute virtual point corresponding to physical display point + locvertices = worldtolens(lenslet, vertexpoints) # transform pupil vertices into local coordinate frame of lens + virtpoint = point(virtualpoint(lens(lenslet), displaypoint)) # compute virtual point corresponding to physical display point for i in 1:N ppoint = locvertices[:,i] - vec = ppoint-virtpoint + vec = ppoint - virtpoint vecdist = ppoint[3] virtdist = virtpoint[3] - scale = abs(virtdist/(vecdist-virtdist)) - planepoint = scale*vec + virtpoint - projectedpoints[1:2,i] = SVector{2,T}(planepoint[1],planepoint[2]) #local lens coordinate frame has z axis aligned with the positive normal to the lens plane. + scale = abs(virtdist / (vecdist - virtdist)) + planepoint = scale * vec + virtpoint + projectedpoints[1:2,i] = SVector{2,T}(planepoint[1], planepoint[2]) # local lens coordinate frame has z axis aligned with the positive normal to the lens plane. end return SMatrix{2,N,T}(projectedpoints) @@ -82,51 +82,54 @@ function convertlazysets(verts::LazySets.VectorIterator) return SMatrix{3,M,T}(temp) end -"""Returns a number between 0 and 1 representing the ratio of the lens radiance to the pupil radiance. Assume lᵣ is the radiance transmitted through the lens from the display point. Some of this radiance, pᵣ, passes through the pupil. The beam energy is the ratio pᵣ/lᵣ.""" -function beamenergy(assy::LensletAssembly{T},displaypoint::AbstractVector{T},pupilpoints::SMatrix{3,N,T}) where{N,T<:Real} + +"""Returns a number between 0 and 1 representing the ratio of the lens radiance to the pupil radiance. Assume lᵣ is the radiance transmitted through the lens from the display point. Some of this radiance, pᵣ, passes through the pupil. The beam energy is the ratio pᵣ/lᵣ. Returns beam energy and the centroid of the intersection polygon to simplify ray casting.""" +function beamenergy(assy::LensletAssembly{T}, displaypoint::AbstractVector{T}, pupilpoints::SMatrix{3,N,T}) where {N,T <: Real} llens::ParaxialLens{T} = lens(assy) - virtpoint = point(virtualpoint(llens,displaypoint)) - projectedpoints = project(assy,displaypoint,pupilpoints) + virtpoint = point(virtualpoint(llens, displaypoint)) + projectedpoints = project(assy, displaypoint, pupilpoints) lensverts = vertices(llens) - rows,cols = size(lensverts) - temp = SMatrix{1,cols,T}(reshape(fill(0,cols),1,cols)) - lensverts3D = vcat(lensverts,temp) - beamlens = SphericalPolygon(lensverts3D,virtpoint,T(1)) + rows, cols = size(lensverts) + temp = SMatrix{1,cols,T}(reshape(fill(0, cols), 1, cols)) + lensverts3D = vcat(lensverts, temp) + beamlens = SphericalPolygon(lensverts3D, virtpoint, T(1)) projpoly = LazySets.VPolygon(projectedpoints) lenspoly = LazySets.VPolygon(lensverts) - intsct = projpoly ∩ lenspoly #this could be slow, especially multithreaded, because it will allocate. Lazysets.vertices returns Vector{SVector}, rather than SMatrix or SVector{SVector}. + intsct = projpoly ∩ lenspoly # this could be slow, especially multithreaded, because it will allocate. Lazysets.vertices returns Vector{SVector}, rather than SMatrix or SVector{SVector}. if isempty(intsct) return T(0) - else - beamintsct = SphericalPolygon(convertlazysets(LazySets.vertices(intsct)),virtpoint,T(1)) - return OpticSim.area(beamintsct)/OpticSim.area(beamlens) + else + verts = convertlazysets(LazySets.vertices(intsct)) + centroid = sum(eachcol(verts)) ./ size(verts)[2] + beamintsct = SphericalPolygon(verts, virtpoint, T(1)) + return OpticSim.area(beamintsct) / OpticSim.area(beamlens), centroid end end """Compute the bounding box in the display plane of the image of the worldpoly on the display plane. Pixels inside this area, conservatively, need to be turned on because some of their rays may pass through the worldpoly""" -function activebox(lenslet::LensletAssembly{T},worldpoly::SVector{N,SVector{3,T}}) where{T<:Real,N} +function activebox(lenslet::LensletAssembly{T}, worldpoly::SVector{N,SVector{3,T}}) where {T <: Real,N} maxx = typemin(T) maxy = typemin(T) minx = typemax(T) miny = typemax(T) lens = lens(lenslet) display = display(lenslet) - + for point in worldpoly - locpoint = transform(lenslet)*point + locpoint = transform(lenslet) * point for lenspoint in vertices(lens) - ray = Ray(locpoint,lenspoint-locpoint) - refractedray,_,_ = processintersection(interface(lens),lenspoint,ray, T(OpticSim.GlassCat.TEMP_REF), T(OpticSim.GlassCat.PRESSURE_REF)) - displaypoint = surfaceintersection(display,refractedray) - maxx = max(maxx,displaypoint[1]) - minx = min(minx,displaypoint[1]) - maxy = max(maxy,displaypoint[2]) - miny = min(miny,displaypoint[2]) + ray = Ray(locpoint, lenspoint - locpoint) + refractedray, _, _ = processintersection(interface(lens), lenspoint, ray, T(OpticSim.GlassCat.TEMP_REF), T(OpticSim.GlassCat.PRESSURE_REF)) + displaypoint = surfaceintersection(display, refractedray) + maxx = max(maxx, displaypoint[1]) + minx = min(minx, displaypoint[1]) + maxy = max(maxy, displaypoint[2]) + miny = min(miny, displaypoint[2]) end end - return minx,miny,maxx,maxy + return minx, miny, maxx, maxy end diff --git a/test/testsets/ParaxialAnalysis.jl b/test/testsets/ParaxialAnalysis.jl index 2fbfd18a1..01981c96a 100644 --- a/test/testsets/ParaxialAnalysis.jl +++ b/test/testsets/ParaxialAnalysis.jl @@ -38,11 +38,9 @@ using Unitful.DefaultSymbols displaypoint = SVector(0.0,0.0,-8.0) #pupil is placed so that only 1/4 of it (approximately) is illuminated by lens beam pupil = Rectangle(1.0,1.0,SVector(0.0,0.0,-1.0),SVector(2.0,2.0,40.0)) - - @test isapprox( - 1/16, - ParaxialAnalysis.beamenergy(lenslet,displaypoint,Geometry.vertices3d(pupil)), - atol = 1e-4) + energy,centroid = ParaxialAnalysis.beamenergy(lenslet,displaypoint,Geometry.vertices3d(pupil)) + @test isapprox(1/16, energy,atol = 1e-4) + @test isapprox([.75,.75,0.0],centroid) end @testset "SphericalPolygon" begin