diff --git a/Documentation/Doxygen/doxygen.bib b/Documentation/Doxygen/doxygen.bib
index 85591545681..c5fe7b121dd 100644
--- a/Documentation/Doxygen/doxygen.bib
+++ b/Documentation/Doxygen/doxygen.bib
@@ -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},
diff --git a/Modules/Core/Common/include/itkSobelOperator.h b/Modules/Core/Common/include/itkSobelOperator.h
index 1fd4f47ad31..161d9280cef 100644
--- a/Modules/Core/Common/include/itkSobelOperator.h
+++ b/Modules/Core/Common/include/itkSobelOperator.h
@@ -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
@@ -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.
  *
@@ -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.
@@ -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;
+  }
+
+  /** 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>));
 
@@ -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
 
diff --git a/Modules/Core/Common/include/itkSobelOperator.hxx b/Modules/Core/Common/include/itkSobelOperator.hxx
index a3acf0384d7..1acff682ece 100644
--- a/Modules/Core/Common/include/itkSobelOperator.hxx
+++ b/Modules/Core/Common/include/itkSobelOperator.hxx
@@ -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
 {
@@ -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());
   }
-  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
 
diff --git a/Modules/Core/Common/test/itkSobelOperatorGTest.cxx b/Modules/Core/Common/test/itkSobelOperatorGTest.cxx
index 05bf24bc6a0..d5f491e813a 100644
--- a/Modules/Core/Common/test/itkSobelOperatorGTest.cxx
+++ b/Modules/Core/Common/test/itkSobelOperatorGTest.cxx
@@ -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];
@@ -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 } });
 }