Fix mesh susceptibility

src/geometries/internal/susceptibility/grid.jl

@@ -7,14 +7,42 @@ import ...BoundingBoxes: BoundingBox, lower
 import ..Base: BaseSusceptibility, single_susceptibility, single_susceptibility_gradient
 
 
-struct SusceptibilityGridElement{N}
-    index :: Int32
-    shift :: Int32
+"""
+A simplified representation of a magnetic susceptibility obstruction, which has been simplified as a point source in N-dimensional space.
+
+There are two versions:
+- [`IsotropicSusceptibilityGridElement`](@ref): for objects with isotropic susceptibilities (and symmetric objects with anisotropic susceptibilities), where the magnetisation is consistently in the direction of the B0 field.
+- [`AnisotropicSusceptibilityGridElement`](@ref): for objects with anisotropic susceptibilities, where the magnetisation can be in any arbitrary orientation.
+"""
+abstract type SusceptibilityGridElement{N} end
+
+
+function (::Type{SusceptibilityGridElement{N}})(position::AbstractVector{<:Number}, radius::Number, susceptibility) where {N}
+    if susceptibility isa Number
+        return IsotropicSusceptibilityGridElement{N}(position, radius, susceptibility)
+    else
+        return AnisotropicSusceptibilityGridElement{N}(position, radius, susceptibility)
+    end
+end
+
+"""
+[`SusceptibilityGridElement`](@ref) with a magnetisation in the direction of the B0 field.
+"""
+struct IsotropicSusceptibilityGridElement{N} <: SusceptibilityGridElement{N}
     position :: SVector{N, Float64}
     radius :: Float64
     susceptibility :: Float64
 end
 
