MGsystem Class Reference

MAIN INTEFFACE: Multigrid schemes for solving least squares approximation to scattered data with B-splines More...

#include <MGsystem.h>

List of all members.

Public Types

typedef boost::shared_ptr<
std::map< int, double > > 		shared_map

Public Member Functions

 MGsystem ()
void initMultilevel (int m0, int n0, int h, boost::shared_ptr< std::vector< double > > U, boost::shared_ptr< std::vector< double > > V, boost::shared_ptr< std::vector< double > > Z, boost::shared_ptr< GenMatrix< UCBspl_real > > x, shared_map weights=shared_map(new std::map< int, double >))
void solveFMG ()
void solveAscend ()
void solveVcycle ()
void relax (int numberOfIterations, int relaxType=3)
UCBspl::SplineSurface getSplineSurface () const
double residual_l2 (bool scaled=false) const
double residual_linf () const
void setDomain (double umin, double vmin, double umax, double vmax)
void getDomain (double &umin, double &vmin, double &umax, double &vmax)
void setSmoothingFactor (double smoothingFac)
bool addPoint (double u, double v, double z, double weight=1.0)
void setWeight (int pointIndex, double weight)
void setNumberOfIterations (int noIter)
void setNumberOfIterationsCoarsest (int noIter)
void cleanup ()
void qweInitX (double value)
void qweSetRHS (const GenMatrix< UCBspl_real > &t_k)

Detailed Description

MAIN INTEFFACE: Multigrid schemes for solving least squares approximation to scattered data with B-splines

MGsystem - The class takes a set of scattered data as input and produces a B-spline surface that approximates the scattered data. The resulting B-spline surface will be smooth and/or accurate depending on the user input. One of the strengths with the methods implemented in this class is, if desired, smooth and "natural" extrapolation of the surface outside the domain of the scattered data or to areas of the domain with no scattered data.

The approximation scheme is a least square fit with a thin plate spline smoothing term. A set of multigrid schemes are implemented to solve the equation system that arises. Different choices of relaxation is implemented in the class LSsystem. The functions solveFMG, solveAscend and solveVcycle (and relax) produces the spline surface from the scattered data. This is done by solving a large linear equation system, but due the multigrid schemes this is extremely fast compared to non-multigrid schemes. Though, due to the fact that iterative equation solvers are used, the equation system is not solved exact. But, there are mechanisms for adjusting the initial solution; see below. The spline surface can be retrieved as a SplineSurface object with getSplineSurface.

No stopping criteria for the equation solvers are provided - the user control is only through the given number of iterations.

Note that the format of the scattered data and the solution vector is exactly the same as in the SINTEF MBA (Multilevel B-spline Approximation) Library. Thus, if one uses the V-cycle scheme or simple relaxation only, which starts at the finest level, a good starting vector can be found by running MBA first. This is also necessary for surfaces with large coefficient grids, which otherwise will converge very slowly with the V-cycle scheme or simple relaxation only.

Examle of use (minimal with default parameters): 

  MGsystem mgsys();
  m0=n0=1; h=7; // Spline surface with 2^h + 3 = 131 coefficients in each direction
  mgsys.initMultilevel(m0,n0,h,xArr,yArr,zArr,PHI);
  mgsys.solveFMG(); // Solve equation system by running the Full Multigrid Scheme

Some hints on usage and setting parameters different from the default.

Choice of solver: The most common multigrid solver is solveFMG, but solveAscend is equally good in many cases and faster when running with default parameters. These two solvers should be used if there is no good starting vector present (see above).

Smoothing factor and number of iterations: The influence of the smoothing factor on the produced B-spline surface is extensively explained in setSmoothingFactor. One should note that when increasing the smoothing factor, more iterations are needed to obtain a correct surface. This can be controlled by the user in different ways: One may just set more iterations at each level with setNumberOfIterations; e.g., 50 or 100 iterations instead of 10 which is the default. But, in most cases this will work better: First run the solver (for example solveAscend) with the default smoothing factor (1.0). Then call setSmoothingFactor and run relax. The last command will run the given number of iterations on the final surface with the new smoothing factor. If you are not satisfied with the result (see Controlling the result below), you may repeat the last two steps with new parameters. Alternatively one may replace relax with solveVcycle.

Controlling the result and visualizing the surface: One way of controlling the result is through residual_l2 and residual_linf. One should note though that even if one or both of the residuals are small, this does not guarantee that the solution of the equation system is accurate (which is due to an ill-conditioned equation system).
The safest way to control the quality of the produced B-spline surface is through visual inspection. The surface can be retrieved at any time with getSplineSurface. This gives a SplineSurface object that can be sampled ( z=f(x,y) ) or piped further to utility functions that prints the surface in visualization format, e.g., VRML.

Extrapolation: The mathematical formulation of the B-spline surface used here (least squares approximation of scattered data with smoothing term) gives overall good extrapolation properties. That is, the surface can be extrapolated beyond the xy-range of the input data in a "natural" way by using the function setDomain before creating the surface. Such "natural" extrapolation of surfaces is important in many applications, for example in geological modelling where horizons must be extrapolated to intersect faults.

Author:
Øyvind Hjelle <Oyvind.Hjelle@math.sintef.no>
See also:
MBA

Member Typedef Documentation

typedef boost::shared_ptr<std::map<int, double> > MGsystem::shared_map

Initialize.

