Adding basic interpolate method to radixmath

radixmath/CMakeLists.txt

@@ -7,19 +7,22 @@ IF("${CMAKE_SYSTEM_NAME}" STREQUAL "Windows")
 ENDIF()
 
 SET(SOURCE
+interpolate.cc
+lognormal.cc
 matrix.cc
 normal.cc
-vector3d.cc
-lognormal.cc
 point3d.cc
 ray.cc
 util.cc
+vector3d.cc
 )
 
 SET(HEADERS
 annealing.hh
 annealing.i.hh
 constants.hh
+interpolate.hh
+lognormal.hh
 matrix.hh
 matrix.i.hh
 normal.hh
@@ -27,7 +30,6 @@ point3d.hh
 ray.hh
 util.hh
 vector3d.hh
-lognormal.hh
 )
 
 #

radixmath/interpolate.cc

+#include "radixmath/interpolate.hh"
+
+#include "radixbug/bug.hh"
+
+#include <cmath>
+
+namespace radix
+{
+std::vector<float> interpolateValues(
+    const std::vector<float> &baseValues,
+    const std::vector<float> &valuesToInterpolate, float missingValue,
+    bool circular)
+{
+  std::vector<float> interpolatedValues = valuesToInterpolate;
+
+  // Ensure pressure and other vector are same length for interpolation
+  if (baseValues.size() != valuesToInterpolate.size())
+  {
+    radix_line("Error! Pressure vector is not same size ("
+               << baseValues.size() << ") as vector to interpolate ("
+               << valuesToInterpolate.size() << ") - can't use to interpolate");
+    return std::vector<float>();
+  }
+
+  // Ensure we have no missing values in the pressure dataset
+  for (float pressure : baseValues)
+  {
+    if (pressure == missingValue)
+    {
+      radix_line("Error! Base vector contains "
+                 << missingValue << " values - can't use to interpolate");
+      radix("  Base vector:");
+      for (float p : baseValues)
+      {
+        radix(" " << p << ";");
+      }
+      radix_line("");
+      return std::vector<float>();
+    }
+  }
+
+  // Write out initial values
+  radix_line("Interpolating " << valuesToInterpolate.size()
+                              << " values; initial:");
+  radix("  ");
+  for (float f : valuesToInterpolate)
+  {
+    radix(f << " ");
+  }
+  radix_line("");
+  radix_line("Base values:");
+  radix("  ");
+  for (float f : baseValues)
+  {
+    radix(f << " ");
+  }
+  radix_line("");
+
+  // Loop through vector first to find first/last non-missing indices
+  size_t firstIndex = valuesToInterpolate.size(), lastIndex = 0;
+  bool foundFirst = false;
+  for (size_t i = 0; i < valuesToInterpolate.size(); ++i)
+  {
+    if (valuesToInterpolate[i] != missingValue)
+    {
+      if (!foundFirst)
+      {
+        firstIndex = i;
+        foundFirst = true;
+      }
+      lastIndex = i;
+    }
+  }
+  // Fill in values before first and after last indices
+  radix_line("  Found first ("
+             << firstIndex << ") and last (" << lastIndex
+             << ") non-missing indices; filling in values before & after...");
+  for (size_t i = 0; i < firstIndex; ++i)
+  {
+    interpolatedValues[i] = valuesToInterpolate[firstIndex];
+  }
+  for (size_t i = lastIndex + 1; i < valuesToInterpolate.size(); ++i)
+  {
+    interpolatedValues[i] = valuesToInterpolate[lastIndex];
+  }
+
+  // Fill in missing values in central parts of vector
+  radix_line("  Filling in missing data in central part of vector...");
+  for (size_t i = firstIndex + 1; i < lastIndex; ++i)
+  {
+    // Search for a missing value
+    size_t lastGood = i - 1;
+    if (valuesToInterpolate[i] == missingValue)
+    {
+      // Get the next good value
+      while (valuesToInterpolate[i] == missingValue)
+      {
+        i++;
+      }
+      size_t nextGood = i;
+      // Interpolate between the two good values
+      for (size_t j = lastGood + 1; j < nextGood; ++j)
+      {
+        float lastGoodValue = valuesToInterpolate[lastGood],
+              nextGoodValue = valuesToInterpolate[nextGood];
+
+        if (circular && fabs(lastGoodValue - nextGoodValue) > 180.0)
+        {
+          radix_line(
+              "    Circular interpolation with distance > 180 "
+              "degrees: performing correction");
+          if (lastGoodValue < nextGoodValue)
+          {
+            lastGoodValue += 360.0;
+          }
+          else
+          {
+            nextGoodValue += 360.0;
+          }
+        }
+
+        interpolatedValues[j] =
+            lastGoodValue + ((nextGoodValue - lastGoodValue) *
+                             ((baseValues[j] - baseValues[lastGood]) /
+                              (baseValues[nextGood] - baseValues[lastGood])));
+
+        if (circular)
+        {
+          interpolatedValues[j] = fmod(interpolatedValues[j], 360.0);
+        }
+      }
+    }
+  }
+
+  return interpolatedValues;
+}
+
+}  // namespace radix

radixmath/interpolate.hh

+/*
+ * File:    interpolate.hh
+ */
+
+#ifndef RADIX_RADIXMATH_INTEPOLATE_HH_
+#define RADIX_RADIXMATH_INTEPOLATE_HH_
+
+#include "radixcore/visibility.hh"
+
+#include <vector>
+
+namespace radix
+{
+std::vector<float> RADIX_PUBLIC
+interpolateValues(const std::vector<float> &baseValues,
+                  const std::vector<float> &valuesToInterpolate,
+                  float missingValue = -9999.f, bool circular = false);
+
+}  // namespace radix
+
+#endif

