Help for NSYTTILT

PURPOSE:

NSYTTILT computes the tilt an InSight instrument would have if placed at
the given point.  The minimum and maximum tilts are provided, as well as a
status band, and the program accounts for clock angle rotation and sinkage.
It will also compute the delta-tilt between the SEIS and WTS instruments.

The input for NSYTTILT is an XYZ image such as that created by MARSXYZ.
This can be a single 3-band file (preferred) or three 1-band files.

The output of NSYTTILT is a 3-band file of type REAL.  The first band is
a status (flag) band indicating whether or not the result is within limits.
The second and third bands are the minimum and maximum tilts, respectively.
NSYTTILT additionally outputs instrument-specific surface normal and Z value
files, which are needed by NSYTROUGH.

EXECUTION:

nsyttilt inp=data.xyz out=data.tfs uix_out=data.uis zix_out=data.zis -seis
nsyttilt inp=data.xyz out=data.tfw uix_out=data.uiw zix_out=data.ziw -wts
nsyttilt inp=data.xyz out=data.tfh uix_out=data.uih zix_out=data.zih -hp3
nsyttilt inp=data.xyz out=data.tds -seis -delta

METHOD:

At each pixel in the output file, NSYTTILT attempts to "place" the given
instrument at that point.  The point represents the grapple on the instrument,
so the XY position of the point is what would be communicated to the IDS
team to actually place the instrument there on the surface.  The XY position
is derived from the XYZ input file for the given pixel.  If there is no XYZ
position for the pixel, the pixel is skipped.

Note that the program works internally in Site frame.  This aligns all the tilts
with the gravity vector, giving true tilt relative to gravity. It will convert 
the given XYZ to this frame automatically.  While other frames could be used 
(e.g. -rover for Lander frame), the tilts would then be measured relative to the
lander, which is generally not what it desired.

Instrument Feet
---------------
The size and spacing of the feet for each instrument is given in the
nsyt_foot_dimensions.h file that is compiled into the program.  This should
be read from a config file, but time constraints may prevent this from
happening.

The program represents the feet as a set of circles - 3 each for SEIS and WTS,
4 for HP3.  The size of these circles matches the size of the feet as they
interact with the ground.  Although one should consult the include file for
specific values, at the time this help was written, the circle radii are 3cm
for SEIS, 4.2665cm for WTS, and 5cm for HP3.

The Z location of the foot is defined as the highest point within the foot
patch, the idea being that the foot will "perch" on the highest point available.
However, this is sucseptible to noise spikes in the data.  Therefore, a
statistical analysis (mean and standard deviation, or sigma) of the Z points
within the foot patch is done, and all points exceeding some number of sigma
away from the mean are excluded from the analysis.  The threshold is determined
by the FILTER_SCALE parameter, which is then multiplied by sigma.

Tether Effects
--------------
The feet are initially set assuming the tether points directly down the
-X axis.  However, the actual instruments while being deployed will tend to
rotate so the tether points back toward the spacecraft.  This is accounted
for by rotating the instrument (i.e. feet position) based on the tether pivot
points.

For SEIS, the tether emanates from a single tether box position.  This
position is hardcoded at x=0.145, y=-0.101 in the lander frame.  The tether
points radially from this point to the instrument deployment location, and
the SEIS feet are rotated to match.

For HP3, the algorithm is more complicated, because the tether can be routed
different ways based on where the instrument is being deployed.  The algorithm
comes from Won Kin in Section 347:

  HP3 has two "virtual" tether anchor points.
   - HP3 tether guard post1 --- P1 = (X1, Y1) = (0.083, -0.435)
   - HP3 mount --- P2 = (X2, Y2) = (-1.279, 0.149)
   - Line equation between two anchor points:
     y = ((Y2-Y1)/(X2-X1)) * (x - X1) + Y1 = mx + b
   - Given the current tether position (X0, Y0):
     - If Y0 > m * X0 + b, use HP3 tether guard post1 P1 as the "virtual" anchor
       point.
     - Otherwise, use HP3 mount P2 as the "virtual" anchor point.
 You can determine the tether direction by assuming that the tether passes
 through the "virtual" anchor point.

For WTS, there is no tether.  However, the arm moves slowly enough that the
WTS is likely to rotate with the arm, meaning that there is a most-likely
rotation value that acts as if there was a tether.

Additionally, there are radial and cross-radial offsets between the WTS and 
SEIS, because the grapples are not coaxial in the nominal placement.  The WTS 
position is thus offset by an amount determined in the WTS_OFF parameter. These
offsets are applied after WTS rotations (instrument clocking and rotation due to
tether) are done. All WTS rotations (such as for clocking) are done centered 
around the nominal(SEIS) location. This means that WTS clocking rotation 
(see below) also introduces a small translation in the WTS.

