Implemented nearest neighbour interpolation for grids

CMakeLists.txt

@@ -88,6 +88,7 @@ easi_add_test (5_function)
 easi_add_test (26_function)
 easi_add_test (33_layered_constant)
 easi_add_test (101_asagi)
+easi_add_test (101_asagi_nearest)
 easi_add_test (f_16_scec)
 easi_add_test (f_120_sumatra)
 easi_add_test (supplied_parameters)

doc/maps.rst

@@ -138,6 +138,7 @@ Looks up values using ASAGI (with trilinear interpolation).
     file: <string>
     parameters: [<parameter>,<parameter>,...]
     var: <string>
+    interpolation: (nearest|linear)
 
 :Domain:
   *inherited*
@@ -152,6 +153,8 @@ Looks up values using ASAGI (with trilinear interpolation).
   the NetCDF file
 :var:
   The NetCDF variable which holds the data (default: data)
+:interpolation:
+  Choose between nearest neighbour and linear interpolation (default: linear)
 
 SCECFile
 --------
@@ -163,6 +166,7 @@ http://scecdata.usc.edu/cvws/download/tpv16/TPV16\_17\_Description\_v03.pdf).
 
     !SCECFile
     file: <string>
+    interpolation: (nearest|linear)
 
 :Domain:
   *inherited*, must be 2D
@@ -171,7 +175,9 @@ http://scecdata.usc.edu/cvws/download/tpv16/TPV16\_17\_Description\_v03.pdf).
 :Example:
   `example <https://github.com/SeisSol/easi/blob/master/examples/f_16_scec.yaml#L8>`_
 :file:
-  Path to a SCEC stress file 
+  Path to a SCEC stress file
+:interpolation:
+  Choose between nearest neighbour and linear interpolation (default: linear)
 
 EvalModel
 ---------

examples/101_asagi_nearest.yaml

@@ -0,0 +1,12 @@
+!AffineMap
+matrix:
+  x: [1.0, 0.0, 0.0]
+  z: [0.0, 1.0, 0.0]
+translation:
+  x: 0.0
+  z: 0.0
+components: !ASAGI
+  file: ../examples/asagi_example.nc
+  parameters: [rho, mu, lambda]
+  var: data
+  interpolation: nearest

include/easi/component/ASAGI.h

@@ -64,6 +64,7 @@ class ASAGI : public Grid<ASAGI> {
 
   void setGrid(std::set<std::string> const& in, std::vector<std::string> const& parameters, asagi::Grid* grid, unsigned numberOfThreads);
 
+  void getNearestNeighbour(Slice<double> const& x, double* buffer);
   void getNeighbours(Slice<double> const& x, double* weights, double* buffer);
   unsigned permutation(unsigned index) const { return m_permutation[index]; }
 
@@ -135,6 +136,18 @@ void ASAGI::setGrid(std::set<std::string> const& in, std::vector<std::string> co
   }
 }
 
