Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix bugs related with SpinRange and flow #488

Merged
merged 5 commits into from
Sep 27, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions KomaMRIBase/src/motion/motionlist/Motion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions KomaMRIBase/src/motion/motionlist/MotionList.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
10 changes: 5 additions & 5 deletions KomaMRIBase/src/motion/motionlist/SpinSpan.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
5 changes: 2 additions & 3 deletions KomaMRIBase/test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -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(
Expand Down
9 changes: 9 additions & 0 deletions KomaMRICore/src/simulation/Flow.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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])
Expand All @@ -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)
Expand Down
26 changes: 5 additions & 21 deletions KomaMRIPlots/src/ui/DisplayFunctions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1062,16 +1062,8 @@
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
Expand Down Expand Up @@ -1360,16 +1352,8 @@
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]

Check warning on line 1356 in KomaMRIPlots/src/ui/DisplayFunctions.jl

View check run for this annotation

Codecov / codecov/patch

KomaMRIPlots/src/ui/DisplayFunctions.jl#L1355-L1356

Added lines #L1355 - L1356 were not covered by tests
end

if length(obj) > max_spins
Expand Down Expand Up @@ -1448,7 +1432,7 @@
l.width = width
end
if view_2d
h = scatter(
h = scattergl(
x=obj.x*1e2,
y=obj.y*1e2,
mode="markers",
Expand Down
Binary file removed examples/2.phantoms/artery.phantom
Binary file not shown.
Loading