Help for FARENC

 PURPOSE:
 "farenc" is a Vicar program to determine the center of a target body
 presenting part of its limb in a picture.  The orientation
 and lighting geometry of the target body can be made constraints.
 The program has 4 main modes. These are:
 
 1. Determining constraints from the SEDR or from parameters.
    Constraints determine the shape,size,and orientation of the limb.
 2. Locating candidate limb points.
 3. Fitting these points to a conic.
 4. Generating GEOM parameters, 
    Updating the SPICE/SEDR OM matrix, and
    Writing the planet center to a file (if the input is a perspective
    projection with a map3 label)

 Three limb fitting functions are available:

              1. An unconstrained circle
              2. A constrained circle
              3. A constrained ellipse

LOCAL AND REMOTE SPICE ACCESS:

FARENC accesses SPICE data either locally or via the MIPS SPICE server.
Currently (Jan 2002) SPICE data is available for Cassini, GLL, and
VGR (at Jupiter only).  For images without SPICE data, parameters may be used
to specify the target geometry.

SPICEMODE specifies whether local or remote SPICE access is to be used.
If defaulted, SPICEMODE is set to the value of the logical name (or
environmental variable) DEFAULTSPICE.


PARAMETERS FOR RETRIEVING THE INITIAL CAMERA POINTING:

Initial camera pointing data is first retrieved from predict C kernels or
from MIPS C kernels.  The following optional parameters permit the user to
specify where this initial pointing is retrieved:

CKNAME and CKID are alternative ways to specify the C kernel from which camera
pointing is to be retrieved.  For example, CKNAME=FARE or CKID=M904 specifies
that the camera pointing should be retrieved from the file MIPS_FARENC.CK.
The CKID is the unique kernel ID assigned to each C kernel.  When CKID is
specified, it overrides the CKNAME parameter.  A complete list
of kernel IDs is located in the ASCII file assigned the logical name (or
environmental variable) KERNELDB.


PARAMETERS FOR STORING THE IMPROVED CAMERA POINTING:

The following optional parameters are used to store provenance information along
with the improved (C-Smithed) camera pointing computed by the program.  This
provenance information is stored in the (C kernel) segment identifier for the
camera pointing data.  In cases where there are several versions of camera
pointing for an image in a given C kernel, this provenance information can
later be used to retrieve a specific instance of camera pointing from the
kernel.

PURPOSE identifies the purpose for creating the camera pointing.
REQNUM identifies the request number associated with the camera pointing.
CDATE specifies the date and time the camera pointing was created.
GROUPID identifies the group which created the camera pointing.
INSTITUTE identifies the facility which created the camera pointing.

See the level 2 help (via the TAE tutor mode) for further details.

