From 540f10b13c4785f9d9169b9c345ee35f0a5989d9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Onur=20=C3=9Clgen?= Date: Mon, 15 Jan 2024 17:34:19 +0000 Subject: [PATCH] Add symmetrise velocity fields regression test #92 --- niftyreg_build_version.txt | 2 +- reg-test/CMakeLists.txt | 1 + ...reg_test_regr_symmetriseVelocityFields.cpp | 158 ++++++++++++++++++ 3 files changed, 160 insertions(+), 1 deletion(-) create mode 100644 reg-test/reg_test_regr_symmetriseVelocityFields.cpp diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index f1386578..e45b99e9 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -383 +384 diff --git a/reg-test/CMakeLists.txt b/reg-test/CMakeLists.txt index b08293d5..1e0304ab 100755 --- a/reg-test/CMakeLists.txt +++ b/reg-test/CMakeLists.txt @@ -130,6 +130,7 @@ if(USE_CUDA) set(EXEC_LIST reg_test_regr_kernelConvolution ${EXEC_LIST}) set(EXEC_LIST reg_test_regr_lts ${EXEC_LIST}) set(EXEC_LIST reg_test_regr_measure ${EXEC_LIST}) + set(EXEC_LIST reg_test_regr_symmetriseVelocityFields ${EXEC_LIST}) endif(USE_CUDA) foreach(EXEC ${EXEC_LIST}) diff --git a/reg-test/reg_test_regr_symmetriseVelocityFields.cpp b/reg-test/reg_test_regr_symmetriseVelocityFields.cpp new file mode 100644 index 00000000..d7149814 --- /dev/null +++ b/reg-test/reg_test_regr_symmetriseVelocityFields.cpp @@ -0,0 +1,158 @@ +#include "reg_test_common.h" +#include "CudaF3dContent.h" + +/** + * Symmetrise velocity fields regression test to ensure the CPU and CUDA versions yield the same output +**/ + +class SymmetriseVelocityFieldsTest { +protected: + using TestData = std::tuple; + using TestCase = std::tuple; + + inline static vector testCases; + +public: + SymmetriseVelocityFieldsTest() { + if (!testCases.empty()) + return; + + // Create a random number generator + std::mt19937 gen(0); + std::uniform_real_distribution distr(-1, 1); + + // Create 2D and 3D reference images + constexpr NiftiImage::dim_t dimSize = 4; + NiftiImage reference2d({ dimSize, dimSize }, NIFTI_TYPE_FLOAT32); + NiftiImage reference3d({ dimSize, dimSize, dimSize }, NIFTI_TYPE_FLOAT32); + + // Create 2D and 3D control point grids + NiftiImage controlPointGrid2d = CreateControlPointGrid(reference2d); + NiftiImage controlPointGridBw2d = CreateControlPointGrid(reference2d); + NiftiImage controlPointGrid3d = CreateControlPointGrid(reference3d); + NiftiImage controlPointGridBw3d = CreateControlPointGrid(reference3d); + + // Add random values to the control point grid coefficients + // No += or + operator for RNifti::NiftiImageData:Element + // so reverting to old school for now + float *cpp2dPtr = static_cast(controlPointGrid2d->data); + float *cpp2dBwPtr = static_cast(controlPointGridBw2d->data); + float *cpp3dPtr = static_cast(controlPointGrid3d->data); + float *cpp3dBwPtr = static_cast(controlPointGridBw3d->data); + for (size_t i = 0; i < controlPointGrid2d.nVoxels(); ++i) { + cpp2dPtr[i] += distr(gen); + cpp2dBwPtr[i] += distr(gen); + } + for (size_t i = 0; i < controlPointGrid3d.nVoxels(); ++i) { + cpp3dPtr[i] += distr(gen); + cpp3dBwPtr[i] += distr(gen); + } + + // Create the affine matrices and fill them with random values + std::array matrices{}; + for (int i = 0; i < matrices.size(); ++i) + for (int j = 0; j < 4; ++j) + for (int k = 0; k < 4; ++k) + matrices[i].m[j][k] = j == k ? distr(gen) : 0; + + // Add the test data + vector testData; + testData.emplace_back(TestData( + "2D", + std::move(reference2d), + std::move(controlPointGrid2d), + std::move(controlPointGridBw2d) + )); + testData.emplace_back(TestData( + "3D", + std::move(reference3d), + std::move(controlPointGrid3d), + std::move(controlPointGridBw3d) + )); + + // Create the platforms + Platform platformCpu(PlatformType::Cpu); + Platform platformCuda(PlatformType::Cuda); + + for (auto&& testData : testData) { + // Make a copy of the test data + auto [testName, reference, controlPointGrid, controlPointGridBw] = testData; + + // Set the affine matrices + controlPointGrid->sform_code = 0; + controlPointGrid->qto_xyz = matrices[0]; + controlPointGridBw->sform_code = 1; + controlPointGridBw->sto_xyz = matrices[1]; + + // Create images + NiftiImage referenceCpu(reference), referenceCuda(reference); + NiftiImage cppCpu(controlPointGrid), cppCuda(controlPointGrid); + NiftiImage cppBwCpu(controlPointGrid), cppBwCuda(controlPointGrid); + + // Create the content + unique_ptr contentCpu{ new F3dContent(referenceCpu, referenceCpu, cppCpu) }; + unique_ptr contentBwCpu{ new F3dContent(referenceCpu, referenceCpu, cppBwCpu) }; + unique_ptr contentCuda{ new CudaF3dContent(referenceCuda, referenceCuda, cppCuda) }; + unique_ptr contentBwCuda{ new CudaF3dContent(referenceCuda, referenceCuda, cppBwCuda) }; + + // Create the computes + unique_ptr computeCpu{ platformCpu.CreateCompute(*contentCpu) }; + unique_ptr computeCuda{ platformCuda.CreateCompute(*contentCuda) }; + + // Symmetrise the velocity fields + computeCpu->SymmetriseVelocityFields(*contentBwCpu); + computeCuda->SymmetriseVelocityFields(*contentBwCuda); + + // Get the results of CUDA since CPU results are already inplace + contentCuda->GetControlPointGrid(); + contentBwCuda->GetControlPointGrid(); + + // Save for testing + testCases.push_back({ testName, std::move(cppCpu), std::move(cppBwCpu), std::move(cppCuda), std::move(cppBwCuda) }); + } + } +}; + +TEST_CASE_METHOD(SymmetriseVelocityFieldsTest, "Regression Symmetrise Velocity Fields", "[regression]") { + // Loop over all generated test cases + for (auto&& testCase : testCases) { + // Retrieve test information + auto&& [sectionName, cppCpu, cppBwCpu, cppCuda, cppBwCuda] = testCase; + + SECTION(sectionName) { + NR_COUT << "\n**************** Section " << sectionName << " ****************" << std::endl; + + // Increase the precision for the output + NR_COUT << std::fixed << std::setprecision(10); + + // Check the results + const auto cppCpuPtr = cppCpu.data(); + const auto cppBwCpuPtr = cppBwCpu.data(); + const auto cppCudaPtr = cppCuda.data(); + const auto cppBwCudaPtr = cppBwCuda.data(); + for (size_t i = 0; i < cppCpu.nVoxels(); i++) { + const float cppCpuVal = cppCpuPtr[i]; + const float cppCudaVal = cppCudaPtr[i]; + const float diff = abs(cppCpuVal - cppCudaVal); + if (diff > 0) { + NR_COUT << "[i]=" << i; + NR_COUT << " | diff=" << diff; + NR_COUT << " | CPU=" << cppCpuVal; + NR_COUT << " | CUDA=" << cppCudaVal << std::endl; + } + REQUIRE(diff == 0); + // Check the results of the backwards + const float cppBwCpuVal = cppBwCpuPtr[i]; + const float cppBwCudaVal = cppBwCudaPtr[i]; + const float diffBw = abs(cppBwCpuVal - cppBwCudaVal); + if (diffBw > 0) { + NR_COUT << "[i]=" << i; + NR_COUT << " | diffBw=" << diffBw; + NR_COUT << " | CPU=" << cppBwCpuVal; + NR_COUT << " | CUDA=" << cppBwCudaVal << std::endl; + } + REQUIRE(diffBw == 0); + } + } + } +}