class Map2D: public TH2D

 Overloads TH2D::Add, has the same basic behaviour,
 except that zero plus anything is zero.

 Overloads TH2D::Add, has the same basic behaviour,
 except that zero plus anything is zero.

 Translates the map image by a distance (dx,dy) in map coordinates,
 and shifts all pixels up or down by vertical distance dz in
 height units.  Image regions shifted outside the range of the
 map are lost.  Image regions shift into the range of the map
 from outside are zero.  Original contents are overwritten.

 Same as Shift(), except a new object is created to hold
 the results and returned, and the original object is not
 touched.  The caller must delete the returned object.

 Rotates the map image by Euler angles phi,theta,psi about the center
 of the map surface.  Angles may be in degrees (degrees=1) or
 radians (degrees=0).  The angle phi determines the direction of
 the theta rotation axis in the xy plane, and does not generate
 an independent rotation.  Map rotations only make sense for small
 theta.  The theta rotation is implemented as a rotation about the
 x axis by angle theta*cos(phi) followed by a rotation about the y
 axis by angle theta*sin(phi), with errors entering at the level of
 theta^3.  Rotations in the xy plane (angle psi) may be of any size,
 and may result in clipping of the image.  Regions of the map that
 are rotated outside the map range are lost.  Original contents of
 the map are overwritten.

 Rotates the map image about the x axis by angle thetarad in radians.
 The rotation axis is located at coordinate y0 and height z0.
 The map must remain single-valued after the rotation, which places
 an upper limit on the magnitude of theta that depends on the
 height profile of the map relative to z0.  An error message is
 printed if this constraint is violated.  A new Map2D object is
 created to hold the result, and the original object is untouched.
 The caller must delete the returned object.

 Rotates the map image about the y axis by angle thetarad in radians.
 The rotation axis is located at coordinate x0 and height z0.
 The map must remain single-valued after the rotation, which places
 an upper limit on the magnitude of theta that depends on the
 height profile of the map relative to z0.  An error message is
 printed if this constraint is violated.  A new Map2D object is
 created to hold the result, and the original object is untouched.
 The caller must delete the returned object.

 Rotates the map image about the z axis by angle thetarad in radians.
 The rotation axis is located at (x0,y0) in map coordinates.  There
 are no assumed bounds on phirad.  Regions of the image that move
 outside the map range are clipped and lost.  Regions of the image
 that move into the range from outside are zero.  A new Map2D object
 is created to hold the result, and the original object is untouched.
 The caller must delete the returned object.

 Rescales the map image along all three axes by independent scale
 factors.  The origin of the rescaling is taken to be the center
 of the map in x and y, and z=0.  The original map is overwritten.
 Note: the special case of sx=-1,sy=1,sz=-1 or sx=1,sy=-1,sz=-1
 can be used to invert the surface, effectively flipping it over.

 Creates a new map with bounds set by corner coordinates (xlow,ylow)
 and (xhigh,yhigh), then copies this map into the new one.  The
 original map (*this) is not modified.  The name of the new map is
 the name of this map with the suffix "_rN", where N is an integer
 starting with 1 and counting forward.

 Rescales the map image along all three axes by independent scale
 factors.  The origin of the rescaling is (x0,y0) in map coordinates,
 and z0 in height units.  A new object is created to hold the result
 and the original object is untouched.  The caller must delete the
 returned object.

 Computes the height correlation with another map, useful for
 aligning two images of the same terrain.  The contrast
 argument serves as a weight factor on terms in the sum
 where one or the other of the two maps is zero.  This allows
 control over how zero regions are treated, such that
	contrast=0 : ignore zero regions
	contrast=1 : treat them like plateaus with z=0
   contrast=BIG : boundaries dominate, maps are B/W
 Agreement between their coordinate ranges is not required,
 but the two maps must have the same dimensions.

 Copies the map from the source Map2D object into this object,
 overwriting this object's data but preserving its dimensions and
 axis limits.  Matte pixels in the source map are ignored, making
 this method an effective way to overlay two non-overlapping images.
 No new objects are created.

 Overlays the input map on the output, and masks (overwrites with matte)
 any cells that are zero in the input map.  Cells in the input map
 that are non-zero are overwritten in *this with the constant "value"
 unless value=0, in which case those pixels are preserved.  If argument
 "neg" is non-zero, the input map is treated as a negative.  Both input
 and output maps must have the same dimensions.  No new objects are
 created.

 Scans through the map and finds any pixels less than the supplied
 threshold value, replacing them with the value zero.

 Searches the map for regions.  Pixels belong to the same region if
 an unbroken path of contiguous pixels can be found that joins them.
 Two pixels are contiguous if they share a common edge; touching at the
 corners does not count.  Regions up to maxsize pixels are erased
 (overwritten with matte) and a count of the remaining islands of greater
 than maxsize pixels is returned.  If called with maxsize=0 then the
 total number of islands is returned, and the map is not modified.

 Searches the map for localized spikes that appear in regions of
 smoothly varying height, and replaces any spikes that are greater
 than maxheight above or below the local average surface height with
 the local neighborhood average.  The size of the neighborhood is
 (2*regionsize+1)x(2*regionsize+1) in pixels.  This operation is not
 reversible.

 Finds the center of gravity of the non-zero portion of the map
 and applies a shift in x,y to locate the center at the midpoint.

 Locates the left-most and the right-most pixels in the map with
 non-zero values and applies an in-plane rotation such that those
 two pixels are in the same row (y-value) of the map.

 Finds the maximum and minimum values of the map, excluding cells
 with zero value, and sets the Max and Min to the corresponding values.

 Forms a closed path around the outermost periphery of the map, and
 computes the total curl of a 2D vector field within that path using
 Stoke's theorem.  The two components of the vector field within the
 plane of the map are accepted as arguments Vx,Vy.  The result is
 returned in the units of V times coordinate length.  The two maps
 given as arguments must have the same dimensions and bounds in x and
 y.  The object being used to compute the curl is overwritten with a
 map filled with zeros everywhere except on the contour around which
 the Stokes integral was performed, with each cell on the contour
 containing the running integral value up to that point.  If optional
 argument silhouette is provided by the caller, it is overwritten upon
 return with a pointer to a new Map2D object with the same shape as
 Vx,Vy whose contents represent an edge-enhanced silhouette of the
 overlaid input maps, with bin values defined as follows.
        0    : the bin is outside the contour,
      (0,1)  : the bin is inside the contour,
     [1,inf) : the bin is on the contour.
 It is up to the user to delete the silhouette object when it is no
 longer needed.

 Computes the curl of the 2D vector field represented by the components
 Vx,Vy and stores the result in the *this map.  The value stored in
 each cell of the output map represents the integral of the curl within
 that cell, so that the sum over all cells in the map should equal the
 line integral of the field around a closed loop enclosing the non-zero
 field region, according to Stoke's theorem.  The complementary method
 TotalCurl() is provided to compute the same result using the loop line
 integral.  The return value is the surface integral of the curl.

 Computes a correction to a 2D vector field defined on a plane such
 that the resulting field has zero curl.  This essentially replaces
 one of the two degrees of freedom of the field with the constraint.
 The argument "comp" tells which of the two components is replaced:
      comp = 0: return a replacement for Vx;
      comp = 1; return a replacement for Vy.
 The read-only map arguments Vx,Vy are not modified.

 Algorithm:
 There are an infinite number of replacements for Vx,Vy such that
 Del x (Vx,Vy) = 0, all of which are related to one another through
 gauge transformations.  The purpose here is to modify the field in
 as minimal a fashion as possible to eliminate the curl.  To simplify
 the problem, consider only the subset of solutions that leaves Vx [Vy]
 unchanged, as selected by the user.  This reduces the solution to
 finding a set of independent columns [rows] that zeros the curl within
 that column [row] subject to some criterion for minimal modification.
 In 1D the zero-curl constrain reduces to a first-order differences
 equation (first-order inhomogenous differential equation in 1D written
 in finite-element form) which has just one overall undetermined offset.
  I chose the offset to minimize the sum of the squared differences
 from the corresponding input map.  This uniquely determines the result.

 Computes the divergence of the vector field V whose x,y components
 are passed in as the maps Vx,Vy.  The contents of map *this are
 overwritten with the result.  The return value is the total
 divergence integrated over the output map, in units of V times
 coordinate length.  The two maps given as arguments must have the
 same dimensions and bounds in x and y.

 Computes the gradient of the scalar field V in map *this and
 stores the x,y components of the gradient in output maps Vx,Vy.
 Arguments Vx,Vy must point to maps which have been created by the
 user at input, otherwise the corresponding gradient component is not
 computed.  At output these maps are overwritten by the computed
 gradient components.  If both pointer arguments are zero at input
 the routine does nothing and returns 0, otherwise it returns the
 total number of pixels in the map *this for which the gradient vector
 has been computed.

 Solves the Poisson equation in two dimensions with a source given by
 input map rho, subject to the boundary conditions contained in this
 object at entry.  If input argument "mask" points to an existing object,
 that object is interpreted as a mask for the pixels in map *this that
 contain boundary values, otherwise all non-zero pixels in map *this
 at input are considered to specify boundary values.  The algorithm
 will make use of as many (or few) boundary values as the user provides,
 although the solution may not satisfy all of them exactly if the
 problem is over-determined, i.e. too many values are provided, or they
 are clustered in too compact a region of the map.  In any case, the
 solution will attempt to match all of the provided data points as
 closely as possible, subject to the constraints of Poisson's equation
 and the resolution of the map.  All cells in the map *this are over
 written with the solution except those specifying boundary conditions.

 Return value:
 The return value is the rank of the matrix that expresses LaPlace's
 equation for the given boundary conditions.  The rank expresses the
 effective number of independent boundary values that were used in
 attempting the satisfy the stated boundary conditions.

 Algorithm:
 First the Green's function for the 2D Poisson equation with a point
 (single pixel) source is evaluated for each non-zero pixel in rho,
 and the results are superimposed.  At this point, the solution obeys
 Poisson's equation but with arbitrary boundary conditions. The map
 is now subtracted from map *this that specify the boundary values
 for the problem, and the resulting map is handed to LaplaceSolve()
 to find the homogeneous part of the solution.  Adding the homogenous
 and inhomogenous parts yields the full solution, which is returned
 in this object.

 Timing:
 The time required to complete this function increases quickly with
 the dimensions of the map.  The following figures are for a single
 Intel Core2 quad-core processor at 3.4GHz running a single thread.

   t = (56 + 5.3e-8 Ndim^4.2 + 4.9e-9 Ndim^4.2) seconds
         ^    ^               ^
         |    |               |
         |    |                \.
         |    |                 \--- compute the homogeneous part
         |    |
         |     \--- compute the Greens function integral
         |
          \--- compute the Greens function (35x35 discrete kernel)

 Example: for a 1000 x 1000 pixel map, t = 229,000 seconds (64 hr)

 Iterates the Poisson equation in two dimensions with a source given by
 input map rho, subject to the boundary conditions contained in this
 object at entry.  If input argument "mask" points to an existing object,
 that object is interpreted as a mask for the pixels in map *this that
 contain boundary values, otherwise all non-zero pixels in map *this
 at input are considered to specify boundary values.  The algorithm
 will make use of as many (or few) boundary values as the user provides,
 although the solution may not satisfy all of them exactly if the
 problem is over-determined, i.e. too many values are provided, or they
 are clustered in too compact a region of the map.  In any case, the
 solution will attempt to match all of the provided data points as
 closely as possible, subject to the constraints of Poisson's equation
 and the resolution of the map.  All cells in the map *this are over
 written with the solution except those specifying boundary conditions.

 Return value:
 The largest correction that was applied to any cell in the map in the
 present iteration, or zero if no changes were made.

 Algorithm:
 The following formula is applied to all cells except those with rho=0
 and the ones that specify boundary conditions.  Boundary condition
 cells are identified by having a non-zero value in the corresponding
 cell of the mask map, or if mask=0, by having non-zero values in the
 *this map at entry.

           / (n)    (n)   \    2   / (n)    (n)   \   2            2  2
          | u   +  u       | dy + | u   +  u       | dx  - rho   dx dy
  (n+1)    \ i-1,j  i+1,j /        \ i,j-1  i,j+1 /           i,j
 u     =   -------------------------------------------------------------
  i,j                        /  2      2 \                            .
                          2 | dx  +  dy   |
                             \           /                            .

 Timing:

  The Green's function for the 2D Poisson equation can be written most
  compactly as follows,

                             1        |      2        2 |
             G(x,y;x ,y ) = ----  log |(x-x ) + (y-y )  |
                    0  0    4 pi      |    0        0   |

 where x0,y0 is the location of the source.  In two dimensions it is
 not possible to set the boundary conditions for G->0 at x0,y0->infinity
 so instead one choses that G=0 at some finite radius R.  The expression
 above has been written with R=1, but adding a constant term can be used
 to shift R to an arbitrary radius, except for the singular values 0 and
 infinity.  The result is loaded into this map, overwriting any contents
 that were present at input, but leaving the dimensions and bounds of
 the original map unchanged.  The point x0,y0 is assumed to lie somewhere
 within the bounds of the map.

 This is the correct expression for the continuum case, but it diverges
 as x,y -> x0,y0, and has large discretization errors for pixels nearby.
 For this reason, I take a NxN pixel square region centered around x0,y0
 and replace the computed value of G for these pixels with the solution
 to the discrete NxN matrix Poisson equation on a square grid, subject
 to the boundary condition that the cells around the outer perimeter of
 the NxN pixel block match onto the above continuum formula.  By chosing
 N to be sufficiently large, the residual discretization errors are kept
 to a negligible level.

 Solves the Laplace equation in two dimensions, subject to the boundary
 conditions contained in this object at entry.  If input argument "mask"
 points to an existing object, that object is interpreted as a mask for
 the pixels in map *this that contain boundary values, otherwise all
 non-zero pixels in map *this at input are considered to specify
 boundary values.  The algorithm will make use of as many (or few)
 boundary values as the user provides, although the solution may not
 satisfy all of them exactly if the problem is over-determined, i.e.
 too many values are provided, or they are clustered in too compact a
 region of the map.  In any case, the solution will attempt to match all
 of the provided data points as closely as possible, subject to the
 constraints of the Laplace equation and the resolution of the map.
  All cells in the map *this are over written with the solution.

 Algorithm:
 The problem is solved in the basis of solutions to Laplace's equation
 in polar coordinates.  Here f_n is a stackable column vector of two
 basis functions

                             /                              .
                         m  |  cos(m phi)
            f (x,y) =   r  <
             m              |  sin(m phi)
                             \                              .

 There is a second set of solutions that are singular at the origin,
 but these are excluded because the interior region of the map is
 assumed to include the origin, taken to lie at the geometric center
 of the map.  Note that for n=0 the second element of f_n is identically
 zero.  To form a column vector of basis functions F_i, imagine a stack
 of f_m column vectors starting with m=0 at the top, and increasing
 downward until the number N of elements in F equals at least M, the
 number the pixels (x_j,y_j), j=0..M that are provided as boundary
 values in the input map.  This results in a large MxN matrix A_ji
 such that
                   A     = F ( x , y )
                    j,i     i   j   j

 Solving Laplace's equation now amounts to solving the linear equation

                   A    C  = z
                    j,i  i    j

 for unknown coefficients C_i in terms of initial values z_j for pixel j.
 In general A is not a square matrix, and has many singular values.  I
 solve this equation using singular value decomposition, and identify
 the rank of the matrix A, which is passed as the return value from this
 method.  The results of the singular value decomposition are used to
 filter the given values z_j into the unknown coefficients C_i, which
 are then used to fill the map with the optimal solution to the boundary
 value problem.

 The user should pay careful attention to this value, as it may
 indicate that the problem is ill-posed.  In general, one expects that
 the rank will be close to the number of boundary values provided.  It
 will never be greater.  If the rank is significantly less than the
 number of boundary values provided, it may indicate that the boundary
 values over-determine the solution, or it may simply mean that fewer
 basis functions are needed to describe the solution than the number
 of boundary values would naively suggest.

 Scans over the map and forms a histogram of the cell heights,
 omitting the cells that contain zero.