+void ASAGI::getNearestNeighbour(Slice<double> const& x, double* buffer) {
+  double pos[MaxDimensions];
+  float* bufferSP = reinterpret_cast<float*>(buffer);
+  for (unsigned d = 0; d < m_grid->getDimensions(); ++d) {
+    pos[d] = x(d);
+  }
+  m_grid->getBuf(bufferSP, pos);
+  for (int j = m_numValues-1; j >= 0; --j) {
+    buffer[j] = static_cast<double>(bufferSP[j]);
+  }
+}
+
 void ASAGI::getNeighbours(Slice<double> const& x, double* weights, double* buffer) {
   double lowPos[MaxDimensions];
   for (unsigned d = 0; d < m_grid->getDimensions(); ++d) {

include/easi/component/Grid.h

@@ -45,10 +45,22 @@ namespace easi {
 template<typename Derived>
 class Grid : public Map {
 public:
+  enum InterpolationType {
+    Nearest,
+    Linear
+  };
+
   virtual ~Grid() {}
 
   virtual Matrix<double> map(Matrix<double>& x);
 
+  void setInterpolationType(std::string const& interpolationType);
+  void setInterpolationType(enum InterpolationType interpolationType) { m_interpolationType = interpolationType; }
+
+  void getNearestNeighbour(Slice<double> const& x, double* buffer) {
+    static_cast<Derived*>(this)->getNearestNeighbour(x, buffer);
+  }
+
   void getNeighbours(Slice<double> const& x, double* weights, double* buffer) {
     static_cast<Derived*>(this)->getNeighbours(x, weights, buffer);
   }
@@ -58,6 +70,9 @@ class Grid : public Map {
 
 protected:
   virtual unsigned numberOfThreads() const = 0;
+
+private:
+  enum InterpolationType m_interpolationType = Linear;
 };
 
 template<typename GridImpl>
@@ -73,15 +88,19 @@ Matrix<double> Grid<GridImpl>::map(Matrix<double>& x) {
     #pragma omp for
 #endif
     for (unsigned i = 0; i < x.rows(); ++i) {
-      getNeighbours(x.rowSlice(i), weights, neighbours);
-
-      // linear interpolation
-      for (int d = static_cast<int>(dimDomain())-1; d >= 0; --d) {
-        for (int p = 0; p < (1 << d); ++p) {
-          for (int v = 0; v < static_cast<int>(dimCodomain()); ++v) {
-            neighbours[p*dimCodomain() + v] = neighbours[p*dimCodomain() + v] * (1.0-weights[d]) + neighbours[((1 << d) + p)*dimCodomain() + v] * weights[d];
+      if (m_interpolationType == Linear) {
+        getNeighbours(x.rowSlice(i), weights, neighbours);
+
+        // linear interpolation
+        for (int d = static_cast<int>(dimDomain())-1; d >= 0; --d) {
+          for (int p = 0; p < (1 << d); ++p) {
+            for (int v = 0; v < static_cast<int>(dimCodomain()); ++v) {
+              neighbours[p*dimCodomain() + v] = neighbours[p*dimCodomain() + v] * (1.0-weights[d]) + neighbours[((1 << d) + p)*dimCodomain() + v] * weights[d];
+            }
           }
         }
+      } else {
+        getNearestNeighbour(x.rowSlice(i), neighbours);
       }
 
       for (int v = 0; v < static_cast<int>(dimCodomain()); ++v) {
@@ -95,6 +114,21 @@ Matrix<double> Grid<GridImpl>::map(Matrix<double>& x) {
 
   return y;
 }
+
+template<typename GridImpl>
+void Grid<GridImpl>::setInterpolationType(std::string const& interpolationType) {
+  std::string iType = interpolationType;
+  std::transform(iType.begin(), iType.end(), iType.begin(), ::tolower);
+  if (iType == "nearest") {
+    setInterpolationType(Nearest);
+  } else if (iType == "linear") {
+    setInterpolationType(Linear);
+  } else {
+    std::stringstream ss;
+    ss << "Invalid interpolation type " << interpolationType << ".";
+    throw std::invalid_argument(ss.str());
+  }
+}
 }
 
 

include/easi/component/SCECFile.h

@@ -53,6 +53,9 @@ class SCECFile : public Grid<SCECFile> {
 
   void setMap(std::set<std::string> const& in, std::string const& fileName);
 
+  void getNearestNeighbour(Slice<double> const& x, double* buffer) {
+    m_grid->getNearestNeighbour(x, buffer);
+  }
   void getNeighbours(Slice<double> const& x, double* weights, double* buffer) {
     m_grid->getNeighbours(x, weights, buffer);
   }

include/easi/parser/YAMLComponentParsers.h

@@ -208,24 +208,34 @@ void parse<PolynomialMap>(PolynomialMap* component, YAML::Node const& node, std:
   parse<Map>(component, node, in, parser);
 }
 
+template<typename GridImpl>
+void parseGrid(Grid<GridImpl>* component, YAML::Node const& node, std::set<std::string> const& in, YAMLAbstractParser* parser) {
+  if (node["interpolation"]) {
+    std::string interpolation = node["interpolation"].as<std::string>();
+    component->setInterpolationType(interpolation);
+  }
+
+  parse<Map>(component, node, in, parser);
+}
+
 template<>
 void parse<SCECFile>(SCECFile* component, YAML::Node const& node, std::set<std::string> const& in, YAMLAbstractParser* parser) {
   checkType(node, "file", {YAML::NodeType::Scalar});
 
   std::string fileName = node["file"].as<std::string>();
   component->setMap(in, fileName);
 
-  parse<Map>(component, node, in, parser);
+  parseGrid<SCECFile>(component, node, in, parser);
 }
 
 #ifdef USE_ASAGI
 template<>
 void parse<ASAGI>(ASAGI* component, YAML::Node const& node, std::set<std::string> const& in, YAMLAbstractParser* parser) {
   checkType(node, "file", {YAML::NodeType::Scalar});
   checkType(node, "parameters", {YAML::NodeType::Sequence});
-  
+
   auto parameters = node["parameters"].as<std::vector<std::string>>();
-  
+
   std::string varName = "data";
   if (node["var"]) {
     varName = node["var"].as<std::string>();
@@ -236,7 +246,7 @@ void parse<ASAGI>(ASAGI* component, YAML::Node const& node, std::set<std::string
 
   component->setGrid(in, parameters, grid, parser->asagiReader()->numberOfThreads());
 
-  parse<Map>(component, node, in, parser);
+  parseGrid<ASAGI>(component, node, in, parser);
 }
 #endif
 

include/easi/util/RegularGrid.h

@@ -54,6 +54,7 @@ class RegularGrid {
   void setVolume(double const* min, double const* max);
   double* operator()(unsigned const* index);
 
+  void getNearestNeighbour(Slice<double> const& x, double* buffer);
   void getNeighbours(Slice<double> const& x, double* weights, double* buffer);
 
 private:
@@ -98,6 +99,27 @@ double* RegularGrid::operator()(unsigned const* index) {
   return m_values + m_numValues * idx;
 }
 
+void RegularGrid::getNearestNeighbour(Slice<double> const& x, double* buffer) {
+  assert(x.size() == m_dimensions);
+
+  unsigned idx[MaxDimensions];
+  for (unsigned d = 0; d < m_dimensions; ++d) {
+    if (x(d) < m_min[d]) {
+      idx[d] = 0;
+    } else if (x(d) >= m_max[d]) {
+      idx[d] = m_num[d]-1;
+    } else {
+      double xn = (x(d) - m_min[d]) / m_delta[d];
+      idx[d] = std::round(xn);
+    }
+  }
+
+  double* values = operator()(idx);
+  for (unsigned v = 0; v < m_numValues; ++v) {
+    buffer[v] = values[v];
+  }
+}
+
 void RegularGrid::getNeighbours(Slice<double> const& x, double* weights, double* buffer) {
   assert(x.size() == m_dimensions);
 

tests/101_asagi_nearest.cpp

@@ -0,0 +1,80 @@
+/**
+ * @file
+ * This file is part of SeisSol.
+ *
+ * @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.)
+ *
+ * @section LICENSE
+ * Copyright (c) 2017, SeisSol Group
+ * All rights reserved.
+ * 
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 
+ * 1. Redistributions of source code must retain the above copyright notice,
+ *    this list of conditions and the following disclaimer.
+ * 
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ * 
+ * 3. Neither the name of the copyright holder nor the names of its
+ *    contributors may be used to endorse or promote products derived from this
+ *    software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ *
+ * @section DESCRIPTION
+ **/
+
+#include <mpi.h>
+#include "easitest.h"
+
+int main(int argc, char** argv) {
+  assert(argc == 2);
+
+  MPI_Init(&argc, &argv);
+
+  easi::Query query = createQuery<3>({
+      {1, {  2.0,  -1.0, 0.0}},
+      {1, {  3.0,   0.0, 0.0}},
+      {1, {  2.5,  -0.5, 0.0}},
+      {1, { 2.75, -0.25, 0.0}},
+      {1, {  2.6, -0.25, 0.0}}
+    });
+  auto material = elasticModel(argv[1], query);
+
+  assert(equal(material[0].rho,      1.0));
+  assert(equal(material[0].mu,      -1.0));
+  assert(equal(material[0].lambda,  10.0));
+
+  assert(equal(material[1].rho,      6.0));
+  assert(equal(material[1].mu,      -6.0));
+  assert(equal(material[1].lambda,  15.0));
+
+  assert(equal(material[2].rho,      5.0));
+  assert(equal(material[2].mu,      -5.0));
+  assert(equal(material[2].lambda,  14.0));
+
+  assert(equal(material[3].rho,      6.0));
+  assert(equal(material[3].mu,      -6.0));
+  assert(equal(material[3].lambda,  15.0));
+
+  assert(equal(material[4].rho,      5.0));
+  assert(equal(material[4].mu,      -5.0));
+  assert(equal(material[4].lambda,  14.0));
+
+  MPI_Finalize();
+
+  return 0;
+}