Skip to content

Commit

Permalink
Add regression tests for kernel convolution #92
Browse files Browse the repository at this point in the history
  • Loading branch information
onurulgen committed Oct 9, 2023
1 parent 47c4a84 commit 9f24fa1
Show file tree
Hide file tree
Showing 8 changed files with 193 additions and 23 deletions.
2 changes: 1 addition & 1 deletion niftyreg_build_version.txt
Original file line number Diff line number Diff line change
@@ -1 +1 @@
336
337
18 changes: 9 additions & 9 deletions reg-lib/cpu/_reg_tools.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -834,7 +834,7 @@ void reg_tools_kernelConvolution(nifti_image *image,
const int& kernelType,
const int *mask,
const bool *timePoints,
const bool *axis) {
const bool *axes) {
if (image->nx > 2048 || image->ny > 2048 || image->nz > 2048)
NR_FATAL_ERROR("This function does not support images with dimensions larger than 2048");

Expand Down Expand Up @@ -867,7 +867,7 @@ void reg_tools_kernelConvolution(nifti_image *image,
}
// Loop over the x, y and z dimensions
for (int n = 0; n < 3; n++) {
if (axis[n] && image->dim[n] > 1) {
if (axes[n] && image->dim[n] > 1) {
double temp;
if (sigma[t] > 0) temp = sigma[t] / image->pixdim[n + 1]; // mm to voxel
else temp = fabs(sigma[t]); // voxel-based if negative value
Expand Down Expand Up @@ -1308,18 +1308,18 @@ void reg_tools_kernelConvolution(nifti_image *image,
const int& kernelType,
const int *mask,
const bool *timePoints,
const bool *axis) {
const bool *axes) {
if (image->datatype != NIFTI_TYPE_FLOAT32 && image->datatype != NIFTI_TYPE_FLOAT64)
NR_FATAL_ERROR("The image is expected to be of floating precision type");

if (image->nt <= 0) image->nt = image->dim[4] = 1;
if (image->nu <= 0) image->nu = image->dim[5] = 1;

bool axisToSmooth[3];
if (axis == nullptr) {
// All axis are smoothed by default
axisToSmooth[0] = axisToSmooth[1] = axisToSmooth[2] = true;
} else for (int i = 0; i < 3; i++) axisToSmooth[i] = axis[i];
bool axesToSmooth[3];
if (axes == nullptr) {
// All axes are smoothed by default
axesToSmooth[0] = axesToSmooth[1] = axesToSmooth[2] = true;
} else for (int i = 0; i < 3; i++) axesToSmooth[i] = axes[i];

const int activeTimePointCount = image->nt * image->nu;
unique_ptr<bool[]> activeTimePoints{ new bool[activeTimePointCount] };
Expand All @@ -1336,7 +1336,7 @@ void reg_tools_kernelConvolution(nifti_image *image,

std::visit([&](auto&& imgDataType) {
using ImgDataType = std::decay_t<decltype(imgDataType)>;
reg_tools_kernelConvolution<ImgDataType>(image, sigma, kernelType, mask, activeTimePoints.get(), axisToSmooth);
reg_tools_kernelConvolution<ImgDataType>(image, sigma, kernelType, mask, activeTimePoints.get(), axesToSmooth);
}, NiftiImage::getFloatingDataType(image));
}
/* *************************************************************** */
Expand Down
8 changes: 4 additions & 4 deletions reg-lib/cpu/_reg_tools.h
Original file line number Diff line number Diff line change
Expand Up @@ -89,15 +89,15 @@ void reg_getRealImageSpacing(nifti_image *image,
* @param mask An integer mask over which the smoothing should occur.
* @param timePoints Boolean array to specify which time points have to be
* smoothed. The array follow the dim array of the nifti header.
* @param axis Boolean array to specify which axis have to be
* @param axes Boolean array to specify which axes have to be
* smoothed. The array follow the dim array of the nifti header.
*/
void reg_tools_kernelConvolution(nifti_image *image,
const float *sigma,
const int& kernelType,
const int *mask = nullptr,
const bool *timePoints = nullptr,
const bool *axis = nullptr);
const bool *axes = nullptr);
/* *************************************************************** */
/** @brief Smooth a label image using a Gaussian kernel
* @param image Image to be smoothed
Expand All @@ -120,13 +120,13 @@ void reg_tools_labelKernelConvolution(nifti_image *image,
* @param type The image is first smoothed using a Gaussian
* kernel of 0.7 voxel standard deviation before being downsample
* if type is set to true.
* @param axis Boolean array to specify which axis have to be
* @param axes Boolean array to specify which axes have to be
* downsampled. The array follow the dim array of the nifti header.
*/
template <class PrecisionType>
void reg_downsampleImage(nifti_image *image,
int type,
bool *axis);
bool *axes);
/* *************************************************************** */
/** @brief Returns the maximal euclidean distance from a
* deformation field image
Expand Down
14 changes: 7 additions & 7 deletions reg-lib/cuda/CudaKernelConvolution.cu
Original file line number Diff line number Diff line change
Expand Up @@ -6,15 +6,15 @@ void NiftyReg::Cuda::KernelConvolution(const nifti_image *image,
const float *sigma,
const int kernelType,
const bool *timePoints,
const bool *axis) {
const bool *axes) {
if (image->nx > 2048 || image->ny > 2048 || image->nz > 2048)
NR_FATAL_ERROR("This function does not support images with dimensions larger than 2048");

bool axisToSmooth[3];
if (axis == nullptr) {
// All axis are smoothed by default
axisToSmooth[0] = axisToSmooth[1] = axisToSmooth[2] = true;
} else for (int i = 0; i < 3; i++) axisToSmooth[i] = axis[i];
bool axesToSmooth[3];
if (axes == nullptr) {
// All axes are smoothed by default
axesToSmooth[0] = axesToSmooth[1] = axesToSmooth[2] = true;
} else for (int i = 0; i < 3; i++) axesToSmooth[i] = axes[i];

const auto activeTimePointCount = std::min(image->nt * image->nu, 4);
bool activeTimePoints[4]{}; // 4 is the maximum number of time points
Expand Down Expand Up @@ -49,7 +49,7 @@ void NiftyReg::Cuda::KernelConvolution(const nifti_image *image,

// Loop over the x, y and z dimensions
for (int n = 0; n < 3; n++) {
if (!axisToSmooth[n] || image->dim[n] <= 1) continue;
if (!axesToSmooth[n] || image->dim[n] <= 1) continue;

double temp;
if (sigma[t] > 0) temp = sigma[t] / image->pixdim[n + 1]; // mm to voxel
Expand Down
4 changes: 2 additions & 2 deletions reg-lib/cuda/CudaKernelConvolution.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,15 +13,15 @@ namespace NiftyReg::Cuda {
* @param kernelType Type of kernel to use.
* @param timePoints Boolean array to specify which time points have to be
* smoothed. The array follow the dim array of the nifti header.
* @param axis Boolean array to specify which axis have to be
* @param axes Boolean array to specify which axes have to be
* smoothed. The array follow the dim array of the nifti header.
*/
void KernelConvolution(const nifti_image *image,
float4 *imageCuda,
const float *sigma,
const int kernelType,
const bool *timePoints = nullptr,
const bool *axis = nullptr);
const bool *axes = nullptr);
/* *************************************************************** */
}
/* *************************************************************** */
1 change: 1 addition & 0 deletions reg-test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,7 @@ set(EXEC_LIST reg_test_voxelCentricToNodeCentric ${EXEC_LIST})
if(USE_CUDA)
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})
set(EXEC_LIST reg_test_regr_lts ${EXEC_LIST})
set(EXEC_LIST reg_test_regr_measure ${EXEC_LIST})
endif(USE_CUDA)
Expand Down
1 change: 1 addition & 0 deletions reg-test/reg_test_common.h
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include <array>
#include <random>
#include <iomanip>
#include <numeric>
#include <catch2/catch_test_macros.hpp>
#include "_reg_lncc.h"
#include "_reg_localTrans.h"
Expand Down
168 changes: 168 additions & 0 deletions reg-test/reg_test_regr_kernelConvolution.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,168 @@
#include "reg_test_common.h"
#include "CudaContent.h"
#include "CudaKernelConvolution.hpp"

