Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 0 additions & 6 deletions Documentation/Doxygen/doxygen.bib
Original file line number Diff line number Diff line change
Expand Up @@ -1127,12 +1127,6 @@ @article{sled1998
doi = {10.1109/42.668698},
url = {https://doi.org/10.1109/42.668698}
}
@techreport{sobel1995,
title = {An Isotropic 3x3x3 Volume Gradient Operator},
author = {Irwin Sobel},
year = 1995,
institution = {Hewlett-Packard Laboratories}
}
@inbook{soille2004,
title = {Geodesic Transformations},
author = {Soille, Pierre},
Expand Down
36 changes: 29 additions & 7 deletions Modules/Core/Common/include/itkSobelOperator.h
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,7 @@ namespace itk
* The current implementation of the Sobel operator is for 2 and 3 dimensions only.
* The ND version is planned for future releases.
*
* The extension to 3D was described in \cite sobel1995.
*
* The Sobel operator in 3D has the kernel
* The Sobel operator in 3D has the following kernel, when using legacy coefficients:
*
* \verbatim
* -1 -3 -1 0 0 0 1 3 1
Expand All @@ -73,6 +71,16 @@ namespace itk
* x-1 x x+1
* \endverbatim
*
* And it has the following kernel, when it is not using legacy coefficients:
*
* \verbatim
* -1 -2 -1 0 0 0 1 2 1
* -2 -4 -2 0 0 0 2 4 2
* -1 -2 -1 0 0 0 1 2 1
*
* x-1 x x+1
* \endverbatim
*
* The \c x kernel is just rotated as required to obtain the kernel in the
* \c y and \c z directions.
*
Expand Down Expand Up @@ -104,10 +112,6 @@ class ITK_TEMPLATE_EXPORT SobelOperator : public NeighborhoodOperator<TPixel, VD

itkOverrideGetNameOfClassMacro(SobelOperator);

static_assert(
VDimension == 2 || VDimension == 3,
"The ND version of the Sobel operator has not been implemented. Currently only 2D and 3D versions are available.");

/** Creates the operator with length only in the specified direction.
* The radius of the operator will be 0 except along the axis on which
* the operator will work.
Expand All @@ -124,6 +128,21 @@ class ITK_TEMPLATE_EXPORT SobelOperator : public NeighborhoodOperator<TPixel, VD
* \sa CreateDirectional \sa Fill */
// virtual void CreateToRadius(const unsigned long);

/** Allows specifying whether or not `GenerateCoefficients()` should return the legacy coordinate values, compatible
* with ITK <= 5.4. */
void
UseLegacyCoefficients(const bool useLegacyCoefficients)
{
m_UseLegacyCoefficients = useLegacyCoefficients;
}
Comment on lines +133 to +137
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that itkSetMacro(UseLegacyCoefficients, bool) does not compile, because itk::SobelOperator is not derived from itk::Object, and it does not support this->Modified().

The macro (itkSetMacro) would make it SetUseLegacyCoefficients(bool _arg) instead. But I think UseLegacyCoefficients(const bool useLegacyCoefficients) is more readable. OK?

What about the other member function, IsUsingLegacyCoefficients()? Is that clear enough?


/** Tells whether or not `GenerateCoefficients()` returns the legacy coordinate values, compatible with ITK <= 5.4. */
bool
IsUsingLegacyCoefficients() const
{
return m_UseLegacyCoefficients;
}

protected:
itkConceptMacro(SignedOutputPixelType, (Concept::Signed<typename NumericTraits<TPixel>::ValueType>));

Expand All @@ -137,6 +156,9 @@ class ITK_TEMPLATE_EXPORT SobelOperator : public NeighborhoodOperator<TPixel, VD
/** Arranges coefficients spatially in the memory buffer. */
void
Fill(const CoefficientVector & coeff) override;

private:
bool m_UseLegacyCoefficients{ true };
};
} // namespace itk

Expand Down
96 changes: 76 additions & 20 deletions Modules/Core/Common/include/itkSobelOperator.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,10 @@
#define itkSobelOperator_hxx

#include "itkIndexRange.h"
#include "itkIntTypes.h" // For size_t.
#include "itkMath.h"
#include "itkObject.h"
#include <array>

