Added nanoflann.hpp.

src/Makefile.am

@@ -113,6 +113,7 @@ libcdo_la_SOURCES =            \
                namelist.cc      \
                namelist.h      \
                namelist_parser.cc     \
+               nanoflann.hpp    \
                normal.cc        \
                nth_element.cc   \
                nth_element.h   \

src/Makefile.in

@@ -622,12 +622,12 @@ libcdo_la_SOURCES = argument.h argument.cc array.h array.cc cdo_int.h \
 	grid_search.cc grid_search.h listarray.cc listarray.h list.cc \
 	list.h listbuf.cc listbuf.h merge_sort2.cc merge_sort2.h \
 	modules.cc modules.h namelist.cc namelist.h namelist_parser.cc \
-	normal.cc nth_element.cc nth_element.h operator_help.h \
-	par_io.cc par_io.h parse_literal.cc percentiles_hist.cc \
-	percentiles_hist.h percentiles.cc percentiles.h pipe.cc pipe.h \
-	pmlist.cc pmlist.h sellist.cc sellist.h \
-	pragma_omp_atomic_update.h printinfo.h process.cc process.h \
-	pstream.cc pstream.h pstream_write.h pstream_int.h \
+	nanoflann.hpp normal.cc nth_element.cc nth_element.h \
+	operator_help.h par_io.cc par_io.h parse_literal.cc \
+	percentiles_hist.cc percentiles_hist.h percentiles.cc \
+	percentiles.h pipe.cc pipe.h pmlist.cc pmlist.h sellist.cc \
+	sellist.h pragma_omp_atomic_update.h printinfo.h process.cc \
+	process.h pstream.cc pstream.h pstream_write.h pstream_int.h \
 	pthread_debug.cc pthread_debug.h readline.cc realtime.cc \
 	remap.h remaplib.cc remapsort.cc remap_scrip_io.cc \
 	remap_search_reg2d.cc remap_search_latbins.cc \

src/grid_search.cc

@@ -68,10 +68,11 @@ 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;
+  if      ( strcmp(methodstr, "kdtree")    == 0 ) gridsearch_method_nn = GS_KDTREE;
+  else if ( strcmp(methodstr, "nanoflann") == 0 ) gridsearch_method_nn = GS_NANOFLANN;
+  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
     cdoAbort("gridsearch method %s not available!", methodstr);
 }
@@ -167,6 +168,41 @@ kdTree_t *gs_create_kdtree(size_t n, const double *restrict lons, const double *
   return kdt;
 }
 
+static
+nfTree_t *gs_create_nanoflann(size_t n, const double *restrict lons, const double *restrict lats, struct gridsearch *gs)
+{
+  PointCloud<float> *pointcloud = new PointCloud<float>();
+  if ( cdoVerbose ) printf("nanoflann init 3D: n=%zu  nthreads=%d\n", n, ompNumThreads);
+
+  float min[3], max[3];
+  min[0] = min[1] = min[2] =  1e9;
+  max[0] = max[1] = max[2] = -1e9;
+  // Generating  Point Cloud
+  float restrict point[3];
+  pointcloud->pts.resize(n);
+#if defined(HAVE_OPENMP4)
+#pragma omp simd
+#endif
+  for ( size_t i = 0; i < n; i++ ) 
+    {
+      LLtoXYZ_f(lons[i], lats[i], point);
+      pointcloud->pts[i].x = point[0];
+      pointcloud->pts[i].y = point[1];
+      pointcloud->pts[i].z = point[2];
+      for ( unsigned j = 0; j < 3; ++j )
+        {
+          min[j] = point[j] < min[j] ? point[j] : min[j];
+          max[j] = point[j] > max[j] ? point[j] : max[j];
+        }
+    }
+
+  // construct a kd-tree index:
+  nfTree_t *kdt = new nfTree_t(3 /*dim*/, *pointcloud, nanoflann::KDTreeSingleIndexAdaptorParams(50 /* max leaf */) );
+  kdt->buildIndex();
+
+  return kdt;
+}
+
 static
 kdTree_t *gs_create_kdsph(size_t n, const double *restrict lons, const double *restrict lats)
 {
@@ -460,10 +496,11 @@ struct gridsearch *gridsearch_create_nn(size_t n, const double *restrict lons, c
   gs->n = n;
   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);
+  if      ( gs->method_nn == GS_KDTREE    ) gs->kdt  = gs_create_kdtree(n, lons, lats);
+  else if ( gs->method_nn == GS_NANOFLANN ) gs->nft  = gs_create_nanoflann(n, lons, lats, gs);
+  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);
 
   gs->search_radius = cdo_default_search_radius();
 
@@ -534,6 +571,36 @@ size_t gs_nearest_kdtree(kdTree_t *kdt, double lon, double lat, double *prange)
   return index;
 }
 
