Merge pull request #32 from csiro-coasts/dev

Dev

conf/configure

@@ -741,6 +741,7 @@ enable_swan
 enable_ecology_custom
 enable_do_timing
 enable_mkl
+enable_kdtree
 with_mkl
 with_gsl
 enable_da
@@ -1385,6 +1386,7 @@ Optional Features:
 --enable-ecology-custom=standard|anm|... - The ecological model custom version
   --enable-do-timing      enables timing profile
   --enable-mkl      enables building with MKL
+  --enable-kdtree         Enable KDTREE feature
   --enable-da         enables building with the DA library
   --disable-openmp        do not use OpenMP
   --enable-omp              enables building with OpenMP
@@ -5553,6 +5555,17 @@ if test "$ENABLE_MKL"; then
 
 fi
 
+# Check whether --enable-kdtree was given.
+if test "${enable_kdtree+set}" = set; then :
+  enableval=$enable_kdtree; if test "$enableval" = "yes"; then
+                  { $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: KDTREE help for delaunay walk is enabled, results may vary." >&5
+$as_echo "$as_me: WARNING: KDTREE help for delaunay walk is enabled, results may vary." >&2;}
+                 $as_echo "#define USE_KDTREE 1" >>confdefs.h
+
+               fi
+fi
+
+
 
 # Check whether --with-mkl was given.
 if test "${with_mkl+set}" = set; then :

conf/configure.in

@@ -334,6 +334,16 @@ if test "$ENABLE_MKL"; then
    AC_DEFINE(HAVE_MKL)
 fi
 
+dnl ###########################################################################
+dnl #### kdtree vs delaunay
+dnl ###########################################################################
+AC_ARG_ENABLE([kdtree],
+              [AS_HELP_STRING([--enable-kdtree], [Enable KDTREE feature])],
+              [if test "$enableval" = "yes"; then
+                  AC_MSG_WARN([KDTREE help for delaunay walk is enabled, results may vary.])
+                 AC_DEFINE([USE_KDTREE])
+               fi])
+
 dnl ###########################################################################
 dnl #### Intel Math Kernel Library (MKL)
 dnl ###########################################################################

conf/ems_version.h

@@ -1,4 +1,4 @@
-#define EMS_VERSION "v1.5.2"
+#define EMS_VERSION "v1.5.3"
 
                 #define EMS_IS_RELEASE 1
 

lib/include/delaunay.h

@@ -12,7 +12,7 @@
  *  reserved. See the license file for disclaimer and full
  *  use/redistribution conditions.
  *
- *  $Id: delaunay.h 6922 2021-10-08 04:40:51Z her127 $
+ *  $Id: delaunay.h 7533 2024-04-24 04:03:08Z tho861 $
  *
  */
 
@@ -21,6 +21,12 @@
 
 #include "grid_utils.h"
 #include "istack.h"
+#include "ems_conf.h" // added so we can switch the kdtree on/off
+
+#ifdef USE_KDTREE
+    #include "nanoflann_wrappers.h" // wrapper for c++ -> c 
+#endif
+
 
 #if !defined(_DELAUNAY_STRUCT)
 #define _DELAUNAY_STRUCT
@@ -59,6 +65,10 @@ typedef struct {
                                  * triangles i-th point belongs to */
     int** point_triangles;      /* point_triangles[i][j] is index of j-th
                                  * triangle i-th point belongs to */
+    #ifdef USE_KDTREE
+        point* centroids;               /* store the centroids for the tree*/
+        KDTreeContext *rt;              /* kdtree (called it rt since we originally implemented an rtree)*/
+    #endif
 
     int nedges;
     int* edges;                 /* n-th edge is formed by points[edges[n*2]]

lib/include/ems.h

@@ -12,7 +12,7 @@
  *  reserved. See the license file for disclaimer and full
  *  use/redistribution conditions.
  *
- *  $Id: ems.h 7477 2023-12-22 01:31:36Z riz008 $
+ *  $Id: ems.h 7577 2024-05-31 00:36:37Z riz008 $
  */
 
 #ifndef	_EMS_H
@@ -98,7 +98,7 @@ int strncasecmp(const char *s1, const char *s2, int n);
 
 /* Release verions and getters */
 #define EMSLIB_MAJOR_VERSION 1
-#define EMSLIB_MINOR_VERSION 3 
+#define EMSLIB_MINOR_VERSION 4
 #define EMSLIB_PATCH_VERSION 1
 
 int get_emslib_major_vers(void);

lib/include/ems_conf.h.in

@@ -129,3 +129,6 @@
 
 /* GSL */
 #undef HAVE_GSL
+
+/* Define to enable KDTREE feature. */
+#undef USE_KDTREE

lib/include/nanoflann_wrappers.h

@@ -0,0 +1,24 @@
+// nanoflann_wrappers.h
+
+#include <stdio.h> // Include for printf
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+// typedef struct {
+//   double x;
+//   double y;
+//   double z;
+//   double *v;
+// } point;
+
+typedef struct KDTreeContext KDTreeContext;
+void* create_kdtree(const point* points, size_t N, size_t max_leaf_size);
+int find_nearest_vertex(void* context, const float query_pt[3]); 
+void delete_kdtree(void* context);
+
+#ifdef __cplusplus
+}
+#endif

lib/interp/delaunay/delaunay.c

@@ -13,7 +13,7 @@
  *  reserved. See the license file for disclaimer and full
  *  use/redistribution conditions.
  *
- *  $Id: delaunay.c 6924 2021-10-08 04:41:47Z her127 $
+ *  $Id: delaunay.c 7533 2024-04-24 04:03:08Z tho861 $
  *
  */
 
@@ -31,6 +31,10 @@
 #include "emsmath.h"
 #include "emsalloc.h"
 
+#ifdef USE_KDTREE
+    #include "errfn.h"
+#endif
+
 //
 int verbose = 0;
 
@@ -132,6 +136,10 @@ delaunay* delaunay_create()
     d->nflags = 0;
     d->nflagsallocated = 0;
     d->flagids = NULL;
+    #if (defined USE_KDTREE) || (defined USE_RTREE)
+        d->centroids = NULL;
+        d->rt = NULL;
+    #endif
 
     return d;
 }
@@ -365,6 +373,32 @@ static void tio2delaunay_v(struct triangulateio* tio_out, struct triangulateio*
     }
 }
 
+
+void output_triangles(delaunay* d, const char* base_filename) {
+    char filename[1024];
+    snprintf(filename, sizeof(filename), "%s_%p.txt", base_filename, (void*)d);
+
+    FILE* file = fopen(filename, "w");
+    if (!file) {
+        perror("Failed to open triangles data file");
+        return;
+    }
+
+    for (int i = 0; i < d->ntriangles; ++i) {
+        triangle* t = &d->triangles[i];
+        // Output format: TriangleID, x1, y1, x2, y2, x3, y3
+        fprintf(file, "%d, %f, %f, %f, %f, %f, %f\n", i,
+                d->points[t->vids[0]].x, d->points[t->vids[0]].y,
+                d->points[t->vids[1]].x, d->points[t->vids[1]].y,
+                d->points[t->vids[2]].x, d->points[t->vids[2]].y);
+    }
+
+    fclose(file);
+    printf("Triangle data written to %s\n", filename);
+}
+
+
+
 /* Builds Delaunay triangulation of the given array of points.
  *
  * @param np Number of points
@@ -439,6 +473,21 @@ delaunay* delaunay_build(int np, point points[], int ns, int segments[], int nh,
     tio_destroy(&tio_in);
     tio_destroy(&tio_out);
 
+    #ifdef USE_KDTREE
+        // make centroids for kdtree:
+        d->centroids = malloc(d->ntriangles * sizeof(point));
+        for (int i = 0; i < d->ntriangles; ++i) {
+            triangle* t = &d->triangles[i];
+            double centroidX = (d->points[t->vids[0]].x + d->points[t->vids[1]].x + d->points[t->vids[2]].x) / 3.0;
+            double centroidY = (d->points[t->vids[0]].y + d->points[t->vids[1]].y + d->points[t->vids[2]].y) / 3.0;
+            d->centroids[i].x = centroidX;
+            d->centroids[i].y = centroidY;
+            // output_triangles(d,"triangles.txt");
+        }
+        warn("building kdtree\n");
+        d->rt = create_kdtree(d->centroids, d->ntriangles, 10); // last arg is max leaf size, perhaps this could be in the prm/tran
+    #endif
+
     return d;
 }
 
@@ -477,6 +526,12 @@ void delaunay_destroy(delaunay* d)
         istack_destroy(d->t_out);
     if (d->flagids != NULL)
         free(d->flagids);
+    #if (defined USE_KDTREE)
+        if (d->rt != NULL)
+            delete_kdtree(d->rt);
+        if (d->centroids != NULL)
+            free(d->centroids);
+    #endif
     free(d);
 }
 
@@ -502,6 +557,12 @@ int delaunay_xytoi(delaunay* d, point* p, int id)
     if (p->x < d->xmin || p->x > d->xmax || p->y < d->ymin || p->y > d->ymax)
         return -1;
 
+    #ifdef USE_KDTREE
+        id = find_nearest_vertex(d->rt, (float[]){p->x, p->y});
+        if (id < 0) warn("kdtree could not find nearest centroid to  %.12f %.12f\n", p->x, p->y);
+        // printf("delaunay_xytoi : kdtree result = %i ", id);
+    #endif
+
     if (id < 0 || id > d->ntriangles)
         id = 0;
     t = &d->triangles[id];
@@ -518,6 +579,7 @@ int delaunay_xytoi(delaunay* d, point* p, int id)
         }
     } while (i < 3);
 
+    // printf("delaunay result = %i \n", id);
     return id;
 }
 
@@ -534,6 +596,12 @@ int delaunay_xytoi_ng(delaunay* d, point* p, int id)
     if (p->x < d->xmin || p->x > d->xmax || p->y < d->ymin || p->y > d->ymax)
         return id;
 
+    #ifdef USE_KDTREE
+        id = find_nearest_vertex(d->rt, (float[]){p->x, p->y});
+        if (id < 0) warn("kdtree could not find nearest centroid to  %.12f %.12f\n", p->x, p->y);
+        // printf("delaunay_xytoi_ng : kdtree result = %i ", id);
+    #endif
+
     if (id < 0 || id > d->ntriangles)
         id = 0;
     t = &d->triangles[id];
@@ -554,6 +622,7 @@ int delaunay_xytoi_ng(delaunay* d, point* p, int id)
         }
     } while (i < 3);
 
+    // printf("delaunay_xytoi : kdtree result = %i ", id);
     return id;
 }
 
@@ -591,6 +660,13 @@ int delaunay_xytoi_lag(delaunay* d, point* p, int id)
 
     if(fabs(p->x - d->points[t->vids[0]].x) < 1e-5 &&
        fabs(p->y - d->points[t->vids[0]].y) < 1e-5) return(id);
+
+    #ifdef USE_KDTREE
+        // // The walk goes into an infinite loop if seeded with an id so not implementing here.
+        // id = find_nearest_vertex(d->rt, (float[]){p->x, p->y});
+        // printf("delaunay_xytoi_lag : kdtree result = %i ", id);
+    #endif
+
     do {
         for (i = 0; i < 3; ++i) {
 	  int i1 = (i + 1) % 3, i2;
@@ -626,7 +702,7 @@ int delaunay_xytoi_lag(delaunay* d, point* p, int id)
             }
         }
     } while (i < 3);
-
+    // printf("delaunay_xytoi_lag : del = %i  \n", id);
     return id;
 }
 
@@ -646,6 +722,13 @@ int delaunay_xytoi_lago(delaunay* d, point* p, int id)
 
     if(fabs(p->x - d->points[t->vids[0]].x) < 1e-5 &&
        fabs(p->y - d->points[t->vids[0]].y) < 1e-5) return(id);
+
+    #ifdef USE_KDTREE
+        // // This function delaunay_xytoi_lago doesn't seem to be called anyway so not including
+        // id = find_nearest_vertex(d->rt, (float[]){p->x, p->y});
+        // printf("delaunay_xytoi_lago : kdtree result = %i ", id);
+    #endif
+
     do {
         for (i = 0; i < 3; ++i) {
             int i1 = (i + 1) % 3;
@@ -681,6 +764,7 @@ int delaunay_xytoi_lago(delaunay* d, point* p, int id)
         }
     } while (i < 3);
 
+    // printf(" delaunay_xytoi_lago : delaunay result = %i  ", id);
     return id;
 }