namespace itk
{
Expand Down Expand Up @@ -60,32 +63,85 @@ template <typename TPixel, unsigned int VDimension, typename TAllocator>
auto
SobelOperator<TPixel, VDimension, TAllocator>::GenerateCoefficients() -> CoefficientVector
{
const unsigned int direction = this->GetDirection();

if constexpr (VDimension == 2)
if (const unsigned int direction{ this->GetDirection() }; direction < VDimension)
{
switch (direction)
if constexpr (VDimension == 3)
{
case 0:
return { -1, 0, 1, -2, 0, 2, -1, 0, 1 };
case 1:
return { -1, -2, -1, 0, 0, 0, 1, 2, 1 };
if (m_UseLegacyCoefficients)
{
switch (direction)
{
case 0:
return { -1, 0, 1, -3, 0, 3, -1, 0, 1, -3, 0, 3, -6, 0, 6, -3, 0, 3, -1, 0, 1, -3, 0, 3, -1, 0, 1 };
case 1:
return { -1, -3, -1, 0, 0, 0, 1, 3, 1, -3, -6, -3, 0, 0, 0, 3, 6, 3, -1, -3, -1, 0, 0, 0, 1, 3, 1 };
case 2:
return { -1, -3, -1, -3, -6, -3, -1, -3, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 1, 3, 6, 3, 1, 3, 1 };
}
}
}

// Create N-dimensional Sobel kernels (one per axis).
// Each kernel is size 3^N, stored in row-major order with the last dimension varying fastest.
// Standard Sobel definition: derivative = [-1, 0, 1], smoothing = [1, 2, 1].
// Kernel for axis a is: K_a(x) = d[x_a] * Product_{i != a} s[x_i], with x_i being element of {-1,0,1}.
static constexpr auto sobelKernels = [] {
constexpr int derivative[3] = { -1, 0, 1 };
constexpr int smoothing[3] = { 1, 2, 1 };

constexpr size_t total = Math::UnsignedPower(3, VDimension);

// Allocate one kernel per axis
std::array<std::array<int, total>, VDimension> kernels{};

// Iterate over all N-D indices in {0,1,2}^N (mapping to {-1,0,1})
for (size_t linear = 0; linear < total; ++linear)
{
// Convert linear index -> base-3 digits (idx[0..N-1])
const auto idx = [linear] {
std::array<int, VDimension> result{};
size_t temp = linear;
for (size_t i = 0; i < VDimension; ++i)
{
result[VDimension - 1 - i] = static_cast<int>(temp % 3);
temp /= 3;
}
return result;
}();

// For each axis, compute Sobel value at this point
for (size_t axis = 0; axis < VDimension; ++axis)
{
int val = derivative[idx[axis]];
if (val == 0)
{ // derivative center => whole product is 0
kernels[axis][linear] = 0;
}
else
{
for (size_t i = 0; i < VDimension; ++i)
{
if (i != axis)
{
val *= smoothing[idx[i]];
}
}
kernels[axis][linear] = val;
}
}
}
return kernels;
}();

// Note: It appears that `sobelKernels` is in the reverse order!
const auto & kernel = sobelKernels[VDimension - 1 - direction];
return CoefficientVector(kernel.cbegin(), kernel.cend());
Comment on lines +137 to +138
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not simply kernel = sobelKernels[direction]? I still don't understand why the kernels returned by the AI generated function from #5702 (comment) are in the reverse order, compared by ITK 🤷

}
if constexpr (VDimension == 3)
else
{
switch (direction)
{
case 0:
return { -1, 0, 1, -3, 0, 3, -1, 0, 1, -3, 0, 3, -6, 0, 6, -3, 0, 3, -1, 0, 1, -3, 0, 3, -1, 0, 1 };
case 1:
return { -1, -3, -1, 0, 0, 0, 1, 3, 1, -3, -6, -3, 0, 0, 0, 3, 6, 3, -1, -3, -1, 0, 0, 0, 1, 3, 1 };
case 2:
return { -1, -3, -1, -3, -6, -3, -1, -3, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 1, 3, 6, 3, 1, 3, 1 };
}
itkExceptionMacro("The direction value (" << direction << ") should be less than the dimensionality (" << VDimension
<< ").");
}
itkExceptionMacro("The direction value (" << direction << ") should be less than the dimensionality (" << VDimension
<< ").");
}
} // namespace itk

Expand Down
41 changes: 37 additions & 4 deletions Modules/Core/Common/test/itkSobelOperatorGTest.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -32,17 +32,22 @@ using ExpectedKernelType = std::vector<int>;

