mrtransform: add -oversample option

cmd/mrtransform.cpp

@@ -123,6 +123,14 @@ void usage ()
         "set the interpolation method to use when reslicing (choices: nearest, linear, cubic, sinc. Default: cubic).")
     + Argument ("method").type_choice (interp_choices)
 
+    + Option ("oversample",
+        "set the amount of over-sampling (in the target space) to perform when regridding. This is particularly "
+        "relevant when downsamping a high-resolution image to a low-resolution image, to avoid aliasing artefacts. "
+        "This can consist of a single integer, or a comma-separated list of 3 integers if different oversampling "
+        "factors are desired along the different axes. Default is determined from ratio of voxel dimensions (disabled "
+        "for nearest-neighbour interpolation).")
+    + Argument ("factor").type_sequence_int()
+
     + OptionGroup ("Non-linear transformation options")
 
     // TODO point users to a documentation page describing the warp field format
@@ -170,19 +178,20 @@ void usage ()
       "Use NaN as the out of bounds value (Default: 0.0)");
 }
 
-void apply_warp (Image<float>& input, Image<float>& output, Image<default_type>& warp, const int interp, const float out_of_bounds_value) {
+void apply_warp (Image<float>& input, Image<float>& output, Image<default_type>& warp,
+  const int interp, const float out_of_bounds_value, const vector<int>& oversample) {
   switch (interp) {
   case 0:
-    Filter::warp<Interp::Nearest> (input, output, warp, out_of_bounds_value);
+    Filter::warp<Interp::Nearest> (input, output, warp, out_of_bounds_value, oversample);
     break;
   case 1:
-    Filter::warp<Interp::Linear> (input, output, warp, out_of_bounds_value);
+    Filter::warp<Interp::Linear> (input, output, warp, out_of_bounds_value, oversample);
     break;
   case 2:
-    Filter::warp<Interp::Cubic> (input, output, warp, out_of_bounds_value);
+    Filter::warp<Interp::Cubic> (input, output, warp, out_of_bounds_value, oversample);
     break;
   case 3:
-    Filter::warp<Interp::Sinc> (input, output, warp, out_of_bounds_value);
+    Filter::warp<Interp::Sinc> (input, output, warp, out_of_bounds_value, oversample);
     break;
   default:
     assert (0);
@@ -441,6 +450,27 @@ void run ()
       WARN ("interpolator choice ignored since the input image will not be regridded");
   }
 
+  // over-sampling
+  vector<int> oversample = Adapter::AutoOverSample;
+  opt = get_options ("oversample");
+  if (opt.size()) {
+    if (!template_header.valid() && !warp)
+      throw Exception ("-oversample option applies only to regridding using the template option or to non-linear transformations");
+    oversample = opt[0][0];
+    if (oversample.size() == 1)
+      oversample.resize (3, oversample[0]);
+    else if (oversample.size() != 3)
+      throw Exception ("-oversample option requires either a single integer, or a comma-separated list of 3 integers");
+    for (const auto x : oversample)
+      if (x < 1)
+        throw Exception ("-oversample factors must be positive integers");
+  }
+  else if (interp == 0)
+    // default for nearest-neighbour is no oversampling
+    oversample = { 1, 1, 1 };
+
+
+
   // Out of bounds value
   float out_of_bounds_value = 0.0;
   opt = get_options ("nan");
@@ -478,16 +508,16 @@ void run ()
 
     switch (interp) {
       case 0:
-        Filter::reslice<Interp::Nearest> (input, output, linear_transform, Adapter::AutoOverSample, out_of_bounds_value);
+        Filter::reslice<Interp::Nearest> (input, output, linear_transform, oversample, out_of_bounds_value);
         break;
       case 1:
-        Filter::reslice<Interp::Linear> (input, output, linear_transform, Adapter::AutoOverSample, out_of_bounds_value);
+        Filter::reslice<Interp::Linear> (input, output, linear_transform, oversample, out_of_bounds_value);
         break;
       case 2:
-        Filter::reslice<Interp::Cubic> (input, output, linear_transform, Adapter::AutoOverSample, out_of_bounds_value);
+        Filter::reslice<Interp::Cubic> (input, output, linear_transform, oversample, out_of_bounds_value);
         break;
       case 3:
-        Filter::reslice<Interp::Sinc> (input, output, linear_transform, Adapter::AutoOverSample, out_of_bounds_value);
+        Filter::reslice<Interp::Sinc> (input, output, linear_transform, oversample, out_of_bounds_value);
         break;
       default:
         assert (0);
@@ -523,21 +553,21 @@ void run ()
       } else {
         warp_deform = Registration::Warp::compute_full_deformation (warp, template_header, from);
       }
-      apply_warp (input, output, warp_deform, interp, out_of_bounds_value);
+      apply_warp (input, output, warp_deform, interp, out_of_bounds_value, oversample);
       if (fod_reorientation)
         Registration::Transform::reorient_warp ("reorienting", output, warp_deform, directions_cartesian.transpose(), modulate);
 
     // Compose and apply input linear and 4D deformation field
     } else if (warp.ndim() == 4 && linear) {
       auto warp_composed = Image<default_type>::scratch (warp);
       Registration::Warp::compose_linear_deformation (linear_transform, warp, warp_composed);
-      apply_warp (input, output, warp_composed, interp, out_of_bounds_value);
+      apply_warp (input, output, warp_composed, interp, out_of_bounds_value, oversample);
       if (fod_reorientation)
         Registration::Transform::reorient_warp ("reorienting", output, warp_composed, directions_cartesian.transpose(), modulate);
 
     // Apply 4D deformation field only
     } else {
-      apply_warp (input, output, warp, interp, out_of_bounds_value);
+      apply_warp (input, output, warp, interp, out_of_bounds_value, oversample);
       if (fod_reorientation)
         Registration::Transform::reorient_warp ("reorienting", output, warp, directions_cartesian.transpose(), modulate);
     }

core/filter/resize.h

@@ -131,7 +131,8 @@ namespace MR
           {
             switch (interp_type) {
             case 0:
-              reslice <Interp::Nearest> (input, output);
+              // Prevent use of oversampling when using nearest-neighbour interpolation
+              reslice <Interp::Nearest> (input, output, Adapter::NoTransform, { 1, 1, 1 });
               break;
             case 1:
               reslice <Interp::Linear> (input, output);

core/filter/warp.h

@@ -63,7 +63,8 @@ namespace MR
           ImageTypeSource& source,
           ImageTypeDestination& destination,
           WarpType& warp,
-          const typename ImageTypeDestination::value_type value_when_out_of_bounds = Interpolator<ImageTypeSource>::default_out_of_bounds_value())
+          const typename ImageTypeDestination::value_type value_when_out_of_bounds = Interpolator<ImageTypeSource>::default_out_of_bounds_value(),
+          vector<int> oversample = Adapter::AutoOverSample)
       {
 
         // reslice warp onto destination grid
@@ -76,7 +77,8 @@ namespace MR
            header.size(3) = 3;
            Stride::set (header, Stride::contiguous_along_axis (3));
            auto warp_resliced = Image<typename WarpType::value_type>::scratch (header);
-           reslice<Interp::Cubic> (warp, warp_resliced);
+           transform_type identity_transform;
+           reslice<Interp::Cubic> (warp, warp_resliced, identity_transform, oversample);
            Adapter::Warp<Interpolator, ImageTypeSource, Image<typename WarpType::value_type> > interp (source, warp_resliced, value_when_out_of_bounds);
 
            if (destination.ndim() == 4)