Skip to content

Commit 970bba7

Browse files
dzenanzdzenanz
authored
Merge pull request #6242 from hjmjohnson/ingest-HigherOrderAccurateGradient
ENH: Ingest ITKHigherOrderAccurateGradient into Modules/Filtering
2 parents 59176b1 + b2712b1 commit 970bba7

18 files changed

Lines changed: 1233 additions & 54 deletions

‎Modules/Filtering/HigherOrderAccurateGradient/CMakeLists.txt‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
project(HigherOrderAccurateGradient)
2+
3+
itk_module_impl()

‎Modules/Filtering/HigherOrderAccurateGradient/README.md‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
# HigherOrderAccurateGradient
2+
3+
In-tree ITK module providing higher-order accurate numerical
4+
derivative and gradient image filters. The flagship classes
5+
`itk::HigherOrderAccurateDerivativeImageFilter`,
6+
`itk::HigherOrderAccurateDerivativeOperator`, and
7+
`itk::HigherOrderAccurateGradientImageFilter` compute centered finite
8+
differences of arbitrary even order on N-dimensional scalar images.
9+
10+
## Origin
11+
12+
Ingested from the standalone remote module
13+
[**InsightSoftwareConsortium/ITKHigherOrderAccurateGradient**](https://github.com/InsightSoftwareConsortium/ITKHigherOrderAccurateGradient)
14+
on 2026-05-08, following the v4 ingestion guidelines defined in
15+
InsightSoftwareConsortium/ITK#6204. The upstream repository will be
16+
archived read-only after this PR merges; it remains reachable at the
17+
URL above for historical reference (notably the `doc/` directory,
18+
which was intentionally left in the upstream archive).
19+
20+
## References
21+
22+
- McCormick M.M., Liu X., Jomier J., Marion C., Ibanez L.
23+
*Higher Order Accurate Derivative and Gradient Calculation in ITK.*
24+
The Insight Journal. January-December. 2014.
25+
<https://hdl.handle.net/10380/3231>

‎Modules/Filtering/HigherOrderAccurateGradient/include/itkHigherOrderAccurateDerivativeImageFilter.h‎

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
/*=========================================================================
2+
*
3+
* Copyright NumFOCUS
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* https://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
#ifndef itkHigherOrderAccurateDerivativeImageFilter_h
19+
#define itkHigherOrderAccurateDerivativeImageFilter_h
20+
21+
#include "itkImageToImageFilter.h"
22+
23+
namespace itk
24+
{
25+
26+
/** \class HigherOrderAccurateDerivativeImageFilter
27+
* \brief Computes the higher order accurate directional derivative of an image.
28+
* The directional derivative at each pixel location is computed by convolution
29+
* with a higher order accurate derivative operator of user-specified order.
30+
*
31+
* SetOrder() specifies the order of the derivative.
32+
*
33+
* SetDirection() specifies the direction of the derivative with respect to the
34+
* coordinate axes of the image.
35+
*
36+
* To specify the order of accuracy, use SetOrderOfAccuracy(). The
37+
* approximation will be accurate to two times the OrderOfAccuracy in terms of
38+
* Taylor series terms.
39+
*
40+
* \sa Image
41+
* \sa Neighborhood
42+
* \sa NeighborhoodOperator
43+
* \sa NeighborhoodIterator
44+
*
45+
* \ingroup ImageFeatureExtraction
46+
* \ingroup HigherOrderAccurateGradient
47+
*/
48+
template <typename TInputImage, typename TOutputImage>
49+
class HigherOrderAccurateDerivativeImageFilter : public ImageToImageFilter<TInputImage, TOutputImage>
50+
{
51+
public:
52+
ITK_DISALLOW_COPY_AND_MOVE(HigherOrderAccurateDerivativeImageFilter);
53+
54+
/** Standard class type alias. */
55+
using Self = HigherOrderAccurateDerivativeImageFilter;
56+
using Superclass = ImageToImageFilter<TInputImage, TOutputImage>;
57+
using Pointer = SmartPointer<Self>;
58+
using ConstPointer = SmartPointer<const Self>;
59+
60+
/** Extract some information from the image types. Dimensionality
61+
* of the two images is assumed to be the same. */
62+
using OutputPixelType = typename TOutputImage::PixelType;
63+
using OutputInternalPixelType = typename TOutputImage::InternalPixelType;
64+
using InputPixelType = typename TInputImage::PixelType;
65+
using InputInternalPixelType = typename TInputImage::InternalPixelType;
66+
67+
/** Extract some information from the image types. Dimensionality
68+
* of the two images is assumed to be the same. */
69+
static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension;
70+
71+
/** Image type alias support. */
72+
using InputImageType = TInputImage;
73+
using OutputImageType = TOutputImage;
74+
75+
/** Method for creation through the object factory. */
76+
itkNewMacro(Self);
77+
78+
/** Run-time type information (and related methods). */
79+
itkTypeMacro(HigherOrderAccurateDerivativeImageFilter, ImageToImageFilter);
80+
81+
/** The output pixel type must be signed. */
82+
#ifdef ITK_USE_CONCEPT_CHECKING
83+
/** Begin concept checking */
84+
itkConceptMacro(SignedOutputPixelType, (Concept::Signed<OutputPixelType>));
85+
/** End concept checking */
86+
#endif
87+
88+
/** Standard get/set macros for filter parameters. */
89+
itkSetMacro(Order, unsigned int);
90+
itkGetConstMacro(Order, unsigned int);
91+
itkSetMacro(OrderOfAccuracy, unsigned int);
92+
itkGetConstMacro(OrderOfAccuracy, unsigned int);
93+
itkSetMacro(Direction, unsigned int);
94+
itkGetConstMacro(Direction, unsigned int);
95+
96+
/** Use the image spacing information in calculations. Use this option if you
97+
* want derivatives in physical space. Default is UseImageSpacingOn. */
98+
void
99+
SetUseImageSpacingOn()
100+
{
101+
this->SetUseImageSpacing(true);
102+
}
103+
104+
/** Ignore the image spacing. Use this option if you want derivatives in
105+
isotropic pixel space. Default is UseImageSpacingOn. */
106+
void
107+
SetUseImageSpacingOff()
108+
{
109+
this->SetUseImageSpacing(false);
110+
}
111+
112+
/** Set/Get whether or not the filter will use the spacing of the input
113+
image in its calculations */
114+
itkSetMacro(UseImageSpacing, bool);
115+
itkGetConstMacro(UseImageSpacing, bool);
116+
117+
protected:
118+
HigherOrderAccurateDerivativeImageFilter() = default;
119+
120+
~HigherOrderAccurateDerivativeImageFilter() override = default;
121+
void
122+
PrintSelf(std::ostream & os, Indent indent) const override;
123+
124+
/** HigherOrderAccurateDerivativeImageFilter needs a larger input requested region than
125+
* the output requested region (larger in the direction of the
126+
* derivative). As such, HigherOrderAccurateDerivativeImageFilter needs to provide an
127+
* implementation for GenerateInputRequestedRegion() in order to
128+
* inform the pipeline execution model.
129+
*
130+
* \sa ImageToImageFilter::GenerateInputRequestedRegion() */
131+
void
132+
GenerateInputRequestedRegion() override;
133+
134+
/** Standard pipeline method. While this class does not implement a
135+
* ThreadedGenerateData(), its GenerateData() delegates all
136+
* calculations to an NeighborhoodOperatorImageFilter. Since the
137+
* NeighborhoodOperatorImageFilter is multithreaded, this filter is
138+
* multithreaded by default. */
139+
void
140+
GenerateData() override;
141+
142+
private:
143+
/** The order of the derivative. */
144+
unsigned int m_Order{ 1 };
145+
146+
/** Order of accuracy. */
147+
unsigned int m_OrderOfAccuracy{ 2 };
148+
149+
/** The direction of the derivative. */
150+
unsigned int m_Direction{ 0 };
151+
152+
bool m_UseImageSpacing{ true };
153+
};
154+
155+
} // end namespace itk
156+
157+
#ifndef ITK_MANUAL_INSTANTIATION
158+
# include "itkHigherOrderAccurateDerivativeImageFilter.hxx"
159+
#endif
160+
161+
#endif

‎Modules/Filtering/HigherOrderAccurateGradient/include/itkHigherOrderAccurateDerivativeImageFilter.hxx‎

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
/*=========================================================================
2+
*
3+
* Copyright NumFOCUS
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* https://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
#ifndef itkHigherOrderAccurateDerivativeImageFilter_hxx
19+
#define itkHigherOrderAccurateDerivativeImageFilter_hxx
20+
21+
#include "itkNumericTraits.h"
22+
#include "itkNeighborhoodOperatorImageFilter.h"
23+
#include "itkHigherOrderAccurateDerivativeOperator.h"
24+
#include "itkProgressAccumulator.h"
25+
26+
namespace itk
27+
{
28+
29+
template <typename TInputImage, typename TOutputImage>
30+
void
31+
HigherOrderAccurateDerivativeImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion()
32+
{
33+
// call the superclass' implementation of this method. this should
34+
// copy the output requested region to the input requested region
35+
Superclass::GenerateInputRequestedRegion();
36+
37+
// get pointers to the input and output
38+
typename Superclass::InputImagePointer inputPtr = const_cast<InputImageType *>(this->GetInput());
39+
40+
if (!inputPtr)
41+
{
42+
return;
43+
}
44+
45+
// Build an operator so that we can determine the kernel size
46+
HigherOrderAccurateDerivativeOperator<OutputPixelType, ImageDimension> oper;
47+
oper.SetDirection(m_Direction);
48+
oper.SetOrder(m_Order);
49+
oper.SetOrderOfAccuracy(m_OrderOfAccuracy);
50+
oper.CreateDirectional();
51+
52+
// get a copy of the input requested region (should equal the output
53+
// requested region)
54+
typename TInputImage::RegionType inputRequestedRegion;
55+
inputRequestedRegion = inputPtr->GetRequestedRegion();
56+
57+
// pad the input requested region by the operator radius
58+
inputRequestedRegion.PadByRadius(oper.GetRadius());
59+
60+
// crop the input requested region at the input's largest possible region
61+
if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
62+
{
63+
inputPtr->SetRequestedRegion(inputRequestedRegion);
64+
return;
65+
}
66+
else
67+
{
68+
// Couldn't crop the region (requested region is outside the largest
69+
// possible region). Throw an exception.
70+
71+
// store what we tried to request (prior to trying to crop)
72+
inputPtr->SetRequestedRegion(inputRequestedRegion);
73+
74+
// build an exception
75+
InvalidRequestedRegionError e(__FILE__, __LINE__);
76+
e.SetLocation(ITK_LOCATION);
77+
e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
78+
e.SetDataObject(inputPtr);
79+
throw e;
80+
}
81+
}
82+
83+
84+
template <typename TInputImage, typename TOutputImage>
85+
void
86+
HigherOrderAccurateDerivativeImageFilter<TInputImage, TOutputImage>::GenerateData()
87+
{
88+
ZeroFluxNeumannBoundaryCondition<TInputImage> nbc;
89+
90+
// Define the operator value type so that we can filter integral
91+
// images and have the proper operator defined.
92+
using OperatorValueType = typename NumericTraits<OutputPixelType>::RealType;
93+
94+
// Filter
95+
HigherOrderAccurateDerivativeOperator<OperatorValueType, ImageDimension> oper;
96+
oper.SetDirection(m_Direction);
97+
oper.SetOrder(m_Order);
98+
oper.SetOrderOfAccuracy(m_OrderOfAccuracy);
99+
oper.CreateDirectional();
100+
oper.FlipAxes();
101+
102+
if (m_UseImageSpacing)
103+
{
104+
if (this->GetInput()->GetSpacing()[m_Direction] == 0.0)
105+
{
106+
itkExceptionMacro(<< "Image spacing cannot be zero.");
107+
}
108+
else
109+
{
110+
oper.ScaleCoefficients(1.0 / this->GetInput()->GetSpacing()[m_Direction]);
111+
}
112+
}
113+
114+
typename NeighborhoodOperatorImageFilter<InputImageType, OutputImageType, OperatorValueType>::Pointer filter =
115+
NeighborhoodOperatorImageFilter<InputImageType, OutputImageType, OperatorValueType>::New();
116+
117+
// Create a process accumulator for tracking the progress of this minipipeline
118+
ProgressAccumulator::Pointer progress = ProgressAccumulator::New();
119+
progress->SetMiniPipelineFilter(this);
120+
121+
// Register the filter with the with progress accumulator using
122+
// equal weight proportion
123+
progress->RegisterInternalFilter(filter, 1.0f);
124+
125+
filter->OverrideBoundaryCondition(&nbc);
126+
127+
//
128+
// Set up the mini-pipline
129+
//
130+
filter->SetOperator(oper);
131+
filter->SetInput(this->GetInput());
132+
133+
// Graft this filter's output to the mini-pipeline. this sets up
134+
// the mini-pipeline to write to this filter's output and copies
135+
// region ivars and meta-data
136+
filter->GraftOutput(this->GetOutput());
137+
138+
// Execute the mini-pipeline.
139+
filter->Update();
140+
141+
// Graft the output of the mini-pipeline back onto the filter's output,
142+
// this copies back the region ivars and meta-data.
143+
this->GraftOutput(filter->GetOutput());
144+
}
145+
146+
147+
template <typename TInputImage, typename TOutputImage>
148+
void
149+
HigherOrderAccurateDerivativeImageFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream & os, Indent indent) const
150+
{
151+
Superclass::PrintSelf(os, indent);
152+
153+
os << indent << "Order: " << m_Order << std::endl;
154+
os << indent << "OrderOfAccuracy: " << m_OrderOfAccuracy << std::endl;
155+
os << indent << "Direction: " << m_Direction << std::endl;
156+
os << indent << "UseImageSpacing: " << m_UseImageSpacing << std::endl;
157+
}
158+
159+
} // end namespace itk
160+
161+
#endif

0 commit comments

Comments
 (0)