Help for PHOTFUNC

PURPOSE:

PHOTFUNC corrects an image for spatial variations in apparent brightness caused
by illumination and viewing geometry.  This variation, known as the photometric
function, is removed to enable comparison of albedo differences on the target
surface.  In addition, the brightness correction enables the generation of
mosaics free from brightness differences along the edges of mosaic elements.

EXECUTION:

    PHOTFUNC INP=(INIMG[,GEOMA][,GRID][,CMAP]) [OUT=OUTIMG]  user-parameters...

  where: INIMG and OUTIMG are the input and output images,
         GEOMA is a geometric correction parameter file,
         GRID is a PHOTFIT grid file,
         CMAP is a classification map (see STATS, FASTCLAS),
	 [] brackets denote optional values.

Note: The grid input is currently deactivated (PHOTFIT is obsolete.
      See PHOTFIT2 and PHOTOM)

OPERATION:  PHOTFUNC can be used on virtually any image of a planet or
satellite, provided that the projection and lighting geometry are known.
The image may be a raw flight image (image space), one that has been corrected
for geometric distortions (object space), or one that has been map projected.
The image may be arbitrarily large and of any data format.

The projection and lighting geometry is normally obtained by via the SPICE
server or from the image label.  Alternatively, these may be provided via
parameters.

Because the computations can be fairly expensive, the function is computed
exactly for a grid of points spaced m samples and n lines apart (see SINC,
LINC, and INCR parameters).  For intermediate points within the grid, the
function is computed using bilinear interpolation.  Note that with increasingly
faster computers, the need for the interpolation is rapidly diminishing.  Use
the 'NOINTERP keyword to force exact computations at each pixel.

PHOTFUNC can also operate in a second mode, in which no photometric 
correction is performed, but information relating to the illumination 
geometry and photometry is calculated and printed (see NOCORREC and PRINT
parameters).

Further information is given under the following headings below:

  THE PHOTOMETRIC FUNCTION
  PROJECTION AND LIGHTING GEOMETRY
  GEOMETRIC DISTORTIONS
  COMPUTING THE OM MATRIX
  COMPUTING THE LIGHTING ANGLES
  INTERPOLATION ALGORITHM
  CLASSIFICATION MAP OPTION
  EXAMPLES

THE PHOTOMETRIC FUNCTION:  A model of the intrinsic reflectivity of a surface
is called a Photometric Function F:
		
	F(I,E,G,w,...) = B/J
where:
	B	is the observed reflected surface brightness,
	J	is the brightness incident on the surface,
	I	is the angle of incidence (i.e., of the sun),
	E	is the angle of emission (i.e., of the observer),
	G	is the phase, 
	w	is the albedo of surface particles, and
	...	denotes any other quantities which may contribute to the 
                functional relationship of B/J.
(See also HELP PHOTFIT2.)

PHOTFUNC uses such a model to correct the observed brightness to the
normalized brightness, B0, which is the brightness that would be 
observed with normal (90-degree) illumination and viewing angle
(I = E = G = 0).  Thus:

	B0/J = F(0,0,0,w,...)

and the correction factor f = B/B0 can be computed by dividing the above 
equations.  Then the apparent brightness B in the input picture is divided
by f in order to obtain an estimate of the normalized brightness B0 in the 
output picture. 

