gridsearch: added method GS_KDSPH.

src/grid_search.cc

@@ -69,6 +69,7 @@ float distance(const float *restrict a, const float *restrict b)
 void gridsearch_set_method(const char *methodstr)
 {
   if      ( strcmp(methodstr, "kdtree")  == 0 ) gridsearch_method_nn = GS_KDTREE;
+  else if ( strcmp(methodstr, "kdsph")   == 0 ) gridsearch_method_nn = GS_KDSPH;
   else if ( strcmp(methodstr, "nearpt3") == 0 ) gridsearch_method_nn = GS_NEARPT3;
   else if ( strcmp(methodstr, "full")    == 0 ) gridsearch_method_nn = GS_FULL;
   else
@@ -156,6 +157,39 @@ struct kdNode *gs_create_kdtree(size_t n, const double *restrict lons, const dou
 }
 
 
+struct kdNode *gs_create_kdsph(size_t n, const double *restrict lons, const double *restrict lats)
+{
+  struct kd_point *pointlist = (struct kd_point *) Malloc(n * sizeof(struct kd_point)); // kd_point contains 3d point
+  // see  example_cartesian.c
+  if ( cdoVerbose ) printf("kdtree lib spherical init: n=%zu  nthreads=%d\n", n, ompNumThreads);
+  kdata_t min[2], max[2];
+  min[0] = min[1] =  1e9;
+  max[0] = max[1] = -1e9;
+  kdata_t *restrict point;
+#if defined(HAVE_OPENMP4)
+#pragma omp simd
+#endif
+  for ( size_t i = 0; i < n; i++ ) 
+    {
+      point = pointlist[i].point;
+      point[0] = lons[i];
+      point[1] = lats[i];
+      point[2] = 0; // dummy
+      for ( size_t j = 0; j < 2; ++j )
+        {
+          min[j] = point[j] < min[j] ? point[j] : min[j];
+          max[j] = point[j] > max[j] ? point[j] : max[j];
+        }
+      pointlist[i].index = i;
+    }
+
+  struct kdNode *kdt = kd_sph_buildTree(pointlist, n, min, max, ompNumThreads);
+  if ( pointlist ) Free(pointlist);
+
+  return kdt;
+}
+
+
 void gs_destroy_nearpt3(struct gsNear *near)
 {
   if ( near )
@@ -274,6 +308,7 @@ struct gridsearch *gridsearch_create_nn(size_t n, const double *restrict lons, c
   if ( n == 0 ) return gs;
 
   if      ( gs->method_nn == GS_KDTREE  ) gs->kdt  = gs_create_kdtree(n, lons, lats);
+  else if ( gs->method_nn == GS_KDSPH   ) gs->kdt  = gs_create_kdsph(n, lons, lats);
   else if ( gs->method_nn == GS_NEARPT3 ) gs->near = gs_create_nearpt3(n, lons, lats);
   else if ( gs->method_nn == GS_FULL    ) gs->full = gs_create_full(n, lons, lats);
 
@@ -341,6 +376,27 @@ kdNode *gs_nearest_kdtree(kdNode *kdt, double lon, double lat, double *prange)
 }
 
 
+kdNode *gs_nearest_kdsph(kdNode *kdt, double lon, double lat, double *prange)
+{
+  if ( kdt == NULL ) return NULL;
+  
+  float range0 = gs_set_range(prange);
+  kdata_t range = KDATA_SCALE(range0);
+
+  kdata_t point[2];
+  point[0] = lon;
+  point[1] = lat;
+
+  kdNode *node = kd_nearest(kdt, point, &range, 3);
+
+  float frange = KDATA_INVSCALE(range);
+  if ( !(frange < range0) ) node = NULL;
+  if ( prange ) *prange = frange;
+
+  return node;
+}
+
+
 unsigned gs_nearest_nearpt3(struct gsNear *near, double lon, double lat, double *prange)
 {
   size_t index = GS_NOT_FOUND;
@@ -429,6 +485,11 @@ size_t gridsearch_nearest(struct gridsearch *gs, double lon, double lat, double
           kdNode *node = gs_nearest_kdtree(gs->kdt, lon, lat, prange);
           if ( node ) index = (int) node->index;
         }
+      else if ( gs->method_nn == GS_KDSPH )
+        {
+          kdNode *node = gs_nearest_kdsph(gs->kdt, lon, lat, prange);
+          if ( node ) index = (int) node->index;
+        }
       else if ( gs->method_nn == GS_NEARPT3 )
         {
           index = gs_nearest_nearpt3(gs->near, lon, lat, prange);

src/grid_search.h

@@ -9,7 +9,7 @@
 #define GS_NOT_FOUND  ULONG_MAX
 
 
-enum T_GRIDSEARCH_METHOD_NN  {GS_FULL=1, GS_KDTREE, GS_NEARPT3};
+enum T_GRIDSEARCH_METHOD_NN  {GS_FULL=1, GS_KDTREE, GS_KDSPH, GS_NEARPT3};
 
 struct gsFull {
   unsigned n;

src/kdtreelib/kdtree.h

@@ -164,10 +164,8 @@ struct kd_thread_data *kd_buildArg(struct kd_point *points,
 struct kdNode *kd_buildTree(struct kd_point *points, unsigned long nPoints,
                             kdata_t *min, kdata_t *max, int dim, int max_threads);
 void *kd_doBuildTree(void *threadarg);
-struct kdNode *kd_sph_buildTree(struct kd_point *points,
-                                unsigned long nPoints,
-                                kdata_t *min, kdata_t *max,
-                                int max_threads);
+struct kdNode *kd_sph_buildTree(struct kd_point *points, unsigned long nPoints,
+                                kdata_t *min, kdata_t *max, int max_threads);
 void *kd_sph_doBuildTree(void *threadarg);
 
 /* Functions for range searches 

src/kdtreelib/kdtree_cartesian.cc

@@ -24,10 +24,10 @@ kdata_t square(const kdata_t x)
 static constexpr
 kdata_t kd_dist_sq(const kdata_t *restrict a, const kdata_t *restrict b)
 {
-  return (square((a[0]-b[0]))+square((a[1]-b[1]))+square((a[2]-b[2])));
+  return square((a[0]-b[0]))+square((a[1]-b[1]))+square((a[2]-b[2]));
 }
 
-inline
+
 kdata_t kd_min(kdata_t x, kdata_t y)
 {
     return x < y ? x : y;