Parameters:
m0,n0 (>=1): The initial size of the spline space in the hierarchical construction. If the rectangular domain is a square, m0=n0=1 is recommended. If the rectangular domain in the y-direction is twice of that the x-direction, m0=1, n0=2 is recommended. In general, if the rectangular domain in the y-direction is k times the length in the x-direction, m0=1, n0=k is recommended.
h,: Number of levels in the hierarchical construction. If, e.g., m0=n0=1 and h=8, The resulting spline surface has a coefficient grid of size 2^h + 3 = 259 in each direction.
U,: Array with x-values of scattered data
V,: Array with y-values of scattered data
Z,: Array with z-values of scattered data
x,: The solution vector, i.e., the tensor product coefficients (formatted as a matrix).
If the solution vector is allocated to the correct size and solveVcycle or relax is used, the solution vector must also be initialized. (For example to the average of the z-values, or to zero).
weights,: (Not required). Each point can be given a weight of importance. The default weight for each point is 1.0.
Example: Assume that points with index 17 and 80 (starts at zero) should have weights 5.0 and 2.5 respectively. The weight argument is then initialized thus: 
  boost::shared_ptr<std::map<int,double> > weights(new std::map<int,double>);
  (*weights.get())[17] = 5.0;
  (*weights.get())[80] = 2.5;

Constructor & Destructor Documentation

MGsystem::MGsystem (  )  [inline]

Default parameters should be used through the API.

Member Function Documentation

void MGsystem::solveFMG (  ) 

Full MultiGrid scheme (starting from the finest level).

void MGsystem::solveAscend (  ) 

Simple ascend method from the coarsest to the finest grid (nested iteration)

void MGsystem::solveVcycle (  ) 

One full V-cycle (starting from the finest level). Run MBA first to find a good starting vector.

void MGsystem::relax ( int  numberOfIterations,
int  relaxType = 3 
)

Run the given number iterations at the finest level. This can be done, e.g., to improve the solution from one of the other solvers, or after a new smoothing parameter has been set with setSmoothingFactor. See also the detailed description of the class in the header.

Parameters:
relaxType 1=Gauss Seidel
2=Jacobi
3=Conjugate Gradients (CG) (recommended)
4=Preconditioned Conjugate gradients

UCBspl::SplineSurface MGsystem::getSplineSurface (  )  const

Retrieve the spline surface

double MGsystem::residual_l2 ( bool  scaled = false  )  const

The l_2 norm of the residual ||b - Ax||_2 after solving the equation system, possibly scaled with grid cell size if indicated. A low value indicates a more exact solution of the equation system than a large value.

double MGsystem::residual_linf (  )  const

The l_inf norm of the residual ||b - Ax||_inf after solving the equation system. A low value indicates a more exact solution of the equation system than a large value.

void MGsystem::setDomain ( double  umin,
double  vmin,
double  umax,
double  vmax 
)

Set the domain over which the surface is to be defined. The default is the xy-range of the scattered data. If used, this must be done before creating the surface.

Note:
This function can only be used to expand the domain beyond the xy-range of the scattered data, i.e., for extrapolation (see above). It is the users responsibility to check that no scattered data falls outside the domain. (use std::min_element and std::max_element to find range of data for std::vector)

void MGsystem::getDomain ( double &  umin,
double &  vmin,
double &  umax,
double &  vmax 
)

Get the domain over which the surface is defined

void MGsystem::setSmoothingFactor ( double  smoothingFac  ) 

Control smoothness of the resulting spline surface. A positive value should be given; default is 1.0.

This is a useful function for controlling smoothness versus accuracy of the resulting spline surface. More specifically, the given smoothing factor indicates how strong the thin plate spline energy should influence on the result. A high value gives a surface with "low energy" which is smooth, but it does not necessarily approximate the scattered data well. A low value of the smoothing factor gives a surface which is a good approximation to the scattered data in a least square sense. Thus the accuracy will dominate over smoothness such that the surface may be rough and not pleasant looking since it is forced to approximate all the scattered data well (even noise and outliers are approximated well).

Theoretically, when the smoothing factor increases, the surface converges to a least squares plane through the scattered data.

One should experiment with this function on typical data. For terrain modelling one may try a smoothing factor in the range 0.2 - 1000.0, but lower/higher values may also be useful.

Note:
A value of zero or close to zero may cause break down of the linear equation system in many cases when the scattered data are not uniformly distributed in the domain, and/or the tensor product grid is large compared to the number of input data points.
See also detailed description in the header.

bool MGsystem::addPoint ( double  u,
double  v,
double  z,
double  weight = 1.0 
)

Add a point to the data set. This can be done after a surface has been calculated from the initial data set as in the example above, and if one want to add more scattered data points and update the B-spline surface accordingly. One of the solve...-functions or relax() must be run afterwards to take the new point into account. More than one point can be added before solve... and/or relax() is run again. The given point is also added to the dataset. A weight of importance different from 1.0 can be given.

Return values:
bool false if (u,v) is outside the domain.

void MGsystem::setWeight ( int  pointIndex,
double  weight 
)

Set weight of importance of the point with the given index (starting from 0). Default weight is 1.0

void MGsystem::setNumberOfIterations ( int  noIter  )  [inline]

Set number of iterations at each level other than the coarsest. Default is 10

void MGsystem::setNumberOfIterationsCoarsest ( int  noIter  )  [inline]

Set number of iterations at coarsest level. Default is 50

void MGsystem::cleanup (  ) 

Clean-up array structures and reduce memory usage

The documentation for this class was generated from the following files:
Generated on Wed Nov 28 12:27:29 2007 for LSMG by  doxygen 1.5.1