From d87f266dffdbcef00407d260eab20f4738a98825 Mon Sep 17 00:00:00 2001
From: stnava <stnava@gmail.com>
Date: Sat, 25 Nov 2023 11:03:18 -0500
Subject: [PATCH] PERF: point-based distances in neighborhood search

---
 Utilities/itkSurfaceImageCurvature.h   | 13 +++++++------
 Utilities/itkSurfaceImageCurvature.hxx | 13 +++++++------
 2 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/Utilities/itkSurfaceImageCurvature.h b/Utilities/itkSurfaceImageCurvature.h
index 9fe4cae06..97cdbb9a4 100644
--- a/Utilities/itkSurfaceImageCurvature.h
+++ b/Utilities/itkSurfaceImageCurvature.h
@@ -54,6 +54,7 @@ class SurfaceImageCurvature final : public SurfaceCurvatureBase<TSurface, 3>
   };
   typedef Image<PixelType, itkGetStaticConstMacro(ImageDimension)> ImageType;
   typedef typename ImageType::IndexType                            IndexType;
+  typedef typename ImageType::SpacingType                          SpacingType;
   typedef typename ImageType::SizeType                             SizeType;
   typedef ImageRegionIteratorWithIndex<ImageType>                  ImageIteratorType;
   /** Image dimension. */
@@ -89,7 +90,7 @@ class SurfaceImageCurvature final : public SurfaceCurvatureBase<TSurface, 3>
    * mean shift algorithm to find the best neighborhood.
    */
   void
-  FindNeighborhood(unsigned int numMeanShifts = 0) override;
+  FindNeighborhood(unsigned int numMeanShifts = 2) override;
 
   void
   FindEuclideanNeighborhood(PointType p);
@@ -173,11 +174,10 @@ class SurfaceImageCurvature final : public SurfaceCurvatureBase<TSurface, 3>
   CurvatureAtIndex(IndexType index)
   {
     PointType p;
-
-    for (unsigned int k = 0; k < ImageDimension; k++)
-    {
-      p[k] = (RealType)index[k];
-    }
+    this->m_FunctionImage->TransformIndexToPhysicalPoint(index, p);
+//    for (unsigned int k = 0; k < ImageDimension; k++) {
+//      p[k] = (RealType)index[k];
+//    }
     this->SetOrigin(p);
     this->EstimateFrameFromGradient(index);
     this->FindNeighborhood();
@@ -272,6 +272,7 @@ class SurfaceImageCurvature final : public SurfaceCurvatureBase<TSurface, 3>
   float                                    m_Area;
   RealType                                 m_MinSpacing;
   typename VectorInterpolatorType::Pointer m_Vinterp;
+  SpacingType                              m_Spacing;
 };
 } // namespace itk
 
diff --git a/Utilities/itkSurfaceImageCurvature.hxx b/Utilities/itkSurfaceImageCurvature.hxx
index 4a6104daf..cff4b1dd9 100644
--- a/Utilities/itkSurfaceImageCurvature.hxx
+++ b/Utilities/itkSurfaceImageCurvature.hxx
@@ -166,7 +166,7 @@ SurfaceImageCurvature<TSurface>::FindEuclideanNeighborhood(
   this->m_FunctionImage->TransformPhysicalPointToIndex(tempp, oindex);
   for (unsigned int i = 0; i < ImageDimension; i++)
   {
-    rad[i] = (long)(m_NeighborhoodRadius);
+    rad[i] = (long)(m_NeighborhoodRadius/this->m_Spacing[i]+0.5);
   }
   m_ti.SetLocation(oindex);
   unsigned int temp = 0;
@@ -190,7 +190,8 @@ SurfaceImageCurvature<TSurface>::FindEuclideanNeighborhood(
             isorigin = false;
           }
           p[k] = ipt[k];
-          RealType delt = static_cast<RealType>(oindex[k]) - static_cast<RealType>(index[k]);
+          RealType delt = rootpoint[k] - ipt[k];
+//          RealType delt = static_cast<RealType>(oindex[k]) - static_cast<RealType>(index[k]);
           dist += delt * delt;
         }
         dist = sqrt(dist);
@@ -419,6 +420,7 @@ void
 SurfaceImageCurvature<TSurface>::EstimateNormalsFromGradient()
 {
   typename ImageType::Pointer image = GetInput();
+  this->m_Spacing = image->GetSpacing();
 
   if (!image)
   {
@@ -748,9 +750,9 @@ SurfaceImageCurvature<TSurface>::ComputeSurfaceArea()
 
   RealType area = 0.0;
   this->m_TotalArea = 0.0;
+  unsigned long ct = 0;
 
   ti.GoToBegin();
-  unsigned int ct = 0;
   while (!ti.IsAtEnd())
   {
     index = ti.GetIndex();
@@ -847,6 +849,7 @@ SurfaceImageCurvature<TSurface>::IntegrateFunctionOverSurface(bool norm)
 
   this->m_ti.Initialize(rad, this->GetInput(), image->GetLargestPossibleRegion());
   this->m_ti2.Initialize(rad2, this->GetInput(), image->GetLargestPossibleRegion());
+  unsigned long ct = 0;
 
   IndexType index;
 
@@ -855,8 +858,6 @@ SurfaceImageCurvature<TSurface>::IntegrateFunctionOverSurface(bool norm)
   //  std::cout << " begin integrate ";
 
   ti.GoToBegin();
-  unsigned int ct = 0;
-  //  std::cout << " begin while " << std::endl;
   while (!ti.IsAtEnd())
   {
     index = ti.GetIndex();
@@ -1039,9 +1040,9 @@ SurfaceImageCurvature<TSurface>::ComputeFrameOverDomain(unsigned int which)
   unsigned int  ct = 1;
   unsigned long ct2 = 0;
   RealType      kpix = 0;
-
   double thresh = 0.0;
 
+
   ti.GoToBegin();
   while (!ti.IsAtEnd())
   {