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

Optimize ArbitraryMotion (continuation) #408

Merged
merged 115 commits into from
Jun 27, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
115 commits
Select commit Hold shift + click to select a range
5120dee
First commit
pvillacorta May 31, 2024
87adbe0
Fix bugs and pass tests
pvillacorta May 31, 2024
3506d03
Merge branch 'master' into optimize-arbitrary
pvillacorta Jun 2, 2024
60e5743
Merge branch 'master' into optimize-arbitrary
pvillacorta Jun 4, 2024
9b2a9c2
Simplify adapt functions
pvillacorta Jun 4, 2024
1ee2bb8
ExplicitArbitraryMotion, `sort_motions!` -> `initialize_motion`
pvillacorta Jun 12, 2024
b70ca80
Minor change in `plot_image`
pvillacorta Jun 12, 2024
a387d22
Merge branch 'master' of https://github.com/pvillacorta/KomaMRI.jl
pvillacorta Jun 12, 2024
97eb2a9
Merge branch 'master' into optimize-arbitrary-parallel
pvillacorta Jun 12, 2024
7c62afa
Delete ExplicitArbitraryMotion, custom `GriddedInterpolation` function
pvillacorta Jun 14, 2024
83e0527
Merge branch 'JuliaHealth:master' into master
pvillacorta Jun 14, 2024
d282f9f
Improve adapt a bit also for SimpleMotion. Use functor
pvillacorta Jun 14, 2024
b82fa19
Remove comments
pvillacorta Jun 14, 2024
352ebbb
Simplify getindex for ArbitraryMotion
pvillacorta Jun 14, 2024
c8e5beb
Merge branch 'master' of https://github.com/pvillacorta/KomaMRI.jl
pvillacorta Jun 15, 2024
8414916
Merge branch 'master' into optimize-arbitrary-view
pvillacorta Jun 15, 2024
49bf6a8
Solve typo
pvillacorta Jun 15, 2024
4d41071
Delete adapt for ArbitraryMotion
pvillacorta Jun 15, 2024
c8a754b
Improve code coverage
pvillacorta Jun 15, 2024
f044ede
Not exporting unnecessary things
pvillacorta Jun 17, 2024
ed0ad5e
Add testset macro
pvillacorta Jun 17, 2024
9ef5d81
Add specific identifiers
pvillacorta Jun 17, 2024
4010c2d
Add path for avoiding errors
pvillacorta Jun 17, 2024
596379e
Add aux functions to tests
pvillacorta Jun 17, 2024
c9f0980
Fix typo
pvillacorta Jun 17, 2024
368cec2
Remove `collect`
pvillacorta Jun 17, 2024
5050c46
Simplify range
pvillacorta Jun 17, 2024
9ab9abb
Fix error from previous commit
pvillacorta Jun 17, 2024
b8a2b83
Interpolation creation with ranges
pvillacorta Jun 17, 2024
b4b2307
Evaluate using UnitRange for spin id
pvillacorta Jun 17, 2024
1733ba6
Change functions names
pvillacorta Jun 17, 2024
ed7afd2
Fix test for ArbitraryMotion
pvillacorta Jun 17, 2024
cfeddb2
Simplify vcat for ArbitraryMotion
pvillacorta Jun 17, 2024
47ccfc0
Recover `return`
pvillacorta Jun 17, 2024
e10474b
Simplify ArbitraryMotion read and write
pvillacorta Jun 17, 2024
58756f9
Remove K from tests
pvillacorta Jun 17, 2024
8dfaf1d
Improve `unit_time` functions
pvillacorta Jun 18, 2024
6134911
Continuation of last commit
pvillacorta Jun 18, 2024
c50da52
Add tests when t_start= t_end
pvillacorta Jun 18, 2024
a58485b
Correct typos in SimpleMotionTypes docstrings
pvillacorta Jun 18, 2024
115815f
Delete `maxlog=1` and use `scattergl`
pvillacorta Jun 18, 2024
eeaf970
Merge branch 'master' into optimize-arbitrary-view
pvillacorta Jun 18, 2024
2bce2e4
Correct `view` bug
pvillacorta Jun 18, 2024
872d16a
Dispatch on Val(Ns)
pvillacorta Jun 18, 2024
7172839
one -> oneunit
pvillacorta Jun 18, 2024
5bc95ef
`fields` -> `phantom_fields`
pvillacorta Jun 18, 2024
57149eb
Merge branch 'master' into fix-plot-phantom
pvillacorta Jun 18, 2024
e040a0f
Resolve functor for SimpleMotion
pvillacorta Jun 19, 2024
079a84c
Simplify methods of `proces_times`
pvillacorta Jun 19, 2024
c9a4dbb
Custom `range` function for keeping type matching
pvillacorta Jun 19, 2024
3a0862e
Remove TupleTools dependency, initialize_motion -> sort_motion (only …
pvillacorta Jun 19, 2024
f1186c3
Solve method redefinition problem
pvillacorta Jun 19, 2024
37cdfae
Merge branch 'master' into fix-plot-phantom
cncastillo Jun 19, 2024
4df277b
Fix rangeT, but use LinRange
pvillacorta Jun 19, 2024
ba92e75
Solve bug
pvillacorta Jun 19, 2024
10a4aa3
Change process_times in plot_phantom_map
pvillacorta Jun 19, 2024
3d45b55
Solve bug
pvillacorta Jun 19, 2024
0310cab
Merge branch 'master' into optimize-arbitrary-view
pvillacorta Jun 19, 2024
1315be2
Reduce allocations in `get_spin_coords` for SimpleMotion
pvillacorta Jun 20, 2024
de1da32
Remove `rangeT`
pvillacorta Jun 21, 2024
513ccf9
Merge branch 'master' into optimize-arbitrary-view
pvillacorta Jun 21, 2024
4107c1c
Remove :skipci for Simple and Arbitrary Motion
pvillacorta Jun 21, 2024
95a9fca
Try vectors instead of ranges
pvillacorta Jun 21, 2024
58df3e6
Renges with correct eltype
pvillacorta Jun 21, 2024
ba4b468
More testing...
pvillacorta Jun 21, 2024
3e16df7
More...
pvillacorta Jun 21, 2024
87c7b4f
t in GPU, id in CPU
pvillacorta Jun 21, 2024
75918fe
Try with vectors
pvillacorta Jun 21, 2024
fbac368
Fix bug
pvillacorta Jun 21, 2024
e7b5861
Another bug
pvillacorta Jun 21, 2024
fa97756
Let's try with `LinRange`
pvillacorta Jun 21, 2024
9288e04
`copyto!`
pvillacorta Jun 21, 2024
fe98989
Debug
pvillacorta Jun 21, 2024
d12f81d
Revert last change
pvillacorta Jun 21, 2024
db425e6
Merge branch 'master' into optimize-arbitrary-view
cncastillo Jun 21, 2024
409a719
Add :skipci to everything except simple and arbitrary
pvillacorta Jun 22, 2024
8085239
Debugging prints
pvillacorta Jun 22, 2024
ce91eb8
Remove @supress temporarily
pvillacorta Jun 22, 2024
81160e0
Try with grid for interpolate solving
pvillacorta Jun 22, 2024
eadf2fa
Try Interpolate directly in test
pvillacorta Jun 22, 2024
7112b64
More debugging
pvillacorta Jun 22, 2024
aad423a
Solve type conflict
pvillacorta Jun 22, 2024
781c051
Change conversion order
pvillacorta Jun 22, 2024
d6f5da7
Try with 2D since 1D works
pvillacorta Jun 22, 2024
eac3e87
Try with using `gpu` function
pvillacorta Jun 22, 2024
f1e3dee
Check if the problem commes from `NoInterp`
pvillacorta Jun 22, 2024
723b418
Try with ranges again
pvillacorta Jun 22, 2024
0348a48
Back again to vectors
pvillacorta Jun 22, 2024
206fc37
Final test with vectors and Gridded in both dimensions
pvillacorta Jun 22, 2024
f7683a5
Remove comments and recover :skipci for all testitems
pvillacorta Jun 22, 2024
c7f4a3b
Recover `@suppress`
pvillacorta Jun 22, 2024
10b814e
Merge branch 'master' into optimize-arbitrary-view
pvillacorta Jun 24, 2024
0f4a85b
Merge branch 'master' into fix-plot-phantom
pvillacorta Jun 24, 2024
431c164
Merge branch 'fix-plot-phantom' into optimize-arbitrary-view
pvillacorta Jun 24, 2024
cfa1165
Keep `scatter` for now
pvillacorta Jun 24, 2024
cbb1710
Test the exact desired situation (but with Gridded in both dimensions)
pvillacorta Jun 26, 2024
f082e99
Solve bug, retry...
pvillacorta Jun 26, 2024
f336659
Check if it passes now
pvillacorta Jun 26, 2024
ade4138
Try with the simplest situation to reproduce
pvillacorta Jun 26, 2024
b9bfbfe
Try again
pvillacorta Jun 26, 2024
88e7451
Keep type matching
pvillacorta Jun 26, 2024
2ea2d98
Try vectors with the simplest reproducible situation
pvillacorta Jun 26, 2024
f06cf21
Try ranges with the simplest reproducible situation
pvillacorta Jun 26, 2024
f9908d7
Try again but forcing range step to be float32
pvillacorta Jun 26, 2024
529a250
Back to vectors, mergeable
pvillacorta Jun 26, 2024
250f4f2
`:skipci` also for ArbitraryMotion (for now)
pvillacorta Jun 26, 2024
1e8d8d7
Merge branch 'master' into optimize-arbitrary-view
pvillacorta Jun 26, 2024
ba82913
Update KomaMRICore/test/runtests.jl
pvillacorta Jun 27, 2024
5c6fe49
_utils.jl -> sorting.jl
pvillacorta Jun 27, 2024
639d9e4
Use `.=` instead of `copy!` for test in Metal
pvillacorta Jun 27, 2024
df33c54
All `skipci` except motions
pvillacorta Jun 27, 2024
4eaac2a
Skipci for all except arbitrarymotion
pvillacorta Jun 27, 2024
c26bd4b
Back to `copyto!`
pvillacorta Jun 27, 2024
c5a49ae
Independent test sets for every simple motion
pvillacorta Jun 27, 2024
5dcc8b8
Leave `skipci`s as they were
pvillacorta Jun 27, 2024
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/KomaMRIBase.jl
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@ include("datatypes/Phantom.jl")
include("datatypes/simulation/DiscreteSequence.jl")
include("timing/TimeStepCalculation.jl")
include("timing/TrapezoidalIntegration.jl")
include("timing/UnitTime.jl")

# Main
export γ # gyro-magnetic ratio [Hz/T]
Expand All @@ -51,8 +52,7 @@ export NoMotion, SimpleMotion, ArbitraryMotion
export SimpleMotionType
export Translation, Rotation, HeartBeat
export PeriodicTranslation, PeriodicRotation, PeriodicHeartBeat
export get_spin_coords, sort_motions!
export LinearInterpolator
export get_spin_coords
# Secondary
export get_kspace, rotx, roty, rotz
# Additionals
Expand Down
21 changes: 10 additions & 11 deletions KomaMRIBase/src/datatypes/Phantom.jl
Original file line number Diff line number Diff line change
Expand Up @@ -69,10 +69,7 @@ Base.getindex(x::Phantom, i::Integer) = x[i:i]
"""Compare two phantoms"""
Base.:(==)(obj1::Phantom, obj2::Phantom) = reduce(
&,
[
getfield(obj1, field) == getfield(obj2, field) for
field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))
],
[getfield(obj1, field) == getfield(obj2, field) for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))],
)
Base.:(≈)(obj1::Phantom, obj2::Phantom) = reduce(&, [getfield(obj1, field) ≈ getfield(obj2, field) for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))])
Base.:(==)(m1::MotionModel, m2::MotionModel) = false
Expand All @@ -88,7 +85,13 @@ Base.getindex(obj::Phantom, p::Union{AbstractRange,AbstractVector,Colon}) = begi
end

