diff --git a/Project.toml b/Project.toml index 6c10fc23e..699a0643f 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "OpticSim" uuid = "24114763-4efb-45e7-af0e-cde916beb153" authors = ["Brian Guenter", "Charlie Hewitt", "Ran Gal", "Alfred Wong"] repository = "https://github.com/microsoft/OpticSim.jl" -version = "0.4.1" +version = "0.4.2" [deps] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" diff --git a/deps/generate.jl b/deps/generate.jl index f364f4229..68125f1bb 100644 --- a/deps/generate.jl +++ b/deps/generate.jl @@ -198,18 +198,15 @@ function catalog_to_modstring(start_id::Integer, catalogname::AbstractString, ca modstrings = [ "module $catalogname", "using ..GlassCat: Glass, GlassID, AGF", - # "using StaticArrays: SVector", "export $(join(keys(catalog), ", "))", "" ] for (glassname, glassinfo) in catalog - # skip docstrings for CI builds to avoid 'missing docstring' warnings in makedocs - docstring = isCI ? "" : glassinfo_to_docstring(glassinfo, id, catalogname, glassname) argstring = glassinfo_to_argstring(glassinfo, id) - append!(modstrings, [docstring, "const $glassname = Glass($argstring)", ""]) + push!(modstrings, "const $glassname = Glass($argstring)") id += 1 end - append!(modstrings, ["end #module", ""]) # last "" is for \n at EOF + append!(modstrings, ["end # module", ""]) # last "" is for \n at EOF return id, join(modstrings, "\n") end diff --git a/docs/make.jl b/docs/make.jl index e155d7b28..2d3817928 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -21,7 +21,7 @@ # SOFTWARE using Documenter -using OpticSim, OpticSim.Geometry, OpticSim.Emitters +using OpticSim makedocs( sitename = "OpticSim.jl", @@ -49,6 +49,7 @@ makedocs( "Glass Functions" => "glasscat.md", "Optimization" => "optimization.md", "Cloud Execution" => "cloud.md", + "Notebook utilities" => "notebooksutils.md", "Reference" => "ref.md", "Roadmap" => "roadmap.md" ], diff --git a/docs/src/basic_types.md b/docs/src/basic_types.md index b1aa1f466..6187dde54 100644 --- a/docs/src/basic_types.md +++ b/docs/src/basic_types.md @@ -11,12 +11,12 @@ using OpticSim, OpticSim.Geometry Representing a 3D vector. ```@docs -Geometry.Vec3 -Geometry.unitX3 -Geometry.unitY3 -Geometry.unitZ3 -Geometry.zero3 -Geometry.one3 +OpticSim.Geometry.Vec3 +OpticSim.Geometry.unitX3 +OpticSim.Geometry.unitY3 +OpticSim.Geometry.unitZ3 +OpticSim.Geometry.zero3 +OpticSim.Geometry.one3 ``` ## Vec4 @@ -24,13 +24,13 @@ Geometry.one3 Representing a 4D vector ```@docs -Geometry.Vec4 -Geometry.unitX4 -Geometry.unitY4 -Geometry.unitZ4 -Geometry.unitW4 -Geometry.zero4 -Geometry.one4 +OpticSim.Geometry.Vec4 +OpticSim.Geometry.unitX4 +OpticSim.Geometry.unitY4 +OpticSim.Geometry.unitZ4 +OpticSim.Geometry.unitW4 +OpticSim.Geometry.zero4 +OpticSim.Geometry.one4 ``` ## Transform @@ -39,17 +39,16 @@ Representing a general 3D transform (4x4 matrix). Currently only used as a rigid Transforms are used to position surfaces within the CSG tree, position emitters in 3D, etc. ```@docs -Geometry.Transform -Geometry.identitytransform -Geometry.rotationX -Geometry.rotationY -Geometry.rotationZ -Geometry.rotation -Geometry.rotationd -Geometry.rotate -Geometry.translation -Geometry.rotmat -Geometry.rotmatd -Geometry.rotmatbetween +OpticSim.Geometry.Transform +OpticSim.Geometry.identitytransform +OpticSim.Geometry.rotationX +OpticSim.Geometry.rotationY +OpticSim.Geometry.rotationZ +OpticSim.Geometry.rotation +OpticSim.Geometry.rotationd +OpticSim.Geometry.rotate +OpticSim.Geometry.translation +OpticSim.Geometry.rotmat +OpticSim.Geometry.rotmatd +OpticSim.Geometry.rotmatbetween ``` - diff --git a/docs/src/csg.md b/docs/src/csg.md index b8170cd7c..1d3c1aba8 100644 --- a/docs/src/csg.md +++ b/docs/src/csg.md @@ -12,7 +12,7 @@ The code for this in our system would look this this: ```@example using OpticSim # hide cyl = Cylinder(0.7) -cyl_cross = csgunion(csgunion(leaf(cyl), leaf(cyl, rotationd(90, 0, 0))), leaf(cyl, rotationd(0, 90, 0))) +cyl_cross = csgunion(csgunion(leaf(cyl), leaf(cyl, Geometry.rotationd(90, 0, 0))), leaf(cyl, Geometry.rotationd(0, 90, 0))) cube = Cuboid(1.0, 1.0, 1.0) sph = Sphere(1.3) diff --git a/docs/src/emitters_new.md b/docs/src/emitters_new.md index ff43fc982..509104795 100644 --- a/docs/src/emitters_new.md +++ b/docs/src/emitters_new.md @@ -21,7 +21,7 @@ The [`OpticSim`](index.html) package comes with various implementations of each * [`Emitters.Origins.Point`](@ref) - a single point * [`Emitters.Origins.RectUniform`](@ref) - a uniformly sampled rectangle with user defined number of samples * [`Emitters.Origins.RectGrid`](@ref) - a rectangle sampled in a grid fashion - * [`Emitters.Origins.Hexapolar`](@ref) - a circle (or an ellipse) sampled in an hexapolar fasion (rings) + * [`Emitters.Origins.Hexapolar`](@ref) - a circle (or an ellipse) sampled in an hexapolar fashion (rings) - Rays Directions Distribution - the interface **length** returns the number of samples, and **generate** returns the n'th sample. * [`Emitters.Directions.Constant`](@ref) * [`Emitters.Directions.RectGrid`](@ref) @@ -30,7 +30,7 @@ The [`OpticSim`](index.html) package comes with various implementations of each ## [Examples of Basic Emitters](@id basic_emitters) -**Note**: All the examples on this page assumes the followin statement was executed: +**Note**: All of the examples on this page assume that the following statement was executed: ```@example using OpticSim, OpticSim.Geometry, OpticSim.Emitters @@ -209,7 +209,7 @@ P = AngularPower.Lambertian() O = Origins.RectGrid(1.0, 1.0, 3, 3) D = Directions.HexapolarCone(deg2rad(5.0), 3) -# construct the srouce. in this example a "pixel" source will contain only one source as we are simulating a "b/w" display. +# construct the source. in this example a "pixel" source will contain only one source as we are simulating a "b/w" display. # for RGB displays we can combine 3 sources to simulate "a pixel". Tr = Transform(Vec3(0.5, 0.5, 0.0)) source1 = Sources.Source(Tr, S, O, D, P) @@ -232,7 +232,7 @@ end Tr = Transform(Vec3(0.0, 0.0, 0.0)) my_display = Sources.CompositeSource(Tr, pixels) -Vis.draw(my_display) # render the display - nothing bu the origins primities +Vis.draw(my_display) # render the display - nothing but the origins primitives rays = AbstractArray{OpticalRay{Float64, 3}}(collect(my_display)) # collect the rays generated by the display Vis.draw!(rays) # render the rays @@ -261,6 +261,7 @@ A display example using composite sources. Emitters.Spectrum.Uniform Emitters.Spectrum.DeltaFunction Emitters.Spectrum.Measured +Emitters.Spectrum.spectrumpower ``` ## [Angular Power Distribution](@id angular_power_distribution) diff --git a/docs/src/examples.md b/docs/src/examples.md index df15345dd..21db87df5 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -37,51 +37,6 @@ nothing # hide ![Cooke triplet visualization](assets/cooke.png) -## Zoom Lens - -```@example -using OpticSim - -using DataFrames -function zoom_lens(pos = 1) - if pos == 0 - stop = 2.89 - zoom = 9.48 - dist = 4.46970613 - elseif pos == 1 - stop = 3.99 - zoom = 4.48 - dist = 21.21 - else - stop = 4.90 - zoom = 2.00 - dist = 43.81 - end - return AxisymmetricOpticalSystem{Float64}( - DataFrame(Surface = [:Object, :Stop, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, :Image], - Radius = [Inf64, Inf64, -1.6202203499676E+01, -4.8875855327468E+01, 1.5666614444619E+01, -4.2955326460481E+01, 1.0869565217391E+02, 2.3623907394283E+01, -1.6059097478722E+01, -4.2553191489362E+02, -3.5435861091425E+01, -1.4146272457208E+01, -2.5125628140704E+02, -2.2502250225023E+01, -1.0583130489999E+01, -4.4444444444444E+01, Inf64], - Aspherics = [missing, missing, missing, missing, missing, [(4, 1.03860000000E-04), (6, 1.42090000000E-07), (8, -8.84950000000E-09), (10, 1.24770000000E-10), (12, -1.03670000000E-12), (14, 3.65560000000E-15)], missing, missing, [(4, 4.27210000000E-05), (6, 1.24840000000E-07), (8, 9.70790000000E-09), (10, -1.84440000000E-10), (12, 1.86440000000E-12), (14, -7.79750000000E-15)], [(4, 1.13390000000E-04), (6, 4.81650000000E-07), (8, 1.87780000000E-08), (10, -5.75710000000E-10), (12, 8.99940000000E-12), (14, -4.67680000000E-14)], missing, missing, missing, missing, missing, missing, missing], - Thickness = [Inf64, 0.0, 5.18, 0.10, 4.40, 0.16, 1.0, 4.96, zoom, 4.04, 1.35, 1.0, 2.80, 3.0, 1.22, dist, missing], - Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.Air, OpticSim.GlassCat.OHARA.S_LAH66, OpticSim.GlassCat.Air, OpticSim.GlassCat.NIKON.LLF6, OpticSim.GlassCat.Air, OpticSim.GlassCat.OHARA.S_TIH6, OpticSim.GlassCat.OHARA.S_FSL5, OpticSim.GlassCat.Air, OpticSim.GlassCat.OHARA.S_FSL5, OpticSim.GlassCat.Air, OpticSim.GlassCat.OHARA.S_LAL8, OpticSim.GlassCat.OHARA.S_FSL5, OpticSim.GlassCat.Air, OpticSim.GlassCat.OHARA.S_LAH66, OpticSim.GlassCat.Air, missing], - SemiDiameter = [Inf64, stop, 3.85433218451, 3.85433218451, 4.36304692871, 4.36304692871, 4.72505505439, 4.72505505439, 4.72505505439, 4.45240784026, 4.45240784026, 4.50974054117, 4.50974054117, 4.50974054117, 4.76271114409, 4.76271114409, 15.0])) -end -@show zoom_lens(0) -Vis.drawtracerays(zoom_lens(0), test = true, trackallrays = true, numdivisions = 50, resolution = (1200, 600)) -Vis.make2dy() # hide -Vis.save("assets/zoom0.png") # hide -Vis.drawtracerays(zoom_lens(1), test = true, trackallrays = true, numdivisions = 50, resolution = (1200, 600)) -Vis.make2dy() # hide -Vis.save("assets/zoom1.png") # hide -Vis.drawtracerays(zoom_lens(2), test = true, trackallrays = true, numdivisions = 50, resolution = (1200, 600)) -Vis.make2dy() # hide -Vis.save("assets/zoom2.png") # hide -nothing # hide -``` - -![Zoom position 1 visualization](assets/zoom0.png) -![Zoom position 2 visualization](assets/zoom1.png) -![Zoom position 3 visualization](assets/zoom2.png) - ## Schmidt Cassegrain Telescope ```@example diff --git a/docs/src/notebooksutils.md b/docs/src/notebooksutils.md new file mode 100644 index 000000000..e9a78f4df --- /dev/null +++ b/docs/src/notebooksutils.md @@ -0,0 +1,14 @@ +# Notebook utilities + +```@meta +CurrentModule = NotebooksUtils +``` + +[TODO] + +```@docs +run_sample +SetBackend +run +InitNotebook +``` diff --git a/docs/src/ref.md b/docs/src/ref.md index 85d4d79b2..028d09221 100644 --- a/docs/src/ref.md +++ b/docs/src/ref.md @@ -14,6 +14,11 @@ Pages = ["ref.md"] Modules = [OpticSim] ``` +## Geometry +```@autodocs +Modules = [OpticSim.Geometry] +``` + ## Zernike ```@autodocs diff --git a/src/Geometry/Geometry.jl b/src/Geometry/Geometry.jl index 1d7ac5f53..332eafbc6 100644 --- a/src/Geometry/Geometry.jl +++ b/src/Geometry/Geometry.jl @@ -26,6 +26,9 @@ module Geometry include("Transform.jl") end # module Geometry +export Geometry + + using .Geometry include("Ray.jl") diff --git a/src/GlassCat/data/CARGILLE.jl b/src/GlassCat/data/CARGILLE.jl index 31bfbed72..bf500c49d 100644 --- a/src/GlassCat/data/CARGILLE.jl +++ b/src/GlassCat/data/CARGILLE.jl @@ -24,49 +24,10 @@ module CARGILLE using ..GlassCat: Glass, GlassID, OTHER -using StaticArrays: SVector -""" CARGILLE.OG0608 -``` -ID: OTHER:1 -RI @ 587nm: 1.457518 -Abbe Number: 57.18978 -ΔPgF: 0.008 -TCE (÷1e-6): 800.0 -Density: 0.878g/m³ -Valid wavelengths: 0.32μm to 1.55μm -Reference Temp: 25°C -``` -""" -const OG0608 = Glass(GlassID(OTHER, 1), -2, 1.4451400, 0.0043176, -1.80659e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.0009083144750540808, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.008, -1.0, -1.0, 800.0, -1.0, 0, -1.0, [SVector(0.32, 0.03, 10.0), SVector(0.365, 0.16, 100.0), SVector(0.4047, 0.40, 100.0), SVector(0.480, 0.71, 100.0), SVector(0.4861, 0.72, 100.0), SVector(0.5461, 0.80, 100.0), SVector(0.5893, 0.90, 100.0), SVector(0.6328, 0.92, 100.0), SVector(0.6439, 0.95, 100.0), SVector(0.6563, 0.96, 100.0), SVector(0.6943, 0.99, 100.0), SVector(0.840, 0.99, 100.0), SVector(0.10648, 0.74, 100.0), SVector(0.1300, 0.39, 100.0), SVector(0.1550, 0.16, 100.0)], 1.457518, -1.0, -1.0, 0, 57.18978, 0, 0.878, -1) - -""" CARGILLE.OG0607 -``` -ID: OTHER:2 -RI @ 587nm: 1.457587 -Abbe Number: 57.19833 -ΔPgF: 0.008 -TCE (÷1e-6): 700.0 -Density: 0.878g/m³ -Valid wavelengths: 0.32μm to 1.55μm -Reference Temp: 25°C -``` -""" -const OG0607 = Glass(GlassID(OTHER, 2), -2, 1.44503, 0.0044096, -2.85878e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.0009083144750540808, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.008, -1.0, -1.0, 700.0, -1.0, 0, -1.0, [SVector(0.32, 0.15, 10.0), SVector(0.365, 0.12, 100.0), SVector(0.4047, 0.42, 100.0), SVector(0.480, 0.78, 100.0), SVector(0.4861, 0.79, 100.0), SVector(0.5461, 0.86, 100.0), SVector(0.5893, 0.90, 100.0), SVector(0.6328, 0.92, 100.0), SVector(0.6439, 0.90, 100.0), SVector(0.6563, 0.92, 100.0), SVector(0.6943, 0.98, 100.0), SVector(0.840, 0.99, 100.0), SVector(0.10648, 0.61, 100.0), SVector(0.1300, 0.39, 100.0), SVector(0.1550, 0.11, 100.0)], 1.457587, -1.0, -1.0, 0, 57.19833, 0, 0.878, -1) - -""" CARGILLE.OG081160 -``` -ID: OTHER:3 -RI @ 587nm: 1.515549 -Abbe Number: 36.82493 -ΔPgF: 0.014 -TCE (÷1e-6): 700.0 -Density: 1.11g/m³ -Valid wavelengths: 0.32μm to 1.55μm -Reference Temp: 25°C -``` -""" -const OG081160 = Glass(GlassID(OTHER, 3), -2, 1.49614, 0.00692199, -8.07052e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.000885983052189022, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.014, -1.0, -1.0, 700.0, -1.0, 0, -1.0, [SVector(0.32, 0.04, 100.0), SVector(0.365, 0.13, 100.0), SVector(0.4047, 0.26, 100.0), SVector(0.480, 0.48, 100.0), SVector(0.4861, 0.49, 100.0), SVector(0.5461, 0.60, 100.0), SVector(0.5893, 0.68, 100.0), SVector(0.6328, 0.71, 100.0), SVector(0.6439, 0.73, 100.0), SVector(0.6563, 0.74, 100.0), SVector(0.6943, 0.76, 100.0), SVector(0.840, 0.83, 100.0), SVector(0.10648, 0.86, 100.0), SVector(0.1300, 0.89, 100.0), SVector(0.1550, 0.90, 100.0)], 1.515549, -1.0, -1.0, 0, 36.82493, 0, 1.11, -1) +const OG0608 = Glass(GlassID(OTHER, 1), -2, 1.4451400, 0.0043176, -1.80659e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.0009083144750540808, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.008, -1.0, -1.0, 800.0, -1.0, 0, -1.0, [(0.32, 0.03, 10.0), (0.365, 0.16, 100.0), (0.4047, 0.40, 100.0), (0.480, 0.71, 100.0), (0.4861, 0.72, 100.0), (0.5461, 0.80, 100.0), (0.5893, 0.90, 100.0), (0.6328, 0.92, 100.0), (0.6439, 0.95, 100.0), (0.6563, 0.96, 100.0), (0.6943, 0.99, 100.0), (0.840, 0.99, 100.0), (0.10648, 0.74, 100.0), (0.1300, 0.39, 100.0), (0.1550, 0.16, 100.0)], 1.457518, -1.0, -1.0, 0, 57.18978, 0, 0.878, -1) +const OG0607 = Glass(GlassID(OTHER, 2), -2, 1.44503, 0.0044096, -2.85878e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.0009083144750540808, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.008, -1.0, -1.0, 700.0, -1.0, 0, -1.0, [(0.32, 0.15, 10.0), (0.365, 0.12, 100.0), (0.4047, 0.42, 100.0), (0.480, 0.78, 100.0), (0.4861, 0.79, 100.0), (0.5461, 0.86, 100.0), (0.5893, 0.90, 100.0), (0.6328, 0.92, 100.0), (0.6439, 0.90, 100.0), (0.6563, 0.92, 100.0), (0.6943, 0.98, 100.0), (0.840, 0.99, 100.0), (0.10648, 0.61, 100.0), (0.1300, 0.39, 100.0), (0.1550, 0.11, 100.0)], 1.457587, -1.0, -1.0, 0, 57.19833, 0, 0.878, -1) +const OG081160 = Glass(GlassID(OTHER, 3), -2, 1.49614, 0.00692199, -8.07052e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.000885983052189022, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.014, -1.0, -1.0, 700.0, -1.0, 0, -1.0, [(0.32, 0.04, 100.0), (0.365, 0.13, 100.0), (0.4047, 0.26, 100.0), (0.480, 0.48, 100.0), (0.4861, 0.49, 100.0), (0.5461, 0.60, 100.0), (0.5893, 0.68, 100.0), (0.6328, 0.71, 100.0), (0.6439, 0.73, 100.0), (0.6563, 0.74, 100.0), (0.6943, 0.76, 100.0), (0.840, 0.83, 100.0), (0.10648, 0.86, 100.0), (0.1300, 0.89, 100.0), (0.1550, 0.90, 100.0)], 1.515549, -1.0, -1.0, 0, 36.82493, 0, 1.11, -1) end diff --git a/src/NotebooksUtils/NotebooksUtils.jl b/src/NotebooksUtils/NotebooksUtils.jl index a7a3e6652..5759e8035 100644 --- a/src/NotebooksUtils/NotebooksUtils.jl +++ b/src/NotebooksUtils/NotebooksUtils.jl @@ -71,7 +71,7 @@ function run(; port=nothing, path=nothing, sysimage_file=nothing, auto_detect_sy local file_path = $path local sysimage_file_path = $sysimage_file - if (file_path==nothing) + if (file_path === nothing) @info "Launching Pluto" else file_path = fix_path(file_path) diff --git a/src/OpticSim.jl b/src/OpticSim.jl index 924799ad1..10c53387d 100644 --- a/src/OpticSim.jl +++ b/src/OpticSim.jl @@ -44,12 +44,7 @@ include("constants.jl") include("utilities.jl") include("Geometry/Geometry.jl") include("Optical/Optical.jl") -include("Visualization/Visualization.jl") - -# TODO: RG: the following include is seperated from the Optical include due to the new approach of putting the visualazation -# part with the logic. We need to come up with a better approch to how and where visualization code should reside. -include("Optical/Emitters.jl") # defines the Emitters module - +include("Vis/Vis.jl") include("Examples/Examples.jl") include("Optimization/Optimization.jl") include("Cloud/Cloud.jl") diff --git a/src/Optical/Emitter2.jl b/src/Optical/Emitter2.jl deleted file mode 100644 index 4d66aa2e3..000000000 --- a/src/Optical/Emitter2.jl +++ /dev/null @@ -1,282 +0,0 @@ -# MIT License - -# Copyright (c) Microsoft Corporation. - -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: - -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. - -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE - -#= -Emitters are defined by Pixels and SpatialLayouts. A emitter has a spectrum, and an optical power distribution over the hemisphere. These are -intrinsic physical properties of the emitter. -=# -#Ray distribution types -using Distributions - -""" -Each AbstractSpectrum type defines a spectrumSample function which returns a uniformly sampled point from the spectrum of the light source and the power at that wavelength. -""" -abstract type AbstractSpectrum{T<:Real} end -""" Every ray direction distribution must implement the functions direction(a::Source)""" -abstract type AbstractRayDirectionDistribution{T<:Real} end -abstract type AbstractAngularPowerDistribution{T<:Real} end -abstract type AbstractRayOriginDistribution{T<:Real} end - -""" -Create an emitter by choosing a type for the spectral distribution, ray origin distribution, ray direction distribution, and angular light power distribution. -For example, to create an emitter with uniform spectrum, rays originating on a regular grid, ray directions distributed around a vector v up to a maximum angle θmax, with Lambertian power distribution call the Source construction with these types and arguments: - -a = Source{UniformSpectrum,GridOrigin, ConeDistribution, Lambertian}(UniformSpectrum(),GridOrigin(),ConeDistribution(θmax),Lambertian()) -""" -struct Source{T<:Real,S<:AbstractSpectrum{T},O<:AbstractRayOriginDistribution{T},D<:AbstractRayDirectionDistribution{T},P<:AbstractAngularPowerDistribution{T}} - spectrum::S - rayorigin::O - raydirection::D - power::P - - Source(spectrum::S, origins::O, directions::D, power::P, ::Type{T} = Float64) where {S<:AbstractSpectrum,O<:AbstractRayOriginDistribution,D<:AbstractRayDirectionDistribution,P<:AbstractAngularPowerDistribution,T<:Real} = new{T,S,O,D,P}(spectrum, origins, directions, power) -end -export Source - -const UNIFORMSHORT = 0.450 #um -const UNIFORMLONG = 0.680 #um - -""" -flat spectrum from 450nm to 680nm -""" -struct UniformSpectrum{T} <: AbstractSpectrum{T} end -export UniformSpectrum -"""returns a tuple (power,wavelength) """ -spectrumsample(::Source{T,UniformSpectrum{T},O,D,P}) where {O,D,P,T} = (T(1), rand(Distributions.Uniform(UNIFORMSHORT, UNIFORMLONG))) - -struct DeltaFunctionSpectrum{T<:Real} <: AbstractSpectrum{T} - λ::T -end -export DeltaFunctionSpectrum -spectrumsample(a::Source{T,DeltaFunctionSpectrum{T},O,D,P}) where {O,D,P,T} = (T(1), a.spectrum.λ) - -""" -Use measured spectrum to compute emitter power. Create spectrum by reading CSV files. - -Evaluate spectrum at arbitrary wavelength with [`spectrumpower`](@ref) -""" -struct MeasuredSpectrum{T<:Real} <: AbstractSpectrum{T} - lowwavelength::Int - highwavelength::Int - wavelengthstep::Int - powersamples::Vector{T} - - function MeasuredSpectrum(samples::DataFrame) - colnames = names(samples) - @assert "Wavelength" in colnames - @assert "Power" in colnames - wavelengths = samples[:Wavelength] - @assert eltype(wavelengths) <: Integer - - power::Vector{T} where {T<:Real} = samples[:Power] # no missing values allowed and must be real numbers - maxpower = maximum(power) - T = eltype(power) - p = sortperm(wavelengths) - wavelengths = wavelengths[p] - power = power[p] ./ maxpower #make sure power values have the same sorted permutation as wavelengths normalized with maximum power equal to 1 - λmin = wavelengths[p[1]] - λmax = wavelengths[p[end]] - step = wavelengths[2] - wavelengths[1] - - for i in 3:length(wavelengths) - @assert wavelengths[i] - wavelengths[i - 1] == step # no missing values allowed and step size must be consistent - end - - return new{T}(λmin, λmax, step, power) - end -end -export MeasuredSpectrum - -"""expects wavelength in nm not um""" -function spectrumpower(spectrum::MeasuredSpectrum{T}, λ::T)::Union{Nothing,T} where {T<:Real} - λ = λ #convert from um to nm - if λ < spectrum.lowwavelength || λ > spectrum.highwavelength - return nothing - end - - lowindex = floor(Int, (λ - spectrum.lowwavelength)) ÷ spectrum.wavelengthstep + 1 - - if lowindex == length(spectrum.powersamples) - return T(spectrum.powersamples[end]) - else - highindex = lowindex + 1 - α = mod(λ, spectrum.wavelengthstep) / spectrum.wavelengthstep - - return T((1 - α) * spectrum.powersamples[lowindex] + α * spectrum.powersamples[highindex]) - end -end - -# Every spectrum type must define a function spectrumsample which randomly selects a wavelength and returns the power at that wavelength. For delta function spectra the same wavelength will be returned every time. - -function spectrumsample(a::Source{T,MeasuredSpectrum{T},O,D,P}) where {O,D,P,T} - spectrum = a.spectrum::MeasuredSpectrum{T} - λ = rand(Distributions.Uniform{T}(T(spectrum.lowwavelength), T(spectrum.highwavelength))) - power = spectrumpower(spectrum, λ) - if power === nothing #added this condition because the compiler was allocating if just returned values directly. - return (nothing, nothing) - else - return (power, λ / T(1000)) #convert from nm to um - end -end - -"""stores the state of ray generation. Keeps track of ray number so you can have a PointOrigin which emits multiple rays. This can be used in the future to do phase summations on the image plane, since only point sources are perfectly spatially coherent.""" -struct RayState - originnumber::Int - directionnumber::Int - - RayState() = new(1, 1) - RayState(originnum::Int, directionnum::Int) = new(originnum, directionnum) -end - -struct PointOrigin{T<:Real} <: AbstractRayOriginDistribution{T} - position::SVector{3,T} -end -export PointOrigin - -numsamples(a::PointOrigin{T}) where {T} = 1 - -origin(a::Source{T,S,PointOrigin{T},D,P}, ::Int) where {S,D,P,T} = a.rayorigin.position - -struct UniformRandomOrigin{T} <: AbstractRayOriginDistribution{T} - width::T - height::T - numsamples::Int -end -export UniformRandomOrigin - -numsamples(a::UniformRandomOrigin{T}) where {T} = a.numsamples - -function origin(a::Source{T,S,UniformRandomOrigin{T},D,P}, ::Int) where {S,D,P,T} - origin = a.origin - x = rand(Distributions.Uniform(-origin.width / 2.0, origin.width / 2.0)) - y = rand(Distributions.Uniform(-origin.height / 2.0, origin.height / 2.0)) - z = T(0) -end - -struct GridOrigin{T} <: AbstractRayOriginDistribution{T} - width::T - height::T - usamples::Int - vsamples::Int - ustep::T - vstep::T - - GridOrigin(width::T, height::T, usamples::Int, vsamples::Int) where {T<:Real} = new{T}(width, height, usamples, vsamples, width / (usamples - 1), height / (vsamples - 1)) -end -export GridOrigin - -numsamples(a::GridOrigin{T}) where {T<:Real} = a.usamples * a.vsamples - -function origin(a::Source{T,S,GridOrigin{T},D,P}, originnumber::Int) where {S,D,P,T} - @assert originnumber <= numsamples(a.rayorigin) - origin = a.rayorigin - - u = originnumber ÷ origin.vsamples - v = mod(originnumber, origin.usamples) - return SVector{3,T}(T(-origin.width / 2.0 + origin.ustep * u), T(-origin.height / 2.0 + origin.vstep * v), T(0)) -end - -struct HexapolarOrigin{T} <: AbstractRayOriginDistribution{T} - position::SVector{3,T} - direction::SVector{3,T} - uvec::SVector{3,T} - vvec::SVector{3,T} - halfsizeu::T - halfsizev::T - nrings::Int - function HexapolarOrigin{T}(nrings::Int, halfsizeu::T, halfsizev::T; position::SVector{3,T} = SVector{3,T}(0.0, 0.0, 0.0), direction::SVector{3,T} = SVector{3,T}(0.0, 0.0, -1.0), rotationvec::SVector{3,T} = SVector{3,T}(0.0, 1.0, 0.0)) where {T<:Real} - new{T}(position, getuvvecs(direction, rotationvec)..., halfsizeu, halfsizev, nrings) - end -end - -struct ConstantRayDirection{T} <: AbstractRayDirectionDistribution{T} - direction::SVector{3,T} -end -export ConstantRayDirection - -numsamples(a::ConstantRayDirection{T}) where {T} = 1 - -direction(a::Source{T,S,O,ConstantRayDirection{T},P}, directionnumber::Int) where {S,O,P,T} = a.raydirection.direction - -struct ConeDistribution{T} <: AbstractRayDirectionDistribution{T} - θmax::T - numsamples::Int -end -export ConeDistribution - -numsamples(a::ConeDistribution{T}) where {T} = a.numsamples - -""" -Generates a unit vector pointing somewhere within the cone with half angle `θmax` around `direction` which is the normal to the emitter surface. -""" -function direction(a::Source{T,S,O,ConeDistribution{T},P}, directionnumber::Int) where {S,O,P,T<:Real} - direction = SVector{3,T}(T(0), T(0), T(1)) - θmax = a.raydirection.θmax - uvec = SVector{3,T}(T(1), T(0), T(0)) - vvec = SVector{3,T}(T(0), T(1), T(0)) - - ϕ = rand(T) * 2π - θ = NaNsafeacos(one(T) + rand(T) * (cos(θmax) - 1)) - return normalize(sin(θ) * (cos(ϕ) * uvec + sin(ϕ) * vvec) + cos(θ) * direction) -end - - - -Base.length(a::Source{T,S,O,D,P}) where {T,S,O,D,P} = numsamples(a.rayorigin) * numsamples(a.raydirection) - -"""raysample is the function that can couple origin and direction generation if necessary. The default function couples them in a simple way but more complex coupling should be possible. For each origin numsamples(a.raydirection) direction samples are taken with identical origin. Then the origin number is incremented. This repeats till all rays have been generated. The origin and direction functions receive an integer indicating the origin or direction number so regular patterns such as rectangular and hexapolar grids can be generated properly.""" -function raysample(a::Source{T,S,O,D,P}, r::RayState) where {T,S,O,D,P} - newray = Ray(origin(a, r.originnumber), direction(a, r.directionnumber)) - dirnum = r.directionnumber + 1 - if dirnum > numsamples(a.raydirection) - return (newray, RayState(r.originnumber + 1, 1)) #reset counter for direction number and move on to next origin - else - return (newray, RayState(r.originnumber, dirnum)) - end - -end - -struct Lambertian{T} <: AbstractAngularPowerDistribution{T} end -export Lambertian - -directionpower(::Source{T,S,O,D,Lambertian{T}}, ::SVector{3,T}) where {S,O,D,T} = T(1) - -function Base.iterate(a::Source)::Tuple{OpticalRay,RayState} - state = RayState() - generateray(a, state) -end - -function Base.iterate(a::Source, b::RayState)::Union{Nothing,Tuple{OpticalRay,RayState}} - if b.originnumber > numsamples(a.rayorigin) - return nothing - else - return generateray(a, b) - end -end - -# """generates an optical ray with origin and spectrum distribution dictated by OriginDistribution and Spectrum type. Ray is generated in the emitter local coordinate frame.""" -function generateray(a::Source, state::RayState, sourcenumber = 0) - spectralpower, λ = spectrumsample(a) - ray, newstate = raysample(a, state) - optray = OpticalRay(ray, spectralpower * directionpower(a, direction(ray)), λ, sourcenum = sourcenumber) - return (optray, newstate) -end diff --git a/src/Optical/Emitters.jl b/src/Optical/Emitters.jl deleted file mode 100644 index 941762643..000000000 --- a/src/Optical/Emitters.jl +++ /dev/null @@ -1,1108 +0,0 @@ -# MIT License - -# Copyright (c) Microsoft Corporation. - -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: - -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. - -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE - -#= -Emitters are defined by Pixels and SpatialLayouts. An emitter has a spectrum, and an optical power distribution over the hemisphere. -These are intrinsic physical properties of the emitter. -=# - -module Emitters - -using OpticSim, OpticSim.Geometry -using LinearAlgebra - -struct MissingImplementationException <: Exception - msg::String -end - -# defining name placeholders to override in nested modules -generate() = 0 -export generate -visual_size() = 0 -export visual_size -apply() = 0 -export apply - -# define some utility functions -right(t::OpticSim.Geometry.Transform{T}) where {T<:Real} = normalize(Vec3(t[1,1], t[2,1], t[3,1])) -up(t::OpticSim.Geometry.Transform{T}) where {T<:Real} = normalize(Vec3(t[1,2], t[2,2], t[3,2])) -forward(t::OpticSim.Geometry.Transform{T}) where {T<:Real} = normalize(Vec3(t[1,3], t[2,3], t[3,3])) -OpticSim.origin(t::OpticSim.Geometry.Transform{T}) where {T<:Real} = Vec3(t[1,4], t[2,4], t[3,4]) -export right, up, forward - -#------------------------------------------------------------------------------ -# Spectrum -#------------------------------------------------------------------------------ -#region Spectrum -module Spectrum - -using ..OpticSim -using DataFrames -using Distributions - -using ..Emitters - -abstract type AbstractSpectrum{T<:Real} end - -const UNIFORMSHORT = 0.450 #um -const UNIFORMLONG = 0.680 #um - -#--------------------------------------- -# Uniform Spectrum -#--------------------------------------- -""" -Encapsulates a flat spectrum range which is sampled uniformly. Unless stated diferrently, the range used will be 450nm to 680nm. - -```julia -Uniform(low::T, high::T) where {T<:Real} -Uniform(::Type{T} = Float64) where {T<:Real} -``` -""" -struct Uniform{T} <: AbstractSpectrum{T} - low_end::T - high_end::T - - # user defined range of spectrum - function Uniform(low::T, high::T) where {T<:Real} - return new{T}(low, high) - end - - # with no specific range we will use the constants' values - function Uniform(::Type{T} = Float64) where {T<:Real} - return new{T}(UNIFORMSHORT, UNIFORMLONG) - end -end -export Uniform - -Emitters.generate(s::Uniform{T}) where {T} = (T(1), rand(Distributions.Uniform(s.low_end, s.high_end))) - -#--------------------------------------- -# Delta Function Spectrum -#--------------------------------------- -""" -Encapsulates a constant spectrum. - -```julia -DeltaFunction{T<:Real} -``` -""" -struct DeltaFunction{T<:Real} <: AbstractSpectrum{T} - λ::T -end -export DeltaFunction - -Emitters.generate(s::DeltaFunction{T}) where {T} = (T(1), s.λ) - -#--------------------------------------- -# Measured Spectrum -#--------------------------------------- - -""" -Encapsulates a measured spectrum to compute emitter power. Create spectrum by reading CSV files. -Evaluate spectrum at arbitrary wavelength with [`spectrumpower`](@ref) (**more technical details coming soon**) - -```julia -Measured(samples::DataFrame) -``` -""" -struct Measured{T<:Real} <: AbstractSpectrum{T} - low_wave_length::Int - high_wave_length::Int - wave_length_step::Int - power_samples::Vector{T} - - function Measured(samples::DataFrame) - colnames = names(samples) - @assert "Wavelength" in colnames - @assert "Power" in colnames - wavelengths = samples[!, :Wavelength] - @assert eltype(wavelengths) <: Integer - - power::Vector{T} where {T<:Real} = samples[!, :Power] # no missing values allowed and must be real numbers - maxpower = maximum(power) - T = eltype(power) - p = sortperm(wavelengths) - wavelengths = wavelengths[p] - power = power[p] ./ maxpower #make sure power values have the same sorted permutation as wavelengths normalized with maximum power equal to 1 - λmin = wavelengths[p[1]] - λmax = wavelengths[p[end]] - step = wavelengths[2] - wavelengths[1] - - for i in 3:length(wavelengths) - @assert wavelengths[i] - wavelengths[i - 1] == step # no missing values allowed and step size must be consistent - end - - return new{T}(λmin, λmax, step, power) - end -end -export Measured - -"""expects wavelength in nm not um""" -function spectrumpower(spectrum::Measured{T}, λ::T)::Union{Nothing,T} where {T<:Real} - λ = λ #convert from um to nm - if λ < spectrum.low_wave_length || λ > spectrum.high_wave_length - return nothing - end - - lowindex = floor(Int, (λ - spectrum.low_wave_length)) ÷ spectrum.wave_length_step + 1 - - if lowindex == length(spectrum.power_samples) - return T(spectrum.power_samples[end]) - else - highindex = lowindex + 1 - α = mod(λ, spectrum.wave_length_step) / spectrum.wave_length_step - - return T((1 - α) * spectrum.power_samples[lowindex] + α * spectrum.power_samples[highindex]) - end -end - - -function Emitters.generate(s::Measured{T}) where {T} - spectrum = s - λ = rand(Distributions.Uniform{T}(T(spectrum.low_wave_length), T(spectrum.high_wave_length))) - power = spectrumpower(spectrum, λ) - if power === nothing #added this condition because the compiler was allocating if just returned values directly. - return (nothing, nothing) - else - return (power, λ / T(1000)) #convert from nm to um - end -end - - -end # module Spectrum -#endregion Spectrum - -#------------------------------------------------------------------------------ -# Directions Distribution -#------------------------------------------------------------------------------ -#region Direction Distribution -module Directions - -using OpticSim, OpticSim.Geometry -using LinearAlgebra -using ..Emitters - -abstract type AbstractDirectionDistribution{T<:Real} end - -#--------------------------------------- -# Direction Distrubution Common Utilities -#--------------------------------------- - -export generate - -Base.iterate(a::AbstractDirectionDistribution, state = 1) = state > length(a) ? nothing : (generate(a, state - 1), state + 1) -Base.getindex(a::AbstractDirectionDistribution, index) = generate(a, index) -Base.firstindex(a::AbstractDirectionDistribution) = 0 -Base.lastindex(a::AbstractDirectionDistribution) = length(a) - 1 -Base.copy(a::AbstractDirectionDistribution) = a # most don't have any heap allocated stuff so don't really need copying - - -#--------------------------------------- -# Constant Ray Direction -#--------------------------------------- -""" -Encapsulates a single ray direction, where the default direction is unitZ3 [0, 0, 1]. - -```julia -Constant(direction::Vec3{T}) where {T<:Real} -Constant(::Type{T} = Float64) where {T<:Real} -``` -""" -struct Constant{T} <: AbstractDirectionDistribution{T} - direction::Vec3{T} - - function Constant(direction::Vec3{T}) where {T<:Real} - return new{T}(direction) - end - - function Constant(::Type{T} = Float64) where {T<:Real} - direction = unitZ3(T) - return new{T}(direction) - end -end -export Constant - -Base.length(d::Constant{T}) where {T} = 1 -function Emitters.generate(d::Constant{T}, n::Int) where {T<:Real} - return d.direction -end - -#--------------------------------------- -# Grid Ray Distribution -#--------------------------------------- -""" -Encapsulates a single ray direction, where the default direction is unitZ3 [0, 0, 1]. - -```julia -Constant(direction::Vec3{T}) where {T<:Real} -Constant(::Type{T} = Float64) where {T<:Real} -``` -""" -struct RectGrid{T} <: AbstractDirectionDistribution{T} - direction::Vec3{T} - halfangleu::T - halfanglev::T - nraysu::Int - nraysv::Int - uvec::Vec3{T} - vvec::Vec3{T} - - function RectGrid(direction::Vec3{T}, halfangleu::T, halfanglev::T, numraysu::Int, numraysv::Int) where {T<:Real} - (uvec, vvec) = get_uv_vectors(direction) - return new{T}(direction, halfangleu, halfanglev, numraysu, numraysv, uvec, vvec) - end - - function RectGrid(halfangleu::T, halfanglev::T, numraysu::Int, numraysv::Int) where {T<:Real} - direction, uvec, vvec = (unitZ3(T), unitX3(T), unitY3(T)) - return new{T}(direction, halfangleu, halfanglev, numraysu, numraysv, uvec, vvec) - end -end -export RectGrid - -Base.length(d::RectGrid{T}) where {T} = d.nraysu * d.nraysv -function Emitters.generate(d::RectGrid{T}, n::Int) where {T<:Real} - direction = d.direction - uvec = d.uvec - vvec = d.vvec - - # distributing evenly across the area of the rectangle which subtends the given angle (*not* evenly across the angle range) - dindex = mod(n, d.nraysu * d.nraysv) - v = d.nraysv == 1 ? zero(T) : 2 * Int(floor(dindex / d.nraysu)) / (d.nraysv - 1) - 1.0 - u = d.nraysu == 1 ? zero(T) : 2 * mod(dindex, d.nraysu) / (d.nraysu - 1) - 1.0 - θu = atan(u * tan(d.halfangleu) / 2) * d.halfangleu - θv = atan(v * tan(d.halfanglev) / 2) * d.halfanglev - dir = cos(θv) * (cos(θu) * direction + sin(θu) * uvec) + sin(θv) * vvec - return dir -end - - -#--------------------------------------- -# Cone Ray Distribution -#--------------------------------------- -""" -Encapsulates a numsamples rays sampled uniformlly from a cone with max angle θmax. - -```julia -UniformCone(direction::Vec3{T}, θmax::T, numsamples::Int) where {T<:Real} -UniformCone(θmax::T, numsamples::Int) where {T<:Real} -``` -""" -struct UniformCone{T} <: AbstractDirectionDistribution{T} - direction::Vec3{T} - θmax::T - numsamples::Int - uvec::Vec3{T} - vvec::Vec3{T} - - function UniformCone(direction::Vec3{T}, θmax::T, numsamples::Int) where {T<:Real} - (uvec, vvec) = get_uv_vectors(direction) - return new{T}(direction, θmax, numsamples, uvec, vvec) - end - - function UniformCone(θmax::T, numsamples::Int) where {T<:Real} - direction, uvec, vvec = (unitZ3(T), unitX3(T), unitY3(T)) - return new{T}(direction, θmax, numsamples, uvec, vvec) - end -end -export UniformCone - -Base.length(d::UniformCone{T}) where {T} = d.numsamples -function Emitters.generate(d::UniformCone{T}, n::Int) where {T<:Real} - direction = d.direction - θmax = d.θmax - uvec = d.uvec - vvec = d.vvec - - ϕ = rand(T) * 2π - θ = acos(clamp(one(T) + rand(T) * (cos(θmax) - 1), -one(T), one(T))) - return normalize(sin(θ) * (cos(ϕ) * uvec + sin(ϕ) * vvec) + cos(θ) * direction) -end - -#---------------------------------------------------------------------------- -# Cone Ray Hexapolar Distribution -#----------------------------------------------------------------------------- -""" -Rays are generated by sampling a cone with θmax angle in an hexapolar fasion. The number of rays depends on the requested rings and is computed using the following formula: -`1 + round(Int, (nrings * (nrings + 1) / 2) * 6)` - -```julia -HexapolarCone(direction::Vec3{T}, θmax::T, nrings::Int) where {T<:Real} -HexapolarCone(θmax::T, nrings::Int = 3) where {T<:Real} -``` -""" -struct HexapolarCone{T} <: AbstractDirectionDistribution{T} - direction::Vec3{T} - θmax::T - nrings::Int - uvec::Vec3{T} - vvec::Vec3{T} - - function HexapolarCone(direction::Vec3{T}, θmax::T, nrings::Int) where {T<:Real} - (uvec, vvec) = get_orthogonal_vectors(direction) - return new{T}(direction, θmax, nrings, uvec, vvec) - end - - # assume canonical directions - function HexapolarCone(θmax::T, nrings::Int = 3) where {T<:Real} - direction, uvec, vvec = (unitZ3(T), unitX3(T), unitY3(T)) - return new{T}(direction, θmax, nrings, uvec, vvec) - end -end -export HexapolarCone - - -Base.length(d::HexapolarCone{T}) where {T} = 1 + round(Int, (d.nrings * (d.nrings + 1) / 2) * 6) -function Emitters.generate(d::HexapolarCone{T}, n::Int) where {T<:Real} - dir = d.direction - θmax = d.θmax - uvec = d.uvec - vvec = d.vvec - - n = mod(n, length(d)) - if n == 0 - return normalize(dir) - else - t = 1 - ringi = 1 - for i in 1:(d.nrings) - t += 6 * i - if n < t - ringi = i - break - end - end - ρ = ringi / d.nrings - pind = n - (t - 6 * ringi) - - ϕ = (pind / (6 * ringi)) * 2π - # elevation calculated as ring fraction multipled by max angle - θ = acos(clamp(one(T) + (cos(ρ * θmax) - 1), -one(T), one(T))) - return normalize(sin(θ) * (cos(ϕ) * uvec + sin(ϕ) * vvec) + cos(θ) * dir) - end -end - -end # module Directions -#endregion Direction Distribution - -#------------------------------------------------------------------------------ -# Origins Distribution -#------------------------------------------------------------------------------ -#region Origins -module Origins - -using OpticSim, OpticSim.Geometry -using LinearAlgebra -using Distributions -using ..Emitters - - -abstract type AbstractOriginDistribution{T<:Real} end - -#--------------------------------------- -# Origin Distrubution Common Utilities -#--------------------------------------- - -Base.iterate(a::AbstractOriginDistribution, state = 1) = state > length(a) ? nothing : (generate(a, state - 1), state + 1) -Base.getindex(a::AbstractOriginDistribution, index) = generateray(a, index) -Base.firstindex(a::AbstractOriginDistribution) = 0 -Base.lastindex(a::AbstractOriginDistribution) = length(a) - 1 -Base.copy(a::AbstractOriginDistribution) = a # most don't have any heap allocated stuff so don't really need copying - -export generate - -#-------------------------------------- -# Point Origin -#-------------------------------------- -""" -Encapsulates a single point origin. - -```julia -Point(position::Vec3{T}) where {T<:Real} -Point(x::T, y::T, z::T) where {T<:Real} -Point(::Type{T} = Float64) where {T<:Real} -``` -""" -struct Point{T} <: AbstractOriginDistribution{T} - origin::Vec3{T} - - function Point(position::Vec3{T}) where {T<:Real} - return new{T}(position) - end - - function Point(x::T, y::T, z::T) where {T<:Real} - return new{T}(Vec3{T}(x, y, z)) - end - - function Point(::Type{T} = Float64) where {T<:Real} - return new{T}(zero(Vec3)) - end - -end -export Point - -Base.length(o::Point{T}) where {T} = 1 -Emitters.visual_size(o::Point{T}) where {T} = 1 - -function Emitters.generate(o::Point{T}, n::Int) where {T<:Real} - return o.origin -end - -#------------------------------------- -# Random Rectangle Origin -#------------------------------------- -""" -Encapsulates a uniformly sampled rectangle with user defined number of samples. - -```julia -RectUniform(width::T, height::T, count::Int) where {T<:Real} -``` -""" -struct RectUniform{T} <: AbstractOriginDistribution{T} - width::T - height::T - samples_count::Int - - function RectUniform(width::T, height::T, count::Int) where {T<:Real} - return new{T}(width, height, count) - end -end -export RectUniform - - -Base.length(o::RectUniform{T}) where {T} = o.samples_count -Emitters.visual_size(o::RectUniform{T}) where {T} = max(o.width, o.height) - -# generate origin on the agrid -function Emitters.generate(o::RectUniform{T}, n::Int) where {T<:Real} - n = mod(n, length(o)) - u = rand(Distributions.Uniform(T(-1), T(1))) - v = rand(Distributions.Uniform(T(-1), T(1))) - return zero(Vec3{T}) + ((o.width / 2) * u * unitX3(T)) + ((o.height/2) * v * unitY3(T)) -end - - -#------------------------------------- -# Grid Rectangle Origin -#------------------------------------- -""" -Encapsulates a rectangle sampled in a grid fashion. - -```julia -RectGrid(width::T, height::T, usamples::Int, vsamples::Int) where {T<:Real} -``` -""" -struct RectGrid{T} <: AbstractOriginDistribution{T} - width::T - height::T - usamples::Int - vsamples::Int - ustep::T - vstep::T - - function RectGrid(width::T, height::T, usamples::Int, vsamples::Int) where {T<:Real} - return new{T}(width, height, usamples, vsamples, width / (usamples - 1), height / (vsamples - 1)) - end -end -export RectGrid - -Base.length(o::RectGrid{T}) where {T} = o.usamples * o.vsamples -Emitters.visual_size(o::RectGrid{T}) where {T} = max(o.width, o.height) - -# generate origin on the agrid -function Emitters.generate(o::RectGrid{T}, n::Int) where {T<:Real} - n = mod(n, length(o)) - v = o.vsamples == 1 ? zero(T) : 2 * Int(floor(n / o.usamples)) / (o.vsamples - 1) - 1.0 - u = o.usamples == 1 ? zero(T) : 2 * mod(n, o.usamples) / (o.usamples - 1) - 1.0 - return zeros(Vec3{T}) + ((o.width / 2) * u * unitX3(T)) + ((o.height/2) * v * unitY3(T)) -end - - -#------------------------------------- -# Hexapolar Origin -#------------------------------------- -""" -Encapsulates an ellipse (or a circle where halfsizeu=halfsizev) sampled in an hexapolar fasion (rings) - -```julia -Hexapolar(nrings::Int, halfsizeu::T, halfsizev::T) where {T<:Real} -``` -""" -struct Hexapolar{T} <: AbstractOriginDistribution{T} - halfsizeu::T - halfsizev::T - nrings::Int - - function Hexapolar(nrings::Int, halfsizeu::T, halfsizev::T) where {T<:Real} - return new{T}(halfsizeu, halfsizev, nrings) - end -end -export Hexapolar - -Base.length(o::Hexapolar{T}) where {T} = 1 + round(Int, (o.nrings * (o.nrings + 1) / 2) * 6) -Emitters.visual_size(o::Hexapolar{T}) where {T} = max(o.halfsizeu*2, o.halfsizev*2) - -function Emitters.generate(o::Hexapolar{T}, n::Int) where {T<:Real} - n = mod(n, length(o)) - if n == 0 - return zeros(Vec3{T}) - else - t = 1 - ringi = 1 - for i in 1:(o.nrings) - t += 6 * i - if n < t - ringi = i - break - end - end - ρ = ringi / o.nrings - pind = n - (t - 6 * ringi) - ϕ = (pind / (6 * ringi)) * 2π - u = cos(ϕ) * o.halfsizeu - v = sin(ϕ) * o.halfsizev - return zeros(Vec3{T}) + ρ * (u * unitX3(T) + v * unitY3(T)) - end -end - - - -end # module Origins -#endregion Origins - - -#------------------------------------------------------------------------------ -# Angular Power Distribution -#------------------------------------------------------------------------------ -#region Angular Power Distribution -module AngularPower - -using OpticSim, OpticSim.Geometry -using LinearAlgebra -using ..Emitters - -abstract type AbstractAngularPowerDistribution{T<:Real} end - -#--------------------------------------- -# Lambertian -#--------------------------------------- -""" -Ray power is unaffected by angle. -""" -struct Lambertian{T} <: AbstractAngularPowerDistribution{T} - Lambertian(::Type{T} = Float64) where {T<:Real} = new{T}() -end - -function Emitters.apply(d::Lambertian{T}, Tr::Transform{T}, power::T, ray::OpticSim.Ray{T,3}) where {T} - return power -end -export Lambertian - -#--------------------------------------- -# Cosine Power Distribution -#--------------------------------------- -""" - Cosine(cosine_exp::T = one(T)) where {T<:Real} - -Cosine power distribution. Ray power is calculated by: - -`power = power * (cosine_angle ^ cosine_exp)` -""" -struct Cosine{T} <: AbstractAngularPowerDistribution{T} - cosine_exp::T - - function Cosine(cosine_exp::T = one(T)) where {T<:Real} - new{T}(cosine_exp) - end -end -export Cosine - -# returning ray power -function Emitters.apply(d::Cosine{T}, Tr::Transform{T}, power::T, ray::OpticSim.Ray{T,3}) where {T} - cosanglebetween = dot(OpticSim.direction(ray), forward(Tr)) - power = power * cosanglebetween ^ d.cosine_exp - return power -end - -#--------------------------------------- -# Gaussian Power Distribution -#--------------------------------------- -""" - Gaussian(gaussianu::T, gaussianv::T) where {T<:Real} - -GGaussian power distribution. Ray power is calculated by: - -`power = power * exp(-(gaussianu * l^2 + gaussianv * m^2))` -where l and m are the cos_angles between the two axes respectivly. -""" -struct Gaussian{T} <: AbstractAngularPowerDistribution{T} - gaussianu::T - gaussianv::T - - function Gaussian(gaussianu::T, gaussianv::T) where {T<:Real} - new{T}(gaussianu, gaussianv) - end -end -export Gaussian - -# returning ray power -function Emitters.apply(d::Gaussian{T}, Tr::Transform{T}, power::T, ray::OpticSim.Ray{T,3}) where {T} - - l = dot(OpticSim.direction(ray), right(Tr)) - m = dot(OpticSim.direction(ray), up(Tr)) - power = power * exp(-(d.gaussianu * l^2 + d.gaussianv * m^2)) - - return power -end - - -end # module Angular Power -#endregion Angular Power Distribution - -#------------------------------------------------------------------------------ -# Sources -#------------------------------------------------------------------------------ -#region Sources -module Sources - -using OpticSim, OpticSim.Geometry -using LinearAlgebra -using Distributions -using ..Emitters -using ..Spectrum -using ..Directions -using ..Origins -using ..AngularPower - -abstract type AbstractSource{T<:Real} <: OpticSim.OpticalRayGenerator{T} end - -Base.getindex(s::AbstractSource, index) = generate(s, index) -Base.firstindex(s::AbstractSource) = 0 -Base.lastindex(s::AbstractSource) = length(s) - 1 -Base.copy(a::AbstractSource) = a # most don't have any heap allocated stuff so don't really need copying - - -#--------------------------------------- -# Source -#--------------------------------------- -""" -This data-type represents the basic emitter (Source), which is a combination of a Spectrum, Angular Power Distribution, Origins and Directions distibution and a 3D Transform. - -```julia -Source(::Type{T} = Float64; - transform::Tr = Transform(), - spectrum::S = Spectrum.Uniform(), - origins::O = Origins.Point(), - directions::D = Directions.Constant(), - power::P = AngularPower.Lambertian(), - sourcenum::Int = 0) where { - Tr<:Transform, - S<:Spectrum.AbstractSpectrum, - O<:Origins.AbstractOriginDistribution, - D<:Directions.AbstractDirectionDistribution, - P<:AngularPower.AbstractAngularPowerDistribution, - T<:Real} -Source(transform::Tr, spectrum::S, origins::O, directions::D, power::P, ::Type{T} = Float64; sourcenum::Int = 0) where { - Tr<:Transform, - S<:Spectrum.AbstractSpectrum, - O<:Origins.AbstractOriginDistribution, - D<:Directions.AbstractDirectionDistribution, - P<:AngularPower.AbstractAngularPowerDistribution, - T<:Real} - -``` -""" -struct Source{ T<:Real, - Tr<:Transform{T}, - S<:Spectrum.AbstractSpectrum{T}, - O<:Origins.AbstractOriginDistribution{T}, - D<:Directions.AbstractDirectionDistribution{T}, - P<:AngularPower.AbstractAngularPowerDistribution{T}} <: AbstractSource{T} - transform::Tr - spectrum::S - origins::O - directions::D - power_distribution::P - sourcenum::Int - - function Source(::Type{T} = Float64; - transform::Tr = Transform(), - spectrum::S = Spectrum.Uniform(), - origins::O = Origins.Point(), - directions::D = Directions.Constant(), - power::P = AngularPower.Lambertian(), - sourcenum::Int = 0) where { - Tr<:Transform, - S<:Spectrum.AbstractSpectrum, - O<:Origins.AbstractOriginDistribution, - D<:Directions.AbstractDirectionDistribution, - P<:AngularPower.AbstractAngularPowerDistribution, - T<:Real} - new{T, Tr, S, O, D, P}(transform, spectrum, origins, directions, power, sourcenum) - end - - function Source(transform::Tr, spectrum::S, origins::O, directions::D, power::P, ::Type{T} = Float64; sourcenum::Int = 0) where { - Tr<:Transform, - S<:Spectrum.AbstractSpectrum, - O<:Origins.AbstractOriginDistribution, - D<:Directions.AbstractDirectionDistribution, - P<:AngularPower.AbstractAngularPowerDistribution, - T<:Real} - new{T, Tr, S, O, D, P}(transform, spectrum, origins, directions, power, sourcenum) - end -end -export Source - -# used to not generate new origin points if we can use last points - mainly to keep consistency of origin generation when randomness is involved -struct SourceGenerationState{T<:Real} - n::Int - last_index_used::Int - last_point_generated::Vec3{T} -end - -Base.length(s::Source{T}) where {T} = length(s.origins) * length(s.directions) - -function Base.iterate(s::Source{T}) where {T<:Real} - state = SourceGenerationState(1, -1, zeros(Vec3{T})) - return iterate(s, state) -end -function Base.iterate(s::Source{T}, state::SourceGenerationState{T}) where {T<:Real} - if (state.n > length(s)) - return nothing - end - - return generate(s, state) -end - - -function Emitters.generate(s::Source{T}, n::Int) where {T<:Real} - return generate(s, SourceGenerationState(n+1, -1, zeros(Vec3{T})))[1] -end - -function Emitters.generate(s::Source{T}, state::SourceGenerationState{T} = SourceGenerationState(-1, -1, zeros(Vec3{T}))) where {T<:Real} - - # @info "New State Generation" - n::Int = state.n - 1 - origin_index = floor(Int, (n / length(s.directions))) - direction_index = mod(n, length(s.directions)) - - if (origin_index == state.last_index_used) - o = state.last_point_generated - else - o = generate(s.origins, origin_index) - o = s.transform * o - end - - # o = generate(s.origins, origin_index) - d = generate(s.directions, direction_index) - - # rotate the direction according to the orientation of the ray-originator - d = rotation(s.transform) * d - - ray = OpticSim.Ray(o, d) - power, wavelength = generate(s.spectrum) - - power = apply(s.power_distribution, s.transform, power, ray) - - # return (OpticSim.Ray(o, d), SourceGenerationState(state.n+1, origin_index, o)) - return (OpticSim.OpticalRay(ray, power, wavelength; sourcenum=s.sourcenum), SourceGenerationState(state.n+1, origin_index, o)) -end - -#--------------------------------------- -# Composite Source -#--------------------------------------- -""" -This data-type represents the composit emitter (Source) which is constructed with a list of basic or composite emitters and a 3D Transform. - -```julia -CompositeSource(transform::Transform{T}, sources::Array, ::Type{T} = Float64) where {T<:Real} -``` -""" -struct CompositeSource{T} <: AbstractSource{T} - transform::Transform{T} - sources::Array - - uniform_length::Int - total_length::Int - start_indexes::Vector{Int} - - function CompositeSource(transform::Transform{T}, sources::Array, ::Type{T} = Float64) where {T<:Real} - - lens = [length(src) for src in sources] - size1 = findmin(lens)[1] - size2 = findmax(lens)[1] - if (size1 == size2) - uniform_length = size1 - start_indexes = Vector{Int}() - total_length = length(sources) * uniform_length - # @info "Uniform Length $size1 total=$total_length" - else - uniform_length = -1 - - start_indexes = Vector{Int}() - start = 0 - for s in sources - push!(start_indexes, start) - start = start + length(s) - end - push!(start_indexes, start) - # @info start_indexes - total_length = start - end - - new{T}(transform, sources, uniform_length, total_length, start_indexes) - end -end -export CompositeSource - -Base.length(s::CompositeSource{T}) where {T} = s.total_length - -function Base.iterate(s::CompositeSource{T}) where {T<:Real} - state = SourceGenerationState(1, -1, zeros(Vec3{T})) - return iterate(s, state) -end -function Base.iterate(s::CompositeSource{T}, state::SourceGenerationState{T}) where {T<:Real} - if (state.n > length(s)) - return nothing - end - - return generate(s, state) -end - -function Emitters.generate(s::CompositeSource{T}, n::Int) where {T<:Real} - return generate(s, SourceGenerationState(n+1, -1, zeros(Vec3{T})))[1] -end - -function Emitters.generate(s::CompositeSource{T}, state::SourceGenerationState{T} = SourceGenerationState(-1, -1, zeros(Vec3{T}))) where {T<:Real} - - # @info "New Composite State Generation" - n::Int = state.n - 1 - - if (s.uniform_length != -1) - source_index = floor(Int, (n / s.uniform_length)) - ray_index = mod(n, s.uniform_length) - # @info "New Composite State Generation: source=$source_index ray=$ray_index ($(length(s.sources)))" - else - n = mod(n, s.total_length) - - ## perform a binary search to find out which source can generate the spesific ray - low = 1 - high = length(s.start_indexes) - 1 # last cell in this vector is a dummy cell containing the index after the last ray - source_index = -1 - while (low <= high) - mid = floor(Int, (low + high) / 2) - - if (n < s.start_indexes[mid]) - high = mid - 1 - elseif (n >= s.start_indexes[mid+1]) - low = mid + 1 - else - source_index = mid - break; - end - end - @assert source_index != -1 - # @info "Bin Search n=$n source_index=$source_index" - source_index = source_index - 1 - ray_index = n - s.start_indexes[source_index+1] - end - - ray, new_state = generate(s.sources[source_index+1], SourceGenerationState(ray_index+1, state.last_index_used, state.last_point_generated)) - ray = s.transform * ray - return (ray, SourceGenerationState(state.n+1, new_state.last_index_used, new_state.last_point_generated)) -end - - -end # module Sources -#endregion Sources - -#------------------------------------------------------------------------------ -# Visualization -#------------------------------------------------------------------------------ -#region Visualization -module Visualization - -using OpticSim, OpticSim.Geometry -using LinearAlgebra -using Distributions -using StaticArrays - -import Makie -import Makie.AbstractPlotting -import Makie.AbstractPlotting.MakieLayout - -using ..Emitters -using ..Spectrum -using ..Directions -using ..Origins -using ..AngularPower -using ..Sources - -const ARRROW_LENGTH = 0.5 -const ARRROW_SIZE = 0.01 -const MARKER_SIZE = 15 - - -#------------------------------------- -# draw debug information - local axes and positions -#------------------------------------- -function maybe_draw_debug_info(scene::MakieLayout.LScene, o::Origins.AbstractOriginDistribution; transform::Geometry.Transform = Transform(), debug::Bool=false, kwargs...) where {T<:Real} - - dir = forward(transform) - uv = SVector{3}(right(transform)) - vv = SVector{3}(up(transform)) - pos = origin(transform) - - if (debug) - # this is a stupid hack to force makie to render in 3d - for some scenes, makie decide with no apperent reason to show in 2d instead of 3d - AbstractPlotting.scatter!(scene, [pos[1], pos[1]+0.1], [pos[2], pos[2]+0.1], [pos[3], pos[3]+0.1], color=:red, markersize=0) - - # draw the origin and normal of the surface - Makie.scatter!(scene, pos, color=:blue, markersize = MARKER_SIZE * visual_size(o)) - - # normal - arrow_start = pos - arrow_end = dir * ARRROW_LENGTH * visual_size(o) - Makie.arrows!(scene.scene, [AbstractPlotting.Point3f0(arrow_start)], [AbstractPlotting.Point3f0(arrow_end)], arrowsize=ARRROW_SIZE * visual_size(o), arrowcolor=:blue) - arrow_end = uv * 0.5 * ARRROW_LENGTH * visual_size(o) - Makie.arrows!(scene.scene, [AbstractPlotting.Point3f0(arrow_start)], [AbstractPlotting.Point3f0(arrow_end)], arrowsize= 0.5 * ARRROW_SIZE * visual_size(o), arrowcolor=:red) - arrow_end = vv * 0.5 * ARRROW_LENGTH * visual_size(o) - Makie.arrows!(scene.scene, [AbstractPlotting.Point3f0(arrow_start)], [AbstractPlotting.Point3f0(arrow_end)], arrowsize= 0.5 * ARRROW_SIZE * visual_size(o), arrowcolor=:green) - - # draw all the samples origins - positions = map(x -> transform*x, collect(o)) - positions = collect(AbstractPlotting.Point3f0, positions) - Makie.scatter!(scene, positions, color=:green, markersize = MARKER_SIZE * visual_size(o)) - - # positions = collect(AbstractPlotting.Point3f0, o) - # Makie.scatter!(scene, positions, color=:green, markersize = MARKER_SIZE * visual_size(o)) - end - -end - - -#------------------------------------- -# draw point origin -#------------------------------------- -# function OpticSim.Vis.draw!(scene::MakieLayout.LScene, o::Origins.Point; transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real} -function OpticSim.Vis.draw!(scene::MakieLayout.LScene, o::Origins.Point; transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real} - maybe_draw_debug_info(scene, o; transform=transform, kwargs...) -end - -#------------------------------------- -# draw RectGrid and RectUniform origins -#------------------------------------- -function OpticSim.Vis.draw!(scene::MakieLayout.LScene, o::Union{Origins.RectGrid, Origins.RectUniform}; transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real} - dir = forward(transform) - uv = SVector{3}(right(transform)) - vv = SVector{3}(up(transform)) - pos = origin(transform) - - # @info "RECT: transform $(pos)" - - plane = OpticSim.Plane(dir, pos) - rect = OpticSim.Rectangle(plane, o.width / 2, o.height / 2, uv, vv) - - OpticSim.Vis.draw!(scene, rect; kwargs...) - - maybe_draw_debug_info(scene, o; transform=transform, kwargs...) -end - - -#------------------------------------- -# draw hexapolar origin -#------------------------------------- -function OpticSim.Vis.draw!(scene::MakieLayout.LScene, o::Origins.Hexapolar; transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real} - dir = forward(transform) - uv = SVector{3}(right(transform)) - vv = SVector{3}(up(transform)) - pos = origin(transform) - - plane = OpticSim.Plane(dir, pos) - ellipse = OpticSim.Ellipse(plane, o.halfsizeu, o.halfsizev, uv, vv) - - OpticSim.Vis.draw!(scene, ellipse; kwargs...) - - maybe_draw_debug_info(scene, o; transform=transform, kwargs...) -end - -#------------------------------------- -# draw source -#------------------------------------- -function OpticSim.Vis.draw!(scene::MakieLayout.LScene, s::Sources.Source{T}; parent_transform::Geometry.Transform = Transform(), debug::Bool=false, kwargs...) where {T<:Real} - - OpticSim.Vis.draw!(scene, s.origins; transform=parent_transform * s.transform, debug=debug, kwargs...) - - if (debug) - m = zeros(T, length(s), 7) - for (index, optical_ray) in enumerate(s) - ray = OpticSim.ray(optical_ray) - ray = parent_transform * ray - m[index, 1:7] = [ray.origin... ray.direction... OpticSim.power(optical_ray)] - end - - m[:, 4:6] .*= m[:, 7] * ARRROW_LENGTH * visual_size(s.origins) - - # Makie.arrows!(scene, [Makie.Point3f0(origin(ray))], [Makie.Point3f0(rayscale * direction(ray))]; kwargs..., arrowsize = min(0.05, rayscale * 0.05), arrowcolor = color, linecolor = color, linewidth = 2) - color = :yellow - Makie.arrows!(scene, m[:,1], m[:,2], m[:,3], m[:,4], m[:,5], m[:,6]; kwargs..., arrowcolor = color, linecolor = color, linewidth = 2, arrowsize=ARRROW_SIZE * visual_size(s.origins)) - end - - # for ray in o - # OpticSim.Vis.draw!(scene, ray) - # end -end - -#------------------------------------- -# draw optical rays -#------------------------------------- -function OpticSim.Vis.draw!(scene::MakieLayout.LScene, rays::AbstractVector{OpticSim.OpticalRay{T, 3}}; kwargs...) where {T<:Real} - m = zeros(T, length(rays)*2, 3) - for (index, optical_ray) in enumerate(rays) - ray = OpticSim.ray(optical_ray) - m[(index-1)*2+1, 1:3] = [origin(optical_ray)...] - m[(index-1)*2+2, 1:3] = [(OpticSim.origin(optical_ray) + OpticSim.direction(optical_ray) * 1 * OpticSim.power(optical_ray))... ] - end - - color = :green - Makie.linesegments!(scene, m[:,1], m[:,2], m[:,3]; kwargs..., color = color, linewidth = 2, ) -end - -#------------------------------------- -# draw composite source -#------------------------------------- -function OpticSim.Vis.draw!(scene::MakieLayout.LScene, s::Sources.CompositeSource{T}; parent_transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real} - for source in s.sources - OpticSim.Vis.draw!(scene, source; parent_transform=parent_transform*s.transform, kwargs...) - end -end - - -end # module Visualization -#endregion Visualization - - -# exporting stuff from the Emitters module -export Geometry, Spectrum, Directions, Origins, AngularPower, Sources - -end # module Emitters - -export Emitters diff --git a/src/Optical/Emitters/AngularPower.jl b/src/Optical/Emitters/AngularPower.jl new file mode 100644 index 000000000..22fcf2fe3 --- /dev/null +++ b/src/Optical/Emitters/AngularPower.jl @@ -0,0 +1,65 @@ +module AngularPower +export Lambertian, Cosine, Gaussian + +using ....OpticSim +using ...Emitters +using ...Geometry +using LinearAlgebra + +abstract type AbstractAngularPowerDistribution{T<:Real} end + +""" + Lambertian{T} <: AbstractAngularPowerDistribution{T} + +Ray power is unaffected by angle. +""" +struct Lambertian{T} <: AbstractAngularPowerDistribution{T} + Lambertian(::Type{T} = Float64) where {T<:Real} = new{T}() +end + +Emitters.apply(::Lambertian, ::Transform, power::Real, ::OpticSim.Ray{<:Real,3}) = power + +""" + Cosine{T} <: AbstractAngularPowerDistribution{T} + +Cosine power distribution. Ray power is calculated by: + +`power = power * (cosine_angle ^ cosine_exp)` +""" +struct Cosine{T} <: AbstractAngularPowerDistribution{T} + cosine_exp::T + + function Cosine(cosine_exp::T = one(T)) where {T<:Real} + new{T}(cosine_exp) + end +end + +function Emitters.apply(d::Cosine, tr::Transform, power::Real, ray::OpticSim.Ray{<:Real,3}) + cosanglebetween = dot(OpticSim.direction(ray), forward(tr)) + return power * cosanglebetween ^ d.cosine_exp +end + +""" + Gaussian{T} <: AbstractAngularPowerDistribution{T} + +GGaussian power distribution. Ray power is calculated by: + +`power = power * exp(-(gaussianu * l^2 + gaussianv * m^2))` +where l and m are the cos_angles between the two axes respectivly. +""" +struct Gaussian{T} <: AbstractAngularPowerDistribution{T} + gaussianu::T + gaussianv::T + + function Gaussian(gaussianu::T, gaussianv::T) where {T<:Real} + new{T}(gaussianu, gaussianv) + end +end + +function Emitters.apply(d::Gaussian, tr::Transform, power::Real, ray::OpticSim.Ray{<:Real,3}) + l = dot(OpticSim.direction(ray), right(tr)) + m = dot(OpticSim.direction(ray), up(tr)) + return power * exp(-(d.gaussianu * l^2 + d.gaussianv * m^2)) +end + +end # module Angular Power diff --git a/src/Optical/Emitters/Directions.jl b/src/Optical/Emitters/Directions.jl new file mode 100644 index 000000000..ac7874440 --- /dev/null +++ b/src/Optical/Emitters/Directions.jl @@ -0,0 +1,189 @@ +module Directions +export Constant, RectGrid, UniformCone, HexapolarCone + +using ....OpticSim +using ...Emitters +using ...Geometry +using LinearAlgebra + +abstract type AbstractDirectionDistribution{T<:Real} end + +Base.iterate(a::AbstractDirectionDistribution, state = 1) = state > length(a) ? nothing : (generate(a, state - 1), state + 1) +Base.getindex(a::AbstractDirectionDistribution, index) = generate(a, index) +Base.firstindex(a::AbstractDirectionDistribution) = 0 +Base.lastindex(a::AbstractDirectionDistribution) = length(a) - 1 +Base.copy(a::AbstractDirectionDistribution) = a # most don't have any heap allocated stuff so don't really need copying + +""" + Constant{T} <: AbstractDirectionDistribution{T} + +Encapsulates a single ray direction, where the default direction is unitZ3 [0, 0, 1]. + +```julia +Constant(direction::Vec3{T}) where {T<:Real} +Constant(::Type{T} = Float64) where {T<:Real} +``` +""" +struct Constant{T} <: AbstractDirectionDistribution{T} + direction::Vec3{T} + + function Constant(direction::Vec3{T}) where {T<:Real} + return new{T}(direction) + end + + function Constant(::Type{T} = Float64) where {T<:Real} + direction = unitZ3(T) + return new{T}(direction) + end +end + +Base.length(d::Constant) = 1 +Emitters.generate(d::Constant, ::Integer) = d.direction + +""" + RectGrid{T} <: AbstractDirectionDistribution{T} + +Encapsulates a single ray direction, where the default direction is unitZ3 [0, 0, 1]. + +```julia +Constant(direction::Vec3{T}) where {T<:Real} +Constant(::Type{T} = Float64) where {T<:Real} +``` +""" +struct RectGrid{T} <: AbstractDirectionDistribution{T} + direction::Vec3{T} + halfangleu::T + halfanglev::T + nraysu::Integer + nraysv::Integer + uvec::Vec3{T} + vvec::Vec3{T} + + function RectGrid(direction::Vec3{T}, halfangleu::T, halfanglev::T, numraysu::Integer, numraysv::Integer) where {T<:Real} + (uvec, vvec) = get_uv_vectors(direction) + return new{T}(direction, halfangleu, halfanglev, numraysu, numraysv, uvec, vvec) + end + + function RectGrid(halfangleu::T, halfanglev::T, numraysu::Integer, numraysv::Integer) where {T<:Real} + direction, uvec, vvec = (unitZ3(T), unitX3(T), unitY3(T)) + return new{T}(direction, halfangleu, halfanglev, numraysu, numraysv, uvec, vvec) + end +end + +Base.length(d::RectGrid) = d.nraysu * d.nraysv +function Emitters.generate(d::RectGrid{T}, n::Integer) where {T<:Real} + direction = d.direction + uvec = d.uvec + vvec = d.vvec + + # distributing evenly across the area of the rectangle which subtends the given angle (*not* evenly across the angle range) + dindex = mod(n, d.nraysu * d.nraysv) + v = d.nraysv == 1 ? zero(T) : 2 * Integer(floor(dindex / d.nraysu)) / (d.nraysv - 1) - 1.0 + u = d.nraysu == 1 ? zero(T) : 2 * mod(dindex, d.nraysu) / (d.nraysu - 1) - 1.0 + θu = atan(u * tan(d.halfangleu) / 2) * d.halfangleu + θv = atan(v * tan(d.halfanglev) / 2) * d.halfanglev + dir = cos(θv) * (cos(θu) * direction + sin(θu) * uvec) + sin(θv) * vvec + return dir +end + +""" + UniformCone{T} <: AbstractDirectionDistribution{T} + +Encapsulates `numsamples` rays sampled uniformly from a cone with max angle θmax. + +```julia +UniformCone(direction::Vec3{T}, θmax::T, numsamples::Integer) where {T<:Real} +UniformCone(θmax::T, numsamples::Integer) where {T<:Real} +``` +""" +struct UniformCone{T} <: AbstractDirectionDistribution{T} + direction::Vec3{T} + θmax::T + numsamples::Integer + uvec::Vec3{T} + vvec::Vec3{T} + + function UniformCone(direction::Vec3{T}, θmax::T, numsamples::Integer) where {T<:Real} + (uvec, vvec) = get_uv_vectors(direction) + return new{T}(direction, θmax, numsamples, uvec, vvec) + end + + function UniformCone(θmax::T, numsamples::Integer) where {T<:Real} + direction, uvec, vvec = (unitZ3(T), unitX3(T), unitY3(T)) + return new{T}(direction, θmax, numsamples, uvec, vvec) + end +end + +Base.length(d::UniformCone) = d.numsamples +function Emitters.generate(d::UniformCone{T}, ::Integer) where {T<:Real} + direction = d.direction + θmax = d.θmax + uvec = d.uvec + vvec = d.vvec + + ϕ = rand(T) * 2π + θ = acos(clamp(one(T) + rand(T) * (cos(θmax) - 1), -one(T), one(T))) + return normalize(sin(θ) * (cos(ϕ) * uvec + sin(ϕ) * vvec) + cos(θ) * direction) +end + +""" + HexapolarCone{T} <: AbstractDirectionDistribution{T} + +Rays are generated by sampling a cone with θmax angle in an hexapolar fashion. The number of rays depends on the requested rings and is computed using the following formula: +`1 + round(Integer, (nrings * (nrings + 1) / 2) * 6)` + +```julia +HexapolarCone(direction::Vec3{T}, θmax::T, nrings::Integer) where {T<:Real} +HexapolarCone(θmax::T, nrings::Integer = 3) where {T<:Real} +``` +""" +struct HexapolarCone{T} <: AbstractDirectionDistribution{T} + direction::Vec3{T} + θmax::T + nrings::Integer + uvec::Vec3{T} + vvec::Vec3{T} + + function HexapolarCone(direction::Vec3{T}, θmax::T, nrings::Integer) where {T<:Real} + (uvec, vvec) = get_orthogonal_vectors(direction) + return new{T}(direction, θmax, nrings, uvec, vvec) + end + + # assume canonical directions + function HexapolarCone(θmax::T, nrings::Integer = 3) where {T<:Real} + direction, uvec, vvec = (unitZ3(T), unitX3(T), unitY3(T)) + return new{T}(direction, θmax, nrings, uvec, vvec) + end +end + +Base.length(d::HexapolarCone) = 1 + round(Integer, (d.nrings * (d.nrings + 1) / 2) * 6) +function Emitters.generate(d::HexapolarCone{T}, n::Integer) where {T<:Real} + dir = d.direction + θmax = d.θmax + uvec = d.uvec + vvec = d.vvec + + n = mod(n, length(d)) + if n == 0 + return normalize(dir) + else + t = 1 + ringi = 1 + for i in 1:(d.nrings) + t += 6 * i + if n < t + ringi = i + break + end + end + ρ = ringi / d.nrings + pind = n - (t - 6 * ringi) + + ϕ = (pind / (6 * ringi)) * 2π + # elevation calculated as ring fraction multipled by max angle + θ = acos(clamp(one(T) + (cos(ρ * θmax) - 1), -one(T), one(T))) + return normalize(sin(θ) * (cos(ϕ) * uvec + sin(ϕ) * vvec) + cos(θ) * dir) + end +end + +end # module Directions diff --git a/src/Optical/Emitters/Emitters.jl b/src/Optical/Emitters/Emitters.jl new file mode 100644 index 000000000..0302326ee --- /dev/null +++ b/src/Optical/Emitters/Emitters.jl @@ -0,0 +1,51 @@ +#= +Emitters are defined by Pixels and SpatialLayouts. An emitter has a spectrum, and an optical power distribution over the hemisphere. +These are intrinsic physical properties of the emitter. +=# + +module Emitters +export generate, visual_size, apply +export right, up, forward +export Spectrum, Directions, Origins, AngularPower, Sources + +using ...OpticSim, ..Geometry +using LinearAlgebra + +# define some utility functions +right(t::OpticSim.Geometry.Transform{<:Real}) = normalize(Vec3(t[1,1], t[2,1], t[3,1])) +up(t::OpticSim.Geometry.Transform{<:Real}) = normalize(Vec3(t[1,2], t[2,2], t[3,2])) +forward(t::OpticSim.Geometry.Transform{<:Real}) = normalize(Vec3(t[1,3], t[2,3], t[3,3])) +OpticSim.origin(t::OpticSim.Geometry.Transform{<:Real}) = Vec3(t[1,4], t[2,4], t[3,4]) + +# defining name placeholders to override in nested modules + +""" + generate(???) + +[TODO] +""" +generate() = 0 + +""" + visual_size(???) + +[TODO] +""" +visual_size() = 0 + +""" + apply(???) + +[TODO] Returns ray power +""" +apply() = 0 + +include("Spectrum.jl") +include("Directions.jl") +include("Origins.jl") +include("AngularPower.jl") +include("Sources.jl") + +end # module Emitters + +export Emitters diff --git a/src/Optical/Emitters/Origins.jl b/src/Optical/Emitters/Origins.jl new file mode 100644 index 000000000..43f9bd459 --- /dev/null +++ b/src/Optical/Emitters/Origins.jl @@ -0,0 +1,157 @@ +module Origins +export Point, RectUniform, RectGrid, Hexapolar + +using ....OpticSim +using ...Emitters +using ...Geometry +using LinearAlgebra +using Distributions + +abstract type AbstractOriginDistribution{T<:Real} end + +Base.iterate(a::AbstractOriginDistribution, state = 1) = state > length(a) ? nothing : (generate(a, state - 1), state + 1) +Base.getindex(a::AbstractOriginDistribution, index) = generateray(a, index) +Base.firstindex(a::AbstractOriginDistribution) = 0 +Base.lastindex(a::AbstractOriginDistribution) = length(a) - 1 +Base.copy(a::AbstractOriginDistribution) = a # most don't have any heap allocated stuff so don't really need copying + +""" + Point{T} <: AbstractOriginDistribution{T} + +Encapsulates a single point origin. + +```julia +Point(position::Vec3{T}) where {T<:Real} +Point(x::T, y::T, z::T) where {T<:Real} +Point(::Type{T} = Float64) where {T<:Real} +``` +""" +struct Point{T} <: AbstractOriginDistribution{T} + origin::Vec3{T} + + function Point(position::Vec3{T}) where {T<:Real} + return new{T}(position) + end + + function Point(x::T, y::T, z::T) where {T<:Real} + return new{T}(Vec3{T}(x, y, z)) + end + + function Point(::Type{T} = Float64) where {T<:Real} + return new{T}(zero(Vec3)) + end +end + +Base.length(o::Point) = 1 +Emitters.visual_size(o::Point) = 1 +Emitters.generate(o::Point, ::Integer) = o.origin + +""" + RectUniform{T} <: AbstractOriginDistribution{T} + +Encapsulates a uniformly sampled rectangle with user defined number of samples. + +```julia +RectUniform(width::T, height::T, count::Integer) where {T<:Real} +``` +""" +struct RectUniform{T} <: AbstractOriginDistribution{T} + width::T + height::T + samples_count::Integer + + function RectUniform(width::T, height::T, count::Integer) where {T<:Real} + return new{T}(width, height, count) + end +end + +Base.length(o::RectUniform) = o.samples_count +Emitters.visual_size(o::RectUniform) = max(o.width, o.height) + +# generate origin on the agrid +function Emitters.generate(o::RectUniform{T}, n::Integer) where {T<:Real} + n = mod(n, length(o)) + u = rand(Distributions.Uniform(-one(T), one(T))) + v = rand(Distributions.Uniform(-one(T), one(T))) + return zero(Vec3{T}) + ((o.width / 2) * u * unitX3(T)) + ((o.height/2) * v * unitY3(T)) +end + +""" + RectGrid{T} <: AbstractOriginDistribution{T} + +Encapsulates a rectangle sampled in a grid fashion. + +```julia +RectGrid(width::T, height::T, usamples::Integer, vsamples::Integer) where {T<:Real} +``` +""" +struct RectGrid{T} <: AbstractOriginDistribution{T} + width::T + height::T + usamples::Integer + vsamples::Integer + ustep::T + vstep::T + + function RectGrid(width::T, height::T, usamples::Integer, vsamples::Integer) where {T<:Real} + return new{T}(width, height, usamples, vsamples, width / (usamples - 1), height / (vsamples - 1)) + end +end + +Base.length(o::RectGrid) = o.usamples * o.vsamples +Emitters.visual_size(o::RectGrid) = max(o.width, o.height) + +# generate origin on the agrid +function Emitters.generate(o::RectGrid{T}, n::Integer) where {T<:Real} + n = mod(n, length(o)) + v = o.vsamples == 1 ? zero(T) : 2 * Integer(floor(n / o.usamples)) / (o.vsamples - 1) - 1.0 + u = o.usamples == 1 ? zero(T) : 2 * mod(n, o.usamples) / (o.usamples - 1) - 1.0 + return zeros(Vec3{T}) + ((o.width / 2) * u * unitX3(T)) + ((o.height/2) * v * unitY3(T)) +end + +""" + Hexapolar{T} <: AbstractOriginDistribution{T} + +Encapsulates an ellipse (or a circle where halfsizeu=halfsizev) sampled in an hexapolar fashion (rings). + +```julia +Hexapolar(nrings::Integer, halfsizeu::T, halfsizev::T) where {T<:Real} +``` +""" +struct Hexapolar{T} <: AbstractOriginDistribution{T} + halfsizeu::T + halfsizev::T + nrings::Integer + + function Hexapolar(nrings::Integer, halfsizeu::T, halfsizev::T) where {T<:Real} + return new{T}(halfsizeu, halfsizev, nrings) + end +end + +Base.length(o::Hexapolar) = 1 + round(Integer, (o.nrings * (o.nrings + 1) / 2) * 6) +Emitters.visual_size(o::Hexapolar) = max(o.halfsizeu*2, o.halfsizev*2) + +function Emitters.generate(o::Hexapolar{T}, n::Integer) where {T<:Real} + n = mod(n, length(o)) + if n == 0 + return zeros(Vec3{T}) + else + t = 1 + ringi = 1 + for i in 1:(o.nrings) + t += 6 * i + if n < t + ringi = i + break + end + end + ρ = ringi / o.nrings + pind = n - (t - 6 * ringi) + ϕ = (pind / (6 * ringi)) * 2π + u = cos(ϕ) * o.halfsizeu + v = sin(ϕ) * o.halfsizev + return zeros(Vec3{T}) + ρ * (u * unitX3(T) + v * unitY3(T)) + end +end + +end # module Origins diff --git a/src/Optical/Emitters/Sources.jl b/src/Optical/Emitters/Sources.jl new file mode 100644 index 000000000..e27749eda --- /dev/null +++ b/src/Optical/Emitters/Sources.jl @@ -0,0 +1,237 @@ +module Sources +export Source, CompositeSource + +using ....OpticSim +using ...Emitters +using ...Geometry +using ..Spectrum +using ..Directions +using ..Origins +using ..AngularPower +using LinearAlgebra +using Distributions + +abstract type AbstractSource{T} <: OpticSim.OpticalRayGenerator{T} end + +Base.getindex(s::AbstractSource, index) = generate(s, index) +Base.firstindex(s::AbstractSource) = 0 +Base.lastindex(s::AbstractSource) = length(s) - 1 +Base.copy(a::AbstractSource) = a # most don't have any heap allocated stuff so don't really need copying + +""" + Source{T<:Real, Tr<:Transform{T}, S<:Spectrum.AbstractSpectrum{T}, O<:Origins.AbstractOriginDistribution{T}, D<:Directions.AbstractDirectionDistribution{T}, P<:AngularPower.AbstractAngularPowerDistribution{T}} <: AbstractSource{T} + +This data-type represents the basic emitter (Source), which is a combination of a Spectrum, Angular Power Distribution, Origins and Directions distibution and a 3D Transform. + +```julia +Source(::Type{T} = Float64; + transform::Tr = Transform(), + spectrum::S = Spectrum.Uniform(), + origins::O = Origins.Point(), + directions::D = Directions.Constant(), + power::P = AngularPower.Lambertian(), + sourcenum::Integer = 0) where { + Tr<:Transform, + S<:Spectrum.AbstractSpectrum, + O<:Origins.AbstractOriginDistribution, + D<:Directions.AbstractDirectionDistribution, + P<:AngularPower.AbstractAngularPowerDistribution, + T<:Real} + +Source(transform::Tr, spectrum::S, origins::O, directions::D, power::P, ::Type{T} = Float64; sourcenum::Integer = 0) where { + Tr<:Transform, + S<:Spectrum.AbstractSpectrum, + O<:Origins.AbstractOriginDistribution, + D<:Directions.AbstractDirectionDistribution, + P<:AngularPower.AbstractAngularPowerDistribution, + T<:Real} +``` +""" +struct Source{ + T<:Real, + Tr<:Transform{T}, + S<:Spectrum.AbstractSpectrum{T}, + O<:Origins.AbstractOriginDistribution{T}, + D<:Directions.AbstractDirectionDistribution{T}, + P<:AngularPower.AbstractAngularPowerDistribution{T} +} <: AbstractSource{T} + transform::Tr + spectrum::S + origins::O + directions::D + power_distribution::P + sourcenum::Integer + + function Source( + ::Type{T} = Float64; + transform::Tr = Transform(), + spectrum::S = Spectrum.Uniform(), + origins::O = Origins.Point(), + directions::D = Directions.Constant(), + power::P = AngularPower.Lambertian(), + sourcenum::Integer = 0 + ) where { + Tr<:Transform, + S<:Spectrum.AbstractSpectrum, + O<:Origins.AbstractOriginDistribution, + D<:Directions.AbstractDirectionDistribution, + P<:AngularPower.AbstractAngularPowerDistribution, + T<:Real + } + new{T, Tr, S, O, D, P}(transform, spectrum, origins, directions, power, sourcenum) + end + + function Source( + transform::Tr, + spectrum::S, + origins::O, + directions::D, + power::P, + ::Type{T} = Float64; + sourcenum::Integer = 0 + ) where { + Tr<:Transform, + S<:Spectrum.AbstractSpectrum, + O<:Origins.AbstractOriginDistribution, + D<:Directions.AbstractDirectionDistribution, + P<:AngularPower.AbstractAngularPowerDistribution, + T<:Real + } + new{T, Tr, S, O, D, P}(transform, spectrum, origins, directions, power, sourcenum) + end +end + +# used to not generate new origin points if we can use last points - mainly to keep consistency of origin generation when randomness is involved +struct SourceGenerationState{T<:Real} + n::Integer + last_index_used::Integer + last_point_generated::Vec3{T} +end + +Base.length(s::Source) = length(s.origins) * length(s.directions) +Base.iterate(s::Source{T}) where {T<:Real} = iterate(s, SourceGenerationState(1, -1, zeros(Vec3{T}))) +Base.iterate(s::Source, state::SourceGenerationState) = state.n > length(s) ? nothing : generate(s, state) + +function Emitters.generate(s::Source{T}, n::Integer) where {T<:Real} + return generate(s, SourceGenerationState(n+1, -1, zeros(Vec3{T})))[1] +end + +function Emitters.generate(s::Source{T}, state::SourceGenerationState{T} = SourceGenerationState(-1, -1, zeros(Vec3{T}))) where {T<:Real} + # @info "New State Generation" + n::Integer = state.n - 1 + origin_index = floor(Integer, (n / length(s.directions))) + direction_index = mod(n, length(s.directions)) + + if (origin_index == state.last_index_used) + o = state.last_point_generated + else + o = generate(s.origins, origin_index) + o = s.transform * o + end + + # o = generate(s.origins, origin_index) + d = generate(s.directions, direction_index) + + # rotate the direction according to the orientation of the ray-originator + d = rotation(s.transform) * d + + ray = OpticSim.Ray(o, d) + power, wavelength = generate(s.spectrum) + + power = apply(s.power_distribution, s.transform, power, ray) + + # return (OpticSim.Ray(o, d), SourceGenerationState(state.n+1, origin_index, o)) + return (OpticSim.OpticalRay(ray, power, wavelength; sourcenum=s.sourcenum), SourceGenerationState(state.n+1, origin_index, o)) +end + +""" + CompositeSource{T} <: AbstractSource{T} + +This data-type represents the composite emitter (Source) which is constructed with a list of basic or composite emitters and a 3D Transform. + +```julia +CompositeSource(transform::Transform{T}, sources::Array) where {T<:Real} +``` +""" +struct CompositeSource{T} <: AbstractSource{T} + transform::Transform{T} + sources::Array + + uniform_length::Integer + total_length::Integer + start_indexes::Vector{Integer} + + function CompositeSource(transform::Transform{T}, sources::Array) where {T<:Real} + lens = [length(src) for src in sources] + size1 = findmin(lens)[1] + size2 = findmax(lens)[1] + if (size1 == size2) + uniform_length = size1 + start_indexes = Vector{Integer}() + total_length = length(sources) * uniform_length + # @info "Uniform Length $size1 total=$total_length" + else + uniform_length = -1 + + start_indexes = Vector{Integer}() + start = 0 + for s in sources + push!(start_indexes, start) + start = start + length(s) + end + push!(start_indexes, start) + # @info start_indexes + total_length = start + end + + new{T}(transform, sources, uniform_length, total_length, start_indexes) + end +end + +Base.length(s::CompositeSource) = s.total_length +Base.iterate(s::CompositeSource{T}) where {T<:Real} = iterate(s, SourceGenerationState(1, -1, zeros(Vec3{T}))) +Base.iterate(s::CompositeSource, state::SourceGenerationState) = state.n > length(s) ? nothing : generate(s, state) + +function Emitters.generate(s::CompositeSource{T}, n::Integer) where {T<:Real} + return generate(s, SourceGenerationState(n+1, -1, zeros(Vec3{T})))[1] +end + +function Emitters.generate(s::CompositeSource{T}, state::SourceGenerationState{T} = SourceGenerationState(-1, -1, zeros(Vec3{T}))) where {T<:Real} + # @info "New Composite State Generation" + n::Integer = state.n - 1 + + if (s.uniform_length != -1) + source_index = floor(Integer, (n / s.uniform_length)) + ray_index = mod(n, s.uniform_length) + # @info "New Composite State Generation: source=$source_index ray=$ray_index ($(length(s.sources)))" + else + n = mod(n, s.total_length) + + ## perform a binary search to find out which source can generate the spesific ray + low = 1 + high = length(s.start_indexes) - 1 # last cell in this vector is a dummy cell containing the index after the last ray + source_index = -1 + while (low <= high) + mid = floor(Integer, (low + high) / 2) + + if (n < s.start_indexes[mid]) + high = mid - 1 + elseif (n >= s.start_indexes[mid+1]) + low = mid + 1 + else + source_index = mid + break; + end + end + @assert source_index != -1 + # @info "Bin Search n=$n source_index=$source_index" + source_index = source_index - 1 + ray_index = n - s.start_indexes[source_index+1] + end + + ray, new_state = generate(s.sources[source_index+1], SourceGenerationState(ray_index+1, state.last_index_used, state.last_point_generated)) + ray = s.transform * ray + return (ray, SourceGenerationState(state.n+1, new_state.last_index_used, new_state.last_point_generated)) +end + +end # module Sources diff --git a/src/Optical/Emitters/Spectrum.jl b/src/Optical/Emitters/Spectrum.jl new file mode 100644 index 000000000..2f8a8d7d0 --- /dev/null +++ b/src/Optical/Emitters/Spectrum.jl @@ -0,0 +1,125 @@ +module Spectrum +export Uniform, DeltaFunction, Measured + +using ....OpticSim +using ...Emitters +using DataFrames +using Distributions + +const UNIFORMSHORT = 0.450 #um +const UNIFORMLONG = 0.680 #um + +abstract type AbstractSpectrum{T<:Real} end + +""" + Uniform{T} <: AbstractSpectrum{T} + +Encapsulates a flat spectrum range which is sampled uniformly. Unless stated diferrently, the range used will be 450nm to 680nm. + +```julia +Uniform(low::T, high::T) where {T<:Real} +Uniform(::Type{T} = Float64) where {T<:Real} +``` +""" +struct Uniform{T} <: AbstractSpectrum{T} + low_end::T + high_end::T + + # user defined range of spectrum + function Uniform(low::T, high::T) where {T<:Real} + return new{T}(low, high) + end + + # with no specific range we will use the constants' values + function Uniform(::Type{T} = Float64) where {T<:Real} + return new{T}(UNIFORMSHORT, UNIFORMLONG) + end +end + +Emitters.generate(s::Uniform{T}) where {T<:Real} = (one(T), rand(Distributions.Uniform(s.low_end, s.high_end))) + +""" + DeltaFunction{T} <: AbstractSpectrum{T} + +Encapsulates a constant spectrum. + +```julia +DeltaFunction{T<:Real} +``` +""" +struct DeltaFunction{T} <: AbstractSpectrum{T} + λ::T +end + +Emitters.generate(s::DeltaFunction{T}) where {T<:Real} = (one(T), s.λ) + +""" + Measured{T} <: AbstractSpectrum{T} + +Encapsulates a measured spectrum to compute emitter power. Create spectrum by reading CSV files. +Evaluate spectrum at arbitrary wavelength with [`spectrumpower`](@ref) (**more technical details coming soon**) + +```julia +Measured(samples::DataFrame) +``` +""" +struct Measured{T} <: AbstractSpectrum{T} + low_wave_length::Integer + high_wave_length::Integer + wave_length_step::Integer + power_samples::Vector{T} + + function Measured(samples::DataFrame) + colnames = names(samples) + @assert "Wavelength" in colnames + @assert "Power" in colnames + wavelengths = samples[!, :Wavelength] + @assert eltype(wavelengths) <: Integer + + power::Vector{T} where {T<:Real} = samples[!, :Power] # no missing values allowed and must be real numbers + maxpower = maximum(power) + T = eltype(power) + p = sortperm(wavelengths) + wavelengths = wavelengths[p] + power = power[p] ./ maxpower #make sure power values have the same sorted permutation as wavelengths normalized with maximum power equal to 1 + λmin = wavelengths[p[1]] + λmax = wavelengths[p[end]] + step = wavelengths[2] - wavelengths[1] + + for i in 3:length(wavelengths) + @assert wavelengths[i] - wavelengths[i - 1] == step # no missing values allowed and step size must be consistent + end + + return new{T}(λmin, λmax, step, power) + end +end + +function Emitters.generate(spectrum::Measured{T}) where {T<:Real} + λ = rand(Distributions.Uniform(convert(T, spectrum.low_wave_length), convert(T, spectrum.high_wave_length))) + power = spectrumpower(spectrum, λ) + if power === nothing #added this condition because the compiler was allocating if just returned values directly. + return (nothing, nothing) + else + return (power, λ / convert(T, 1000)) #convert from nm to um + end +end + +"""expects wavelength in nm not um""" +function spectrumpower(spectrum::Measured{T}, λ::T) where {T<:Real} + if λ < spectrum.low_wave_length || λ > spectrum.high_wave_length + return nothing + end + + lowindex = floor(Integer, (λ - spectrum.low_wave_length)) ÷ spectrum.wave_length_step + 1 + + if lowindex == length(spectrum.power_samples) + return convert(T, spectrum.power_samples[end]) + else + highindex = lowindex + 1 + α = mod(λ, spectrum.wave_length_step) / spectrum.wave_length_step + + return convert(T, (1 - α) * spectrum.power_samples[lowindex] + α * spectrum.power_samples[highindex]) + end +end + +end # module Spectrum diff --git a/src/Optical/Optical.jl b/src/Optical/Optical.jl index eaf31199d..1581839a3 100644 --- a/src/Optical/Optical.jl +++ b/src/Optical/Optical.jl @@ -24,7 +24,7 @@ include("OpticalRay.jl") include("Emitter.jl") -include("Emitter2.jl") +include("Emitters/Emitters.jl") include("Fresnel.jl") include("Paraxial.jl") include("Grating.jl") diff --git a/src/Vis/Emitters.jl b/src/Vis/Emitters.jl new file mode 100644 index 000000000..7a04bd076 --- /dev/null +++ b/src/Vis/Emitters.jl @@ -0,0 +1,154 @@ +using ..OpticSim, .Geometry +using LinearAlgebra +using Distributions +using StaticArrays + +import Makie +import Makie.AbstractPlotting +import Makie.AbstractPlotting.MakieLayout + +using .Emitters +using .Emitters.Spectrum +using .Emitters.Directions +using .Emitters.Origins +using .Emitters.AngularPower +using .Emitters.Sources + +const ARRROW_LENGTH = 0.5 +const ARRROW_SIZE = 0.01 +const MARKER_SIZE = 15 + + +#------------------------------------- +# draw debug information - local axes and positions +#------------------------------------- +function maybe_draw_debug_info(scene::MakieLayout.LScene, o::Origins.AbstractOriginDistribution; transform::Geometry.Transform = Transform(), debug::Bool=false, kwargs...) where {T<:Real} + + dir = forward(transform) + uv = SVector{3}(right(transform)) + vv = SVector{3}(up(transform)) + pos = origin(transform) + + if (debug) + # this is a stupid hack to force makie to render in 3d - for some scenes, makie decide with no apperent reason to show in 2d instead of 3d + AbstractPlotting.scatter!(scene, [pos[1], pos[1]+0.1], [pos[2], pos[2]+0.1], [pos[3], pos[3]+0.1], color=:red, markersize=0) + + # draw the origin and normal of the surface + Makie.scatter!(scene, pos, color=:blue, markersize = MARKER_SIZE * visual_size(o)) + + # normal + arrow_start = pos + arrow_end = dir * ARRROW_LENGTH * visual_size(o) + Makie.arrows!(scene.scene, [AbstractPlotting.Point3f0(arrow_start)], [AbstractPlotting.Point3f0(arrow_end)], arrowsize=ARRROW_SIZE * visual_size(o), arrowcolor=:blue) + arrow_end = uv * 0.5 * ARRROW_LENGTH * visual_size(o) + Makie.arrows!(scene.scene, [AbstractPlotting.Point3f0(arrow_start)], [AbstractPlotting.Point3f0(arrow_end)], arrowsize= 0.5 * ARRROW_SIZE * visual_size(o), arrowcolor=:red) + arrow_end = vv * 0.5 * ARRROW_LENGTH * visual_size(o) + Makie.arrows!(scene.scene, [AbstractPlotting.Point3f0(arrow_start)], [AbstractPlotting.Point3f0(arrow_end)], arrowsize= 0.5 * ARRROW_SIZE * visual_size(o), arrowcolor=:green) + + # draw all the samples origins + positions = map(x -> transform*x, collect(o)) + positions = collect(AbstractPlotting.Point3f0, positions) + Makie.scatter!(scene, positions, color=:green, markersize = MARKER_SIZE * visual_size(o)) + + # positions = collect(AbstractPlotting.Point3f0, o) + # Makie.scatter!(scene, positions, color=:green, markersize = MARKER_SIZE * visual_size(o)) + end + +end + + +#------------------------------------- +# draw point origin +#------------------------------------- +# function OpticSim.Vis.draw!(scene::MakieLayout.LScene, o::Origins.Point; transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real} +function OpticSim.Vis.draw!(scene::MakieLayout.LScene, o::Origins.Point; transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real} + maybe_draw_debug_info(scene, o; transform=transform, kwargs...) +end + +#------------------------------------- +# draw RectGrid and RectUniform origins +#------------------------------------- +function OpticSim.Vis.draw!(scene::MakieLayout.LScene, o::Union{Origins.RectGrid, Origins.RectUniform}; transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real} + dir = forward(transform) + uv = SVector{3}(right(transform)) + vv = SVector{3}(up(transform)) + pos = origin(transform) + + # @info "RECT: transform $(pos)" + + plane = OpticSim.Plane(dir, pos) + rect = OpticSim.Rectangle(plane, o.width / 2, o.height / 2, uv, vv) + + OpticSim.Vis.draw!(scene, rect; kwargs...) + + maybe_draw_debug_info(scene, o; transform=transform, kwargs...) +end + + +#------------------------------------- +# draw hexapolar origin +#------------------------------------- +function OpticSim.Vis.draw!(scene::MakieLayout.LScene, o::Origins.Hexapolar; transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real} + dir = forward(transform) + uv = SVector{3}(right(transform)) + vv = SVector{3}(up(transform)) + pos = origin(transform) + + plane = OpticSim.Plane(dir, pos) + ellipse = OpticSim.Ellipse(plane, o.halfsizeu, o.halfsizev, uv, vv) + + OpticSim.Vis.draw!(scene, ellipse; kwargs...) + + maybe_draw_debug_info(scene, o; transform=transform, kwargs...) +end + +#------------------------------------- +# draw source +#------------------------------------- +function OpticSim.Vis.draw!(scene::MakieLayout.LScene, s::Sources.Source{T}; parent_transform::Geometry.Transform = Transform(), debug::Bool=false, kwargs...) where {T<:Real} + + OpticSim.Vis.draw!(scene, s.origins; transform=parent_transform * s.transform, debug=debug, kwargs...) + + if (debug) + m = zeros(T, length(s), 7) + for (index, optical_ray) in enumerate(s) + ray = OpticSim.ray(optical_ray) + ray = parent_transform * ray + m[index, 1:7] = [ray.origin... ray.direction... OpticSim.power(optical_ray)] + end + + m[:, 4:6] .*= m[:, 7] * ARRROW_LENGTH * visual_size(s.origins) + + # Makie.arrows!(scene, [Makie.Point3f0(origin(ray))], [Makie.Point3f0(rayscale * direction(ray))]; kwargs..., arrowsize = min(0.05, rayscale * 0.05), arrowcolor = color, linecolor = color, linewidth = 2) + color = :yellow + Makie.arrows!(scene, m[:,1], m[:,2], m[:,3], m[:,4], m[:,5], m[:,6]; kwargs..., arrowcolor = color, linecolor = color, linewidth = 2, arrowsize=ARRROW_SIZE * visual_size(s.origins)) + end + + # for ray in o + # OpticSim.Vis.draw!(scene, ray) + # end +end + +#------------------------------------- +# draw optical rays +#------------------------------------- +function OpticSim.Vis.draw!(scene::MakieLayout.LScene, rays::AbstractVector{OpticSim.OpticalRay{T, 3}}; kwargs...) where {T<:Real} + m = zeros(T, length(rays)*2, 3) + for (index, optical_ray) in enumerate(rays) + ray = OpticSim.ray(optical_ray) + m[(index-1)*2+1, 1:3] = [origin(optical_ray)...] + m[(index-1)*2+2, 1:3] = [(OpticSim.origin(optical_ray) + OpticSim.direction(optical_ray) * 1 * OpticSim.power(optical_ray))... ] + end + + color = :green + Makie.linesegments!(scene, m[:,1], m[:,2], m[:,3]; kwargs..., color = color, linewidth = 2, ) +end + +#------------------------------------- +# draw composite source +#------------------------------------- +function OpticSim.Vis.draw!(scene::MakieLayout.LScene, s::Sources.CompositeSource{T}; parent_transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real} + for source in s.sources + OpticSim.Vis.draw!(scene, source; parent_transform=parent_transform*s.transform, kwargs...) + end +end diff --git a/src/Vis/Vis.jl b/src/Vis/Vis.jl new file mode 100644 index 000000000..df54f4e63 --- /dev/null +++ b/src/Vis/Vis.jl @@ -0,0 +1,12 @@ +module Vis + +# If using precompiled system image (which we always are) you have to run AbstractPlotting.__init__() after loading Makie +# during precompilation, the display stack gets shuffled around such that the Makie display does not take priority. +# See https://discourse.julialang.org/t/makie-doesnt-display-plot-when-using-a-custom-julia-sysimage/38515. +__init__() = AbstractPlotting.__init__() + +include("Visualization.jl") +include("Emitters.jl") + +end # module Vis +export Vis diff --git a/src/Visualization/Visualization.jl b/src/Vis/Visualization.jl similarity index 99% rename from src/Visualization/Visualization.jl rename to src/Vis/Visualization.jl index 6c77fdd58..14838b4cf 100644 --- a/src/Visualization/Visualization.jl +++ b/src/Vis/Visualization.jl @@ -20,11 +20,9 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE -module Vis - using ..OpticSim using ..OpticSim: euclideancontrolpoints, evalcsg, vertex, makiemesh, detector, centroid, origingen, lower, upper, intervals, α -using ..Geometry +using ..OpticSim.Geometry using Unitful using ImageView @@ -40,11 +38,6 @@ import Plots import Luxor using FileIO -# If using precompiled system image (which we always are) you have to run AbstractPlotting.__init__() after loading Makie -# during precompilation, the display stack gets shuffled around such that the Makie display does not take priority. -# See https://discourse.julialang.org/t/makie-doesnt-display-plot-when-using-a-custom-julia-sysimage/38515. -__init__() = AbstractPlotting.__init__() - ############################################################################# function drawcurve!(canvassize::Int, curve::Spline{P,S,N,M}, numpoints::Int; linewidth = 0.5, curvecolor = RGB(1, 0, 1), controlpointcolor = RGB(0, 0, 0), controlpointsize = 5, controlpolygoncolor = RGB(0, 0, 0)) where {P,S,N,M} @@ -1055,6 +1048,3 @@ function eyebox_eval_planar(assembly::LensAssembly{T}, raygen::OpticalRayGenerat Vis.show(out_image) end end - -end # module Vis -export Vis