diff --git a/.gitignore b/.gitignore index a45d4bb..a31e38d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,9 +1,12 @@ +# Editors +.vscode + # Files generated by invoking Julia with --code-coverage *.jl.cov *.jl.*.cov *.obj - +examples/**.png # Files generated by invoking Julia with --track-allocation *.jl.mem @@ -25,3 +28,14 @@ docs/site/ # committed for packages, but should be committed for applications that require a static # environment. Manifest.toml + +# Misc. +*.tmp +*.temp + +*sync-conflict* +src/debug_script.jl +examples/personal +examples/*.png + + diff --git a/.gitignore.orig b/.gitignore.orig new file mode 100644 index 0000000..6000d9f --- /dev/null +++ b/.gitignore.orig @@ -0,0 +1,54 @@ +# Editors +.vscode + +# Files generated by invoking Julia with --code-coverage +*.jl.cov +*.jl.*.cov + +*.obj +it*.png +loss*.png +change*.png +delta*.png +error*.png +img*.png +orig*.png + + +# Files generated by invoking Julia with --track-allocation +*.jl.mem + +# System-specific files and directories generated by the BinaryProvider and BinDeps packages +# They contain absolute paths specific to the host computer, and so should not be committed +deps/deps.jl +deps/build.log +deps/downloads/ +deps/usr/ +deps/src/ + +# Build artifacts for creating documentation generated by the Documenter package +docs/build* +docs/site/ + +# File generated by Pkg, the package manager, based on a corresponding Project.toml +# It records a fixed state of all packages used by the project. As such, it should not be +# committed for packages, but should be committed for applications that require a static +# environment. +Manifest.toml + +# Misc. +*.tmp +*.temp + +*sync-conflict* +src/debug_script.jl +examples/personal +examples/*.png + + +<<<<<<< HEAD + +src/debug_script.jl +examples/personal +======= +>>>>>>> 9473c949cf07d0e9431c75b57947a2304ac33a97 diff --git a/Project.toml b/Project.toml index 61a1bcf..ad64858 100644 --- a/Project.toml +++ b/Project.toml @@ -1,15 +1,30 @@ name = "CausticsEngineering" uuid = "e17b973c-3eae-4e30-9cf2-ee583e3dfd43" authors = ["Matt Ferraro", "Emmanuel Rialland"] -version = "0.1.0" +version = "0.2.0" [deps] +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +Debugger = "31a5f54b-26ea-5ae9-a837-f05ce5417438" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" +FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" +GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" +LibGit2 = "76f85450-5226-5b5a-8eaa-529ad045b433" +MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118" +MeshViz = "9ecf9c4f-6e5a-4b5e-83ae-06f2c7a661d8" +Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa" +OhMyREPL = "5fb14364-9ced-5910-84b2-373655c76a03" +PkgTemplates = "14b8a8f1-9102-5b29-a752-f990bacb7fe1" +Plotly = "58dd65bb-95f3-509e-9936-c39a10fdeae7" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" +Weave = "44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9" [compat] julia = "1" diff --git a/docs/make.jl b/docs/make.jl index f3664cc..b54d21f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,7 +1,7 @@ push!(LOAD_PATH, "../src/") using CausticsEngineering -using Documenter +using Documenter, DocStringExtensions # # HTML docs diff --git a/examples/loss_it1.png b/examples/loss_it1.png deleted file mode 100644 index 2a0709b..0000000 Binary files a/examples/loss_it1.png and /dev/null differ diff --git a/examples/loss_it2.png b/examples/loss_it2.png deleted file mode 100644 index 3188ba7..0000000 Binary files a/examples/loss_it2.png and /dev/null differ diff --git a/examples/loss_it3.png b/examples/loss_it3.png deleted file mode 100644 index f5b0012..0000000 Binary files a/examples/loss_it3.png and /dev/null differ diff --git a/run.jl b/run.jl index bf33ef8..8a6cd3d 100644 --- a/run.jl +++ b/run.jl @@ -1,7 +1,21 @@ +# Work in a temporary environment using Pkg -Pkg.activate(".") +Pkg.activate(; temp = true) -using Images, CausticsEngineering +# Speed up by avoiding updating the repository when adding packages +Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true -image = Images.load("./examples/cat_posing.jpg") # Check current working directory with pwd() +# Add useful package +Pkg.add([ + "Revise", "Images" +]) + +Pkg.develop(path = @__DIR__) + +using Revise, Images + +using CausticsEngineering + +# Check current working directory with pwd() +image = Images.load("./examples/personal/caricature.jpg") engineer_caustics(image); diff --git a/src/CausticsEngineering.jl b/src/CausticsEngineering.jl index 1bed38b..e132126 100644 --- a/src/CausticsEngineering.jl +++ b/src/CausticsEngineering.jl @@ -1,13 +1,71 @@ module CausticsEngineering using DocStringExtensions +using Random, Images -using Images, Plots +using Plots gr() +using Meshes, FileIO, MeshIO + + +include("../examples/personal/original.jl") + + +include("parameters.jl") include("utilities.jl") +include("math_utilities.jl") +include("mesh.jl") +include("io.jl") +include("plots.jl") + include("create_mesh.jl") -export main, engineer_caustics + +""" +$(SIGNATURES) +""" +function main() + @assert size(ARGS) == (1,) "Intented usage is: julia create_mesh.jl image.png" + + img = Images.load(ARGS[1]) + return engineer_caustics(img) end + +export + # Constants + Meters_Per_Pixel, + Bottom_Offset, + Top_Offset, + + # Types + Vertex3D, + FieldVertex3D, + FaceMesh, + + # math + average, + average_absolute, + laplacian, + ∇, + + # Mesh + create_solid, + create_solid_as_dict, + get_lens_pixels_area, + + # Plotting + plot_as_quiver, + + # I/O + save_obj!, + + # Main procedures + propagate_poisson!, + march_mesh!, + engineer_caustics, + main + + +end # End module diff --git a/src/create_mesh.jl b/src/create_mesh.jl index ad6c5b6..fb5bff7 100644 --- a/src/create_mesh.jl +++ b/src/create_mesh.jl @@ -1,957 +1,356 @@ -# Currently unused. -const grid_definition = 512 - """ $(SIGNATURES) -This func returns a square mesh, centered on zero, with (width * height) nodes -""" -function squareMesh(width::Int, height::Int) - nodeList = Vector{Point3D}(undef, height * width) - nodeArray = Matrix{Point3D}(undef, width, height) - count = 1 - midpoint = width / 2 - for y in 1:height, x in 1:width - newPoint = Point3D(x, y, 0, x, y) - nodeList[count] = newPoint - nodeArray[x, y] = newPoint - count += 1 - end - - triangles = Vector{Triangle}(undef, (width - 1) * (height - 1) * 2) - count = 1 - for y = 1:(height - 1) - for x = 1:(width - 1) - # here x and y establish the column of squares we're in - index_ul = (y - 1) * width + x - index_ur = index_ul + 1 - - index_ll = y * width + x - index_lr = index_ll + 1 - - triangles[count] = Triangle(index_ul, index_ll, index_ur) - count += 1 - triangles[count] = Triangle(index_lr, index_ur, index_ll) - count += 1 - end - end - - return Mesh(nodeList, nodeArray, triangles, width, height) -end - - """ -$(SIGNATURES) +function engineer_caustics(source_image) + imageBW = Float64.(Gray.(source_image)) -Warning: not guaranteed to work in 3D? -""" -function centroid(mesh::Mesh, index::Int) - triangle = mesh.triangles[index] - p1 = mesh.nodes[triangle.pt1] - p2 = mesh.nodes[triangle.pt2] - p3 = mesh.nodes[triangle.pt3] - return centroid(p1, p2, p3) -end + # Set global size parameters + height, width = size(imageBW) + global N_Pixel_Height = height + global N_Pixel_Width = width + global Meters_Per_Pixel = Caustics_Long_Side / N_Pixel_Height + println("Image size: $(N_Pixel_Height) x $(N_Pixel_Width)") + # mesh is the same size as the image with an extra row/column to have coordinates to + # cover each image corner with a triangle. + mesh = FaceMesh(N_Pixel_Height, N_Pixel_Width) -""" -$(SIGNATURES) + # Create a pre-allocated solver + solve_velocity_potential! = solver_potential!(mesh) -Given 3 points and 3 velocities, calculate the `t` required to bring the area of that triangle to zero -""" -function findT( - p1::Point3D, - p2::Point3D, - p3::Point3D, - dp1::Point3D, - dp2::Point3D, - dp3::Point3D, -) - x1 = p2.x - p1.x - y1 = p2.y - p1.y - - x2 = p3.x - p1.x - y2 = p3.y - p1.y - - u1 = dp2.x - dp1.x - v1 = dp2.y - dp1.y - - u2 = dp3.x - dp1.x - v2 = dp3.y - dp1.y - - a = u1 * v2 - u2 * v1 - b = x1 * v1 + y2 * u1 - x2 * v1 - y1 * u2 - c = x1 * y2 - x2 * y1 - - if a != 0 - quotient = b^2 - 4a * c - if quotient >= 0 - d = sqrt(quotient) - return (-b - d) / 2a, (-b + d) / 2a - else - return -123.0, -123.0 - end - else + # imageBW is normalised to have 1 unit per pixel on average. + imageBW /= average(imageBW) - # cool, there just isn't any dependence on t^2, but there is still on t! - return -c / b, -c / b - end -end + marginal_change = nothing + max_update = Inf + counter = 0 + while (abs(max_update) > 1e-6 && counter < 4) + counter += 1 + print( + """ + ================================================================================================================ + STARTING ITERATION $(counter): + starting ϕ field = $(field_summary(mesh.corners.ϕ)) -""" -$(SIGNATURES) - -This function saves the mesh object in obj format. -""" -function saveObj!(mesh::Mesh, filename::String; scale=1.0, scalez=1.0, reverse=false, flipxy=false) - # This function saves the mesh object in stl format - open(filename, "w") do io - for vertex in mesh.nodes - if flipxy - println(io, "v ", vertex.y * scale, " ", vertex.x * scale, " ", vertex.z * scalez) - else - println(io, "v ", vertex.x * scale, " ", vertex.y * scale, " ", vertex.z * scalez) - end - end - - for face in mesh.triangles - if reverse - println(io, "f ", face.pt3, " ", face.pt2, " ", face.pt1) - else - println(io, "f ", face.pt1, " ", face.pt2, " ", face.pt3) - end - end - - println(io, "dims ", mesh.width, " ", mesh.height) - end -end - - -""" -$(SIGNATURES) -""" -function Obj2Mesh(filename) - lines = readlines(filename) - - vertexLines = [l for l in lines if startswith(l, "v")] - nodeList = Vector{Point3D}(undef, size(vertexLines)) - count = 1 - - for line in vertexLines - elements = split(line, " ") - x = parse(Float64, elements[2]) - y = parse(Float64, elements[3]) - z = parse(Float64, elements[4]) * 10 - pt = Point3D(x, y, z, 0, 0) - nodeList[count] = pt - count += 1 - end - - faceLines = [l for l in lines if startswith(l, "f")] - triangles = Vector{Triangle}(undef, size(faceLines)) - for line in faceLines - elements = split(line, " ") - triangle = Triangle( - parse(Int64, elements[2]), - parse(Int64, elements[3]), - parse(Int64, elements[4]), + """, ) - end - dimsLines = [l for l in lines if startswith(l, "dims")] - elements = split(dimsLines[1], " ") + ε, max_update = solve_velocity_potential!(imageBW, "it$(counter)") - return Mesh(nodeList, triangles, parse(Int64, elements[2]), parse(Int64, elements[3])) -end + print( + """ + RESULT AT ITERATION $(counter): + Luminosity error = $(field_summary(ε)) + Vertical move max update = $(max_update) -""" -$(SIGNATURES) -""" -function ∇(f::Matrix{Float64}) - width, height = size(f) + $(field_summary(mesh.corners.vr, "∇u")) + $(field_summary(mesh.corners.vc, "∇v")) + $(field_summary(mesh.corners.r - mesh.corners.rows_numbers, "total mesh changes on row")) + $(field_summary(mesh.corners.c - mesh.corners.cols_numbers, "total mesh changes on col")) - ∇fᵤ = zeros(Float64, width, height) # the right edge will be filled with zeros - ∇fᵥ = zeros(Float64, width, height) # the buttom edge will be filled with zeros + end ϕ field = $(field_summary(mesh.corners.ϕ)) + ================================================================================================================ - for x in 1:width, y in 1:height - ∇fᵤ[x, y] = (x == width ? 0 : f[x + 1, y] - f[x, y]) - ∇fᵥ[x, y] = (y == height ? 0 : f[x, y + 1] - f[x, y]) + """, + ) + plot_as_quiver(mesh, n_steps = 50, scale = height / 10, max_length = height / 20) end - return ∇fᵤ, ∇fᵥ -end + println("\nSTARTING HORIZONTAL ITERATION ---") + max_update = solve_height_potential(mesh, imageBW; f = Focal_Length) + println("\t Horizontal max update = $(max_update)") -""" -$(SIGNATURES) -""" -function getPixelArea(mesh::Mesh) - # A Mesh is a grid of 3D points. The X and Y coordinates are not necessarily aligned or square - # The Z coordinate represents the value. brightness is just proportional to area. - # pixelAreas = Matrix{Float64}(undef, mesh.width - 1, mesh.height - 1) - pixelAreas = zeros(mesh.width - 1, mesh.height - 1) - - for x in 1:mesh.width - 1, y in 1:mesh.height - 1 - upperLeft = mesh.nodeArray[x, y] - upperRight = mesh.nodeArray[x + 1, y] - - lowerLeft = mesh.nodeArray[x, y + 1] - lowerRight = mesh.nodeArray[x + 1, y + 1] - - #= - *------* - | / | - | / | - | / | - | / | - *------* =# - pixelAreas[x, y] = - triangle_area(lowerLeft, upperRight, upperLeft) + - triangle_area(lowerLeft, lowerRight, upperRight) - end + # Move the around a nil average. + mesh.corners.ϕ .-= average(mesh.corners.ϕ) - return pixelAreas + return mesh, imageBW end -""" -$(SIGNATURES) -""" -function relax!(matrix::Matrix{Float64}, D::Matrix{Float64}) - # This function implements successive over relaxation for a matrix and its associated error matrix - # There is a hardcoded assumption of Neumann boundary conditions--that the derivative across the - # boundary must be zero in all cases. See: - # https://math.stackexchange.com/questions/3790299/how-to-iteratively-solve-poissons-equation-with-no-boundary-conditions - # sz = size(matrix) - # width = sz[1] - # height = sz[2] - width, height = size(matrix) - # ω = 2 / (1 + π / width) - ω = 1.99 - # println("OMEGA $(ω)") - max_update = 0 - for y = 1:height - for x = 1:width - val = matrix[x, y] - - if x == 1 && y == 1 - # Top left corner - val_down = matrix[x, y + 1] - val_right = matrix[x + 1, y] - delta = ω / 2 * (val_down + val_right - 2 * val - D[x, y]) - if abs(delta) > max_update - max_update = abs(delta) - end - matrix[x, y] += delta - elseif x == 1 && y == height - # Bottom left corner - val_up = matrix[x, y - 1] - val_right = matrix[x + 1, y] - delta = ω / 2 * (val_up + val_right - 2 * val - D[x, y]) - if abs(delta) > max_update - max_update = abs(delta) - end - matrix[x, y] += delta - elseif x == width && y == 1 - # Top right corner - val_down = matrix[x, y + 1] - val_left = matrix[x - 1, y] - delta = ω / 2 * (val_down + val_left - 2 * val - D[x, y]) - if abs(delta) > max_update - max_update = abs(delta) - end - matrix[x, y] += delta - elseif x == width && y == height - # Bottom right corner - val_up = matrix[x, y - 1] - val_left = matrix[x - 1, y] - delta = ω / 2 * (val_up + val_left - 2 * val - D[x, y]) - if abs(delta) > max_update - max_update = abs(delta) - end - matrix[x, y] += delta - - elseif x == 1 - # Along the left edge, but not the top or buttom corner - val_up = matrix[x, y - 1] - val_down = matrix[x, y + 1] - val_right = matrix[x + 1, y] - delta = ω / 3 * (val_up + val_down + val_right - 3 * val - D[x, y]) - if abs(delta) > max_update - max_update = abs(delta) - end - matrix[x, y] += delta - elseif x == width - # Along the right edge, but not the top or buttom corner - val_up = matrix[x, y - 1] - val_down = matrix[x, y + 1] - val_left = matrix[x - 1, y] - delta = ω / 3 * (val_up + val_down + val_left - 3 * val - D[x, y]) - if abs(delta) > max_update - max_update = abs(delta) - end - matrix[x, y] += delta - elseif y == 1 - # Along the top edge, but not the left or right corner - val_down = matrix[x, y + 1] - val_left = matrix[x - 1, y] - val_right = matrix[x + 1, y] - delta = ω / 3 * (val_down + val_left + val_right - 3 * val - D[x, y]) - if abs(delta) > max_update - max_update = abs(delta) - end - matrix[x, y] += delta - elseif y == height - # Along the bottom edge, but not the left or right corner - val_up = matrix[x, y - 1] - val_left = matrix[x - 1, y] - val_right = matrix[x + 1, y] - delta = ω / 3 * (val_up + val_left + val_right - 3 * val - D[x, y]) - if abs(delta) > max_update - max_update = abs(delta) - end - matrix[x, y] += delta - else - # The normal case, in the middle of the mesh! - val_up = matrix[x, y - 1] - val_down = matrix[x, y + 1] - val_left = matrix[x - 1, y] - val_right = matrix[x + 1, y] - - # The new way - # ∇x₁ = - - - # The old way - delta = ω / 4 * (val_up + val_down + val_left + val_right - 4 * val - D[x, y]) - if abs(delta) > max_update - max_update = abs(delta) - end - matrix[x, y] += delta - end - # node.z = .25 * (node_up.z + node_down.z + node_left.z + node_right.z) # simple averaging - # node.z += ω/4 * (node_up.z + node_down.z + node_left.z + node_right.z - 4 * node.z) - - # matrix[x, y] += ω/4 * (val_up + val_down + val_left + val_right - 4 * val - D[x, y]) - end - end - - max_update - - # for y = 1:height - # for x = 1:width - # val = matrix[x, y] - # end - # end -end - """ $(SIGNATURES) -This function will take a `grid_definition x grid_definition` matrix and returns a `grid_definition x grid_definition` mesh. - -It currently takes the size of the matrix passed as argument. +ϕ is the _velocity potential_. The velocity represents the direction towards which the mesh vertices should move. +Zones of high luminosity should be covered with more triangles and should attract more vertices. """ -function matrix_to_mesh(matrix::Matrix{Float64}) - w, h = size(matrix) - - retval = squareMesh(w, h) - for x in 1:w, y in 1:h - index = (y - 1) * retval.width + x - node = retval.nodes[index] - node.z = matrix[x, y] - end +function solver_potential!(mesh) - return retval -end + # Start with a new clean field + (height, width) = size(mesh) + ε = zeros(Float64, height + 1, width + 1) + lens_pixels_area = zeros(Float64, height, width) + # Preallocate the Poisson solver + propagate_poisson! = poisson_propagator!(height + 1, width + 1) -""" -$(SIGNATURES) -""" -function marchMesh!(mesh::Mesh, ϕ::Matrix{Float64}) - ∇ϕᵤ, ∇ϕᵥ = ∇(ϕ) + return (image, prefix) -> begin - imgWidth, imgHeight = size(ϕ) # should be grid_definitionxgrid_definition + # Reset at each call + fill!(ε, 0.0) + fill!(lens_pixels_area, 0.0) - # For each point in the mesh we need to figure out its velocity - velocities = Matrix{Point3D}(undef, mesh.width, mesh.height) - for x in 1:mesh.width, y in 1:mesh.height - # XY coordinates in the mesh ARE XY coordinates in the image. The mesh just needs an extra row and column - # at the bottom right edge so that the triangles can be closed - - if x == mesh.width - u = 0 - else - u = (y == mesh.height ? ∇ϕᵤ[x, y - 1] : ∇ϕᵤ[x, y]) - end + # Get the area of each individual pixel as stretch/shrunk on the lens. Area = energy. + # _FENCES_SIZED_ + # Illumination only depends on the position of the corners, not their heights. ϕ is not relevant. + lens_pixels_area .= get_lens_pixels_area(mesh) - if y == mesh.height - v = 0 - else - v = (x == mesh.width ? ∇ϕᵥ[x - 1, y] : ∇ϕᵥ[x, y]) - end + # Positive error => the triangle needs to shrink (less light). Enforce nil average error. + ε[1:end-1, 1:end-1] = Float64.(lens_pixels_area - image) + ε .-= average(ε) - velocities[x, y] = Point3D(-u, -v, 0, 0, 0) - end - - min_t = 10000 - triangleCount = 1 - for triangle in mesh.triangles - p1 = mesh.nodes[triangle.pt1] - p2 = mesh.nodes[triangle.pt2] - p3 = mesh.nodes[triangle.pt3] - - v1 = velocities[p1.ix, p1.iy] - v2 = velocities[p2.ix, p2.iy] - v3 = velocities[p3.ix, p3.iy] - - t1, t2 = findT(p1, p2, p3, v1, v2, v3) - - if 0 < t1 < min_t - min_t = t1 - end + println(""" + solve_velocity_potential! before loop: + $(field_summary(lens_pixels_area, "Pixel area")) + $(field_summary(ε, "luminosity error")) + """) - if 0 < t2 < min_t - min_t = t2 - end - triangleCount += 1 - end - - println("Overall min_t:", min_t) - δ = min_t / 2 - - for point in mesh.nodes - v = velocities[point.ix, point.iy] - point.x = v.x * δ + point.x - point.y = v.y * δ + point.y - end - - # saveObj(mesh, "gateau.obj") -end - - - """ -$(SIGNATURES) -""" -function quantifyLoss!(D, suffix, img) - println("Loss:") - println("Minimum: $(minimum(D))") - println("Maximum: $(maximum(D))") - println("SUM: $(sum(D))") - - blue = zeros(size(D)) - blue[D .> 0] = D[D .> 0] - red = zeros(size(D)) - red[D .< 0] = -D[D .< 0] - green = zeros(size(D)) - - println(size(blue)) - println(size(red)) - println(size(green)) - - rgbImg = RGB.(red, green, blue)' - save("./examples/loss_$(suffix).png", map(clamp01nan, rgbImg)) -end - - -""" -$(SIGNATURES) -""" -function oneIteration(meshy, img, suffix) - # Remember meshy is (will be) `grid_definition x grid_definition` just like the image - # `grid_definition x grid_definition`, so LJ is `grid_definition x grid_definition`. - LJ = getPixelArea(meshy) - D = Float64.(LJ - img) - # Shift D to ensure its sum is zero - D .-= sum(D) / (512 * 512) - - # Save the loss image as a png - println(minimum(D)) - println(maximum(D)) - quantifyLoss!(D, suffix, img) - - # These lines are just for plotting the quiver of L, which I needed for the blog post - # ∇Lᵤ, ∇Lᵥ = ∇(D) - # plotVAsQuiver(∇Lᵤ, ∇Lᵥ, stride=10, scale=10, max_length=200) - # println("okay") - # return - - width, height = size(img) - - ϕ = zeros(width, height) - - println("Building Phi") - for i = 1:10000 - max_update = relax!(ϕ, D) - - if isnan(max_update) - println("MAX UPDATE WAS NaN. CANNOT BUILD PHI") - return - end - - if i % 500 == 0 - println(max_update) - end - - if max_update < 0.00001 - println("Convergence reached at step $(i) with max_update of $(max_update)") - break - end - end - - # saveObj!( - # matrix_to_mesh(ϕ * 0.02), - # "./examples/phi_$(suffix).obj", - # reverse=false, - # flipxy=true, - # ) - # plotAsQuiver(ϕ * -1.0, stride=30, scale=1.0, max_length=200, flipxy=true, reversex=false, reversey=false) - # saveObj(matrix_to_mesh(D * 10), "D_$(suffix).obj") - - # Now we need to march the x,y locations in our mesh according to this gradient! - marchMesh!(meshy, ϕ) -# saveObj!(meshy, "./examples/mesh_$(suffix).obj", flipxy=true) -end + # Save the loss image as a png. + save_plot_scalar_field!(ε, "error_$(prefix)") + # Save the generated image as a png. + luminosity_ratio = sum(image) / sum(lens_pixels_area) + save( + "./examples/img_$(prefix).png", + Gray.(clamp.(lens_pixels_area * luminosity_ratio, 0.0, 1.0)), + ) -""" -$(SIGNATURES) -""" -function setHeights!(mesh, heights, heightScale=1.0, heightOffset=10) - width, height = size(heights) - for y = 1:height - for x = 1:width - mesh.nodeArray[x, y].z = heights[x, y] * heightScale + heightOffset - if x == 100 && y == 100 - println("Example heights: $(heights[x, y]) and $(heights[x, y] * heightScale) and $(heights[x, y] * heightScale + heightOffset)") + # Start with a clean, flat potential field. + mesh.corners.ϕ = zeros(Float64, height + 1, width + 1) + count = 0 + max_update = Inf + ∇²ϕ_est = δ = nothing + while (1e-6 < max_update && count < 10_000) + count += 1 + max_update, ∇²ϕ_est, δ = propagate_poisson!(mesh.corners.ϕ, ε) + + # Save the loss image as a png. + if count % 1_000 == 1 + println(""" + Iteration $(count) - max_update = $(round(max_update, sigdigits=4)) + $(field_summary(mesh.corners.ϕ, "ϕ")) + $(field_summary(∇²ϕ_est, "∇²ϕ_est")) + $(field_summary(δ, "δ")) + """) + save_plot_scalar_field!(∇²ϕ_est, "∇²ϕ_est_$(prefix)") + save_plot_scalar_field!(δ, "delta_$(prefix)") + save_plot_scalar_field!(mesh.corners.ϕ, "change_mesh.corners.ϕ_$(prefix)") end end - end - # get the side edge - for y = 1:height - mesh.nodeArray[width + 1, y].z = mesh.nodeArray[width, y].z - end + println( + """ + Finished iteration at count $(count) - max_update = $(round(max_update, sigdigits=4)) + $(field_summary(mesh.corners.ϕ, "ϕ")) + $(field_summary(∇²ϕ_est, "∇²ϕ_est")) + $(field_summary(δ, "δ")) + """, + ) - # get the bottom edge - for x = 1:width + 1 - mesh.nodeArray[x, height + 1].z = mesh.nodeArray[x, height].z + # Now we need to march the mesh row,col corner locations according to this gradient. + δ = march_mesh!(mesh) + return ε, max_update end - - - # # get the pesky corner! -# mesh.nodeArray[width + 1, height + 1].z = mesh.nodeArray[width, height].z end """ $(SIGNATURES) -""" -function solidify(inputMesh, offset=100) - width = inputMesh.width - height = inputMesh.height - totalNodes = width * height * 2 - nodeList = Vector{Point3D}(undef, totalNodes) - nodeArrayTop = Matrix{Point3D}(undef, width, height) - nodeArrayBottom = Matrix{Point3D}(undef, width, height) - - # imagine a 4x4 image. 4 * 2 + 2 * 2 = 12 - numEdgeNodes = width * 2 + (height - 2) * 2 - - numTrianglesTop = (width - 1) * (height - 1) * 2 - numTrianglesBottom = numTrianglesTop - numTrianglesEdges = numEdgeNodes * 2 - - totalTriangles = numTrianglesBottom + numTrianglesTop + numTrianglesEdges - - println("Specs: $(width) $(height) $(totalNodes) $(numEdgeNodes) $(numTrianglesBottom) $(totalTriangles)") - - # Build the bottom surface - count = 1 - for y = 1:height - for x = 1:width - newPoint = Point3D(x, y, -offset, x, y) - nodeList[count] = newPoint - nodeArrayBottom[x, y] = newPoint - count += 1 - end - end - # Copy in the top surface - for y = 1:height - for x = 1:width - node = inputMesh.nodeArray[x, y] - copiedPoint = Point3D(node.x, node.y, node.z, node.ix, node.iy) - if node.ix != x - println("OH NO POINTS NOT MATCHED $(x) vs $(node.ix)") - end - if node.iy != y - println("OH NO POINTS NOT MATCHED $(y) vs $(node.iy)") - end +ϕ is the height map. - nodeList[count] = copiedPoint - nodeArrayTop[x, y] = copiedPoint - count += 1 - end - end - - println("We now have $(count - 1) valid nodes") - - triangles = Vector{Triangle}(undef, totalTriangles) - # Build the triangles for the bottom surface - count = 1 - for y = 1:(height - 1) - for x = 1:(width - 1) - # here x and y establish the column of squares we're in - index_ul = (y - 1) * width + x - index_ur = index_ul + 1 - - index_ll = y * width + x - index_lr = index_ll + 1 +`march_mesh!` flexes the mesh +""" +function march_mesh!(mesh::FaceMesh) - triangles[count] = Triangle(index_ul, index_ll, index_ur) - count += 1 - triangles[count] = Triangle(index_lr, index_ur, index_ll) - count += 1 - end - end + # Calculate the gradient of the velocity potential and reverse its direction. + mesh.corners.vr, mesh.corners.vc = ∇(mesh.corners.ϕ) + mesh.corners.vr .*= -1.0 + mesh.corners.vc .*= -1.0 - println("We've filled up $(count - 1) triangles") - if count != numTrianglesBottom + 1 - println("Hmm aren't count and triangles bottom equal? $(count) vs $(numTrianglesBottom + 1)") - end + # Clean up the mesh borders + reset_border_values!(mesh.corners) - # Build the triangles for the top surface - for y = 1:(height - 1) - for x = 1:(width - 1) - # here x and y establish the column of squares we're in - index_ul = (y - 1) * width + x + totalNodes / 2 - index_ur = index_ul + 1 + # For each point in the mesh we need to figure out its velocity + # However all the nodes located at a border will never move + # I.e. velocity (Vx, Vy) = (0, 0) and the square of acrylate will remain of the same size. + # Get the time, at that velocity, for the area of the triangle to be nil. + # We are only interested by positive values to only move in the direction of the (opposite) gradient + height, width = size(mesh) + list_triangles = vcat( + [triangle3D(mesh, row, col, :top) for row ∈ 1:height, col ∈ 1:width], + [triangle3D(mesh, row, col, :bottom) for row ∈ 1:height, col ∈ 1:width], + ) - index_ll = y * width + x + totalNodes / 2 - index_lr = index_ll + 1 + list_maximum_t = [t for t ∈ find_maximum_t.(list_triangles) if !isnothing(t) && t > 0.0] - triangles[count] = Triangle(index_ul, index_ur, index_ll) - count += 1 - triangles[count] = Triangle(index_lr, index_ll, index_ur) - count += 1 - end - end + δ = minimum(list_maximum_t) / 2.0 + mesh.corners.r .-= δ * mesh.corners.vr + mesh.corners.c .-= δ * mesh.corners.vc - println("We've filled up $(count - 1) triangles") - - # Build the triangles to close the mesh - x = 1 - for y = 1:(height - 1) - ll = (y - 1) * width + x - ul = ll + totalNodes / 2 - lr = y * width + x - ur = lr + totalNodes / 2 - triangles[count] = Triangle(ll, ul, ur) - count += 1 - triangles[count] = Triangle(ur, lr, ll) - count += 1 - end + println(""" - x = width - for y = 1:(height - 1) - ll = (y - 1) * width + x - ul = ll + totalNodes / 2 - lr = y * width + x - ur = lr + totalNodes / 2 - triangles[count] = Triangle(ll, ur, ul) - count += 1 - triangles[count] = Triangle(ur, ll, lr) - count += 1 - end + March mesh with correction_ratio δ = $(δ) + $(field_summary(mesh.corners.vr, "∇u")) + $(field_summary(mesh.corners.vc, "∇v")) - y = 1 - for x = 2:width - ll = (y - 1) * width + x - ul = ll + totalNodes / 2 - lr = (y - 1) * width + (x - 1) - ur = lr + totalNodes / 2 - triangles[count] = Triangle(ll, ul, ur) - count += 1 - triangles[count] = Triangle(ur, lr, ll) - count += 1 - end - - y = height - for x = 2:width - ll = (y - 1) * width + x - ul = ll + totalNodes / 2 - lr = (y - 1) * width + (x - 1) - ur = lr + totalNodes / 2 - triangles[count] = Triangle(ll, ur, ul) - count += 1 - triangles[count] = Triangle(ur, ll, lr) - count += 1 - end + """) -Mesh(nodeList, nodeArrayBottom, triangles, width, height) + return δ end + """ $(SIGNATURES) """ -function findSurface(mesh, image, f, imgWidth) - width, height = size(image) - - # imgWidth = .1 # m - # f = 1.0 # m - H = f - metersPerPixel = imgWidth / width - println(metersPerPixel) +function solve_height_potential(mesh::FaceMesh, image; f = Focal_Length) + # height indexes rows, width indexes columns + height, width = size(mesh) - # η = 1.49 - n₂ = 1 - n₁ = 1.49 - Nx = zeros(width + 1, height + 1) - Ny = zeros(width + 1, height + 1) + # Preallocate the Poisson solver + propagate_poisson! = poisson_propagator!(height + 1, width + 1) - for j = 1:height - for i = 1:width - node = mesh.nodeArray[i, j] - dx = (node.ix - node.x) * metersPerPixel - dy = (node.iy - node.y) * metersPerPixel - little_h = node.z * metersPerPixel - H_minus_h = H - little_h - dz = H_minus_h + # Coordinates on the caustics = simple rectangular values in corners. + # Coordinates on the lens face comes from corner field. + # d_row/d_col are _POSTS-SIZED_ + d_row = mesh.corners.rows_numbers - mesh.corners.r + d_col = mesh.corners.cols_numbers - mesh.corners.c + # The focal length is in meters. H is in corners. + H = f / Meters_Per_Pixel - # k = η * sqrt(dx * dx + dy * dy + H_minus_h * H_minus_h) - H_minus_h - # Nx[i, j] = 1/k * dx - # Ny[i, j] = 1/k * dy - Ny[i, j] = tan(atan(dy / dz) / (n₁ - 1)) - Nx[i, j] = tan(atan(dx / dz) / (n₁ - 1)) + # true_H is the real distance beteen the surface of the lens and the projection + # H is the initial distance from the surface of the carved face to the projection screen. + # H is constant and does not account for the change in heights due to carving. + # Note: Higher location means closer to the caustics => negative sign + # true_H is _POSTS_SIZED_ + true_H = zeros(Float64, height + 1, width + 1) + true_H = -mesh.corners.ϕ .+ H + # Normals. + # N_row/N_col are _POSTS_SIZED_ + N_row = zeros(Float64, size(true_H)) + N_col = zeros(Float64, size(true_H)) + N_row = tan.(atan.(d_row ./ true_H) / (n₁ - n₂)) + N_col = tan.(atan.(d_col ./ true_H) / (n₁ - n₂)) - end - end + # @. N_row[1:end, 1:end] = + # tan(atan(d_row[1:end, 1:end] / true_H[1:end, 1:end]) / (n₁ - n₂)) + # @. N_col[1:end, 1:end] = + # tan(atan(d_col[1:end, 1:end] / true_H[1:end, 1:end]) / (n₁ - n₂)) - divergence = zeros(width, height) # We need to find the divergence of the Vector field described by Nx and Ny - - for j = 1:height - for i = 1:width - δx = (Nx[i + 1, j] - Nx[i, j]) - δy = (Ny[i, j + 1] - Ny[i, j]) - divergence[i, j] = δx + δy - end - end - println("Have all the divergences") - println("Divergence sum: $(sum(divergence))") - divergence .-= sum(divergence) / (width * height) - - h = zeros(width, height) - max_update = 0 - for i = 1:10000 - max_update = relax!(h, divergence) - - if i % 500 == 0 - println(max_update) - end - if max_update < 0.00001 - println("Convergence reached at step $(i) with max_update of $(max_update)") + # div_row/div_col are _FENCES_SIZED_ + div_row = zeros(Float64, size(true_H)) + div_col = zeros(Float64, size(true_H)) + div_row[1:end-1, 1:end-1] = N_row[2:end, 1:end-1] .- N_row[1:end-1, 1:end-1] + div_col[1:end-1, 1:end-1] = N_col[1:end-1, 2:end] .- N_col[1:end-1, 1:end-1] + + # Normalised divergence + # divergence_direction is _FENCES_SIZED_ + divergence_direction = div_row + div_col + # divergence_direction .-= average(divergence_direction) + + max_update = Inf + old_update = Inf + for i ∈ 1:1e6 + old_update = max_update + max_update, _, _ = propagate_poisson!(mesh.corners.ϕ, -divergence_direction) + + i % 1_000 == 1 && + println("Convergence horiz. Δ at counter = $(i) max update = $(max_update)") + + if abs(max_update) <= 1e-6 + println( + "Convergence horiz. Δ stopped at counter = $(i) with max_update of $(max_update)", + ) break end end - # println("HEIGHT STATS") - # h .-= sum(h) / (width * height) - # println(minimum(h)) - # println(maximum(h)) - # println(sum(h)) - # println("-----------") - - # saveObj!(matrix_to_mesh(h * 10), "examples/heightmap.obj") - h, metersPerPixel + return max_update end - -""" -$(SIGNATURES) -""" -function testSquareMesh!() - mesh = squareMesh(100, 50) - - println(mesh.nodeArray[1, 1]) - println(mesh.nodes[1]) - mesh.nodeArray[1, 1].x = 8 - println(mesh.nodeArray[1, 1]) - println(mesh.nodes[1]) - mesh.nodes[1].y += 12 -println(mesh.nodeArray[1, 1]) - println(mesh.nodes[1]) - -end """ $(SIGNATURES) -""" -function testSolidify!() - println("Testing solidification") - width = 100 - height = 100 - origMesh = squareMesh(width, height) - - for y in 1:height, x in 1:width - x2 = (x - width / 2) / width - y2 = (y - height / 2) / height - value = x2 * x2 + y2 * y2 - origMesh.nodeArray[x, y].z = 15 - value * 25 - end - - saveObj!(origMesh, "./examples/testSolidify.obj") - solidMesh = solidify(origMesh, 0) - saveObj!(solidMesh, "./examples/testSolidify2.obj") -end +This function implements successive over relaxation for a matrix and its associated error matrix +There is a hardcoded assumption of Neumann boundary conditions -- that the derivative across the +boundary must be zero in all cases. See: +https://math.stackexchange.com/questions/3790299/how-to-iteratively-solve-poissons-equation-with-no-boundary-conditions -""" -$(SIGNATURES) -""" -function plotAsQuiver( - g; - stride=4, - scale=300, - max_length=2, - flipxy=false, - reversey=false, - reversex=false, -) - h, w = size(g) - xs = Float64[] - ys = Float64[] - us = Float64[] - vs = Float64[] - - for x in 1:stride:w, y in 1:stride:h - reversex ? push!(xs, x) : push!(xs, -x) - reversey ? push!(ys, -y) : push!(ys, y) - - p1 = g[y, x] - u = (g[y, x + 1] - g[y, x]) * scale - v = (g[y + 1, x] - g[y, x]) * scale - - u = -u - - reversey && (v = -v) - reversex && (u = -u) - - # println(u, v) - u >= 0 ? push!(us, min(u, max_length)) : push!(us, max(u, -max_length)) - v >= 0 ? push!(vs, min(v, max_length)) : push!(vs, max(v, -max_length)) - end - - q = - flipxy ? quiver(ys, xs, quiver=(vs, us), aspect_ratio=:equal) : - quiver(xs, ys, quiver=(us, vs), aspect_ratio=:equal) +- ϕ is the potential to be solved subject to the Poisson equation and is the same size as the corners (_POSTS_SIZED_) +- target is ∇ϕ and is of the size in pixels (_FENCES_SIZED_) - display(q) - readline() -end +In order of how to think about the flow of the program for the velocity potential, + - if δ > 0, the Laplacian is too high at that point (∇²ϕ is too low). + - To decrease the Laplacian ∇²ϕ, note that ∇² is a divergence (of a gradient) that represents net flow going outside. + - Decrease the Laplacian + - => decrease a divergence + - => decrease flow towards outside + - The flow is a gradient of the potential ϕ. + - => decrease flow towards outside + - => either decrease the values of the potential around that point, or increase the value of the potential at that point. + - => Increase the value of ϕ at that point means a positive change when δ > 0. + - δ > 0 => correction of the SAME sign. """ -$(SIGNATURES) -""" -function plotVAsQuiver(vx, vy; stride=4, scale=300, max_length=2) - h, w = size(vx) +function poisson_propagator!(height, width) - xs = Float64[] - ys = Float64[] - us = Float64[] - vs = Float64[] + ∇²ϕ_est = zeros(Float64, height, width) + δ = zeros(Float64, height, width) - for x in 1:stride:w, y in 1:stride:h - push!(xs, x) - push!(ys, h - y) + # Return a closure over preallocaed arrays + return (ϕ::Matrix{Float64}, ε::Matrix{Float64}) -> begin + ω = 1.99 - u = max(vx[x, y], 0.001) - v = max(vy[x, y], 0.001) + # Reset at each call + fill!(∇²ϕ_est, 0.0) + fill!(δ, 0.0) - push!(us, u) - push!(vs, v) - # println(u, ": ", v) - end + for row = 1:height, col = 1:width + val_up = row == 1 ? 0.0 : ϕ[row-1, col] + val_down = row == height ? 0.0 : ϕ[row+1, col] + val_left = col == 1 ? 0.0 : ϕ[row, col-1] + val_right = col == width ? 0.0 : ϕ[row, col+1] - # readline() - q = quiver(xs, ys, quiver=(us, vs), aspect_ratio=:equal) - display(q) - readline() -end + ∇²ϕ_est[row, col] = + val_up + val_down + val_left + val_right - 4 * ϕ[row, col] + δ[row, col] = ω / 4.0 * (∇²ϕ_est[row, col] - ε[row, col]) + ϕ[row, col] += δ[row, col] + end + return maximum(abs.(δ)), ∇²ϕ_est, δ + # Problem with the in-place updates if using the followng. + # # Laplacian + # # ∇²ϕ_est represent the estimate of the Laplacian of ϕ (second-order value) + # ∇²ϕ_est = laplacian(ϕ) -""" -$(SIGNATURES) -""" -function engineer_caustics(img) - img = Gray.(img) - img2 = permutedims(img) * 1.0 - width, height = size(img2) - - # meshy is the same size as the image - meshy = squareMesh(width + 1, height + 1) - - # We need to boost the brightness of the image so that its sum and the sum of the area are equal - mesh_sum = width * height - image_sum = sum(img2) - boost_ratio = mesh_sum / image_sum - - # img3 is `grid_definition x grid_definition` and is normalised to the same (sort of) _energy_ as the - # original image. - img3 = img2 .* boost_ratio - - oneIteration(meshy, img3, "it1") - oneIteration(meshy, img3, "it2") - oneIteration(meshy, img3, "it3") - oneIteration(meshy, img3, "it4") - # oneIteration(meshy, img3, "it5") - # oneIteration(meshy, img3, "it6") - - artifactSize = 0.1 # meters - focalLength = 0.2 # meters - h, metersPerPixel = findSurface(meshy, img3, focalLength, artifactSize) - - setHeights!(meshy, h, 1.0, 10) - - solidMesh = solidify(meshy) - saveObj!( - solidMesh, - "./examples/original_image.obj", - scale=1 / 512 * artifactSize, - scalez=1 / 512.0 * artifactSize, - ) - return meshy, img3 -end + # # In the case of the velocity potential, δ is the difference between the divergence of the luminosity error + # # and the current one for the current field ϕ (divergence of its gradient). + # # High luminosity error ∇²ϕ means that the triangles are too bright = too large. + # # The Laplacian Lϕ has to converge towards the ∇²ϕ. Difference to calculate speed of the descent. + # # Correction ratio to avoid setting the triangle to nothing + # δ = ω / 4.0 * (∇²ϕ_est - ε) -""" -$(SIGNATURES) -""" -function main() - @assert size(ARGS) == (1,) "Intented usage is: julia create_mesh.jl image.png" + # # In the case of the velocity potential, if δ > 0, the Laplacian is too high at that point (∇²ϕ is too low). + # ϕ[1:end-1, 1:end-1] .+= δ - img = load(ARGS[1]) - return engineer_caustics(img) + # return ∇²ϕ_est, δ, maximum(abs.(δ)) + end end - - -""" -$(SIGNATURES) -""" -main() diff --git a/src/debug_script.jl b/src/debug_script.jl new file mode 100644 index 0000000..cf32b11 --- /dev/null +++ b/src/debug_script.jl @@ -0,0 +1,7 @@ +using Revise, Debugger + +using Images +image = Images.load("./examples/cat_posing.jpg"); # Check current working directory with pwd() + +using CausticsEngineering +mesh, image_caustics = engineer_caustics(image); diff --git a/src/io.jl b/src/io.jl new file mode 100644 index 0000000..cc49642 --- /dev/null +++ b/src/io.jl @@ -0,0 +1,112 @@ + +# Utility function to name a vertex +vertex_name(face::String, row::Int64, col::Int64)::String = "Vertex_$(name)_$(row)_$(col)" + +# Utility functions to create the names +triangle_name(face::String, row::Int64, col::Int64, side::Union{Symbol,Symbol})::String = + "Tri_$(face)_$(row)_$(col)_$(side)" + + + +""" +$(SIGNATURES) + +""" +function convert(::Meshes.SimpleMesh, mesh::FaceMesh) + height, width = size(mesh) + + points = Meshes.Point3.([mesh.corners[ci] for ci ∈ CartesianIndices(mesh.corners)]) + + top_connections = + [connect(mesh.toptriangles[ci]) for ci ∈ CartesianIndices(mesh.corners)] + bot_connections = + [connect(mesh.bottriangles[ci]) for ci ∈ CartesianIndices(mesh.corners)] + + return SimpleMesh(points, vcat(top_connections, bot_connections)) +end + + +""" +$(SIGNATURES) + +""" +function save_obj!( + triangle_index, + vertex_index::AbstractMatrix{Float64}, + filename::String; + scale = 1.0, + scaleh = 1.0, + reverse = false, + flipxy = false, +) + + open(filename, "w") do io + println(io, "solid engineered_caustics") + + xs = vertex_index[:, 1] + ys = vertex_index[:, 2] + zs = vertex_index[:, 3] + for i ∈ 1:length(xs) + if flipxy + println(io, "v $(ys[i] * scale) $(xs[i][1] * scale) $(zs[i] * scaleh)") + else + println(io, "v $(xs[i] * scale) $(ys[i][1] * scale) $(zs[i] * scaleh)") + end + end + + t1 = triangle_index[1] + t2 = triangle_index[2] + t3 = triangle_index[3] + for i ∈ 1:length(t1) + println(io, "f $(t1[i]) $(t2[i]) $(t3[i])") + end + + println(io, "endsolid engineered_caustics") + end +end + + + +""" +$(SIGNATURES) + +TO REFACTOR. +""" +function stl_2_mesh(filename) + lines = readlines(filename) + + vertexLines = [l for l in lines if startswith(l, "v")] + nodeList = Vector{Vertex3D}(undef, size(vertexLines)) + + count = 0 + for line in vertexLines + count += 1 + elements = split(line, " ") + r = parse(Float64, elements[2]) + c = parse(Float64, elements[3]) + h = parse(Float64, elements[4]) * 10 + pt = Vertex3D(r, c, h, 0, 0) + nodeList[count] = pt + end + + faceLines = [l for l in lines if startswith(l, "f")] + triangles = Vector{Triangle}(undef, size(faceLines)) + for line in faceLines + elements = split(line, " ") + triangles[line] = Triangle( + parse(Int64, elements[2]), + parse(Int64, elements[3]), + parse(Int64, elements[4]), + ) + end + + dimsLines = [l for l in lines if startswith(l, "dims")] + elements = split(dimsLines[1], " ") + + return RectangleMesh( + nodeList, + triangles, + parse(Int64, elements[2]), + parse(Int64, elements[3]), + ) +end diff --git a/src/math_utilities.jl b/src/math_utilities.jl new file mode 100644 index 0000000..863fbda --- /dev/null +++ b/src/math_utilities.jl @@ -0,0 +1,77 @@ + +""" +$(SIGNATURES) +""" +average(m::AbstractMatrix) = (sum(m) / length(m)) + + +""" +$(SIGNATURES) +""" +average_absolute(m::AbstractMatrix) = (sum(abs.(m)) / length(m)) + + +""" +$(SIGNATURES) +""" +test_average_absolute(m::AbstractMatrix) = (average_absolute(m)) < 1e6 + + +""" +$(SIGNATURES) +""" +function smallest_positive(x1::Float64, x2::Float64) + (x1 > 0.0 && x2 < 0.0) && return x1 + (x1 < 0.0 && x2 > 0.0) && return x2 + (x1 > 0.0 && x2 > 0.0) && return min(x1, x2) + return nothing +end + + +""" +$(SIGNATURES) + +Approximates the gradient of a scalar field. +""" +function ∇(ϕ::AbstractMatrix{Float64}) + height, width = size(ϕ) + + # x, y only represents the first and the second variables. No reference to graphical representations. + ∇ϕx = zeros(Float64, height, width) # divergence on the right edge will be filled with zeros + ∇ϕy = zeros(Float64, height, width) # divergence on bottom edge will be filled with zeros + + ∇ϕx[begin:end-1, :] = ϕ[begin+1:end, :] - ϕ[begin:end-1, :] + ∇ϕy[:, begin:end-1] = ϕ[:, begin+1:end] - ϕ[:, begin:end-1] + + fill_borders!(∇ϕx, 0.0) + fill_borders!(∇ϕy, 0.0) + + return ∇ϕx, ∇ϕy +end + + +""" +$(SIGNATURES) + +Calculate the second-order Laplace operator by convolution of a kernel. +""" +function laplacian(ϕ::AbstractMatrix{Float64}) + height, width = size(ϕ) + + # Embed matrix within a larger matrix for better vectorization and avoid duplicated code + # Padded matrix adds 1 row/col after the size of ϕ. + # ϕ is inserted in padded matrix within 1:1+height+1 x 1:1+width+1. + # The rest of the padded matrix (its borders) are set at 0.0. + padded_ϕ = zeros(Float64, 1 + height + 1, 1 + width + 1) + padded_ϕ[begin+1:end-1, begin+1:end-1] .= ϕ[1:end, 1:end] + + # Convolution = up + down + left + right + ∇²ϕ = + padded_ϕ[begin+1-1:end-1-1, begin+1+0:end-1+0] + + padded_ϕ[begin+1+1:end-1+1, begin+1+0:end-1+0] + + padded_ϕ[begin+1+0:end-1+0, begin+1-1:end-1-1] + + padded_ϕ[begin+1+0:end-1+0, begin+1+1:end-1+1] - + 4.0 * padded_ϕ[begin+1:end-1, begin+1:end-1] + + return ∇²ϕ[1:end-1, 1:end-1] +end diff --git a/src/mesh.jl b/src/mesh.jl new file mode 100644 index 0000000..12e8e98 --- /dev/null +++ b/src/mesh.jl @@ -0,0 +1,273 @@ +""" +$(SIGNATURES) + +Given 3 points and their velocities, calculate the time `t` required to bring the area of that triangle to zero +""" +function find_maximum_t(p1::Vertex3D, p2::Vertex3D, p3::Vertex3D) + # Three points A, B and C, with coordinates (x, y) + # The area of a triangle is 1/2 * [ Ax (By - Cy) + Bx (Cy - Ay) + Cx (Ay - By)] + # where each point of the triangle is where it will be after time t + # i.e. a point goes from P to P+tV where V is the velocity of that point. + # Notation is here B -> B + t_vB + + # To make the calculation simpler, everything is translated so that A is at the + # origin of the plane and its velocity is nil. + Br = p2.r - p1.r + Bc = p2.c - p1.c + Cr = p3.r - p1.r + Cc = p3.c - p1.c + + t_vBr = p2.vr - p1.vr + t_vBc = p2.vc - p1.vc + t_vCr = p3.vr - p1.vr + t_vCc = p3.vc - p1.vc + + # After this, given that Ar = Ac = t_vAr = t_vAc = 0, the area is nil iff + # (Br + t_vBr) (Cc + t_vCc ) - (Cr + t_vCr) (Bc + t_vBc) = 0. + # After expansion and reshuffling to have a quadratic equation where t + # is the variable, the coefficients of that equation are: + a = t_vCc * t_vBr - t_vBc * t_vCr + b = -Bc * t_vCr - Cr * t_vBc + Br * t_vCc + Cc * t_vBr + c = Br * Cc - Cr * Bc + + # if a = 0, this is just a linear equation. + if a == 0 && b != 0 + return smallest_positive(-c / b, c / b) + else + discriminant = b^2 - 4a * c + + # If there is a solution + if discriminant >= 0 + d = sqrt(discriminant) + return smallest_positive((-b - d) / 2a, (-b + d) / 2a) + end + end + # There can be no solution if, after translation, B abd C move in parallel direction. + # C will never end up on the line AB. + # Very unlikely with Float64. + # If nothing found, no move + return nothing +end + +find_maximum_t(p::Tuple{Vertex3D,Vertex3D,Vertex3D}) = find_maximum_t(p[1], p[2], p[3]) + + +""" +$(SIGNATURES) + +Create a mash of the solid block to be carved. +Imagining a mid-plane where the heightmap is carved, start with a block whose lower plane is +bottom_distance away and the highr plane is top_distance above. +The distances are measured in meters and converted to pixels when necessary. + +The final mesh is a vector of vertices each with an ID number and a position. + +""" +function create_solid_as_dict( + mesh::FaceMesh; + bottom_distance = Bottom_Offset, + top_distance = Top_Offset, +) + + height, width = size(mesh) + bottom_offset = bottom_distance / Meters_Per_Pixel + top_offset = top_distance / Meters_Per_Pixel + + # To do this, we use several dictionaries before saving + # Dictionary 1: Triangle to 3 strings each naming a triangle vertex. + # Dictionary 2: Named vertex to its x, y, z coordinates. + # Dictionary 3: Named vertex to a unique numerical index. + triangle_dict = Dict{String,Vector{String}}() + vertex_dict = Dict{String,Tuple{Float64,Float64,Float64}}() + + mesh_height, mesh_width = size(mesh.corners.ϕ) + for row ∈ 1:mesh_height-1, col ∈ 1:mesh_width-1 + + # Top face, add vertices, top and bottom triangles + for side ∈ [:top, :bottom] + side_name = side == :top ? "Top" : "Bottom" + t = side == :top ? mesh.toptriangles[row, col] : mesh.bottriangles[row, col] + + # Add the vertices 1 by 1 + vertices_names = ["", "", ""] + for i ∈ 1:3 + v_row, v_col = t[i] + r = mesh.corners.r[v_row, v_col] + c = mesh.corners.c[v_row, v_col] + ϕ = mesh.corners.ϕ[v_row, v_col] + name = vertex_name("Top", row, col) + vertices_names[i] = name + vertex_dict[name] = (r, c, ϕ) + end + + # Add the triangle + triangle_dict[triangle_name(side_name, row, col, side)] = vertices_names + end + + # Bottom face + v_coord_top = [(row, col), (row + 1, col), (row, col + 1)] + vertices_names_top = ["", "", ""] + + v_coord_bot = [(row, col + 1), (row + 1, col), (row + 1, col + 1)] + vertices_names_bot = ["", "", ""] + + # Add the vertices 1 by one + for i ∈ 1:3 + (v_row, v_col) = v_coord_top[i] + name = vertex_name("Bottom", v_row, v_col) + vertices_names_top[i] = name + vertex_dict[name] = + (Float64(v_row), Float64(v_col), Float64(-bottom_offset * 0)) + + (v_row, v_col) = v_coord_bot[i] + name = vertex_name("Bottom", v_row, v_col) + vertices_names_bot[i] = name + vertex_dict[name] = (Float64(v_row), Float64(v_col), Float64(-bottom_offset)) + end + + # Add the triangles + triangle_dict[triangle_name("Bottom", row, col, :top)] = vertices_names_top + triangle_dict[triangle_name("Bottom", row, col, :bottom)] = vertices_names_bot + end + + + ### + ### Side meshes + ### + # Build triangles to create side meshes and close the mesh + # The mesh is made of a single couple top/bottom triangles joining the bottom face + # and the top face. + # All the vertices already exist + + # (Looking from above), west/east sides. + for (col, face_name) ∈ [(1, "West"), (width, "East")] + for row = 1:height-1 + + v_coord_top = [(:top, row, col), (:bottom, row, col), (:top, row + 1, col)] + vertices_names_top = ["", "", ""] + for i ∈ 1:3 + vertices_names_top[i] = vertex_name( + (v_coord_top[i][1] == :top ? "Top" : "Bottom"), + v_coord_top[i][2], + v_coord_top[i][3], + ) + end + + v_coord_bot = + [(:top, row + 1, col), (:bottom, row + 1, col), (:bottom, row, col)] + vertices_names_bot = ["", "", ""] + for i ∈ 1:3 + vertices_names_bot[i] = vertex_name( + (v_coord_bot[i][1] == :top ? "Top" : "Bottom"), + v_coord_bot[i][2], + v_coord_bot[i][3], + ) + end + + # Add the triangles + triangle_dict[triangle_name(face_name, row, col, :top)] = vertices_names_top + triangle_dict[triangle_name(face_name, row, col, :bottom)] = vertices_names_bot + end + end + + # (Looking from above), north/south sides. + for (row, face_name) ∈ [(1, "North"), (height, "South")] + for col ∈ 1:width-1 + + v_coord_top = [(:top, row, col), (:bottom, row, col), (:top, row, col + 1)] + vertices_names_top = ["", "", ""] + for i ∈ 1:3 + vertices_names_top[i] = vertex_name( + (v_coord_top[i][1] == :top ? "Top" : "Bottom"), + v_coord_top[i][2], + v_coord_top[i][3], + ) + end + + v_coord_bot = [(:top, row, col + 1), (:bottom, row, col + 1), (:top, row, col)] + vertices_names_bot = ["", "", ""] + for i ∈ 1:3 + vertices_names_bot[i] = vertex_name( + (v_coord_bot[i][1] == :top ? "Top" : "Bottom"), + v_coord_bot[i][2], + v_coord_bot[i][3], + ) + end + + + # Add the triangles + triangle_dict[triangle_name(face_name, row, col, :top)] = vertices_names_top + triangle_dict[triangle_name(face_name, row, col, :bottom)] = vertices_names_bot + end + end + + return triangle_dict, vertex_dict +end + + +""" +$(SIGNATURES) + +""" +function create_solid( + mesh::FaceMesh; + bottom_distance = Bottom_Offset, + top_distance = Top_Offset, +) + + triangle_dict, vertex_dict = create_solid_as_dict( + mesh; + bottom_distance = Bottom_Offset, + top_distance = Top_Offset, + ) + + _, _, triangle_index, vertex_index = create_solid( + triangle_dict, + vertex_dict; + bottom_distance = bottom_distance, + top_distance = top_distance, + ) + + return triangle_dict, vertex_dict, triangle_index, vertex_index +end + + + +""" +$(SIGNATURES) + +""" +function create_solid( + triangle_dict::Dict{String,Vector{String}}, + vertex_dict::Dict{String,Tuple{Float64,Float64,Float64}}; + bottom_distance = Bottom_Offset, + top_distance = Top_Offset, +) + + index_dict = Dict{String,Int64}() + index_of_vertex = 0 + for name ∈ keys(vertex_dict) + index_of_vertex += 1 + index_dict[name] = index_of_vertex + end + + vertex_index = Vector{Tuple{Float64,Float64,Float64}}(undef, length(vertex_dict)) + # For each vertex using its name and coordinates + for (name, v) ∈ vertex_dict + # Determine the index of the vertex + i = index_dict[name] + x = v[1] + y = v[2] + z = v[3] + vertex_index[i] = (x, y, z) + end + + # The connections are listed as a tuple of 3 arrays. + # First array contains the first point of a triangle, and so on + triangle_index = [ + (index_dict[v[1]], index_dict[v[2]], index_dict[v[3]]) for v ∈ values(triangle_dict) + ] + + return triangle_dict, vertex_dict, triangle_index, vertex_index + +end diff --git a/src/parameters.jl b/src/parameters.jl new file mode 100644 index 0000000..da60a3b --- /dev/null +++ b/src/parameters.jl @@ -0,0 +1,42 @@ +# Constants + +# Physical parameters + +## Materials +const n₁ = 1.49 # acrylate +const n₂ = 1.00 # air + +## The sides of the block to be carved +const Block_Side = 0.15 # meters +const Block_Height = Block_Side +const Block_Width = Block_Side + +## Depth parameters. If we imagine a reference place, the top face (facing the caustics) is at +## Top_Offset above it. The flat face facing the light source is Bottom below it. +const Top_Offset = 0.01 # 1 centimeters +const Bottom_Offset = 0.02 # 2 centimeters + +## Optics +const Focal_Length = 1.0 # meters + +const Caustics_Long_Side = 0.20 # m +global Caustics_Height = Caustics_Long_Side +global Caustics_Width = Caustics_Long_Side + +# Caustics picture +global N_Pixel_Side = 512 +global N_Pixel_Height = N_Pixel_Side +global N_Pixel_Width = N_Pixel_Side + +global Meters_Per_Pixel = Caustics_Height / N_Pixel_Height + + +# calculation +const ω = 1.99 +# const ω = 2 / (1 + π / sqrt(N_Pixel_Height * N_Pixel_Width)) +const N_Iterations_Convergence = 10_000 # CHECK: What is a reasonable number of iterations? + +# # Global allocations to avoid excessive re-allocations. USEFUL ??? +# global container = zeros(Float64, 1_024, 1_024) +# global divergence_intensity = zeros(Float64, 1_024, 1_024) +# global target_map = zeros(Float64, 1_024, 1_024) diff --git a/src/plots.jl b/src/plots.jl new file mode 100644 index 0000000..b3bfa82 --- /dev/null +++ b/src/plots.jl @@ -0,0 +1,201 @@ + +""" +$(SIGNATURES) +""" +function plot_as_quiver( + mesh::FaceMesh; + n_steps = 20, + scale = 30, + max_length = N_Pixel_Width / n_steps, + fliprc = false, + reverser = false, + reversec = false, +) + + height, width = size(mesh) + stride = Int64(ceil(max(height, width) / n_steps)) + + row_s = Float64[] + col_s = Float64[] + v_rows = Float64[] + v_cols = Float64[] + + # The gradient vectors reflect change of luminosity:outward flow = bigger cells = more luminosity + # Plots are cosmetically look better when showing the light focusing to a specific location. + # Those gradients are the opposite direction. + ϕ = -mesh.corners.ϕ + + v_length = sqrt.(mesh.corners.vr .^ 2 + mesh.corners.vr .^ 2) + v_max = maximum(v_length) + + mat_vr = (reverser ? -1 : 1) .* mesh.corners.vr * scale / v_max + mat_vc = (reversec ? -1 : 1) .* mesh.corners.vc * scale / v_max + + mat_vr = clamp.(mat_vr, -max_length, max_length) + mat_vc = clamp.(mat_vc, -max_length, max_length) + + for row = 1:stride:height, col = 1:stride:width + reverser ? push!(row_s, row) : push!(row_s, -row) + reversec ? push!(col_s, -col) : push!(col_s, col) + + push!(v_rows, mat_vr[row, col]) + push!(v_cols, mat_vc[row, col]) + end + + display( + fliprc ? quiver(row_s, col_s, quiver = (v_rows, v_cols)) : + quiver(col_s, row_s, quiver = (v_cols, v_rows)), + ) + return nothing +end + + +""" +$(SIGNATURES) +""" +function plot_as_quiver( + ϕ; + n_steps = 20, + scale = 300, + max_length = N_Pixel_Width / n_steps, + fliprc = false, + reverser = false, + reversec = false, +) + + height, width = size(ϕ) + stride = Int64(ceil(max(height, width) / n_steps)) + + row_s = Float64[] + col_s = Float64[] + Δrow = Float64[] + Δcol = Float64[] + + for row = 1:stride:height, col = 1:stride:width + dr = (-(ϕ[row, col+1] - ϕ[row, col]) * scale) + dc = (ϕ[row+1, col] - ϕ[row, col]) * scale + + reverser ? push!(row_s, col) : push!(row_s, width + 1 - col) + reversec ? push!(col_s, height + 1 - row) : push!(col_s, row) + reverser && (dr = -dr) + reversec && (dc = -dc) + + push!(Δrow, clamp(dr, -maxlength, max_length)) + push!(Δcol, clamp(dc, -maxlength, max_length)) + end + + # By default, we flip to match the original image. + display( + fliprc ? quiver(row_s, col_s, quiver = (Δrow, Δcol)) : + quiver(col_s, row_s, quiver = (Δcol, Δrow)), + ) + + return nothing +end + + +""" +$(SIGNATURES) +""" +function plot_velocities_as_quiver( + vr, + vc; + n_steps = 20, + scale = 300, + max_length = N_Pixel_Width / n_steps, +) + + height, width = size(vr) + stride = Int64(ceil(max(height, width) / n_steps)) + + row_s = Float64[] + col_s = Float64[] + Δrow = Float64[] + Δcol = Float64[] + + dr = max.(vr, 0.001) + dc = max.(vc, 0.001) + + for row = 1:stride:height, col = 1:stride:width + push!(row_s, height - row) + push!(col_s, col) + + push!(Δrow, dr[row, col]) + push!(Δcol, dc[row, col]) + end + + # readline() + display(quiver(row_s, col_s, quiver = (Δrow, Δcol), aspect_ratio = :equal)) +end + + + +""" +$(SIGNATURES) +""" +function save_plot_scalar_field!(scalar_field, filename) + blue = zeros(Float64, size(scalar_field)) + red = zeros(Float64, size(scalar_field)) + green = zeros(Float64, size(scalar_field)) + + height, width = size(scalar_field) + + for r ∈ 1:height, c ∈ 1:width + field = scalar_field[r, c] + if field < -1.0 + # Green if luminosity error < -1.0 + red[r, c] = 0.0 + blue[r, c] = 0.0 + green[r, c] = 1.0 + + elseif -1.0 <= field < 0.0 + # Red shade if negative + red[r, c] = -field + blue[r, c] = 0.0 + green[r, c] = 0.0 + + elseif 0.0 <= field < 1.0 + # Blue shade if positive + red[r, c] = 0.0 + blue[r, c] = field + green[r, c] = 0.0 + + elseif 1.0 <= field + # White if > 1 + red[r, c] = 1.0 + blue[r, c] = 1.0 + green[r, c] = 1.0 + end + end + + rgbImg = RGB.(red, green, blue) + save("./examples/$(filename).png", rgbImg) +end + + +""" +$(SIGNATURES) +""" +function save_plot_scalar_field(scalar_field, prefix, img) + normalised_D_max = scalar_field + normalised_D_min = scalar_field + # normalised_D_max = scalar_field ./ maximum(scalar_field) + # normalised_D_min = scalar_field ./ minimum(scalar_field) + + blue = zeros(size(scalar_field)) + red = zeros(size(scalar_field)) + green = zeros(size(scalar_field)) + + blue[scalar_field.>0] = normalised_D_max[scalar_field.>0] + red[scalar_field.<0] = -normalised_D_min[scalar_field.<0] + + rgbImg = map(clamp01nan, RGB.(red, green, blue)) + save("./examples/$(prefix)_loss.png", map(clamp01nan, rgbImg)) + + println("Saving output image:") + println(typeof(img)) + E = Gray.(scalar_field) + println(typeof(E)) + outputImg = img - E + save("./examples/$(prefix)_actual.png", outputImg) +end diff --git a/src/scratch2.jl b/src/scratch2.jl new file mode 100644 index 0000000..0951fb2 --- /dev/null +++ b/src/scratch2.jl @@ -0,0 +1,29 @@ +struct P2D + x::AbstractMatrix{Real} + y::AbstractMatrix{Real} + + function P2D(dim1, dim2) + x = Matrix{Float64}(undef, dim1, dim2) + fill!(x, 0.) + y = Matrix{Float64}(undef, dim1, dim2) + fill!(y, 0.) + + return new(x, y) + end +end + + + +m = P2D(2, 3) + +m.x +m.y + +m.x .+= 2.0 + +function f1(v::P2D, a) + v.x .+= a +end + +f1(m, 3.) + diff --git a/src/scratchpad.jl b/src/scratchpad.jl index ea8b4e5..049175c 100644 --- a/src/scratchpad.jl +++ b/src/scratchpad.jl @@ -1,8 +1,528 @@ -using Images, CausticsEngineering +using Revise, Debugger +using Images, Plots +gr() +using CausticsEngineering + +image = Images.load("./examples/personal/slashdot.jpg"); # Check current working directory with pwd() + +image = Images.load("./examples/cat_posing.jpg"); # Check current working directory with pwd() +image = Images.load("./examples/personal/goose.jpg"); # Check current working directory with pwd() + +image = Images.load("./examples/personal/salvador_dali_1.jpg"); # Check current working directory with pwd() +image = Images.load("./examples/personal/salvador_dali_2.jpg"); # Check current working directory with pwd() + +image = Images.load("./examples/personal/statue_of_liberty_1.jpg"); # Check current working directory with pwd() +image = Images.load("./examples/personal/statue_of_liberty_2.jpg"); # Check current working directory with pwd() +image = Images.load("./examples/personal/statue_of_liberty_3.jpg"); # Check current working directory with pwd() + +image = Images.load("./examples/personal/caricature.jpg"); # Check current working directory with pwd() +image = Images.load("./examples/personal/portrait.jpg"); # Check current working directory with pwd() +image = Images.load("./examples/personal/bilal.jpg"); # Check current working directory with pwd() + +mesh, imageBW = engineer_caustics(image); + +imageBW = Float64.(Gray.(image)); +Gray.(imageBW) +sum(imageBW) +sum(imageBW) / length(imageBW) + +CausticsEngineering.average(imageBW) + + +new_img = CausticsEngineering.get_lens_pixels_area(mesh); +luminosity_ratio = sum(imageBW) / sum(new_img) +Gray.(new_img * luminosity_ratio) + + +ε = Float64.(lens_pixels_area - image) +ε .-= average(ε) + +# Save the loss image as a png +save_plot_scalar_field!(ε, "loss_" * prefix) + + + +mesh, imageBW = CausticsEngineering.original_engineer_caustics(image); +imageBW = Float64.(Gray.(image)) +Gray.(imageBW) + +new_img = CausticsEngineering.get_lens_pixels_area(mesh) +luminosity_ratio = sum(imageBW) / sum(new_img) +Gray.(new_img * luminosity_ratio) + + + + +imageBW = Float64.(Gray.(image)); +average_energy_caustics = CausticsEngineering.average(imageBW) +imageBW = imageBW / average_energy_caustics; + +(sum(imageBW) / length(imageBW) - 1.0) * 1e12 + + +# Check a few values to make sure they make sense +mesh.corners.r +mesh.corners.vr +minimum(mesh.corners.vr) +maximum(mesh.corners.vr) + +mesh.corners.c +mesh.corners.vc +minimum(mesh.corners.vc) +maximum(mesh.corners.vc) + +mesh.corners.ϕ +minimum(mesh.corners.ϕ) +maximum(mesh.corners.ϕ) +mean_ϕ = sum(mesh.corners.ϕ) / length(mesh.corners.ϕ) + +# Plot the last vector field +p = plot_as_quiver(mesh, n_steps = 60, scale = 5.0, max_length = 20) + + + +# Generate dictionaries and arrays containing the vertices and triangles. +triangle_dict, vertex_dict, triangle_index, vertex_index = + create_solid(mesh; bottom_distance = Bottom_Offset, top_distance = Top_Offset) + + +# Can be done in 2 steps +# +# triangle_dict, vertex_dict = create_solid_as_dict( +# mesh; +# bottom_distance = Bottom_Offset, +# top_distance = Top_Offset +# ) +# +# # Generate array or coordinates and indices usable to create traditional mesh objects +# _, _, triangle_index, vertex_index = create_solid( +# triangle_dict, +# vertex_dict; +# bottom_distance = Bottom_Offset, +# top_distance = Top_Offset, +# ) + +# Area of all the distorted pixels +am = CausticsEngineering.get_area_corners(mesh) +maximum(am) +minimum(am) + +# Save as an obj file +save_obj!( + triangle_index, + vertex_index, + "./examples/result_cat_mesh.obj", + scale = Float64(Meters_Per_Pixel), + scaleh = Float64(Meters_Per_Pixel), +) + + +# Convert to a mesh structure +using Meshes + +# Zone to plot +row_min = 100.0 +row_max = 200.0 +col_min = 100.0 +col_max = 200.0 + +function clip_trangle_index( + triangle_index, + vertex_index; + row_min = 100.0, + row_max = 200.0, + col_min = 100.0, + col_max = 200.0, +) + + function is_in_zone(t)::Bool + p1, p2, p3 = t + v1_r, v1_c = vertex_index[p1] + v2_r, v2_c = vertex_index[p2] + v3_r, v3_c = vertex_index[p3] + + return row_min <= v1_r <= row_max && + col_min <= v1_c <= col_max && + row_min <= v2_r <= row_max && + col_min <= v2_c <= col_max && + row_min <= v3_r <= row_max && + col_min <= v3_c <= col_max + end + + return [t for t ∈ triangle_index if is_in_zone(t)] +end + +# Select relevant triangles +triangle_index3D = clip_trangle_index(triangle_index, vertex_index) +connections = Meshes.connect.(triangle_index3D) + +# 3D Mesh +vertices3D = Meshes.Point3.(vertex_index) +cat_mesh3D = Meshes.SimpleMesh(vertices3D, connections) + +# 2D Mesh +vertex_index2D = [(v[1], v[2]) for v ∈ vertex_index] +vertices2D = Meshes.Point2.(vertex_index2D) +cat_mesh2D = Meshes.SimpleMesh(vertices2D, connections) + + +# NOT WORKING... +using MeshViz + +# Explicitly import the preferred backend (depending on use case) +import GLMakie +MeshViz.viz(cat_mesh3D) +MeshViz.viz(cat_mesh2D) + + + + + + + + + +imageBW = Float64.(Gray.(image)) + +# Set global size parameters +height, width = size(imageBW) +global N_Pixel_Height = height +global N_Pixel_Width = width +global Meters_Per_Pixel = Caustics_Long_Side / N_Pixel_Height +println("Image size: $(N_Pixel_Height) x $(N_Pixel_Width)") + +# mesh is the same size as the image with an extra row/column to have coordinates to +# cover each image corner with a triangle. +mesh = FaceMesh(N_Pixel_Height, N_Pixel_Width) + +# The total energy going through the lens is equal to the amount of energy on the caustics +total_energy_lens = N_Pixel_Height * N_Pixel_Width * 1.0 # 1 unit of energy per pixel +total_energy_caustics = sum(imageBW) +average_energy_per_pixel = average(imageBW) + +# imageBW is normalised to the same (sort of) _energy_ as the original image. +imageBW = imageBW / average_energy_per_pixel + +lens_pixels_area = get_lens_pixels_area(mesh) + +# Positive error => the triangle needs to shrink (less light) +error_luminosity = Float64.(lens_pixels_area - imageBW) +error_luminosity = error_luminosity .- average(error_luminosity) + +# Save the loss image as a png +println( + """ +Luminosity: + Max/min pixel areas: $(maximum(lens_pixels_area)) / $(minimum(lens_pixels_area)) + Average pixel areas: $(average_absolute(lens_pixels_area)) + Average pixel areas: $(average(lens_pixels_area)) + + Max/min luminosity error: $(maximum(error_luminosity)) / $(minimum(error_luminosity)) + Average abs error: $(average_absolute(error_luminosity)) + Average error: $(average(error_luminosity)) + """, +) + + + +#################################################################################################################### +## STEPPING + +using Revise, Debugger, Images, Plots; +gr() + +using CausticsEngineering + + +######################## Load image +image = Images.load("./examples/personal/slashdot.jpg"); # Check current working directory with pwd() +imageBW = Float64.(Gray.(image)); +imageBW /= average(imageBW); +sum(imageBW) + +h, w = size(imageBW) + + +######################## Create mesh +mesh = CausticsEngineering.FaceMesh(h, w); +CausticsEngineering.field_summary(mesh.corners.ϕ) + +origmesh = CausticsEngineering.OriginalMesh(w + 1, h + 1); +origϕ = zeros(w, h); +CausticsEngineering.field_summary(origϕ) + + +######################## Lens luminosity +area_distorted_corners = CausticsEngineering.get_lens_pixels_area(mesh); +CausticsEngineering.field_summary(area_distorted_corners) + +origarea_distorted_corners = CausticsEngineering.getPixelArea(origmesh); +CausticsEngineering.field_summary(origarea_distorted_corners) + +CausticsEngineering.field_summary(area_distorted_corners - origarea_distorted_corners') + + +######################## Laplacian +Lϕ = CausticsEngineering.laplacian(mesh.corners.ϕ); +CausticsEngineering.field_summary(Lϕ) + +function orig_laplacian(ϕ) + width, height = size(ϕ) + Lϕ = zeros(Float64, size(ϕ)) + + for x = 1:width, y = 1:height + val_up = (y == 1) ? 0.0 : ϕ[x, y-1] + val_down = (y == height) ? 0.0 : ϕ[x, y+1] + val_left = (x == 1) ? 0.0 : ϕ[max(1, x - 1), y] + val_right = (x == width) ? 0.0 : ϕ[x+1, y] + + Lϕ[x, y] = val_up + val_down + val_left + val_right - 4 * ϕ[x, y] + + # if Lϕ[x, y] != 0.0 + # println(val_up, val_down, val_left, val_right, ϕ[x, y], x, y) + # break + # end + end + return Lϕ +end + +origLϕ = orig_laplacian(origϕ); +CausticsEngineering.field_summary(origLϕ) + + +######################## Error luminosity +error_luminosity = zeros(Float64, h + 1, w + 1); +error_luminosity[1:end-1, 1:end-1] = Float64.(area_distorted_corners - imageBW); +CausticsEngineering.field_summary(error_luminosity) + +origerror_luminosity = zeros(Float64, w + 1, h + 1); +origerror_luminosity[1:end-1, 1:end-1] = Float64.(origarea_distorted_corners - imageBW'); +CausticsEngineering.field_summary(origerror_luminosity) + +CausticsEngineering.field_summary(error_luminosity - origerror_luminosity') + + +######################## Normalised error luminosity +error_luminosity = error_luminosity .- average(error_luminosity); +CausticsEngineering.field_summary(error_luminosity) + +origerror_luminosity = origerror_luminosity .- average(origerror_luminosity); +CausticsEngineering.field_summary(origerror_luminosity) + +CausticsEngineering.field_summary(error_luminosity - origerror_luminosity') + + +######################## Distance to Laplacian & propagate +max_update, Lϕ, δ = CausticsEngineering.propagate_poisson!(mesh.corners.ϕ, error_luminosity); +CausticsEngineering.field_summary(mesh.corners.ϕ) +CausticsEngineering.field_summary(Lϕ) +CausticsEngineering.field_summary(δ) + +max_update, origLϕ, origδ = + CausticsEngineering.orig_propagate_poisson!(origϕ, origerror_luminosity); +CausticsEngineering.field_summary(origϕ) +CausticsEngineering.field_summary(origLϕ) +CausticsEngineering.field_summary(origδ) + +CausticsEngineering.field_summary(mesh.corners.ϕ[1:end-1, 1:end-1] - origϕ') +CausticsEngineering.field_summary(mesh.corners.ϕ[2:end, 1:end-1] - origϕ') +CausticsEngineering.field_summary(mesh.corners.ϕ[1:end-1, 2:end] - origϕ') +CausticsEngineering.field_summary(mesh.corners.ϕ[2:end, 2:end] - origϕ') +m = + ( + mesh.corners.ϕ[1:end-1, 1:end-1] + + mesh.corners.ϕ[2:end, 1:end-1] + + mesh.corners.ϕ[1:end-1, 2:end] + + mesh.corners.ϕ[2:end, 2:end] + ) / 4.0; +CausticsEngineering.field_summary(m - origϕ') + + + +δ = Lϕ - error_luminosity; +CausticsEngineering.field_summary(δ) + +δ .*= 0.50; +ϕ = mesh.corners.ϕ[1:end-1, 1:end-1]; +ϕ .+= δ; +mesh.corners.ϕ[1:end-1, 1:end-1] .= ϕ .- average(ϕ); +CausticsEngineering.field_summary(mesh.corners.ϕ) + + + + +mesh = CausticsEngineering.FaceMesh(height, width); +CausticsEngineering.field_summary( + CausticsEngineering.laplacian(mesh.corners.ϕ) - error_luminosity, +) + + + +origϕ = zeros(w, h); + + + + + + + +mesh = CausticsEngineering.FaceMesh(height, width); +new_update = 10_000 +old_update = 2 * new_update +for i ∈ 1:1_000 + old_update = new_update + new_update, Lϕ, δ = + CausticsEngineering.propagate_poisson!(mesh.corners.ϕ, error_luminosity) + i % 50 == 0 && println(""" + $(i) => $(max_update) + $(CausticsEngineering.field_summary(mesh.corners.ϕ)) + $(CausticsEngineering.field_summary(Lϕ)) + $(CausticsEngineering.field_summary(δ)) + """) + if new_update > old_update || i == 5_000 + new_update, Lϕ, δ = + CausticsEngineering.orig_propagate_poisson!(origϕ, origerror_luminosity) + println(""" + $(i) => $(new_update) + $(CausticsEngineering.field_summary(mesh.corners.ϕ)) + $(CausticsEngineering.field_summary(Lϕ)) + $(CausticsEngineering.field_summary(δ)) + """) + break + end +end + +origϕ = zeros(width, height); +for i ∈ 1:5_000 + max_update, origLϕ, origδ = + CausticsEngineering.orig_propagate_poisson!(origϕ, origerror_luminosity) + i % 500 == 0 && println(""" + $(i) => $(max_update) + $(CausticsEngineering.field_summary(origϕ)) + $(CausticsEngineering.field_summary(origLϕ)) + $(CausticsEngineering.field_summary(origδ)) + """) + if max_update <= 1e-3 + max_update, origLϕ, origδ = + CausticsEngineering.orig_propagate_poisson!(origϕ, origerror_luminosity) + println(CausticsEngineering.field_summary(origϕ)) + println(CausticsEngineering.field_summary(origLϕ)) + println(CausticsEngineering.field_summary(origδ)) + break + end +end + + + + + + + + + + +CausticsEngineering.save_plot_scalar_field!(error_luminosity, "trace_1", imageBW) + +max_update = Inf +old_update = Inf + +any(isnan.(error_luminosity)) + +any(isnan.(mesh.corners.r)) +any(isnan.(mesh.corners.c)) +any(isnan.(mesh.corners.vr)) +any(isnan.(mesh.corners.vc)) + +any(isnan.(mesh.corners.ϕ)) + +i = 0 +begin + i += 1 + old_update = max_update + max_update = CausticsEngineering.propagate_poisson!(mesh.corners.ϕ, error_luminosity) + # i % 100 == 0 && println(max_update) + # (abs(old_update - max_update) <= 1e-6 || abs(max_update) <= 1e-4) && break +end + +max_update +old_update + + +f() + + +∇ϕᵤ, ∇ϕᵥ = CausticsEngineering.∇(mesh.corners.ϕ) + +sum(∇ϕᵤ) +average(∇ϕᵤ) +average_absolute(∇ϕᵤ) + + +height, width = size(imageBW) +mesh.corners.vr .= -∇ϕᵤ +mesh.corners.vc .= -∇ϕᵥ +list_triangles = vcat( + [ + CausticsEngineering.triangle3D(mesh, row, col, :top) for row ∈ 1:height, + col ∈ 1:width + ], + [ + CausticsEngineering.triangle3D(mesh, row, col, :bottom) for row ∈ 1:height, + col ∈ 1:width + ], +) + +list_maximum_t = [ + time for time ∈ CausticsEngineering.find_maximum_t.(list_triangles) if + !isnothing(time) && time > 0.0 +] +minimum(list_maximum_t) +maximum(list_maximum_t) + +δ = minimum(list_maximum_t) / 2.0 + +mesh.corners.r += δ * ∇ϕᵤ +mesh.corners.c += δ * ∇ϕᵥ -image = Images.load("./examples/cat_posing.jpg") # Check current working directory with pwd() -mesh, image_caustics = engineer_caustics(image); nothing; -image_caustics' + +function my_laplacian(ϕ::AbstractMatrix{Float64}) + height, width = size(ϕ) + + # Embed matrix within a larger matrix for better vectorization and avoid duplicated code + # Padded matrix adds 1 row/col after the size of ϕ. + # ϕ is inserted in padded matrix within 1:1+height+1 x 1:1+width+1. + # The rest of the padded matrix (its borders) are set at 0.0. + padded_ϕ = zeros(Float64, 1 + height + 1, 1 + width + 1) + padded_ϕ[begin+1:end-1, begin+1:end-1] .= ϕ[1:end, 1:end] + + # Convolution = up + down + left + right + ∇²ϕ = + padded_ϕ[begin+1-1:end-1-1, begin+1:end-1] + + padded_ϕ[begin+1+1:end-1+1, begin+1:end-1] + + padded_ϕ[begin+1:end-1, begin+1-1:end-1-1] + + padded_ϕ[begin+1:end-1, begin+1+1:end-1+1] - + 4.0 * padded_ϕ[begin+1:end-1, begin+1:end-1] + + return ∇²ϕ[1:end-1, 1:end-1] +end + + + + +m = [[1.0 2.0 3.0 2.0]; [4.0 5.0 6.0 5.0]; [7.0 8.0 9.0 8.0]; [10.0 11.0 12.0 11.0]] +h, w = size(m) +padded_m = zeros(Float64, 1 + h + 1, 1 + w + 1) +padded_m[begin+1:end-1, begin+1:end-1] .= m[1:end, 1:end] +padded_m + +padded_m[begin+1-1:end-1-1, begin+1:end-1] + +padded_m[begin+1+1:end-1+1, begin+1:end-1] + +padded_m[begin+1:end-1, begin+1-1:end-1-1] + +padded_m[begin+1:end-1, begin+1+1:end-1+1] - 4.0 * padded_m[begin+1:end-1, begin+1:end-1] + + +sum(CausticsEngineering.laplacian(m)) +sum(orig_laplacian(m)) diff --git a/src/test_functions.jl b/src/test_functions.jl new file mode 100644 index 0000000..275836a --- /dev/null +++ b/src/test_functions.jl @@ -0,0 +1,128 @@ +using CausticsEngineering +using Test + +""" +$(SIGNATURES) + +""" +function testSquareMesh!() + mesh = create_mesh(100, 50) + + println(mesh.nodeArray[1, 1]) + println(mesh.nodes[1]) + + mesh.nodeArray[1, 1].x = 8 + println(mesh.nodeArray[1, 1]) + println(mesh.nodes[1]) + + mesh.nodes[1].y += 12 + println(mesh.nodeArray[1, 1]) + println(mesh.nodes[1]) +end + + +""" +$(SIGNATURES) + +TO REFACTOR. +""" +function testSolidify!() + println("Testing solidification") + width = 100 + height = 100 + origMesh = create_mesh(width, height) + + for y = 1:height, x = 1:width + x2 = (x - width / 2) / width + y2 = (y - height / 2) / height + value = x2 * x2 + y2 * y2 + origMesh.nodeArray[x, y].z = 15 - value * 25 + end + + save_stl!(origMesh, "./examples/testSolidify.obj") + solidMesh = solidify(origMesh, 0) + save_stl!(solidMesh, "./examples/testSolidify2.obj") +end + + + + + +function testAreaTriangle(x1, y1, x2, y2, x3, y3) + a = sqrt((x2 - x1)^2 + (y2 - y1)^2) + b = sqrt((x3 - x2)^2 + (y3 - y2)^2) + c = sqrt((x1 - x3)^2 + (y1 - y3)^2) + s = (a + b + c) / 2.0 + + surface_sq = (a + b + c) * (-a + b + c) * (a - b + c) * (a + b - c) / 16.0 + return surface_sq <= 1e-100 ? 0.0 : surface_sq +end + +function doTestAreaTriangle(N = 10_000) + X1 = 10_000 .* rand(N) + Y1 = 10_000 .* rand(N) + X2 = 10_000 .* rand(N) + Y2 = 10_000 .* rand(N) + X3 = 10_000 .* rand(N) + Y3 = 10_000 .* rand(N) + + for x1 ∈ X1, y1 ∈ Y1, x2 ∈ X2, y2 ∈ Y2, x3 ∈ X3, y3 ∈ Y3 + + _ = testAreaTriangle(x1, y1, x2, y2, x3, y3) > 0.0 + end +end + +@testset "triangle area" begin + for x1 ∈ 10_000 .* rand(10_000), + y1 ∈ 10_000 .* rand(10_000), + x2 ∈ 10_000 .* rand(10_000), + y2 ∈ 10_000 .* rand(10_000), + x3 ∈ 10_000 .* rand(10_000), + y3 ∈ 10_000 .* rand(10_000) + + @test testAreaTriangle(x1, y1, x2, y2, x3, y3) >= 0.0 + end +end + +@testset "More triangle area" begin + @test area(0, 1, 1, 0, 0, 0) +end + + +v1 = CausticsEngineering.Vertex3D(0.0, 1.0, 0.0, 0.0, 0.0) +v2 = CausticsEngineering.Vertex3D(1.0, 0.0, 0.0, 0.0, 0.0) +v3 = CausticsEngineering.Vertex3D(0.0, 0.0, 0.0, 0.0, 0.0) +CausticsEngineering.area(v1, v2, v3) + +v1 = CausticsEngineering.Vertex3D(0.0, 2.0, 0.0, 0.0, 0.0) +v2 = CausticsEngineering.Vertex3D(2.0, 0.0, 0.0, 0.0, 0.0) +v3 = CausticsEngineering.Vertex3D(0.0, 0.0, 0.0, 0.0, 0.0) +CausticsEngineering.area(v1, v2, v3) + +fm = CausticsEngineering.FaceMesh(512, 512); +fm.toptriangles[1, 1] +CausticsEngineering.triangle3D(fm, 1, 1, :top) +CausticsEngineering.area(CausticsEngineering.triangle3D(fm, 1, 1, :top)...) +CausticsEngineering.area(CausticsEngineering.triangle3D(fm, 1, 1, :bottom)...) + +CausticsEngineering.area(CausticsEngineering.triangle3D(fm, 126, 13, :top)...) +CausticsEngineering.area(CausticsEngineering.triangle3D(fm, 126, 13, :bottom)...) + + + +v1 = [7250.768146661309, 4849.994680605491] +v2 = [1725.2665442814496, 8169.783790398901] +v3 = [6560.1725812433, 5264.912964907599] + +a = sqrt(sum((v2 - v1) .^ 2)) +b = sqrt(sum((v3 - v2) .^ 2)) +c = sqrt(sum((v1 - v3) .^ 2)) +s = (a + b + c) / 2.0 + +s * (s - a) * (s - b) * (s - c) +(a + b + c) * (-a + b + c) * (a - b + c) * (a + b - c) / 16.0 + +b + c +a + +10^8 * (b + c - a) diff --git a/src/utilities.jl b/src/utilities.jl index 0cc97f4..65af945 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -1,96 +1,374 @@ """ $(TYPEDEF) + +Coordinates as row (vertical), columns (horizontal). Represent the data of each topleft corner of a pixel. +vr, vc: velocity vector of that corner. +Warning: posts/fences: One more corner than number of pixels. """ -mutable struct Point3D - # Coordinates - x::Float64 - y::Float64 - z::Float64 +struct Vertex3D + r::Float64 + c::Float64 + ϕ::Float64 + + vr::Float64 + vc::Float64 - # Velocity at those coordinates - ix::Int - iy::Int + Vertex3D(r, c, ϕ, vr, vc) = new(r, c, ϕ, vr, vc) + Vertex3D() = new(0.0, 0.0, 0.0, 0.0, 0.0) end +""" +$(SIGNATURES) +""" +function dist(p1::Tuple{Float64,Float64,Float64}, p2::Tuple{Float64,Float64,Float64}) + dr = p2[1] - p1[1] + dc = p2[2] - p1[2] + return √(dr^2 + dc^2) +end """ $(SIGNATURES) """ -function dist(p1::Point3D, p2::Point3D) - dx = p2.x - p1.x - dy = p2.y - p1.y - return sqrt(dx * dx + dy * dy) +function dist(p1::Vertex3D, p2::Vertex3D) + dr = p2.r - p1.r + dc = p2.c - p1.c + return √(dr^2 + dc^2) end """ $(SIGNATURES) - -A midpoint is the average between two points """ -midpoint(p1::Point3D, p2::Point3D) = - Point3D(0.5p1.x + 0.5p2.x, 0.5p1.y + 0.5p2.y, 0.5p1.z + 0.5p2.z, 0, 0) +function area(v1::Vertex3D, v2::Vertex3D, v3::Vertex3D) + a = dist(v1, v2) + b = dist(v2, v3) + c = dist(v3, v1) + s = (a + b + c) / 2.0 + + # surface_sq = (a + b + c) * (- a + b + c) * (a - b + c) * (a + b - c) / 16.0 + surface_sq = s * (s - a) * (s - b) * (s - c) + return surface_sq <= 1e-100 ? 0.0 : √(surface_sq) +end """ $(SIGNATURES) -Barycentre of three points. +Centroid of three points. +""" +centroid(p1::Vertex3D, p2::Vertex3D, p3::Vertex3D) = + Vertex3D((p1.r + p2.r + p3.r) / 3.0, (p1.c + p2.c + p3.c) / 3.0, 0.0, 0.0, 0.0) + + +""" +$(SIGNATURES) """ -centroid(p1::Point3D, p2::Point3D, p3::Point3D) = Point3D( - (p1.x + p2.x + p3.x) / 3.0, - (p1.y + p2.y + p3.y) / 3.0, - (p1.z + p2.z + p3.z) / 3.0, - 0, - 0, -) +function fill_borders!(mat::AbstractMatrix{T}, val::T) where {T} + mat[begin:end, begin] .= val + mat[begin:end, end] .= val + mat[begin, begin:end] .= val + mat[end, begin:end] .= val +end """ $(TYPEDEF) + +## Coordinates + +Origin is top left, going right and down .`x` goes horizontal and follows columns. x is within 1 and width +`y` goes vertical and follows rows. y is within 1 and height. + +# Velocity + +Velocity vector along `row` and `col` (velocity along h not used). + +Implementation is a struc of arrays. Easier to vectorise. + """ -struct Triangle - pt1::Int64 - pt2::Int64 - pt3::Int64 +mutable struct FieldVertex3D + size::Tuple{Int,Int} + + r::AbstractMatrix{Float64} + c::AbstractMatrix{Float64} + vr::AbstractMatrix{Float64} + vc::AbstractMatrix{Float64} - Triangle() = new(0, 0, 0) - Triangle(pt1, pt2, pt3) = new(pt1, pt2, pt3) + # Velocity potential at each corner + ϕ::AbstractMatrix{Float64} + + rows_numbers::AbstractMatrix{Float64} + cols_numbers::AbstractMatrix{Float64} end +# Information about each corner surrounding a pixel. Posts&fences warning!!! +function FieldVertex3D(height, width) + # Two matrices comtaining the row and column number of each element of the matrix + rows_numbers = repeat(Float64.(1:height+1), 1, width + 1) + cols_numbers = repeat(Float64.(1:width+1)', height + 1, 1) + + mr = rows_numbers + mc = cols_numbers + mvr = zeros(Float64, height + 1, width + 1) + mvc = zeros(Float64, height + 1, width + 1) + + mϕ = zeros(Float64, height + 1, width + 1) + + fv = FieldVertex3D((height, width), mr, mc, mvr, mvc, mϕ, rows_numbers, cols_numbers) + reset_border_values!(fv) + return fv +end + + +""" +$(SIGNATURES) + +Reset the values of the borders corners: all positions at integer values, all velocities at 0. +""" +function reset_border_values!(corners::FieldVertex3D) + corners.r[1, :] .= corners.rows_numbers[1, :] + corners.r[end, :] .= corners.rows_numbers[end, :] + + corners.r[:, 1] .= corners.rows_numbers[:, 1] + corners.r[:, end] .= corners.rows_numbers[:, end] + + corners.c[1, :] .= corners.cols_numbers[1, :] + corners.c[end, :] .= corners.cols_numbers[end, :] + + corners.c[:, 1] .= corners.cols_numbers[:, 1] + corners.c[:, end] .= corners.cols_numbers[:, end] + + fill_borders!(corners.vr, 0.0) + fill_borders!(corners.vc, 0.0) +end + + +Base.size(fv::FieldVertex3D) = fv.size + + +""" +$(SIGNATURES) +""" +function Vertex3D(fv::FieldVertex3D, row, col) + height, width = size(fv) + + if (1 <= row <= height + 1) && (1 <= col <= width + 1) + Vertex3D( + fv.r[row, col], + fv.c[row, col], + fv.vr[row, col], + fv.vc[row, col], + fv.ϕ[row, col], + ) + else + missing + end +end + + """ $(TYPEDEF) + +The algorithm works by allocating a rectangle from the exit mesh (the face where the light exits the block +of acrylate). Initially, each rectangle representsa single pixel on that face and faces s single pixel on +the caustics picture. + +Each rectangle is split into 2 trangles: one upper-left, the other bottom-right like: + +``` +*------* +| / | +| / | +| / | +| / | +*------* +``` +Although the corners of each triangle match corners of the rectangles, it is sometimes easier to think in terms +of triangles rather than corners coming from different rectangles. +""" +struct FaceMesh + size::Tuple{Int,Int} + corners::FieldVertex3D + + toptriangles::AbstractMatrix{Tuple{Tuple{Int,Int},Tuple{Int,Int},Tuple{Int,Int}}} + bottriangles::AbstractMatrix{Tuple{Tuple{Int,Int},Tuple{Int,Int},Tuple{Int,Int}}} + + FaceMesh(s, tl, tt, bt) = new(s, tl, tt, bt) + + function FaceMesh(height::Int, width::Int) + m_corner = FieldVertex3D(height, width) + t_triangles = Matrix{Tuple{Tuple{Int,Int},Tuple{Int,Int},Tuple{Int,Int}}}( + undef, + height, + width, + ) + b_triangles = Matrix{Tuple{Tuple{Int,Int},Tuple{Int,Int},Tuple{Int,Int}}}( + undef, + height, + width, + ) + + for row ∈ 1:height, col ∈ 1:width + t_triangles[row, col] = ((row, col), (row + 1, col), (row, col + 1)) + b_triangles[row, col] = ((row, col + 1), (row + 1, col), (row + 1, col + 1)) + end + + return FaceMesh((height, width), m_corner, t_triangles, b_triangles) + end +end + + +""" +$(SIGNATURES) + +Return the dimension of the mesh as the number of rectangles. It does not return the number of corner points. """ -struct Mesh - nodes::Vector{Point3D} - nodeArray::Matrix{Point3D} - triangles::Vector{Triangle} - width::Int - height::Int +Base.size(mesh::FaceMesh) = mesh.size + + +""" +$(SIGNATURES) + +Converts a triangle as a triplet of references to mesh vertices to a triplet of 3D coordinates. +""" +function triangle3D(mesh::FaceMesh, row::Int, col::Int, side = Union{:top,:bottom}) + height, width = size(mesh) + max_distance_squared = (height + 1)^2 + (width + 1)^2 + + if 1 <= row <= height && 1 <= col <= width + if side == :top + t = mesh.toptriangles[row, col] + elseif side == :bottom + t = mesh.bottriangles[row, col] + end + + t1_row, t1_col = t[1] + p1 = Vertex3D(mesh.corners, t1_row, t1_col) + + t2_row, t2_col = t[2] + p2 = Vertex3D(mesh.corners, t2_row, t2_col) + + t3_row, t3_col = t[3] + p3 = Vertex3D(mesh.corners, t3_row, t3_col) + + return (p1, p2, p3) + else + # return (missing, missing, missing) + @assert false "Coordinates out of bounds at $(row), $(col)" + end end +triangle3D(mesh::FaceMesh, ci::CartesianIndex{2}, side = Union{:top,:bottom}) = + triangle3D(mesh, ci[1], ci[2], side) + """ $(SIGNATURES) + +Centroid of a specific mesh triangle. """ -function triangle_area(mesh::Mesh, index::Int) - triangle = mesh.triangles[index] - pt1 = mesh.nodes[triangle.pt1] - pt2 = mesh.nodes[triangle.pt2] - pt3 = mesh.nodes[triangle.pt3] +function centroid(mesh::FaceMesh, ci::CartesianIndex; side = Union{:top,:bottom}) + if side == :top + return centroid(top_triangle3D(mesh, ci)...) + elseif side == :bottom + return centroid(bot_triangle3D(mesh, ci)...) + else + return missing + end +end + +""" +$(SIGNATURES) - return triangle_area(pt1, pt2, pt3) +Centroid of a specific mesh triangle. +""" +function centroid(mesh::FaceMesh, row::Int, col::Int; side = Union{:top,:bottom}) + if side == :top + return centroid(top_triangle3D(mesh, row, col)...) + elseif side == :bottom + return centroid(bot_triangle3D(mesh, row, col)...) + else + return missing + end end +""" +$(SIGNATURES) +""" +function area(mesh::FaceMesh, ci::CartesianIndex; side = Union{:top,:bottom}) + if side == :top + return area(top_triangle3D(mesh, ci)...) + elseif side == :bottom + return area(top_triangle3D(mesh, ci)...) + else + return missing + end +end """ $(SIGNATURES) """ -function triangle_area(p1::Point3D, p2::Point3D, p3::Point3D) - a = dist(p1, p2) - b = dist(p2, p3) - c = dist(p3, p1) - s = (a + b + c) / 2.0 +function area(mesh::FaceMesh, row::Int, col::Int; side = Union{:top,:bottom}) + if side == :top + return area(top_triangle3D(mesh, row, col)...) + elseif side == :bottom + return area(top_triangle3D(mesh, row, col)...) + else + return missing + end +end - return sqrt(s * (s - a) * (s - b) * (s - c)) +""" +$(SIGNATURES) +""" +area(t::Tuple{Vertex3D,Vertex3D,Vertex3D}) = area(t...) + + +""" +$(SIGNATURES) + +A Mesh is a collection of triangles. The brightness flowing through a given triangle is just proportional to its +area in the r, r plane. h is ignored. + +The function returns a matrix with the quantity of light coming from each 'rectangle' around a corner. That 'rectangle' +has been shifted and flexed around. +""" +function get_lens_pixels_area(mesh::FaceMesh) + height, width = size(mesh) + + top_tri_area = + area.([triangle3D(mesh, row, col, :top) for row ∈ 1:height, col ∈ 1:width]) + bot_tri_area = + area.([triangle3D(mesh, row, col, :bottom) for row ∈ 1:height, col ∈ 1:width]) + + total_area = top_tri_area + bot_tri_area + + average_energy_per_pixel = average(total_area) + luminosity_pixels = total_area / average_energy_per_pixel + + @assert abs(sum(luminosity_pixels) - height * width) < 1.0 "Total lens luminosity is incorrect" + + return luminosity_pixels end + + +""" +$(SIGNATURES) + +""" +function field_summary(field, fieldname) + s = round(sum(field), sigdigits = 6) + max = round(maximum(field), sigdigits = 4) + min = round(minimum(field), sigdigits = 4) + avg = round(average(field), sigdigits = 4) + abs = round(average_absolute(field), sigdigits = 4) + + return """$(fieldname) (Sum/Max/Min/Avg/Avg abs.): $(s) / $(max) / $(min) / $(avg) / $(abs)""" +end + + +""" +$(SIGNATURES) + +""" +field_summary(field) = field_summary(field, "") diff --git a/test/runtests.jl b/test/runtests.jl index 406e8cb..0c068d9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,4 @@ using CausticsEngineering using Test -@testset "CausticsEngineering.jl" begin - # Write your tests here. -end +include("triangles.jl") diff --git a/test/triangles.jl b/test/triangles.jl new file mode 100644 index 0000000..e7aebf8 --- /dev/null +++ b/test/triangles.jl @@ -0,0 +1,42 @@ + + +@testset "Find maximum time to nil triangle" begin + t3 = ( + CausticsEngineering.Vertex3D(5.0, 5.0, 0.0, 0.0, 0.0), + CausticsEngineering.Vertex3D(6.0, 5.0, 0.0, 0.0, 0.0), + CausticsEngineering.Vertex3D(5.5, 5.5, 0.0, 0.0, -3.0), + ) + @test CausticsEngineering.find_maximum_t(t3)[1] ≈ 0.5 / 3.0 atol = 0.01 + + + t3 = ( + CausticsEngineering.Vertex3D(5.0, 5.0, 0.0, 0.0, 0.0), + CausticsEngineering.Vertex3D(6.0, 5.0, 0.0, 0.0, 0.0), + CausticsEngineering.Vertex3D(5.5, 5.5, 0.0, 0.7, -3.0), + ) + @test CausticsEngineering.find_maximum_t(t3)[1] ≈ 0.5 / 3.0 atol = 0.01 + + + t3 = ( + CausticsEngineering.Vertex3D(5.0, 5.0, 0.0, 0.0, 0.0), + CausticsEngineering.Vertex3D(6.0, 5.0, 0.0, 0.0, 0.0), + CausticsEngineering.Vertex3D(5.5, 5.5, 0.0, -10.0, -3.0), + ) + @test CausticsEngineering.find_maximum_t(t3)[1] ≈ 0.5 / 3.0 atol = 0.01 + + + t3 = ( + CausticsEngineering.Vertex3D(5.0, 5.0, 0.01, 0.1, -2.0), + CausticsEngineering.Vertex3D(6.0, 5.0, 0.01, 0.1, 2.0), + CausticsEngineering.Vertex3D(5.5, 5.5, 0.01, 0.5, -3.0), + ) + @test CausticsEngineering.find_maximum_t(t3)[1] ≈ 0.1540 atol = 0.01 + + t3 = ( + CausticsEngineering.Vertex3D(5.0, 5.0, 0.01, 0.1, -2.0), + CausticsEngineering.Vertex3D(6.0, 5.0, 0.01, 0.1, 2.0), + CausticsEngineering.Vertex3D(5.5, 5.5, 0.01, 0.5, -3.0), + ) + @test CausticsEngineering.find_maximum_t(t3)[1] ≈ 0.1540 atol = 0.01 + +end