BUG: Fix InverseDisplacementFieldImageFilter subsampling geometry

Modules/Filtering/DisplacementField/include/itkInverseDisplacementFieldImageFilter.h

@@ -111,6 +111,7 @@ class ITK_TEMPLATE_EXPORT InverseDisplacementFieldImageFilter : public ImageToIm
   /** Image spacing type alias */
   using SpacingType = typename TOutputImage::SpacingType;
   using OriginPointType = typename TOutputImage::PointType;
+  using DirectionType = typename TOutputImage::DirectionType;
 
   /** Get/Set the coordinate transformation.
    * Set the KernelBase spline used for resampling the displacement grid.
@@ -143,6 +144,14 @@ class ITK_TEMPLATE_EXPORT InverseDisplacementFieldImageFilter : public ImageToIm
   /** Get the output image origin. */
   itkGetConstReferenceMacro(OutputOrigin, OriginPointType);
 
+  /** Set/Get the output image direction cosines. Unless explicitly set, the
+   * output inherits the direction of the input displacement field. */
+  /** @ITKStartGrouping */
+  virtual void
+  SetOutputDirection(const DirectionType & direction);
+  itkGetConstReferenceMacro(OutputDirection, DirectionType);
+  /** @ITKEndGrouping */
+
   /** Set/Get the factor used for subsampling the input displacement field.  A
    * large value in this factor will produce a fast computation of the inverse
    * field but with low precision. A small value of this factor will produce a
@@ -172,6 +181,15 @@ class ITK_TEMPLATE_EXPORT InverseDisplacementFieldImageFilter : public ImageToIm
   [[nodiscard]] ModifiedTimeType
   GetMTime() const override;
 
+  /** Verify that SubsamplingFactor is non-zero. */
+  void
+  VerifyPreconditions() const override;
+
+  /** Verify that the input is large enough to produce at least one landmark
+   * per axis after subsampling. */
+  void
+  VerifyInputInformation() const override;
+
   itkConceptMacro(OutputHasNumericTraitsCheck, (Concept::HasNumericTraits<OutputPixelComponentType>));
 
 protected:
@@ -199,6 +217,8 @@ class ITK_TEMPLATE_EXPORT InverseDisplacementFieldImageFilter : public ImageToIm
                                                   // use
   SpacingType     m_OutputSpacing{};              // output image spacing
   OriginPointType m_OutputOrigin{};               // output image origin
+  DirectionType   m_OutputDirection{};            // output image direction cosines
+  bool            m_OutputDirectionSpecified{ false };
 
   unsigned int m_SubsamplingFactor{}; // factor to subsample the
                                       // input field.

Modules/Filtering/DisplacementField/include/itkInverseDisplacementFieldImageFilter.hxx