+"""
+[`SusceptibilityGridElement`](@ref) with an arbitrarily oriented magnetisation.
+"""
+struct AnisotropicSusceptibilityGridElement{N} <: SusceptibilityGridElement{N}
+    position :: SVector{N, Float64}
+    radius :: Float64
+    susceptibility :: SVector{N, Float64}
+end
+
 
 BoundingBox(element::SusceptibilityGridElement{N}) where {N} = BoundingBox{N}(
     element.position .- radius,
@@ -29,28 +57,33 @@ The sources for which the field contribution should still be added are stored in
 These are stored of vectors of [`SuscetibilityGridElement`](@ref), where each element
 contains all the information to compute the dipole field approximation.
 """
-abstract type SusceptibilityGrid{N, O, K} end
+abstract type SusceptibilityGrid end
+
+const IndexRepeat = @NamedTuple{index::Int32, shift::Int32}
 
 """
 Specialised version of [`SusceptibilityGrid`](@ref) for repeating geometries.
 """
-struct SusceptibilityGridRepeat{N, O, K} <: SusceptibilityGrid{N, O, K}
+struct SusceptibilityGridRepeat{N, O, E, K} <: SusceptibilityGrid
     inv_resolution :: SVector{N, Float64}
     rotation :: SMatrix{N, 3, Float64, K}
 
     off_resonance :: Array{Float64, N}
     half_repeats :: SVector{N, Float64}
 
-    indices :: Array{Vector{SusceptibilityGridElement{N}}, N}
+    indices :: Array{Vector{Tuple{E, IndexRepeat}}, N}
     sources :: Vector{O}
     shifts :: Vector{SVector{N, Float64}}
     B0_field :: SVector{N, Float64}
 end
 
+
+const IndexNoRepeat = @NamedTuple{index::Int32}
+
 """
 Specialised version of [`SusceptibilityGrid`](@ref) for non-repeating geometries.
 """
-struct SusceptibilityGridNoRepeat{N, O, K} <: SusceptibilityGrid{N, O, K}
+struct SusceptibilityGridNoRepeat{N, O, E, K} <: SusceptibilityGrid
     inv_resolution :: SVector{N, Float64}
     rotation :: SMatrix{N, 3, Float64, K}
 
@@ -58,11 +91,11 @@ struct SusceptibilityGridNoRepeat{N, O, K} <: SusceptibilityGrid{N, O, K}
     off_resonance :: Array{Float64, N}
 
     bounding_box_indices :: BoundingBox{N}
-    indices :: Array{Vector{SusceptibilityGridElement{N}}, N}
+    indices :: Array{Vector{Tuple{E, IndexNoRepeat}}, N}
     sources :: Vector{O}
     shifts :: Vector{SVector{N, Float64}}
     B0_field :: SVector{N, Float64}
-    total_susceptibility :: Float64
+    total_susceptibility :: SVector{N, Float64}
     center_susceptibility :: SVector{N, Float64}
 end
 
@@ -108,13 +141,13 @@ function susceptibility_off_resonance(grid::SusceptibilityGrid, position::SVecto
         return field
     end
 
-    for element in grid.indices[coord_indices...]
-        if grid isa SusceptibilityGridNoRepeat || iszero(element.shift)
+    for (element, index) in grid.indices[coord_indices...]
+        if grid isa SusceptibilityGridNoRepeat || iszero(index.shift)
             shifted = normed
         else
-            shifted = normed .- grid.shifts[element.shift]
+            shifted = normed .- grid.shifts[index.shift]
         end
-        field += element_susceptibility(element, grid, shifted, inside)
+        field += element_susceptibility(element, index.index, grid, shifted, inside)
     end
     return field
 end
@@ -128,9 +161,13 @@ susceptibility_off_resonance(fixed::FixedSusceptibility{0}, position::SVector{3,
 
 """
     dipole_approximation(susceptibility, offset, distance, B0_field)
+    dipole_approximation(magnetisation, offset, distance, B0_field)
 
 Computes the shift in magnetic field due to a shift of susceptibility which is at the given offset (with given pre-computed distance) along the B0 field.
 
+The `susceptibility` is given as a scalar value and is assumed to generate a field in the z-direction.
+The `magnetisation` is a vector representing a magnetisation in an arbitrary direction (as required for anisotropic susceptibilities).
+
 Equations from [schenck96_role_magnet_suscep_magnet_reson_imagin](@cite).
 """
 function dipole_approximation(susceptiblity::Float64, offset::SVector{3, Float64}, dist::Float64, B0_field::SVector{3, Float64})
@@ -145,6 +182,24 @@ function dipole_approximation(susceptiblity::Float64, offset::SVector{3, Float64
     return susceptiblity / (4π * dist^5) * (3 * offset_B0^2 - dist^2)
 end
 
+function dipole_approximation(magnetisation::SVector{3, Float64}, offset::SVector{3, Float64}, dist::Float64, B0_field::SVector{3, Float64})
+    offset_mag = (
+        magnetisation[1] * offset[1] +
+        magnetisation[2] * offset[2] +
+        magnetisation[3] * offset[3]
+    )
+    offset_b0 = (
+        B0_field[1] * offset[1] +
+        B0_field[2] * offset[2] +
+        B0_field[3] * offset[3]
+    )
+    mag_b0 = (
+        B0_field[1] * magnetisation[1] +
+        B0_field[2] * magnetisation[2] +
+        B0_field[3] * magnetisation[3]
+    )
+    return 1 / (4π * dist^5) * (3 * offset_b0 * offset_mag - mag_b0 * dist^2)
+end
 
 average_field(x, y, nsteps) = mean([1/sqrt((d+x)^2 + y^2)^3 for d in range(-0.5, 0.5, nsteps * 2 + 1)[2:2:end-1]])
 
@@ -166,6 +221,22 @@ function dipole_approximation(susceptiblity::Float64, offset::SVector{2, Float64
     return susceptiblity / (2π * dist^4) * (2 * offset_B0^2 - dist^2 * inplane_B0_sqr)
 end
 
+function dipole_approximation(magnetisation::SVector{2, Float64}, offset::SVector{2, Float64}, dist::Float64, B0_field::SVector{2, Float64})
+    offset_mag = (
+        magnetisation[1] * offset[1] +
+        magnetisation[2] * offset[2]
+    )
+    offset_b0 = (
+        B0_field[1] * offset[1] +
+        B0_field[2] * offset[2]
+    )
+    mag_b0 = (
+        B0_field[1] * magnetisation[1] +
+        B0_field[2] * magnetisation[2]
+    )
+    return 1 / (2π * dist^4) * (2 * offset_b0 * offset_mag - mag_b0 * dist^2)
+end
+
 function dipole_approximation_repeat(susceptiblity::Float64, offset::SVector{N, Float64}, B0_field::SVector{N, Float64}, repeats::SVector{N, Float64}) where{N}
     half_repeats = repeats ./ 2
     normed = @. mod(offset + half_repeats, repeats) - half_repeats
@@ -185,7 +256,7 @@ Computes the off-resonance field contribution from a [`SuscetibilityGridElement`
 For a `position` within twice the `source.radius` of `source.position`, this will call [`dipole_approximation`](@ref).
 For any closer `position` [`single_susceptibility`](@ref) will be called on the appropriate element in `grid.sources`.
 """
-function element_susceptibility(element::SusceptibilityGridElement, grid::SusceptibilityGrid, position::AbstractVector, stuck_inside::Union{Nothing, Bool})
+function element_susceptibility(element::SusceptibilityGridElement, index::Int32, grid::SusceptibilityGrid, position::AbstractVector, stuck_inside::Union{Nothing, Bool})
     offset = position - element.position
     dist = norm(offset)
 
@@ -193,7 +264,7 @@ function element_susceptibility(element::SusceptibilityGridElement, grid::Suscep
         return dipole_approximation(element.susceptibility, offset, dist, grid.B0_field)
     else
         return single_susceptibility(
-            grid.sources[element.index],
+            grid.sources[index],
             offset,
             dist,
             stuck_inside,