cumsum(trj_l) for speed with long trajectories

src/BackProjection.jl

@@ -30,17 +30,17 @@ function calculateBackProjection(data::AbstractVector{<:AbstractArray{cT}}, trj:
     Ncoil = size(data[1], 2)
     xbp = Array{cT}(undef, img_shape..., Ncoef, Ncoil)
 
-    trj_l = [size(trj[it], 2) for it in eachindex(trj)]
-    data_temp = Vector{cT}(undef, sum(trj_l))
+    trj_idx = cumsum([size(trj[it],2) for it in eachindex(trj)])
+    data_temp = Vector{cT}(undef, trj_idx[end])
 
     img_idx = CartesianIndices(img_shape)
     verbose && println("calculating backprojection...")
     flush(stdout)
     for icoef ∈ axes(U, 2)
         t = @elapsed for icoil ∈ axes(data[1], 2)
             @simd for it in eachindex(data)
-                idx1 = sum(trj_l[1:it-1]) + 1
-                idx2 = sum(trj_l[1:it])
+                idx1 = (it == 1) ? 1 : trj_idx[it-1] + 1
+                idx2 = trj_idx[it]
                 @views data_temp[idx1:idx2] .= data[it][:, icoil] .* conj(U[it, icoef])
             end
             applyDensityCompensation!(data_temp, trj_v; density_compensation)
@@ -64,16 +64,16 @@ function calculateBackProjection(data::AbstractVector{<:AbstractMatrix{cT}}, trj
     xbp = zeros(cT, img_shape..., Ncoef)
     xtmp = Array{cT}(undef, img_shape)
 
-    trj_l = [size(trj[it], 2) for it in eachindex(trj)]
-    data_temp = Vector{cT}(undef, sum(trj_l))
+    trj_idx = cumsum([size(trj[it],2) for it in eachindex(trj)])
+    data_temp = Vector{cT}(undef, trj_idx[end])
     img_idx = CartesianIndices(img_shape)
     verbose && println("calculating backprojection...")
     flush(stdout)
     for icoef ∈ axes(U, 2)
         t = @elapsed for icoil ∈ eachindex(cmaps)
             @simd for it in eachindex(data)
-                idx1 = sum(trj_l[1:it-1]) + 1
-                idx2 = sum(trj_l[1:it])
+                idx1 = (it == 1) ? 1 : trj_idx[it-1] + 1
+                idx2 = trj_idx[it]
                 @views data_temp[idx1:idx2] .= data[it][:, icoil] .* conj(U[it, icoef])
             end
             applyDensityCompensation!(data_temp, trj_v; density_compensation)

src/NFFTNormalOpBasisFunc.jl

@@ -105,21 +105,21 @@ function calculateToeplitzKernelBasis(img_shape_os, trj::AbstractVector{<:Abstra
     λ2 = similar(λ)
     Λ  = Array{Complex{T}}(undef, Ncoeff, Ncoeff, length(kmask_indcs))
 
-    trj_l = [size(trj[it],2) for it in eachindex(trj)]
-    S  = Vector{Complex{T}}(undef, sum(trj_l))
+    trj_idx = cumsum([size(trj[it],2) for it in eachindex(trj)])
+    S  = Vector{Complex{T}}(undef, trj_idx[end])
 
     fftplan  = plan_fft(λ; flags = FFTW.MEASURE, num_threads=Threads.nthreads())
     nfftplan = plan_nfft(reduce(hcat, trj), img_shape_os; precompute = TENSOR, blocking = true, fftflags = FFTW.MEASURE, m=5, σ=2)
 
-    # Evaluating only the upper triangular matrix assumes that the PSF from the rightmost voxel to the leftmost voxel is the adjoint of the PSF in the opposite direction. 
-    # For the outmost voxel, this is not correct, but the resulting images are virtually identical in our test cases. 
-    # To avoid this error, remove the `if ic2 >= ic1` and the `Λ[ic2,ic1,it] = conj.(λ[kmask_indcs[it]])` statements at the cost of computation time.  
+    # Evaluating only the upper triangular matrix assumes that the PSF from the rightmost voxel to the leftmost voxel is the adjoint of the PSF in the opposite direction.
+    # For the outmost voxel, this is not correct, but the resulting images are virtually identical in our test cases.
+    # To avoid this error, remove the `if ic2 >= ic1` and the `Λ[ic2,ic1,it] = conj.(λ[kmask_indcs[it]])` statements at the cost of computation time.
     for ic2 ∈ axes(Λ, 2), ic1 ∈ axes(Λ, 1)
         if ic2 >= ic1 # eval. only upper triangular matrix
             t = @elapsed begin
                 @simd for it ∈ axes(U,1)
-                    idx1 = sum(trj_l[1:it-1]) + 1
-                    idx2 = sum(trj_l[1:it])
+                    idx1 = (it == 1) ? 1 : trj_idx[it-1] + 1
+                    idx2 = trj_idx[it]
                     @inbounds S[idx1:idx2] .= conj(U[it,ic1]) * U[it,ic2]
                 end
 
@@ -128,7 +128,7 @@ function calculateToeplitzKernelBasis(img_shape_os, trj::AbstractVector{<:Abstra
                 mul!(λ, fftplan, λ2)
 
                 Threads.@threads for it ∈ eachindex(kmask_indcs)
-                    @inbounds Λ[ic2,ic1,it] = conj.(λ[kmask_indcs[it]]) 
+                    @inbounds Λ[ic2,ic1,it] = conj.(λ[kmask_indcs[it]])
                     @inbounds Λ[ic1,ic2,it] =       λ[kmask_indcs[it]]
                 end
             end