template <unsigned int VDimension>
void
CheckKernelCoordinates(const std::array<ExpectedKernelType, VDimension> & expectedKernels)
CheckKernelCoordinates(const bool useLegacyCoordinates,
const std::array<ExpectedKernelType, VDimension> & expectedKernels)
{
using PixelType = float;
for (unsigned int direction{}; direction < VDimension; ++direction)
{
SCOPED_TRACE("`VDimension` = " + std::to_string(VDimension) + ", `direction` = " + std::to_string(direction));
SCOPED_TRACE("`VDimension` = " + std::to_string(VDimension) + ", `useLegacyCoordinates` = " +
(useLegacyCoordinates ? "true" : "false") + ", `direction` = " + std::to_string(direction));

itk::SobelOperator<PixelType, VDimension> sobelOperator;
sobelOperator.SetDirection(direction);
EXPECT_EQ(sobelOperator.GetDirection(), direction);

sobelOperator.UseLegacyCoefficients(useLegacyCoordinates);
EXPECT_EQ(sobelOperator.IsUsingLegacyCoefficients(), useLegacyCoordinates);

sobelOperator.CreateToRadius(itk::Size<VDimension>::Filled(1));

const ExpectedKernelType & expectedKernel = expectedKernels[direction];
Expand Down Expand Up @@ -71,10 +76,38 @@ TEST(SobelOperator, ExerciseBasicObjectMethods)
// Checks that the coordinates of the kernels have the expected values.
TEST(SobelOperator, CheckKernelCoordinates)
{
CheckKernelCoordinates<2>(
{ ExpectedKernelType{ -1, 0, 1, -2, 0, 2, -1, 0, 1 }, ExpectedKernelType{ -1, -2, -1, 0, 0, 0, 1, 2, 1 } });
for (const bool useLegacyCoordinates : { false, true })
{
CheckKernelCoordinates<1>(useLegacyCoordinates, { ExpectedKernelType{ -1, 0, 1 } });
CheckKernelCoordinates<2>(
useLegacyCoordinates,
{ ExpectedKernelType{ -1, 0, 1, -2, 0, 2, -1, 0, 1 }, ExpectedKernelType{ -1, -2, -1, 0, 0, 0, 1, 2, 1 } });
CheckKernelCoordinates<4>(
useLegacyCoordinates,
{ ExpectedKernelType{ -1, 0, 1, -2, 0, 2, -1, 0, 1, -2, 0, 2, -4, 0, 4, -2, 0, 2, -1, 0, 1, -2, 0, 2, -1, 0, 1,
-2, 0, 2, -4, 0, 4, -2, 0, 2, -4, 0, 4, -8, 0, 8, -4, 0, 4, -2, 0, 2, -4, 0, 4, -2, 0, 2,
-1, 0, 1, -2, 0, 2, -1, 0, 1, -2, 0, 2, -4, 0, 4, -2, 0, 2, -1, 0, 1, -2, 0, 2, -1, 0, 1 },
ExpectedKernelType{ -1, -2, -1, 0, 0, 0, 1, 2, 1, -2, -4, -2, 0, 0, 0, 2, 4, 2, -1, -2, -1, 0, 0, 0, 1, 2, 1,
-2, -4, -2, 0, 0, 0, 2, 4, 2, -4, -8, -4, 0, 0, 0, 4, 8, 4, -2, -4, -2, 0, 0, 0, 2, 4, 2,
-1, -2, -1, 0, 0, 0, 1, 2, 1, -2, -4, -2, 0, 0, 0, 2, 4, 2, -1, -2, -1, 0, 0, 0, 1, 2, 1 },
ExpectedKernelType{ -1, -2, -1, -2, -4, -2, -1, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 2, 4, 2, 1, 2, 1,
-2, -4, -2, -4, -8, -4, -2, -4, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 4, 2, 4, 8, 4, 2, 4, 2,
-1, -2, -1, -2, -4, -2, -1, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 2, 4, 2, 1, 2, 1 },
ExpectedKernelType{ -1, -2, -1, -2, -4, -2, -1, -2, -1, -2, -4, -2, -4, -8, -4, -2, -4, -2, -1, -2, -1,
-2, -4, -2, -1, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 2, 4, 2, 1, 2, 1,
2, 4, 2, 4, 8, 4, 2, 4, 2, 1, 2, 1, 2, 4, 2, 1, 2, 1 } });
}

// For 3D, the legacy and the non-legacy coordinates are different.
CheckKernelCoordinates<3>(
true,
{ ExpectedKernelType{ -1, 0, 1, -3, 0, 3, -1, 0, 1, -3, 0, 3, -6, 0, 6, -3, 0, 3, -1, 0, 1, -3, 0, 3, -1, 0, 1 },
ExpectedKernelType{ -1, -3, -1, 0, 0, 0, 1, 3, 1, -3, -6, -3, 0, 0, 0, 3, 6, 3, -1, -3, -1, 0, 0, 0, 1, 3, 1 },
ExpectedKernelType{ -1, -3, -1, -3, -6, -3, -1, -3, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 1, 3, 6, 3, 1, 3, 1 } });
CheckKernelCoordinates<3>(
false,
{ ExpectedKernelType{ -1, 0, 1, -2, 0, 2, -1, 0, 1, -2, 0, 2, -4, 0, 4, -2, 0, 2, -1, 0, 1, -2, 0, 2, -1, 0, 1 },
ExpectedKernelType{ -1, -2, -1, 0, 0, 0, 1, 2, 1, -2, -4, -2, 0, 0, 0, 2, 4, 2, -1, -2, -1, 0, 0, 0, 1, 2, 1 },
ExpectedKernelType{ -1, -2, -1, -2, -4, -2, -1, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 2, 4, 2, 1, 2, 1 } });
}
Loading