Examples:  'LOCAL CKNAME=NAIF specifies that SPICE data be retrieved from
          local kernels using camera pointing from predicts or AACS telemetry.
          Improved camera pointing will later be stored in the local C kernel
          specific to this program.

           'REMOTE INSTITUTE=DLR USERID=TYR retrieves SPICE data remotely
          (from MIPS) via the SPICE server.  When improved camera pointing is
          stored (at MIPS), provenance data regarding the facility (DLR) and
          user (Thomas Roatsch) who created the data is included.


 "farenc" determines limb points by computing a value called
 activity.  Activity is defined as the sum of the absolute
 values of the differences between the DN's of opposite
 sides of a 3 by 3 pixel box, i.e.:

          --- --- ---
         | 1 | 2 | 3 |
          --- --- ---
         | 4 | 5 | 6 |    |DN1 - DN9| + |DN3 - DN7| = Activity
          --- --- ---
         | 7 | 8 | 9 |
          --- --- ---
 Points having the highest activity lie along the limb of a
 target body.  If ACTI is the activity threshold and DNTH is the DN
 threshold ( As applied to DN in the above figure) then the
 tests required for a point to be accepted as a limb candidate are:

       1. DN5 >= DNTH
       2. |DN2 - DN8| + |DN4 - DN6| >= ACTI
       3. |DN1 - DN9| + |DN3 - DN7| >= ACTI
       4. MIN (DN1,DN2,DN3,DN4,DN6,DN7,DN8,DN9) <= BELOW
       5. The direction of the planet center obtained from:
         atan2(dn1+dn2+dn3-dn7-dn8-dn9,dn3+dn6+dn9-dn1-dn4-dn7)
          must be within 90 degrees of the sun illumination 
          direction.
 For each line and column, two limb candidate points are chosen.
 These points must be the two largest activity values in that
 dimension which:

       1. are more than or equal to DIST pixels apart.
       2. form a ratio of smallest/largest activity at least
           equal to HEIGHT.

 Points are then removed which fail to have sufficient numbers
 of neighbors, implying that they are isolated events.

 Limb fitting is performed in two stages. First an initial fit is
 made using simulated annealing. Then a least squares fit is performed
 which permits iterative point rejection until all the points
 lie within TOLERANCE.
 
 1. Simulated annealing. This is a complex method which treats
    the problem like a thermodynamic cooling system and whose
    object is to cool the system into the lowest energy state.
    The lowest state corresponds to the minimum residual fit
    error for the limb points.

 2. Iterative least squares fit.
    Each cycle rejects values which depart from the
    mean error by more than +/-SIGMA * THETA, until all values fall
    within computed radius +/-TOLE.

 If an oblate spheroid is imaged from a distance less than about
 20 diameters and if the oblateness is greater than about 5%,
 the center of the target body may not coincide acceptably with the
 center of the projected ellipse.  An approximate solution to
 the problem ( assuming the projected shape is an ellipse )
 which provides an additive offset to the computed ellipse
 center to give the true target body center is provided by "farenc".
 This correction is not made to the ellipse center, but is
 only printed along with the ellipse center.  A target body whose
 oblateness is 10%, viewed from two diameters at a subspace-
 craft latitude of +/-45 degrees will suffer a target body center
 discrepancy error from that of the ellipse center of about
 1% of the ellipse major axis.

 Unless data specific to the constrained ellipse or circle
 has been provided, "farenc" will compute only the general
 circle.  Unless 'NOSEDR is specified,
 data for the constrained fit will be obtained from the SEDR
 or SPICE.  User specified parameters will over-
 ride the equivalent SEDR values.

 If a reseau or GEOMA data has been provided, the ellipse
 and circle centers will be computed in object space by
 first transforming all of the limb points to object space.
 The principal output file is a dataset containing GEOMA parameters.
 (This file must always be the LAST output.)  This data set can be
 input to GEOMA to geometrically correct and translate the target body
 such that the center of the target body will lie at NL/2.0, NS/2.0.
 This data set could also be input to other programs as a
 geometric correction data set.  These other programs such
 as "photfunc", "photcat", "photcal" and "map2" must have "farenc"
 modes for determining the OM matrix (the camera to picture
 body rotation matrix).  The object space target body center
 will be NL/2.0,NS/2.0 by default.
EXECUTION:

 "farenc" may be executed using the following format:

    farenc INP=(IN,GEOM) OUT=(CURVE,POINTS,G) SIZE=(SL,SS,NL,NS) PARAMS

  where PARAMS represents parameters such as those that follow.

SUGGESTED EXECUTION FORMAT:
   farenc inp=(in,geom) out=g 'auto

