Replace MiniQHull

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"

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",

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

examples/driver.jl

@@ -18,19 +18,20 @@ 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, = @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);

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,
@@ -50,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
@@ -165,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")

src/geometry/convexhull.jl

@@ -1,109 +1,160 @@
 """
-    elements, points = convexhull(points)
+    convexhull2(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`).
+Convex hull of 2D points (gift wrapping).
+https://en.wikipedia.org/wiki/Gift_wrapping_algorithm
 """
-function convexhull(points)
+function convexhull2(p)
+    function swap!(points, i, j)
+        tmp = points[:, i]
+        points[:, i] = points[:, j]
+        points[:, j] = tmp
+    end
 
-    # Sizes
-    dim, npoint = size(points)
+    isclock(u, v) = u[2]v[1] - u[1]v[2] ≥ 0
 
-    # Compute Delaunay triangulation
-    D = Int.(delaunay(points))
+    npoint = size(p, 2)
 
-    # Sort each column of d
-    for i = 1:size(D, 2)
-        sort!(@view D[:, i])
-    end
+    # First point with lowest x-coordinate
+    swap!(p, 1, argmin(p[1, :]))
 
-    # 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
+    # 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
-        else
-            # A[i] is a copy of A[i-1], but A[i+1] is a new element. Perform next iteration
-            lookup = true
+        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
 
-    # Keep the last element if it is different from the one before
-    ikeep[end] = lookup
+    n = size(p, 2)
+    e = [(1:n)'; (2:n)' 1]
 
-    # Keep unique elements
-    A = A[ikeep]
+    e, p
+end
 
-    # 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
+"""
+    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
 
-    E = mapreduce(Vector, hcat, A)
+    pointsfirst = points[tet]
+    prepend!(deleteat!(points, tet), pointsfirst)
 
-    # Important: `unique` preserves order of edges in 2D case
-    ikeep = unique(E)
+    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
 
-    # Inverse point map
-    ikeep_inv = fill(0, npoint)
-    for i = 1:length(ikeep)
-        ikeep_inv[ikeep[i]] = i
+    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
 
-    # Renumber points from 1 to nkeep
-    boundary_elements = ikeep_inv[E]
-    boundary_points = points[:, ikeep]
+    faces = reduce(hcat, faces)
+    points = reduce(hcat, points)
+
+    faces, points
+end
 
-    boundary_elements, boundary_points
+"""
+    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)

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)"

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

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