"""Separate object spins in a sub-group (lightweigth)."""
Base.view(obj::Phantom, p::Union{AbstractRange,AbstractVector,Colon}) = @views obj[p]
Base.view(obj::Phantom, p::Union{AbstractRange,AbstractVector,Colon}) = begin
fields = []
for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))
push!(fields, (field, @view(getfield(obj, field)[p])))
end
return Phantom(; name=obj.name, fields...)
end

"""Addition of phantoms"""
+(obj1::Phantom, obj2::Phantom) = begin
Expand Down Expand Up @@ -117,10 +120,6 @@ function get_dims(obj::Phantom)
return dims
end

function sort_motions!(motion::MotionModel)
return nothing
end

"""
obj = heart_phantom(...)

Expand Down Expand Up @@ -178,7 +177,7 @@ function heart_phantom(
Dλ1=Dλ1[ρ .!= 0],
Dλ2=Dλ2[ρ .!= 0],
Dθ=Dθ[ρ .!= 0],
motion=SimpleMotion([
motion=SimpleMotion(
PeriodicHeartBeat(;
period=period,
asymmetry=asymmetry,
Expand All @@ -189,7 +188,7 @@ function heart_phantom(
PeriodicRotation(;
period=period, asymmetry=asymmetry, yaw=rotation_angle, pitch=0.0, roll=0.0
),
]),
),
)
return phantom
end
Expand Down
153 changes: 70 additions & 83 deletions KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,13 +2,24 @@
# Interpolator{T,Degree,ETPType},
# Degree = Linear,Cubic....
# ETPType = Periodic, Flat...
const LinearInterpolator = Interpolations.Extrapolation{
T,
1,
Interpolations.GriddedInterpolation{T,1,V,Gridded{Linear{Throw{OnGrid}}},Tuple{V}},
Gridded{Linear{Throw{OnGrid}}},
Periodic{Nothing},
} where {T<:Real,V<:AbstractVector{T}}

const Interpolator1D = Interpolations.GriddedInterpolation{
T,1,V,Itp,K
} where {
T<:Real,
V<:AbstractArray{T},
Itp<:Interpolations.Gridded{Linear{Throw{OnGrid}}},
K<:Tuple{AbstractVector{T}},
}

const Interpolator2D = Interpolations.GriddedInterpolation{
T,2,V,Itp,K
} where {
T<:Real,
V<:AbstractArray{T},
Itp<:Interpolations.Gridded{Linear{Throw{OnGrid}}},
K<:Tuple{AbstractVector{T}, AbstractVector{T}},
}

"""
motion = ArbitraryMotion(period_durations, dx, dy, dz)
Expand Down Expand Up @@ -43,106 +54,82 @@
)
```
"""
struct ArbitraryMotion{T<:Real,V<:AbstractVector{T}} <: MotionModel{T}
period_durations::Vector{T}
dx::Array{T,2}
dy::Array{T,2}
dz::Array{T,2}
ux::Vector{LinearInterpolator{T,V}}
uy::Vector{LinearInterpolator{T,V}}
uz::Vector{LinearInterpolator{T,V}}
end

function ArbitraryMotion(
period_durations::AbstractVector{T},
dx::AbstractArray{T,2},
dy::AbstractArray{T,2},
dz::AbstractArray{T,2},
) where {T<:Real}
@warn "Note that ArbitraryMotion is under development so it is not optimized so far" maxlog = 1
Ns = size(dx)[1]
num_pieces = size(dx)[2] + 1
limits = times(period_durations, num_pieces)

#! format: off
Δ = zeros(Ns,length(limits),4)
Δ[:,:,1] = hcat(repeat(hcat(zeros(Ns,1),dx),1,length(period_durations)),zeros(Ns,1))
Δ[:,:,2] = hcat(repeat(hcat(zeros(Ns,1),dy),1,length(period_durations)),zeros(Ns,1))
Δ[:,:,3] = hcat(repeat(hcat(zeros(Ns,1),dz),1,length(period_durations)),zeros(Ns,1))

etpx = [extrapolate(interpolate((limits,), Δ[i,:,1], Gridded(Linear())), Periodic()) for i in 1:Ns]
etpy = [extrapolate(interpolate((limits,), Δ[i,:,2], Gridded(Linear())), Periodic()) for i in 1:Ns]
etpz = [extrapolate(interpolate((limits,), Δ[i,:,3], Gridded(Linear())), Periodic()) for i in 1:Ns]
#! format: on

return ArbitraryMotion(period_durations, dx, dy, dz, etpx, etpy, etpz)
struct ArbitraryMotion{T} <: MotionModel{T}
t_start::T
t_end::T
dx::AbstractArray{T}
dy::AbstractArray{T}
dz::AbstractArray{T}
end

function Base.getindex(
motion::ArbitraryMotion, p::Union{AbstractRange,AbstractVector,Colon}
)
fields = []
for field in fieldnames(ArbitraryMotion)
if field in (:dx, :dy, :dz)
push!(fields, getfield(motion, field)[p, :])
elseif field in (:ux, :uy, :uz)
push!(fields, getfield(motion, field)[p])
else
push!(fields, getfield(motion, field))
end
end
return ArbitraryMotion(fields...)
return ArbitraryMotion(motion.t_start, motion.t_end, motion.dx[p,:], motion.dy[p,:], motion.dz[p,:])
pvillacorta marked this conversation as resolved.
Show resolved Hide resolved
end
function Base.view(

Check warning on line 70 in KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl

View check run for this annotation

Codecov / codecov/patch

KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl#L70

Added line #L70 was not covered by tests
motion::ArbitraryMotion, p::Union{AbstractRange,AbstractVector,Colon}
)
return ArbitraryMotion(motion.t_start, motion.t_end, @view(motion.dx[p,:]), @view(motion.dy[p,:]), @view(motion.dz[p,:]))

Check warning on line 73 in KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl

View check run for this annotation

Codecov / codecov/patch

KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl#L73

Added line #L73 was not covered by tests
end

Base.:(==)(m1::ArbitraryMotion, m2::ArbitraryMotion) = reduce(&, [getfield(m1, field) == getfield(m2, field) for field in fieldnames(ArbitraryMotion)])
Base.:(≈)(m1::ArbitraryMotion, m2::ArbitraryMotion) = reduce(&, [getfield(m1, field) ≈ getfield(m2, field) for field in fieldnames(ArbitraryMotion)])

function Base.vcat(m1::ArbitraryMotion, m2::ArbitraryMotion)
fields = []
@assert m1.period_durations == m2.period_durations "period_durations of both ArbitraryMotions must be the same"
for field in
Iterators.filter(x -> !(x == :period_durations), fieldnames(ArbitraryMotion))
push!(fields, [getfield(m1, field); getfield(m2, field)])
end
return ArbitraryMotion(m1.period_durations, fields...)
@assert (m1.t_start == m2.t_start) && (m1.t_end == m2.t_end) "Starting and ending times must be the same"
return ArbitraryMotion(m1.t_start, m1.t_end, [m1.dx; m2.dx], [m1.dy; m2.dy], [m1.dz; m2.dz])
end

"""
limits = times(obj.motion)
"""
function times(motion::ArbitraryMotion)
period_durations = motion.period_durations
num_pieces = size(motion.dx)[2] + 1
return times(period_durations, num_pieces)
return range(motion.t_start, motion.t_end, length=size(motion.dx, 2))
end

function GriddedInterpolation(nodes, A, ITP)
return Interpolations.GriddedInterpolation{eltype(A), length(nodes), typeof(A), typeof(ITP), typeof(nodes)}(nodes, A, ITP)
end

function interpolate(motion::ArbitraryMotion{T}, Ns::Val{1}) where {T<:Real}
_, Nt = size(motion.dx)
t = similar(motion.dx, Nt); copyto!(t, collect(range(zero(T), oneunit(T), Nt)))
cncastillo marked this conversation as resolved.
Show resolved Hide resolved
itpx = GriddedInterpolation((t, ), motion.dx[:], Gridded(Linear()))
itpy = GriddedInterpolation((t, ), motion.dy[:], Gridded(Linear()))
itpz = GriddedInterpolation((t, ), motion.dz[:], Gridded(Linear()))
return itpx, itpy, itpz
end

function interpolate(motion::ArbitraryMotion{T}, Ns::Val) where {T<:Real}
Ns, Nt = size(motion.dx)
id = similar(motion.dx, Ns); copyto!(id, collect(range(oneunit(T), T(Ns), Ns)))
t = similar(motion.dx, Nt); copyto!(t, collect(range(zero(T), oneunit(T), Nt)))
Comment on lines +106 to +107
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
id = similar(motion.dx, Ns); copyto!(id, collect(range(oneunit(T), T(Ns), Ns)))
t = similar(motion.dx, Nt); copyto!(t, collect(range(zero(T), oneunit(T), Nt)))
id = similar(motion.dx, Ns); id .= range(oneunit(T), T(Ns), Ns)
t = similar(motion.dx, Nt); t .= range(zero(T), oneunit(T), Nt)

Is id .= 1:Ns not enough?
Are the copyto!(collect( necessary? They are a little bit ugly.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like metal and oneAPI don't have this operation defined for ranges or something like that.
When assigning wit =. and a range, the following error appears (see this):

ERROR: LoadError: InvalidIRError: compiling MethodInstance for (::Metal.var"#broadcast_linear#202")(::Metal.MtlDeviceVector{Float32, 1}, ::Base.Broadcast.Broadcasted{Metal.MtlArrayStyle{1, Metal.MTL.Private}, Tuple{Base.OneTo{Int64}}, typeof(identity), Tuple{Base.Broadcast.Extruded{StepRangeLen{Float32, Float64, Float64, Int64}, Tuple{Bool}, Tuple{Int64}}}}) resulted in invalid LLVM IR

From here, it looks like =. is defined for assigning scalars, but I think it's not the case for ranges.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I think submitting a few issues to Metal/OneAPI would be good. I don't know what the MWE is, @rkierulf says these work

julia> using Metal

julia> A = ones(Float32, 4)
4-element Vector{Float32}:
 1.0
 1.0
 1.0
 1.0

julia> B = convert(MtlArray,A)
4-element MtlVector{Float32, Private}:
 1.0
 1.0
 1.0
 1.0

julia> B .= range(1,4) # Error 1: id .= range(oneunit(T), T(Ns), Ns)
4-element MtlVector{Float32, Private}:
 1.0
 2.0
 3.0
 4.0

julia> B .+ B' # Error 2: xt, yt, zt = x .+ 0*t, y .+ 0*t, z .+ 0*t
4×4 MtlMatrix{Float32, Private}:
 2.0  3.0  4.0  5.0
 3.0  4.0  5.0  6.0
 4.0  5.0  6.0  7.0
 5.0  6.0  7.0  8.0
  • Error 1: id .= range(oneunit(T), T(Ns), Ns)
  • Error 2: xt, yt, zt = x .+ 0*t, y .+ 0*t, z .+ 0*t

itpx = GriddedInterpolation((id, t), motion.dx, Gridded(Linear()))
itpy = GriddedInterpolation((id, t), motion.dy, Gridded(Linear()))
itpz = GriddedInterpolation((id, t), motion.dz, Gridded(Linear()))
return itpx, itpy, itpz
end

function resample(itpx::Interpolator1D{T}, itpy::Interpolator1D{T}, itpz::Interpolator1D{T}, t::AbstractArray{T}) where {T<:Real}
return itpx.(t), itpy.(t), itpz.(t)
end

function times(period_durations::AbstractVector, num_pieces::Int)
# Pre-allocating memory
limits = zeros(eltype(period_durations), num_pieces * length(period_durations) + 1)

idx = 1
for i in 1:length(period_durations)
segment_increment = period_durations[i] / num_pieces
cumulative_sum = limits[idx] # Start from the last computed value in limits
for j in 1:num_pieces
cumulative_sum += segment_increment
limits[idx + 1] = cumulative_sum
idx += 1
end
end
return limits
function resample(itpx::Interpolator2D{T}, itpy::Interpolator2D{T}, itpz::Interpolator2D{T}, t::AbstractArray{T}) where {T<:Real}
Ns = size(itpx.coefs, 1)
id = similar(itpx.coefs, Ns)
copyto!(id, collect(range(oneunit(T), T(Ns), Ns)))
cncastillo marked this conversation as resolved.
Show resolved Hide resolved
return itpx.(id, t), itpy.(id, t), itpz.(id, t)
end

# TODO: Calculate interpolation functions "on the fly"
function get_spin_coords(
motion::ArbitraryMotion{T},
x::AbstractVector{T},
y::AbstractVector{T},
z::AbstractVector{T},
t::AbstractArray{T},
t::AbstractArray{T}
) where {T<:Real}
xt = x .+ reduce(vcat, [etp.(t) for etp in motion.ux])
yt = y .+ reduce(vcat, [etp.(t) for etp in motion.uy])
zt = z .+ reduce(vcat, [etp.(t) for etp in motion.uz])
return xt, yt, zt
end
motion_functions = interpolate(motion, Val(size(x,1)))
ux, uy, uz = resample(motion_functions..., unit_time(t, motion.t_start, motion.t_end))
return x .+ ux, y .+ uy, z .+ uz
end
3 changes: 2 additions & 1 deletion KomaMRIBase/src/datatypes/phantom/motion/NoMotion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ x = x
struct NoMotion{T<:Real} <: MotionModel{T} end

Base.getindex(motion::NoMotion, p::Union{AbstractRange,AbstractVector,Colon}) = motion
Base.view(motion::NoMotion, p::Union{AbstractRange,AbstractVector,Colon}) = motion

Base.:(==)(m1::NoMotion, m2::NoMotion) = true
Base.:(≈)(m1::NoMotion, m2::NoMotion) = true
Expand All @@ -18,7 +19,7 @@ function get_spin_coords(
x::AbstractVector{T},
y::AbstractVector{T},
z::AbstractVector{T},
t::AbstractArray{T},
t::AbstractArray{T}
) where {T<:Real}
return x, y, z
end
Expand Down
Loading