Help for GEOMA


GEOMA is a VICAR applications program which makes geometric changes in
pictures.  It can be used for correcting geometric distortion, or performing
other types of geometric transformations.  The geometric changes
(transformations) are defined as a mapping from a set of locations (called
tiepoints or control points) in the input image to a corresponding set of
locations in the output image.  Between the tiepoints, the mapping is defined
by bilinear interpolation.  GEOMA is similar in function to VICAR program

The following TAE command line formats show the most common usages:
      GEOMA INP=(a,a2) OUT=b SIZE=(1,1,nl,ns) optional parameters
      GEOMA INP=(a,a2) OUT=b  NL=nl NS=ns optional parameters
      GEOMA (a,a2) b (1,1,nl,ns) optional parameters
      GEOMA (a,a2) b optional parameters
      GEOMA INP=a OUT=b SIZE=(1,1,nl,ns) optional parameters
      GEOMA INP=a OUT=b  NL=nl NS=ns optional parameters
      GEOMA a b (1,1,nl,ns) optional parameters
      GEOMA a b optional parameters
      GEOMA INP=a OUT=b SIZE=(1,1,nl,ns) PARMS=c optional parameters
      GEOMA INP=a OUT=b  NL=nl NS=ns PARMS=c optional parameters
      GEOMA a b (1,1,nl,ns) PARMS=c optional parameters
      GEOMA a b PARMS=c optional parameters

       Here 'a' represents the input image file name,
       'a2' represents an IBIS input file containing tiepoint information,
       'b' represents the output image file name,
       and 'c' represents the parameter file name.

       Note: GEOMA no longer supports a non-IBIS second input file.  
             Old non-IBIS 'GEOMA' files (which were often used to specify
             the geometric transformation needed to correct for camera
             distortion) can be converted on the VAX or Alpha using the program

In both GEOMA and LGEOM the picture to be geometrically transformed 
is subdivided into a number of four-sided areas (quadrilaterals).
The user must specify to the program the exact mapping from input picture to
transformed picture for each of the four corner points of each quadrilateral.
Within the quadrilaterals, the programs use a bilinear interpolation to map
the interior points from input picture to transformed picture. The LGEOM 
program restricts the specification of the mapping of the quadrilateral
corners (control points) such that, in the transformed (output) picture, they
form a rectangular grid, (i.e. the quadrilaterals are rectangles in the
output image.)  The GEOMA program does not impose this restriction, and 
allows more arbitrary specification of the control point mapping. In fact, 
two of the four control points forming the corners of any given quadrilateral
may coincide, causing the quadrilateral to degenerate into a triangle (a 
quadrilateral having one side of zero length). The potential user is WARNED 
that GEOMA can produce some undesirable side effects due to the basic 
difference from the LGEOM program.  This is explained further in the 
OPERATION section.

As with LGEOM, GEOMA will accept its control parameters as the second input
file or as a parameter file instead of or in addition to the conventional 
VICAR parameters.

The size of the output image is determined by the number of lines and number 
of samples in the SIZE field if the SIZE field is entered by the user.  If the
SIZE field is not entered, the output file is the same size as the input file.

The data type of the input image may either be byte or halfword data. The data
type is obtained from the VICAR label of the input image or from the FORMAT
parameter. The output image has the same data format  (byte or halfword) as the
input image. 