Note: If the input image has a map3 perspective label placed there by
      the "perslab" program "farenc" will ignore the sedr/spice and
      will write the planet center to an output file called: FARENC.POS
      This file is a Vicar image, has 1 line & 2 samples with the
      real values line,sample. See test file for example using
      Space Telescope project imagery.

 The following parameters perform editing operations on the
 initial set of limb points in order to reject spurious points
 from the final least squares fit:

   	SIGACT  SIGMA  TOLERANCE   SUNANGLE  PRINT   MAXNUM

 For additional information type help followed by the parameter name.

 The following parameters control the nature of the geometric
 transformation parameters ouput in data set "G".

       GEOM              RESCALE

 For additional information type HELP followed by the parameter name.

 The following parameters allow the constrained ellipse to be
 used provided that either of the two sets of parameters which follow
 are completely specified.  In the first case the center of the
 projected ellipse will be computed but not the correction from
 ellipse center to oblate spheroid target body center ( they can disagree ).
 If SEDR data are available, the values for SET #II will be
 automatically computed. If data from the SEDR are not desired or available
 the 'NOSEDR keyword should be used. The default is to read the
 SEDR (Image Catalog).

 SET #I
    NORANGLE or ANGLE       SMAA              SMIA

 SET #II
          NOSEDR
    NORANGLE or ANGLE FOCL or FOCAL SCALE  RADIUS
    REQUATOR    RPOLE FAR   RMAG  SLAT   LATI  RING
    LATI  LONG    SOLAR  PLANET

 For additional information type HELP followed by the parameter name.
 Program Operation Under Annealing Step.

 Each of the three functions (general circle, constrained circle &
 constrained ellipse) is allowed 6 attempts to bring N % of the
 limb points within 1 pixel of mean radial error. If all six
 attempts fail the smallest residual case is selected.
 The residual returned & printed will be either:
 1. The mean radial error if N % of points are within 1 pixel of 
    the fit.
 2. The mean error + (smaa/5)/ (# points selected within 50 
    pixels radial error)
 For each attempted fit N becomes: 80,70,60,50,40,and 30 %.   

 The annealing cost or objective function to be minimized is:
 (if n= total # pts, m=# points within smaa/5 radial error,
  k=# points within 3 pixels radial error, sumdr=sum radial
  residuals)

  if(k > N% of n)then cost=sumdr/k
  else if(m > 0)then cost=sumdr/m + (smaa/5)/m
  else cost=sumdr/n + 1

  The cost function is designed to punish for poor fits
  yet reward for including more points. It's a tradeoff.

 PRECISION: 
  Some differences in intermediate results such as shown below might
occur between successive runs of FARENC:
*****************
No  convergence,residual=  22.380320

Constrained annealing circle fit object space:
Line of center=  282.28452
Samp of center=  410.81024
Radial fit mean residual=  22.374784
******
No  convergence,residual=  22.374847

Constrained annealing circle fit object space:
Line of center=  282.29230
Samp of center=  410.80789
Radial fit mean residual=  22.374775
******************
  The differences above come from the first stage of the planet center
computation (with subroutine METROPOLIS).  This stage uses random numbers
in searching for the planet center location which best agrees with the points
found on the planet limb.  Because of the random numbers, the results of the
first stage vary slightly from run to run.
  The second stage of the planet center computation, the least squares
solution, refines the results of the first stage and produces line and sample
coordinates that agree well from run to run and machine to machine.  The
"Sigma of Radius Errors" and the "Largest Radius Error" often will vary
slightly for different machine types.
  In the test log for FARENC, there can be slight differences in output from
program GEOMA for different machine types.  Likewise, there are slight
differences in output for different machine types when FARENC is run on an
image written by program GEOMA.  This is due to round-off building up over
multiple processing steps.

Also, updates to the SPICE kernels may cause differences in the section showing
FINAL DATA FROM SEDR and PARAMETERS.

