Merge pull request #488 from JuliaHealth/fix-spinspan

Fix bugs related with `SpinRange` and flow

KomaMRIBase/src/motion/motionlist/Motion.jl

@@ -164,11 +164,11 @@ Base.:(≈)(m1::Motion, m2::Motion)  = (typeof(m1) == typeof(m2)) & reduce(&, [g
 """ Motion sub-group """
 function Base.getindex(m::Motion, p)
     idx, spin_range = m.spins[p]
-    return Motion(m.action[idx], m.time, spin_range)
+    return idx !== nothing ? Motion(m.action[idx], m.time, spin_range) : nothing
 end
 function Base.view(m::Motion, p)
     idx, spin_range = @view(m.spins[p])
-    return Motion(@view(m.action[idx]), m.time, spin_range)
+    return idx !== nothing ? Motion(@view(m.action[idx]), m.time, spin_range) : nothing
 end
 
 # Auxiliary functions

KomaMRIBase/src/motion/motionlist/MotionList.jl

@@ -38,14 +38,14 @@ MotionList(motions...) = length([motions]) > 0 ? MotionList([motions...]) : @err
 function Base.getindex(mv::MotionList{T}, p) where {T<:Real}
     motion_array_aux = Motion{T}[]
     for m in mv.motions
-        push!(motion_array_aux, m[p])
+        m[p] !== nothing ? push!(motion_array_aux, m[p]) : nothing
     end
     return length(motion_array_aux) > 0 ? MotionList(motion_array_aux) : NoMotion{T}()
 end
 function Base.view(mv::MotionList{T}, p) where {T<:Real}
     motion_array_aux = Motion{T}[]
     for m in mv.motions
-        push!(motion_array_aux, @view(m[p]))
+        @view(m[p]) !== nothing ? push!(motion_array_aux, @view(m[p])) : nothing
     end
     return length(motion_array_aux) > 0 ? MotionList(motion_array_aux) : NoMotion{T}()
 end

KomaMRIBase/src/motion/motionlist/SpinSpan.jl

@@ -58,12 +58,12 @@ end
 # Functions
 function Base.getindex(spins::SpinRange, p)
     idx = intersect_idx(spins.range, p)
-    return intersect_idx(p, spins.range), SpinRange(idx)
-end
-function Base.view(spins::SpinRange, p)
-    idx = intersect_idx(spins.range, p)
-    return intersect_idx(p, spins.range), SpinRange(idx)
+    l = length(idx)
+    intersect  = l >= 1 ? intersect_idx(p, spins.range) : nothing
+    spin_range = l >= 2 ? SpinRange(idx) : (l == 1 ? SpinRange(idx[1]:idx[1]) : nothing)
+    return intersect, spin_range
 end
+Base.view(spins::SpinRange, p) = spins[p]
 Base.:(==)(sr1::SpinRange, sr2::SpinRange) = sr1.range == sr2.range
 Base.length(sr::SpinRange) = length(sr.range)
 get_indexing_range(spins::SpinRange) = spins.range

KomaMRIBase/test/runtests.jl

@@ -555,14 +555,14 @@ end
 
     simplemotion = MotionList(
         Translate(0.05, 0.05, 0.0, Periodic(period=0.5, asymmetry=0.5)),
-        Rotate(0.0, 0.0, 90.0, TimeRange(t_start=0.05, t_end=0.5))
+        Rotate(0.0, 0.0, 90.0, TimeRange(t_start=0.05, t_end=0.5), SpinRange(1:3))
     )
 
     Ns = length(obj1)
     Nt = 3
     t_start = 0.0
     t_end = 1.0
-    arbitrarymotion = MotionList(Path(0.01 .* rand(Ns, Nt), 0.01 .* rand(Ns, Nt), 0.01 .* rand(Ns, Nt), TimeRange(t_start, t_end)))
+    arbitrarymotion = MotionList(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(
@@ -602,7 +602,6 @@ end
     obs1.motion = arbitrarymotion
     obs2.motion = arbitrarymotion[rng]
     @test obs1[rng] == obs2
-    # @test @view(obs1[rng]) == obs2
 
     # Test addition of phantoms
     oba = Phantom(

KomaMRICore/src/simulation/Flow.jl

@@ -30,6 +30,7 @@ function outflow_spin_reset!(
    # Get spin state range affected by the spin span
    idx = KomaMRIBase.get_indexing_range(spin_span)
    spin_state_matrix = @view(spin_state_matrix[idx, :])
+   replace_by = replace_view(replace_by, idx)
    # Obtain mask
    mask = get_mask(action.spin_reset, ts)
    # Modify spin state: reset and replace by initial value
@@ -53,6 +54,7 @@ function outflow_spin_reset!(
    # Get spin state range affected by the spin span
    idx = KomaMRIBase.get_indexing_range(spin_span)
    M = @view(M[idx])
+   replace_by = replace_view(replace_by, idx)
    # Obtain mask
    mask = get_mask(action.spin_reset, ts)
    mask = @view(mask[:, end])
@@ -72,6 +74,13 @@ function init_time(t, seq_t, add_t0)
    return t
 end
 
+function replace_view(replace_by::AbstractArray, idx)
+   return @view(replace_by[idx])
+end
+function replace_view(replace_by, idx)
+   return replace_by
+end
+
 function get_mask(spin_reset, t::Real)
    itp  = KomaMRIBase.interpolate(spin_reset, KomaMRIBase.Gridded(KomaMRIBase.Constant{KomaMRIBase.Previous}()), Val(size(spin_reset, 1)), t)
    return KomaMRIBase.resample(itp, t)

KomaMRIPlots/src/ui/DisplayFunctions.jl

@@ -1062,16 +1062,8 @@ function plot_phantom_map(
     end
 
     function decimate_uniform_phantom(obj, num_points::Int)
-        dimx, dimy, dimz = KomaMRIBase.get_dims(obj)
-        ss = Int(ceil((length(obj) / num_points)^(1 / sum(KomaMRIBase.get_dims(obj)))))
-        ssx = dimx ? ss : 1
-        ssy = dimy ? ss : 1
-        ssz = dimz ? ss : 1
-        ix = sortperm(obj.x)[1:ssx:end]
-        iy = sortperm(obj.y)[1:ssy:end]
-        iz = sortperm(obj.z)[1:ssz:end]
-        idx = intersect(ix, iy, iz)
-        return obj[idx]
+        ss = Int(ceil(length(obj) / num_points))
+        return obj[1:ss:end]
     end
 
     if length(obj) > max_spins
@@ -1360,16 +1352,8 @@ function plot_phantom_map(
     kwargs...,
 )
     function decimate_uniform_phantom(obj, num_points::Int)
-        dimx, dimy, dimz = KomaMRIBase.get_dims(obj)
-        ss = Int(ceil((length(obj) / num_points)^(1 / sum(KomaMRIBase.get_dims(obj)))))
-        ssx = dimx ? ss : 1
-        ssy = dimy ? ss : 1
-        ssz = dimz ? ss : 1
-        ix = sortperm(obj.x)[1:ssx:end]
-        iy = sortperm(obj.y)[1:ssy:end]
-        iz = sortperm(obj.z)[1:ssz:end]
-        idx = intersect(ix, iy, iz)
-        return obj[idx]
+        ss = Int(ceil(length(obj) / num_points))
+        return obj[1:ss:end]
     end
 
     if length(obj) > max_spins
@@ -1448,7 +1432,7 @@ function plot_phantom_map(
         l.width = width
     end
 	if view_2d
-        h = scatter( 
+        h = scattergl( 
             x=obj.x*1e2,
             y=obj.y*1e2,
             mode="markers",