update normal calculation

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "98fe4994922ae0460c4ee5f09975a0b15c7fcff3"
+project_hash = "7deb69ebd2bc655c0098fd7f47a66961056e436c"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -241,6 +241,12 @@ version = "0.18.20"
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
+[[deps.Dbus_jll]]
+deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "fc173b380865f70627d7dd1190dc2fce6cc105af"
+uuid = "ee1fde0b-3d02-5ea6-8484-8dfef6360eab"
+version = "1.14.10+0"
+
 [[deps.DelimitedFiles]]
 deps = ["Mmap"]
 git-tree-sha1 = "9e2f36d3c96a820c678f2f1f1782582fcf685bae"
@@ -360,7 +366,7 @@ uuid = "559328eb-81f9-559d-9380-de523a88c83c"
 version = "1.0.14+0"
 
 [[deps.GLFW_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll", "xkbcommon_jll"]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll", "libdecor_jll", "xkbcommon_jll"]
 git-tree-sha1 = "3f74912a156096bd8fdbef211eff66ab446e7297"
 uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
 version = "3.4.0+0"
@@ -828,6 +834,12 @@ git-tree-sha1 = "fa7fd067dca76cadd880f1ca937b4f387975a9f5"
 uuid = "d3a379c0-f9a3-5b72-a4c0-6bf4d2e8af0f"
 version = "2.16.0+0"
 
+[[deps.Loess]]
+deps = ["Distances", "LinearAlgebra", "Statistics", "StatsAPI"]
+git-tree-sha1 = "a113a8be4c6d0c64e217b472fb6e61c760eb4022"
+uuid = "4345ca2d-374a-55d4-8d30-97f9976e7612"
+version = "0.6.3"
+
 [[deps.LogExpFunctions]]
 deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
 git-tree-sha1 = "a2d09619db4e765091ee5c6ffe8872849de0feea"
@@ -1052,6 +1064,12 @@ git-tree-sha1 = "0fac6313486baae819364c52b4f483450a9d793f"
 uuid = "5432bcbf-9aad-5242-b902-cca2824c8663"
 version = "0.5.12"
 
+[[deps.Pango_jll]]
+deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "FriBidi_jll", "Glib_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "cb5a2ab6763464ae0f19c86c56c63d4a2b0f5bda"
+uuid = "36c8627f-9965-5494-a995-c6b170f724f3"
+version = "1.52.2+0"
+
 [[deps.Parameters]]
 deps = ["OrderedCollections", "UnPack"]
 git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
@@ -1756,6 +1774,12 @@ deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
 version = "5.8.0+1"
 
+[[deps.libdecor_jll]]
+deps = ["Artifacts", "Dbus_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pango_jll", "Wayland_jll", "xkbcommon_jll"]
+git-tree-sha1 = "9bf7903af251d2050b467f76bdbe57ce541f7f4f"
+uuid = "1183f4f0-6f2a-5f1a-908b-139f9cdfea6f"
+version = "0.2.2+0"
+
 [[deps.libevdev_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "141fe65dc3efabb0b1d5ba74e91f6ad26f84cc22"

Project.toml

@@ -4,13 +4,17 @@ authors = ["Jessica Kent", "Ted Marku"]
 version = "1.0.0-DEV"
 
 [deps]
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 ImageSegmentation = "80713f31-8817-5129-9cf8-209ff8fb23e1"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Loess = "4345ca2d-374a-55d4-8d30-97f9976e7612"
 MultipleScattering = "80a8ab25-5750-5d93-a6d7-4adc97cdd5fb"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]

docs/examples/a-maze-ing_example.jl

@@ -0,0 +1,73 @@
+using RayTracing
+using Plots
+
+# we can use the image_to_domain function to convert an image to a domain
+domain = image_to_domain("dev/RayTracing.jl/docs/examples/maze.png")
+
+# this step ensures that the boundaries are in the right place
+
+domain = Domain(size(domain.domain,1), size(domain.domain, 2), iszero.(domain.domain), -domain.boundary_normals)
+
+# plot domain!
+heatmap(domain.domain)
+
+dvdx          = 0.0;
+dvdz          = 0.0;
+base_velocity = 1.0;        
+
+#define a velocity gradient.
+
+velocity_grads = VelocityGradients(dvdx,dvdz,base_velocity); 
+
+# inialise the fields
+fields = initialise_fields(domain, grad_to_field(domain, velocity_grads), zeros(512,512))
+
+# slowness map
+
+heatmap(fields.slowness_field)
+
+#define Parameters
+
+t_max       = 1000;
+no_rays     = 100;
+angle_range = [0, 360];
+params      = RayTracing.initialise_parameters(t_max,no_rays,angle_range,fields);
+
+# initial ray parameters
+
+x_0       = 256.0;  
+z_0       = 215.0;
+position  = [x_0 , z_0];
+amp_0     = 1;
+theta     = params.angles[1];
+u_initial = RayTracing.interpolate_field(fields.slowness_field, x_0,z_0);
+
+# initials raays and ray data
+ray = Ray([513.0, 256.0], rand(2), 1.0, rand(2))
+RayTracing.in_domain(ray, domain)
+results = Vector{RayTracing.RaySimulationResult}(undef, size(params.angles,1))
+no_samples  = 1000
+
+for (i,theta) in enumerate(params.angles)
+    ICs      = InitialConditions(position,amp_0,theta,u_initial)
+    ray_data = initialise_ray_data(domain, params,no_samples)
+    sim = RaySimulation(ICs,domain,params,fields,ray_data)
+    results[i] = calculate_ray(sim, ReflectionType())
+end
+
+source_data = [r.ray_data for r in results]
+
+plot_ray(domain,source_data)
+scatter!([x_0],[z_0])
+my_palette = palette(:Oranges, 100)
+anim = @animate for i in 1:10:1000
+    x = [s.position[1:i,1] for s in source_data]
+    z = [s.position[1:i,2] for s in source_data]
+    heatmap(domain.domain, c=:blues, leg=false)
+    for i in 1:length(x)
+        plot!(x[i], z[i],yflip=true, xlims=(0,512), ylims=(0,512), palette=my_palette)
+    end
+    scatter!([x_0],[z_0])
+end
+
+gif(anim,"dev/RayTracing.jl/docs/imgs/a-maze-ing.gif", fps = 7)

src/RayTracing.jl

@@ -41,6 +41,7 @@ using MultipleScattering
 using LinearAlgebra
 using Images, ImageSegmentation
 using Statistics
+using FFTW
 
 include("types.jl")
 include("utils.jl")

src/boundary_logic.jl

@@ -75,4 +75,14 @@ function scatter_logic(domain::Domain, fields::DomainFields, ray::Ray, params::P
     end
 
     return rays
-end
+end
+
+function in_boundary(normal::Vector, direction::Vector)
+    is_in_boundary = false
+
+    if dot(normal, direction) == 0.0
+        is_in_boundary = true
+    end
+
+    return is_in_boundary
+end

src/geometry.jl

@@ -1,7 +1,6 @@
 struct Boundary
-    indices::Vector{CartesianIndex{2}}
+    boundary_points::AbstractArray
     id::Int
-    normals::Array{Float64, 3}
 end
 
 function identify_regions(filename::String, threshold::Float64)
@@ -14,15 +13,36 @@ function identify_regions(filename::String, threshold::Float64)
     return id_map, domain_ids
 end
 
-function get_region(domain::Matrix{Float64}, id::Int)
-    inds = findall(iszero.(domain .- id), domain)
+function smooth_domain(domain::Matrix)
 
-    return Region(inds)
+    x_dims = size(domain,2)
+    z_dims = size(domain,1)
+    area = x_dims*z_dims
+
+    α = 2pi*0.5^2/(9*(var(domain)*0.01*area))
+    sample_rate = Int(round(sqrt(area)/2))
+
+    gaussians1 = [[exp((-(norm([(i-k)/z_dims,j/x_dims]))^2)/2α) for i in 1:z_dims, j in 1:x_dims] for k in 0:sample_rate:z_dims]
+
+    gaussians2 = [[exp((-(norm([(i-k)/z_dims,(j-sample_rate)/x_dims]))^2)/2α) for i in 1:z_dims, j in 1:x_dims] for k in 0:sample_rate:z_dims]
+
+    gaussians3 = [[exp((-(norm([(i-k)/z_dims,(j-2*sample_rate)/x_dims]))^2)/2α) for i in 1:z_dims, j in 1:x_dims] for k in 0:sample_rate:z_dims]
+
+    gaussians = sum(gaussians1)+sum(gaussians2)+sum(gaussians3)
+
+    dom_spectrum = fft(domain)
+
+    pointwise = [dom_spectrum[i,j]*gaussians[i,j] for i in 1:z_dims, j in 1:x_dims]
+
+    smoothed = ifft(pointwise)
+
+    return abs.(smoothed)
 end
 
 function get_normals(boundaries::Matrix)
     nz,nx = size(boundaries)
     boundary_normals = zeros(Float64,nz,nx,2)
+    boundaries = smooth_domain(boundaries)
     for i = 2:nx-1, j = 2:nz-1
         grad_z = 0; grad_x = 0
         for ii = (i-1):(i+1)
@@ -43,19 +63,47 @@ function get_normals(boundaries::Matrix)
     return boundary_normals
 end
 
-function get_boundaries(domain::Matrix{Float64})
-
-    domain_field = domain.domain
-    normals = domain.boundary_normals
-    boundaries_field = (abs2.(normals[:,:,1])+abs2.(normals[:,:,2]))/2
-    boundary_inds = findall(!iszero, boundaries_field)
+function center_of_mass(matrix::AbstractArray)
+    total_mass = sum(matrix)
+    xs = 1:size(matrix,2)
+    zs = 1:size(matrix,1)
+    weighted_positions = [matrix[i,j]*[j,i] for i in zs, j in xs]
+
+    center = sum(weighted_positions)/total_mass
+
+    return center
+end
+
+
+function is_closed(matrix::AbstractArray)
+    closed = false
+    center = Int.(round.(center_of_mass(matrix)))
+    x_slice = matrix[:,center[1]]
+    z_slice = matrix[center[2],:]
+    first_zero_x = findfirst(iszero, x_slice)
+    last_one_x = findlast(isone, x_slice)
+    first_zero_z = findfirst(iszero, z_slice)
+    last_one_z = findlast(isone, z_slice)
 
-    for inds in boundary_inds
-        boundaries_field[inds] *= domain_field[inds]
+    if first_zero_x + (size(matrix,2)-last_one_x) -1 != size(matrix,2) && first_zero_z + (size(matrix,1)-last_one_z) -1 != size(matrix,1)
+        closed = true
     end
+    return closed
+end
 
+function march(direction::Vector, center::Vector, matrix::AbstractArray)
+    position = center + direction
+    at_boundary = false
 
+    while at_boundary == false
+        submatrix, nodex, nodez = get_submatrix(matrix, position[2], position[1])
 
+        if sum(submatrix) != 9
+            at_boundary = true
+        else
+            position += direction
+        end
+    end
 
-
+    return position
 end

src/ray_physics.jl

@@ -1,9 +1,9 @@
 # ray physics
 
 function snells_law(domain::Domain, ray::Ray, media::Vector{Acoustic{Float64, 2}})
-    position = Int.(round.(ray.position))
-    normal = domain.boundary_normals[position[2], position[1],:]
-    tangent = [-normal[2], normal[1]]
+
+    normal = [interpolate_field(domain.boundary_normals[:,:,1], ray.position[1], ray.position[2]), interpolate_field(domain.boundary_normals[:,:,2], ray.position[1], ray.position[2])]
+    tangent = [normal[2], -normal[1]]
     cs = [real(m.c) for m in media]
 
     #calculate angle from normal
@@ -23,9 +23,8 @@ end
 function snells_law(domain::Domain, ray::Ray)
     #Only returns reflected direction
 
-    position = Int.(round.(ray.position))
-    normal = domain.boundary_normals[position[2], position[1],:]
-    tangent = [-normal[2], normal[1]]
+    normal = [interpolate_field(domain.boundary_normals[:,:,1], ray.position[1], ray.position[2]), interpolate_field(domain.boundary_normals[:,:,2], ray.position[1], ray.position[2])]
+    tangent = [normal[2], -normal[1]]
 
     #calculate angle from normal
     cosθ1 = dot(normal, ray.direction)
@@ -61,11 +60,11 @@ function reflect_ray(domain::Domain,params::Parameters,fields::DomainFields,ray:
     x = ray.position[1];
     z = ray.position[2]
 
+
 
     v2 = snells_law(domain, ray)
 
     dx = v2[1]; dz = v2[2];
-    position += v2;
     u = interpolate_field(fields.slowness_field, x, z)
     sx = (dx/dt)*u^2; sz = (dz/dt)*u^2;
 

src/ray_tracing.jl

@@ -159,6 +159,7 @@ function calculate_ray(sim::RaySimulation, boundary_type::ReflectionType)
     steps = size(params.time, 1)
     distance_travelled = 0.0
     amps = ray_data.amplitudes
+    reflected = false
     for i = 2:steps-1
         is_in_domain = in_domain(ray, domain)
 
@@ -167,17 +168,27 @@ function calculate_ray(sim::RaySimulation, boundary_type::ReflectionType)
             nodex = nodes[1]
             nodez = nodes[2]
             amps[nodez,nodex] += ray.amp
-            at_boundary = boundary_detect(domain, ray)
+            at_boundary = boundary_detect(domain, ray; transmitted = reflected)
+
+            nodes = Int.(round.(ray.position))
+            nodex = nodes[1]
+            nodez = nodes[2]
             if at_boundary==true
                 nodes = Int.(round.(ray.position))
                 nodex = nodes[1]
                 nodez = nodes[2]
                 ray = reflect_ray(domain,params,fields,ray)
-            else
-                old_position = ray.position
-                ray = advance_ray(ICs,params,fields,ray, distance_travelled)
-                distance_travelled += norm(ray.position - old_position)
+                reflected = true
+                println("position:", ray.position)
+                println("direction:", ray.direction)
+                println("normal:", domain.boundary_normals[nodez, nodex, :])
             end
+
+            old_position = ray.position
+            ray = advance_ray(ICs,params,fields,ray, distance_travelled)
+            distance_travelled += norm(ray.position - old_position)
+            reflected = false
+
         else
             # Ray has left the Domain
             break