radixmath/tests/CMakeLists.txt

 INCLUDE(GoogleTest)
 
+ADD_GOOGLE_TEST(tstInterpolate.cc NP 1)
+ADD_GOOGLE_TEST(tstLognormal.cc NP 1)
 ADD_GOOGLE_TEST(tstMatrix.cc NP 1)
 ADD_GOOGLE_TEST(tstUtil.cc NP 1)
-ADD_GOOGLE_TEST(tstLognormal.cc NP 1)

radixmath/tests/tstInterpolate.cc

+#include "gtest/gtest.h"
+
+#include "radixmath/interpolate.hh"
+
+#include <vector>
+
+using namespace radix;
+
+TEST(Interpolate, InternalInterpolation)
+{
+  float missingValue = -9999.f;
+  float tolerance    = 0.0001;
+
+  std::vector<float> baseValues1{10.f, 20.f, 30.f, 40.f, 50.f, 60.f, 70.f};
+  std::vector<float> interpValues1{1.f, missingValue, 3.f, missingValue,
+                                   6.f, missingValue, 1.f};
+  std::vector<float> expectValues1{1.f, 2.f, 3.f, 4.5, 6.f, 3.5, 1.f};
+  std::vector<float> testValues1 =
+      interpolateValues(baseValues1, interpValues1);
+  ASSERT_EQ(expectValues1.size(), testValues1.size());
+  for (int i = 0; i < testValues1.size(); ++i)
+  {
+    EXPECT_NEAR(expectValues1[i], testValues1[i], expectValues1[i] * tolerance);
+  }
+
+  std::vector<float> baseValues2{1.f, 2.f, 4.f, 10.f, 22.f, 32.f, 52.f, 62.f};
+  std::vector<float> interpValues2{2.f,          missingValue, 32.f,
+                                   missingValue, 44.f,         missingValue,
+                                   missingValue, 84.f};
+  std::vector<float> expectValues2{2.f,  12.f, 32.f, 36.f,
+                                   44.f, 54.f, 74.f, 84.f};
+  std::vector<float> testValues2 =
+      interpolateValues(baseValues2, interpValues2);
+  ASSERT_EQ(expectValues2.size(), testValues2.size());
+  for (int i = 0; i < testValues2.size(); ++i)
+  {
+    EXPECT_NEAR(expectValues2[i], testValues2[i], expectValues2[i] * tolerance);
+  }
+
+  missingValue = 3e12;
+  std::vector<float> baseValues3{100.f, 50.f, 25.f, 12.5};
+  std::vector<float> interpValues3{16.f, missingValue, missingValue, 2.f};
+  std::vector<float> expectValues3{16.f, 8.f, 4.f, 2.f};
+  std::vector<float> testValues3 =
+      interpolateValues(baseValues3, interpValues3, missingValue);
+  ASSERT_EQ(expectValues3.size(), testValues3.size());
+  for (int i = 0; i < testValues3.size(); ++i)
+  {
+    EXPECT_NEAR(expectValues3[i], testValues3[i], expectValues3[i] * tolerance);
+  }
+}
+
+TEST(Interpolate, CircularInterpolation)
+{
+  float missingValue = 713.3;
+  float tolerance    = 0.0001;
+
+  std::vector<float> baseValues1{1.f, 2.f, 3.f, 4.f, 5.f, 6.f};
+  std::vector<float> interpValues1{90.f,  missingValue, 270.f,
+                                   350.f, missingValue, 30.f};
+  std::vector<float> expectValues1{90.f, 180.f, 270.f, 350.f, 10.f, 30.f};
+  std::vector<float> testValues1 =
+      interpolateValues(baseValues1, interpValues1, missingValue, true);
+  ASSERT_EQ(expectValues1.size(), testValues1.size());
+  for (int i = 0; i < testValues1.size(); ++i)
+  {
+    EXPECT_NEAR(expectValues1[i], testValues1[i], expectValues1[i] * tolerance);
+  }
+}
+
+TEST(Interpolate, FillInBeginning)
+{
+  float missingValue = -9999.f;
+  float tolerance    = 0.0001;
+
+  std::vector<float> baseValues1{12.f, 23.f, 100.f, 100.f};
+  std::vector<float> interpValues1{missingValue, missingValue, 103.23, 200.2};
+  std::vector<float> expectValues1{103.23, 103.23, 103.23, 200.2};
+  std::vector<float> testValues1 =
+      interpolateValues(baseValues1, interpValues1, missingValue, true);
+  ASSERT_EQ(expectValues1.size(), testValues1.size());
+  for (int i = 0; i < testValues1.size(); ++i)
+  {
+    EXPECT_NEAR(expectValues1[i], testValues1[i], expectValues1[i] * tolerance);
+  }
+}
+
+TEST(Interpolate, FillInEnd)
+{
+  float missingValue = -9999.f;
+  float tolerance    = 0.0001;
+
+  std::vector<float> baseValues1{12.f, 23.f, 100.f, 100.f};
+  std::vector<float> interpValues1{103.23, 200.2, missingValue, missingValue};
+  std::vector<float> expectValues1{103.23, 200.2, 200.2, 200.2};
+  std::vector<float> testValues1 =
+      interpolateValues(baseValues1, interpValues1, missingValue, true);
+  ASSERT_EQ(expectValues1.size(), testValues1.size());
+  for (int i = 0; i < testValues1.size(); ++i)
+  {
+    EXPECT_NEAR(expectValues1[i], testValues1[i], expectValues1[i] * tolerance);
+  }
+}