The tiepoint grid must be specified as a set of parameters in the
TAE command line or in the parameter file or in a second input file.
The tiepoint grid is specified using the parameters NAH, NAV, and TIEPOINT.

   NAH=nah  NAV=nav
     TIEPOINT=(nl1,ns1,ol1,os1   nl2,ns2,ol2,os2...

where nah is the number of grid cells horizontally (across the top in the
sample direction), nav is the number of grid cells vertically in the output
image space, the point pairs (nli,nsi,oli,osi) are line-sample coordinates in
the output (new) and input (old) spaces respectively. The number of pairs k
must equal (nah+1)*(nav+1). 

The RESTART parameter has been removed since it was deemed not maintainable.
(See FR 16750 of February 1986.)
Applications of GEOMA:

GEOMA and LGEOM are both VICAR programs for performing general geometric
transformations.  Both of them require a grid of tiepoints for defining the
transformation.  In LGEOM, the grid must be rectangular (i.e. the tiepoints
are arranged in straight rows and columns.)  In GEOMA the grid does not need
to be rectangular. This allows GEOMA to be more flexible or general purpose.

GEOMA and LGEOM are almost always used following a VICAR program that
generates the control information (tiepoints) in either a parameter
file or an input file.  This keeps the user from having to enter the
tiepoints manually.  This also makes most of the details regarding the
tiepoint grid transparent to the user.  GEOMA is often used following
VICAR program RESLOC and sometimes following VICAR program PICREG.

GEOMA's greater flexibility is not often needed.  GEOMA is generally 
slower than LGEOM, and it is less commonly used than LGEOM.  The
most common use of GEOMA is in conjunction with RESLOC.  In this
application GEOMA is used to remove geometric distortion in Voyager
pictures based on the reseau locations.   RESLOC writes the
control information in a file which is used as the second input file
for GEOMA.  Sometimes VICAR program FARENC is used between RESLOC 
and GEOMA to add a translation to the geometric transformation when
dealing with a limb of a planet or satelite.  GEOMA is used in
this application because the reseau do not lie in a rectangular grid.

In general, program LGEOM is recommended over GEOMA because of speed unless the
greater flexibility of GEOMA is needed.  For instance, it is more common to use
LGEOM than GEOMA following PICREG.  Where GEOMA is used, it is common to use a
grid of degenerate quadrilaterals (triangles) as shown below in the diagrams in
the OPERATION section.

Users can also use VICAR program MGEOM as an alternative to GEOMA when the
output tiepoints lie in a rectangular grid.  The advantages of MGEOM are
listed in the help for MGEOM.

1.    RESLOC (IMG,TF) (R,G) 
      GEOMA  (IMG,G) IMGB (1,1,800,800)

In this example GEOMA is used to remove geometric distortion in a Voyager
image named IMG based on the reseau locations.   RESLOC writes the
control information in file G which is used as the second input file
for GEOMA.  GEOMA produces an 800 by 800 output image named IMGB.
2.   PICREG (J,J2) (T,G)
       ... numerous steps to interactively select tiepoints
       NHOR=10 NVER=10 
     GEOMA INP=J OUT=JB PARMS=G (1,1,800,800)

In this example GEOMA is used to register image J to image J2 using
tiepoints selected interactively in program PICREG.  The tiepoints
selected by the user are converted into a rectangular grid with 10+1
rows and columns using a first-order surface fit (TPFORM=3).  Each
rectangle is divided into two triangles to make the final grid of
degenerate quadrilaterals (triangles) for GEOMA.  The control
information is passed to GEOMA via parameter file G.  GEOMA produces
the registered image JB.
3.   GEOMA INP=T3 OUT=T4 SIZE=(1,1,20,20) NAH=1 NAV=1          TIEPOINT=(1,1,1,1  1,20,1,15,  20,1,10,1  20,1,10,1)

In this example GEOMA is used to perform an affine (first order)
transformation defined by three tiepoints.  (1,1) is mapped to (1,1).
(1,15) is mapped to (1,20).  (10,1) is mapped to (20,1).  The rest of
the mapping is defined by linear interpolation.
4.  The last example is the test procedure for GEOMA.  This is
    a complete example that could be run by the user and that 
    demonstrates uses of the possible parameters.

    gen tgen1 10 10 linc=5 sinc=10
    geoma inp=tgen1 out=tgen2 size=(1,1,20,20) nah=1 nav=1 tiepoint=(1,1,1,1,     1,20,1,10,20,1,10,1,20,20,10,10)
    list tgen2
    gen tgen3 10 20 'half sinc=-1 linc=2
    geoma inp=tgen3 out=tgen4 size=(1,1,20,20) 'half nah=1 nav=1     tiepoint=(1,1,1,1,1,20,1,20,20,1,10,1,20,20,10,20)
    list tgen4
    gen tgen1 10 10 ival=10 sinc=2 linc=5
    geoma inp=tgen1 out=tgen6 size=(1,1,100,100) nah=1 nav=1 tiepoint=(1,1,1,1,     1,100,1,10,100,1,10,1,100,100,10,10) 'nointerp
    list tgen6
    gen tgen1 10 10
    geoma inp=tgen1 out=tgen7 size=(1,1,10,10) nah=1 nav=1 tiepoint=(1,1,10,10,     1,10,10,1,10,1,1,10,10,10,1,1) notify=20
    list tgen7
    gen tgen1 10 10 ival=10 sinc=2 linc=5
    geoma inp=(VGR:F1636832.FIC,VGR:F1636832.IGPR) out=tgen5
    label-list tgen5

GEOMA computes the DN value of each pixel in the output picture as follows.
First, decide within which quadrilateral area in the output picture the pixel
lies. Let the coordinates of the four control points defining that
quadrilateral in the output picture be called (x(j),y(j),j=1,2,3,4), and in the
input picture let them be called (x'(j), y'(j),j=1,2,3,4). (x may be considered
the line coordinate and y the sample, although they could as well be reversed
for this discussion.) The values of all 16 of these numbers are part of the
control parameters for GEOMA. Define a transformation by 
		x' = ax + by + cxy + d                                (1)
		y' = ex + fy + gxy + h
where the values of the coefficients a,b,..,h are determined by requiring that
		x'(j) = ax(j) + by(j) + cx(j)y(j) + d                 (2)
		y'(j) = ex(j) + fy(j) + gx(j)y(j) + h
		for j = 1,2,3,4

This is the condition that the defined transformation exactly maps
the control point coordinates as specified by the parameters.

The transformation is used to transform the line-sample coordinates of
the output picture pixel being processed into its corresponding coordinates in
the input picture. In general, although the output coordinates are always
integer, the input coordinates are not, so the input coordinates do not
correspond to any one pixel but instead fall in between the input pixels. If
intensity interpolation is NOT to be performed, as specified by the
NOINTERP parameter, then the DN value assigned to the output pixel
is that of the pixel in the input picture nearest the transformed coordinates.
The use of this mode considerably reduces the CPU time required by GEOMA. The
quality of the result is in general degraded but still may be quite
acceptable, depending on the nature of the transformation.

Assuming interpolation is to be performed to compute the DN value of
the output picture, the four pixels nearest to the transformed coordinates in
the input picture are selected.

The interpolation formula is
		DN = px' + qy' + rx'y' + s 			(3)
where x' and y' are the transformed coordinates of the output picture pixel,
DN is the data number to be assigned to the output pixel, and p,q,r and s are
constants chosen such that
		DN(k) = px'(k) + qy'(k) + rx'(k)y'(k) + s       (4)
		for k = 1,2,3,4
where x'(k), y'(k) and DN(k), k = 1,2,3,4 are the coordinates and DN values of
the four selected input pixels nearest to the transformed coordinates.

In order for GEOMA to work properly, the tiepoint coordinates in the
output picture must satisfy certain constraints. They must be organizable into
a rectangular matrix of points. That is, each tiepoint must belong to one row
and one column of points, each row must have the same number of points as every
other row, and each column must have the same number as every other column. It
should be possible to number the rows sequentially from top to bottom, and the
columns form left to right. Quadrilaterals used by the program for performing
the geometric transformation are formed by connecting each tiepoint to the ad-
jacent tiepoints in the same column and in the same row. The quadrilaterals
should all have the property that they are not concave (no interior angle of
any quadrilateral shoud be greater than 180 degrees). In addition, each tie-
point should be below any tiepoint in its column and in any preceding row, and
above any tiepoint in its column in any subsequent row. Similarly, each tie-
point should be to the right of any tiepoint in its row and in a preceding
column, and to the left of any tiepoint in its row and in a subsequent column.
The number of areas vertically is the number of rows of tiepoints less 1, and
the number of areas horizontally is the number of columns of tiepoints less 1.
The following caution is issued to potential users of GEOMA. Although the
transformation in equation (1) is very well behaved within the quadrilateral
in which it applies, there is no assurance that it is continuous across the
boundary between adjacent quadrilaterals unless the boundary is precisely
vertical or horizontal. The degree of discontinuity depends on the details
of the tiepoint specifications. In some cases the discontinuity may be so
small that the transformation picture has no visible defect, but in others
it may be quite visible.

The user has two choices to minimize the discontinuity. He may
carefully choose the transformation such that discontinuities are minimized,
or he may specify the tiepoints in such a way that the quadrilaterals
degenerate into triangles. A quadrilateral degenerates into a triangle if two
adjacent tiepoints defining it have the same coordinates. By suitable choice of
tiepoints, the user can specify the coordinates in such a way that all
quadrilaterals degenerate into triangles.  This, in fact, is commonly done.
(See programs RESLOC and PICREG.)  The figures below show sample arrays
of tiepoints satisfying this condition in output-image views.  (The rows
and columns are straight below because it was easier to draw, but they are not 
required to be straight.)
           1,2       3,4        5
            +---------+---------             | \       | \       |
            |   \     |   \     |
            |     \   |     \   |
            |       \7|8      \ |
          6 +---------+---------+ 9,10
            |       / |       / |
            |     /   |     /   |
            |   /     |   /     |
            | /       | /       |
            +---------+---------           11,12     13,14       15

   Sample Triangular Tiepoint Array: 4 areas horizontally, 2 areas vertically.
           1,2       3,4        5,6      7,8       
            +---------+---------+---------             | \       | \       | \       | \
            |   \     |   \     |   \     |   \
            |     \   |     \   |     \   |     \
            |       \ |10     \ |12     \ |14     \
          9 +---------+---------+---------+---------+ 16
            |       / |11     / |13     / |15     /
            |     /   |     /   |     /   |     /  
            |   /     |   /     |   /     |   /    
            | /       | /       | /       | /      
            +---------+---------+---------           17,18     19,20     21,22     23,24

   Sample Triangular Tiepoint Array: 7 areas horizontally, 2 areas vertically.
When an area is a triangle, the cross terms in equation (1) are
eliminated and the transformation becomes
		x' = ax + by + d				(5)
		y' = ex + fy + h
where the values of the coefficients are determined by requiring that
		x'(j) = ax(j) + by(j) + d			(6)
		y'(j) = ex(j) + fy(j) + h
		for j = 1,2,3,4
The values of j designate the three distinct tiepoints defining the triangle.
Because equation (5) is linear, the transformation is guaranteed to be
continuous at the triangle boundaries and discontinuities in the output
picture cannot occur.

There is a question about how to transform picture samples which fall
outside all the quadrilaterals defined by parameters. One answer is to avoid
this situation by specifying the control point locations so that every sample
in the output picture falls within some quadrilateral. This requires that some
control points fall exactly on or outside the border of the output picture.
The program is not affected in any way by using control points outside the

If some picture samples do fall outside the defined quadrilaterals, the
program will process each such sample by assigning it to a nearby quadrila-
teral. This is done by "extending" the boundaries of all "edge" quadrilaterals
to the picture borders. The details of the precise extrapolation algorithm
will not be given. The user should understand that the algorithm can lead to
discontinuities in the processed picture in the regions outside the defined
quadrilaterals, even when the control areas are triangles.
1. The input and output images must be byte or halfword data.
2. There are essentially no restrictions in GEOMA on picture size or 
   number of tiepoints.  (The input image may not have >32767 lines.)
3. The maximum number of tiepoints that can be passed via a 
   parameter file is currently 2047.  
4. The maximum number of tiepoints that can be specified as 
   part of the TAE command line is 125.
5. The maximum number of tiepoints that can be passed via an IBIS tiepoint
   file is set to 10,000.  (This could be increased, but you would probably
   want to increase the size of BIGBUF proportionally.)

Under certain circumstances, such as 90 degree rotations, GEOMA will be much
slower than LGEOM.  (An "LGEOMA" could be written, using the internal algorithm
of LGEOM, which would be much faster in those circumstances. However, such a
program does not currently exist.)   For small rotations GEOMA will be faster
than LGEOM.  As with most VICAR programs, the elapsed time for execution will
be shortened if the input and output files are on separate disk drives. 

Under VMS the execution time and the number of page faults can be affected by
the working set size that was set by the system manager.  If a large
number of page faults are observed, the working set size should be increased.
Alternatively, it may be possible to reduce the NBBUF value in the program. 


This portable version of GEOMA will generally produce identical output on
all of the MIPS-supported platforms.  This output will frequently not be
identical to the output of the unported version of GEOMA, but the differences
will be trivial and should be improvements.

This section is required by the MSTP SRD for all programs that
do not meet the precision requirements given in the SRD.  The two requirements
not met by GEOMA are:
1) Integer data output by a ported program must be the same as that produced
   by the unported program.
