Merge pull request #520 from JuliaHealth/easy-motion

Solve bug in `SpinRanges` when indexing a phantom with a BitVector
JuliaHealth · Dec 10, 2024 · 6d4deb8 · 6d4deb8
2 parents 96fb134 + 00eee39
commit 6d4deb8
Show file tree

Hide file tree

Showing 3 changed files with 14 additions and 2 deletions.
diff --git a/KomaMRIBase/src/motion/Motion.jl b/KomaMRIBase/src/motion/Motion.jl
@@ -191,7 +191,7 @@ end
 function get_spin_coords(
     m::Motion{T}, x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T}, t
 ) where {T<:Real}
-    ux, uy, uz = x .+ 0*t, y .+ 0*t, z .+ 0*t # Buffers for displacements
+    ux, uy, uz = x .* (0*t), y .* (0*t), z .* (0*t) # Buffers for displacements
     t_unit = unit_time(t, m.time)
     idx = get_indexing_range(m.spins)
     displacement_x!(@view(ux[idx, :]), m.action, @view(x[idx]), @view(y[idx]), @view(z[idx]), t_unit)

diff --git a/KomaMRIBase/src/motion/SpinSpan.jl b/KomaMRIBase/src/motion/SpinSpan.jl
@@ -69,3 +69,5 @@ Base.length(sr::SpinRange) = length(sr.range)
 get_indexing_range(spins::SpinRange) = spins.range
 expand(sr::SpinRange, Ns::Int) = sr
 intersect_idx(a, b) = findall(x -> x in a, b)
+intersect_idx(a, b::BitVector) = findall(x -> x in a, findall(x->x==true, b))
+intersect_idx(a::BitVector, b) = findall(x -> x in findall(x->x==true, a), b)
diff --git a/KomaMRIBase/test/runtests.jl b/KomaMRIBase/test/runtests.jl
@@ -400,7 +400,6 @@ end
     t_end = 1.0
     arbitrarymotion = Path(0.01 .* rand(Ns, Nt), 0.01 .* rand(Ns, Nt), 0.01 .* rand(Ns, Nt), TimeRange(t_start, t_end), SpinRange(2:2:4))
 
-    # Test phantom subset
     obs1 = Phantom(
         name,
         x,
@@ -416,6 +415,8 @@ end
         Dθ,
         simplemotion
     )
+
+    # Test phantom subset (simple range)
     rng = 1:2:5
     obs2 = Phantom(
         name,
@@ -457,6 +458,15 @@ end
     )
     @test obs1 + obs2 == oba
 
+    # Test phantom subset (BitVector range)
+    obs3 = copy(obs1)
+    obs4 = copy(obs1)
+    rng = obs3.x .> 0
+    obs3.motion = Translate(5e-4, 6e-4, 7e-4, TimeRange(0.0, 1.0), SpinRange(rng))
+    obs4.motion = Translate(5e-4, 6e-4, 7e-4, TimeRange(0.0, 1.0), SpinRange(1:length(obs4)))
+    @test obs3[rng] == obs4[rng]
+    @test obs3[rng].motion == obs4.motion[rng]
+
     # Test scalar multiplication of a phantom
     c = 7
     obc = Phantom(name=name, x=x, y=y, z=z, ρ=c*ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw, Dλ1=Dλ1, Dλ2=Dλ2, Dθ=Dθ)