From 64d7a35a391de270a7aeb28e5ca3d5664ec7f0d1 Mon Sep 17 00:00:00 2001 From: Ali Aghaeifar <11705505+aghaeifar@users.noreply.github.com> Date: Thu, 26 Dec 2024 21:03:59 +0100 Subject: [PATCH] adding unit test - clean up and ... (#20) * handler for sim module * update * adding test kernel * adding test kernel to workflow * clean up * clean up --- .github/workflows/build_cmake.yml | 3 + CMakeLists.txt | 5 +- containers/create_dockers.sh | 2 +- src/config/config_generator.h | 4 - src/config/handler.h | 5 +- src/dwi/handler.cpp | 1 - src/dwi/handler.h | 5 +- src/dwi/pgse.h | 6 +- src/phantom/handler.h | 4 + src/phantom/phantom_base.h | 2 +- src/phantom/phantom_cylinder.cpp | 2 +- src/phantom/phantom_cylinder.h | 4 +- src/phantom/phantom_sphere.cpp | 2 +- src/phantom/phantom_sphere.h | 4 +- src/sim/CMakeLists.txt | 1 + src/sim/device_helper.cu | 4 - src/sim/device_helper.cuh | 3 - src/sim/h5_helper.h | 4 +- src/sim/handler.cu | 18 ++ src/sim/handler.cuh | 22 +++ src/sim/kernels.cu | 10 +- src/sim/kernels.cuh | 2 +- src/sim/monte_carlo.cu | 2 - src/sim/monte_carlo.cuh | 2 - src/sim/simulation_parameters.cuh | 4 +- src/{spinwalk.cu => spinwalk.cpp} | 33 ++-- tests/CMakeLists.txt | 2 +- tests/create_config.sh | 77 -------- tests/create_phantom.sh | 34 ---- tests/plot_results.m | 300 ------------------------------ tests/test_kernel.cpp | 83 +++++++++ tests/test_phantom.cpp | 29 +-- 32 files changed, 182 insertions(+), 497 deletions(-) create mode 100644 src/sim/handler.cu create mode 100644 src/sim/handler.cuh rename src/{spinwalk.cu => spinwalk.cpp} (86%) delete mode 100755 tests/create_config.sh delete mode 100755 tests/create_phantom.sh delete mode 100644 tests/plot_results.m create mode 100644 tests/test_kernel.cpp diff --git a/.github/workflows/build_cmake.yml b/.github/workflows/build_cmake.yml index 33f4d74..c39cf97 100644 --- a/.github/workflows/build_cmake.yml +++ b/.github/workflows/build_cmake.yml @@ -71,3 +71,6 @@ jobs: - name: test config creation run: ${{github.workspace}}/spinwalk_test --run_test=test_config_creation + + - name: test kernel + run: ${{github.workspace}}/spinwalk_test --run_test=test_kernel diff --git a/CMakeLists.txt b/CMakeLists.txt index 84256d4..df45977 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -51,7 +51,7 @@ add_subdirectory(src/sim) # Add the sources files set(SOURCE_MODULES ${SUBCOMMAND_PHANTOM} ${SUBCOMMAND_SIM} ${SUBCOMMAND_DWI} ${SUBCOMMAND_CONFIG}) -set(SOURCE_MAIN ./src/spinwalk.cu) +set(SOURCE_MAIN ./src/spinwalk.cpp) # change extension of the files if CUDA is not found if(NOT CMAKE_CUDA_COMPILER) @@ -67,9 +67,6 @@ if(NOT CMAKE_CUDA_COMPILER) list(APPEND RENAMED_SOURCES ${NEW_FILE}) endforeach() set(SOURCE_MODULES ${RENAMED_SOURCES}) - - file(COPY_FILE ./src/spinwalk.cu ./src/spinwalk.cpp) - set(SOURCE_MAIN ./src/spinwalk.cpp) endif() # Add the executable diff --git a/containers/create_dockers.sh b/containers/create_dockers.sh index 6371b77..12fbf04 100755 --- a/containers/create_dockers.sh +++ b/containers/create_dockers.sh @@ -6,7 +6,7 @@ pushd "$SCRIPT_DIR" > /dev/null # docker system prune -a --volumes --force docker system prune --force -CUDA_VERSIONS=("12.0.0" "12.2.0" "12.6.3") +CUDA_VERSIONS=("12.0.0" "12.6.3") for VERSION in "${CUDA_VERSIONS[@]}"; do echo "Building Docker image for CUDA version $VERSION..." diff --git a/src/config/config_generator.h b/src/config/config_generator.h index 12dadea..15d313b 100644 --- a/src/config/config_generator.h +++ b/src/config/config_generator.h @@ -1,7 +1,3 @@ -// -------------------------------------------------------------------------------------------- -// generate default configuration file -// -------------------------------------------------------------------------------------------- - #include #include diff --git a/src/config/handler.h b/src/config/handler.h index f949c32..86bdf30 100644 --- a/src/config/handler.h +++ b/src/config/handler.h @@ -1,4 +1,5 @@ - +#ifndef CONFIG_HANDLER_H +#define CONFIG_HANDLER_H #include #include @@ -18,3 +19,5 @@ namespace config { static bool execute(const execute_args& args); }; } + +#endif // CONFIG_HANDLER_H \ No newline at end of file diff --git a/src/dwi/handler.cpp b/src/dwi/handler.cpp index ebfa7a5..15180a9 100644 --- a/src/dwi/handler.cpp +++ b/src/dwi/handler.cpp @@ -1,6 +1,5 @@ #include - #include "handler.h" #include "pgse.h" diff --git a/src/dwi/handler.h b/src/dwi/handler.h index 8cfa491..a4ce13e 100644 --- a/src/dwi/handler.h +++ b/src/dwi/handler.h @@ -1,4 +1,5 @@ - +#ifndef DWI_HANDLER_H +#define DWI_HANDLER_H #include #include @@ -19,3 +20,5 @@ namespace dMRI { static bool execute(const execute_args& args); }; } + +#endif // DWI_HANDLER_H \ No newline at end of file diff --git a/src/dwi/pgse.h b/src/dwi/pgse.h index 2676b3c..be823bd 100644 --- a/src/dwi/pgse.h +++ b/src/dwi/pgse.h @@ -1,4 +1,6 @@ +#ifndef PGSE_H +#define PGSE_H #include #include @@ -25,4 +27,6 @@ class pgse std::vector b_value; }; -} \ No newline at end of file +} + +#endif // PGSE_H diff --git a/src/phantom/handler.h b/src/phantom/handler.h index f327cbe..59a9e8c 100644 --- a/src/phantom/handler.h +++ b/src/phantom/handler.h @@ -1,4 +1,6 @@ +#ifndef PHANTOM_HANDLER_H +#define PHANTOM_HANDLER_H #include #include @@ -23,3 +25,5 @@ namespace phantom { static bool execute(const execute_args& args); }; } + +#endif // PHANTOM_HANDLER_H diff --git a/src/phantom/phantom_base.h b/src/phantom/phantom_base.h index fcdac77..89937bb 100644 --- a/src/phantom/phantom_base.h +++ b/src/phantom/phantom_base.h @@ -5,7 +5,7 @@ * * Author : Ali Aghaeifar * Date : 10.02.2023 - * Descrip : simulating BOLD in microvascular network + * Descrip : * -------------------------------------------------------------------------- */ #ifndef PHANTOM_BASE_H diff --git a/src/phantom/phantom_cylinder.cpp b/src/phantom/phantom_cylinder.cpp index 2398173..d18b796 100644 --- a/src/phantom/phantom_cylinder.cpp +++ b/src/phantom/phantom_cylinder.cpp @@ -1,7 +1,7 @@ /* -------------------------------------------------------------------------- * Project: SpinWalk - * File: arbitrary_gradient.h + * File: phantom_cylinder.cpp * * Author : Ali Aghaeifar * Date : 15.05.2024 diff --git a/src/phantom/phantom_cylinder.h b/src/phantom/phantom_cylinder.h index 813e23d..eb06aba 100644 --- a/src/phantom/phantom_cylinder.h +++ b/src/phantom/phantom_cylinder.h @@ -1,11 +1,11 @@ /* -------------------------------------------------------------------------- * Project: SpinWalk - * File: file_utils.cpp + * File: phantom_cylinder.h * * Author : Ali Aghaeifar * Date : 10.02.2023 - * Descrip : simulating BOLD in microvascular network + * Descrip : * -------------------------------------------------------------------------- */ #ifndef CYLINDER_H diff --git a/src/phantom/phantom_sphere.cpp b/src/phantom/phantom_sphere.cpp index bb8c3b2..1d6234b 100644 --- a/src/phantom/phantom_sphere.cpp +++ b/src/phantom/phantom_sphere.cpp @@ -1,7 +1,7 @@ /* -------------------------------------------------------------------------- * Project: SpinWalk - * File: arbitrary_gradient.h + * File: phantom_sphere.cpp * * Author : Ali Aghaeifar * Date : 15.05.2024 diff --git a/src/phantom/phantom_sphere.h b/src/phantom/phantom_sphere.h index 1b605d8..d461768 100644 --- a/src/phantom/phantom_sphere.h +++ b/src/phantom/phantom_sphere.h @@ -1,11 +1,11 @@ /* -------------------------------------------------------------------------- * Project: SpinWalk - * File: file_utils.cpp + * File: phantom_sphere.cpp * * Author : Ali Aghaeifar * Date : 10.02.2023 - * Descrip : simulating BOLD in microvascular network + * Descrip : * -------------------------------------------------------------------------- */ #ifndef SPHERE_H diff --git a/src/sim/CMakeLists.txt b/src/sim/CMakeLists.txt index 803c341..78bf913 100644 --- a/src/sim/CMakeLists.txt +++ b/src/sim/CMakeLists.txt @@ -1,6 +1,7 @@ set(SUBCOMMAND_SIM ${CMAKE_CURRENT_SOURCE_DIR}/monte_carlo.cu ${CMAKE_CURRENT_SOURCE_DIR}/kernels.cu + ${CMAKE_CURRENT_SOURCE_DIR}/handler.cu ${CMAKE_CURRENT_SOURCE_DIR}/device_helper.cu ${CMAKE_CURRENT_SOURCE_DIR}/config_reader.cpp PARENT_SCOPE) \ No newline at end of file diff --git a/src/sim/device_helper.cu b/src/sim/device_helper.cu index 892bf0b..ef1af5f 100644 --- a/src/sim/device_helper.cu +++ b/src/sim/device_helper.cu @@ -1,7 +1,4 @@ -#ifndef DEVICE_HELPER_CUH -#define DEVICE_HELPER_CUH - #include #include #include "device_helper.cuh" @@ -104,4 +101,3 @@ bool check_memory_size(size_t required_size_MB) } } // namespace sim -#endif // DEVICE_HELPER_CUH diff --git a/src/sim/device_helper.cuh b/src/sim/device_helper.cuh index 2e3601f..529b82f 100644 --- a/src/sim/device_helper.cuh +++ b/src/sim/device_helper.cuh @@ -3,9 +3,6 @@ #define DEVICE_HELPER_CUH #include -//--------------------------------------------------------------------------------------------- -// check for CUDA and GPU device -//--------------------------------------------------------------------------------------------- namespace sim { diff --git a/src/sim/h5_helper.h b/src/sim/h5_helper.h index 13e5464..f8a3a7b 100644 --- a/src/sim/h5_helper.h +++ b/src/sim/h5_helper.h @@ -1,11 +1,11 @@ /* -------------------------------------------------------------------------- * Project: SpinWalk - * File: file_utils.h + * File: h5_helper.h * * Author : Ali Aghaeifar * Date : 10.02.2023 - * Descrip : simulating BOLD in microvascular network + * Descrip : * -------------------------------------------------------------------------- */ #ifndef __H5_HELPER_H__ diff --git a/src/sim/handler.cu b/src/sim/handler.cu new file mode 100644 index 0000000..5dcb2d7 --- /dev/null +++ b/src/sim/handler.cu @@ -0,0 +1,18 @@ + +#include +#include +#include "handler.cuh" +#include "monte_carlo.cuh" + +namespace sim { + bool handler::execute(const execute_args& args) { + sim::monte_carlo mc(args.use_cpu, args.device_id); + for(const auto& config_file : args.config_files){ + std::cout << "<" << std::filesystem::path(config_file).filename().string() << ">\n"; + if(mc.run(config_file) == false) + return 1; + } + std::cout << "Simulation completed successfully. See the log file\n"; + return true; + } +} \ No newline at end of file diff --git a/src/sim/handler.cuh b/src/sim/handler.cuh new file mode 100644 index 0000000..f34dea4 --- /dev/null +++ b/src/sim/handler.cuh @@ -0,0 +1,22 @@ + +#ifndef SIM_HANDLER_H +#define SIM_HANDLER_H + +#include +#include +#include + +namespace sim { + struct execute_args { + bool use_cpu = false; + uint32_t device_id = 0; + std::vector config_files; // Default value for name + }; + + class handler { + public: + static bool execute(const execute_args& args); + }; +} + +#endif // SIM_HANDLER_H \ No newline at end of file diff --git a/src/sim/kernels.cu b/src/sim/kernels.cu index c463ed8..54f5744 100644 --- a/src/sim/kernels.cu +++ b/src/sim/kernels.cu @@ -1,10 +1,10 @@ /* -------------------------------------------------------------------------- * Project: SpinWalk - * File: kernels.cuh + * File: kernels.cu * * Author : Ali Aghaeifar * Date : 10.02.2023 - * Descrip : simulating BOLD in microvascular network + * Descrip : * -------------------------------------------------------------------------- */ #include @@ -20,12 +20,8 @@ #include #include #endif -//--------------------------------------------------------------------------------------------- -// -//--------------------------------------------------------------------------------------------- -#define ABS(x) ((x) < 0 ? -(x) : (x)) - +#define ABS(x) ((x) < 0 ? -(x) : (x)) namespace sim { diff --git a/src/sim/kernels.cuh b/src/sim/kernels.cuh index 13a466f..8efdf7b 100644 --- a/src/sim/kernels.cuh +++ b/src/sim/kernels.cuh @@ -4,7 +4,7 @@ * * Author : Ali Aghaeifar * Date : 10.02.2023 - * Descrip : simulating BOLD in microvascular network + * Descrip : * -------------------------------------------------------------------------- */ diff --git a/src/sim/monte_carlo.cu b/src/sim/monte_carlo.cu index 56d2662..4f02c13 100644 --- a/src/sim/monte_carlo.cu +++ b/src/sim/monte_carlo.cu @@ -56,8 +56,6 @@ monte_carlo::monte_carlo(bool gpu_disabled, int32_t device_id) monte_carlo::~monte_carlo() { - // if(param) delete param; - // if(config) delete config; } void monte_carlo::allocate_memory() diff --git a/src/sim/monte_carlo.cuh b/src/sim/monte_carlo.cuh index 02e3d8a..9055299 100644 --- a/src/sim/monte_carlo.cuh +++ b/src/sim/monte_carlo.cuh @@ -6,8 +6,6 @@ #include "config_reader.h" #include "simulation_parameters.cuh" -// struct simulation_parameters; - namespace sim { // class config_reader; diff --git a/src/sim/simulation_parameters.cuh b/src/sim/simulation_parameters.cuh index fc5f09d..67e3d1b 100644 --- a/src/sim/simulation_parameters.cuh +++ b/src/sim/simulation_parameters.cuh @@ -1,10 +1,10 @@ /* -------------------------------------------------------------------------- * Project: Microvascular - * File: miscellaneous.h + * File: simulation_parameters.cuh * * Author : Ali Aghaeifar * Date : 10.02.2023 - * Descrip : simulating BOLD in microvascular network + * Descrip : * -------------------------------------------------------------------------- */ #ifndef __SIMULATION_PARAMETERS_H__ diff --git a/src/spinwalk.cu b/src/spinwalk.cpp similarity index 86% rename from src/spinwalk.cu rename to src/spinwalk.cpp index 771b5fa..0b18c36 100644 --- a/src/spinwalk.cu +++ b/src/spinwalk.cpp @@ -1,15 +1,11 @@ /* -------------------------------------------------------------------------- * Project: SpinWalk - * File: spinwalk.cu + * File: spinwalk.cpp * * Author : Ali Aghaeifar * Date : 10.02.2023 - * Descrip : Monte Carlo simulation of the spin dynamics in the presence of off-resonance fields. + * Descrip : * -------------------------------------------------------------------------- */ - -// compile(lin) : nvcc ./src/spinwalk.cu ./src/kernels.cu -I ./include/ -Xptxas -v -O3 -arch=compute_75 -code=sm_75 -Xcompiler -fopenmp -o spinwalk -// compile(win) : nvcc ./src/spinwalk.cu ./src/kernels.cu -I ./include/ -Xptxas -v -O3 -arch=compute_86 -code=sm_86 -Xcompiler /openmp -std=c++17 -o spinwalk - // standard libraries #include #include @@ -24,10 +20,10 @@ // custom headers #include "definitions.h" -#include "sim/monte_carlo.cuh" #include "phantom/handler.h" #include "dwi/handler.h" #include "config/handler.h" +#include "sim/handler.cuh" #include "sim/device_helper.cuh" namespace bl = boost::log; @@ -40,7 +36,7 @@ int main(int argc, char * argv[]) float arg_radius=50.f, arg_fov=1000.f, arg_dchi=0.11e-6, arg_oxy_level=0.75, arg_ori=90.f, arg_vol_fra=4.f; uint32_t arg_res=500, device_id=0, TE_us=1000, timestep_us=10; int32_t arg_seed = -1; - std::vector config_files, phantom_files; + std::vector config_files, phantom_files; std::vector bdelta; std::vector bvector; std::vector bvalue; @@ -127,20 +123,13 @@ int main(int argc, char * argv[]) return 1; } - // ========== loop over configs and simulate ========== - if (config_files.size() == 0) - return 0; - - std::cout << "Running simulation for " << config_files.size() << " config(s)..." << "\n\n"; - - sim::monte_carlo mc(use_cpu, device_id); - for(const auto& config_file : config_files){ - std::cout << "<" << std::filesystem::path(config_file).filename().string() << ">\n"; - if(mc.run(config_file) == false){ + // ========== Monte-Carlo simulation ========== + if(subcommand_sim->parsed()){ + if (sim::handler::execute({.use_cpu=use_cpu, .device_id=device_id, .config_files=config_files}) == false){ std::cout << ERR_MSG << "Simulation failed. See the log file " << log_filename <<", Aborting...!" << "\n"; - return 1; - } - } - std::cout << "Simulation completed successfully. See the log file " << log_filename << "\n"; + return 1; + } + } + return 0; } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f6790ea..0a966fc 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -23,7 +23,7 @@ endif() # include_directories(${CMAKE_SOURCE_DIR}/include ${CMAKE_BINARY_DIR}/include ${Boost_INCLUDE_DIR} ${CMAKE_CURRENT_BINARY_DIR}) include_directories(${Boost_INCLUDE_DIRS}) -add_executable(${project_test} test_phantom.cpp test_config.cpp) +add_executable(${project_test} test_phantom.cpp test_config.cpp test_kernel.cpp) target_link_libraries(${project_test} ${Boost_LIBRARIES} SpinWalk_lib) set_target_properties(${project_test} PROPERTIES OUTPUT_NAME "spinwalk_test") diff --git a/tests/create_config.sh b/tests/create_config.sh deleted file mode 100755 index f57df23..0000000 --- a/tests/create_config.sh +++ /dev/null @@ -1,77 +0,0 @@ -#!/bin/bash - -# Get the directory where the script is located -SCRIPT_DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" - -# Change to the script directory -pushd "$SCRIPT_DIR" > /dev/null - -spinwalk config -s GRE --TE 20000 -t 50 -p ./phantoms/r8_Y0.78_vf4_ori90_fov600_res600.h5 ./phantoms/r8_Y0.85_vf4_ori90_fov600_res600.h5 -o ./gre.ini -spinwalk config -s SE --TE 20000 -t 50 -p ./phantoms/r8_Y0.78_vf4_ori90_fov600_res600.h5 ./phantoms/r8_Y0.85_vf4_ori90_fov600_res600.h5 -o ./se.ini - -# replace the default config file with the scales.ini as the parent config -sed -i 's/default_config.ini/scales.ini/g' gre.ini -sed -i 's/default_config.ini/scales.ini/g' se.ini - - -# Define the output file path -output_file="scales.ini" - -# Write to the file -cat > "$output_file" < /dev/null - -# Define the output directory -output_dir="./phantoms" -mkdir -p "$output_dir" - -vol_frac=4 -oxy_level_act=0.85 -oxy_level_rest=0.78 -dChi=0.00000011 -orientation=90 -resolution=600 -fov=600 -radius=8 - -# Define the output file -output_file_rest="${output_dir}/r${radius}_Y${oxy_level_rest}_vf${vol_frac}_ori${orientation}_fov${fov}_res${resolution}.h5" -# Call the command with the variable parameter and redirect the output -spinwalk phantom -c -r "$radius" -v "$vol_frac" -f "$fov" -z "$resolution" -d "$dChi" -y "$oxy_level_rest" -n "$orientation" -e 0 -o "$output_file_rest" - -# Define the output file -output_file_act="${output_dir}/r${radius}_Y${oxy_level_act}_vf${vol_frac}_ori${orientation}_fov${fov}_res${resolution}.h5" -# Call the command with the variable parameter and redirect the output -spinwalk phantom -c -r "$radius" -v "$vol_frac" -f "$fov" -z "$resolution" -d "$dChi" -y "$oxy_level_act" -n "$orientation" -e 0 -o "$output_file_act" - -# check phantoms are there -echo -e "\nList of existing phantoms in phantom folder..." -ls -l --block-size=M ./phantoms \ No newline at end of file diff --git a/tests/plot_results.m b/tests/plot_results.m deleted file mode 100644 index 8f42a6f..0000000 --- a/tests/plot_results.m +++ /dev/null @@ -1,300 +0,0 @@ -clc -clear - -fname{1}{1} = '../../outputs/gre_phantom_0.h5'; -fname{1}{2} = '../../outputs/gre_phantom_1.h5'; -fname{2}{1} = '../../outputs/se_phantom_0.h5'; -fname{2}{2} = '../../outputs/se_phantom_1.h5'; -% fname{3}{1} = '../../outputs/ssfp_phantom_0.h5'; -% fname{3}{2} = '../../outputs/ssfp_phantom_1.h5'; - -dim_vessel_size = 1; -dim_spin = 2; -dim_echo = 3; -dim_xyz = 4; -rad_ref_um = 53.367; - -tissue_type = 0; % 0 = extra-vascular, 1 = intra-vascular, [0,1] = combined - -signal_magnitude = cell(numel(fname), numel(fname{1})); -for seq = 1:numel(fname) - for i=1:numel(fname{seq}) - m1 = permute(h5read(fname{seq}{i}, '/M'), 4:-1:1); - scales = permute(h5read(fname{seq}{i}, '/scales'), 4:-1:1); - T = permute(h5read(fname{seq}{i}, '/T'), 4:-1:1); - - m1_t = zeros(numel(scales), 1); - for s=1:numel(scales) - ind = ismember(T(s,:), tissue_type); - m1_f = m1(s, ind(:), end, :); - m1_f = squeeze(m1_f); - m1_t(s) = abs(complex(mean(m1_f(:,1)), mean(m1_f(:,2)))); - end - signal_magnitude{seq, i} = m1_t; - end -end - -clf -vessel_radius = rad_ref_um * scales; -for seq = 1:numel(fname) - relative_signal = 100 * (1 - signal_magnitude{seq, 1} ./ signal_magnitude{seq, 2}); - h = semilogx(vessel_radius, relative_signal); xlabel('Vessel radius (um)'); ylabel('BOLD Signal %'); - hold on; - % ylim([0, 7]) -end -legend('GRE', 'SE', 'SSFP') - - -%% grase sequence - -% fname{4}{1} = '../../outputs/grase_m1_0.dat'; -% fname{4}{2} = '../../outputs/grase_m1_1.dat'; -TR = 0.2; -EcoSpc = 0.005; -% subplot(2,1,2) -% relative_signal = signal_magnitude{seq}(:,:,:,:,2) - signal_magnitude{seq}(:,:,:,:,1); -% relative_signal = 100 * (1 - signal_magnitude{seq}(:,:,:,:,1)./ signal_magnitude{seq}(:,:,:,:,2)); -% relative_signal = squeeze(relative_signal); -% imagesc(relative_signal); -% axis image -% set(gca,'YDir','normal') -% -% ind = floor(linspace(1, numel(scales), 6)); -% set(gca,'XTick', ind, 'XTickLabel', round(vessel_radius(ind)), 'color','none'); box off; -% xlabel('Vessl radius (um)'); -% -% ind = floor(linspace(1, TR/EcoSpc, 8)); -% set(gca,'YTick', ind, 'YTickLabel', EcoSpc * ind * 1e3, 'color','none'); -% ylabel('Echo Time (ms)'); -% -% % title('BOLD Signal') -% colormap(inferno); -% colorbar - -%% stimulated echo -clc -clear -close all -fname{1} = '../../outputs/ste_phantom_0.h5'; -fname{2} = '../../outputs/ste_phantom_1.h5'; - -dim_vessel_size = 1; -dim_spin = 2; -dim_echo = 3; -dim_xyz = 4; -rad_ref_um = 53.367; - -tissue_type = 0; % 0 = extra-vascular, 1 = intra-vascular, [0,1] = combined - -spins_xyz = []; -signal_magnitude = cell(numel(fname), 1); -for i=1:numel(fname) - m1 = permute(h5read(fname{i}, '/M'), 4:-1:1); - scales = permute(h5read(fname{i}, '/scales'), 4:-1:1); - T = permute(h5read(fname{i}, '/T'), 4:-1:1); - m1_t = zeros(numel(scales), size(m1, dim_echo)-1); - for s=1:numel(scales) - ind = ismember(T(s,:,end-1), tissue_type); - m1_f = m1(s, ind(:), 1:end-1, :); - m1_f = squeeze(m1_f); - m1_t(s, :) = abs(complex(sum(m1_f(:,:,1)), sum(m1_f(:,:,2)))); - end - signal_magnitude{i} = m1_t; -end -vessel_radius = rad_ref_um * scales; - -relative_signal = 100 * (1 - signal_magnitude{1} ./ signal_magnitude{2}); -difference_signal = signal_magnitude{2} - signal_magnitude{1}; - -figure(1) -subplot(1,2,1) -cla -semilogx(vessel_radius, squeeze(relative_signal(:,1))); -hold on -semilogx(vessel_radius, squeeze(relative_signal(:,2))); xlabel('Vessel radius (um)'); ylabel('100 * (1 - S_r_e_s_t / S_a_c_t )'); -hold off -legend('SE', 'STE') - -subplot(1,2,2) -cla -semilogx(vessel_radius, squeeze(difference_signal(:,1))); -hold on -semilogx(vessel_radius, squeeze(difference_signal(:,2))); xlabel('Vessel radius (um)'); ylabel('S_a_c_t - S_r_e_s_t'); -hold off -legend('SE', 'STE') - -figure(2) -cla -plot(squeeze(signal_magnitude{1}(:,2))) -hold on -plot(squeeze(signal_magnitude{2}(:,2))) - -plot(squeeze(signal_magnitude{1}(:,1))) -plot(squeeze(signal_magnitude{2}(:,1))) -hold off -xlabel('Vessel radius (um)'); ylabel('Magnitude'); -legend('STE Rest', 'STE Act', 'SE Rest', 'SE Act') - -%% linear gradient calculations -clc -fov = [0.0e-3, 0.0e-3, 0.0e-3 - 0.9e-3, 0.9e-3, 0.9e-3]; -phase_range = 2*pi; % gradient will create this phase range in FoV -l = diff(fov); -time_step = 50e-6; -gamma = 267515315; % rad/s.T -if range(l) ~= 0 - error('fov must be isotropic'); -end - -G = phase_range / (gamma * l(1) * time_step); % T/m -disp(['Gradient strength = ' num2str(G * 1000) ' mT/m']) -% generate initial XYZ0 positions and save to h5 -x = write_xyz0(fov, 101*101*101, '/DATA/aaghaeifar/Nextcloud/Projects/microvascular/field_maps/xyz0_allzero.h5'); - -%% linear gradient plot -clear -fname = '/DATA/aaghaeifar/Nextcloud/Projects/microvascular/outputs/gradient_phantom_allzero.h5'; -m_xyz = permute(h5read(fname, '/M'), 4:-1:1); -m_xy = squeeze(complex(m_xyz(1,:,end-1,1), m_xyz(1,:,end-1,2))); -m_xy = double(m_xy); -sz = nthroot(numel(m_xy), 3); -m_xy = reshape(m_xy, [sz, sz, sz]) ; -phase = angle(m_xy .* exp(-i*pi)); % shift the range to [-pi pi] -close all -vin(phase) % 3D plot - -%% random-walk trajectory -clear -clc; -filename = '../../outputs/gre_trajectory_phantom_0.h5'; -xyz_all = permute(h5read(filename, '/XYZ'), 4:-1:1); -size(xyz_all) -xyz_all = xyz_all * 1e6; -n_spins = size(xyz_all, 3); -for s=round(linspace(1, n_spins, 10)) - xyz = xyz_all(:,:,s,2); - plot3(xyz(1,:), xyz(2,:), xyz(3,:)) - hold on; -end -hold off - -% nonzeros(xyz_all) - -%% Calculate Gradient Strength for certain b-value in diffusion simulations -clc -b = 100:100:5000; % s/mm^2 -gamma = 267.5153151e6; -d = 50e-6; % gradient duration -D = 30e-3; % distance between gradients - -G = b * 1e6 ./ (gamma * gamma * d * d * (D-d/3)); -G = sqrt(G); % T/m -G = G * 1000; % mT/m -disp(G') - -for i = 1:length(b) - % Open text file for writing - filename = sprintf('/DATA/aaghaeifar/Nextcloud/Projects/microvascular/SpinWalk/config/dwi/dwi_x_%d.ini', b(i)); - fid = fopen(filename, 'w'); - % Write to text file - fprintf(fid, 'SEQ_NAME = dwi_x_%d\n', b(i)); - fprintf(fid, '[PARENT]\n'); - fprintf(fid, 'PARENT_CONFIG = dwi_baseline.ini\n\n'); - fprintf(fid, '[SCAN_PARAMETERS]\n'); - fprintf(fid, '; Apply gradient in mT/m for each axis. Gradients are active for one DWELL_TIME\n'); - fprintf(fid, 'GRADIENT_XYZ[0] = %.5f 0.0 0.0\n', G(i)); - fprintf(fid, 'GRADIENT_XYZ[1] = %.5f 0.0 0.0\n', G(i)); - % Close text file - fclose(fid); -end - - -%% diffusion -clc -clear - -folder = '../../outputs/dwi'; -base = 'dwi_base_free_diffusion_model.h5'; - -b = 100:100:5000; - -m1 = permute(h5read(fullfile(folder, base), '/M'), 4:-1:1); -s0 = abs(complex(sum(m1(1,:,end-1,1)), sum(m1(1,:,end-1,2)))); -s_dwi = zeros(numel(b), 1); -for i=1:numel(b) - f = fullfile(folder, ['dwi_x_' num2str(b(i)) '_free_diffusion_model.h5' ]); - m1 = permute(h5read(f, '/M'), 4:-1:1); - s_dwi(i) = abs(complex(sum(m1(1,:,end-1,1)), sum(m1(1,:,end-1,2)))); -end -s = [s0; s_dwi]; -b = [0, b]; - -[xData, yData] = prepareCurveData( double(b(:)), double(s(:)) ); - -% Virtual b-value -b = 100:100:5000; -f = fullfile(folder, ['dwi_x_' num2str(b(1)) '_free_diffusion_model.h5' ]); -m1 = permute(h5read(f, '/M'), 4:-1:1); -xyz = permute(h5read(f, '/XYZ'), 4:-1:1); -[theta, rho] = cart2pol(m1(1,:,end-1,1), m1(1,:,end-1,2)); -G_scale = sqrt(b / b(1)); -for i=1:numel(b) - % sig = squeeze(rho) .* exp(1j*(squeeze(theta) * G_scale(i))); - % s_dwi(i) = abs(sum(sig)); - [x, y] = pol2cart(squeeze(theta) * G_scale(i), squeeze(rho)); - s_dwi(i) = abs(complex(sum(x), sum(y))); -end -s2 = [s0; s_dwi]; -b2 = [0, b]; - -[xData_v, yData_v] = prepareCurveData( double(b2(:)), double(s2(:)) ); - -% Set up fittype and options. -ft = fittype( 'exp1' ); -opts = fitoptions( 'Method', 'NonlinearLeastSquares' ); -opts.Algorithm = 'Levenberg-Marquardt'; -opts.DiffMinChange = 1e-09; -opts.Display = 'Off'; -opts.MaxFunEvals = 6000; -opts.MaxIter = 4000; -opts.StartPoint = [988912.916299669 -0.000461964926600495]; -opts.TolFun = 1e-09; -opts.TolX = 1e-09; - -% Fit model to data. -[fitresult, gof] = fit( xData, yData, ft, opts ); -[fitresult_v, gof_v] = fit( xData_v, yData_v, ft, opts ); - -% Plot fit with data. -plot( fitresult, xData, yData); hold on; -plot(fitresult_v, xData_v, yData_v); hold off; -legend('S vs. b', 'dwi', 'S_v vs. b', 'dwi_v','Location', 'NorthEast', 'Interpreter', 'none' ); -% Label axes -xlabel( 'b', 'Interpreter', 'none' ); -ylabel( 'S', 'Interpreter', 'none' ); -grid on -disp(['Diff. Const = ' num2str(abs(fitresult.b) * 1e-6) ' m/s']); -disp(['Diff. Const v = ' num2str(abs(fitresult_v.b) * 1e-6) ' m/s']); - - -%% - - - - - - - - - - - - - - - - - - - diff --git a/tests/test_kernel.cpp b/tests/test_kernel.cpp new file mode 100644 index 0000000..8cb5723 --- /dev/null +++ b/tests/test_kernel.cpp @@ -0,0 +1,83 @@ +#include +#include +#include + +#include "../src/sim/kernels.cuh" + +BOOST_AUTO_TEST_SUITE(test_kernel) + +BOOST_AUTO_TEST_CASE(test_sub2ind_3d_row_major) { + int64_t x = 1, y = 2, z = 3; + int64_t len_dim_x = 10, len_dim_y = 10, len_dim_z = 10; + + int64_t expected = x * len_dim_z * len_dim_y + y * len_dim_z + z; + BOOST_CHECK_EQUAL(sim::sub2ind(x, y, z, len_dim_x, len_dim_y, len_dim_z), expected); +} + +BOOST_AUTO_TEST_CASE(test_sub2ind_4d_row_major) { + int64_t x = 1, y = 2, z = 3, o = 4; + int64_t len_dim_x = 10, len_dim_y = 10, len_dim_z = 10, len_dim_o = 10; + + int64_t expected = x * len_dim_o * len_dim_z * len_dim_y + + y * len_dim_z * len_dim_o + + z * len_dim_o + + o; + BOOST_CHECK_EQUAL(sim::sub2ind(x, y, z, o, len_dim_x, len_dim_y, len_dim_z, len_dim_o), expected); +} + +BOOST_AUTO_TEST_CASE(test_xrot) { + float m0[3] = {0.0f, 0.0f, 1.0f}; + float m1[3] = {0.0f, 0.0f, 0.0f}; + float theta = 90.0f; + float sin_theta = std::sin(theta * M_PI / 180.0f); + float cos_theta = std::cos(theta * M_PI / 180.0f); + + sim::xrot(sin_theta, cos_theta, m0, m1); + + BOOST_TEST(std::abs(m1[0] - 0.0f) < 1e-5); + BOOST_TEST(std::abs(m1[1] - (-1.0f)) < 1e-5); + BOOST_TEST(std::abs(m1[2] - 0.0f) < 1e-5); +} + +BOOST_AUTO_TEST_CASE(test_yrot) { + float m0[3] = {1.0f, 0.0f, 0.0f}; + float m1[3] = {0.0f, 0.0f, 0.0f}; + float theta = 90.0f; + float sin_theta = std::sin(theta * M_PI / 180.0f); + float cos_theta = std::cos(theta * M_PI / 180.0f); + + sim::yrot(sin_theta, cos_theta, m0, m1); + + BOOST_TEST(std::abs(m1[0] - 0.0f) < 1e-5); + BOOST_TEST(std::abs(m1[1] - 0.0f) < 1e-5); + BOOST_TEST(std::abs(m1[2] - (-1.0f)) < 1e-5); +} + +BOOST_AUTO_TEST_CASE(test_zrot) { + float m0[3] = {1.0f, 0.0f, 0.0f}; + float m1[3] = {0.0f, 0.0f, 0.0f}; + float theta = 90.0f; + float sin_theta = std::sin(theta * M_PI / 180.0f); + float cos_theta = std::cos(theta * M_PI / 180.0f); + + sim::zrot(sin_theta, cos_theta, m0, m1); + + BOOST_TEST(std::abs(m1[0] - 0.0f) < 1e-5); + BOOST_TEST(std::abs(m1[1] - 1.0f) < 1e-5); + BOOST_TEST(std::abs(m1[2] - 0.0f) < 1e-5); +} + +BOOST_AUTO_TEST_CASE(test_relax) { + float m0[3] = {1.0f, 0.5f, -0.5f}; + float m1[3] = {0.0f, 0.0f, 0.0f}; + float e1 = 0.9f; + float e2 = 0.8f; + + sim::relax(e1, e2, m0, m1); + + BOOST_TEST(std::abs(m1[0] - m0[0] * e2) < 1e-5); + BOOST_TEST(std::abs(m1[1] - m0[1] * e2) < 1e-5); + BOOST_TEST(std::abs(m1[2] - (1.0f + e1 * (m0[2] - 1.0f))) < 1e-5); +} + +BOOST_AUTO_TEST_SUITE_END() \ No newline at end of file diff --git a/tests/test_phantom.cpp b/tests/test_phantom.cpp index 77ae6bd..b727401 100644 --- a/tests/test_phantom.cpp +++ b/tests/test_phantom.cpp @@ -11,31 +11,20 @@ BOOST_AUTO_TEST_CASE(cylinder_creation) { boost::log::core::get()->set_logging_enabled(false); float diff_margin = 2.f; - std::vector volume_fraction = {5.0, 10.0}; - for(const auto &vf : volume_fraction){ - phantom::cylinder cyl(600, 300, 0.11e-6, -1, -20, vf, 0, 0, ""); - BOOST_TEST(cyl.run(false)); - BOOST_TEST(std::abs(cyl.get_actual_volume_fraction() - vf) < diff_margin); - } - - std::vector angles = {5.0, 10.0}; - float vf = 10.0; - for(const auto &a : angles){ - phantom::cylinder cyl(600, 300, 0.11e-6, -1, -20, vf, 0, a, ""); - BOOST_TEST(cyl.run(false)); - BOOST_TEST(std::abs(cyl.get_actual_volume_fraction() - vf) < diff_margin); - } + float vf = 10.0f; + float angles = 10.0f; + phantom::cylinder cyl(600, 300, 0.11e-6, -1, -20, vf, 0, angles, ""); + BOOST_TEST(cyl.run(false)); + BOOST_TEST(std::abs(cyl.get_actual_volume_fraction() - vf) < diff_margin); } BOOST_AUTO_TEST_CASE(sphere_creation) { boost::log::core::get()->set_logging_enabled(false); - std::vector volume_fraction = {5.0, 10.0}; float diff_margin = 2.f; - for(const auto &vf : volume_fraction){ - phantom::sphere sph(600, 300, 0.11e-6, -1, -20, vf, 0, ""); - BOOST_CHECK(sph.run(false)); - BOOST_CHECK(std::abs(sph.get_actual_volume_fraction() - vf) < diff_margin); - } + float vf = 12.0; + phantom::sphere sph(600, 300, 0.11e-6, -1, -20, vf, 0, ""); + BOOST_CHECK(sph.run(false)); + BOOST_CHECK(std::abs(sph.get_actual_volume_fraction() - vf) < diff_margin); }