Skip to content

Commit

Permalink
#92: block match & LTS CPU/CUDA regression tests
Browse files Browse the repository at this point in the history
  • Loading branch information
mmodat committed Jul 14, 2023
1 parent 46bb6c8 commit c425883
Show file tree
Hide file tree
Showing 8 changed files with 408 additions and 19 deletions.
2 changes: 1 addition & 1 deletion niftyreg_build_version.txt
Original file line number Diff line number Diff line change
@@ -1 +1 @@
271
272
32 changes: 31 additions & 1 deletion reg-lib/cpu/_reg_blockMatching.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,36 @@
#include <map>
#include <iostream>
#include <cmath>

_reg_blockMatchingParam::_reg_blockMatchingParam(_reg_blockMatchingParam *in)
{
this->totalBlockNumber=in->totalBlockNumber;
this->dim=in->dim;
this->percent_to_keep=in->percent_to_keep;
this->activeBlockNumber=in->activeBlockNumber;
this->definedActiveBlockNumber=in->definedActiveBlockNumber;
this->stepSize=in->stepSize;
this->voxelCaptureRange=in->voxelCaptureRange;
this->blockNumber[0]=in->blockNumber[0];
this->blockNumber[1]=in->blockNumber[1];
this->blockNumber[2]=in->blockNumber[2];
this->totalBlock = (int *)malloc(this->totalBlockNumber * sizeof(int));
for(int i=0; i<this->totalBlockNumber; ++i)
this->totalBlock[i] = in->totalBlock[i];

this->referencePosition = (float *)malloc(this->activeBlockNumber * this->dim * sizeof(float));
this->warpedPosition = (float *)malloc(this->activeBlockNumber * this->dim * sizeof(float));
for(int i=0; i<this->activeBlockNumber*this->dim ; ++i){
this->referencePosition[i] = in->referencePosition[i];
this->warpedPosition[i] = in->warpedPosition[i];
}
}
_reg_blockMatchingParam::~_reg_blockMatchingParam()
{
if (referencePosition) free(referencePosition);
if (warpedPosition) free(warpedPosition);
if (totalBlock) free(totalBlock);
}
/* *************************************************************** */
template<class DataType>
void _reg_set_active_blocks(nifti_image *referenceImage, _reg_blockMatchingParam *params, int *mask, bool runningOnGPU) {
Expand Down Expand Up @@ -247,7 +277,7 @@ void initialise_block_matching_method(nifti_image * reference,
reg_print_msg_debug(text)
#endif
//params->activeBlock = (int *)malloc(params->activeBlockNumber * sizeof(int));
params->referencePosition = (float *)malloc(params->activeBlockNumber * params->dim * sizeof(float));
params->referencePosition = (float *)malloc(params->activeBlockNumber * params->dim * sizeof(float));
params->warpedPosition = (float *)malloc(params->activeBlockNumber * params->dim * sizeof(float));

#ifndef NDEBUG
Expand Down
10 changes: 4 additions & 6 deletions reg-lib/cpu/_reg_blockMatching.h
Original file line number Diff line number Diff line change
Expand Up @@ -69,12 +69,10 @@ struct _reg_blockMatchingParam
stepSize(0)
{}

~_reg_blockMatchingParam()
{
if (referencePosition) free(referencePosition);
if (warpedPosition) free(warpedPosition);
if (totalBlock) free(totalBlock);
}
// Perform a deep copy
_reg_blockMatchingParam(_reg_blockMatchingParam *);

~_reg_blockMatchingParam();
};
/* *************************************************************** */
/** @brief This function initialise a _reg_blockMatchingParam structure
Expand Down
22 changes: 12 additions & 10 deletions reg-lib/cuda/CudaAladinContent.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ CudaAladinContent::~CudaAladinContent() {
void CudaAladinContent::InitVars() {
referenceImageArray_d = nullptr;
floatingImageArray_d = nullptr;
transformationMatrix_d = nullptr;
warpedImageArray_d = nullptr;
deformationFieldArray_d = nullptr;
referencePosition_d = nullptr;
Expand Down Expand Up @@ -404,31 +405,33 @@ int* CudaAladinContent::GetFloatingDims() {
}
/* *************************************************************** */
void CudaAladinContent::FreeCuPtrs() {
if (transformationMatrix != nullptr)
if (transformationMatrix_d != nullptr)
cudaCommon_free(transformationMatrix_d);

if (reference != nullptr) {
if (referenceImageArray_d != nullptr)
cudaCommon_free(referenceImageArray_d);
if (referenceMat_d != nullptr)
cudaCommon_free(referenceMat_d);
}

if (floating != nullptr) {
if (floatingImageArray_d != nullptr)
cudaCommon_free(floatingImageArray_d);
if (floIJKMat_d != nullptr)
cudaCommon_free(floIJKMat_d);
}

if (warped != nullptr)
if (warpedImageArray_d != nullptr)
cudaCommon_free(warpedImageArray_d);

if (deformationField != nullptr)
if (deformationFieldArray_d != nullptr)
cudaCommon_free(deformationFieldArray_d);

if (referenceMask != nullptr)
if (mask_d != nullptr)
cudaCommon_free(mask_d);

if (blockMatchingParams != nullptr) {
if (totalBlock_d != nullptr)
cudaCommon_free(totalBlock_d);
if (referencePosition_d != nullptr)
cudaCommon_free(referencePosition_d);
if (warpedPosition_d != nullptr)
cudaCommon_free(warpedPosition_d);
/*
cudaCommon_free(AR_d);
Expand All @@ -438,7 +441,6 @@ void CudaAladinContent::FreeCuPtrs() {
cudaCommon_free(lengths_d);
cudaCommon_free(newWarpedPos_d);
*/
}
}
/* *************************************************************** */
bool CudaAladinContent::IsCurrentComputationDoubleCapable() {
Expand Down
3 changes: 2 additions & 1 deletion reg-lib/cuda/_reg_common_cuda.cu
Original file line number Diff line number Diff line change
Expand Up @@ -577,7 +577,8 @@ void cudaCommon_free(cudaArray *cuArray_d) {
/* *************************************************************** */
template <class DataType>
void cudaCommon_free(DataType *array_d) {
NR_CUDA_SAFE_CALL(cudaFree(array_d));
if (array_d != nullptr)
NR_CUDA_SAFE_CALL(cudaFree(array_d));
}
template void cudaCommon_free<int>(int*);
template void cudaCommon_free<float>(float*);
Expand Down
5 changes: 5 additions & 0 deletions reg-test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,11 @@ set(EXEC_LIST reg_test_interpolation ${EXEC_LIST})
set(EXEC_LIST reg_test_lncc ${EXEC_LIST})
set(EXEC_LIST reg_test_normaliseGradient ${EXEC_LIST})
set(EXEC_LIST reg_test_voxelCentricToNodeCentric ${EXEC_LIST})
if(USE_CUDA)
set(EXEC_LIST reg_test_regr_blockMatching ${EXEC_LIST})
set(EXEC_LIST reg_test_regr_lts ${EXEC_LIST})
endif(USE_CUDA)


foreach(EXEC ${EXEC_LIST})
add_executable(${EXEC} ${EXEC}.cpp)
Expand Down
175 changes: 175 additions & 0 deletions reg-test/reg_test_regr_blockMatching.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,175 @@
#include "reg_test_common.h"
#include "_reg_blockMatching.h"
#include "CpuBlockMatchingKernel.h"
#include "CudaBlockMatchingKernel.h"

/*
This test file contains a regression test to ensure the CPU and GPU version yield the same output
*/

class BMTest {
/*
Class to call the block matching function
*/
protected:
using TestData = std::tuple<std::string, NiftiImage, NiftiImage, int>;
using TestCase = std::tuple<std::string, _reg_blockMatchingParam *, _reg_blockMatchingParam *>;
inline static vector<TestCase> testCases;
NiftiImage reference2d;
NiftiImage floating2d;
NiftiImage reference3d;
NiftiImage floating3d;
public:
~BMTest() {
std::cout << "Calling destructor" << std::endl;
}
BMTest() {
std::cout << "Calling constructor" << std::endl;
if (!testCases.empty())
return;

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

// Create a reference and floating 2D images
NiftiImage::dim_t size = 64;
vector<NiftiImage::dim_t> dim{ size, size };
this->reference2d = NiftiImage(dim, NIFTI_TYPE_FLOAT32);
this->floating2d = NiftiImage(dim, NIFTI_TYPE_FLOAT32);

// Create a reference 3D image
dim.push_back(size);
this->reference3d = NiftiImage(dim, NIFTI_TYPE_FLOAT32);
this->floating3d = NiftiImage(dim, NIFTI_TYPE_FLOAT32);

// Fill images with random values
float *ref2dPtr = static_cast<float *>(reference2d->data);
float *flo2dPtr = static_cast<float *>(floating2d->data);
for (int y = 0; y < reference2d->ny; ++y)
for (int x = 0; x < reference2d->nx; ++x) {
*ref2dPtr++ = distr(gen);
*flo2dPtr++ = distr(gen);
}

// Fill images with random values
float *ref3dPtr = static_cast<float *>(reference3d->data);
float *flo3dPtr = static_cast<float *>(floating3d->data);
for (int z = 0; z < reference3d->nz; ++z)
for (int y = 0; y < reference3d->ny; ++y)
for (int x = 0; x < reference3d->nx; ++x) {
*ref3dPtr++ = distr(gen);
*flo3dPtr++ = distr(gen);
}


// Create the data container for the regression test
vector<TestData> testData;
for(int b=50; b<=100; b+=50){
testData.emplace_back(TestData(
"BlockMatching 2D block " + std::to_string(b),
std::move(NiftiImage(this->reference2d)),
std::move(NiftiImage(this->floating2d)),
b
));
testData.emplace_back(TestData(
"BlockMatching 3D block " + std::to_string(b),
std::move(NiftiImage(this->reference3d)),
std::move(NiftiImage(this->floating3d)),
b
));
}

for (auto&& data : testData) {
unique_ptr<Platform> platformCPU{ new Platform(PlatformType::Cpu) };
unique_ptr<Platform> platformCUDA{ new Platform(PlatformType::Cuda) };
// Make a copy of the test data
auto&& [testName, reference, floating, block] = data;
// Create content creator
unique_ptr<AladinContentCreator> contentCreatorCPU{
dynamic_cast<AladinContentCreator*>(platformCPU->CreateContentCreator(ContentType::Aladin))
};
unique_ptr<AladinContentCreator> contentCreatorCUDA{
dynamic_cast<AladinContentCreator*>(platformCUDA->CreateContentCreator(ContentType::Aladin))
};
// Create the contents
unique_ptr<AladinContent> contentCPU{ contentCreatorCPU->Create(
NiftiImage(reference).disown(),
NiftiImage(floating).disown(),
nullptr,
nullptr,
sizeof(float),
100,
block,
1
)};
unique_ptr<AladinContent> contentCUDA{ contentCreatorCUDA->Create(
NiftiImage(reference).disown(),
NiftiImage(floating).disown(),
nullptr,
nullptr,
sizeof(float),
100,
block,
1
)};
// Initialise the warped image
contentCPU->SetWarped(NiftiImage(floating).disown());
contentCUDA->SetWarped(NiftiImage(floating).disown());
// Initialise the block matching
std::unique_ptr<Kernel> kernelCPU = nullptr;
kernelCPU.reset(platformCPU->CreateKernel(BlockMatchingKernel::GetName(), contentCPU.get()));
std::unique_ptr<Kernel> kernelCUDA = nullptr;
kernelCUDA.reset(platformCUDA->CreateKernel(BlockMatchingKernel::GetName(), contentCUDA.get()));

// run the computation
kernelCPU->template castTo<CpuBlockMatchingKernel>()->Calculate();
kernelCUDA->template castTo<CudaBlockMatchingKernel>()->Calculate();

// Retrieve the information
_reg_blockMatchingParam *blockMatchingParamsCPU = new _reg_blockMatchingParam(contentCPU->GetBlockMatchingParams());
_reg_blockMatchingParam *blockMatchingParamsCUDA = new _reg_blockMatchingParam(contentCUDA->GetBlockMatchingParams());

testCases.push_back({
testName,
blockMatchingParamsCPU,
blockMatchingParamsCUDA
});
contentCPU.reset();
contentCUDA.reset();
}
}
};

TEST_CASE_METHOD(BMTest, "Regression BlockMatching", "[regression]") {
// Loop over all generated test cases
for (auto&& testCase : this->testCases) {
// Retrieve test information
auto&& [testName, blockMatchingParamsCPU, blockMatchingParamsCUDA] = testCase;

SECTION(testName) {

// Ensure both approaches retreive the same number of voxel
REQUIRE(blockMatchingParamsCPU->activeBlockNumber==blockMatchingParamsCUDA->activeBlockNumber);

// Loop over the block and ensure all values are identical
for(int b=0; b<blockMatchingParamsCPU->activeBlockNumber*blockMatchingParamsCPU->dim; ++b){
float delta = blockMatchingParamsCPU->referencePosition[b] - blockMatchingParamsCUDA->referencePosition[b];
if(fabs(delta) > EPS){
std::cout << "HERE " << delta << std::endl;
std::cout.flush();
}
REQUIRE(fabs(delta) < EPS);
delta = blockMatchingParamsCPU->warpedPosition[b] - blockMatchingParamsCUDA->warpedPosition[b];
if(fabs(delta) > EPS){
std::cout << "HERE " << delta << std::endl;
std::cout.flush();
}
REQUIRE(fabs(delta) < EPS);
}
delete blockMatchingParamsCPU;
delete blockMatchingParamsCUDA;
}
}
}
Loading

0 comments on commit c425883

Please sign in to comment.