PROGRAM HISTORY:

 WRITTEN BY: J. J. Lorre         December 1, 1976
 CONVERTED TO VAX BY:  Joel Mosher,  23 Aug. 1983
 COGNIZANT PROGRAMMER: J. J. Lorre
 REVISIONS:
  1/25/81 -JAM-UPDATE VGR SEDR WITH OM MATRIX AND RS VECTOR
  1/22/82 -JAM- REPLACE WATONA WITH MWATNA
  9/29/82 -JAM- PUT BACK IN MESSAGE FOR NUMBER OF POINTS USED
  10/6/82 -JAM- PUT IN CHANGE TO RESET 'OFFSET' WHEN AUTO IS USED
  10/9/82 -JAM- PUT IN NEW FIX FOR OFFSET USING # PIXELS ON PLANET
  10/15/82 -JAM- PUT IN CODE TO SAVE GEOMA PARAMETERS
  10/25/82 -JAM- REPLACE CALL TO SEDR79 WITH GETCDR
  10/31/82 -JAM- FIX CONVERSION BACK TO GEN. ELLIPSE IN IMAGE SPACE
  2/21/83 -JAM- REPLACE CALLS TO R4SORT WITH RBSORT
  8/17/83 -JAM- REPLACE CALLS TO RBSORT WITH R4SORT
  5/10/84 -CCA- DEL LABEL IO & SYS000 TO COMPUTE PLANET CENTER AND RADIUS
  11 Feb. 1985  - L.W.Kamp -  Get radii from PBDATA, fixed bug in IG
                array size, removed parm INTERACT, added parm FRACTION,
                cleaned up code.
  4/12/85 -JAM- CONVERT TO VICAR2
  5/9/85  -JAM- PUT IN 'NOSEDR KEYWORD
  5/12/85 -jam- make "nogeom" the default
  10/24/85 -CCA- FIXED OVERFLOW BUG IN LSQING, COMMENTED OUT NEW AUTO CODE
  86-2-11 -LWK- replaced hard-coded radii by call to PBDATA; also
                used PBDATA for TARGET parameter;
                fixed overflow of IG;  removed parm INTERACT;
                removed subs CURSOR, WATONA, PRNT2, RESO;
                removed DSRNSE, NOGEOM(always=1); changed COUNT of INP to 2;
                revised treatment of parameter defaults;
                restored JJL's AUTO code, controlled by new parm FRACTION.
 9/29/87 JJL 1. Introduced Metropolis annealing function.
             2. Introduced terminator rejection based on Sobel edge
                direction relative to the sun.
             3. Changed excess point truncation to alternate point selection.
 01 Oct 1987 -  J. J. Lorre - extensive modifications:
                 1. Selection of points only on solar illuminated limb
                    by computing the angle of the local edge.
                 2. Introduction of annealing code to avoid
                    instabilities in least squares routines when
                    bad data is present. 'OLD keyword.
 5/9/88  -FFM-  Incorporate 'SEDRSRC keyword
 10 Jun 1988 -  G.Yagi - fix determination of planet from ERT.
 10 Jun 1988 -  G.Yagi - fix test of ABLE77V2 indicator.
 1/10/90  JJL   Converted to project independent subroutines.
 2/1/93  jjl  Placed annealing and least squares methods in sequence.
 15 Aug 93 jjl  Accepts map3 perspective input label.
 15 APR 94 CCA  Add output to TCL variables PCL, PCS
 10 Jul 95  ams  (CRI) Made portable for UNIX
 11 Sep 95  gmy Changed GETSPICE call to use remote SPICE server (FR 85898)
 27 Mar 96  SP  Changed parameters and method for obtaining SPICE data according
                to new method specified by Gary Yagi using GETSPICE2/PUTSPICE2
                instead of GETSPICE/PUTSPICE.  Old parameter SEDRSRC is replaced
                by parameter CKNAME.       (Note that
                the ported version does not support a second input file.)
                Also added some VALID range checks in PDF for extreme value
                handling.
 12 jun 96  SP   Add code to distinguish between OBJECT SPACE and IMAGE SPACE.
                 (IMAGE space is raw image.  Object space is corrected for
                 camera geometric distortion.)  This corrects an error for GLL
                 where coordinates of limb points not converted to OBJECT
                 SPACE for calculations of planet center.
 22 Jul 96  SP   Upgraded to handle IBIS tiepoint files and other mission
                 besides GLL (allow second input file).  Put back in CCA's
                 option for outputting planet center to TCL variables.
                 JJL made a major fix to correct first stage (METROPOLIS)
                 computation of planet center.
 18 Aug 96 lwk  Fixed bugs in the case of Point Perspective input label;
		added support for NIMS labels (subsolar point)
 20 Dec 01 GMY  Add Cassini capability (CR 7953-MIPS).  Major changes to
                test script.
 20 Jul 05 lwk  Added another check in LSQP to prevent infinite loop when a matrix
                element becomes zero

PARAMETERS: