From db09c2f159a9b4a6cd15db454eba67e8063847b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Onur=20=C3=9Clgen?= Date: Fri, 26 Jan 2024 12:14:42 +0000 Subject: [PATCH] Refactor affine deformation field unit test #92 --- niftyreg_build_version.txt | 2 +- reg-test/reg_test_affineDeformationField.cpp | 375 ++++++++++--------- 2 files changed, 190 insertions(+), 187 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index 0ca45a09..e537bfeb 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -394 +395 diff --git a/reg-test/reg_test_affineDeformationField.cpp b/reg-test/reg_test_affineDeformationField.cpp index dd39cf4e..f38ce164 100644 --- a/reg-test/reg_test_affineDeformationField.cpp +++ b/reg-test/reg_test_affineDeformationField.cpp @@ -1,186 +1,189 @@ -#include "reg_test_common.h" - -/* - This test file contains the following unit tests: - test function: creation of a deformation field from an affine matrix - In 2D and 3D - identity - translation - affine -*/ - - -typedef std::tuple TestData; -typedef std::tuple, unique_ptr> ContentDesc; - -TEST_CASE("Affine Deformation Field", "[unit]") { - // Create a reference 2D image - int dim[8] = { 2, 2, 2, 1, 1, 1, 1, 1 }; - nifti_image *reference2d = nifti_make_new_nim(dim, NIFTI_TYPE_FLOAT32, true); - reg_checkAndCorrectDimension(reference2d); - - // Create a reference 3D image - dim[0] = 3; - dim[3] = 2; - nifti_image *reference3d = nifti_make_new_nim(dim, NIFTI_TYPE_FLOAT32, true); - reg_checkAndCorrectDimension(reference3d); - - // Generate the different test cases - vector testCases; - - // Identity use case - 2D - mat44 identity; - reg_mat44_eye(&identity); - // Test order [0,0] [1,0] [0,1] [1,1] - float identityResult2x[4] = { 0, 1, 0, 1 }; - float identityResult2y[4] = { 0, 0, 1, 1 }; - testCases.emplace_back(TestData( - "identity 2D", - reference2d, - &identity, - identityResult2x, - identityResult2y, - nullptr - )); - - // Identity use case - 3D - // Test order [0,0,0] [1,0,0] [0,1,0] [1,1,0],[0,0,1] [1,0,1] [0,1,1] [1,1,1] - float identityResult3x[8] = { 0, 1, 0, 1, 0, 1, 0, 1 }; - float identityResult3y[8] = { 0, 0, 1, 1, 0, 0, 1, 1 }; - float identityResult3z[8] = { 0, 0, 0, 0, 1, 1, 1, 1 }; - testCases.emplace_back(TestData( - "identity 3D", - reference3d, - &identity, - identityResult3x, - identityResult3y, - identityResult3z - )); - - // Translation - 2D - mat44 translation; - reg_mat44_eye(&translation); - translation.m[0][3] = -0.5; - translation.m[1][3] = 1.5; - translation.m[2][3] = 0.75; - // Test order [0,0] [1,0] [0,1] [1,1] - float translationResult2x[4] = { -0.5, .5, -0.5, .5 }; - float translationResult2y[4] = { 1.5, 1.5, 2.5, 2.5 }; - testCases.emplace_back(TestData( - "translation 2D", - reference2d, - &translation, - translationResult2x, - translationResult2y, - nullptr - )); - - // Translation - 3D - // Test order [0,0,0] [1,0,0] [0,1,0] [1,1,0],[0,0,1] [1,0,1] [0,1,1] [1,1,1] - float translationResult3x[8] = { -0.5, .5, -0.5, .5, -0.5, .5, -0.5, .5 }; - float translationResult3y[8] = { 1.5, 1.5, 2.5, 2.5, 1.5, 1.5, 2.5, 2.5 }; - float translationResult3z[8] = { .75, .75, .75, .75, 1.75, 1.75, 1.75, 1.75 }; - testCases.emplace_back(TestData( - "translation 3D", - reference3d, - &translation, - translationResult3x, - translationResult3y, - translationResult3z - )); - - // Full affine - 2D - // Test order [0,0] [1,0] [0,1] [1,1] - mat44 affine; - reg_mat44_eye(&affine); - affine.m[0][3] = -0.5; - affine.m[1][3] = 1.5; - affine.m[2][3] = 0.75; - for (int i = 0; i < 4; ++i) { - for (int j = 0; j < 4; ++j) { - affine.m[i][j] += ((static_cast(rand()) / RAND_MAX) - 0.5f) / 10.f; - } - } - float affineResult2x[4]; - float affineResult2y[4]; - for (int i = 0; i < 4; ++i) { - auto x = identityResult2x[i]; - auto y = identityResult2y[i]; - affineResult2x[i] = affine.m[0][3] + affine.m[0][0] * x + affine.m[0][1] * y; - affineResult2y[i] = affine.m[1][3] + affine.m[1][0] * x + affine.m[1][1] * y; - - } - testCases.emplace_back(TestData( - "full affine 2D", - reference2d, - &affine, - affineResult2x, - affineResult2y, - nullptr - )); - - // Full affine - 3D - // Test order [0,0,0] [1,0,0] [0,1,0] [1,1,0],[0,0,1] [1,0,1] [0,1,1] [1,1,1] - float affineResult3x[8]; - float affineResult3y[8]; - float affineResult3z[8]; - for (int i = 0; i < 8; ++i) { - auto x = identityResult3x[i]; - auto y = identityResult3y[i]; - auto z = identityResult3z[i]; - affineResult3x[i] = affine.m[0][3] + affine.m[0][0] * x + affine.m[0][1] * y + affine.m[0][2] * z; - affineResult3y[i] = affine.m[1][3] + affine.m[1][0] * x + affine.m[1][1] * y + affine.m[1][2] * z; - affineResult3z[i] = affine.m[2][3] + affine.m[2][0] * x + affine.m[2][1] * y + affine.m[2][2] * z; - } - testCases.emplace_back(TestData( - "affine 3D", - reference3d, - &affine, - affineResult3x, - affineResult3y, - affineResult3z - )); - - // Loop over all generated test cases - for (auto&& testCase : testCases) { - // Retrieve test information - auto&& [testName, reference, testMat, testResX, testResY, testResZ] = testCase; - - // Accumulate all required contents with a vector - vector contentDescs; - for (auto&& platformType : PlatformTypes) { - unique_ptr platform{ new Platform(platformType) }; - unique_ptr contentCreator{ dynamic_cast(platform->CreateContentCreator(ContentType::Aladin)) }; - unique_ptr content{ contentCreator->Create(reference, reference, nullptr, testMat, sizeof(float)) }; - contentDescs.push_back({ std::move(content), std::move(platform) }); - } - // Loop over all possibles contents for each test - for (auto&& contentDesc : contentDescs) { - auto&& [content, platform] = contentDesc; - const std::string sectionName = testName + " " + platform->GetName(); - SECTION(sectionName) { - NR_COUT << "\n**************** Section " << sectionName << " ****************" << std::endl; - - // Do the calculation - unique_ptr affineDeformKernel{ platform->CreateKernel(AffineDeformationFieldKernel::GetName(), content.get()) }; - affineDeformKernel->castTo()->Calculate(); - - // Check all values - nifti_image *defField = content->GetDeformationField(); - auto defFieldPtrX = static_cast(defField->data); - const size_t voxelNumber = NiftiImage::calcVoxelNumber(defField, 3); - auto defFieldPtrY = &defFieldPtrX[voxelNumber]; - auto defFieldPtrZ = &defFieldPtrY[voxelNumber]; - for (size_t i = 0; i < voxelNumber; ++i) { - REQUIRE(fabs(defFieldPtrX[i] - testResX[i]) < EPS); - REQUIRE(fabs(defFieldPtrY[i] - testResY[i]) < EPS); - if (testResZ) - REQUIRE(fabs(defFieldPtrZ[i] - testResZ[i]) < EPS); - } - } - } - } - // Clean up - nifti_image_free(reference2d); - nifti_image_free(reference3d); -} +#include "reg_test_common.h" + +/* + This test file contains the following unit tests: + test function: creation of a deformation field from an affine matrix + In 2D and 3D + Identity + Translation + Affine +*/ + +struct float3 { + float x, y, z; + + std::string to_string() const { + return "(" + std::to_string(x) + ", " + std::to_string(y) + ", " + std::to_string(z) + ")"; + } +}; + +class AffineDeformationFieldTest { +protected: + using TestData = std::tuple>; + using TestCase = std::tuple>; + + inline static vector testCases; + +public: + AffineDeformationFieldTest() { + if (!testCases.empty()) + return; + + // Create reference images + constexpr NiftiImage::dim_t size = 2; + NiftiImage reference2d({ size, size }, NIFTI_TYPE_FLOAT32); + NiftiImage reference3d({ size, size, size }, NIFTI_TYPE_FLOAT32); + + // Data container for the test data + vector testData; + + // Identity use case - 2D + mat44 identity; + reg_mat44_eye(&identity); + // Test order [0,0] [1,0] [0,1] [1,1] + vector identityResult2d{ { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 1, 1, 0 } }; + testData.emplace_back(TestData( + "2D Identity", + reference2d, + identity, + identityResult2d + )); + + // Identity use case - 3D + // Test order [0,0,0] [1,0,0] [0,1,0] [1,1,0],[0,0,1] [1,0,1] [0,1,1] [1,1,1] + vector identityResult3d{ { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 1, 1, 0 }, { 0, 0, 1 }, { 1, 0, 1 }, { 0, 1, 1 }, { 1, 1, 1 } }; + testData.emplace_back(TestData( + "3D Identity", + reference3d, + identity, + identityResult3d + )); + + // Translation - 2D + mat44 translation; + reg_mat44_eye(&translation); + translation.m[0][3] = -0.5; + translation.m[1][3] = 1.5; + translation.m[2][3] = 0.75; + // Test order [0,0] [1,0] [0,1] [1,1] + vector translationResult2d{ { -0.5f, 1.5f, 0 }, { 0.5f, 1.5f, 0 }, { -0.5f, 2.5f, 0 }, { 0.5f, 2.5f, 0 } }; + testData.emplace_back(TestData( + "2D Translation", + reference2d, + translation, + std::move(translationResult2d) + )); + + // Translation - 3D + // Test order [0,0,0] [1,0,0] [0,1,0] [1,1,0],[0,0,1] [1,0,1] [0,1,1] [1,1,1] + vector translationResult3d{ { -0.5f, 1.5f, 0.75f }, { 0.5f, 1.5f, 0.75f }, + { -0.5f, 2.5f, 0.75f }, { 0.5f, 2.5f, 0.75f }, + { -0.5f, 1.5f, 1.75f }, { 0.5f, 1.5f, 1.75f }, + { -0.5f, 2.5f, 1.75f }, { 0.5f, 2.5f, 1.75f } }; + testData.emplace_back(TestData( + "3D Translation", + reference3d, + translation, + std::move(translationResult3d) + )); + + // Full affine - 2D + // Test order [0,0] [1,0] [0,1] [1,1] + mat44 affine; + reg_mat44_eye(&affine); + affine.m[0][3] = -0.5; + affine.m[1][3] = 1.5; + affine.m[2][3] = 0.75; + for (int i = 0; i < 4; ++i) + for (int j = 0; j < 4; ++j) + affine.m[i][j] += ((static_cast(rand()) / RAND_MAX) - 0.5f) / 10.f; + vector affineResult2d(4); + for (int i = 0; i < 4; ++i) { + double x = identityResult2d[i].x; + double y = identityResult2d[i].y; + affineResult2d[i].x = static_cast(affine.m[0][3] + affine.m[0][0] * x + affine.m[0][1] * y); + affineResult2d[i].y = static_cast(affine.m[1][3] + affine.m[1][0] * x + affine.m[1][1] * y); + + } + testData.emplace_back(TestData( + "2D Affine", + reference2d, + affine, + std::move(affineResult2d) + )); + + // Full affine - 3D + // Test order [0,0,0] [1,0,0] [0,1,0] [1,1,0],[0,0,1] [1,0,1] [0,1,1] [1,1,1] + vector affineResult3d(8); + for (int i = 0; i < 8; ++i) { + double x = identityResult3d[i].x; + double y = identityResult3d[i].y; + double z = identityResult3d[i].z; + affineResult3d[i].x = static_cast(affine.m[0][3] + affine.m[0][0] * x + affine.m[0][1] * y + affine.m[0][2] * z); + affineResult3d[i].y = static_cast(affine.m[1][3] + affine.m[1][0] * x + affine.m[1][1] * y + affine.m[1][2] * z); + affineResult3d[i].z = static_cast(affine.m[2][3] + affine.m[2][0] * x + affine.m[2][1] * y + affine.m[2][2] * z); + } + testData.emplace_back(TestData( + "3D Affine", + reference3d, + affine, + std::move(affineResult3d) + )); + + for (auto&& testData : testData) { + for (auto&& platformType : PlatformTypes) { + // Make a copy of the test data + auto [testName, reference, transMat, expRes] = testData; + + // Create the platform + unique_ptr platform{ new Platform(platformType) }; + testName += " "s + platform->GetName(); + + // Create the content for Aladin + unique_ptr aladinContentCreator{ dynamic_cast(platform->CreateContentCreator(ContentType::Aladin)) }; + unique_ptr aladinContent{ aladinContentCreator->Create(reference, reference, nullptr, &transMat, sizeof(float)) }; + + // Do the calculation for Aladin + unique_ptr affineDeformKernel{ platform->CreateKernel(AffineDeformationFieldKernel::GetName(), aladinContent.get()) }; + affineDeformKernel->castTo()->Calculate(); + + // Get the result + NiftiImage defField(aladinContent->GetDeformationField(), NiftiImage::Copy::Image); + + // Save for testing + testCases.push_back({ testName + " - Aladin", std::move(defField), std::move(expRes) }); + } + } + } +}; + +TEST_CASE_METHOD(AffineDeformationFieldTest, "Affine Deformation Field", "[unit]") { + // Loop over all possibles contents for each test + for (auto&& testCase : testCases) { + auto&& [testName, defField, expected] = testCase; + SECTION(testName) { + NR_COUT << "\n**************** Section " << testName << " ****************" << std::endl; + + // Increase the precision for the output + NR_COUT << std::fixed << std::setprecision(10); + + // Check all values + const bool is3d = defField->nz > 1; + const size_t voxelNumber = defField.nVoxelsPerVolume(); + const auto defFieldPtrX = defField.data(0); + const auto defFieldPtrY = defField.data(1); + const auto defFieldPtrZ = defField.data(2); + for (auto i = 0; i < voxelNumber; i++) { + float3 result{ static_cast(defFieldPtrX[i]), static_cast(defFieldPtrY[i]), is3d ? defFieldPtrZ[i] : 0.f }; + float3 diff{ abs(result.x - expected[i].x), abs(result.y - expected[i].y), abs(result.z - expected[i].z) }; + if (diff.x > 0 || diff.y > 0 || diff.z > 0) { + NR_COUT << "[i]=" << i; + NR_COUT << " | diff=" << diff.to_string(); + NR_COUT << " | Result=" << result.to_string(); + NR_COUT << " | Expected=" << expected[i].to_string() << std::endl; + } + REQUIRE((diff.x == 0 && diff.y == 0 && diff.z == 0)); + } + } + } +}