From b8066a6730a0f4d1b51c4c6924b6d8ba06582468 Mon Sep 17 00:00:00 2001 From: stijn Date: Tue, 16 Apr 2024 12:18:08 +0200 Subject: [PATCH] Make jacobian functions accept/returns signal as 3D instead 2D array --- CompasToolkit.jl/examples/jacobian.jl | 2 +- CompasToolkit.jl/examples/jacobian_hermitian.jl | 6 ++++-- CompasToolkit.jl/src/CompasToolkit.jl | 10 +++++----- julia-bindings/src/bindings.cpp | 8 ++++---- src/jacobian/hermitian.cu | 14 +++++++------- src/jacobian/product.h | 4 ++-- src/jacobian/regular.cu | 8 ++++---- src/simulate/phase_encoding.cu | 17 +++++++++-------- src/simulate/phase_encoding.h | 7 +++---- src/simulate/phase_encoding_kernels.cuh | 4 ++-- 10 files changed, 41 insertions(+), 39 deletions(-) diff --git a/CompasToolkit.jl/examples/jacobian.jl b/CompasToolkit.jl/examples/jacobian.jl index ddba662..6ba4df1 100644 --- a/CompasToolkit.jl/examples/jacobian.jl +++ b/CompasToolkit.jl/examples/jacobian.jl @@ -137,4 +137,4 @@ Jv = CompasToolkit.compute_jacobian( v ) -print_equals_check(Jv_ref, transpose(collect(Jv))) \ No newline at end of file +print_equals_check(Jv_ref, transpose(reshape(collect(Jv), :, ncoils))) \ No newline at end of file diff --git a/CompasToolkit.jl/examples/jacobian_hermitian.jl b/CompasToolkit.jl/examples/jacobian_hermitian.jl index b11914e..04cf0d1 100644 --- a/CompasToolkit.jl/examples/jacobian_hermitian.jl +++ b/CompasToolkit.jl/examples/jacobian_hermitian.jl @@ -135,8 +135,10 @@ echos = generate_echos(N, pssfp_ref) ∂echos = generate_delta_echos(N, pssfp_ref) Random.seed!(1337) -v = rand(ComplexF32, nsamples(trajectory_ref), ncoils) -v_ref = map(SVector{ncoils}, eachcol(v)...) +v = rand(ComplexF32, trajectory_ref.nsamplesperreadout, trajectory_ref.nreadouts, ncoils) + +v_ref = reshape(nsamples, nsamples(trajectory_ref), ncoils) +v_ref = map(SVector{ncoils}, eachcol(v_ref)...) Jᴴv_ref = compute_Jᴴv(gpu(echos), gpu(∂echos), gpu(parameters_ref), gpu(coil_sensitivities_ref), gpu(trajectory_ref), gpu(v_ref)) Jᴴv_ref = reduce(hcat, collect(Jᴴv_ref)) # Vector{Svector} -> Matrix diff --git a/CompasToolkit.jl/src/CompasToolkit.jl b/CompasToolkit.jl/src/CompasToolkit.jl index c28b44d..dcf741a 100644 --- a/CompasToolkit.jl/src/CompasToolkit.jl +++ b/CompasToolkit.jl/src/CompasToolkit.jl @@ -546,7 +546,7 @@ function compute_jacobian( trajectory::Trajectory, coils::AbstractMatrix, v::AbstractMatrix -) +)::CompasArray{ComplexF32, 3} context = get_context() ncoils = size(coils, 2) nreadouts::Int64 = trajectory.nreadouts @@ -575,7 +575,7 @@ function compute_jacobian( pointer(v)::Ptr{Cvoid} )::Ptr{Cvoid} - return CompasArray{ComplexF32, 2}(context, Jv_ptr, (ncoils, nreadouts * samples_per_readout)) + return CompasArray{ComplexF32, 3}(context, Jv_ptr, (ncoils, nreadouts, samples_per_readout)) end function compute_jacobian_hermitian( @@ -584,8 +584,8 @@ function compute_jacobian_hermitian( parameters::TissueParameters, trajectory::Trajectory, coils::AbstractMatrix, - v::AbstractMatrix -) + v::AbstractArray{<:Any,3} +)::CompasArray{ComplexF32, 2} context = get_context() ncoils = size(coils, 2) nreadouts::Int64 = trajectory.nreadouts @@ -596,7 +596,7 @@ function compute_jacobian_hermitian( 𝜕echos_T1 = convert_array(𝜕echos.T1, ComplexF32, nvoxels, nreadouts) 𝜕echos_T2 = convert_array(𝜕echos.T2, ComplexF32, nvoxels, nreadouts) coils = convert_array(coils, Float32, nvoxels, ncoils) - v = convert_array(v, ComplexF32, nreadouts * samples_per_readout, ncoils) + v = convert_array(v, ComplexF32, samples_per_readout, nreadouts, ncoils) Jᴴv_ptr = @ccall LIBRARY.compas_compute_jacobian_hermitian( pointer(context)::Ptr{Cvoid}, diff --git a/julia-bindings/src/bindings.cpp b/julia-bindings/src/bindings.cpp index ad4bc01..f379fda 100644 --- a/julia-bindings/src/bindings.cpp +++ b/julia-bindings/src/bindings.cpp @@ -4,10 +4,10 @@ #include "parameters/tissue.h" #include "sequences/pssfp.h" #include "simulate/derivative.h" +#include "simulate/phase_encoding.h" #include "simulate/residual.h" #include "simulate/sequence.h" #include "simulate/signal.h" -#include "simulate/phase_encoding.h" #include "trajectories/cartesian.h" #include "trajectories/spiral.h" #include "trajectories/trajectory.h" @@ -239,7 +239,7 @@ extern "C" kmm::ArrayBase* compas_magnetization_to_signal_spiral( }); } -extern "C" compas::Array* compas_compute_jacobian( +extern "C" compas::Array* compas_compute_jacobian( const compas::CudaContext* context, int ncoils, const compas::Array* echos, @@ -271,7 +271,7 @@ extern "C" compas::Array* compas_compute_jacobian( *coils, *vector); - return new compas::Array(Jv); + return new compas::Array(Jv); }); } @@ -288,7 +288,7 @@ extern "C" compas::Array* compas_compute_jacobian_hermitian( float delta_t, const compas::Array* k_start, cfloat delta_k, - const compas::Array* vector) { + const compas::Array* vector) { return catch_exceptions([&] { auto trajectory = compas::CartesianTrajectory { nreadouts, diff --git a/src/jacobian/hermitian.cu b/src/jacobian/hermitian.cu index 0258dba..7fcbb95 100644 --- a/src/jacobian/hermitian.cu +++ b/src/jacobian/hermitian.cu @@ -21,7 +21,7 @@ static __device__ void expand_readout_and_accumulate_mhv( TissueVoxel p, CartesianTrajectoryView trajectory, int readout, - cuda_view vector) { + cuda_view vector) { auto ns = trajectory.samples_per_readout; auto delta_t = trajectory.delta_t; auto delta_k = trajectory.delta_k; @@ -47,12 +47,11 @@ static __device__ void expand_readout_and_accumulate_mhv( #pragma unroll 8 for (int sample = 0; sample < ns; sample++) { - index_t t = readout * ns + sample; - // accumulate dot product in mHv #pragma unroll for (int icoil = 0; icoil < ncoils; icoil++) { - auto v = vector[icoil][t]; + auto v = vector[icoil][readout][sample]; + mHv[icoil] = add_mul_conj(mHv[icoil], v, ms); dmHv[icoil][0] = add_mul_conj(dmHv[icoil][0], v, dms[0]); dmHv[icoil][1] = add_mul_conj(dmHv[icoil][1], v, dms[1]); @@ -73,7 +72,7 @@ __global__ void jacobian_hermitian_product( TissueParametersView parameters, CartesianTrajectoryView trajectory, cuda_view coil_sensitivities, - cuda_view vector) { + cuda_view vector) { auto voxel = index_t(blockIdx.x * blockDim.x + threadIdx.x); if (voxel < parameters.nvoxels) { @@ -130,7 +129,7 @@ Array compute_jacobian_hermitian( TissueParameters parameters, CartesianTrajectory trajectory, Array coil_sensitivities, - Array vector) { + Array vector) { int ns = trajectory.samples_per_readout; int nreadouts = trajectory.nreadouts; int ncoils = coil_sensitivities.size(0); @@ -145,7 +144,8 @@ Array compute_jacobian_hermitian( COMPAS_ASSERT(coil_sensitivities.size(0) == ncoils); COMPAS_ASSERT(coil_sensitivities.size(1) == nvoxels); COMPAS_ASSERT(vector.size(0) == ncoils); - COMPAS_ASSERT(vector.size(1) == nreadouts * ns); + COMPAS_ASSERT(vector.size(1) == nreadouts); + COMPAS_ASSERT(vector.size(2) == ns); // four reconstruction parameters: T1, T2, rho_x, rho_y auto JHv = Array(4, nvoxels); diff --git a/src/jacobian/product.h b/src/jacobian/product.h index fb02508..7938d5f 100644 --- a/src/jacobian/product.h +++ b/src/jacobian/product.h @@ -5,7 +5,7 @@ namespace compas { -Array compute_jacobian( +Array compute_jacobian( const CudaContext& ctx, Array echos, Array delta_echos_T1, @@ -23,6 +23,6 @@ Array compute_jacobian_hermitian( TissueParameters parameters, CartesianTrajectory trajectory, Array coil_sensitivities, - Array vector); + Array vector); } // namespace compas \ No newline at end of file diff --git a/src/jacobian/regular.cu b/src/jacobian/regular.cu index 596c174..f395848 100644 --- a/src/jacobian/regular.cu +++ b/src/jacobian/regular.cu @@ -42,7 +42,7 @@ __global__ void delta_to_sample_exponent( template __launch_bounds__(256, 16) __global__ void jacobian_product( - cuda_view_mut Jv, + cuda_view_mut Jv, cuda_view echos, cuda_view delta_echos_T1, cuda_view delta_echos_T2, @@ -109,14 +109,14 @@ __launch_bounds__(256, 16) __global__ void jacobian_product( } if (lane_id == 0) { - Jv[icoil][t] = value; + Jv[icoil][r][s] = value; } } } } } // namespace kernels -Array compute_jacobian( +Array compute_jacobian( const CudaContext& ctx, Array echos, Array delta_echos_T1, @@ -142,7 +142,7 @@ Array compute_jacobian( COMPAS_ASSERT(vector.size(0) == 4); // four reconstruction parameters: T1, T2, rho_x, rho_y COMPAS_ASSERT(vector.size(1) == nvoxels); - auto Jv = Array(ncoils, nreadouts * ns); + auto Jv = Array(ncoils, nreadouts, ns); auto E = Array(ns, nvoxels); auto dEdT2 = Array(ns, nvoxels); diff --git a/src/simulate/phase_encoding.cu b/src/simulate/phase_encoding.cu index a26132f..eeff239 100644 --- a/src/simulate/phase_encoding.cu +++ b/src/simulate/phase_encoding.cu @@ -16,15 +16,16 @@ Array phase_encoding( dim3 block_dim = {32, 4}; dim3 grid_dim = {div_ceil(uint(nvoxels), block_dim.x), div_ceil(uint(nreadouts), block_dim.y)}; - ctx.submit_kernel(grid_dim, - block_dim, - kernels::phase_encoding, - write(ph_en_echos), - echos, - parameters, - trajectory); + ctx.submit_kernel( + grid_dim, + block_dim, + kernels::phase_encoding, + write(ph_en_echos), + echos, + parameters, + trajectory); return ph_en_echos; } -} // compas +} // namespace compas diff --git a/src/simulate/phase_encoding.h b/src/simulate/phase_encoding.h index 1ec6fd3..85ef948 100644 --- a/src/simulate/phase_encoding.h +++ b/src/simulate/phase_encoding.h @@ -1,7 +1,7 @@ #pragma once -#include "core/context.h" #include "core/complex_type.h" +#include "core/context.h" #include "parameters/tissue.h" #include "trajectories/cartesian.h" @@ -14,7 +14,6 @@ Array phase_encoding( const CudaContext& ctx, const Array& echos, const TissueParameters& parameters, - const CartesianTrajectory& trajectory - ); + const CartesianTrajectory& trajectory); -} // compas \ No newline at end of file +} // namespace compas \ No newline at end of file diff --git a/src/simulate/phase_encoding_kernels.cuh b/src/simulate/phase_encoding_kernels.cuh index 3148b89..285b599 100644 --- a/src/simulate/phase_encoding_kernels.cuh +++ b/src/simulate/phase_encoding_kernels.cuh @@ -23,5 +23,5 @@ __global__ void phase_encoding( } } -} // kernels -} // compas \ No newline at end of file +} // namespace kernels +} // namespace compas \ No newline at end of file