Instrument Clocking
-------------------
Although the tether controls the rotation of the instrument while hanging on
the arm to some extent, the instrument can still rotate within a certain range.
This rotation is called "clocking".  

Zero clocking is defined as the rotation resulting from the tether calculation.

The program analyzes the tilt throughout this range of clocking angles.
Because it is computationally intractable to compute every possible angle,
we settle for checking a discrete set of clock angles.  The set of angles
checked are controlled by the CLOCK_RANGE and CLOCK_STEP parameters.  Generally
these should be set to be centered around 0 clock angle (and to include 0
clock angle), but this is not mandatory.

Settling of Feet
----------------
The Z location of each foot is defined as the highest point found within the
foot circle (after removing outliers) - the theory being, the foot will be
perched on the highest point in the area.  However, if the ground is sandy,
the feet are likely to "settle" into the sand by some amount.  For this reason,
the feet are "settled" by an amount specified by the SINKAGE parameter.  This
amount is added to each foot location (because +Z is down, adding the value
lowers the foot location).  Tilt analysis is then done using all combinations
of settled and unsettled feet.  This should pick up the worst-case scenario on
both ends of the spectrum - no settling and maximum settling.

Tilt Computation
----------------
For SEIS and WTS, tilt is computed by taking the Z value of each foot,
constructing vectors between them, and using the cross-product to determine
the normal of the plane defined by those points.  As noted above, this is done
for each combination of settled and unsettled feet, and the min and max tilts
across all combinations are gathered.

HP3 has 4 feet.  The 4 foot points are unlikely to lie in a single plane;
it it possible for the HP3 to sit unstably on any of its 3 feet.  For that
reason, the above computation is performed for all combinations of 3 of the 4
feet.  Thus the best (min) and worst (max) cases are determined no matter which
feet are actually stably on the ground.  This means that each combination of
3 feet is also analyzed for each combination of settled and unsettled feet.
In reality, the compliance of the ground is likely to allow all 4 feet to sit
on the ground, but this algorithm ensures we look at the worst case.

The min and max tilts are computed as above for each clock value, and then
the minimum and maximum across all clocks are computed to determine the final
output value.

Delta-Tilt Computation
----------------------
The delta-tilt is a comparison between the tilt of the SEIS and the tilt
of the WTS that covers it.  It is computed by first computing tilt as above
for the SEIS.  At each SEIS clock point, the tilt of the WTS is computed.
This also follows the above tilt algorithm, except the parameters
DELTA_WTS_RANGE and DELTA_WTS_TILT are used to determine the clock range.

At each combination of SEIS and WTS tilts, the delta tilt is computed, as
the arccos of the dot product between the two normal vectors.  The min and
max delta tilt is then computed across all combinations of clocks for both
instruments.

It is important to note that this is not finding the min/max tilt of each
instrument and then computing the angle between them.  It is computing the
angle between each and every combination of clock angles (subject to the
step size and range parameters), and determining the min/max of those.

Output File Format and State Band
---------------------------------
The output file consists of 3 bands in floating-point format:

Band 1: State band (contains goodness flag)
Band 2: Minimum tilt or delta tilt
Band 3: Maximum tilt or delta tilt

The min/max tilt/delta-tilt are expressed in degrees.

The State band can contain the following values:

0 = no value
1 = n/a
2 = n/a
3 = exceeds tilt threshold
4 = n/a
5 = tilt within limits (good)

As with all the InSight placement products, 5 indicates everything is within
limits, with other values indicating what parameter is out of limits.

Output Instrument Surface Normal, Z, and XYZ
--------------------------------------------
There are four optional output files, UIX_OUT, ZIX_OUT, UPX_OUT, and XPX_OUT.  
UIX_OUT and ZIX_OUT files contain intermediate results that are intended to be 
supplied as inputs to the NSYTROUGH program. UPX_OUT and XPX_OUT files are 
similar to UIX_OUT and ZIX_OUT except that they are in general measured in 
Lander frame instead of Site frame. Also, XPX_OUT contains all XYZ coordinates
instead of just the Z band. 

The UIX_OUT file is a 3-band float unit-vector normal file in the same format
as produced by the MARSUVW program.  However, the normal is calculated in a
special way.  It is the normal not of the surface, but of the instrument
itself, when located at the given pixel.

Examination of the tilt algorithm shows that the instrument normal is not
unique - the combinations of clock angles, sinkage, and multiple feet lead
to a variety of possible surface normals.  However, for the purpose of the
NSYTROUGH program, we need one "most likely" normal.  Therefore we choose
the one closest to the center of the clock range, with no sinkage.  We do
still look at all combinations of feet, and take the average normal across
those combinations (applies only to HP3).  Note that no sinkage and full
sinkage result in the same normal, since all feet sink the same amount.

The ZIX_OUT file is a 1-band float file that contains the Z position of the
instrument if it was placed at that location.  The Z position is the coordinate
of the point directly below the grapple, in the plane of the feet.  Similar to
UIX_OUT, this takes the center of the clock range.  Note that it uses the
nominal feet, NOT adjusted for settling.  There is a parameter in NSYTROUGH to
include settling if that is desired.

Because the surface normal is an average, using it with the foot locations of
the HP3 creates more than one solution.  In all cases the foot position is
used to fix the Z value of a plane that is normal to the instrument normal,
and then the Z position of the center of the instrument (grapple point) is
computed.  The final Z value is the lowest (highest Z since +Z is down) of
these computed positions - again, modeling the worst case (since that means
more interaction with the terrain).

Parallel Processing
-------------------
This program requires a lot of computational resources - especially delta
tilt.  The number of combinations of positions that must be checked can get
quite large.  For this reason, the 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 each 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.

Implementation Notes
--------------------
The instrument shapes should be in a data file but for now are compiled in
to the program, via nsyt_instruments(_h).com.

Clocking is simulated by stepping to different clock values and checking each.
There is a tradeoff between fidelity and execution time.  The smaller the clock
steps, the more accurate the result but the slower it runs (dramatically so).

Once an instrument foot point is determined, we must find the pixel that is
closest to that point as a starting point.  This is implemented by using a
K-D tree as implemented in the "nanoflann" package.  The POINT_EPSILON parameter
determines how far away this point can be from the ideal point (default .004m),
and the LEAF_MAX_SIZE determines how big the leaves of the tree can be before
splitting them.

From there, a region growing algorithm is used to determine all the XYZ points
that belong to the foot.  The region spirals out one pixel at a time until a
complete box (one pixel wide) is found with no hits within the confines of the
foot circle.  This region growing is bounded by the FOOT_WINDOW parameter; once
the box half-width exceeds the FOOT_WINDOW, the growing stops.  This is for
efficiency (to prevent spiraling out to the entire image in pathological
cases), but in practice the window should be set so the entire foot is within
it.

After region growing terminates, the mean and standard deviation of the
points in the region are computed.  Outliers (more than FILTER_SCALE * sigma
away from the mean in Z) are rejected.  Then the highest point (smallest Z)
amongst the remaining points is retained as the position of the foot.  The min
and max tilt are then computed, using a combinations library (combinations.hpp)
to pick up all combinations of feet.  All sets of settled/unsettled are
determined using a binary bit mask counting from 0 to 7 (one bit per foot).

Note that for efficiency, the foot Z locations are cached for reuse as they
are computed, since settling and clocking mean each foot location is used
many, many times.  The K-D tree must be searched for each foot location,
but once it is determined, the Z of the foot at that spot is computed only
once.  There are two separate caches for delta-tilt, since the SEIS and WTS
feet are not the same size.

NOTE:  Holes or gaps in the XYZ image are generally ignored.  The central
position of each foot has to match the given XYZ coordinate (within
POINT_EPSILON), so holes in the middle of a foot, if big enough, could cause
the location to be rejected.  But no analysis is done of other holes in
the area the foot "should" cover.  Holes represent unknowns, so normally the
presence of holes would be a red flag for placement.  However, the feet are
generally small enough that significant unknown terrain within one foot patch
is unlikely.  Also, we look for the maximum within the foot patch, so holes
would have to also be the highest terrain around to cause a problem.  However,
it does open the possibility of single pixels representing an entire foot
patch.

TBD Items
---------
TBD!!!!:  Improve handling of holes by computing how many pixels the foot
should cover at the given distance and have a percentage threshold to indicate
the maximum number of holes (i.e. minimum coverage of pixels) to allow for a
given location.  Perhaps also a threshold for the absolute minimum number of
pixels in a foot patch.  Throw out locations that do not meet this criteria.
This is complicated by the potential lack of a camera model (if an ortho XYZ
is input, which is a common use case) - a test would have to be made whether
a camera model was available, and assume a constant iFOV if not.

TBD!!!!:  Output labels are not being set correctly.


HISTORY:
  2015-06    S.Myint - Initial version of tilt and delta tilt programs
  2015-10    rgd - Nearly complete rewrite; merged tilt and delta tilt
  2015-11    rgd - Added clocking; added normal and Z outputs; fixed many
  		   other issues; wrote help.
  2016-01    rgd - Added outlier rejection, and caching of foot points
  2018-01    Steven Lu - Implement radial and cross-radial offsets for WTS.
  2018-07    Steven Lu - Added UPX_OUT and XPX_OUT.
  2019-12-23 Walt Bunch - Initialized some variables; cleaned up -Wall warnings.
                          Added unit test.
  2020-05-26 Walt Bunch - Replaced sprintfs.

COGNIZANT PROGRAMMER: B. Deen

PARAMETERS: