CC @ntustison

Thanks @PranjalSahu for kicking this off — point-to-plane ICP is a real gap in ITK and the structural choices in this PR (per-point normals via PixelType, integration with the v4 metric framework, multithreaded accumulation) are the right shape. The PR has been dormant for ~3 years without a primary driver, so we're closing it on grounds of inactivity and consolidating findings here so a future contributor can pick this up cold.

Findings from a 2026-05-10 algorithmic review

The current implementation has three correctness issues that mean it needs a rewrite (not a patch) before merging. Documenting them so the next attempt doesn't repeat them.

localDerivative = closestPoint - point;
for (int i = 0; i < pixel.Size(); ++i)
{
  localDerivative[i] = localDerivative[i] * pixel[i];   // element-wise — wrong
}

MeasureType e = 0;
for (unsigned int d = 0; d < PointDimension; ++d) e += (point[d] - closestPoint[d]) * pixel[d];
for (unsigned int d = 0; d < PointDimension; ++d) localDerivative[d] = -e * pixel[d];

threadLocalTransformDerivative[par] += new_jacobian[par] * (pointDerivative[0] + pointDerivative[1]);

Reference math (canonical)

• Loss: E(R, t) = Σᵢ (nᵢ · (R·sᵢ + t − qᵢ))²
• Per-point Jacobian (small-angle linearization, parameters [α, β, γ, tₓ, t_y, t_z]): Jᵢ = [sᵢ × nᵢ, nᵢ] — rotation block first, translation block second
• Normal equations: H Δx = −b where H = Σᵢ JᵢᵀJᵢ (6×6), b = Σᵢ Jᵢᵀ rᵢ (6×1), rᵢ = nᵢ · (sᵢ − qᵢ)
• Update reconstruction: rebuild R from (α, β, γ) using full ZYX-Euler trig per iteration (matches PCL constructTransformationMatrix and Open3D TransformVector6dToMatrix4d)

Reference: K.-L. Low, "Linear Least-Squares Optimization for Point-to-Plane ICP Surface Registration", UNC TR04-004 (2004). Cross-checked against:

• PCL registration/include/pcl/registration/impl/transformation_estimation_point_to_plane_lls.hpp
• Open3D cpp/open3d/pipelines/registration/TransformationEstimationPointToPlane.cpp

VTK's vtkIterativeClosestPointTransform is point-to-point only — confirms the gap this PR was filling is real.

Hints for a future driver

A working draft for a fresh-start takeover lived briefly at ~/src/ITK-p2plane during this review (now discarded) — concrete starting points:

1. Don't touch the base class. Override only GetLocalNeighborhoodValue and GetLocalNeighborhoodValueAndDerivative in the new metric class. The chain-rule logic in PointSetToPointSetMetricWithIndexv4::CalculateValueAndDerivative already does the right thing once the per-point gradient is correct.

2. Per-point implementation:

   // GetLocalNeighborhoodValueAndDerivative:
   const auto closestPoint = this->m_MovingTransformedPointsLocator->FindClosestPoint(point);
   MeasureType e = 0;
   for (unsigned int d = 0; d < PointDimension; ++d) e += (point[d] - closestPoint[d]) * pixel[d];
   measure = e;  // signed plane distance — base class accumulates squared via per-point²
   for (unsigned int d = 0; d < PointDimension; ++d) localDerivative[d] = -e * pixel[d];

3. NumPy reference + ground-truth tests: a stand-alone NumPy implementation of the canonical algorithm (Low 2004 / PCL / Open3D) was developed during the review and verified against a sphere-surface fixture — converges in 13 iterations, recovers (R, t) to ~1e-5 RMS error. Worth committing alongside any future C++ implementation as itkPointToPlanePointSetMetricNumPyReference.py so future contributors can re-derive expected fixture values rather than copy-pasting.

4. Test fixtures:

   • Sphere-surface (well-conditioned, all 6 DOF observable) — use as primary recovery test
   • Flat plane — should fail gracefully (singular H matrix; in-plane translation and normal-axis rotation are unobservable)
   • Convergence-rate comparison vs EuclideanDistancePointSetToPointSetMetricv4 on identical input — point-to-plane should converge in roughly half the iterations for surface-like data

5. PixelType convention: PointSet<CovariantVector<float, N>, N> is the natural fit for storing per-point unit normals; the sphere fixture trivially generates these as the position vectors.

6. Adjacent ecosystem: the original author's other point-cloud-registration PRs (ENH: Add RANSAC as ITK remote module #3733 RANSAC remote, ENH: Add FPFH as remote module #3734 FPFH remote, ENH: Add DistanceThreshold parameter in EuclideanDistance Metricv4 #3728 EuclideanDistance DistanceThreshold) form a coherent feature → matching → fine-alignment toolkit. Landing point-to-plane ICP completes that pipeline; it would be worth coordinating documentation/examples with those modules.

Discourse motivation

Multiple Discourse threads have asked for surface-to-surface ICP capability that point-to-plane addresses (faster convergence than point-to-point on smooth surfaces):

• t/4044 (2021-04-17) — "(coarse) point set registration in python"
• t/3010 (2020-04-27) — "VTK IterativeClosestPointTransform to ITK AffineTransform"
• t/2214 (2019-09-03) — "Non-rigid registration of two surface meshes with non-corresponding vertices"
• t/3742 (2021-01-04) — "Point cloud / depth map / mesh registration SimpleITK/ITK Python"
• t/7410 (2025-01-28) — "Registration of MRI and Ultrasound volumes" (often uses surface ICP as a step)

The functional motivation is solid; this PR is closed only because it stalled without a champion. Reopening or starting fresh against a new ITK base is encouraged when someone has the bandwidth to drive it through review.

Closing on grounds of inactivity (3+ years, WIP, no primary driver). Thanks again @PranjalSahu — the structural pattern, the FPFH/RANSAC remote modules, and the DistanceThreshold enhancement remain in active use.

I have access to AI agents now and using the review comments I have created some changes. I will check them more for correctness and will create a new PR.