Added interpolation to new base values function

radixmath/interpolate.cc

@@ -144,7 +144,7 @@ std::vector<float> interpolateValues(
     }
   }
 
-  // Write out initial values
+  // Write out final values
   radix_line("Interpolation complete. Final values:");
   radix("  ");
   for (float f : interpolatedValues)
@@ -156,19 +156,119 @@ std::vector<float> interpolateValues(
   return interpolatedValues;
 }
 
+/**
+ * @brief interpolateToOtherBaseValues
+ * @param baseValues
+ * @param newBaseValues
+ * @param valuesToInterpolate
+ * @param circular
+ * @param missingValues
+ * @return
+ */
 std::vector<float> interpolateToOtherBaseValues(
     const std::vector<float> &baseValues,
     const std::vector<float> &newBaseValues,
     const std::vector<float> &valuesToInterpolate, const bool circular,
     const float missingValues)
 {
+  std::vector<float> finalInterpolatedValues(newBaseValues.size());
+
+  radix_line("Interpolating values to new base values");
+  radix("  Initial base values: ");
+  radix_block(for (float f : baseValues) { radix(f << " "); });
+  radix_line("");
+  radix("  Initial interp values: ");
+  radix_block(for (float f : valuesToInterpolate) { radix(f << " "); });
+  radix_line("");
+
   // First interpolate across the existing data to ensure there are no missing
+  radix_line("  Initial interpolation to remove missing values in input data");
   std::vector<float> initialInterpolatedValues = interpolateValues(
       baseValues, valuesToInterpolate, circular, missingValues);
+  // If this fails, we have to fail this one too
+  if (initialInterpolatedValues.size() == 0)
+  {
+    radix_line("Initial interpolation failed - returning empty vector");
+    return std::vector<float>();
+  }
+  radix_line("  Initial interpolation complete");
 
-  // TODO the rest
+  // Iterate over the new base vector to calculate values at these positions
+  radix_line("  Interpolating to new base values");
+  for (size_t i = 0; i < newBaseValues.size(); ++i)
+  {
+    // Fill in values with a new base lower than the lowest original
+    if (newBaseValues[i] < baseValues.front())
+    {
+      radix_line("    Base value below range: filling with lowest");
+      finalInterpolatedValues[i] = initialInterpolatedValues.front();
+    }
+    // Fill in values with a new base higher than the highest original
+    else if (newBaseValues[i] > baseValues.back())
+    {
+      radix_line("    Base value above range: filling with highest");
+      finalInterpolatedValues[i] = initialInterpolatedValues.back();
+    }
+    // Fill in central values
+    else
+    {
+      radix_line("    Filling in internal value");
+      // Get the base values immediately above and below this
+      size_t lowerBaseIndex = 0, higherBaseIndex = 0;
+      for (size_t j = 0; j < baseValues.size(); ++j)
+      {
+        if (baseValues[j] > newBaseValues[i])
+        {
+          higherBaseIndex = j;
+          break;
+        }
+        else
+        {
+          lowerBaseIndex = j;
+        }
+      }
+
+      // Interpolate between these values
+      float lastGoodValue = initialInterpolatedValues[lowerBaseIndex],
+            nextGoodValue = initialInterpolatedValues[higherBaseIndex];
+
+      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;
+        }
+      }
+
+      finalInterpolatedValues[i] =
+          lastGoodValue +
+          ((nextGoodValue - lastGoodValue) *
+           ((newBaseValues[i] - baseValues[lowerBaseIndex]) /
+            (baseValues[higherBaseIndex] - baseValues[lowerBaseIndex])));
+
+      if (circular)
+      {
+        finalInterpolatedValues[i] = fmod(finalInterpolatedValues[i], 360.0);
+      }
+    }
+  }
+
+  radix_line("Interpolation to new base values complete");
+  radix("  Final base values: ");
+  radix_block(for (float f : newBaseValues) { radix(f << " "); });
+  radix_line("");
+  radix("  Final interp values: ");
+  radix_block(for (float f : finalInterpolatedValues) { radix(f << " "); });
+  radix_line("");
 
-  return std::vector<float>();
+  return finalInterpolatedValues;
 }
 
 }  // namespace radix

radixmath/tests/tstInterpolate.cc

@@ -107,10 +107,10 @@ TEST(Interpolate, NewBaseValues)
   float missingValue = -9999.f;
   float tolerance    = 0.0001;
 
-  std::vector<float> baseValues1{12.f, 23.f, 100.f, 100.f};
-  std::vector<float> newBaseValues1{10.f, 20.f, 30.f, 40.f, 50.f, 60.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> baseValues1{1000.f, 2000.f, 3000.f};
+  std::vector<float> newBaseValues1{500.f, 1500.f, 2500.f, 3500.f};
+  std::vector<float> interpValues1{10.5, 20.5, 30.5};
+  std::vector<float> expectValues1{10.5, 15.5, 25.5, 30.5};
   std::vector<float> testValues1 =
       interpolateToOtherBaseValues(baseValues1, newBaseValues1, interpValues1);
   ASSERT_EQ(expectValues1.size(), testValues1.size());
@@ -118,4 +118,29 @@ TEST(Interpolate, NewBaseValues)
   {
     EXPECT_NEAR(expectValues1[i], testValues1[i], expectValues1[i] * tolerance);
   }
+
+  std::vector<float> interpValuesMissing{10.5, missingValue, 30.5};
+  std::vector<float> testValuesMissing = interpolateToOtherBaseValues(
+      baseValues1, newBaseValues1, interpValuesMissing);
+  ASSERT_EQ(expectValues1.size(), testValuesMissing.size());
+  for (int i = 0; i < testValuesMissing.size(); ++i)
+  {
+    EXPECT_NEAR(expectValues1[i], testValuesMissing[i],
+                expectValues1[i] * tolerance);
+  }
+
+  // Uncomment if negative gradient base values required
+  //  std::vector<float> baseValues2{3000.f, 2000.f, 1000.f};
+  //  std::vector<float> newBaseValues2{3500.f, 2500.f, 1500.f, 500.f};
+  //  std::vector<float> interpValues2{30.5, 20.5, 10.5};
+  //  std::vector<float> expectValues2{30.5, 25.5, 15.5, 10.5};
+  //  std::vector<float> testValues2 =
+  //      interpolateToOtherBaseValues(baseValues2, newBaseValues2,
+  //      interpValues2);
+  //  ASSERT_EQ(expectValues2.size(), testValues2.size());
+  //  for (int i = 0; i < testValues2.size(); ++i)
+  //  {
+  //    EXPECT_NEAR(expectValues2[i], testValues2[i], expectValues2[i] *
+  //    tolerance);
+  //  }
 }