+static
+size_t gs_nearest_nanoflann(nfTree_t *nft, double lon, double lat, double *prange)
+{
+  size_t index = GS_NOT_FOUND;
+  if ( nft == NULL ) return index;
+  
+  float range0 = gs_set_range(prange);
+  float range = range0;
+
+  float query_pt[3];
+  LLtoXYZ_f(lon, lat, query_pt);
+
+  const size_t num_results = 1;
+  size_t ret_index;
+  float out_dist_sqr;
+  nanoflann::KNNResultSet<float> resultSet(num_results);
+  resultSet.init(&ret_index, &out_dist_sqr);
+  nft->findNeighbors(resultSet, query_pt, nanoflann::SearchParams(10));
+
+  index = ret_index;
+  *prange = out_dist_sqr;
+  //float frange = range;
+  //if ( !(frange < range0) ) node = NULL;
+  //if ( prange ) *prange = frange;
+
+  //if ( node ) index = node->index;
+
+  return index;
+}
+
 static
 size_t gs_nearest_kdsph(kdTree_t *kdt, double lon, double lat, double *prange)
 {
@@ -643,10 +710,11 @@ size_t gridsearch_nearest(struct gridsearch *gs, double lon, double lat, double
   if ( gs )
     {
       // clang-format off
-      if      ( gs->method_nn == GS_KDTREE )  index = gs_nearest_kdtree(gs->kdt, lon, lat, prange);
-      else if ( gs->method_nn == GS_KDSPH )   index = gs_nearest_kdsph(gs->kdt, lon, lat, prange);
-      else if ( gs->method_nn == GS_NEARPT3 ) index = gs_nearest_nearpt3(gs->near, lon, lat, prange);
-      else if ( gs->method_nn == GS_FULL )    index = gs_nearest_full(gs->full, lon, lat, prange);
+      if      ( gs->method_nn == GS_KDTREE )    index = gs_nearest_kdtree(gs->kdt, lon, lat, prange);
+      else if ( gs->method_nn == GS_NANOFLANN ) index = gs_nearest_nanoflann(gs->nft, lon, lat, prange);
+      else if ( gs->method_nn == GS_KDSPH )     index = gs_nearest_kdsph(gs->kdt, lon, lat, prange);
+      else if ( gs->method_nn == GS_NEARPT3 )   index = gs_nearest_nearpt3(gs->near, lon, lat, prange);
+      else if ( gs->method_nn == GS_FULL )      index = gs_nearest_full(gs->full, lon, lat, prange);
       else cdoAbort("gridsearch_nearest::method_nn undefined!");
       // clang-format on
     }

src/grid_search.h

-#ifndef _GRID_SEARCH_H_
-#define _GRID_SEARCH_H_
+#ifndef  GRID_SEARCH_H_
+#define  GRID_SEARCH_H_
 
 #include <stdbool.h>
 #include "kdtreelib/kdtree.h"
 #include "nearpt3c.h"
+#include "nanoflann.hpp"
 
 #define GS_NOT_FOUND  SIZE_MAX
 
+enum T_GRIDSEARCH_METHOD_NN  {GS_FULL=1, GS_NANOFLANN, GS_KDTREE, GS_KDSPH, GS_NEARPT3};
 
-enum T_GRIDSEARCH_METHOD_NN  {GS_FULL=1, GS_KDTREE, GS_KDSPH, GS_NEARPT3};
+template <typename T>
+struct PointCloud
+{
+  struct Point
+  {
+    T  x,y,z;
+  };
+
+  std::vector<Point>  pts;
+
+  // Must return the number of data points
+  inline size_t kdtree_get_point_count() const { return pts.size(); }
+
+  // Returns the dim'th component of the idx'th point in the class:
+  // Since this is inlined and the "dim" argument is typically an immediate value, the
+  //  "if/else's" are actually solved at compile time.
+  inline T kdtree_get_pt(const size_t idx, int dim) const
+  {
+    if (dim == 0) return pts[idx].x;
+    else if (dim == 1) return pts[idx].y;
+    else return pts[idx].z;
+  }
+
+  // Optional bounding-box computation: return false to default to a standard bbox computation loop.
+  //   Return true if the BBOX was already computed by the class and returned in "bb" so it can be avoided to redo it again.
+  //   Look at bb.size() to find out the expected dimensionality (e.g. 2 or 3 for point clouds)
+  template <class BBOX>
+  bool kdtree_get_bbox(BBOX& /* bb */) const { return false; }
+};
+
+typedef nanoflann::KDTreeSingleIndexAdaptor<
+    nanoflann::L2_Simple_Adaptor<float, PointCloud<float> > ,
+    PointCloud<float>,
+    3 /* dim */
+    > nfTree_t;
 
 struct gsFull {
   size_t n;
@@ -35,6 +71,7 @@ struct gridsearch {
 
   struct gsNear *near;
   struct kdTree *kdt;
+  nfTree_t *nft;
   struct gsFull *full;
 
   double search_radius;

src/text.h

-#ifndef _TEXT_H
-#define _TEXT_H
+#ifndef  TEXT_H
+#define  TEXT_H
 
 #include <stdio.h>
 
-#define RESET		0
-#define BRIGHT 		1
-#define DIM		2
-#define UNDERLINE 	4
-#define BLINK		5
-#define REVERSE		7
-#define HIDDEN		8
-
-#define BLACK 		0
-#define RED		1
-#define GREEN		2
-#define YELLOW		3
-#define BLUE		4
-#define MAGENTA		5
-#define CYAN		6
-#define	WHITE		7
-
+enum text_mode {RESET=0, BRIGHT=1, DIM=2, UNDERLINE=4, BLINK=5, REVERSE=7, HIDDEN=8};
+enum text_color {BLACK=0, RED=1, GREEN=2, YELLOW=3, BLUE=4, MAGENTA=5, CYAN=6, WHITE=7};
 
 #define COLOR_STDOUT (stdout_is_tty && CDO_Color)
 #define COLOR_STDERR (stderr_is_tty && CDO_Color)
@@ -27,4 +12,4 @@
 void set_text_color(FILE *fp, int attr, int fg);
 void reset_text_color(FILE *fp);
 
-#endif  /* _TEXT_H */
+#endif  /* TEXT_H */