/**
* Kernel convolution regression test to ensure the CPU and CUDA versions yield the same output
**/

class KernelConvolutionTest {
protected:
using TestData = std::tuple<std::string, NiftiImage&, vector<float>, int, bool*, bool*>;
using TestCase = std::tuple<std::string, NiftiImage, NiftiImage>;

inline static vector<TestCase> testCases;

public:
KernelConvolutionTest() {
if (!testCases.empty())
return;

// Create a random number generator
std::mt19937 gen(0);
std::uniform_real_distribution<float> distr(0, 1);

// Create images
constexpr int imageCount = 8;
constexpr NiftiImage::dim_t size = 16;
vector<NiftiImage::dim_t> dims[imageCount]{ { size, size },
{ size, size, 1, 1, 2 },
{ size, size, 1, 1, 3 },
{ size, size, 1, 2, 2 },
{ size, size, size },
{ size, size, size, 2, 1 },
{ size, size, size, 3, 1 },
{ size, size, size, 2, 2 } };
NiftiImage images[imageCount];

// Fill images with random values
for (int i = 0; i < imageCount; i++) {
images[i] = NiftiImage(dims[i], NIFTI_TYPE_FLOAT32);
auto imagePtr = images[i].data();
for (size_t j = 0; j < images[i].nVoxels(); j++)
imagePtr[j] = distr(gen);
}

// Create a lambda to concatenate strings for std::accumulate
auto strConcat = [](const std::string& str, const auto& val) { return str + " "s + std::to_string(val); };

// Create the data container for the regression test
constexpr int kernelTypeCount = 4;
distr.param(std::uniform_real_distribution<float>::param_type(1, 10)); // Change the range of the distribution
vector<TestData> testData;
for (int i = 0; i < imageCount; i++) {
for (int kernelType = 0; kernelType < kernelTypeCount; kernelType++) {
vector<float> sigmaValues(images[i]->nt * images[i]->nu);
std::generate(sigmaValues.begin(), sigmaValues.end(), [&]() { return distr(gen); });
const std::string sigmaStr = std::accumulate(sigmaValues.begin(), sigmaValues.end(), ""s, strConcat);
const std::string dimsStr = std::accumulate(dims[i].begin(), dims[i].end(), ""s, strConcat);
testData.emplace_back(TestData(
"Kernel: "s + std::to_string(kernelType) + " Sigma:"s + sigmaStr + " Dims:"s + dimsStr,
images[i],
std::move(sigmaValues),
kernelType,
nullptr,
nullptr
));
}
}

// Define time points and axes to smooth
constexpr auto timePointCount = 4;
bool timePoints[timePointCount][4]{ { true, false, false, false },
{ false, true, false, false },
{ false, false, true, false },
{ false, false, false, true } };
bool axes[timePointCount][3]{ { true, false, false },
{ false, true, false },
{ false, false, true },
{ true, true, true } };

// Add the time points and axes to the latest test data
auto latestTestData = testData.end() - timePointCount;
for (int i = 0; i < timePointCount; i++) {
auto&& [testName, image, sigmaValues, kernelType, activeTimePoints, activeAxes] = latestTestData[i];
const std::string timePointsStr = std::accumulate(timePoints[i], timePoints[i] + 4, ""s, strConcat);
const std::string axesStr = std::accumulate(axes[i], axes[i] + 3, ""s, strConcat);
testData.emplace_back(TestData(
testName + " TimePoints:"s + timePointsStr + " Axes:"s + axesStr,
image,
sigmaValues,
kernelType,
timePoints[i],
axes[i]
));
}

// Create the platforms
Platform platformCpu(PlatformType::Cpu);
Platform platformCuda(PlatformType::Cuda);

for (auto&& testData : testData) {
// Get the test data
auto&& [testName, image, sigmaValues, kernelType, activeTimePoints, activeAxes] = testData;

// Create images
NiftiImage imageCpu(image), imageCuda(image);

// Create the contents
unique_ptr<Content> contentCpu{ new Content(
imageCpu,
imageCpu,
nullptr,
nullptr,
sizeof(float)
) };
unique_ptr<CudaContent> contentCuda{ new CudaContent(
imageCuda,
imageCuda,
nullptr,
nullptr,
sizeof(float)
) };

// Use deformation fields to store images
contentCpu->SetDeformationField(imageCpu.disown());
contentCuda->SetDeformationField(imageCuda.disown());

// Compute the kernel convolution for CPU and CUDA
reg_tools_kernelConvolution(contentCpu->GetDeformationField(), sigmaValues.data(), kernelType, nullptr, activeTimePoints, activeAxes);
Cuda::KernelConvolution(contentCuda->Content::GetDeformationField(), contentCuda->GetDeformationFieldCuda(), sigmaValues.data(), kernelType, activeTimePoints, activeAxes);

// Get the images
imageCpu = NiftiImage(contentCpu->GetDeformationField(), NiftiImage::Copy::Image);
imageCuda = NiftiImage(contentCuda->GetDeformationField(), NiftiImage::Copy::Image);

// Save for testing
testCases.push_back({ testName, std::move(imageCpu), std::move(imageCuda) });
}
}
};

TEST_CASE_METHOD(KernelConvolutionTest, "Regression Kernel Convolution", "[regression]") {
// Loop over all generated test cases
for (auto&& testCase : testCases) {
// Retrieve test information
auto&& [testName, imageCpu, imageCuda] = 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 images
const auto imageCpuPtr = imageCpu.data();
const auto imageCudaPtr = imageCuda.data();
for (size_t i = 0; i < imageCpu.nVoxels(); ++i) {
const float cpuVal = imageCpuPtr[i];
const float cudaVal = imageCudaPtr[i];
if (cpuVal != cpuVal && cudaVal != cudaVal) continue; // Skip NaN values
const float diff = fabs(cpuVal - cudaVal);
if (diff > EPS)
NR_COUT << i << " " << cpuVal << " " << cudaVal << std::endl;
REQUIRE(diff < EPS);
}
}
}
}

0 comments on commit 9f24fa1

Please sign in to comment.