Help for MARSXYZFOCUS

PURPOSE:
WATSON, just like MSL MAHLI, computes and returns depth-from-focus products such 
as the 'EDM' depth map along with the corresponding 'EZS' z-stack image.  The 
goal of this program is to use these maps to create an actual mesh, for which an
XYZ image needs to be created first and the run through the usual mesh processes.
The same process applies to MSL MAHLI, except for how searching is done, and 
also may apply to the ACI camera.

MARSXYZFOCUS takes a list of EZS/EDM parent images, typically eight, and 
first extracts their focus values. It then creates a DN to range curve based on
those values, given working distance values specified in meters. This curve is 
used to determine a ray which intersects a plane parallel to the instrument's 
sapphire window. The intersection point is saved as the output XYZ coordinate 
for each pixel.

EXECUTION:
marsxyzfocus inp=edm_image out=xyz_image parents=parentimages.list
 -OR-
marsxyzfocus inp=edm_image out=xyz_image focusvalues=array

where:
 - edm_image is the input EDM grayscale image
 - xyz_image is the output XYZ image
 - parentimages.list is a list of parent images for a given EZS/EDM pair
 - array is a comma-separeted set of focus values for each of the parent images. 
This array is optional, and should be used in case the input images themselves
are unavailable.


METHOD:
Given a set of parent images (typically 2-31) and/or their associated focus 
motor counts, this software first converts these values to range from the front 
surface of the lens, or the surface of the sapphire window for cameras such as
WATSON and ACI. In particular, for M20 the WATSON and ACI toolframes' origins 
are at this surface. The formulas for range conversion were obtained here:  
     https://link.springer.com/article/10.1007/s11214-021-00812-z
     https://zenodo.org/record/5555292#.Yi-qjBDMLZ4

More specifically, the coefficients for the WATSON case were obtained
from data acquired on both Mars and Earth and describe the relation between 
working distance (dw, originally in CENTIMETERS but converted to meters in the
program) and stepper motor count (mopen) when the WATSON 
dust cover is open and elements of the WATSON image are in focus. Those  
coefficients are:

	a = 1.09106×10^6
	b = –332.921
	c = 3.82592×10^(–2)
	d = –1.96922×10^(–6)
	e = 3.84562×10^(–11)
	dw = (a*mopen^(–1) + b + c*mopen + d*mopen^2 + e*mopen^3)^–1

For ACI, the formula is as follows:

	a = 0.005
	b = 20.34
	dw = a*mopen - b

As of now these values are hard-coded, but can be parameterized, preferably with 
something in $MARS_CONFIG_PATH, so it can be adjusted if needed in the future.

NOTE: Just to re-iterate since this is important, the centimeter values are 
converted to meters within the code, to be consistent with the units used in 
XYZ maps.

Next, there is a table which describes the relationship between the grayscale
data value (DN) in the EDM and which image corresponds to the the best focus,
which is a function of the nmber of parent images. For the i-th parent image out
of N in total, the DN value is computed per:
	DN = floor(255-i*255/N)

Now, we fit a curve which allows for the conversion of DN to range, given 
the DN vector as well as the correpsonding motor counts per parent image. 
Piecewise linear interpolation is an acceptable first-order solution, but this 
program performs a more-accurate 2nd-order polynomial fit. The fit is performed
via a least-squares solve using Eigen.

Finally, converting range (in m) to XYZ involves several steps:

1) Range is relative to the front surface of the sapphire window. This is the 
origin of the WATSON/ACI tool frame on M20. On MSL, the origin of the frame is 
where the camera optic axis intersects the plane defined by the top of the 
contact sensor pokers, which is 19mm (0.019m) in front of the lens cover. 
Alternatively, an arbitrary offset can be provided via the "PROJ_ORIGIN"
parameter. Given the projection origin and the 'A' camera orientation vector, 
also computed previously, we can create a plane that is coincident with the 
sapphire and is perpendicular to the optical axis. The origin and 'A' vector are
used in Step 5 below.
	 
2) For each pixel, we find the closest input image, using the DN vector, and 
extract the camera model from that image.  Currently, we have only one camera 
model across all focuses, but in the future we may get focus-dependent camera 
models, so we'll pull from the correct image. Also, camera model interpolation 
between focus values likely won't make a difference, so the closest value should 
suffice.

3) Use the chosen camera model to project a ray into space for the current pixel.

4) Use the computed range to offset the plane from (1) as defined by the origin
and 'A' vector along the plane's normal. 

5) Intersect the projected ray with the offset plane, which defines the XYZ 
coordinate for that pixel and save the XYZ coordinate in an XYZ image.

Finally, the XYZ coordinates for all pixels are saved as an XYZ image, from 
which a mesh can be created using for example 'marsmesh':

     $MARSLIB/marsmesh inp=xyz.VIC out=mesh.obj in_skin=image.png
	  - where image.png corresponds to the input EDM depth map, used for texture

ADDITIONAL NOTES:

The code does not current work well with thumbnails provided as parent images.
Another future improvement could be to handle cases where thumbnails are 
provided for most of the images but one full frame. That frame can be used to 
determine the coordinate system and other functions, but the focus values should
still be picked up from the rest of the thumbs.


COGNIZANT PROGRAMMER: Mauricio Hess-Flores

HISTORY:
  2022-04-01 M. Hess-Flores - Initial marsxyzfocus 

PARAMETERS: