limit maxNeighbors to knnWeights.maxNeighbors()

src/knn_weights.h

@@ -89,9 +89,9 @@ public:
   inline bool
   distance_is_less(double distance, double distx, size_t index, size_t index2)
   {
-    constexpr double cmp_tol = 1.e-12;
+    constexpr double cmpTolerance = 1.e-12;
     // return (distance < distx || (distance <= distx && index < index2));
-    return (distance + cmp_tol) < distx || (index < index2 && std::fabs(distance - distx) < cmp_tol);
+    return (distance + cmpTolerance) < distx || (index < index2 && std::fabs(distance - distx) < cmpTolerance);
   }
 
   inline void
@@ -126,6 +126,12 @@ public:
           }
       }
   }
+
+  inline void
+  store_distance(size_t index, double distance)
+  {
+    store_distance(index, distance, m_maxNeighbors);
+  }
   /*
   inline void
   set_distance(const size_t *indices, const double *distance, size_t numNeighbors)

src/pointsearch_healpix.h

@@ -53,24 +53,24 @@ private:
 
     auto sqrSearchRadius = cdo::sqr(searchRadius);
 
-    double xyz[3];
+    double searchPoint[3];
     size_t numWeights = 0;
     for (size_t i = 0; i < numIndices; ++i)
       {
-        gcLLtoXYZ(lons[i], lats[i], xyz);
+        gcLLtoXYZ(lons[i], lats[i], searchPoint);
         // Find distance to this point
-        double sqrDist = (float) cdo::sqr_distance(queryPoint, xyz);
+        double sqrDistance = (float) cdo::sqr_distance(queryPoint, searchPoint);
         // Store the index and distance if this is one of the smallest so far
-        if (sqrDist <= sqrSearchRadius)
+        if (sqrDistance <= sqrSearchRadius)
           {
             indices[numWeights] = indices[i];
-            distances[numWeights] = std::sqrt(sqrDist);
+            distances[numWeights] = std::sqrt(sqrDistance);
             numWeights++;
           }
       }
 
-    numWeights = std::min(numWeights, knnWeights.maxNeighbors());
-    for (size_t i = 0; i < numWeights; ++i) { knnWeights.store_distance(indices[i], distances[i], numWeights); }
+    auto maxNeighbors = std::min(numWeights, knnWeights.maxNeighbors());
+    for (size_t i = 0; i < numWeights; ++i) { knnWeights.store_distance(indices[i], distances[i], maxNeighbors); }
 
     knnWeights.check_distance();
   }

src/pointsearch_reg2d.cc

@@ -39,18 +39,16 @@ fill_src_indices(bool isCyclic, long nx, long ny, long ii, long jj, long k, size
   return numIndices;
 }
 
-/*
 void
-PointsearchReg2d::compute()
+PointsearchReg2d::compute_point(size_t index, double (&xyz)[3])
 {
-  auto iy = index / nx;
-  auto ix = index - iy * nx;
-
+  auto iy = index / m_nx;
+  auto ix = index - iy * m_nx;
   xyz[0] = m_cosLats[iy] * m_cosLons[ix];
   xyz[1] = m_cosLats[iy] * m_sinLons[ix];
   xyz[2] = m_sinLats[iy];
 }
-*/
+
 void
 PointsearchReg2d::store_distance_reg2d(double plon, double plat, knnWeightsType &knnWeights, size_t nx, size_t numIndices,
                                        size_t *indices, double *distances, double searchRadius)
@@ -59,30 +57,24 @@ PointsearchReg2d::store_distance_reg2d(double plon, double plat, knnWeightsType
   gcLLtoXYZ(plon, plat, queryPoint);
   auto sqrSearchRadius = cdo::sqr(searchRadius);
 
-  double xyz[3];
+  double searchPoint[3];
   size_t numWeights = 0;
   for (size_t i = 0; i < numIndices; ++i)
     {
-      auto index = indices[i];
-      auto iy = index / nx;
-      auto ix = index - iy * nx;
-
-      xyz[0] = m_cosLats[iy] * m_cosLons[ix];
-      xyz[1] = m_cosLats[iy] * m_sinLons[ix];
-      xyz[2] = m_sinLats[iy];
+      compute_point(indices[i], searchPoint);
       // Find distance to this point
-      double sqrDist = (float) cdo::sqr_distance(queryPoint, xyz);
+      double sqrDistance = (float) cdo::sqr_distance(queryPoint, searchPoint);
       // Store the index and distance if this is one of the smallest so far
-      if (sqrDist <= sqrSearchRadius)
+      if (sqrDistance <= sqrSearchRadius)  // knnWeights.store_distance(index, std::sqrt(sqrDist));
         {
-          indices[numWeights] = index;
-          distances[numWeights] = std::sqrt(sqrDist);
+          indices[numWeights] = indices[i];
+          distances[numWeights] = std::sqrt(sqrDistance);
           numWeights++;
         }
     }
 
-  numWeights = std::min(numWeights, knnWeights.maxNeighbors());
-  for (size_t i = 0; i < numWeights; ++i) { knnWeights.store_distance(indices[i], distances[i], numWeights); }
+  auto maxNeighbors = std::min(numWeights, knnWeights.maxNeighbors());
+  for (size_t i = 0; i < numWeights; ++i) { knnWeights.store_distance(indices[i], distances[i], maxNeighbors); }
 
   knnWeights.check_distance();
 }

src/pointsearch_reg2d.h

@@ -75,6 +75,8 @@ private:
   size_t m_ny{ 0 };
   Varray<double> m_cosLats, m_sinLats;  // cosine, sine of grid lats (for distance)
   Varray<double> m_cosLons, m_sinLons;  // cosine, sine of grid lons (for distance)
+
+  void compute_point(size_t index, double (&xyz)[3]);
 };
 
 #endif