GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BLI2DNonUnifObjectGrid.h
Go to the documentation of this file.
1 #ifndef BLI2DNONUNIF_OBJECT_GRID_H_
2 #define BLI2DNONUNIF_OBJECT_GRID_H_
3 
4 namespace genie {
5 
6 //____________________________________________________________________________
7 /*!
8 
9 \brief A class template that performs bilinear interpolation on a
10  non-uniform grid with an implementation similar to that of
11  genie::BLI2DNonUnifGrid.
12 
13 \details The main differences between this class template and
14  genie::BLI2DNonUnifGrid are
15 
16  * Values for the z coordinate can be any arbitrary object Object
17  that implements the member functions operator*(double),
18  operator*(const Object&), and operator+(const Object&)
19 
20  * Rather than C-style arrays, the grid values are accessed via
21  pointers to std::vector objects
22 
23  * Upper and lower bounds on the grid are found using
24  std::lower_bound() rather than a manual linear search
25 
26  * The genie::BLI2DNonUnifGrid object does not take ownership of the
27  grid vectors, which must be stored elsewhere
28 
29 \tparam ZObject Type of the object describing each z coordinate
30 \tparam IndexType Type to use when computing indices in the vectors
31 \tparam XType Type used to represent x coordinates
32 \tparam YType Type used to represent y coordinates
33 
34 \todo Think about how to have this class inherit from the other BLI2D
35  classes
36 
37 \author Steven Gardiner <gardiner \at fnal.gov>
38  Fermi National Accelerator Laboratory
39 
40 \created August 23, 2018
41 
42 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
43  For the full text of the license visit http://copyright.genie-mc.org
44  or see $GENIE/LICENSE
45 */
46 //____________________________________________________________________________
47 
48 template<typename ZObject, typename IndexType = int, typename XType = double,
49  typename YType = double> class BLI2DNonUnifObjectGrid
50 {
51  public:
52 
53  /// \param[in] X Pointer to a vector of x coordinates
54  /// \param[in] Y Pointer to a vector of y coordinates
55  /// \param[in] Z Pointer to a vector of z coordinates
56  /// \param[in] extrapolate Whether to allow bilinear extrapolation (true)
57  /// or to use the grid endpoints (false) when evaluating z values outside
58  /// of the grid
59  BLI2DNonUnifObjectGrid(const std::vector<XType>* X,
60  const std::vector<YType>* Y, const std::vector<ZObject>* Z,
61  bool extrapolate = false) : fX(X), fY(Y), fZ(Z), fExtrapolate(extrapolate)
62  {}
63 
64  /// Retrieve the minimum x value
65  inline XType x_min() const { return fX->front(); }
66 
67  /// Retrieve the maximum x value
68  inline XType x_max() const { return fX->back(); }
69 
70  /// Retrieve the minimum y value
71  inline YType y_min() const { return fY->front(); }
72 
73  /// Retrieve the maximum y value
74  inline YType y_max() const { return fY->back(); }
75 
76  /// Calculates the index in the vector of z coordinates that
77  /// corresponds to a given set of x and y indices
78  /// \param[in] ix Index of the desired grid point on the x axis
79  /// \param[in] iy Index of the desired grid point on the y axis
80  IndexType index_Z(IndexType ix, IndexType iy) const {
81  IndexType num_y_points = fY->size();
82 
83  return (num_y_points * ix) + iy;
84  }
85 
86  /// Uses bilinear interpolation to compute the z coordinate (represented
87  /// by an object of type ZObject) that corresponds to the given x and y
88  /// coordinates
89  ZObject interpolate(double x, double y) const
90  {
91  IndexType ix_lo, ix_hi, iy_lo, iy_hi;
92 
93  // For points outside the grid, evaluate the function at the end points
94  // unless extrapolation is enabled.
95  XType evalx = x;
96  if (!fExtrapolate) {
97  evalx = std::min(x, x_max());
98  evalx = std::max(evalx, x_min());
99  }
100 
101  YType evaly = y;
102  if (!fExtrapolate) {
103  evaly = std::min(y, y_max());
104  evaly = std::max(evaly, y_min());
105  }
106 
107  // Find the indices of the grid points on either side of the
108  // desired x and y values. If the desired point is outside of
109  // the x or y grid limits, get the indices of the two closest
110  // grid points to use for possible extrapolation.
111  get_bound_indices(fX, evalx, ix_lo, ix_hi);
112  get_bound_indices(fY, evaly, iy_lo, iy_hi);
113 
114  // Get the x and y values corresponding to the lower (x1, y1) and
115  // upper (x2, y2) bounds found previously
116  XType x1 = fX->at( ix_lo );
117  XType x2 = fX->at( ix_hi );
118  YType y1 = fY->at( iy_lo );
119  YType y2 = fY->at( iy_hi );
120 
121  // Retrieve the z values corresponding to each of the four locations
122  // that will be used for the bilinear interpolation
123  const ZObject& z11 = fZ->at( this->index_Z(ix_lo, iy_lo) );
124  const ZObject& z21 = fZ->at( this->index_Z(ix_hi, iy_lo) );
125  const ZObject& z12 = fZ->at( this->index_Z(ix_lo, iy_hi) );
126  const ZObject& z22 = fZ->at( this->index_Z(ix_hi, iy_hi) );
127 
128  // Perform the interpolation (first y, then x)
129  ZObject z1 = z11 * (y2-evaly)/(y2-y1) + z12 * (evaly-y1)/(y2-y1);
130  ZObject z2 = z21 * (y2-evaly)/(y2-y1) + z22 * (evaly-y1)/(y2-y1);
131  ZObject z = z1 * (x2-evalx)/(x2-x1) + z2 * (evalx-x1)/(x2-x1);
132 
133  return z;
134  }
135 
136 protected:
137 
138  const std::vector<XType>* fX; ///< Pointer to the vector of x coordinates
139  const std::vector<YType>* fY; ///< Pointer to the vector of y coordinates
140 
141  /// Pointer to the vector of z coordinate objects
142  const std::vector<ZObject>* fZ;
143 
144  /// Whether to allow bilinear extrapolation (true) or to compute z values for
145  /// x and coordinates outside of the grid using the grid endpoints (false)
147 
148  /// Determines the indices for the two gridpoints surrounding a requested
149  /// x or y coordinate. If the x or y coordinate is outside of the grid,
150  /// this function returns the two closest grid points.
151  /// \param[in] vec A vector of grid point coordinates
152  /// \param[in] val The requested x or y coordinate
153  /// \param[out] lower_index The index of the closest grid point less than or
154  /// equal to the requested value, or the lower of the two nearest grid points
155  /// if the value falls outside of the grid
156  /// \param[out] upper_index The index of the closest grid point greater than
157  /// the requested value, or the higher of the two nearest grid points
158  /// if the value falls outside of the grid
159  /// \return true if the requested value is within the grid, or false
160  /// otherwise
161  template <typename Type> bool get_bound_indices(const std::vector<Type>* vec,
162  Type val, int& lower_index, int& upper_index) const
163  {
164  /// \todo Check that the vector contains at least two entries
165 
166  bool within = true;
167 
168  typedef typename std::vector<Type>::const_iterator Iterator;
169  Iterator begin = vec->begin();
170  Iterator end = vec->end();
171 
172  // std::lower_bound returns an iterator to the first element of the
173  // container which is not less than the supplied value
174  Iterator not_less_point = std::lower_bound(begin, end, val);
175 
176  Iterator lower_point;
177 
178  // Check whether the requested grid point is within the grid limits
179  if (not_less_point == begin) {
180  lower_point = begin;
181  // first element of vec > val
182  if (*begin != val) within = false;
183  }
184  else if (not_less_point == end) {
185  // last element of vec < val
186  within = false;
187  lower_point = end - 2;
188  }
189  else {
190  // x is within the grid limits
191  lower_point = not_less_point - 1;
192  }
193 
194  lower_index = std::distance(begin, lower_point);
195  upper_index = lower_index + 1;
196 
197  return within;
198  }
199 
200 }; // template class BLI2DNonUnifObjectGrid
201 
202 } // namespace genie
203 #endif
IndexType index_Z(IndexType ix, IndexType iy) const
YType y_max() const
Retrieve the maximum y value.
const std::vector< ZObject > * fZ
Pointer to the vector of z coordinate objects.
YType y_min() const
Retrieve the minimum y value.
const std::vector< XType > * fX
Pointer to the vector of x coordinates.
XType x_max() const
Retrieve the maximum x value.
A class template that performs bilinear interpolation on a non-uniform grid with an implementation si...
const std::vector< YType > * fY
Pointer to the vector of y coordinates.
XType x_min() const
Retrieve the minimum x value.
bool get_bound_indices(const std::vector< Type > *vec, Type val, int &lower_index, int &upper_index) const
BLI2DNonUnifObjectGrid(const std::vector< XType > *X, const std::vector< YType > *Y, const std::vector< ZObject > *Z, bool extrapolate=false)
ZObject interpolate(double x, double y) const