INTERFACE PlanetGeometry;
(*
  This module provides routines for converting between 2D pixel
  coordinates in a planetary image to 3D coordinates in a
  planet-relative frame.
  
  The module also provides tools for determining the 3D position and
  orientation of the planet, given the 2D position of its limb and a
  few reference features.
  
  Finally, the module provides tools for deducing the sun direction
  from the position of the terminator, and for computing the areal
  illumination and light direction at any point on the planet's
  surface.
  
  The geometric parameter are expressed in terms of three coordinate
  systems:
  
    IMAGE: coordinates for the image plane, in pixels, with origin at the 
    center of the top left pixel of the image, X axis pointing right, Y axis
    pointing down, and Z axis (when relevant) pointing out.  Note that
    this is a ""left-handed" coordinate system.
    
    CAMERA: coordinates for three-space, in km, with origin at the planet's
    center, X axis parallel to the IMAGE X axis, Y axis parallel to IMAGE Y
    but pointing up, and Z axis pointing almost towards the camera.
    Note that this is a right-handed coordinate system.
    
    PLANET: coordinates for three-space, in km, with origin at the 
    planet's center, with the Z axis pointing north,
    the X axis at longitude zero, and the Y axis in the right-handed
    direction.
    
*)

IMPORT LR2, LR3;
IMPORT Wr, Rd, ParseParams;

TYPE
  NAT = CARDINAL;
  LONG = LONGREAL;

TYPE
  T = RECORD
      diskCtr: LR2.T;      (* Center of disk, in pixel coords *)
      diskR: LONG;         (* Radius of disk, pixels *)
      planetER: LONG;      (* Planet's equatorial radius, km *)
      planetPR: LONG;      (* Planet's polar radius, km *)
      camSys: LR3x3.T;     (* Planet-to-camera transformation matrix *)
      camDist: LONG;       (* Distance from camera to planet center, km *)
    END;
    
PROCEDURE FromParams(
    pp: ParseParams.T;
    VAR f: ARRAY OF PlanetFeature.T
  ): T RAISES {ParseParams.Error} =
(*
  Extracts from command line options a description of the geometry
  of a planetary image.
  
  The user must specify the equatorial and polar radii of the planet,
  "planetER" and "planetPR", and the distance "camDist" from the
  camera to the planet's center; both in km.  (At present, the
  procedure requires that "planetER" and "planetPR" be equal, and
  "camDist" be infinite.)  The remaining parameters are computed by
  the program from user-given data.
  
  First, the procedure determines the center "diskCtr" and radius
  "diskR" of the planet's disk in the image.  This data is either
  given directly by the user, or computed from user-given pixel
  coordinates of three points along its limb.
  
  Next, the program determines a 3x3 orthogonal matrix "M" that
  relates the PLANET and CAMERA coordinate systems.  (More precisely,
  if "p" are the PLANET coordiantes of a point, then its CAMERA
  coordinates are "p M".)  The matrix "M" is either given directly by
  the user, or computed from the feature data "f", which must include
  the latitude, longitude, and pixel coordinates of at least three
  independent reference points on the planet's surface.  If these data
  are not available, the program sets "M" to the identity matrix.

  If there are features in "f" whose latitude, longitude,
  or pixel coordinates are undefined, this procedure will
  compute those missing data from the other date ???
  or are either 
  
  The given and computed data are all written to "stdout" in the
  format below.  (The comments in each line, from "#" to "\n", are
  not written):
  
    -disk <diskX> <diskY> <diskR>  # Image disk geometry, pixels.
    -planetRadius <plaER> <plaPR>  # Planet radii, equatorial and polar, km.
    -camDist <camDist>             # Spacecraft-planet distance, km; "0" means +oo.
    -planetAxes                    # Planet axes in image coordsys:
       <m00> <m01> <m02>           #   north axis;
       <m10> <m11> <m12>           #   "greenwhich" (zero longigude) axis;
       <m20> <m21> <m22>           #   "east" axis.
    -sunDir <sunX> <sunY> <sunZ>   # Sun direction in planet coordsys.
    -feature <name>                # This entry repeated zero or more times:
      <imgX> <imgY>                #   feature coords in image;
      <lat> <lon>                  #   feature's latitude and longitude;
      <plX> <plY> <plZ>            #   feature's direction rel. planet center.

  The following parameters are parsed and written, with
  default completion, but are otherwise ignored:

    -maxVal <maxV>                 # Maximum pixel value, as in PBM format.
    -minVal <minV>                 # Minimum significant pixel value.
    -lightVal <lightV>             # Nominal pixel value at subsolar point.
    -darkVal <minV>                # Pixel value for non-illuminated parts.
    -bias <BL>                     # Illumination bias, to displace terminator.
    -crop <uX> <uY> <NX> <NY>      # Useful image is "NX"*"NY", top left at "(uX,uY)".
    -planetR <planetER> <planetPR> # Equatorial and polar radii of planet, km.
    -camDist <camDist>             # Dist of camera from planet center, km; "0" means oo.
   
  TO DO: Extend computations to finite-distance perspective.
  
  TO DO: Extend computations to oblate spheroids.

*)
~

END PlanetGeometry.