2) Integer data output by a ported program must be the same on all of the
   MIPS-supported platforms.
All of the programs that do geometric transformations are going to have problems
meeting one or both of these requirements.   This includes MAP3, LGEOM, MGEOM,
GEOMA, and several others.  GEOMA was modified to increase the use of
Double Precision floating point arithmetic so that the first requirement is met
"virtually always" and the second requirement is met "well enough". 
The definition of the above qouted terms is given towards the end of this 

In the port, GEOMA was modified to retain the full precision returned from
subroutine DGELG.  This was done to eliminate or minimize the differences
produced on different platforms.  The results are somewhat more accurate
than those of the unported GEOMA, but since the output DN values are rounded
to the nearest integer, the output images are essentially the same as for the
unported GEOMA.

In computing the output pixel value (DN), GEOMA usually interpolates between
four input pixel values, but if the input line coordinate is close enough to an
integer value, it interpolates between two input pixel values.  For byte images,
the two pixel interpolation and the four pixel interpolation yield the
same result virtually always once the result is rounded to the nearest
integer.  For halfword images, where the pixel values are often a hundred times
greater, the two pixel interpolation and the four pixel interpolation
occassionally yield slightly different results.

The image f1636832.fic used in the GEOMA test pdf shows what type of differences
are likely to occur (with the ported GEOMA) between machines, and what type of
differences are likely to occur between the ported GEOMA and the unported GEOMA.
The output file from the ported GEOMA is identical on all of the MIPS supported
machines.  The output file from the ported GEOMA is very slightly different
from the output file from the ported GEOMA.  98% of the pixels are identical.
Almost all of the rest differ in value by 1.  The minimum difference is -10.
The maximum difference is 17.  Program DIFPIC shows the AVE VAL OF
DIFFS=-0.125000E-03, or approximately .0001. 

