diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index b4eed3b8..cf7ff50f 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -358 +359 diff --git a/reg-test/CMakeLists.txt b/reg-test/CMakeLists.txt index 941ed995..b08293d5 100755 --- a/reg-test/CMakeLists.txt +++ b/reg-test/CMakeLists.txt @@ -124,6 +124,7 @@ set(EXEC_LIST reg_test_normaliseGradient ${EXEC_LIST}) set(EXEC_LIST reg_test_regr_getDeformationField ${EXEC_LIST}) set(EXEC_LIST reg_test_voxelCentricToNodeCentric ${EXEC_LIST}) if(USE_CUDA) + set(EXEC_LIST reg_test_regr_approxBendingEnergyGradient ${EXEC_LIST}) set(EXEC_LIST reg_test_regr_approxLinearEnergyGradient ${EXEC_LIST}) set(EXEC_LIST reg_test_regr_blockMatching ${EXEC_LIST}) set(EXEC_LIST reg_test_regr_kernelConvolution ${EXEC_LIST}) diff --git a/reg-test/reg_test_regr_approxBendingEnergyGradient.cpp b/reg-test/reg_test_regr_approxBendingEnergyGradient.cpp new file mode 100644 index 00000000..a2a01bdf --- /dev/null +++ b/reg-test/reg_test_regr_approxBendingEnergyGradient.cpp @@ -0,0 +1,154 @@ +#include "reg_test_common.h" +#include "CudaF3dContent.h" + +/** + * Approximate bending energy and approximate bending energy gradient regression tests + * to ensure the CPU and CUDA versions yield the same output +**/ + +class ApproxBendingEnergyGradientTest { +protected: + using TestData = std::tuple; + using TestCase = std::tuple; + + inline static vector testCases; + +public: + ApproxBendingEnergyGradientTest() { + if (!testCases.empty()) + return; + + // Create a random number generator + std::mt19937 gen(0); + std::uniform_real_distribution distr(0, 10); + + // Create 2D reference, floating and control point grid images + constexpr NiftiImage::dim_t size = 4; + vector dim{ size, size }; + NiftiImage reference2d(dim, NIFTI_TYPE_FLOAT32); + NiftiImage floating2d(dim, NIFTI_TYPE_FLOAT32); + NiftiImage controlPointGrid = CreateControlPointGrid(reference2d); + NiftiImage controlPointGrid2d[3]{ controlPointGrid, controlPointGrid, controlPointGrid }; + + // Create 3D reference, floating and control point grid images + dim.push_back(size); + NiftiImage reference3d(dim, NIFTI_TYPE_FLOAT32); + NiftiImage floating3d(dim, NIFTI_TYPE_FLOAT32); + controlPointGrid = CreateControlPointGrid(reference3d); + NiftiImage controlPointGrid3d[3]{ controlPointGrid, controlPointGrid, controlPointGrid }; + + // Fill control point grids with random values + for (int i = 0; i < 3; i++) { + auto controlPointGridPtr = controlPointGrid2d[i].data(); + for (size_t j = 0; j < controlPointGrid2d[i].nVoxels(); j++) + controlPointGridPtr[j] = distr(gen); + controlPointGridPtr = controlPointGrid3d[i].data(); + for (size_t j = 0; j < controlPointGrid3d[i].nVoxels(); j++) + controlPointGridPtr[j] = distr(gen); + } + + // Create the data container for the regression test + vector testData; + for (int i = 0; i < 3; i++) { + const float weight = distr(gen); + testData.emplace_back(TestData( + "2D weight: "s + std::to_string(weight), + reference2d, + floating2d, + controlPointGrid2d[i], + weight + )); + testData.emplace_back(TestData( + "3D weight: "s + std::to_string(weight), + reference3d, + floating3d, + controlPointGrid3d[i], + weight + )); + } + + // Create the platforms + Platform platformCpu(PlatformType::Cpu); + Platform platformCuda(PlatformType::Cuda); + + for (auto&& testData : testData) { + // Get the test data + auto&& [testName, reference, floating, controlPointGrid, weight] = testData; + + // Create images + NiftiImage referenceCpu(reference), referenceCuda(reference); + NiftiImage floatingCpu(floating), floatingCuda(floating); + NiftiImage controlPointGridCpu(controlPointGrid), controlPointGridCuda(controlPointGrid); + + // Create the contents + unique_ptr contentCpu{ new F3dContent( + referenceCpu, + floatingCpu, + controlPointGridCpu, + nullptr, + nullptr, + nullptr, + sizeof(float) + ) }; + unique_ptr contentCuda{ new CudaF3dContent( + referenceCuda, + floatingCuda, + controlPointGridCuda, + nullptr, + nullptr, + nullptr, + sizeof(float) + ) }; + + // Create the computes + unique_ptr computeCpu{ platformCpu.CreateCompute(*contentCpu) }; + unique_ptr computeCuda{ platformCuda.CreateCompute(*contentCuda) }; + + // Compute the approximate bending energy for CPU and CUDA + const double approxBendingEnergyCpu = computeCpu->ApproxBendingEnergy(); + const double approxBendingEnergyCuda = computeCuda->ApproxBendingEnergy(); + + // Compute the approximate bending energy gradient for CPU and CUDA + computeCpu->ApproxBendingEnergyGradient(weight); + computeCuda->ApproxBendingEnergyGradient(weight); + + // Get the transformation gradients + NiftiImage transGradCpu(contentCpu->GetTransformationGradient(), NiftiImage::Copy::Image); + NiftiImage transGradCuda(contentCuda->GetTransformationGradient(), NiftiImage::Copy::Image); + + // Save for testing + testCases.push_back({ testName, approxBendingEnergyCpu, approxBendingEnergyCuda, std::move(transGradCpu), std::move(transGradCuda) }); + } + } +}; + +TEST_CASE_METHOD(ApproxBendingEnergyGradientTest, "Regression Approximate Bending Energy Gradient", "[regression]") { + // Loop over all generated test cases + for (auto&& testCase : testCases) { + // Retrieve test information + auto&& [testName, approxBendingEnergyCpu, approxBendingEnergyCuda, transGradCpu, transGradCuda] = testCase; + + SECTION(testName) { + NR_COUT << "\n**************** Section " << testName << " ****************" << std::endl; + + // Increase the precision for the output + NR_COUT << std::fixed << std::setprecision(10); + + // Check the approximate bending energy values + NR_COUT << "Approx Bending Energy: " << approxBendingEnergyCpu << " " << approxBendingEnergyCuda << std::endl; + REQUIRE(abs(approxBendingEnergyCpu - approxBendingEnergyCuda) < EPS); + + // Check the transformation gradients + const auto transGradCpuPtr = transGradCpu.data(); + const auto transGradCudaPtr = transGradCuda.data(); + for (size_t i = 0; i < transGradCpu.nVoxels(); ++i) { + const float cpuVal = transGradCpuPtr[i]; + const float cudaVal = transGradCudaPtr[i]; + const auto diff = abs(cpuVal - cudaVal); + if (diff > 0) + NR_COUT << i << " " << cpuVal << " " << cudaVal << std::endl; + REQUIRE(diff < EPS); + } + } + } +} diff --git a/reg-test/reg_test_regr_approxLinearEnergyGradient.cpp b/reg-test/reg_test_regr_approxLinearEnergyGradient.cpp index 1cf5b166..530d404b 100644 --- a/reg-test/reg_test_regr_approxLinearEnergyGradient.cpp +++ b/reg-test/reg_test_regr_approxLinearEnergyGradient.cpp @@ -136,7 +136,7 @@ TEST_CASE_METHOD(ApproxLinearEnergyGradientTest, "Regression Approximate Linear // Increase the precision for the output NR_COUT << std::fixed << std::setprecision(10); - // Check the approximate linear energy + // Check the approximate linear energy values NR_COUT << "Approx Linear Energy: " << approxLinearEnergyCpu << " " << approxLinearEnergyCuda << std::endl; REQUIRE(abs(approxLinearEnergyCpu - approxLinearEnergyCuda) < EPS);