Currently, PHOTFUNC supports eight models for F:
  the Minnaert function:  F=B0*(COS(I)**k*COS(E)**(k-1),
    (the Lambert function is equivalent to Minnaert for k=1)
  three Hapke functions (1978 version, Cook modification, and 1984 version),
  the Veverka function, and
  the Mosher function (a variant of the Veverka function).
  the Buratti-Veverka function.
  the Irvine function.
See parameters MINNAERT, HAPKE, COOK, VEVERKA, BURATTI, MOSHER, and IRVINE.

Since the program actually only corrects by the normalized factor f, constant
factors, such as B0 in the first model, drop out and they are not needed
as input to the program.  (Thus Minnaert needs only 1 parameter in PHOTFUNC.)

When the gain introduced by the function is greater than a constant (see
MAXCOR parameter), no correction is applied.  All points which are uncorrected
for any reason (e.g. points off the target) are set to ZERO.

Planets also have a systematic brightness variation with solar phase
angle. Planetary science literature contains the magnitude phase
coefficients for satellites and planets. Burns (1977) "Planetary
Satellites" is one good source. If the geometric effects are removed,
the planet must still be brightened up by a factor obtained from the
magnitude phase coefficient depending on the phase angle. CAUTION:
The coefficients quoted in the literature usually refer to the
integrated disk of the planet. Care must be taken to account for any
phase defect (if the terminator is in the picture) which effectively
reduces the planet's area. By dividing the DNs in the picture by a
factor of (1+COS(G))/2 where G is the phase angle, the phase
defect can be removed. 

Also it may be necessary to determine the magnitude phase coefficient
if the picture was taken at a phase angle not previously seen and/or
different parts of the planet have different coefficients. The
coefficient can change over a range of phase angles. Callisto is a
good example. The leading hemisphere has a different coefficient than
the trailing hemisphere. Callisto's coefficient also changes at about
50 degrees phase. 

PROJECTION AND LIGHTING GEOMETRY

To remove the photometric function from an image, the projection and
lighting geometry for the image must be known.  Under normal circumstances,
this information is automatically retrieved from the image label or by
accessing the SPICE server.  The following description is intended to help
the user to provide any missing information via parameters.

The projection and lighting information is used to:

(1) relate each pixel in the image to a surface point on the target body
(2) compute the incidence, emission, and phase angles of light reflected off
    that surface point.

Both of these problems require a model of the target body.  PHOTFUNC uses a
tri-axial ellipsoid model defined by the three radii of the ellipsoid (See
RADII parameter).

The projection information required depends on whether the input image is a
raw flight image or one that has been transformed to a standard map projection.

For raw images, the following information is required (user parameters in
parenthesis):

  camera constants: focal length, pixel scale, line-samp coordinates of
    optical axis intercept (FOCL, PSCALE, LAXIS, SAXIS)
  geometric distortion parameters: see section on GEOMETRIC DISTORTIONS below.
  spacecraft vector: lat-long coordinates of the spacecraft and it's distance
    from the center of the target body (SLAT, SLON, RMAG, SPACE, RSVECTOR).
  camera pointing direction: transformation matrix from (rotating) target-
    centered coordinates to camera-centered coordinates (OMMATRIX).

For map projected images, the projection geometry is defined by the following:

  type of projection: e.g. orthographic, stereographic, mercator (PROJECT)
  picture scale in pixels/km (SCALE)
  special point line-sample (LINE,SAMPLE)
  special point lat-long (LATITUDE,LONGITUD)
  special parallels (PAR1,PAR2) for Lambert projection only

The special point is different for each projection.  See program MAP3 for a
description of the above parameters.

The lighting geometry for the image is defined by:

  spacecraft vector: same as above
  solar vector: lat-long coordinates of the sun (SOLAR).

The program first determines whether the image has been map projected by
scanning for a map label (see programs MAP3, MAPTRANS).  The required projection
geometry is normally provided in this map label.  However, the map label
currently does not contain the lighting geometry (spacecraft and solar vectors).
The program will attempt to retrieve these vectors from the SPICE server using
the same strategy as for raw flight images (see immediately below).  The user
may avoid SPICE access by supplying these vectors via the SPACE and SOLAR
parameters.

For raw flight images, PHOTFUNC will first scan the image label to determine
the target name, spacecraft ID, camera ID, and Spacecraft-Event-Time (see
parameters TARGET, MISSION, CAMERA, and SCET).  For a "supported" mission (see
parameter MISSION), the camera constants and geometric distortion parameters
are automatically retrieved from built-in tables (via a call to getcamcon).
See also the section on GEOMETRIC DISTORTIONS below.

If SPICE data is available for the mission (e.g. Voyager, Galileo, Cassini),
the projection and lighting geometry is automatically retrieved from the SPICE
server.  If not, that information can be supplied in the image label (see
program PERSLAB) or via parameters.  If the OM matrix is not know, see the
section on COMPUTING THE OM MATRIX below.

GEOMETRIC DISTORTIONS

If the input image is a flight image (i.e. not map projected), the program will
check if it has already been corrected for geometric camera distortions (by
scanning for the GEOMA keyword in the image label).  An uncorrected image is
is said to be in "image space" while a corrected image is said to be in "object
space" (see keywords DIST, IMAGE, OBJECT).  Note that "object space" is
identical to the perspective projection output by program MAP3.

If the image is uncorrected, a model of the geometric distortions must be
available so that it can be accounted for in the projection geometry.  For CCD
cameras (e.g. Galileo, Cassini), the geometric distortions are largely due to
the radial distortions of the optics.  For vidicon cameras (e.g. Voyager and
earlier missions) the distortions are defined by displacements in the reseau
pattern embedded on the camera's face-plate.  For each "supported" mission,
a nominal model of the geometric distortions is retrieved from built-in tables.
For CCD cameras, this model is adequate.  For vidicon cameras, more accurate
results may be obtained by locating the reseau in the image and inputing the
resulting GEOMA parameters (see parameter INP and programs RESLOC and RESLOCVO).

   resloc f1636832.fic (res,geo)
   photfunc (f1636832.fic,geo) a.img HAPKE=(.5,.2,10.,.5,.2) target=IO

A raw Voyager image containing geometric distortions is first input to program
RESLOC to locate the reseau and compute geometric correction parameters (geo).
These parameters are input to PHOTFUNC together with the image (the
image must be first) and are used to account for the geometric distortions.

   resloc f1636832.fic (res,geo)
   farenc (f1636832.fic,geo) (c,p,g) target=IO
   geoma (f1636832.fic,g) b.img
   photfunc b.img a.img HAPKE=(.5,.2,10.,.5,.2) target=IO       sscpt=(400.,400.)

This is similar to the previous example, except that program FARENC is used to
correct the camera pointing and output GEOMA parameters (g) which center the
target in the image and correct for geometric distortion.  GEOMA is then used to
create an object-space image.  The result is input into PHOTFUNC, which
recognizes that the image is in object-space from the GEOMA keyword in the
image label.  The SSCPT parameter is required to indicate that the target
center has been translated to the center of the image by FARENC.  Note that
if an object-space output is desired, then running GEOMA before PHOTFUNC is
more effecient since PHOTFUNC runs faster if it does not need to account for
geometric distortions.

Currently, no nominal models are available for Space Telescope or SIPS.

COMPUTING THE OM MATRIX

If the OM matrix is not available, it may be computed via the (1) far-encounter
or (2) tiepoints algorithms.  Both of these methods require that the remaining
projection information be known.

The far encounter algorithm is used in images in which the target limb is
visible in the image.  VICAR programs FARENC or NAV is used to fit the
limb and obtain the line-sample coordinates of the target center.  Parameters
SSCPT, OSSCPT, or ISSCPT may be used to input these value to PHOTFUNC.  This
algorithm also requires the north angle (orientation of the target's spin axis
measured clockwise from up in the image).  See NORANGLE parameter.
The program inputs these values into the far encounter algorithm (see VICAR
subroutine MOMATI) to compute the OM matrix.

The tiepoints algorithm is used in high-resolution images in which the target
limb is not available.  The image is first paired with a reference image whose
camera pointing is known (see programs MANMATCH, NAV2).  Common features
(tiepoints) are located in both images (by determining their line-sample
coordinates) and the corresponding latitude-longitude coordinates are computed.
These (line,sample) and (lat,long) pairs are input to PHOTFUNC via the
TIEPOINTS parameter.

Note that a common use of the far encounter and tiepoints algorithms is to
improve the nominal camera pointing initially available from the SPICE server.
The programs referred to above normally store their results into SPICE C-kernels
from which they are automatically retrieved by PHOTFUNC.  The FARENC and
TIEPOINTS parameters are therefore necessary only if the SPICE server interface
is circumvented or unavailable.

COMPUTING THE LIGHTING ANGLES

To compute the lighting angles at a given pixel, the programs must transform
that pixel to a surface point on the target, and use the spacecraft and solar
vectors to compute vectors that point from that surface point to the sun and
spacecraft.  In addition, the target body model (polar and equatorial radii) is
used to compute the surface normal at that point.

The incidence angle is the angle between the surface normal and the vector
pointing to the sun.

The emission angle is the angle between the surface normal and the vector
pointing to the spacecraft.

The phase angle is the angle between the vectors pointing to the sun and
spacecraft.  Note that since the geometry is three-dimensional, the phase angle
is normally not equal to the sum of the incidence and emission angles.

INTERPOLATION ALGORITHM:  The photometric function is computed exactly for a
grid of points whose spacing is specified by parameters LINC, SINC, or INCR.
For intermediate points within the grid, the function is computed using
bilinear interpolation.  In general, a smaller grid spacing reduces the error
introduced by the interpolation, but increases the execution time.

If one or more of the four grid points surrounding a given pixel is off the
target body, the photometric function will be computed exactly for every pixel
within the box within the four grid points.  This criteria can be extended to
apply to points within LIMB degrees or the limb or TERM degrees of the
terminator (see LIMB and TERM parameters).  This strategy for avoiding large
interpolation errors near the limb or terminator may not be the best, as shall
be clear below.

The keyword 'NOINTERP causes the photometric function to be computed exactly
at every pixel in the image.  In this mode, the grid (LINC, SINC) is used to
skip efficiently over dark sky regions of the image.

It is instructive to compute the error introduced by the interpolation
algorithm by photometrically correcting a full-disk image with and without
interpolation:

	photfunc f1636832.geo a.img MINN=0.7
	photfunc f1636832.geo b.img MINN=0.7 'NOIN
        f2 (a.img,b.img) c.img func="100*in1/in2"
        fit c.img 'EHIST

The error introduced near the limb of the target is significant.  For high
resolution images (no visible limb) the error introduced by interpolation is
negligible.

Keep in mind that the interpolation algorithm was developed in the '70s when
computers were 1000 times slower and 1000 times more expensive.  In an age of
more powerful computers, the interpolation algorithm seems less relevant.  The
following timing results were obtained on a Sun Ultra 5 workstation:

Timing data using interpolation with a 10x10 grid spacing (the default):
					CPU sec		Wall clock sec
Galileo 800x800 Image Space image:	 1.63		 3.52
Cassini 1024x1024 Image Space image:	 1.84		 4.77
Voyager 800x800 Image Space image:	 3.96		 5.35
Voyager 1000x1000 Object Space image:    1.57		 3.44
500x500 Rectilinear projection:		 1.49		 2.70

Corresponding timing data using the 'NOIN keyword:
					CPU sec		Wall clock sec
Galileo 800x800 Image Space image:	 8.50            9.80
Cassini 1024x1024 Image Space image:	15.33		17.00
Voyager 800x800 Image Space image:	25.22		26.63
Voyager 1000x1000 Object Space image:    4.84		 6.09
500x500 Rectilinear projection:		13.44		15.04

Also, some of the photometric functions are computationally expensive if
calculated for every pixel.  The following are timing results using the 'NOIN
keyword on a 1024x1024 Cassini image space image:
						CPU sec		Wall clock sec
Minnaert:					15.33		17.00
BURATTI:					17.13		19.23
VEVERKA:					15.44		17.20
Old HAPKE function:				15.87		17.44
HAPKE function with COOK modification:		15.86		17.39
HAPKE with one term Henyey-Greenstein:		25.72		27.18
HAPKE with two term Legendre polynomial:	24.93		27.20
IRVINE:						16.99		18.95

CLASSIFICATION MAP OPTION:

It is possible to have PHOTFUNC apply its photometric correction only
to certain points in the image, by supplying a classification map (such
as can be generated by program FASTCLAS) as a supplementary input file, 
and specifying a class by the CLASS parameter.  Then only pixels which 
have the specified DN value in the classification map are corrected.  
This can be useful if different types of terrain in an image are to have 
different photometric parameters, as would be the case for the more 
physically realistic functions, such as those of Hapke.  

The following example illustrates the use of the classification map option.
Note that the example fails because of this programmer's ignorance in
using STATS and FASTCLAS.  However, the method should work in principal.

fit io.org a.img (1,1,450,450) 'byte
fit io.blu b.img (1,1,450,450) 'byte
fit io.vio c.img (1,1,450,450) 'byte

stats (a.img,b.img,c.img) stats.img exclude=0    class1=(188,81,20,20) class2=(369,80,20,20)
fastclas (a.img,b.img,c.img,stats.img) class.img sigma=(20.,10.)

copy io.vio c.img (1,1,450,450)
photfunc (c.img,class.img) out.img IRVINE=(1.14,.118,.0019) CLASS=2     'CYLI SCALE=8. NORTH=0.     LATI=80. LONG=230. LINE=225. SAMP=1.
   SOLAR=(0.541,171.276) SPACE=(-.032,156.474,806030.)  
A 450x450 area of a color triplet is extracted from an IO global mosaic and
convert to byte data (program STATS will only accept byte inputs).  Programs
STATS and FASTCLAS are run to create a classification map (class.img) which is
then input to PHOTFUNC.  Only pixels from Class 2 of the map will be corrected.
All other pixels are output as 0 DN.  Note that the image input to PHOTFUNC
must represent the same image area as that of the classification map.  This is
accomplished by copying the area from the mosaic BEFORE inputing into PHOTFUNC.

This example also illustrates how to enter the projection and lighting geometry
for map-projected images if these are not included in the image label.

As the program is currently structured, a separate run must be made for 
each class, with the output of one run forming the input for the next, 
until all classes have been processed.  (The option to do all classes
except for a specified one is also supported.)  It should be possible 
for a future modification to process all classes in one run.

EXAMPLES:

!--------------------------Galileo image-------------------------------

  photfunc inp=18494401.1 out=a.img MINN=0.7

A Minnaert function is applied to a Galileo image of Callisto.  PHOTFUNC will
retrieve the mission, camera id, Spacecraft-Event-Time, and target (Callisto)
from the flight label. These image identifiers are used to automatically
retrieve the image and lighting geometry via the SPICE server.  If everything
goes well, running PHOTFUNC can be as easy as pi.

!--------------------------Cassini image-------------------------------

photfunc n1354897340.1 a.img BURATTI=(0.5,.6,-.003,.14,.14,1.0)     SSCPT=(512.,512.)

Assume that the input image has been previously processed by program FARENC to
center the target body in the image.  The SSCPT parameter specifies that the
center of the target is at line-sample coordinates (512,512).  This will result
in a re-computation of the OM matrix (see description of far encounter algorithm
above).

!--------------------------VGR image space image-------------------------------

photfunc f1636832.fic a.img HAPKE=(.5,.2,10.,.5,.2) target=IO

A raw Voyager image containing geometric distortions is input to PHOTFUNC.
Since no geometric distortion parameters are provided, nominal values are
retrieved from built-in tables and used to account for the geometric distortions
in the projection geometry.  Note that the TARGET parameter is required since
the target body is not identified in Voyager flight labels.  See the section
on GEOMETRIC DISTORTIONS above.

!------------------What to do when SPICE data is not available-----------------

photfunc f1636832.fic a.img BURATTI=(0.5,.6,-.003,.14,.14,1.0) 'NOSPICE      TARGET=IO SOLAR=(0.541,171.276) SPACE=(-.032,156.474,806030.)      SSCPT=(539.67,601.21) NORANGLE=15.7

Assume that the input image has a Voyager flight label.  The 'NOSPICE keyword
prevents the program from attempting to retrieve data from a non-existent
SPICE server.  Based on the mission (VGR-1) and camera identifiers in the label,
the camera constants and nominal geometric distortion parameters are retrieved
from built-in tables.  Based on the TARGET parameter, the polar and equatorial
radii are retrieved as well (careful, these may not be the latest or greatest).

Assuming that the image label contains no further useful information, the
solar and spacecraft vectors must be input via parameters.  Finally, the
subspacecraft line-sample coordinates (a.k.a target center) and north angle are
input to enable the program to compute the camera pointing (see description of
far encounter algorithm above).

photfunc f1636832.geo a.img VEVERKA=(0.5,-0.01,0.5,0.01) 'NOSPICE
    RADII=(1824.3,1824.3,1815.7)      SOLAR=(0.541,171.276) SPACE=(-.032,156.474,806030.)      TIEPOINT=(381.86,382.64,19.35,229.21,          382.17,498.94,25.83,190.00,          381.98,615.31,32.86,163.58,          498.40,498.49,2.539,179.32,          498.53,615.11, 8.65,156.13)

The above example demonstrates the use of the tiepoints mode to specify the
camera pointing for a high-resolution Voyager image (see description of
tiepoints algorithm above).

!------------------------Unsupported missions---------------------------------

  photfunc inp=a.img out=b.img VEVERKA=(0.5,-0.01,0.5,0.01) 'OBJECT      TARGET=IO FOCL=1500.19 LAXIS=500. SAXIS=500. PSCALE=84.8214      SOLAR=(0.541,171.276) SPACE=(-.032,156.474,806030.)      SSCPT=(539.67,601.21) NORANGLE=15.7

If the image is taken from a spacecraft which is not recognized by PHOTFUNC,
all projection and lighting geometry information must be supplied via
parameters.  Here the target constants are specified by identifying the target
body, and the OM matrix is computed from the target center and north angle
(SSCPT and NORANGLE).

!--------------------------Map projected image---------------------------------

map3 &"path1"f1636832.geo b.img NL=500 NS=500     'RECT TARGET=IO SCALE=10.     LINE=1. SAMP=1. LATI=80. LONG=230.
photfunc b.img a.img IRVINE=(1.14,.118,.0019) TARGET=IO

A Voyager image of IO is mapped to a rectilinear projection via program MAP3.
The resulting projection is input to PHOTFUNC.  Although the projection
information can be retrieved from the "map label" output by MAP3, this label
does not currently contain the lighting geometry.  PHOTFUNC will attempt to
retrieve that information from SPICE.  Since both MAP3 and PHOTFUNC will
access the SPICE server (in this example), and since the Voyager label does
not include the target name, this must be provided via the TARGET parameter.

RESTRICTIONS:

1. The GRID option is currently disabled.

HISTORY

Written By: Joel Mosher			1 AUGUST 1978
Cognizant Programmer: L.W.Kamp
Revision: 1				9 FEBRUARY 1981
Revision: 2     CONVERT TO VICAR2      16 MAY 1985
Revision: 3     Viking additions       11 MAY 1986
Rev. 4:     Mods for new PHOTFIT.  2 APR 1987  -- LWK.
Rev. 5:     Added classificaion map option,  19 Aug. 1987  -- LWK.
Rev. 6:     Use normalized phot. func., 2 Sept. 1987 -- LWK
Rev. 7:     Add 6-parameter Hapke fcn., bug fixes, 8 Oct. 1987 -- lwk
Rev 8       Many bug fixes. Put in exact photometric function at
            every point for limb & terminator.
23 Aug 89  ..GMY..  Delete fiddling FDS for WA frames of simultaneous exposure
15 Aug 93   jjl ability to read map3 perspective label type.
06 Sep 02  GMY: Major revisions to help and test files.
	   Major revisions to projection and lighting geometry determination
	   Add NAIF keyword and delete obsolete SEDR keyword (AR 107288)
           Updated to read IBIS format geometric distortion files (AR 107368)
           Changed NOSEDR keyword to NOSPICE.
	   Deleted GEOMA, NAH, NAV, and START parameters.
	   Deleted LUMS and THRESH parameters.
           Logic controlling line loop and interpolation redesigned.
           Fix bug in farenc algorithm.  Fix to work on floating point images.
	   Fix computation of surface normal vector.
07 Oct 02  Replaced oblate spheroid model with triaxial ellipsoid.  All calls
           to convev have been replaced by calls to mp routines.
           Added 'NOIN parameter.
10 Jan 13 -lwk- fixed CHARACTER continuation lines for new compiler flag on Solaris

PARAMETERS: