MeshKernel
BilinearInterpolationOnGriddedSamples.hpp
1 //---- GPL ---------------------------------------------------------------------
2 //
3 // Copyright (C) Stichting Deltares, 2011-2021.
4 //
5 // This program is free software: you can redistribute it and/or modify
6 // it under the terms of the GNU General Public License as published by
7 // the Free Software Foundation version 3.
8 //
9 // This program is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU General Public License for more details.
13 //
14 // You should have received a copy of the GNU General Public License
15 // along with this program. If not, see <http://www.gnu.org/licenses/>.
16 //
17 // contact: delft3d.support@deltares.nl
18 // Stichting Deltares
19 // P.O. Box 177
20 // 2600 MH Delft, The Netherlands
21 //
22 // All indications and logos of, and references to, "Delft3D" and "Deltares"
23 // are registered trademarks of Stichting Deltares, and remain the property of
24 // Stichting Deltares. All rights reserved.
25 //
26 //------------------------------------------------------------------------------
27 
28 #pragma once
29 
30 #include <MeshKernel/Mesh2D.hpp>
31 #include <MeshKernel/MeshInterpolation.hpp>
32 
33 #include <concepts>
34 #include <span>
35 
36 namespace meshkernel
37 {
38 
40  template <typename T>
41  concept InterpolatableType = std::same_as<T, short> ||
42  std::same_as<T, float> ||
43  std::same_as<T, double> ||
44  std::same_as<T, int>;
45 
47  template <InterpolatableType T>
48  class BilinearInterpolationOnGriddedSamples : public MeshInterpolation
49  {
50  public:
58  BilinearInterpolationOnGriddedSamples(const Mesh2D& mesh,
59  UInt numXCoord,
60  UInt numYCoord,
61  const Point& origin,
62  double cellSize,
63  std::span<T const> values);
64 
70  BilinearInterpolationOnGriddedSamples(const Mesh2D& mesh,
71  std::span<double const> xCoordinates,
72  std::span<double const> yCoordinates,
73  // const std::vector<double>& xCoordinates,
74  // const std::vector<double>& yCoordinates,
75  std::span<T const> values);
76 
78  void Compute() override;
79 
80  private:
84  [[nodiscard]] double Interpolation(const Point& point) const;
85 
89  [[nodiscard]] double GetFractionalNumberOfColumns(const Point& point) const;
90 
94  [[nodiscard]] double GetFractionalNumberOfRows(const Point& point) const;
95 
98  [[nodiscard]] double GetGriddedValue(UInt columnIndex, UInt rowIndex) const;
99 
100  const Mesh2D& m_mesh;
101  UInt m_numXCoord;
102  UInt m_numYCoord;
103  Point m_origin;
104  double m_cellSize = 0.0;
105 
106  std::span<double const> m_xCoordinates;
107  std::span<double const> m_yCoordinates;
108  std::span<T const> m_values;
109  bool m_isCellSizeConstant;
110  };
111 
112  template <InterpolatableType T>
113  BilinearInterpolationOnGriddedSamples<T>::BilinearInterpolationOnGriddedSamples(const Mesh2D& mesh,
114  UInt numXCoord,
115  UInt numYCoord,
116  const Point& origin,
117  double cellSize,
118  std::span<T const> values)
119  : m_mesh(mesh),
120  m_numXCoord(numXCoord),
121  m_numYCoord(numYCoord),
122  m_origin(origin),
123  m_cellSize(cellSize),
124  m_values(values),
125  m_isCellSizeConstant{true}
126  {
127  }
128 
129  template <InterpolatableType T>
130  BilinearInterpolationOnGriddedSamples<T>::BilinearInterpolationOnGriddedSamples(const Mesh2D& mesh,
131  std::span<double const> xCoordinates,
132  std::span<double const> yCoordinates,
133  std::span<T const> values)
134  : m_mesh(mesh),
135  m_numXCoord(static_cast<UInt>(xCoordinates.size())),
136  m_numYCoord(static_cast<UInt>(yCoordinates.size())),
137  m_xCoordinates(xCoordinates),
138  m_yCoordinates(yCoordinates),
139  m_values{values},
140  m_isCellSizeConstant(false)
141  {
142  }
143 
144  template <InterpolatableType T>
145  void BilinearInterpolationOnGriddedSamples<T>::Compute()
146  {
147  const auto numNodes = m_mesh.GetNumNodes();
148  const auto numEdges = m_mesh.GetNumEdges();
149  const auto numFaces = m_mesh.GetNumFaces();
150 
151  m_nodeResults.resize(numNodes);
152  std::ranges::fill(m_nodeResults, constants::missing::doubleValue);
153  for (UInt n = 0; n < numNodes; ++n)
154  {
155  const auto node = m_mesh.Node(n);
156  m_nodeResults[n] = Interpolation(node);
157  }
158 
159  m_edgeResults.resize(numEdges);
160  std::ranges::fill(m_edgeResults, constants::missing::doubleValue);
161  for (UInt e = 0; e < numEdges; ++e)
162  {
163  const auto& [first, second] = m_mesh.GetEdge(e);
164  m_edgeResults[e] = 0.5 * (m_nodeResults[first] + m_nodeResults[second]);
165  }
166 
167  m_faceResults.resize(numFaces, constants::missing::doubleValue);
168  std::ranges::fill(m_faceResults, constants::missing::doubleValue);
169  for (UInt f = 0; f < numFaces; ++f)
170  {
171  m_faceResults[f] = Interpolation(m_mesh.m_facesMassCenters[f]);
172  }
173  }
174 
175  template <InterpolatableType T>
176  double BilinearInterpolationOnGriddedSamples<T>::Interpolation(const Point& point) const
177  {
178 
179  double fractionalColumnIndex = GetFractionalNumberOfColumns(point);
180  double fractionalRowIndex = GetFractionalNumberOfRows(point);
181 
182  double columnIndexTmp;
183  fractionalColumnIndex = std::modf(fractionalColumnIndex, &columnIndexTmp);
184 
185  double rowIndexTmp;
186  fractionalRowIndex = std::modf(fractionalRowIndex, &rowIndexTmp);
187 
188  if (columnIndexTmp < 0 || rowIndexTmp < 0)
189  {
190  return constants::missing::doubleValue;
191  }
192 
193  UInt const columnIndex = static_cast<UInt>(columnIndexTmp);
194  UInt const rowIndex = static_cast<UInt>(rowIndexTmp);
195 
196  if (columnIndex + 1 >= m_numXCoord || rowIndex + 1 >= m_numYCoord)
197  {
198  return constants::missing::doubleValue;
199  }
200 
201  const auto result = fractionalColumnIndex * fractionalRowIndex * GetGriddedValue(columnIndex + 1, rowIndex + 1) +
202  (1.0 - fractionalColumnIndex) * fractionalRowIndex * GetGriddedValue(columnIndex, rowIndex + 1) +
203  (1.0 - fractionalColumnIndex) * (1.0 - fractionalRowIndex) * GetGriddedValue(columnIndex, rowIndex) +
204  fractionalColumnIndex * (1.0 - fractionalRowIndex) * GetGriddedValue(columnIndex + 1, rowIndex);
205  return result;
206  }
207 
208  template <InterpolatableType T>
209  double BilinearInterpolationOnGriddedSamples<T>::GetFractionalNumberOfColumns(const Point& point) const
210  {
211  if (m_isCellSizeConstant)
212  {
213  return (point.x - m_origin.x) / m_cellSize;
214  }
215  double result = constants::missing::doubleValue;
216  if (m_xCoordinates.size() < 2)
217  {
218  return result;
219  }
220  for (UInt i = 0; i < m_xCoordinates.size() - 1; ++i)
221  {
222 
223  if (point.x >= m_xCoordinates[i] && point.x < m_xCoordinates[i + 1])
224  {
225  const double dx = m_xCoordinates[i + 1] - m_xCoordinates[i];
226  result = static_cast<double>(i) + (point.x - m_xCoordinates[i]) / dx;
227  break;
228  }
229  }
230  return result;
231  }
232 
233  template <InterpolatableType T>
234  double BilinearInterpolationOnGriddedSamples<T>::GetFractionalNumberOfRows(const Point& point) const
235  {
236  if (m_isCellSizeConstant)
237  {
238  return (point.y - m_origin.y) / m_cellSize;
239  }
240 
241  double result = constants::missing::doubleValue;
242  if (m_yCoordinates.size() < 2)
243  {
244  return result;
245  }
246 
247  for (UInt i = 0; i < m_yCoordinates.size() - 1; ++i)
248  {
249 
250  if (point.y >= m_yCoordinates[i] && point.y < m_yCoordinates[i + 1])
251  {
252  const double dy = m_yCoordinates[i + 1] - m_yCoordinates[i];
253  result = static_cast<double>(i) + (point.y - m_yCoordinates[i]) / dy;
254  break;
255  }
256  }
257  return result;
258  }
259 
260  template <InterpolatableType T>
261  double BilinearInterpolationOnGriddedSamples<T>::GetGriddedValue(UInt columnIndex, UInt rowIndex) const
262  {
263  const auto index = rowIndex * m_numXCoord + columnIndex;
264  return static_cast<double>(m_values[index]);
265  }
266 
267 } // namespace meshkernel
meshkernel::MeshInterpolation
Interface for interpolation methods.
Definition: MeshInterpolation.hpp:42
meshkernel::Mesh2D
A class derived from Mesh, which describes unstructures 2d meshes.
Definition: Mesh2D.hpp:55
meshkernel::Point
A struct describing a point in a two-dimensional space.
Definition: Point.hpp:40
meshkernel
Contains the logic of the C++ static library.
Definition: AveragingInterpolation.hpp:36
meshkernel::InterpolatableType
concept InterpolatableType
Defines the iterpolatable data types.
Definition: BilinearInterpolationOnGriddedSamples.hpp:41
meshkernel::UInt
std::uint32_t UInt
Integer type used when indexing mesh graph entities.
Definition: Definitions.hpp:38
meshkernel::BilinearInterpolationOnGriddedSamples
A class for performing bilinear interpolation on gridded samples.
Definition: BilinearInterpolationOnGriddedSamples.hpp:48
meshkernel::BilinearInterpolationOnGriddedSamples::Compute
void Compute() override
Compute interpolation.
Definition: BilinearInterpolationOnGriddedSamples.hpp:145
meshkernel::BilinearInterpolationOnGriddedSamples::BilinearInterpolationOnGriddedSamples
BilinearInterpolationOnGriddedSamples(const Mesh2D &mesh, UInt numXCoord, UInt numYCoord, const Point &origin, double cellSize, std::span< T const > values)
Bilinear interpolation with constant cell size (faster because no linear search is performed for each...
Definition: BilinearInterpolationOnGriddedSamples.hpp:113
Generated by   doxygen 1.8.17