Isosurface: Added memory support for 32-bit float data.

ChangeLog

@@ -2,6 +2,10 @@
 
 	* Version 1.9.10 release
 
+2021-01-05  Uwe Schulzweida
+
+	* Isosurface: Added memory support for 32-bit float data.
+
 2020-12-17  Uwe Schulzweida
 
 	* Added warning message if a non-thread-safe NetCDF4/HDF5 library is used

src/Isosurface.cc

@@ -23,24 +23,21 @@
 #include "param_conversion.h"
 #include "interpol.h"
 
+template <typename T>
 static void
-isosurface(double isoval, int nlevels, const Varray<double> &levels, const FieldVector &field3D, Field &field2D)
+isosurface_kernel(double isoval, size_t nmiss, const Varray<double> &levels, int nlevels, size_t gridsize, T missval, const Varray<const T *> &data3D, Varray<T> &data2D)
 {
-  const auto gridsize = gridInqSize(field3D[0].grid);
-  const auto missval = field3D[0].missval;
-  auto &data2D = field2D.vec_d;
-
-  auto nmiss = field3D[0].nmiss;
-  for (int k = 1; k < nlevels; ++k) nmiss += field3D[k].nmiss;
-
+#ifdef _OPENMP
+#pragma omp parallel for default(shared)
+#endif
   for (size_t i = 0; i < gridsize; ++i)
     {
       data2D[i] = missval;
 
       for (int k = 0; k < (nlevels - 1); ++k)
         {
-          const auto val1 = field3D[k].vec_d[i];
-          const auto val2 = field3D[k + 1].vec_d[i];
+          const double val1 = data3D[k][i];
+          const double val2 = data3D[k + 1][i];
 
           if (nmiss)
             {
@@ -59,27 +56,48 @@ isosurface(double isoval, int nlevels, const Varray<double> &levels, const Field
             }
         }
     }
-
-  fieldNumMV(field2D);
 }
 
+
 static void
-layerValueMin(int nlevels, const FieldVector &field3D, Field &field2D)
+isosurface(double isoval, int nlevels, const Varray<double> &levels, const FieldVector &field3D, Field &field2D)
 {
   const auto gridsize = gridInqSize(field3D[0].grid);
   const auto missval = field3D[0].missval;
-  auto &data2D = field2D.vec_d;
 
   auto nmiss = field3D[0].nmiss;
   for (int k = 1; k < nlevels; ++k) nmiss += field3D[k].nmiss;
 
+  if (field3D[0].memType == MemType::Float)
+    {
+      std::vector<const float *> data3D(nlevels);
+      for (int k = 0; k < nlevels; ++k) data3D[k] = field3D[k].vec_f.data();
+      isosurface_kernel(isoval, nmiss, levels, nlevels, gridsize, (float)missval, data3D, field2D.vec_f);
+    }
+  else
+    {
+      std::vector<const double *> data3D(nlevels);
+      for (int k = 0; k < nlevels; ++k) data3D[k] = field3D[k].vec_d.data();
+      isosurface_kernel(isoval, nmiss, levels, nlevels, gridsize, missval, data3D, field2D.vec_d);
+    }
+
+  fieldNumMV(field2D);
+}
+
+template <typename T>
+static void
+layer_value_min_kernel(int nlevels, size_t gridsize, T missval, const Varray<const T *> &data3D, Varray<T> &data2D)
+{
+#ifdef _OPENMP
+#pragma omp parallel for default(shared)
+#endif
   for (size_t i = 0; i < gridsize; ++i)
     {
       data2D[i] = missval;
 
       for (int k = 0; k < nlevels; ++k)
         {
-          const auto val = field3D[k].vec_d[i];
+          const auto val = data3D[k][i];
           if (!DBL_IS_EQUAL(val, missval))
             {
               data2D[i] = val;
@@ -87,27 +105,47 @@ layerValueMin(int nlevels, const FieldVector &field3D, Field &field2D)
             }
         }
     }
-
-  fieldNumMV(field2D);
 }
 
 static void
-layerValueMax(int nlevels, const FieldVector &field3D, Field &field2D)
+layer_value_min(int nlevels, const FieldVector &field3D, Field &field2D)
 {
   const auto gridsize = gridInqSize(field3D[0].grid);
   const auto missval = field3D[0].missval;
-  auto &data2D = field2D.vec_d;
 
   auto nmiss = field3D[0].nmiss;
   for (int k = 1; k < nlevels; ++k) nmiss += field3D[k].nmiss;
 
+  if (field3D[0].memType == MemType::Float)
+    {
+      std::vector<const float *> data3D(nlevels);
+      for (int k = 0; k < nlevels; ++k) data3D[k] = field3D[k].vec_f.data();
+      layer_value_min_kernel(nlevels, gridsize, (float)missval, data3D, field2D.vec_f);
+    }
+  else
+    {
+      std::vector<const double *> data3D(nlevels);
+      for (int k = 0; k < nlevels; ++k) data3D[k] = field3D[k].vec_d.data();
+      layer_value_min_kernel(nlevels, gridsize, missval, data3D, field2D.vec_d);
+    }
+
+  fieldNumMV(field2D);
+}
+
+template <typename T>
+static void
+layer_value_max_kernel(int nlevels, size_t gridsize, T missval, const Varray<const T *> &data3D, Varray<T> &data2D)
+{
+#ifdef _OPENMP
+#pragma omp parallel for default(shared)
+#endif
   for (size_t i = 0; i < gridsize; ++i)
     {
       data2D[i] = missval;
 
       for (int k = nlevels - 1; k >= 0; --k)
         {
-          const auto val = field3D[k].vec_d[i];
+          const auto val = data3D[k][i];
           if (!DBL_IS_EQUAL(val, missval))
             {
               data2D[i] = val;
@@ -115,6 +153,29 @@ layerValueMax(int nlevels, const FieldVector &field3D, Field &field2D)
             }
         }
     }
+}
+
+static void
+layer_value_max(int nlevels, const FieldVector &field3D, Field &field2D)
+{
+  const auto gridsize = gridInqSize(field3D[0].grid);
+  const auto missval = field3D[0].missval;
+
+  auto nmiss = field3D[0].nmiss;
+  for (int k = 1; k < nlevels; ++k) nmiss += field3D[k].nmiss;
+
+  if (field3D[0].memType == MemType::Float)
+    {
+      std::vector<const float *> data3D(nlevels);
+      for (int k = 0; k < nlevels; ++k) data3D[k] = field3D[k].vec_f.data();
+      layer_value_max_kernel(nlevels, gridsize, (float)missval, data3D, field2D.vec_f);
+    }
+  else
+    {
+      std::vector<const double *> data3D(nlevels);
+      for (int k = 0; k < nlevels; ++k) data3D[k] = field3D[k].vec_d.data();
+      layer_value_max_kernel(nlevels, gridsize, missval, data3D, field2D.vec_d);
+    }
 
   fieldNumMV(field2D);
 }
@@ -178,9 +239,9 @@ Isosurface(void *process)
 
   const auto lpositive = !(levels[0] < 0.0 && levels[nlevels - 1] < 0.0);
   const auto lreverse = (levels[0] > levels[nlevels - 1]);
-  auto bottomValueFunc = lreverse ? layerValueMax : layerValueMin;
-  auto topValueFunc = lreverse ? layerValueMin : layerValueMax;
-  if (lpositive && positiveIsDown(zaxisID1)) std::swap(bottomValueFunc, topValueFunc);
+  auto bottom_value_func = lreverse ? layer_value_max : layer_value_min;
+  auto top_value_func = lreverse ? layer_value_min : layer_value_max;
+  if (lpositive && positiveIsDown(zaxisID1)) std::swap(bottom_value_func, top_value_func);
 
   const auto zaxisIDsfc = zaxisFromName("surface");
   for (int i = 0; i < nzaxis; i++)
@@ -191,7 +252,6 @@ Isosurface(void *process)
 
   VarList varList1;
   varListInit(varList1, vlistID1);
-  for (auto &var : varList1) var.memType = MemType::Double;
 
   const auto nvars = vlistNvars(vlistID1);
   std::vector<bool> lvar3D(nvars), vars(nvars);
@@ -232,8 +292,8 @@ Isosurface(void *process)
                   field2.init(varList1[varID]);
                   // clang-format off
                   if      (operatorID == ISOSURFACE)  isosurface(isoval, nlevels, levels, vars1[varID], field2);
-                  else if (operatorID == BOTTOMVALUE) bottomValueFunc(nlevels, vars1[varID], field2);
-                  else if (operatorID == TOPVALUE)    topValueFunc(nlevels, vars1[varID], field2);
+                  else if (operatorID == BOTTOMVALUE) bottom_value_func(nlevels, vars1[varID], field2);
+                  else if (operatorID == TOPVALUE)    top_value_func(nlevels, vars1[varID], field2);
                   // clang-format on
 
                   cdoDefRecord(streamID2, varID, 0);