@@ -35,6 +35,7 @@ InverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::InverseDisplacem
 {
   m_OutputSpacing.Fill(1.0);
   m_OutputOrigin.Fill(0.0);
+  m_OutputDirection.SetIdentity();
   for (unsigned int i = 0; i < ImageDimension; ++i)
   {
     m_Size[i] = 0;
@@ -61,6 +62,8 @@ InverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::PrintSelf(std::o
   os << indent << "Size:              " << m_Size << std::endl;
   os << indent << "OutputSpacing:     " << m_OutputSpacing << std::endl;
   os << indent << "OutputOrigin:      " << m_OutputOrigin << std::endl;
+  os << indent << "OutputDirection:   " << m_OutputDirection << std::endl;
+  itkPrintSelfBooleanMacro(OutputDirectionSpecified);
   os << indent << "KernelTransform:   " << m_KernelTransform.GetPointer() << std::endl;
   os << indent << "SubsamplingFactor: " << m_SubsamplingFactor << std::endl;
 }
@@ -91,6 +94,21 @@ InverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::SetOutputOrigin(
   this->SetOutputOrigin(OriginPointType(origin));
 }
 
+/**
+ * Set the output image direction.
+ */
+template <typename TInputImage, typename TOutputImage>
+void
+InverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::SetOutputDirection(const DirectionType & direction)
+{
+  if (!m_OutputDirectionSpecified || m_OutputDirection != direction)
+  {
+    m_OutputDirection = direction;
+    m_OutputDirectionSpecified = true;
+    this->Modified();
+  }
+}
+
 /**
  * Sub-sample the input displacement field and prepare the KernelBase
  * BSpline
@@ -117,13 +135,13 @@ InverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::PrepareKernelBas
   const InputImageType * inputImage = this->GetInput();
 
   resampler->SetInput(inputImage);
-  resampler->SetOutputOrigin(inputImage->GetOrigin());
   resampler->SetOutputDirection(inputImage->GetDirection());
 
   typename InputImageType::SpacingType spacing = inputImage->GetSpacing();
 
   using InputRegionType = typename InputImageType::RegionType;
   using InputSizeType = typename InputImageType::SizeType;
+  using InputIndexType = typename InputImageType::IndexType;
 
   const InputRegionType region = inputImage->GetLargestPossibleRegion();
 
@@ -135,7 +153,20 @@ InverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::PrepareKernelBas
     spacing[i] *= m_SubsamplingFactor;
   }
 
-  const InputRegionType subsampledRegion(region.GetIndex(), size);
+  // Center the landmark lattice: start at fine index (k-1)/2 past the input
+  // region start so the unsampled margins are balanced while landmarks stay
+  // exactly on input voxel centers. The lattice runs on a zero-based grid so
+  // the origin (physical location of index 0) maps to the first landmark.
+  InputIndexType latticeStart = region.GetIndex();
+  for (unsigned int i = 0; i < ImageDimension; ++i)
+  {
+    latticeStart[i] += static_cast<IndexValueType>((m_SubsamplingFactor - 1) / 2);
+  }
+  typename InputImageType::PointType latticeOrigin;
+  inputImage->TransformIndexToPhysicalPoint(latticeStart, latticeOrigin);
+  resampler->SetOutputOrigin(latticeOrigin);
+
+  const InputRegionType subsampledRegion(InputIndexType{}, size);
 
   resampler->SetSize(size);
   resampler->SetOutputStartIndex(subsampledRegion.GetIndex());
@@ -302,6 +333,39 @@ InverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::GenerateOutputIn
   // Set spacing and origin
   outputPtr->SetSpacing(m_OutputSpacing);
   outputPtr->SetOrigin(m_OutputOrigin);
+  if (m_OutputDirectionSpecified)
+  {
+    outputPtr->SetDirection(m_OutputDirection);
+  }
+}
+
+template <typename TInputImage, typename TOutputImage>
+void
+InverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::VerifyPreconditions() const
+{
+  Superclass::VerifyPreconditions();
+
+  if (m_SubsamplingFactor == 0)
+  {
+    itkExceptionMacro("SubsamplingFactor must be non-zero.");
+  }
+}
+
+template <typename TInputImage, typename TOutputImage>
+void
+InverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::VerifyInputInformation() const
+{
+  Superclass::VerifyInputInformation();
+
+  const typename InputImageType::SizeType inputSize = this->GetInput()->GetLargestPossibleRegion().GetSize();
+  for (unsigned int i = 0; i < ImageDimension; ++i)
+  {
+    if (inputSize[i] < m_SubsamplingFactor)
+    {
+      itkExceptionMacro("Input size " << inputSize << " must be at least SubsamplingFactor (" << m_SubsamplingFactor
+                                      << ") along every axis to produce landmarks.");
+    }
+  }
 }
 
 /**

Modules/Filtering/DisplacementField/test/itkInverseDisplacementFieldImageFilterGTest.cxx

@@ -19,13 +19,28 @@
 #include "itkInverseDisplacementFieldImageFilter.h"
 #include "itkThinPlateSplineKernelTransform.h"
 #include "itkImageRegionConstIterator.h"
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkImageRegionIteratorWithIndex.h"
+#include "itkMath.h"
 #include "itkGTest.h"
+#include <algorithm> // For max.
 
 namespace
 {
 using VectorType = itk::Vector<float, 2>;
 using FieldType = itk::Image<VectorType, 2>;
 using FilterType = itk::InverseDisplacementFieldImageFilter<FieldType, FieldType>;
+
+FieldType::DirectionType
+MakeRotation()
+{
+  FieldType::DirectionType direction;
+  direction(0, 0) = 0.0;
+  direction(0, 1) = -1.0;
+  direction(1, 0) = 1.0;
+  direction(1, 1) = 0.0;
+  return direction;
+}
 } // namespace
 
 // The kernel-spline subsampler must inherit the input direction; the exact inverse of a
@@ -35,12 +50,7 @@ TEST(InverseDisplacementFieldImageFilter, RecoversConstantInverseOnRotatedField)
 {
   auto field = FieldType::New();
   field->SetRegions(FieldType::RegionType{ FieldType::IndexType{}, FieldType::SizeType::Filled(32) });
-  FieldType::DirectionType direction;
-  direction(0, 0) = 0.0;
-  direction(0, 1) = -1.0;
-  direction(1, 0) = 1.0;
-  direction(1, 1) = 0.0;
-  field->SetDirection(direction);
+  field->SetDirection(MakeRotation());
   field->Allocate();
   field->FillBuffer(itk::MakeVector(2.0F, 3.0F));
 
@@ -62,3 +72,167 @@ TEST(InverseDisplacementFieldImageFilter, RecoversConstantInverseOnRotatedField)
     EXPECT_NEAR(it.Get()[1], -3.0, 0.05);
   }
 }
+
+// Unless SetOutputDirection is called, the output inherits the input direction.
+TEST(InverseDisplacementFieldImageFilter, OutputDirectionDefaultsToInputDirection)
+{
+  auto field = FieldType::New();
+  field->SetRegions(FieldType::RegionType{ FieldType::IndexType{}, FieldType::SizeType::Filled(32) });
+  field->SetDirection(MakeRotation());
+  field->Allocate();
+  field->FillBuffer(itk::MakeVector(2.0F, 3.0F));
+
+  auto filter = FilterType::New();
+  filter->SetInput(field);
+  filter->SetKernelTransform(itk::ThinPlateSplineKernelTransform<double, 2>::New());
+  filter->SetSubsamplingFactor(4);
+  filter->SetSize(FieldType::SizeType::Filled(16));
+  filter->SetOutputSpacing(itk::MakeFilled<FieldType::SpacingType>(1.0));
+  filter->SetOutputOrigin(itk::MakePoint(-16.0, 8.0));
+  filter->UpdateOutputInformation();
+
+  EXPECT_EQ(filter->GetOutput()->GetDirection(), MakeRotation());
+}
+
+// An explicitly requested output direction is honored and the inverse is evaluated
+// consistently on that grid: the inverse of a constant displacement is its negation.
+TEST(InverseDisplacementFieldImageFilter, HonorsExplicitOutputDirection)
+{
+  auto field = FieldType::New();
+  field->SetRegions(FieldType::RegionType{ FieldType::IndexType{}, FieldType::SizeType::Filled(32) });
+  field->SetDirection(MakeRotation());
+  field->Allocate();
+  field->FillBuffer(itk::MakeVector(2.0F, 3.0F));
+
+  auto filter = FilterType::New();
+  filter->SetInput(field);
+  filter->SetKernelTransform(itk::ThinPlateSplineKernelTransform<double, 2>::New());
+  filter->SetSubsamplingFactor(4);
+  filter->SetSize(FieldType::SizeType::Filled(16));
+  filter->SetOutputSpacing(itk::MakeFilled<FieldType::SpacingType>(1.0));
+  // Axis-aligned output window inside the subsampled landmark hull.
+  filter->SetOutputOrigin(itk::MakePoint(-16.0, 8.0));
+  FieldType::DirectionType identity;
+  identity.SetIdentity();
+  filter->SetOutputDirection(identity);
+  filter->UpdateLargestPossibleRegion();
+
+  const FieldType * output = filter->GetOutput();
+  EXPECT_EQ(output->GetDirection(), identity);
+  itk::ImageRegionConstIterator<FieldType> it(output, output->GetBufferedRegion());
+  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
+  {
+    EXPECT_NEAR(it.Get()[0], -2.0, 0.05);
+    EXPECT_NEAR(it.Get()[1], -3.0, 0.05);
+  }
+}
+
+// The landmark lattice starts at fine index (k-1)/2 so unsampled margins are balanced;
+// for a nonlinear field this bounds the round-trip error near the high edge, where the
+// uncentered lattice (samples at 0, k, ..., N-k) had to extrapolate over k-1 voxels.
+TEST(InverseDisplacementFieldImageFilter, CenteredSubsamplingBoundsHighEdgeError)
+{
+  auto field = FieldType::New();
+  field->SetRegions(FieldType::RegionType{ FieldType::IndexType{}, FieldType::SizeType::Filled(32) });
+  field->Allocate();
+  for (itk::ImageRegionIteratorWithIndex<FieldType> it(field, field->GetBufferedRegion()); !it.IsAtEnd(); ++it)
+  {
+    const auto index = it.GetIndex();
+    it.Set(itk::MakeVector(0.002F * index[0] * index[0], 0.002F * index[1] * index[1]));
+  }
+
+  auto filter = FilterType::New();
+  filter->SetInput(field);
+  filter->SetKernelTransform(itk::ThinPlateSplineKernelTransform<double, 2>::New());
+  filter->SetSubsamplingFactor(4);
+  filter->SetSize(FieldType::SizeType::Filled(32));
+  filter->SetOutputSpacing(itk::MakeFilled<FieldType::SpacingType>(1.0));
+  filter->SetOutputOrigin(itk::MakePoint(0.0, 0.0));
+  filter->UpdateLargestPossibleRegion();
+
+  const FieldType * output = filter->GetOutput();
+  double            maxError = 0.0;
+  for (itk::ImageRegionConstIteratorWithIndex<FieldType> it(output, output->GetBufferedRegion()); !it.IsAtEnd(); ++it)
+  {
+    FieldType::PointType point;
+    output->TransformIndexToPhysicalPoint(it.GetIndex(), point);
+    // Round trip: apply the inverse displacement, then the analytic forward field.
+    itk::Point<double, 2> back;
+    for (unsigned int i = 0; i < 2; ++i)
+    {
+      back[i] = point[i] + it.Get()[i];
+      back[i] += 0.002 * back[i] * back[i];
+      maxError = std::max(maxError, itk::Math::abs(back[i] - point[i]));
+    }
+  }
+  EXPECT_LT(maxError, 0.06);
+}
+
+// An input smaller than the subsampling factor yields a zero-size landmark lattice; the
+// filter must reject it up front instead of building a spline with no landmarks.
+TEST(InverseDisplacementFieldImageFilter, RejectsInputSmallerThanSubsamplingFactor)
+{
+  auto field = FieldType::New();
+  field->SetRegions(FieldType::RegionType{ FieldType::IndexType{}, FieldType::SizeType::Filled(8) });
+  field->Allocate();
+  field->FillBuffer(itk::MakeVector(0.0F, 0.0F));
+
+  auto filter = FilterType::New();
+  filter->SetInput(field);
+  filter->SetKernelTransform(itk::ThinPlateSplineKernelTransform<double, 2>::New());
+  filter->SetSubsamplingFactor(16);
+  filter->SetSize(FieldType::SizeType::Filled(8));
+  filter->SetOutputSpacing(itk::MakeFilled<FieldType::SpacingType>(1.0));
+  filter->SetOutputOrigin(itk::MakePoint(0.0, 0.0));
+  EXPECT_THROW(filter->UpdateLargestPossibleRegion(), itk::ExceptionObject);
+
+  filter->SetSubsamplingFactor(0);
+  EXPECT_THROW(filter->UpdateLargestPossibleRegion(), itk::ExceptionObject);
+
+  filter->SetSubsamplingFactor(4);
+  EXPECT_NO_THROW(filter->UpdateLargestPossibleRegion());
+}
+
+// The centered lattice must account for a non-zero input region start index: the
+// landmarks are placed relative to the region start, not the absolute index 0, so a
+// cropped/non-zero-start input inverts as accurately as a zero-start one (issue #6582
+// follow-up; greptile P1). The input below occupies physical [0,31] via index [5,36]
+// with a compensating origin, matching CenteredSubsamplingBoundsHighEdgeError's grid.
+TEST(InverseDisplacementFieldImageFilter, NonZeroStartRegionInvertsCorrectly)
+{
+  auto field = FieldType::New();
+  field->SetRegions(FieldType::RegionType{ itk::MakeIndex(5, 5), FieldType::SizeType::Filled(32) });
+  field->SetOrigin(itk::MakePoint(-5.0, -5.0));
+  field->Allocate();
+  for (itk::ImageRegionIteratorWithIndex<FieldType> it(field, field->GetBufferedRegion()); !it.IsAtEnd(); ++it)
+  {
+    FieldType::PointType p;
+    field->TransformIndexToPhysicalPoint(it.GetIndex(), p);
+    it.Set(itk::MakeVector(static_cast<float>(0.002 * p[0] * p[0]), static_cast<float>(0.002 * p[1] * p[1])));
+  }
+
+  auto filter = FilterType::New();
+  filter->SetInput(field);
+  filter->SetKernelTransform(itk::ThinPlateSplineKernelTransform<double, 2>::New());
+  filter->SetSubsamplingFactor(4);
+  filter->SetSize(FieldType::SizeType::Filled(32));
+  filter->SetOutputSpacing(itk::MakeFilled<FieldType::SpacingType>(1.0));
+  filter->SetOutputOrigin(itk::MakePoint(0.0, 0.0));
+  filter->UpdateLargestPossibleRegion();
+
+  const FieldType * output = filter->GetOutput();
+  double            maxError = 0.0;
+  for (itk::ImageRegionConstIteratorWithIndex<FieldType> it(output, output->GetBufferedRegion()); !it.IsAtEnd(); ++it)
+  {
+    FieldType::PointType point;
+    output->TransformIndexToPhysicalPoint(it.GetIndex(), point);
+    itk::Point<double, 2> back;
+    for (unsigned int i = 0; i < 2; ++i)
+    {
+      back[i] = point[i] + it.Get()[i];
+      back[i] += 0.002 * back[i] * back[i];
+      maxError = std::max(maxError, itk::Math::abs(back[i] - point[i]));
+    }
+  }
+  EXPECT_LT(maxError, 0.06);
+}