From caaf14d5064a37e3fa522d926f420940e16eb712 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Syver=20D=C3=B8ving=20Agdestein?= Date: Tue, 15 Mar 2022 19:00:54 +0100 Subject: [PATCH 1/5] Make 2D convex hull without MiniQHull --- src/geometry/convexhull.jl | 72 +++++++++++++++++++++++++++++++++++++- 1 file changed, 71 insertions(+), 1 deletion(-) diff --git a/src/geometry/convexhull.jl b/src/geometry/convexhull.jl index ef333e0..74db664 100644 --- a/src/geometry/convexhull.jl +++ b/src/geometry/convexhull.jl @@ -11,7 +11,12 @@ function convexhull(points) dim, npoint = size(points) # Compute Delaunay triangulation - D = Int.(delaunay(points)) + if dim == 2 + return convexhull2(points) + else + # D = Int.(delaunay3(points)) + D = Int.(delaunay(points)) + end # Sort each column of d for i = 1:size(D, 2) @@ -107,3 +112,68 @@ end # display(pl) # sleep(0.5) # end + + +""" + delaunay3(points) + +3D Delaunay triangulation from Tetgen. +""" +function delaunay3(points) + input = RawTetGenIO{Cdouble}() + input.pointlist = points + tetgen = tetrahedralize(input, "cQ") + + tetgen.tetrahedronlist +end + + +""" + convexhull2(points) + +Convex hull of 2D points (gift wrapping). +https://en.wikipedia.org/wiki/Gift_wrapping_algorithm +""" +function convexhull2(p) + function swap!(points, i, j) + tmp = points[:, i] + points[:, i] = points[:, j] + points[:, j] = tmp + end + + isclock(u, v) = u[2]v[1] - u[1]v[2] ≥ 0 + + npoint = size(p, 2) + + # First point with lowest x-coordinate + swap!(p, 1, argmin(p[1, :])) + + # First u (turn clockwise from here) + u = [0, 1] + + # Iterate through points + i = 1 + for i = 1:npoint-1 + done = true + jnext = i + for j = i+1:npoint + if isclock(u, p[:, j] - p[:, i]) + u = p[:, j] - p[:, i] + jnext = j + done = false + end + end + swap!(p, i + 1, jnext) + u = p[:, i] - p[:, i+1] + + if done || (i ≥ 2 && isclock(u, p[:, i] - p[:, 1])) + p = p[:, 1:i] + break + end + end + + n = size(p, 2) + e = [(1:n)'; (2:n)' 1] + + e, p +end From 13d0d83026f04be6430e8dc0ae9e479158e044d9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Syver=20D=C3=B8ving=20Agdestein?= Date: Tue, 15 Mar 2022 22:52:18 +0100 Subject: [PATCH 2/5] Add convexhull3 and remove MiniQHull --- Project.toml | 2 - src/SpinDoctor.jl | 1 - src/geometry/convexhull.jl | 241 +++++++++++++++++-------------------- 3 files changed, 111 insertions(+), 133 deletions(-) diff --git a/Project.toml b/Project.toml index c34eabb..3aab661 100644 --- a/Project.toml +++ b/Project.toml @@ -10,7 +10,6 @@ Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4" GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" -MiniQhull = "978d7f02-9e05-4691-894f-ae31a51d76ca" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -30,7 +29,6 @@ DiffEqCallbacks = "2" Expokit = "0.2" GLPK = "0.14, 0.15, 1" Makie = "0.15, 0.16" -MiniQhull = "0.3" OrdinaryDiffEq = "5, 6" Polynomials = "2, 3" QuadGK = "2" diff --git a/src/SpinDoctor.jl b/src/SpinDoctor.jl index dbaa5b9..2221b80 100644 --- a/src/SpinDoctor.jl +++ b/src/SpinDoctor.jl @@ -11,7 +11,6 @@ using Expokit: expmv! using GLPK: Optimizer using LinearAlgebra using Makie -using MiniQhull: delaunay using OrdinaryDiffEq: OrdinaryDiffEq using OrdinaryDiffEq: DiffEqArrayOperator, diff --git a/src/geometry/convexhull.jl b/src/geometry/convexhull.jl index 74db664..84aa23c 100644 --- a/src/geometry/convexhull.jl +++ b/src/geometry/convexhull.jl @@ -1,133 +1,3 @@ -""" - elements, points = convexhull(points) - -Compute the convex hull of a set of points `points` (of size `dim * npoint`). Return a -matrix of boundary elements `elements` (of size`dim * nelement`) and a restriction of the -original points to the boundary (size `dim * npoint_keep`). -""" -function convexhull(points) - - # Sizes - dim, npoint = size(points) - - # Compute Delaunay triangulation - if dim == 2 - return convexhull2(points) - else - # D = Int.(delaunay3(points)) - D = Int.(delaunay(points)) - end - - # Sort each column of d - for i = 1:size(D, 2) - sort!(@view D[:, i]) - end - - # Identify boundary elements. In 3D, those are facets not shared by two elements, in 2D - # those are edges not shared by two facets - if dim == 2 - combs = [[1, 2], [1, 3], [2, 3]] - else - combs = [[1, 2, 3], [1, 2, 4], [1, 3, 4], [2, 3, 4]] - end - A = [SVector{dim}(D[comb, i]) for i = 1:size(D, 2) for comb ∈ combs] - sort!(A) - - # Indices of unique elements - ikeep = fill(false, length(A)) - - # Perform lookup for the first element - lookup = true - - # Loop through facets or edges. Each one occurs exactly once or twice, and they are - # ordered - for i ∈ 1:length(A)-1 - if lookup - # A[i] is different than A[i-1] - if A[i] == A[i+1] - # A[i] and A[i+1] is a double element, and should not be kept. Skip next - # iteration - lookup = false - else - # A[i] is unique, and is thus a boundary element - ikeep[i] = true - end - else - # A[i] is a copy of A[i-1], but A[i+1] is a new element. Perform next iteration - lookup = true - end - end - - # Keep the last element if it is different from the one before - ikeep[end] = lookup - - # Keep unique elements - A = A[ikeep] - - # Order edges to form a connected wrapping in 2D case - if dim == 2 - for i = 1:length(A)-1 - # Find next edge from endpoint of previous edge - j = i + findfirst(e -> A[i][2] ∈ e, @view A[i+1:end]) - - # Reorder - tmp = A[i+1] - A[i+1] = A[j] - A[j] = tmp - - # Orient next edge correctly - if A[i+1][2] == A[i][2] - A[i+1] = SVector{2}(A[i+1][[2, 1]]) - end - end - end - - E = mapreduce(Vector, hcat, A) - - # Important: `unique` preserves order of edges in 2D case - ikeep = unique(E) - - # Inverse point map - ikeep_inv = fill(0, npoint) - for i = 1:length(ikeep) - ikeep_inv[ikeep[i]] = i - end - - # Renumber points from 1 to nkeep - boundary_elements = ikeep_inv[E] - boundary_points = points[:, ikeep] - - boundary_elements, boundary_points -end - -# a = rand(2, 100) -# e, p = convexhull(a) - -# pl = scatter(a[1, :], a[2, :]) - -# scatter!(pl, p[1, :], p[2, :], color = :red) - -# for i = 1:size(e, 2) -# plot!(pl, p[1, e[:, i]], p[2, e[:, i]]) -# display(pl) -# sleep(0.5) -# end - - -""" - delaunay3(points) - -3D Delaunay triangulation from Tetgen. -""" -function delaunay3(points) - input = RawTetGenIO{Cdouble}() - input.pointlist = points - tetgen = tetrahedralize(input, "cQ") - - tetgen.tetrahedronlist -end - - """ convexhull2(points) @@ -177,3 +47,114 @@ function convexhull2(p) e, p end + +""" + convexhull3(points) + +Convex hull of 3D points. Based on +https://github.com/rileyborgard/convex-hull-3d/blob/master/convexhull3d.cpp +""" +function convexhull3(points) + n = size(points, 2) + @assert n ≥ 4 + + points = SVector{3}.(eachcol(points)) + ϵ = 1e-10 + + tet = [1] + for i = 2:n + length(tet) < 4 || break + a = length(tet) == 1 && norm(points[tet[1]] - points[i]) > ϵ + b = + length(tet) == 2 && + norm((points[tet[2]] - points[tet[1]]) × (points[i] - points[tet[1]])) > ϵ + c = + length(tet) == 3 && + abs( + (points[i] - points[tet[1]]) ⋅ + cross(points[tet[2]] - points[tet[1]], points[tet[3]] - points[tet[1]]), + ) > ϵ + if a || b || c + push!(tet, i) + end + end + @assert length(tet) == 4 + + pointsfirst = points[tet] + prepend!(deleteat!(points, tet), pointsfirst) + + faces = SVector{3,Int}[] + + dead = trues(n, n) + function addface(a, b, c) + push!(faces, SVector{3,Int}(a, b, c)) + dead[a, b] = dead[b, c] = dead[c, a] = false + end + + addface(1, 2, 3) + addface(1, 3, 2) + + for i = 4:n + faces2 = SVector{3,Int}[] + for f ∈ faces + a, b, c = f[1], f[2], f[3] + n = (points[b] - points[a]) × (points[c] - points[a]) + n /= norm(n) + if (points[i] - points[a]) ⋅ n > ϵ + dead[a, b] = dead[b, c] = dead[c, a] = true + else + push!(faces2, f) + end + end + empty!(faces) + for f ∈ faces2 + a, b, c = f[1], f[2], f[3] + dead[b, a] && addface(b, a, i) + dead[c, b] && addface(c, b, i) + dead[a, c] && addface(a, c, i) + end + append!(faces, faces2) + end + + faces = reduce(hcat, faces) + points = reduce(hcat, points) + + faces, points +end + +""" + elements, points = convexhull(points) + +Compute the convex hull of a set of points `points` (of size `dim * npoint`). Return a +matrix of boundary elements `elements` (of size`dim * nelement`) and a restriction of the +original points to the boundary (size `dim * npoint_keep`). +""" +function convexhull(points) + dim = size(points, 1) + if dim == 2 + e, p = convexhull2(points) + elseif dim == 3 + e, p = convexhull3(points) + else + error("convexhull only implemented for 2D and 3D") + end + e, p +end + +# # 2D +# a = rand(2, 100) +# e, p = convexhull(a) +# pl = scatter(a[1, :], a[2, :]) +# scatter!(pl, p[1, :], p[2, :], color = :red) +# for i = 1:size(e, 2) +# plot!(pl, p[1, e[:, i]], p[2, e[:, i]]) +# display(pl) +# sleep(0.5) +# end + + +# # 3D +# points = randn(3, 1000) +# e, p = convexhull3(points) +# poly(p', e'; color = p[3, :], transparency = true, strokewidth = 1, shading = false) +# scatter!(points) From fab71475bb60fe1f0e780c3018bb6aa5087cf33b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Syver=20D=C3=B8ving=20Agdestein?= Date: Tue, 15 Mar 2022 22:53:55 +0100 Subject: [PATCH 3/5] Add surface triangulation plot --- examples/driver.jl | 5 +++-- src/SpinDoctor.jl | 5 +++-- src/plot/plot_surfaces.jl | 21 +++++++++++++++++++++ 3 files changed, 27 insertions(+), 4 deletions(-) create mode 100644 src/plot/plot_surfaces.jl diff --git a/examples/driver.jl b/examples/driver.jl index df78a38..adaf754 100644 --- a/examples/driver.jl +++ b/examples/driver.jl @@ -23,14 +23,15 @@ set_theme!(theme_black()) # include("setups/spheres.jl") include("setups/neuron.jl") -mesh, = @time create_geometry(setup; recreate = true); +mesh, surfaces, cells = @time create_geometry(setup; recreate = true); model = Model(; mesh, coeffs...); volumes = get_cmpt_volumes(model.mesh) D_avg = 1 / 3 * tr.(model.D)' * volumes / sum(volumes) @info "Number of nodes per compartment:" length.(model.mesh.points) ## Plot mesh -plot_mesh(model.mesh) +plot_surfaces(surfaces, 1:3) +plot_mesh(model.mesh, 1:1) ## Assemble finite element matrices matrices = @time assemble_matrices(model); diff --git a/src/SpinDoctor.jl b/src/SpinDoctor.jl index 2221b80..f32de87 100644 --- a/src/SpinDoctor.jl +++ b/src/SpinDoctor.jl @@ -49,7 +49,7 @@ export savefield export compute_adc_sta export fit_adc export fit_tensors -export plot_mesh, plot_field, plot_hardi +export plot_surfaces, plot_mesh, plot_field, plot_hardi export BTPDE, HADC, Karger export Laplace, MatrixFormalism, AnalyticalLaplace, AnalyticalMatrixFormalism @@ -164,9 +164,10 @@ include("postprocess/fit_tensors.jl") include("postprocess/savefield.jl") # Plot -include("plot/plot_mesh.jl") include("plot/plot_field.jl") include("plot/plot_hardi.jl") +include("plot/plot_mesh.jl") +include("plot/plot_surfaces.jl") # Callbacks include("callbacks/callbacks.jl") diff --git a/src/plot/plot_surfaces.jl b/src/plot/plot_surfaces.jl new file mode 100644 index 0000000..02a285f --- /dev/null +++ b/src/plot/plot_surfaces.jl @@ -0,0 +1,21 @@ +""" + plot_surfaces(femesh, compartments = 1:ncompartment) + +Plot surfaces. +""" +function plot_surfaces(surfaces, boundaries = 1:maximum(surfaces.facetmarkers)) + (; points, facets, facetmarkers) = surfaces + scene = nothing + first = true + for b ∈ boundaries + f = facets[:, facetmarkers .== b] + color = points[3, :] + if first + scene = poly(points', f'; color, strokewidth = 1, shading = false) + first = false + else + poly!(points', f'; color, strokewidth = 1, shading = false) + end + end + scene +end From db93a762c550ac17dfa47c02d32f95f8b78b0aa4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Syver=20D=C3=B8ving=20Agdestein?= Date: Tue, 15 Mar 2022 22:58:32 +0100 Subject: [PATCH 4/5] Add last example with working sphere --- docs/make.jl | 2 +- examples/custom_gradients.jl | 1 - examples/driver.jl | 4 ++-- src/gradients/gradient.jl | 4 +++- 4 files changed, 6 insertions(+), 5 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 9528ea5..8ad72c8 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,7 +7,7 @@ DocMeta.setdocmeta!(SpinDoctor, :DocTestSetup, :(using SpinDoctor); recursive = # Generate examples examples = [ "Solve BTPDE" => "solve_btpde", - # "Custom gradients" => "custom_gradients", + "Custom gradients" => "custom_gradients", "Compare ADCs" => "compare_adcs", "Matrix Formalism" => "matrix_formalism", "High Angular Resolution" => "hardi", diff --git a/examples/custom_gradients.jl b/examples/custom_gradients.jl index d2973b8..7972d53 100644 --- a/examples/custom_gradients.jl +++ b/examples/custom_gradients.jl @@ -63,7 +63,6 @@ R = [cos(φ) sin(φ) 0; -sin(φ) cos(φ) 0; 0 0 1] g⃗(t) = 1.0 * R * [sin(2π * t / TE), sin(20π * t / TE) / 5, cos(2π * t / TE)] gradient = GeneralGradient(; g⃗, TE) - # In order to follow the evolution of the solution during time stepping, we add a # [`Plotter`](@ref) to a list of callbacks. Other available callbacks are [`Printer`](@ref) # for showing time stepping information, and [`VTKWriter`](@ref) for saving the solution time diff --git a/examples/driver.jl b/examples/driver.jl index adaf754..6959c68 100644 --- a/examples/driver.jl +++ b/examples/driver.jl @@ -18,10 +18,10 @@ set_theme!(theme_black()) ## Create model from setup recipe # include("setups/axon.jl") # include("setups/sphere.jl") -# include("setups/plates.jl") +include("setups/plates.jl") # include("setups/cylinders.jl") # include("setups/spheres.jl") -include("setups/neuron.jl") +# include("setups/neuron.jl") mesh, surfaces, cells = @time create_geometry(setup; recreate = true); model = Model(; mesh, coeffs...); diff --git a/src/gradients/gradient.jl b/src/gradients/gradient.jl index 3c70550..4a6114a 100644 --- a/src/gradients/gradient.jl +++ b/src/gradients/gradient.jl @@ -6,7 +6,7 @@ Magnetic field gradient ``\\vec{g}(t) \\in \\mathbb{R}^3``. abstract type AbstractGradient{T} end """ - GeneralGradient(g⃗) + GeneralGradient(g⃗, TE) General gradient sequence ``\\vec{g}(t) \\in \\mathbb{R}^3``. The direction and amplitude may vary in time. @@ -36,3 +36,5 @@ end (g::GeneralGradient)(t) = g.g⃗(t) (g::ScalarGradient)(t) = g.amplitude * g.profile(t) * g.dir + +Base.show(io::IO, g::ScalarGradient) = "$(g.amplitude) * $(g.profile) * $(g.dir)" From d04ab230984f3974dac5eaf8a36a6b65f5a18520 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Syver=20D=C3=B8ving=20Agdestein?= Date: Tue, 15 Mar 2022 23:03:16 +0100 Subject: [PATCH 5/5] Test no longer broken --- test/workflow.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/workflow.jl b/test/workflow.jl index 21b1e8a..14d8070 100644 --- a/test/workflow.jl +++ b/test/workflow.jl @@ -24,8 +24,8 @@ T = Float64 @testset "SphereSetup" begin setup = get_setup(SphereSetup{T}) coeffs = get_coeffs(setup) - @test_broken mesh, = create_geometry(setup; recreate = true) - @test_broken model = Model(; mesh, coeffs...) + mesh, = create_geometry(setup; recreate = true) + model = Model(; mesh, coeffs...) end @testset "NeuronSetup" begin