contour.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 /** \file contour.h
24  \brief File defining \ref o2scl::contour
25 */
26 #ifndef O2SCL_CONTOUR_H
27 #define O2SCL_CONTOUR_H
28 
29 #include <cmath>
30 
31 #include <gsl/gsl_math.h>
32 
33 #include <boost/numeric/ublas/vector.hpp>
34 #include <boost/numeric/ublas/matrix.hpp>
35 
36 #include <o2scl/interp.h>
37 #include <o2scl/uniform_grid.h>
38 
39 #ifndef DOXYGEN_NO_O2NS
40 namespace o2scl {
41 #endif
42 
43  /** \brief A contour line
44 
45  The contour lines generated by the \ref o2scl::contour class are
46  given as objects of this type.
47 
48  \future Make this a subclass of \ref o2scl::contour ?
49  (12/12/16 Maybe not, as it is a generic concept which can
50  be computed outside the \ref o2scl::contour class.)
51  */
52  class contour_line {
53 
54  public:
55 
56  /// The contour level
57  double level;
58  /// The line x coordinates
59  std::vector<double> x;
60  /// The line y coordinates
61  std::vector<double> y;
62 
63  /// Create an empty line
65  }
66 
67  /// Copy constructor
69  level=c.level;
70  x=c.x;
71  y=c.y;
72  }
73 
74  /// Copy constructor with operator=()
76  if (this==&c) return *this;
77  level=c.level;
78  x=c.x;
79  y=c.y;
80  return *this;
81  }
82 
83  };
84 
85  /** \brief Edges for the contour class
86 
87  The edge crossings generated by the \ref contour class are given
88  as objects of this type.
89 
90  The \ref status matrix contains one of four possible values
91  - 0 - empty (no edge)
92  - 1 - edge which has not yet been assigned to a contour
93  - 2 - edge assigned to contour point
94  - 3 - edge which has been designated as a contour endpoint
95 
96  The matrices returned by \ref contour are not square, their
97  size depends on whether or not they contain the "bottom edges"
98  or the "right edges".
99 
100  \future Make this a subclass of \ref o2scl::contour .
101  */
103 
104  public:
105 
110 
111  /// Edge status
112  ubmatrix_int status;
113  /// Edge values
114  ubmatrix values;
115 
116  /// Create an empty object
118  }
119 
120  /// Copy constructor
122  status=ec.status;
123  values=ec.values;
124  }
125 
126  /// Copy constructor with operator=()
128  if (this==&ec) return *this;
129  status=ec.status;
130  values=ec.values;
131  return *this;
132  }
133 
134  };
135 
136  /** \brief Calculate contour lines from a two-dimensional data set
137 
138  \b Basic \b Usage
139 
140  - Specify the data as a two-dimensional square grid of "heights"
141  with set_data().
142  - Specify the contour levels with set_levels().
143  - Compute the contours with calc_contours()
144 
145  The contours are generated as a series of x- and y-coordinates,
146  defining a line. If the contour is closed, then the first and
147  the last set of coordinates will be equal.
148 
149  The convention used by this class is that the first (row) index
150  of the matrix enumerates the x coordinate and that the second
151  (column) index enumerates the y coordinate. See the discussion
152  in the User's guide in the section called \ref rowcol_subsect.
153 
154  The data is copied by set_data(), so changing the data will not
155  change the contours unless set_data() is called again. The
156  functions set_levels() and calc_contours() can be called several
157  times for the same data without calling set_data() again.
158 
159  Note that in order to simplify the algorithm for computing
160  contour lines, the calc_contours() function will sometimes
161  adjust the user-specified contour levels slightly in order to
162  ensure that no contour line passes exactly through any data
163  point on the grid. The contours are adjusted by multiplying the
164  original contour level by 1 plus a small number (\f$ 10^{-8} \f$
165  by default), which is specified in \ref lev_adjust.
166 
167  Linear interpolation is used to decide whether or not a line
168  segment from the grid and a contour cross. This choice is
169  intentional, since (in addition to making the algorithm much
170  simpler) it is the user (and not this contour class) which is
171  likely best able to refine the data. In case a simple refinement
172  scheme is desired, the method regrid_data() is provided which
173  refines the internally stored data for any interpolation type.
174 
175  Since linear interpolation is used, the contour calculation
176  implicitly assumes that there is not more than one intersection
177  of any contour level with any line segment. For contours which
178  do not close inside the region of interest, the results will
179  always end at either the minimum or maximum values of the x or y
180  grid points (no extrapolation is ever done). Note also that the
181  points defining the contour are not necessarily equally spaced.
182  Two neighboring points will never be farther apart than the
183  distance across opposite corners of one cell in the grid.
184 
185  \b The \b Algorithm:
186 
187  This works by viewing the data as defining a square
188  two-dimensional grid. The function calc_contours() exhaustively
189  enumerates every line segment in the grid which involves a level
190  crossing and then organizes the points defined by the
191  intersection of a line segment with a level curve into a full
192  contour line.
193 
194  \future Copy constructor
195 
196  \future Improve the algorithm to ensure that no contour
197  line ends on an internal point. I am not sure the
198  best way to do this, but it could be done recursively
199  just by trying all possible links instead of just using
200  the one that minimizes the distance.
201 
202  \future Rewrite the code which adjusts the contour levels
203  to ensure contours don't go through the data to adjust the
204  internal copy of the data instead? This should be more
205  accurate because we're perturbing one point instead of
206  perturbing the entire line.
207 
208  \future Change nx and ny to size_t?
209 
210  \future It would be nice to have a function which creates a set
211  of closed regions to fill which represent the data. However,
212  this likely requires a completely new algorithm, because it's
213  not easy to simply close the contours already generated by the
214  calc_contours() function. There are, for example, several cases
215  which are difficult to handle, such as filling a region in
216  between several closed contours
217 
218  \comment
219 
220  6/28/12 - Thinking about how the algorithm could work.
221 
222  1) Start at the lowest contour level.
223 
224  2) How to create the enclosure in general if it's closed?
225  Just make sure that lines betwen closed region and boundary
226  don't go through any other contours of the same level.
227 
228  3) If one of the two enclosures would enclose a different
229  contour line with the same level, and choose the other
230  enclosure. If neither enclosure would enclose a contour line
231  with the same level, then only one enclosure encloses all other
232  contour lines, so that is the one to choose.
233 
234  4) Proceed to the next contour line of the same level, and
235  then go to the next contour line of the second lowest
236  contour level.
237 
238  \endcomment
239  */
240  class contour {
241 
242  public:
243 
247 
248  contour();
249 
250  ~contour();
251 
252  /// \name Basic usage
253  //@{
254 
255  /** \brief Set the data
256 
257  The types \c vec_t and \c mat_t can be any types which have \c
258  operator[] and \c operator[][] for array and matrix indexing.
259 
260  Note that this method copies all of the user-specified data to
261  local storage so that changes in the data after calling this
262  function will not be reflected in the contours that are
263  generated.
264  */
265  template<class vec_t, class mat_t>
266  void set_data(size_t sizex, size_t sizey, const vec_t &x_fun,
267  const vec_t &y_fun, const mat_t &udata) {
268 
269  if (sizex<2 || sizey<2) {
270  O2SCL_ERR("Not enough data (must be at least 2x2) in set_data().",
271  exc_einval);
272  }
273 
274  nx=sizex;
275  ny=sizey;
276  xfun.resize(nx);
277  yfun.resize(ny);
278  data.resize(nx,ny);
279  for(int i=0;i<nx;i++) xfun[i]=x_fun[i];
280  for(int i=0;i<ny;i++) yfun[i]=y_fun[i];
281  for(int i=0;i<nx;i++) {
282  for(int j=0;j<ny;j++) {
283  data(i,j)=udata(i,j);
284  }
285  }
286 
287  check_data();
288  return;
289  }
290 
291  /** \brief Set the data
292 
293  The type \c mat_t can be any type which has \c operator[][]
294  for matrix indexing.
295 
296  Note that this method copies all of the user-specified data to
297  local storage so that changes in the data after calling this
298  function will not be reflected in the contours that are
299  generated.
300  */
301  template<class mat_t>
302  void set_data(const uniform_grid<double> &ugx,
303  const uniform_grid<double> &ugy,
304  const mat_t &udata) {
305 
306  size_t sizex=ugx.get_npoints();
307  size_t sizey=ugy.get_npoints();
308 
309  if (sizex<2 || sizey<2) {
310  O2SCL_ERR("Not enough data (must be at least 2x2) in set_data().",
311  exc_einval);
312  }
313 
314  nx=sizex;
315  ny=sizey;
316  xfun.resize(nx);
317  yfun.resize(ny);
318  data.resize(nx,ny);
319  for(int i=0;i<nx;i++) xfun[i]=ugx[i];
320  for(int i=0;i<ny;i++) yfun[i]=ugy[i];
321  for(int i=0;i<nx;i++) {
322  for(int j=0;j<ny;j++) {
323  data(i,j)=udata(i,j);
324  }
325  }
326 
327  check_data();
328  return;
329  }
330 
331  /** \brief Set the contour levels
332 
333  This is separate from the function calc_contours() so that
334  the user can compute the contours for different data sets using
335  the same levels.
336  */
337  template<class vec_t> void set_levels(size_t nlevels, vec_t &ulevels) {
338  nlev=nlevels;
339  levels.resize(nlevels);
340  for(size_t i=0;i<nlevels;i++) {
341  levels[i]=ulevels[i];
342  }
343  levels_set=true;
344  return;
345  }
346 
347  /** \brief Calculate the contours
348 
349  \note There may be zero or more than one contour line for
350  each level, so the size of \c clines is not necessarily
351  equal to the number of levels specified in \ref set_levels().
352  */
353  void calc_contours(std::vector<contour_line> &clines);
354  //@}
355 
356  /// \name Regrid function
357  //@{
358  /** \brief Regrid the data
359 
360  Use interpolation to refine the data set. This can be called
361  before calc_contours() in order to attempt make the contour
362  levels smoother by providing a smaller grid size. If the
363  original number of data points is \f$
364  (\mathrm{nx},\mathrm{ny}) \f$, then the new number of data
365  points is
366  \f[
367  (\mathrm{xfact}~(\mathrm{nx}-1)+1,
368  \mathrm{yfact}~(\mathrm{ny}-1)+1)
369  \f]
370  The parameters \c xfact and \c yfact must both be larger than
371  zero and they cannot both be 1.
372  */
373  void regrid_data(size_t xfact, size_t yfact,
374  size_t interp_type=o2scl::itp_cspline);
375  //@}
376 
377  /// \name Obtain internal data
378  //@{
379  /** \brief Get the data
380 
381  This is useful to see how the data has changed after
382  a call to regrid_data().
383 
384  \future There is probably a better way than returning
385  pointers to the internal data.
386  */
387  void get_data(size_t &sizex, size_t &sizey, ubvector *&x_fun,
388  ubvector *&y_fun, ubmatrix *&udata) {
389  if (nx==0) {
390  O2SCL_ERR("Data not set in get_data().",exc_einval);
391  }
392  sizex=nx;
393  sizey=ny;
394  x_fun=&xfun;
395  y_fun=&yfun;
396  udata=&data;
397  return;
398  }
399 
400  /** \brief Return the edges for each contour level
401 
402  The size of \c y_edges and \c x_edges will both be equal to
403  the number of levels set by \ref set_levels().
404  */
405  void get_edges(std::vector<edge_crossings> &x_edges,
406  std::vector<edge_crossings> &y_edges) {
407  x_edges=xed;
408  y_edges=yed;
409  return;
410  }
411 
412  /// Print out the edges to cout
413  void print_edges_yhoriz(edge_crossings &xedges,
414  edge_crossings &yedges);
415 
416  /// Print out the edges to cout
417  void print_edges_xhoriz(edge_crossings &xedges,
418  edge_crossings &yedges);
419 
420  //@}
421 
422  /** \brief Verbosity parameter (default 0)
423 
424  If verbose is greater than 0, then adjustments to the contour
425  levels will be output and the contour lines will be output as
426  they are built up from the intersections. If verbose is
427  greater than 1, then all edges will be output. If verbose is
428  greater than 2, a keypress will be required after each contour
429  line is constructed.
430  */
431  int verbose;
432 
433  /// (default \f$ 10^{-8} \f$)
434  double lev_adjust;
435 
436  /** \brief If true, debug the functions which determine the
437  next point functions (default false)p
438  */
440 
441  /// \name Edge status
442  //@{
443  static const int empty=0;
444  static const int edge=1;
445  static const int contourp=2;
446  static const int endpoint=3;
447  //@}
448 
449 #ifndef DOXYGEN_INTERNAL
450 
451  protected:
452 
453  /// \name Edge direction
454  //@{
455  static const int dxdir=0;
456  static const int dydir=1;
457  //@}
458 
459  /// \name Edge found or not found
460  //@{
461  static const int efound=1;
462  static const int enot_found=0;
463  //@}
464 
465  /// \name User-specified data
466  //@{
467  int nx, ny;
468  ubvector xfun, yfun;
469  ubmatrix data;
470  //@}
471 
472  /// \name User-specified contour levels
473  //@{
474  int nlev;
475  ubvector levels;
476  bool levels_set;
477  //@}
478 
479  /// Right edge list
480  std::vector<edge_crossings> yed;
481  /// Bottom edge list
482  std::vector<edge_crossings> xed;
483 
484  /// Find next point starting from a point on a right edge
485  int find_next_point_y_direct(int j, int k, int &jnext, int &knext,
486  int &dir_next, int nsw,
487  edge_crossings &xedges,
488  edge_crossings &yedges);
489 
490  /// Find next point starting from a point on a bottom edge
491  int find_next_point_x_direct(int j, int k, int &jnext, int &knext,
492  int &dir_next, int nsw,
493  edge_crossings &xedges,
494  edge_crossings &yedges);
495 
496  /// Find all of the intersections of the edges with the contour level
497  void find_intersections(size_t ilev, double &level,
498  edge_crossings &xedges, edge_crossings &yedges);
499 
500  /// Interpolate all right edge crossings
501  void edges_in_y_direct(double level, interp<ubvector> &si,
502  edge_crossings &yedges);
503 
504  /// Interpolate all bottom edge crossings
505  void edges_in_x_direct(double level, interp<ubvector> &si,
506  edge_crossings &xedges);
507 
508  /// Create a contour line from a starting edge
509  void process_line(int j, int k, int dir, std::vector<double> &x,
510  std::vector<double> &y,
511  bool first, edge_crossings &xedges,
512  edge_crossings &yedges);
513 
514  /// Check to ensure the x- and y-arrays are monotonic
515  void check_data();
516 
517  private:
518 
519  contour(const contour &);
520  contour& operator=(const contour&);
521 
522 #endif
523 
524  };
525 
526 #ifndef DOXYGEN_NO_O2NS
527 }
528 #endif
529 
530 #endif
531 
532 
edge_crossings()
Create an empty object.
Definition: contour.h:117
Edges for the contour class.
Definition: contour.h:102
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
contour_line()
Create an empty line.
Definition: contour.h:64
edge_crossings(const edge_crossings &ec)
Copy constructor.
Definition: contour.h:121
invalid argument supplied by user
Definition: err_hnd.h:59
contour_line(const contour_line &c)
Copy constructor.
Definition: contour.h:68
void set_data(size_t sizex, size_t sizey, const vec_t &x_fun, const vec_t &y_fun, const mat_t &udata)
Set the data.
Definition: contour.h:266
size_t get_npoints() const
Get the number of points in the grid (always get_nbins()+1)
Definition: uniform_grid.h:201
Cubic spline for natural boundary conditions.
Definition: interp.h:72
void get_data(size_t &sizex, size_t &sizey, ubvector *&x_fun, ubvector *&y_fun, ubmatrix *&udata)
Get the data.
Definition: contour.h:387
std::vector< edge_crossings > xed
Bottom edge list.
Definition: contour.h:482
Calculate contour lines from a two-dimensional data set.
Definition: contour.h:240
ubmatrix_int status
Edge status.
Definition: contour.h:112
ubmatrix values
Edge values.
Definition: contour.h:114
std::vector< double > x
The line x coordinates.
Definition: contour.h:59
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
void set_data(const uniform_grid< double > &ugx, const uniform_grid< double > &ugy, const mat_t &udata)
Set the data.
Definition: contour.h:302
A contour line.
Definition: contour.h:52
double lev_adjust
(default )
Definition: contour.h:434
std::vector< edge_crossings > yed
Right edge list.
Definition: contour.h:480
int verbose
Verbosity parameter (default 0)
Definition: contour.h:431
A class representing a uniform linear or logarithmic grid.
Definition: uniform_grid.h:38
double level
The contour level.
Definition: contour.h:57
void set_levels(size_t nlevels, vec_t &ulevels)
Set the contour levels.
Definition: contour.h:337
edge_crossings & operator=(const edge_crossings &ec)
Copy constructor with operator=()
Definition: contour.h:127
contour_line & operator=(const contour_line &c)
Copy constructor with operator=()
Definition: contour.h:75
std::vector< double > y
The line y coordinates.
Definition: contour.h:61
bool debug_next_point
If true, debug the functions which determine the next point functions (default false)p.
Definition: contour.h:439
Interpolation class for general vectors.
Definition: interp.h:1558
void get_edges(std::vector< edge_crossings > &x_edges, std::vector< edge_crossings > &y_edges)
Return the edges for each contour level.
Definition: contour.h:405

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