#148 fix mask size for 4D images

Author: Marc Modat

reg-lib/_reg_aladin.cpp

@@ -93,6 +93,11 @@ void reg_aladin<T>::InitialiseRegistration() {
 
     this->Print();
 
+    if (this->inputReference->nt > 1 || this->inputReference->nu > 1)
+        NR_WARN("The reference image contains more than one volume - only the first volume is used for the optimisation");
+    if (this->inputFloating->nt > 1 || this->inputFloating->nu > 1)
+        NR_WARN("The floating image contains more than one volume - only the first volume is used for the optimisation");
+
     // CREATE THE PYRAMID IMAGES
     this->referencePyramid = vector<NiftiImage>(this->levelsToPerform);
     this->floatingPyramid = vector<NiftiImage>(this->levelsToPerform);

reg-lib/cuda/CudaAladinContent.cpp

@@ -22,10 +22,12 @@ CudaAladinContent::CudaAladinContent(NiftiImage& referenceIn,
 /* *************************************************************** */
 void CudaAladinContent::AllocateMask() {
     if (!referenceMask) return;
+    // The host mask buffer is allocated per volume
+    const size_t voxelNumber = reference.nVoxelsPerVolume();
     int *mask;
-    Cuda::Allocate(&mask, reference->nvox);
+    Cuda::Allocate(&mask, voxelNumber);
     maskCuda.reset(mask);
-    Cuda::TransferNiftiToDevice(maskCuda.get(), referenceMask, reference->nvox);
+    Cuda::TransferNiftiToDevice(maskCuda.get(), referenceMask, voxelNumber);
 }
 /* *************************************************************** */
 void CudaAladinContent::AllocateReferenceMat() {

reg-test/reg_test_blockMatching.cpp

@@ -111,8 +111,35 @@ class BMTest {
             contentResampling3d->GetWarped(),
             mask3d.get()
         ));
+
+        // Create 4D variants of the 3D reference and warped images: volume 0 holds the same data,
+        // while the extra volumes are filled with values that must not influence the matching.
+        // The mask stays allocated per volume, mirroring the backward content of a symmetric
+        // registration with a 4D floating image
+        vector<NiftiImage::dim_t> dim4d{ size, size, size, 3 };
+        NiftiImage reference4d(dim4d, NIFTI_TYPE_FLOAT32);
+        NiftiImage warped4d(dim4d, NIFTI_TYPE_FLOAT32);
+        const size_t volumeSize = reference3d.nVoxels();
+        const auto warped3dPtr = contentResampling3d->GetWarped().data();
+        auto ref4dPtr = reference4d.data();
+        auto warped4dPtr = warped4d.data();
+        for (size_t i = 0; i < volumeSize; ++i) {
+            ref4dPtr[i] = static_cast<float>(ref3dPtr[i]);
+            warped4dPtr[i] = static_cast<float>(warped3dPtr[i]);
+        }
+        for (size_t i = volumeSize; i < reference4d.nVoxels(); ++i) {
+            ref4dPtr[i] = 1000.f * distr(gen);
+            warped4dPtr[i] = 1000.f * distr(gen);
+        }
         contentResampling3d.release();
 
+        testData.emplace_back(TestData(
+            "BlockMatching 4D",
+            std::move(reference4d),
+            std::move(warped4d),
+            mask3d.get()
+        ));
+
         for (auto&& data : testData) {
             // Get the test data
             auto&& [testName, reference, warped, mask] = data;