Help for MARSUVW

PURPOSE:

MARSUVW computes a surface normal for each pixel, given an XYZ image as
input.  This surface normal is a unit vector pointing "out" of the surface
(up, for an in-situ image of the ground).  UVW is a term used to name the
surface normal unit vector components, to avoid confusion with the XYZ
coordinates in Cartesian space of the input.  The U component aligns with
the X axis, V with Y, and W with Z.

The input for MARSUVW is an XYZ image such as that created by the MARSXYZ
program.  It can be a single 3-band file or three 1-band files; see MARSXYZ
for details.  The C point of the camera model is also extracted from the input
image, which is used only to determine which way is "out" of the surface.

The output of MARSUVW is either a single 3-band file, or three 1-band files.
Each band is of type REAL, and the U, V, and W components are in order,
corresponding to the X, Y, and Z axes respectively.  As with XYZ images,
the output image may contain holes (missing data), defined by MISSING_CONSTANT.
There will likely be many more holes than in the input image, since the
normal algorithm requires a patch around each pixel.

EXECUTION:

marsuvw inp=data.xyz out=data.uvw
where:
data.xyz is an input 3-band image of type REAL with the X, Y and Z values
at that pixel in meters (the unit is actually irrelevant).
data.uvw is an output 3-band image of type REAL the U, V, and W components
of the unit vector for each pixel.

Any program producing 3D position could be used; marsxyz is simply an example.

METHOD:

MARSUVW computes, for each pixel, the surface normal for the pixel by
analyzing a patch of neighboring pixels.  First, pixels with X/Y values
outside the box specified by X/Y_CENTER and BOX_RADIUS are rejected and
no normal is computed.  For the survivors, the neighboring pixels are
gathered, in a window specified by the RADIUS parameter.  The neighbor
points are then culled by distance in Cartesian space from the pixel in
question, specified by SEPARATION.

If enough points remain (MIN_POINTS), then a best-fit plane is computed.
The plane fit error (defined as the average distance of each point from
the plane) is compared to the ERROR parameter.  If the error is too big,
pixels with individual errors exceeding the REJECT ratio (times the average
error) are removed, the plane is recomputed, and the error checked again.
The process continues until the error is acceptable, or until fewer than
MIN_POINTS remains.

After this process, if enough points still remain, then the normal to the
computed plane is written to the output.

The plane fit is computed by forming a covariance matrix for the locations
of the candidate points relative to the centroid.  The eigenvectors and
eigenvalues of the matrix are computed.  The eigenvector corresponding
to the smallest eigenvalue indicates the direction in which there is the
least spread of points - i.e. the plane normal.  Some checks on the relative
magnitudes of the eigenvalues are used, since even a fuzzy blob of points
is likely to have one eigenvalue smaller than the others.

In order to determine whether the computed normal or its negation are "out"
of the surface, the normal is dotted with the vector from the camera to the
point.  If the angle is less than 90 degrees (dot product > 0), the vector
is used, otherwise -vector is used.  This is theoretically all that is needed,
but in practice noise in the XYZ values occasionally creates normals that
tip slightly away from the camera, when looking nearly edge-on to a distant
plane (such as near the horizon).  In order to mitigate this effect, a
threshold is applied to the W (Z) component of the vector.  If the component
is larger than the threshold, then the negation of the vector is returned.
This reduces the flipping problem at the expense of not being able to properly
detect nearly horizontal overhangs.  Since such overhangs are exceedingly
rare, the threshold test is usually a good thing.  Note that a value > 1.0
will effectively disable the sign flipping.  The flipping is only performed
if the point is farther away from the camera than FLIP_DIST.

The program can also be run in, for lack of a better name, "slope" mode
(using -slope).  This mode was originally developed to compute surface
normals over a rover-sized patch for determining terrain slope (as opposed
to an arm-sized patch in the original code).  There are three differences in
slope mode:

1) The window size (RADIUS) is dynamically adjusted based on the range to
the pixel (range is usually in meters).  Less than 10 is full size, 10-15
is 75%, 15-20 is 66%, and >20 is 50%.  This helps with efficiency over large
areas.

2) Not every pixel is considered.  Because the SEPARATION is much larger for
the slope case, every 2nd, 3rd etc. pixel may be skipped, again for efficiency.
The distance steps are the same as above: >20 uses every pixel, 15-20 uses
every other pixel, 10-15 uses every third pixel, and <10 uses every 4th pixel.

3) There is no REJECT loop.  The normal is computed once, and if the plane
fit error is greater than ERROR, the point is rejected and no normal is
computed.

It is recommended for slope mode that larger values be used for RADIUS and
SEPARATION, and that a very large value be given for BOX_RADIUS to effectively
disable it.

The actual algorithm and core code for marsuvw was developed by Chris Leger
of Section 348.

Because a certain number of neighboring pixels are needed to compute the
normal, the output will typically be much sparser than the XYZ input image.
Reducing the number of neighbors required (MIN_POINTS) can help with the
sparseness of the image, but the resultant normals are likely to be
less accurate due to noise.

The output will be expressed in the coordinate system specified by the
COORD and COORD_INDEX parameters.


Parallel Processing
-------------------
This program has been parallelized using Open MP (OMP), which is built in
to the g++ compiler.

By default the number of threads used equals the number of cores on the machine
where the program is being run.  Each image line is assigned to a different
core, with "dynamic" scheduling to keep the workload for eeach core similar.

Parallel processing can be disabled via the -OMP_OFF keyword.  The number
of threads can be controlled by setting the OMP_NUM_THREADS environment
variable before running the program.  There are numerous other OMP variables
that can be set; see the OMP documentation.  However, the number of threads
is the only one that is likely to be useful in most cases.


HISTORY:
  2003-04-15 B. Deen - Initial marsuvw based on code from Chris Leger (348).
  2005-02-17 B. Deen - Merge of marsuvw_slope, and addition of "flip" capability.
  2015-04-16 B. Deen - Added CAMERA_CENTER parameter.
  2016-06-09 B. Deen - Parallelization of code.
  2020-05-21 W. Bunch - Replaced sprintf calls; added unit test.

COGNIZANT PROGRAMMER: B. Deen

PARAMETERS: