interp2_neigh.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_INTERP2_NEIGH_H
24 #define O2SCL_INTERP2_NEIGH_H
25 
26 /** \file interp2_neigh.h
27  \brief File defining \ref o2scl::interp2_neigh
28 */
29 
30 #include <iostream>
31 #include <string>
32 #include <cmath>
33 
34 #include <o2scl/err_hnd.h>
35 
36 #ifndef DOXYGEN_NO_O2NS
37 namespace o2scl {
38 #endif
39 
40  /** \brief Nearest-neighbor interpolation in two dimensions
41 
42  This class performs nearest-neighbor interpolation when the data
43  points are not arranged in a specified order (i.e. not on a
44  grid). For a set of data \f$ {x_i,y_i,f_i} \f$, the value of \f$
45  f \f$ is predicted given a new value of x and y. Distances are
46  determined with
47  \f[
48  d_{ij} = \sqrt{\left(\frac{x_i-x_j}{\Delta x}\right)^2 +
49  \left(\frac{y_i-y_j}{\Delta y}\right)^2}
50  \f]
51  The values \f$ \Delta_x \f$ and \f$ \Delta_y \f$ are specified
52  in \ref x_scale and \ref y_scale, respectively. If these values
53  are negative (the default) then they are computed with \f$
54  \Delta x = x_{\mathrm{max}}-x_{\mathrm{min}} \f$ and \f$ \Delta
55  y = y_{\mathrm{max}}-y_{\mathrm{min}} \f$ .
56 
57  This class stores pointers to the data, not a copy. The data can
58  be changed between interpolations without an additional call to
59  \ref set_data(), but the scales may need to be recomputed with
60  \ref compute_scale().
61 
62  The vector type can be any type with a suitably defined \c
63  operator[].
64 
65  \note This class operates by performing a \f$ {\cal O}(N) \f$
66  brute-force search to find the closest points.
67 
68  \future Make a parent class for this and \ref o2scl::interp2_planar.
69 
70  \future Maybe interpm_idw subsumes this functionality and
71  makes this class obsolete? Or is this specialization
72  particularly efficient?
73  */
74  template<class vec_t> class interp2_neigh {
75 
76  protected:
77 
78  /// The scale in the x direction
79  double dx;
80 
81  /// The scale in the y direction
82  double dy;
83 
84  public:
85 
88 
89  interp2_neigh() {
90  data_set=false;
91  x_scale=-1.0;
92  y_scale=-1.0;
93  dx=0.0;
94  dy=0.0;
95  }
96 
97  /// The user-specified x scale (default -1)
98  double x_scale;
99 
100  /// The user-specified y scale (default -1)
101  double y_scale;
102 
103  /// Find scaling
104  void compute_scale() {
105  if (x_scale<0.0) {
106  double minx=(*ux)[0], maxx=(*ux)[0];
107  for(size_t i=1;i<np;i++) {
108  if ((*ux)[i]<minx) minx=(*ux)[i];
109  if ((*ux)[i]>maxx) maxx=(*ux)[i];
110  }
111  dx=maxx-minx;
112  } else {
113  dx=x_scale;
114  }
115  if (y_scale<0.0) {
116  double miny=(*uy)[0], maxy=(*uy)[0];
117  for(size_t i=1;i<np;i++) {
118  if ((*uy)[i]<miny) miny=(*uy)[i];
119  if ((*uy)[i]>maxy) maxy=(*uy)[i];
120  }
121  dy=maxy-miny;
122  } else {
123  dy=y_scale;
124  }
125 
126  if (dx<=0.0 || dy<=0.0) {
127  O2SCL_ERR("No scale in interp2_planar::set_data().",exc_einval);
128  }
129 
130  return;
131  }
132 
133  /** \brief Initialize the data for the neigh interpolation
134 
135  This function will call the error handler if \c n_points
136  is zero.
137  */
138  void set_data(size_t n_points, vec_t &x, vec_t &y, vec_t &f) {
139  if (n_points<1) {
140  O2SCL_ERR2("Must provide at least one point in ",
141  "interp2_neigh::set_data()",exc_efailed);
142  }
143  np=n_points;
144  ux=&x;
145  uy=&y;
146  uf=&f;
147  data_set=true;
148 
149  compute_scale();
150 
151  return;
152  }
153 
154  /** \brief Perform the interpolation
155  */
156  double eval(double x, double y) const {
157  double x1, y1, f;
158  size_t i1;
159  eval_point(x,y,f,i1,x1,y1);
160  return f;
161  }
162 
163  /** \brief Perform the interpolation
164  */
165  double operator()(double x, double y) const {
166  return eval(x,y);
167  }
168 
169  /** \brief Perform the planar interpolation using the first two
170  elements of \c v as input
171  */
172  template<class vec2_t> double operator()(vec2_t &v) const {
173  return eval(v[0],v[1]);
174  }
175 
176  /** \brief Interpolation returning the closest point
177 
178  This function interpolates \c x and \c y into the data
179  returning \c f. It also returns the closest x- and y-values
180  found.
181  */
182  void eval_point(double x, double y, double &f,
183  size_t &i1, double &x1, double &y1) const {
184 
185  if (data_set==false) {
186  O2SCL_ERR("Data not set in interp_planar::eval_points().",
187  exc_einval);
188  }
189 
190  // Exhaustively search the data
191  i1=0;
192  double dist_min=pow((x-(*ux)[i1])/dx,2.0)+pow((y-(*uy)[i1])/dy,2.0);
193  for(size_t index=1;index<np;index++) {
194  double dist=pow((x-(*ux)[index])/dx,2.0)+pow((y-(*uy)[index])/dy,2.0);
195  if (dist<dist_min) {
196  i1=index;
197  dist_min=dist;
198  }
199  }
200 
201  // Return the function value
202 
203  f=(*uf)[i1];
204 
205  return;
206  }
207 
208 #ifndef DOXYGEN_INTERNAL
209 
210  protected:
211 
212  /// The number of points
213  size_t np;
214  /// The x-values
215  vec_t *ux;
216  /// The y-values
217  vec_t *uy;
218  /// The f-values
219  vec_t *uf;
220  /// True if the data has been specified
221  bool data_set;
222 
223 #endif
224 
225  };
226 
227 #ifndef DOXYGEN_NO_O2NS
228 }
229 #endif
230 
231 #endif
232 
233 
234 
Nearest-neighbor interpolation in two dimensions.
Definition: interp2_neigh.h:74
double y_scale
The user-specified y scale (default -1)
vec_t * ux
The x-values.
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
double x_scale
The user-specified x scale (default -1)
Definition: interp2_neigh.h:98
double operator()(double x, double y) const
Perform the interpolation.
vec_t * uf
The f-values.
bool data_set
True if the data has been specified.
invalid argument supplied by user
Definition: err_hnd.h:59
double dx
The scale in the x direction.
Definition: interp2_neigh.h:79
double eval(double x, double y) const
Perform the interpolation.
generic failure
Definition: err_hnd.h:61
double operator()(vec2_t &v) const
Perform the planar interpolation using the first two elements of v as input.
vec_t * uy
The y-values.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
size_t np
The number of points.
void set_data(size_t n_points, vec_t &x, vec_t &y, vec_t &f)
Initialize the data for the neigh interpolation.
double dy
The scale in the y direction.
Definition: interp2_neigh.h:82
static const double x1[5]
Definition: inte_qng_gsl.h:48
void compute_scale()
Find scaling.
void eval_point(double x, double y, double &f, size_t &i1, double &x1, double &y1) const
Interpolation returning the closest point.

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).