Class contour (o2scl)

O2scl : Class List

class contour

Calculate contour lines from a two-dimensional data set.

Basic Usage

The contours are generated as a series of x- and y-coordinates, defining a line. If the contour is closed, then the first and the last set of coordinates will be equal.

The convention used by this class is that the first (row) index of the matrix enumerates the x coordinate and that the second (column) index enumerates the y coordinate. See the discussion in the User’s guide in the section called Rows and columns vs. x and y.

The data is copied by set_data(), so changing the data will not change the contours unless set_data() is called again. The functions set_levels() and calc_contours() can be called several times for the same data without calling set_data() again.

Note that in order to simplify the algorithm for computing contour lines, the calc_contours() function will sometimes adjust the user-specified contour levels slightly in order to ensure that no contour line passes exactly through any data point on the grid. The contours are adjusted by multiplying the original contour level by 1 plus a small number ( \( 10^{-8} \) by default), which is specified in lev_adjust.

Linear interpolation is used to decide whether or not a line segment from the grid and a contour cross. This choice is intentional, since (in addition to making the algorithm much simpler) it is the user (and not this contour class) which is likely best able to refine the data. In case a simple refinement scheme is desired, the method regrid_data() is provided which refines the internally stored data for any interpolation type.

Since linear interpolation is used, the contour calculation implicitly assumes that there is not more than one intersection of any contour level with any line segment. For contours which do not close inside the region of interest, the results will always end at either the minimum or maximum values of the x or y grid points (no extrapolation is ever done). Note also that the points defining the contour are not necessarily equally spaced. Two neighboring points will never be farther apart than the distance across opposite corners of one cell in the grid.

The Algorithm:

This works by viewing the data as defining a square two-dimensional grid. The function calc_contours() exhaustively enumerates every line segment in the grid which involves a level crossing and then organizes the points defined by the intersection of a line segment with a level curve into a full contour line.

Todo

Future: Copy constructor

Future: Improve the algorithm to ensure that no contour line ends on an internal point. I am not sure the best way to do this, but it could be done recursively just by trying all possible links instead of just using the one that minimizes the distance.

Future: Rewrite the code which adjusts the contour levels to ensure contours don’t go through the data to adjust the internal copy of the data instead? This should be more accurate because we’re perturbing one point instead of perturbing the entire line.

Future: Change nx and ny to size_t?

Future: It would be nice to have a function which creates a set of closed regions to fill which represent the data. However, this likely requires a completely new algorithm, because it’s not easy to simply close the contours already generated by the calc_contours() function. There are, for example, several cases which are difficult to handle, such as filling a region in between several closed contours

Obtain internal data

int verbose

Verbosity parameter (default 0)

If verbose is greater than 0, then adjustments to the contour levels will be output and the contour lines will be output as they are built up from the intersections. If verbose is greater than 1, then all edges will be output. If verbose is greater than 2, a keypress will be required after each contour line is constructed.

double lev_adjust

(default \( 10^{-8} \))

bool debug_next_point

If true, debug the functions which determine the next point functions (default false)p.

inline void get_data(size_t &sizex, size_t &sizey, ubvector *&x_fun, ubvector *&y_fun, ubmatrix *&udata)

Get the data.

This is useful to see how the data has changed after a call to regrid_data(). 

\verbatim embed:rst
.. todo:: 

   Future: There is probably a better way than returning
   pointers to the internal data.

inline void get_edges(std::vector<edge_crossings> &x_edges, std::vector<edge_crossings> &y_edges)

Return the edges for each contour level.

The size of y_edges and x_edges will both be equal to the number of levels set by set_levels().

void print_edges_yhoriz(edge_crossings &xedges, edge_crossings &yedges)

Print out the edges to cout.

void print_edges_xhoriz(edge_crossings &xedges, edge_crossings &yedges)

Print out the edges to cout.

Edge status

static const int empty = 0
static const int edge = 1
static const int contourp = 2
static const int endpoint = 3

Edge direction

static const int dxdir = 0
static const int dydir = 1

Edge found or not found

static const int efound = 1
static const int enot_found = 0

User-specified data

int nx
int ny
ubvector xfun
ubvector yfun
ubmatrix data

User-specified contour levels

int nlev
ubvector levels
bool levels_set
std::vector<edge_crossings> yed

Right edge list.

std::vector<edge_crossings> xed

Bottom edge list.

int find_next_point_y_direct(int j, int k, int &jnext, int &knext, int &dir_next, int nsw, edge_crossings &xedges, edge_crossings &yedges)

Find next point starting from a point on a right edge.

int find_next_point_x_direct(int j, int k, int &jnext, int &knext, int &dir_next, int nsw, edge_crossings &xedges, edge_crossings &yedges)

Find next point starting from a point on a bottom edge.

void find_intersections(size_t ilev, double &level, edge_crossings &xedges, edge_crossings &yedges)

Find all of the intersections of the edges with the contour level.

void edges_in_y_direct(double level, interp_vec<ubvector> &si, edge_crossings &yedges)

Interpolate all right edge crossings.

void edges_in_x_direct(double level, interp_vec<ubvector> &si, edge_crossings &xedges)

Interpolate all bottom edge crossings.

void process_line(int j, int k, int dir, std::vector<double> &x, std::vector<double> &y, bool first, edge_crossings &xedges, edge_crossings &yedges)

Create a contour line from a starting edge.

void check_data()

Check to ensure the x- and y-arrays are monotonic.

contour(const contour&)
contour &operator=(const contour&)

Basic usage

template<class vec_t, class mat_t>
inline 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.

The types vec_t and mat_t can be any types which have operator[] and operator[][] for array and matrix indexing.

Note that this method copies all of the user-specified data to local storage so that changes in the data after calling this function will not be reflected in the contours that are generated.

template<class mat_t>
inline void set_data(const uniform_grid<double> &ugx, const uniform_grid<double> &ugy, const mat_t &udata)

Set the data.

The type mat_t can be any type which has operator[][] for matrix indexing.

Note that this method copies all of the user-specified data to local storage so that changes in the data after calling this function will not be reflected in the contours that are generated.

template<class vec_t>
inline void set_levels(size_t nlevels, vec_t &ulevels)

Set the contour levels.

This is separate from the function calc_contours() so that the user can compute the contours for different data sets using the same levels.

void calc_contours(std::vector<contour_line> &clines)

Calculate the contours.

Note

There may be zero or more than one contour line for each level, so the size of clines is not necessarily equal to the number of levels specified in set_levels().

Regrid function

void regrid_data(size_t xfact, size_t yfact, size_t interp_type = o2scl::itp_cspline)

Regrid the data.

Use interpolation to refine the data set. This can be called before calc_contours() in order to attempt make the contour levels smoother by providing a smaller grid size. If the original number of data points is \( (\mathrm{nx},\mathrm{ny}) \), then the new number of data points is

\[ (\mathrm{xfact}~(\mathrm{nx}-1)+1, \mathrm{yfact}~(\mathrm{ny}-1)+1) \]
The parameters xfact and yfact must both be larger than zero and they cannot both be 1.

Public Types

typedef boost::numeric::ublas::vector<double> ubvector
typedef boost::numeric::ublas::vector<int> ubvector_int
typedef boost::numeric::ublas::matrix<double> ubmatrix

Public Functions

contour()
~contour()