Add level set geometry example

examples/level-set-example.jl

@@ -0,0 +1,50 @@
+# # Raytracing a Level set geometry
+# A level set geoemtry is defined by a constraint equations $f(x,y,z)=0$.
+# We will ray trace an example parabaloid geometry in this example as a simple geometric jet model.
+using Krang
+import CairoMakie as GLMk
+using GeometryBasics
+using FileIO
+
+# Lets create a camera with a screen of 20Mx20M at a resolution of 200x200 pixels for a high spin black hole.
+metric = Krang.Kerr(0.9) # Kerr metric with a spin of 0.99
+θo = 89 / 180 * π # Inclination angle of the observer
+ρmax = 20.0 # Horizontal and Vertical extent of the screen
+sze = 200 # Resolution of the screen is sze x sze
+
+camera = Krang.SlowLightIntensityCamera(metric, θo, -ρmax, ρmax, -ρmax, ρmax, sze)
+
+# We will define a parabaloid geometry to the set of all points satisfied by the equation $\cos(θ)=1-(r/r_h)^n$.
+struct Parabaloid{T} <: Krang.AbstractLevelSetGeometry{T}
+    rh::T
+    index::T
+end
+function (geometry::Parabaloid)(x,y,z)
+    r = sqrt(x^2+y^2+z^2)
+    return 1-(r/geometry.rh)^geometry.index*(1-z/r)
+end
+
+# The jet will be emit a constant intensity whose physics we define in the `XMaterial`.
+# [!NOTE] We are ignoring relativistic effects in this example. 
+struct XMaterial <: Krang.AbstractMaterial end
+function (mat::XMaterial)(pix, intersection) where T
+    return 1.0
+end
+
+# We will ray trace the geometry and plot the image.
+parabaloid = Parabaloid(Krang.horizon(metric), 1.0)
+
+fig = GLMk.Figure();
+ax = GLMk.Axis(fig[1, 1], aspect = 1)
+
+
+intersections = raytrace(camera, Krang.Mesh(parabaloid, XMaterial()), res = 1_00)
+
+# And plot the image with GLMakie, 
+
+GLMk.heatmap!(ax, intersections, colormap = :afmhot);
+fig
+
+GLMk.save("geometric_jet.png", fig)
+
+# ![Simple geometric jet model](geometric_jet.png)

examples/neural-net-example.jl

@@ -35,7 +35,7 @@ function (m::ImageModel)(x, ps, st)
     emission_vals = zeros(Float64, 1, sze * sze)
     for n = 0:1
         for i = 1:sze
-            Threads.@threads for j = 1:sze
+            for j = 1:sze
                 pix = pixels[i+(j-1)*sze]
                 α, β = Krang.screen_coordinate(pix)
                 T = typeof(α)
@@ -108,7 +108,7 @@ heatmap!(
     Axis(fig[1, 2], aspect = 1, title = "Image Model (Lensed Emission Model)"),
     received_intensity,
 )
-save("emission_model_and_target_model.png", fig)
+CairoMakie.save("emission_model_and_target_model.png", fig)
 
 # ![image](emission_model_and_target_model.png)
 
@@ -154,7 +154,7 @@ heatmap!(
     received_intensity,
     colormap = :afmhot,
 )
-save("emission_model_and_image_model.png", fig)
+CairoMakie.save("emission_model_and_image_model.png", fig)
 
 # ![image](emission_model_and_image_model.png)
 
@@ -205,30 +205,32 @@ received_intensity, st =
 acc_intensity, st = ((x) -> (reshape(x[1], sze, sze), x[2]))(image_model(pixels, ps, st))
 loss_function(pixels, target_img, ps, st)
 loss_function(pixels, target_img, ps_trained, st_trained)
+
 using Printf
-begin
-    fig = Figure(size = (700, 300))
-    heatmap!(
-        Axis(fig[1, 1], aspect = 1, title = "Target Image"),
-        reshape(target_img, sze, sze),
-    )
-    heatmap!(
-        Axis(
-            fig[1, 2],
-            aspect = 1,
-            title = "Starting State (loss=$(@sprintf("%0.2e", loss_function(pixels, target_img, ps, st)[1])))",
-        ),
-        acc_intensity,
-    )
-    heatmap!(
-        Axis(
-            fig[1, 3],
-            aspect = 1,
-            title = "Fitted State (loss=$(@sprintf("%0.2e", loss_function(pixels, target_img, ps_trained, st_trained)[1])))",
-        ),
-        received_intensity,
-    )
-    save("neural_net_results.png", fig)
-end
+
+fig = Figure(size = (700, 300))
+heatmap!(
+    Axis(fig[1, 1], aspect = 1, title = "Target Image"),
+    reshape(target_img, sze, sze),
+)
+heatmap!(
+    Axis(
+        fig[1, 2],
+        aspect = 1,
+        title = "Starting State (loss=$(@sprintf("%0.2e", loss_function(pixels, target_img, ps, st)[1])))",
+    ),
+    acc_intensity,
+)
+heatmap!(
+    Axis(
+        fig[1, 3],
+        aspect = 1,
+        title = "Fitted State (loss=$(@sprintf("%0.2e", loss_function(pixels, target_img, ps_trained, st_trained)[1])))",
+    ),
+    received_intensity,
+)
+fig
+
+save("neural_net_results.png", fig)
 
 # ![image](neural_net_results.png)