ENH: Ingest ITKHigherOrderAccurateGradient into Modules/Filtering

Modules/Filtering/HigherOrderAccurateGradient/CMakeLists.txt

@@ -0,0 +1,3 @@
+project(HigherOrderAccurateGradient)
+
+itk_module_impl()

Modules/Filtering/HigherOrderAccurateGradient/README.md

@@ -0,0 +1,25 @@
+# HigherOrderAccurateGradient
+
+In-tree ITK module providing higher-order accurate numerical
+derivative and gradient image filters. The flagship classes
+`itk::HigherOrderAccurateDerivativeImageFilter`,
+`itk::HigherOrderAccurateDerivativeOperator`, and
+`itk::HigherOrderAccurateGradientImageFilter` compute centered finite
+differences of arbitrary even order on N-dimensional scalar images.
+
+## Origin
+
+Ingested from the standalone remote module
+[**InsightSoftwareConsortium/ITKHigherOrderAccurateGradient**](https://github.com/InsightSoftwareConsortium/ITKHigherOrderAccurateGradient)
+on 2026-05-08, following the v4 ingestion guidelines defined in
+InsightSoftwareConsortium/ITK#6204. The upstream repository will be
+archived read-only after this PR merges; it remains reachable at the
+URL above for historical reference (notably the `doc/` directory,
+which was intentionally left in the upstream archive).
+
+## References
+
+- McCormick M.M., Liu X., Jomier J., Marion C., Ibanez L.
+  *Higher Order Accurate Derivative and Gradient Calculation in ITK.*
+  The Insight Journal. January-December. 2014.
+  <https://hdl.handle.net/10380/3231>

Modules/Filtering/HigherOrderAccurateGradient/include/itkHigherOrderAccurateDerivativeImageFilter.h

@@ -0,0 +1,161 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         https://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkHigherOrderAccurateDerivativeImageFilter_h
+#define itkHigherOrderAccurateDerivativeImageFilter_h
+
+#include "itkImageToImageFilter.h"
+
+namespace itk
+{
+
+/** \class HigherOrderAccurateDerivativeImageFilter
+ * \brief Computes the higher order accurate directional derivative of an image.
+ * The directional derivative at each pixel location is computed by convolution
+ * with a higher order accurate derivative operator of user-specified order.
+ *
+ * SetOrder() specifies the order of the derivative.
+ *
+ * SetDirection() specifies the direction of the derivative with respect to the
+ * coordinate axes of the image.
+ *
+ * To specify the order of accuracy, use SetOrderOfAccuracy().  The
+ * approximation will be accurate to two times the OrderOfAccuracy in terms of
+ * Taylor series terms.
+ *
+ * \sa Image
+ * \sa Neighborhood
+ * \sa NeighborhoodOperator
+ * \sa NeighborhoodIterator
+ *
+ * \ingroup ImageFeatureExtraction
+ * \ingroup HigherOrderAccurateGradient
+ */
+template <typename TInputImage, typename TOutputImage>
+class HigherOrderAccurateDerivativeImageFilter : public ImageToImageFilter<TInputImage, TOutputImage>
+{
+public:
+  ITK_DISALLOW_COPY_AND_MOVE(HigherOrderAccurateDerivativeImageFilter);
+
+  /** Standard class type alias. */
+  using Self = HigherOrderAccurateDerivativeImageFilter;
+  using Superclass = ImageToImageFilter<TInputImage, TOutputImage>;
+  using Pointer = SmartPointer<Self>;
+  using ConstPointer = SmartPointer<const Self>;
+
+  /** Extract some information from the image types.  Dimensionality
+   * of the two images is assumed to be the same. */
+  using OutputPixelType = typename TOutputImage::PixelType;
+  using OutputInternalPixelType = typename TOutputImage::InternalPixelType;
+  using InputPixelType = typename TInputImage::PixelType;
+  using InputInternalPixelType = typename TInputImage::InternalPixelType;
+
+  /** Extract some information from the image types.  Dimensionality
+   * of the two images is assumed to be the same. */
+  static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension;
+
+  /** Image type alias support. */
+  using InputImageType = TInputImage;
+  using OutputImageType = TOutputImage;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(HigherOrderAccurateDerivativeImageFilter, ImageToImageFilter);
+
+  /** The output pixel type must be signed. */
+#ifdef ITK_USE_CONCEPT_CHECKING
+  /** Begin concept checking */
+  itkConceptMacro(SignedOutputPixelType, (Concept::Signed<OutputPixelType>));
+  /** End concept checking */
+#endif
+
+  /** Standard get/set macros for filter parameters. */
+  itkSetMacro(Order, unsigned int);
+  itkGetConstMacro(Order, unsigned int);
+  itkSetMacro(OrderOfAccuracy, unsigned int);
+  itkGetConstMacro(OrderOfAccuracy, unsigned int);
+  itkSetMacro(Direction, unsigned int);
+  itkGetConstMacro(Direction, unsigned int);
+
+  /** Use the image spacing information in calculations. Use this option if you
+   *  want derivatives in physical space. Default is UseImageSpacingOn. */
+  void
+  SetUseImageSpacingOn()
+  {
+    this->SetUseImageSpacing(true);
+  }
+
+  /** Ignore the image spacing. Use this option if you want derivatives in
+      isotropic pixel space.  Default is UseImageSpacingOn. */
+  void
+  SetUseImageSpacingOff()
+  {
+    this->SetUseImageSpacing(false);
+  }
+
+  /** Set/Get whether or not the filter will use the spacing of the input
+      image in its calculations */
+  itkSetMacro(UseImageSpacing, bool);
+  itkGetConstMacro(UseImageSpacing, bool);
+
+protected:
+  HigherOrderAccurateDerivativeImageFilter() = default;
+
+  ~HigherOrderAccurateDerivativeImageFilter() override = default;
+  void
+  PrintSelf(std::ostream & os, Indent indent) const override;
+
+  /** HigherOrderAccurateDerivativeImageFilter needs a larger input requested region than
+   * the output requested region (larger in the direction of the
+   * derivative).  As such, HigherOrderAccurateDerivativeImageFilter needs to provide an
+   * implementation for GenerateInputRequestedRegion() in order to
+   * inform the pipeline execution model.
+   *
+   * \sa ImageToImageFilter::GenerateInputRequestedRegion() */
+  void
+  GenerateInputRequestedRegion() override;
+
+  /** Standard pipeline method. While this class does not implement a
+   * ThreadedGenerateData(), its GenerateData() delegates all
+   * calculations to an NeighborhoodOperatorImageFilter.  Since the
+   * NeighborhoodOperatorImageFilter is multithreaded, this filter is
+   * multithreaded by default. */
+  void
+  GenerateData() override;
+
+private:
+  /** The order of the derivative. */
+  unsigned int m_Order{ 1 };
+
+  /** Order of accuracy. */
+  unsigned int m_OrderOfAccuracy{ 2 };
+
+  /** The direction of the derivative. */
+  unsigned int m_Direction{ 0 };
+
+  bool m_UseImageSpacing{ true };
+};
+
+} // end namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#  include "itkHigherOrderAccurateDerivativeImageFilter.hxx"
+#endif
+
+#endif

Modules/Filtering/HigherOrderAccurateGradient/include/itkHigherOrderAccurateDerivativeImageFilter.hxx

@@ -0,0 +1,161 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         https://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkHigherOrderAccurateDerivativeImageFilter_hxx
+#define itkHigherOrderAccurateDerivativeImageFilter_hxx
+
+#include "itkNumericTraits.h"
+#include "itkNeighborhoodOperatorImageFilter.h"
+#include "itkHigherOrderAccurateDerivativeOperator.h"
+#include "itkProgressAccumulator.h"
+
+namespace itk
+{
+
+template <typename TInputImage, typename TOutputImage>
+void
+HigherOrderAccurateDerivativeImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion()
+{
+  // call the superclass' implementation of this method. this should
+  // copy the output requested region to the input requested region
+  Superclass::GenerateInputRequestedRegion();
+
+  // get pointers to the input and output
+  typename Superclass::InputImagePointer inputPtr = const_cast<InputImageType *>(this->GetInput());
+
+  if (!inputPtr)
+  {
+    return;
+  }
+
+  // Build an operator so that we can determine the kernel size
+  HigherOrderAccurateDerivativeOperator<OutputPixelType, ImageDimension> oper;
+  oper.SetDirection(m_Direction);
+  oper.SetOrder(m_Order);
+  oper.SetOrderOfAccuracy(m_OrderOfAccuracy);
+  oper.CreateDirectional();
+
+  // get a copy of the input requested region (should equal the output
+  // requested region)
+  typename TInputImage::RegionType inputRequestedRegion;
+  inputRequestedRegion = inputPtr->GetRequestedRegion();
+
+  // pad the input requested region by the operator radius
+  inputRequestedRegion.PadByRadius(oper.GetRadius());
+
+  // crop the input requested region at the input's largest possible region
+  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
+  {
+    inputPtr->SetRequestedRegion(inputRequestedRegion);
+    return;
+  }
+  else
+  {
+    // Couldn't crop the region (requested region is outside the largest
+    // possible region).  Throw an exception.
+
+    // store what we tried to request (prior to trying to crop)
+    inputPtr->SetRequestedRegion(inputRequestedRegion);
+
+    // build an exception
+    InvalidRequestedRegionError e(__FILE__, __LINE__);
+    e.SetLocation(ITK_LOCATION);
+    e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
+    e.SetDataObject(inputPtr);
+    throw e;
+  }
+}
+
+
+template <typename TInputImage, typename TOutputImage>
+void
+HigherOrderAccurateDerivativeImageFilter<TInputImage, TOutputImage>::GenerateData()
+{
+  ZeroFluxNeumannBoundaryCondition<TInputImage> nbc;
+
+  // Define the operator value type so that we can filter integral
+  // images and have the proper operator defined.
+  using OperatorValueType = typename NumericTraits<OutputPixelType>::RealType;
+
+  // Filter
+  HigherOrderAccurateDerivativeOperator<OperatorValueType, ImageDimension> oper;
+  oper.SetDirection(m_Direction);
+  oper.SetOrder(m_Order);
+  oper.SetOrderOfAccuracy(m_OrderOfAccuracy);
+  oper.CreateDirectional();
+  oper.FlipAxes();
+
+  if (m_UseImageSpacing)
+  {
+    if (this->GetInput()->GetSpacing()[m_Direction] == 0.0)
+    {
+      itkExceptionMacro(<< "Image spacing cannot be zero.");
+    }
+    else
+    {
+      oper.ScaleCoefficients(1.0 / this->GetInput()->GetSpacing()[m_Direction]);
+    }
+  }
+
+  typename NeighborhoodOperatorImageFilter<InputImageType, OutputImageType, OperatorValueType>::Pointer filter =
+    NeighborhoodOperatorImageFilter<InputImageType, OutputImageType, OperatorValueType>::New();
+
+  // Create a process accumulator for tracking the progress of this minipipeline
+  ProgressAccumulator::Pointer progress = ProgressAccumulator::New();
+  progress->SetMiniPipelineFilter(this);
+
+  // Register the filter with the with progress accumulator using
+  // equal weight proportion
+  progress->RegisterInternalFilter(filter, 1.0f);
+
+  filter->OverrideBoundaryCondition(&nbc);
+
+  //
+  // Set up the mini-pipline
+  //
+  filter->SetOperator(oper);
+  filter->SetInput(this->GetInput());
+
+  // Graft this filter's output to the mini-pipeline.  this sets up
+  // the mini-pipeline to write to this filter's output and copies
+  // region ivars and meta-data
+  filter->GraftOutput(this->GetOutput());
+
+  // Execute the mini-pipeline.
+  filter->Update();
+
+  // Graft the output of the mini-pipeline back onto the filter's output,
+  // this copies back the region ivars and meta-data.
+  this->GraftOutput(filter->GetOutput());
+}
+
+
+template <typename TInputImage, typename TOutputImage>
+void
+HigherOrderAccurateDerivativeImageFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream & os, Indent indent) const
+{
+  Superclass::PrintSelf(os, indent);
+
+  os << indent << "Order: " << m_Order << std::endl;
+  os << indent << "OrderOfAccuracy: " << m_OrderOfAccuracy << std::endl;
+  os << indent << "Direction: " << m_Direction << std::endl;
+  os << indent << "UseImageSpacing: " << m_UseImageSpacing << std::endl;
+}
+
+} // end namespace itk
+
+#endif