The cause of the differences is that the single precision computations
(unported GEOMA) round differently than the double precision computations
(ported GEOMA), resulting in differences in when two pixel interpolation is
used and when four pixel interpolation is used.  The above 
halfword image has a number of locations where the difference
between adjacent pixels is greater than 10,000.  This combines to produce
a slight but insignificant difference.  This difference is estimated at
less than 1% of the difference between using 'NOINTERP and using regular
interpolation, which is the definition I will use for "well enough".
I believe that such differences will not be humanly observable in an image

To complete this discussion, I want to define the other quoted term mentioned
above. "Virtually always" is always minus the probability of a difference in
output between different MIPS supported machines.  Differences are possible
1) one machine does a two pixel interpolation while another machine does
   a four pixel interpolation.
2) one machine produces a double-precision output pixel value very slightly
   less than halfway between two integers and another machine produces an
   output pixel value equal to or slightly greater than halfway between
   two integers.
3) one machine determines that an output pixel lies outside the geometrically
   transformed boundaries of the input image while another machine decides
   the output pixel lies within the boundaries.

Possibilty 2 is approximately the probability of a double precision
result in the range between 0.0 and 1.0 being close enough to 0.5 to round
differently on two machines.  Double-precision rounding differences start at
about 1.E-14 on the MIPS-supported machines and might grow to a maximum of
1.E-8 going through the computations in GEOMA.  So the above probability is 
less than the probability of the result being in the range 0.49999999 and
0.50000000.  Assuming a fairly uniform distribution, this is probability
is 1.E-8.

Possibility 1 can be shown to have a similar probability.  Possibility 3 has an
even lower probability because the border pixels are such a small fraction of
the total pixels.  However, possibily 3 could cause one machine to generate a
pixel value of 0 while another machine could generate a relatively large pixel 
value from the edge of the input image.  The other possibilities only generate
a difference of 1 (possibility 2) or something small relative to the difference
between adjacent pixels (possibility 1).  In an image display, none of these
differences should be humanly observable.   Based on this discussion, "virtually
always" is defined as no more than 1 different pixel in a set of 1.E8 pixels.
This amounts to no more than one difference in a hundred 1000 by 1000 pixel



Made portable for UNIX ... Jim Turner (CRI)      5 September 1994


  7-95    SP  Added support for an IBIS tiepoint file as a second input.
  6-96    OAM Deleted default value of parameter OUT in geoma.pdf to make 
              program compatible with SAGE (DFR).



Input file name and optional IBIS tiepoint file


Output file name.


Standard VICAR size field.


Number of lines in output.


Number of samples in output.


Input data format. Valid: BYTE,HALF.


Number of areas horizontally. See explanation.


Number of areas vertically. See explanation.


Specifies mapping of control points between output and input pictures.


Specifies that intensity interpolation is not to be perfomed after the mapping of each pixel.


Causes a message to be typed at the user's terminal after completion of each NOTIFY percent of output lines.

See Examples:

Cognizant Programmer: