
#define VERSION "1.0"                      /* don't touch this */
#define VDATE   "November 1, 1991"

/*
Name:
 vs2ris


---------------------------------------------------------------------------
Purpose:
---------------------------------------------------------------------------
 To produce an 8-bit Raster Image Set (RIS8) from two-dimensional HDF Vset
 data or "curvilinear" HDF data.  The HDF Vset file may contain both
 triangular and quadrangular elements.  Data defined on a curvilinear grid
 may also be supplied in the straight HDF file format.  In this case input
 consists of three common HDF input files, each containing a single 
 two-dimensional array (SDS) describing the gridpoint data values, 
 x coordinates, and y coordinates respectively.

                                   -----------
       floating point data        |           |
       (single HDF Vset file ---> |  vs2ris   | ------> RIS8
        or three HDF files)       |           |
                                   -----------


---------------------------------------------------------------------------
Source Availability:
---------------------------------------------------------------------------
 This software was developed and made available by:

    Fred Walsteijn
    Institute for Marine and Atmospheric research
    University of Utrecht
    Princetonplein 5
    3584 CC Utrecht
    The Netherlands

    E-mail:  walsteyn@fys.ruu.nl

 You are free to use this software, however all commercial rights are reserved
 by Fred Walsteijn.  This software may not be sold or distributed for profit,
 or included with other software which is sold or distributed for profit
 without the prior permission of the author Fred Walsteijn.

 HDF and HDF Vset file format and interface developed at the National
 Center for Supercomputing Applications at the University of Illinois at
 Urbana-Champaign.

 FRED WALSTEIJN AND THE UNIVERSITIES OF ILLINOIS AND UTRECHT GIVE NO WARRANTY,
 EXPRESSED OR IMPLIED, FOR THE SOFTWARE AND/OR DOCUMENTATION PROVIDED,
 INCLUDING, WITHOUT LIMITATION, WARRANTY OF MERCHANTABILITY AND WARRANTY OF
 FITNESS FOR A PARTICULAR PURPOSE.


---------------------------------------------------------------------------
Syntax:
---------------------------------------------------------------------------
Read the source of function usage(), or compile this tool and execute it
without arguments.


---------------------------------------------------------------------------
Implementation Notes -1- :
---------------------------------------------------------------------------
The "kernel" of this tool consists of functions generating pixel values
inside a triangle and quadrangle.  A few versions of each of these functions
are given below and partly to users of this tool as well via the command line
options -iso and -isot.  I will briefly discuss the pros and cons of the
supplied functions.
     None of these functions can handle "degenerate" polygons: triangles
with all vertices lying on a straight line, quadrangles with crossing sides.
However, I suppose that you preclude such elements from your calculations too,
so this restriction should not bother you.

 void triangle()
 ---------------
  Generates pixels inside a triangle by an extension of the commonly used
  "computer graphics" polygon-fill routine based on intersections of horizontal
  pixel lines and the triangle sides.
  Advantages:
  - the algorithm is fast.
  Disadvantages:
  - ???

 void iso_triangle()
 -------------------
  Generates pixels inside a triangle using an isoparametric mapping known
  from finite element methods.  If you specify -isot at the command line,
  then "vs2ris" uses this function instead of triangle().  The results should
  be (almost) identical to the raster generated by triangle().
  Advantages:
  - concise code, easy to understand.
  - easy to generalize to 3D.
  Disadvantages:
  - Less efficient than triangle() due to the simple (initial) pixel selection
    strategy.  This strategy becomes very slow for slim triangles oriented at
    an angle of 45 degrees with respect to pixel lines.
  
 void quadrangle()
 -----------------
  Generates pixels inside a quadrangle by decomposing it into two triangles.
  Again usage of the command line switch -isot determines whether "vs2ris"
  uses triangle() or iso_triangle() for these two triangles.  Note
  that the decomposition can be made in two ways if and only if the
  quadrangle is convex.  No special strategy is implemented yet to make the
  best choice in this case.  However, if the quadrangle is not convex, only
  one decomposition is allowed, which of course is determined by function
  quadrangle().
  Advantages:
  - is as robust/fast as the "triangle-function" it uses.  Can handle
    non-convex quadrangles.
  Disadvantages:
  - if the number of pixels used for a single quadrangle is large, the
    resulting raster may show a "kink".  In this case you may conclude that
    the technique is not a very natural one (e.g.: if your data is generated
    by a bilinear finite element method).  However, if there is only a
    small number of pixels in the quadrangle (or your data is smooth) this
    will not be noticable.


 void iso_quadrangle()
 ---------------------
  Generates pixels inside a quadrangle using a bilinear isoparametric
  mapping.  Can be considered as the generalization of the NCSA tool "fp2hdf"
  to non-rectangular data.  It yields results identical to "fp2hdf" if
  applied to rectangular data.
  This function is used instead of quadrangle() if you specify
  the -iso command line option.
  Advantages:
  - gives best results, even if the number of vertices is small.
    Small-scale details and saddle points are treated symmetrically.
    Also, this is the natural technique if your data is generated by a
    bilinear finite element technique.  However the difference with
    function quadrangle() is noticable only if you are using a "large" number
    of pixels for a single quadrangle.
  Disadvantages:
  - cannot handle non-convex quadrangles (a limitation of isoparametric
    mappings in general).  The current implementation recognizes this case
    automatically and will locally use function quadrangle() instead.  Only
    in this case is the usage of the command line flag -isot relevant.
  - highly distorted (but still convex) quadrangles lead to
    somewhat larger CPU times, since more Newton iterations are needed.

 void simple_iso_quadrangle()
 ----------------------------
  Simpler version of iso_quadrangle().  Additional pros & cons:
  Advantages:
  - concise code.
  - easy to generalize to 3D.
  Disadvantages:
  - this simple implementation does the final interior/exterior pixel
    classification in isoparametric coordinates.  This may lead to failure
    of the algorithm in the case of a highly distorted but still convex
    quadrangle, since the Jacobian of the mapping may become singular just
    outside the quadrangle.
  - this version becomes very slow for slim quadrangles oriented at
    an angle of 45 degrees with respect to pixel lines.

---------------------------------------------------------------------------
Implementation Notes -2- :
---------------------------------------------------------------------------
The current "read_vset" can handle both LOCAL_INTTYPE and LOCAL_LONGTYPE
connectivity lists.  For small Vsets one can still use LOCAL_INTTYPE lists;
to read pre version 2.0 Vsets one certainly needs this capability, whereas
for large files LOCAL_LONGTYPE storage is mandatory.  Nowadays simulations
with more than 32767 gridpoints/vertices are quite common.  Such data sets
cannot be stored in a HDF Vset file without the use of LOCAL_LONGTYPE
connectivity lists.

I have used the undocumented but PUBLIC function VSsizeof() to distinguish
LOCAL_INTTYPE from LOCAL_LONGTYPE lists.  Without using this function only
the total vdata element size returned in arg vsize by VSinquire() is
available.  This however includes the size of other fields if present.
Thus the list element size is obtained only if the list were the only field
present in its vdata.  I consider this too restrictive, and therefore have
used VSsizeof() instead, and urge the NCSA developers to include this
function in the documentation.
 

---------------------------------------------------------------------------
Speed comparisons:
---------------------------------------------------------------------------
Speed test on a MacIIsi with coprocessor and built-in colors disabled.
Rectangular input data test cases are included to allow comparison of
"vs2ris" with the standard NCSA HDF tool "fp2hdf".
Timings accurate to one second.  Tools compiled with MPW C 3.2 using
highest optimization ("-opt full").

   -------------------------------------------------------------------------------------
   tool-          SDS1=40*40     CRV1=40*40     CRV2=40*40    SDS2=257*257  SDS2=257*257
   name            RIS=300*300    RIS=424*424    RIS=390*390   RIS=300*300   RIS=600*600
   -------------------------------------------------------------------------------------
   fp2hdf -i         11 secs        --             --            12 secs       44 secs
   vs2ris            16 secs        17 secs        16 secs       84 secs      130 secs  # default
   vs2ris -isot      26 secs        26 secs        34 secs       79 secs      145 secs  # result equiv to default
   vs2ris -iso       61 secs        60 secs        75 secs      104 secs      271 secs  # best results
   -------------------------------------------------------------------------------------

SDS1:  
  rectangular input data. For "vs2ris" two additional (trivial) grid
  coordinate HDF files had to be constructed (the first contained an
  SDS x(i,j)==i, the second an SDS y(i,j)==j).
CRV1: 
  curvilinear data obtained by rotating SDS1 45 degrees (by x(i,j)==i-j,
  and y(i,j)==i+j). As a result the quadrangle sides are no longer aligned
  with pixel lines, and the standard NCSA tools are no longer applicable.
  Execution speeds are almost identical to the SDS case.
CRV2:
  again SDS1, but now on a polar grid. Due to the presence of slim
  gridcells (// 45 degrees with pixel lines) the performance of -isot
  degrades. The -iso option is slower because 3 Newton iterations
  (instead of 2) are needed in approx 60% of the pixels.
SDS2:
  similar to SDS1, but with larger array dimensions. The "vs2ris" overhead
  per gridcell now becomes very noticable.
 (In this case enabling built-in color may increase CPU times of all tools
  by almost 50%.  In the other cases enabling color increased CPU times by 10%.
  No speed penalty is involved on MacII's with a color video card.)

Conclusion:
  The default setting of "vs2ris" is very efficient for small SDSs, and may
  be acceptable for large SDSs. The performance of "vs2ris -iso" is acceptable.

Notes:
- CRV1 & CRV2 image dimensions have been chosen such that the number of
  pixels inside the domain is approx equal to the SDS1 case.  It is therefore
  legal to compare execution times.
- "vs2ris -iso" is a generalization of "fp2hdf" to non-rectangular data, and
  in certain cases yields much better results than the default setting.
- vs2ris has been tested with HDF Vset input also, by storing the curvilinear
  data in a single HDF Vset file. The "vs2ris -vs ..." timings are similar
  to the above.

---------------------------------------------------------------------------
"To do":
---------------------------------------------------------------------------
1. Function quadrangle(): improve function diagonal_choice such that it
   yields an optimal choice if the quadrangle is convex.  Remember, for
   convex quadrangles both diagonals can be used, however if the quadrangle
   is "almost" non-convex one of the two may be preferable.  The optimal
   choice can be found by "balancing/comparison of the sub-triangle areas".
2. HDF Vset routine VSinquire() expects a pointer to int to store the number
   of elements in.  This should be modified into a pointer to int32. Otherwise
   no large Vset files can be read on machines with 16-bit ints.  This one is
   for NCSA.


---------------------------------------------------------------------------
Limitations & Bugs & ...:
---------------------------------------------------------------------------
1. Triangles and quadrangles should not be degenerate. Is hard to fix without
   losing pixels in other cases.
2. The sequence of quadrangle vertex-indices as given in each of the
   (HDF Vset) connectivity list fields should correspond to a monotonic walk
   around the quadrangle perimeter. No "butterfly" check is made.
3. On machines with 16-bit ints it does not seem possible to read large Vset
   files. This is a bug in the HDF Vset lib (?), see "To Do".

If you encounter a bug, report it (and a fix) to:
   Fred Walsteijn
   walsteyn@fys.ruu.nl
 

---------------------------------------------------------------------------
History:
---------------------------------------------------------------------------
June  9, 1991 -- first beta version
Sept 15, 1991 -- last cosmetic changes
Oct  23, 1991 -- parse_cmd() now signals missing mandatory args
Nov   1, 1991 -- declaration of connectivity lists modified: now runs on Cray
             
      (by:  Fred Walsteijn,  walsteyn@fys.ruu.nl)

*/


/* ---------------------------- includes --------------------------- */

#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <math.h>
#include <string.h>
#include "vg.h"

/* ---------------------------- macros --------------------------- */

#undef malloc
#undef free      
       /* ^ the NCSA def of free(p) for the Mac fails if p==NULL */

/* The rounding macros below use a temporary var, which is more
 * efficient since they are invoked with complicated args.
 * Functions using rounding macros must define:  float tmp;
 */

/*  #define ROUND( x )       ( (x)<0?(int32)(x-0.5):(int32)(x+0.5) )   */
#define ROUNDUP( x )     ( tmp=(x), tmp<0? (int32)tmp   :(int32)tmp+1 )
#define ROUNDDOWN( x )   ( tmp=(x), tmp<0? (int32)tmp-1 :(int32)tmp   )

/* The min/max macros are only invoked with simple args: */

#define MIN(a,b)         ( (a)<(b)?(a):(b) )
#define MAX(a,b)         ( (a)>(b)?(a):(b) )
#define MIN3(a,b,c)      ( (a)<(b)?   MIN(a,c)  :  MIN(b,c) )
#define MAX3(a,b,c)      ( (a)>(b)?   MAX(a,c)  :  MAX(b,c) )
#define MIN4(a,b,c,d)    ( (a)<(b)?                                   \
                                   ((a)<(c)?  MIN(a,d)  :  MIN(c,d) ) \
                                  :                                   \
                                   ((b)<(c)?  MIN(b,d)  :  MIN(c,d) ) )
#define MAX4(a,b,c,d)    ( (a)>(b)?                                   \
                                   ((a)>(c)?  MAX(a,d)  :  MAX(c,d) ) \
                                  :                                   \
                                   ((b)>(c)?  MAX(b,d)  :  MAX(c,d) ) )


#define IMAGE(i,j)       *(im->image + (im->vres-1-(j))*im->hres + (i))
        /*
         * In "vs2ris" the indices i,j are always proportional to x,y; 
         * so IMAGE(0,0) corresponds to the pixel in the lower/left corner
         * of the raster, which is schematically:  raster[vres-1][0].
         *    Equivalence:  IMAGE(i,j) == raster[vres-1-j][i].
         *    Dimensions:  IMAGE(hres,vres);  raster[vres][hres]
         * im->image points to the first element of the raster.
         */

/* ---------------------------- types, consts --------------------------- */

#define PIXEL      unsigned char   /* image pixel type */
#define PAL        unsigned char   /* palette entry type */
#define BACKGR_COL 0               /* background color in image */
#define FNAME_MAX  FILENAME_MAX    /* maximum number of chars in filenames */
#define STR_MAX    64              /* maximum number of chars in strings */
#define CMAX       10              /* maximum number of connectivity lists to look for */
#define MIN_RES    3               /* minimum number of pixels in x & y direction */

#define LOC_INT    int             /* the LOCAL_INTTYPE */
#define LOC_LONG   long            /* the LOCAL_LONGTYPE */

enum boolean { FALSE, TRUE };
enum filetype { CRV, VSET };         /* input file(s) type: CRV == 3*HDF, VSET == 1 HDF Vset */
enum elmtype { INTTYPE, LONGTYPE };  /* integral type of HDF Vset connectivity lists:
                                        INTTYPE == LOCAL_INTTYPE, LONGTYPE == LOCAL_LONGTYPE */

/*
 * structure definition for the input data
 */
struct Input {
  enum filetype ftyp;       /* ftyp == CRV then curvilinear grid/data (from three HDF files)
                               ftyp == VSET then vertex set input (from single HDF Vset file) */
  int32 dim0, dim1;         /* curvilinear input dimensions - ncols, nrows; CRV only */
  float32 *x;               /* horizontal coords of vertices */
  float32 *y;               /* vertical coords of vertices */
  float32 *data;            /* vertex data values */
                              /* CRV:  x[dim0][dim1], y[dim0][dim1], data[dim0][dim1]
                               * VSET: x[nv],         y[nv],         data[nv]
                               */

  int32 nv,                 /* number of vertices, VSET only */
        ntri, nquad;        /* number of triangles, quadrangles in (s)plist3, (s)plist4 */
  enum elmtype etyp3,       /* etyp3 == INTTYPE then triangles in splist3; otherwise plist3 */
               etyp4;       /* etyp4 == INTTYPE then quadrangles in splist4; otherwise plist4 */
  LOC_LONG *plist3, *plist4;   /* pointers to LOCAL_LONGTYPE connectivity lists */
  LOC_INT *splist3, *splist4;  /* pointers to LOCAL_INTTYPE connectivity lists */
};


/*
 * structure definition for the output raster image
 */
struct Raster {
  int32 hres;       /* horizontal resolution of the image */
  int32 vres;       /* vertical resolution of the image */
  PIXEL *image;
};


/*
 * structure definition for command line options
 */
struct Options {
  char infile[FNAME_MAX], 
       gridx[FNAME_MAX], 
       gridy[FNAME_MAX];      /* input HDF/HDF Vset file names */
  char outfile[FNAME_MAX];    /* output HDF RIS8 file name */
  char palfile[FNAME_MAX];    /* palette file name */
  enum boolean pal,           /* load palette file */
               domain,        /* spatial (x,y-) domain is specified */
               minmax,        /* datamin and datamax specified */
               connect;       /* connectivity list(s) was specified by user */
  int ncon;                   /* number of connectivity lists specified, VSET only */
  char groupname[STR_MAX],  
       pxname[STR_MAX],
       pyname[STR_MAX],
       connectname[CMAX][STR_MAX],
       scalarname[STR_MAX];   /* vgroup and vdata names to use, VSET only */
};

/* ---------------------------- globals --------------------------- */

enum boolean iso, isot;   /* iso  == TRUE means "use bilinear quadrangles" */
                          /* isot == TRUE means "use iso_triangle() instead of triangle()" */
float32 datamin, datamax, /* data range, may scale or clip */
        xmin, xmax, 
        ymin, ymax;       /* selected domain, may clip true domain */


/* ---------------------------- prototypes --------------------------- */

void quit( char *msg );
void clear( struct Input *in, struct Options *opt, struct Raster *im );
void atexit_free( struct Input *in, struct Raster *im );
void free_mem( void );
void usage( char *tool );
void parse_cmd( int argc, char *argv[], 
                struct Input *in, struct Raster *im, struct Options *opt );
void write_opt( struct Options opt );
void write_in( struct Input in );
void set_pal( char *palfile );
void read_crv( struct Input *in, struct Options *opt );
void read_vset( struct Input *in, struct Options *opt ); 
void gminmax( float32 *xmin, float32 *xmax, float32 *x, int32 dim );
void setup_image( struct Raster *im );
void write_image( char *file, struct Raster im );
void curvi2ris( struct Input in, struct Raster *im );
void vset2ris( struct Input in, struct Raster *im );
void triangle( float32 x[3], float32 y[3], float32 dat[3], struct Raster *im );
void iso_triangle( float32 x[3], float32 y[3], float32 dat[3], struct Raster *im );
void quadrangle( float32 x[4], float32 y[4], float32 dat[4], struct Raster *im );
int diagonal_choice( float32 x[4], float32 y[4] );
int convex( float32 x[4], float32 y[4] );
void simple_iso_quadrangle( float32 x[4], float32 y[4], float32 dat[4], struct Raster *im );
void write_iter( void );
void iso_quadrangle( float32 x[4], float32 y[4], float32 dat[4], struct Raster *im );

/* ---------------------------- HDF 3.10 prototypes --------------------------- */

/*  DF *DFopen( char *name, int access, int ndds ); ---  already oldstyle prototyped in df.h !!! */
int DFclose( DF *dfile );
int DFPgetpal( char *filename, char *palette );
int DFR8setpalette( char *pal );
int DFR8putimage( char *filename, char *image, int32 xdim, int32 ydim, int compress );
int DFSDgetdims( char *filename, int *prank, int32 sizes[], int maxrank );
int DFSDgetmaxmin( float32 *pmax, float32 *pmin );
int DFSDgetdata( char *filename, int rank, int32 maxsizes[], float32 data[] );


/* ---------------------------- HDF Vset 2.0 prototypes --------------------------- */

VGROUP *Vattach( DF *f, int vgid, char *accesstype );  
void Vdetach( VGROUP *vg );
void Vgetname( VGROUP *vg, char *vgname );
int Visvs( VGROUP *vg, int id );
int Vgetid( DF* f, int vgid );
int Vgetnext( VGROUP *vg, int id );
VDATA *VSattach( DF *f, int vsid, char *accesstype );
void VSdetach( VDATA *vs );
int VSread( VDATA *vs, unsigned char *buf, int32 nelt, int interlace );
int VSfexist( VDATA *vs, char *fields );
int VSinquire( VDATA *vs, int *nelt, int *interlace, char *fields, int *eltsize, char *vsname);
int VSsetfields( VDATA *vs, char *fields );
int VSsizeof( VDATA *vs, char *fields );



/*********************************************************************************
 * main: Read "curvilinear" HDF input files, or single HDF Vset file.            *
 *       Calculate raster by interpolation in each triangle and/or quadrangle.   *
 *       Write raster to HDF RIS8 file.                                          *
 *********************************************************************************/

main( int argc, char *argv[] )
{ 
  static struct Input  in;       /* curvilinear or vset input data */
  static struct Raster im;       /* raster image to be generated */
  struct Options opt;            /* command line options */

  clear( &in, &opt, &im );       /* set all fields to 0, NULL, ... */
  atexit_free( &in, &im );       /* remember to free malloc-ed stuff */
  
  parse_cmd( argc, argv, &in, &im, &opt );  
                                 /* get switches, filenames, and output dims */
#ifdef DEBUGGING
  write_opt( opt );
#endif

  setup_image( &im );            /* initialise image storage space. */

  if ( opt.pal == TRUE ) 
    set_pal( opt.palfile );      /* read HDF palette and make it default. */

  if ( in.ftyp == CRV )
  { read_crv( &in, &opt );       /* read 3 HDF files */
    curvi2ris( in, &im );        /* generate the raster from curvilinear data. */
  }
  else if( in.ftyp == VSET )
  { read_vset( &in, &opt );      /* read single HDF Vset file */
    vset2ris( in, &im );         /* generate the raster from vertex set. */
  }
  else
    quit("main: Unknown 'in.ftyp' value.");

  write_image( opt.outfile, im );

#ifdef DEBUGGING
  write_in( in );
  write_iter();
#endif
  return 0;
}

/****************************************************
 *  quit:  Write message and terminate program.     *
 ****************************************************/

void quit( char *msg )
{
  printf("%s\n", msg );
  exit(EXIT_FAILURE);
}

/****************************************************
 *  clear:  Sets all struct fields to 'zero'.       *
 ****************************************************/

void clear( struct Input *in, struct Options *opt, struct Raster *im )
{
  in->ftyp = CRV;
  in->x = in->y = in->data  = NULL;
  in->plist3  = in->plist4  = NULL;
  in->splist3 = in->splist4 = NULL;
  in->etyp3 = in->etyp4 = INTTYPE;
  in->dim0 = in->dim1 = in->ntri = in->nquad = in->nv = 0;

  (void) strcpy( opt->infile,  "" );  /* just in case these are printed (DEBUGGING) but */
  (void) strcpy( opt->gridx,   "" );  /* are not applicable ... */
  (void) strcpy( opt->gridy,   "" );
  (void) strcpy( opt->outfile, "" );
  (void) strcpy( opt->palfile, "" );

  im->image = NULL;
}

/****************************************************************************
 *  atexit_free:  Store pointers to structs which will contain malloc-ed    *
 *                data.  By an atexit() call these data will automatically  *
 *                be free-ed on exit.                                       *
 ****************************************************************************/

   /* 
    * Globals for atexit_free() and free_mem():
    */

   struct Input  *ptr_in;
   struct Raster *ptr_im;

void atexit_free( struct Input *in, struct Raster *im )
{
  extern struct Input  *ptr_in;
  extern struct Raster *ptr_im;

  ptr_in = in;
  ptr_im = im;
  if ( 0 != atexit( free_mem ) )
    quit("atexit_free: atexit() failed.");
}

/**************************************************************************
 *  free_mem:  Free malloc-ed stuff.                                      *
 *             ASSUMES: every pointer has been malloc-ed, or equals NULL. *
 **************************************************************************/

void free_mem(void)
{
  extern struct Input  *ptr_in;
  extern struct Raster *ptr_im;

  free( (void *) ptr_in->data );
  free( (void *) ptr_in->y );
  free( (void *) ptr_in->x );

  free( (void *) ptr_in->splist3 );
  free( (void *) ptr_in->plist3 );
  free( (void *) ptr_in->splist4 );
  free( (void *) ptr_in->plist4 );

  free( (void *) ptr_im->image );
}


/***************************************************************
 *  usage:  Give command line syntax and terminate program.    *
 ***************************************************************/

void usage( char *tool )
{
  printf( "\n\n");
  printf( "This is vs2ris (version: %s, dated: %s).\n", VERSION, VDATE );
  printf( "Generates an HDF RIS8 raster from two-dimensional HDF Vset data\n");
  printf( "or 'curvilinear' HDF data.\n");
  printf( "By:  Fred Walsteijn,  University of Utrecht,  The Netherlands\n");
  printf( "                      (walsteyn@fys.ruu.nl)\n\n");

  printf( "SYNTAX:\n" );
  printf( "  %s  -vs <infile>  \n", tool );
  printf( "          -o <outfile>  -res <hres> <vres>  [OPTIONS ...]  [VS-OPTIONS ...]\n\n" );

  printf( "OR:\n" );
  printf( "  %s  -crv <datafile> <gridxfile> <gridyfile>\n", tool );
  printf( "          -o <outfile>  -res <hres> <vres>  [OPTIONS ...]\n\n" );

  printf( "where:\n" );
  printf( "    -vs <infile>\n" );
  printf( "           Input HDF Vset file containing two-dimensional vertex data.\n" );
  printf( "    -crv <datafile> <gridxfile> <gridyfile>\n" );
  printf( "           Three HDF input files specifying two-dimensional curvilinear data.\n" );
  printf( "           Each file contains an SDS describing the data values, x coordinates,\n" );
  printf( "           and y coordinates of the gridpoints.\n" );
  printf( "    -o <outfile>\n" );
  printf( "           HDF RIS8 output file containing a raster image.\n" );
  printf( "    -res <hres> <vres>\n" );
  printf( "           Resolution of output file. Determines the raster aspect ratio.\n\n" );

  printf( "OPTIONS:\n" );
  printf( "    -iso   Use an iso-parametric bilinear mapping to calculate the color of\n" );
  printf( "           pixels inside all convex quadrangles. Non-convex quadrangles are\n" );
  printf( "           subdivided into two triangles. Default: all quadrangles are\n" );
  printf( "           split into two triangles. (Note: -iso gives best results for\n" );
  printf( "           quadrangles but is slower than default.)\n" );
/*
 *printf( "    -isot  Use an iso-parametric (linear) mapping to calculate the color of\n" );
 *printf( "           pixels inside all triangles. This option should have no effect on\n" );
 *printf( "           the raster. Default: computer graphics 'intersection' approach.\n" );
 * (Printing this option is confusing.)
 */
  printf( "    -domain <xmin> <xmax> <ymin> <ymax>\n" );
  printf( "           Specifies the domain part to be visualized. Allows zooming in/out.\n" );
  printf( "           Values do not affect the output raster aspect ratio, therefore a\n" );
  printf( "           stretched image results if you specify a domain with aspect ratio\n" );
  printf( "           different from the raster.  Default: min/max of coordinates in \n" );
  printf( "           input file(s).\n" );
  printf( "    -range <datamin> <datamax>\n" );
  printf( "           The palette is scaled on this data range. Allows clipping.\n" );
  printf( "           Default: min/max of input file data.\n" );
  printf( "    -p <palfile>\n" );
  printf( "           The palette from this HDF file is written to the output raster file.\n" );
  printf( "           Default: no palette is written to the output file.\n\n" );

  printf( "VS-OPTIONS:\n" );
  printf( "    -group <groupname>\n" );
  printf( "           Name of vgroup in the Vset input file to read the vertex vdata(s)\n" );
  printf( "           from. Default: the first vgroup encountered in the Vset file.\n" );
  printf( "    -px <px_vdataname>\n" );
  printf( "           Name of vdata field containing x coords. Default: 'px'.\n" );
  printf( "    -py <py_vdataname>\n" );
  printf( "           Name of vdata field containing y coords. Default: 'py'.\n" );
  printf( "  { -connect <connect_vdataname> } ...\n" );
  printf( "           Name(s) of vdata field(s) containing a 'triangular' and/or\n" );
  printf( "           'quadrangular' connectivity list. Each field name must be\n" );
  printf( "           preceded by the flag -connect. From the specified connectivity\n" );
  printf( "           lists only the first two of different type encountered in the\n" );
  printf( "           Vset file are read: one list of triangles and/or a list of\n" );
  printf( "           quadrangles, which may contain triangles also (as in PolyView\n" );
  printf( "           the 4th vertex id should be zero). Other field names are skipped.\n" );
  printf( "           Default: -connect plist -connect plist3 -connect plist4\n" );
  printf( "    -scalar <scalar_vdataname>\n" );
  printf( "           Name of vdata field containing scalar values. Default: 'scalar'.\n\n" );
/*
 *printf( "    -triangulate | -quadrangulate\n" );
 * (Would be nice too.  This would enable input of random data, without the
 *  need to supply a connectivity list: vs2ris would construct one by itself.)
 */

  printf( "NOTE: all options may be specified in an arbitrary order.\n\n" );
  exit(EXIT_FAILURE);
}


/*********************************************************
 *  parse_cmd:   Parse the command line. Defines:        *
 *                opt:     ALL fields,                   *
 *                im:      hres & vres fields only,      *
 *                in:      ftyp field only,              *
 *                globals: datamin, datamax,             *
 *                         xmin, xmax, ymin, ymax,       *
 *                         iso, isot.                    *
 *********************************************************/

void parse_cmd( int argc, char *argv[], 
                struct Input *in, struct Raster *im, struct Options *opt )
{
  enum boolean  is_crv, is_vs, outp, res;  /* mandatory switches */
  int    i;                                /* argv-loop index */
  long   hres, vres;                       /* for sscanf() */
  double x0, x1, y0, y1, d0, d1;           /* for sscanf() */

/* 
 *  Option defaults:
 */
  iso = isot = opt->domain = opt->minmax = opt->pal = opt->connect = FALSE;

/*
 *  Default vgroup and vdata names similar (not identical) to PolyView defaults:
 */
  (void) strcpy( opt->groupname, "" );     /* "the 1st vgroup" */
  (void) strcpy( opt->pxname, "px" );
  (void) strcpy( opt->pyname, "py" );
  (void) strcpy( opt->scalarname, "scalar" );
  (void) strcpy( opt->connectname[0], "plist" );
  (void) strcpy( opt->connectname[1], "plist3" );
  (void) strcpy( opt->connectname[2], "plist4" );
  opt->ncon = 0;                           /* init for loop below; final assignment is done later */

/*
 *  Local var init:
 */
  is_crv = is_vs = outp = res = FALSE;


  if (argc<2) usage( argv[0] );

  for ( i=1; i<argc; i++ ) 
  {
    if ( !strcmp(argv[i],"-crv") &&  i<argc-3 )
    { is_crv = TRUE;
      (void) strcpy(opt->infile, argv[++i]);
      (void) strcpy(opt->gridx, argv[++i]);
      (void) strcpy(opt->gridy, argv[++i]);
    }
    else if ( !strcmp(argv[i],"-vs") &&  i<argc-1 )
    { is_vs = TRUE;
      (void) strcpy(opt->infile, argv[++i]);
    }
    else if ( !strcmp(argv[i],"-o") &&  i<argc-1 )
    { outp = TRUE;
      (void) strcpy(opt->outfile, argv[++i]);
    }
    else if ( !strcmp(argv[i],"-iso") )
      iso = TRUE;
    else if ( !strcmp(argv[i],"-isot") )
      isot = TRUE;
    else if ( !strcmp(argv[i],"-domain") &&  i<argc-4 )
    { opt->domain = TRUE;
      if ( sscanf(argv[++i], "%lf", &x0) != 1  ||
           sscanf(argv[++i], "%lf", &x1) != 1  ||
           sscanf(argv[++i], "%lf", &y0) != 1  ||
           sscanf(argv[++i], "%lf", &y1) != 1  ||
           x0 >= x1 || y0 >= y1                 )
        quit("parse_cmd: Error in domain spec.");
      xmin = x0;  xmax = x1;
      ymin = y0;  ymax = y1;
    }
    else if ( !strcmp(argv[i],"-range") &&  i<argc-2 )
    { opt->minmax = TRUE;
      if ( sscanf(argv[++i], "%lf", &d0) != 1  ||
           sscanf(argv[++i], "%lf", &d1) != 1  ||
           d0 >= d1                            )
        quit("parse_cmd: Error in range spec.");
      datamin = d0;
      datamax = d1;
    }
    else if ( !strcmp(argv[i],"-res") &&  i<argc-2 )
    { res = TRUE;
      if ( sscanf(argv[++i], "%ld", &hres) != 1  ||
           sscanf(argv[++i], "%ld", &vres) != 1  ||
           hres < MIN_RES || vres < MIN_RES  )
        quit("parse_cmd: Error in resolution spec.");
      im->hres = hres;
      im->vres = vres;
    }
    else if ( !strcmp(argv[i],"-p") &&  i<argc-1 )
    { opt->pal = TRUE;
      (void) strcpy(opt->palfile, argv[++i]);
    }
    else if (!strcmp(argv[i],"-group") && i<argc-1 )
      (void) strcpy(opt->groupname,argv[++i]);
    else if (!strcmp(argv[i],"-px") && i<argc-1 )
      (void) strcpy(opt->pxname,argv[++i]);
    else if (!strcmp(argv[i],"-py") && i<argc-1 )
      (void) strcpy(opt->pyname,argv[++i]);
    else if (!strcmp(argv[i],"-scalar") && i<argc-1 )
      (void) strcpy(opt->scalarname,argv[++i]);
    else if (!strcmp(argv[i],"-connect") && i<argc-1 && opt->ncon < CMAX )
    { opt->connect = TRUE;
      (void) strcpy(opt->connectname[opt->ncon++],argv[++i]);
    }
    else
    { printf("parse_cmd: Unknown flag(s) encountered.\n");
      usage( argv[0] );
    }
  }

/*
 *  Do final checks on the input:
 */

  if ( outp == FALSE )
    quit("parse_cmd: Output file is mandatory. Use option -o.\n");

  if ( res == FALSE )
    quit("parse_cmd: Output raster resolution is mandatory. Use option -res.\n");


  if ( opt->connect == FALSE ) 
    opt->ncon = 3;                      /* use the three default connectivity lists */
  else if ( opt->ncon > 2 )
    printf("parse_cmd: Warning: at most 2 connectivity lists will be used.\n");


  if ( is_crv == TRUE )                 /* Three straight HDF files specified. */
  { if ( is_vs == TRUE )
      quit("parse_cmd: Use only one of -vs, -crv.\n");
    in->ftyp = CRV;
  }
  else if ( is_vs == TRUE )             /* Single HDF Vset file specified. */
    in->ftyp = VSET;
  else
    quit("parse_cmd: Input file(s) are mandatory. Use one of -vs, -crv.\n");
}





#ifdef DEBUGGING
       /* ^ needed since DEC Ultrix C cannot handle the macros below */

/************************************************************
 * write_opt:  Write contents of 'opt' to stdout.           *
 *             For debugging.                               *
 ************************************************************/

#define pr_bool( arg )   printf(#arg " = %s\n", (arg==FALSE)?"FALSE":((arg==TRUE)?"TRUE":"*TRASH*") )
#define pr_float( arg )  printf(#arg " = %lf\n", (double)arg )
#define pr_int( arg )    printf(#arg " = %ld\n", (long)arg )
#define pr_string( arg ) printf(#arg " = %s\n", arg )
#define pr_ftyp( arg )   printf(#arg " = %s\n", (arg==CRV)?"CRV":((arg==VSET)?"VSET":"*TRASH*") )
#define pr_etyp( arg )   printf(#arg " = %s\n", (arg==INTTYPE)?"INTTYPE":((arg==LONGTYPE)?"LONGTYPE":"*TRASH*") )

void write_opt( struct Options opt )
{
  int i;

  printf("\n\n");
  printf("write_opt: --- contents of 'struct Options opt' ---\n");

  pr_string( opt.infile );
  pr_string( opt.gridx );
  pr_string( opt.gridy );
  pr_string( opt.outfile );
  pr_string( opt.palfile );

  pr_string( opt.groupname );
  pr_string( opt.pxname );
  pr_string( opt.pyname );
  pr_string( opt.scalarname );
  pr_int( opt.ncon );
  for (i=0; i<opt.ncon; i++)
    printf("opt.connectname[%d] = %s\n", i, opt.connectname[i] );

  pr_bool( opt.pal );
  pr_bool( opt.domain );
  pr_bool( opt.minmax );
  pr_bool( opt.connect );

  printf("\n\n");
}

/******************************************************
 * write_in:   Write contents of 'in' to stdout.      *
 *             For debugging.                         *
 ******************************************************/

void write_in( struct Input in )
{
  LOC_LONG i, p1, p2, p3, p4;

  printf("\n\n");
  printf("write_in: --- contents of 'struct Input in' ---\n");

  pr_ftyp( in.ftyp );

  pr_int( in.dim0 );
  pr_int( in.dim1 );

  pr_int( in.nv );
  pr_int( in.ntri );
  pr_int( in.nquad );

  pr_etyp( in.etyp3 );
  pr_etyp( in.etyp4 );

/*  if ( in.ftyp==VSET && in.etyp4==LONGTYPE )
 *  {
 *    printf("ARRAY --- LOC_LONG *in.plist4:\n");
 *    for ( i=0; i<in.nquad; i++ )
 *    { p1 = *in.plist4++; 
 *      p2 = *in.plist4++; 
 *      p3 = *in.plist4++; 
 *      p4 = *in.plist4++; 
 *      printf("%ld,%ld,%ld,%ld\n", (long)p1, (long)p2, (long)p3, (long)p4 );
 *    }
 *  }
 */

  printf("\n\n");
}

#endif /* DEBUGGING */



/***************************************************************
 *  set_pal:  Reads HDF palette and makes it default.          *
 ***************************************************************/

void set_pal( char *palfile )
{
  static PAL palette[768];       /* palette storage */

  if ( DFPgetpal( palfile, (char *)palette) )
    quit("set_pal: Unable to open palette file.");
  if ( DFR8setpalette( (char *)palette ) )
    quit("set_pal: Unable to set default palette.");
}

/*********************************************************************************
 * read_crv:  Reads curvilinear data from three straight HDF file containing     *
 *            one SDS each:  data, gridx, gridy.                                 *
 *********************************************************************************/

void read_crv( struct Input *in, struct Options *opt )
{
  int   rank;
  int32 hdfdims[2];
  
  if ( DFSDgetdims(opt->infile, &rank, hdfdims, 2) ) 
    quit("read_crv: Unable to get dimensions from data file.");

  if( rank != 2)
    quit("read_crv: Invalid data rank in data file.");

  if ((hdfdims[0] < 2) || (hdfdims[1] < 2)) 
    quit("read_crv: Dimension(s) less than '2' in data file.");

  in->dim0 = hdfdims[0];
  in->dim1 = hdfdims[1];

  if ( opt->minmax == FALSE               && 
       !DFSDgetmaxmin(&datamax, &datamin) &&
       datamax > datamin )  
    opt->minmax = TRUE;

  if ( DFSDgetdims(opt->gridx, &rank, hdfdims, 2) ||
       rank != 2              || 
       hdfdims[0] != in->dim0 ||
       hdfdims[1] != in->dim1    ) 
    quit("read_crv: Incompatible gridx file.");

  if ( DFSDgetdims(opt->gridy, &rank, hdfdims, 2) ||
       rank != 2              || 
       hdfdims[0] != in->dim0 ||
       hdfdims[1] != in->dim1    ) 
    quit("read_crv: Incompatible gridy file.");

  if ((in->x    = (float32 *) malloc((size_t) (in->dim0*in->dim1*sizeof(float32)) )) == NULL || 
      (in->y    = (float32 *) malloc((size_t) (in->dim0*in->dim1*sizeof(float32)) )) == NULL || 
      (in->data = (float32 *) malloc((size_t) (in->dim0*in->dim1*sizeof(float32)) )) == NULL)
    quit("read_crv: Cannot allocate mem for grid and/or data.");

  if (DFSDgetdata(opt->infile, 2, hdfdims, in->data) ||
      DFSDgetdata(opt->gridx, 2, hdfdims, in->x)     ||
      DFSDgetdata(opt->gridy, 2, hdfdims, in->y)     ) 
    quit("read_crv: Unable to get SDS from input file(s)." );

  if ( opt->domain == FALSE )
  { gminmax( &xmin, &xmax, in->x, in->dim0*in->dim1 );
    gminmax( &ymin, &ymax, in->y, in->dim0*in->dim1 );
  }
  if ( opt->minmax == FALSE )  
    gminmax( &datamin, &datamax, in->data, in->dim0*in->dim1 );

  if (datamin >= datamax)           quit("read_crv: Data max/min err.");
  if (xmin >= xmax || ymin >= ymax) quit("read_crv: Domain spec err.");
}

/****************************************************
 *  read_vset:  Read HDF Vset file.                 *
 ****************************************************/

void read_vset( struct Input *in, struct Options *opt )
{
  VGROUP *vg;    
  VDATA  *vs;
  DF     *f;
  char   vgname[STR_MAX], vsname[STR_MAX], fields[STR_MAX];
  int    vgid, vsid;
  enum boolean 
         found, pxfound, pyfound, scalarfound, 
         connectfound[CMAX], confound, oneconfound, con3found, con4found;
  int    i, csize, vsize, interlace;
  int    nelements;            /* --- hopefully int==int32 ... --- */
  int32  npx, npy, nscalar;   


  if ( (f = DFopen( opt->infile, DFACC_ALL, 0 )) == NULL)
    quit("read_vset: Cannot open input file.");

  /*
   * Now try to find vgroup with name opt->groupname.  
   * If it equals "" (empty string) then the first vgroup is taken.
   */

  found = FALSE;
  vgid  = -1;

  while ((vgid = Vgetid(f,vgid)) != -1)
  { if ((vg = Vattach(f, vgid, "r")) != NULL)
    { Vgetname(vg, vgname);
      if (strcmp(vgname,opt->groupname) == 0 || strcmp(opt->groupname,"") == 0 )
      { found = TRUE;
        break;
      }
    }
    Vdetach(vg);
  }

  if ( found == FALSE )
    quit("read_vset: Group not found.");

  /*
   * Now scan all the vdatas in the vgroup found.
   * Read the px, py, connect, scalar fields if encountered.
   */

  pxfound = pyfound = scalarfound = 
            confound = oneconfound = con3found = con4found = FALSE;
  for (i=0; i<CMAX; i++)  connectfound[i] = FALSE;
  vsid = -1;

  while ((vsid = Vgetnext(vg, vsid)) != -1 &&
         (pxfound == FALSE || pyfound == FALSE || 
          scalarfound == FALSE || confound == FALSE) )
  {
    if (Visvs(vg,vsid))             /* the vsid is indeed referring to a vdata */
    {
      if ((vs = VSattach( f, vsid, "r" )) != NULL &&
          VSinquire(vs, &nelements, &interlace, fields, &vsize, vsname) == 1 &&
          nelements > 0 )
      {
#ifdef DEBUGGING
        printf("vdata %s has %d elements of size %d, and fields %s\n",
                vsname, nelements, vsize, fields );
#endif
        if ( pxfound == FALSE &&
             VSfexist(vs, opt->pxname) == 1 && 
             VSsetfields(vs, opt->pxname) == 1 )
        {
          if ((in->x = (float32 *) malloc((size_t) (nelements*sizeof(float32)))) == NULL)
            quit("read_vset: Cannot allocate mem for px.");
          VSread(vs, (unsigned char *)(in->x), nelements, NO_INTERLACE );
          pxfound = TRUE;
          npx = nelements;
#ifdef DEBUGGING
          printf("   px read\n" );
#endif
        }

        if ( pyfound == FALSE &&
             VSfexist(vs, opt->pyname) == 1 && 
             VSsetfields(vs, opt->pyname) == 1 )
        {
          if ((in->y = (float32 *) malloc((size_t) (nelements*sizeof(float32)))) == NULL)
            quit("read_vset: Cannot allocate mem for py.");
          VSread(vs, (unsigned char *)(in->y), nelements, NO_INTERLACE );
          pyfound = TRUE;
          npy = nelements;
#ifdef DEBUGGING
          printf("   py read\n" );
#endif
        }

        if ( scalarfound == FALSE &&
             VSfexist(vs, opt->scalarname) == 1 && 
             VSsetfields(vs, opt->scalarname) == 1 )
        {
          if ((in->data = (float32 *) malloc((size_t) (nelements*sizeof(float32)))) == NULL)
            quit("read_vset: Cannot allocate mem for scalar.");
          VSread(vs, (unsigned char *)(in->data), nelements, NO_INTERLACE );
          scalarfound = TRUE;
          nscalar = nelements;
#ifdef DEBUGGING
          printf("   scalar read\n" );
#endif
        }

        for (i=0; i < opt->ncon; i++)
          if ( connectfound[i] == FALSE &&
               VSfexist(vs, opt->connectname[i]) == 1 && 
               VSsetfields(vs, opt->connectname[i]) == 1 &&
               (csize = VSsizeof(vs, opt->connectname[i])) > 0)
          {
#ifdef DEBUGGING
            printf("   attempting conn list named %s\n", opt->connectname[i] );
            printf("              size found %d\n", csize );
#endif
            if ( con3found == FALSE &&
                 csize == 3*sizeof(LOC_INT) )       /* LOCAL_INTTYPE triangle conn list */
            {
              if ((in->splist3 = (LOC_INT *) malloc((size_t) ((int32)nelements*csize))) == NULL)
                quit("read_vset: Cannot allocate mem for connectivity list.");
              VSread(vs, (unsigned char *)(in->splist3), nelements, NO_INTERLACE );
              in->etyp3 = INTTYPE;
              in->ntri  = nelements;
              con3found = TRUE;
#ifdef DEBUGGING
              printf("              it's a tri-INTTYPE\n" );
#endif
            }
            else if ( con4found == FALSE &&
                      csize == 4*sizeof(LOC_INT) )  /* LOCAL_INTTYPE quadrangle conn list */
            {
              if ((in->splist4 = (LOC_INT *) malloc((size_t) ((int32)nelements*csize))) == NULL)
                quit("read_vset: Cannot allocate mem for connectivity list.");
              VSread(vs, (unsigned char *)(in->splist4), nelements, NO_INTERLACE );
              in->etyp4 = INTTYPE;
              in->nquad = nelements;
              con4found = TRUE;
#ifdef DEBUGGING
              printf("              it's a quad-INTTYPE\n" );
#endif
            }
            else if ( con3found == FALSE &&
                      csize == 3*sizeof(LOC_LONG) )  /* LOCAL_LONGTYPE triangle conn list */
            {
              if ((in->plist3 = (LOC_LONG *) malloc((size_t) ((int32)nelements*csize))) == NULL)
                quit("read_vset: Cannot allocate mem for connectivity list.");
              VSread(vs, (unsigned char *)(in->plist3), nelements, NO_INTERLACE );
              in->etyp3 = LONGTYPE;
              in->ntri  = nelements;
              con3found = TRUE;
#ifdef DEBUGGING
              printf("              it's a tri-LONGTYPE\n" );
#endif
            }
            else if ( con4found == FALSE &&
                      csize == 4*sizeof(LOC_LONG) )  /* LOCAL_LONGTYPE quadrangle conn list */
            {
              if ((in->plist4 = (LOC_LONG *) malloc((size_t) ((int32)nelements*csize))) == NULL)
                quit("read_vset: Cannot allocate mem for connectivity list.");
              VSread(vs, (unsigned char *)(in->plist4), nelements, NO_INTERLACE );
              in->etyp4 = LONGTYPE;
              in->nquad = nelements;
              con4found = TRUE;
#ifdef DEBUGGING
              printf("              it's a quad-LONGTYPE\n" );
#endif
            }
            else
              quit("read_vset: Illegal or duplicate connectivity list.");
            connectfound[i] = TRUE;
          }
        confound = TRUE;
        oneconfound = FALSE;
        for (i=0; i < opt->ncon; i++)
        { confound = confound && connectfound[i];
          oneconfound = oneconfound || connectfound[i];
        }
        confound = confound || ( con3found && con4found );
      }
      VSdetach(vs);
    }
  }

  Vdetach(vg);
  DFclose(f);

  if ( pxfound == FALSE )  quit("read_vset: Cannot find px.");
  if ( pyfound == FALSE )  quit("read_vset: Cannot find py.");

  if ( scalarfound == FALSE )  quit("read_vset: Cannot find scalar.");
  if ( oneconfound == FALSE )  quit("read_vset: Cannot find connect.");
  if ( confound == FALSE && opt->connect == TRUE )  
    printf("read_vset: Warning: Couldn't load/find all connect lists you specified.\n");

  if ( npx != npy  ||  npx != nscalar )
    quit("read_vset: (px, py, scalar) have unequal lengths.");
  else
    in->nv = npx;

  if ( opt->domain == FALSE )
  { gminmax( &xmin, &xmax, in->x, in->nv );
    gminmax( &ymin, &ymax, in->y, in->nv );
  }
  if ( opt->minmax == FALSE )  
    gminmax( &datamin, &datamax, in->data, in->nv );

  if (datamin >= datamax)           quit("read_vset: Data max/min err.");
  if (xmin >= xmax || ymin >= ymax) quit("read_vset: Domain spec err.");
}

/*****************************************************************
 *  gminmax:  Determines min and max of float32 array x[dim].    *
 *****************************************************************/

void gminmax( float32 *xmin, float32 *xmax, float32 *x, int32 dim )
{
  for ( *xmin = *x, *xmax = *x; dim>1; dim-- )  
  {
    x++;
    if (*x < *xmin) *xmin = *x;
    else if (*x > *xmax) *xmax = *x;
  }
}

/***************************************************************
 *  setup_image:  Allocates storage for image, and sets        *
 *                background colour.                           *
 ***************************************************************/

void setup_image( struct Raster *im )
{
  int32 i;
  PIXEL *ras;

  if ((im->image = (PIXEL *) malloc((size_t) (im->hres*im->vres*sizeof(PIXEL)))) == NULL)
    quit("setup_image: Cannot allocate mem for image.");

  /* 
   * Set background color:  
   */  

  ras = im->image;

  for ( i=0; i < im->hres*im->vres; i++ )             
    *ras++ = BACKGR_COL;
}

/*******************************************************************
 * write_image:  Writes raster image to RIS8 file with name        *
 *               specified by the arg file[].                      *
 *******************************************************************/

void write_image( char *file, struct Raster im )
{

  if ( DFR8putimage( file, (char *)im.image, im.hres, im.vres, DFTAG_RLE) )
    printf("write_image: Image could not be written.\n");
#if ( defined(MAC) && defined(MPW) )
  else
    fsetfileinfo( file, 'NImg', 'NIMG' );  /* creator resp type of NCSA Image RIS8 file */
#endif
}

/******************************************************************
 * curvi2ris:   Generates raster of data defined on curvilinear   *
 *              grid.                                             *
 ******************************************************************/

#define  DATA(i,j)  *(in.data + (i)*in.dim1 + (j))
#define  X(i,j)     *(in.x    + (i)*in.dim1 + (j))
#define  Y(i,j)     *(in.y    + (i)*in.dim1 + (j))

void curvi2ris( struct Input in, struct Raster *im )
{
  int32   i, j;
  float32 x[4], y[4], dat[4];

  for ( i=0; i<in.dim0-1; i++ )
    for ( j=0; j<in.dim1-1; j++ )   /* plot each gridbox */
    {
      x[0] = X(i,  j  );
      x[1] = X(i+1,j  );
      x[2] = X(i+1,j+1);
      x[3] = X(i,  j+1);

      y[0] = Y(i,  j  );
      y[1] = Y(i+1,j  );
      y[2] = Y(i+1,j+1);
      y[3] = Y(i,  j+1);

      dat[0] = DATA(i,  j  );
      dat[1] = DATA(i+1,j  );
      dat[2] = DATA(i+1,j+1);
      dat[3] = DATA(i,  j+1);

      if ( iso == TRUE )     /* bilinear elements */
        iso_quadrangle( x, y, dat, im );
      else                   /* subdivide into two linear (triangular) elements */
        quadrangle( x, y, dat, im );
    }
}

/******************************************************************
 * vset2ris:   Generates raster of data defined on vertex set.    *
 ******************************************************************/

void vset2ris( struct Input in, struct Raster *im )
{  
  int32   i, j, ptr;
  float32 x[4], y[4], dat[4];


  /*
   *  First the triangles.
   */

  for (i=0; i < in.ntri; i++)
  {
    for (j=0; j<3; j++)
    {
      ptr = (in.etyp3==INTTYPE)? in.splist3[3*i+j] - 1
                               : in.plist3[3*i+j] - 1;
      if ( ptr < 0  ||  ptr > in.nv-1 )
        quit("vset2ris: Connectivity list element out of range.");
      x[j]   = in.x[ptr];
      y[j]   = in.y[ptr];
      dat[j] = in.data[ptr];
    }

    if ( isot == TRUE )  
      iso_triangle( x, y, dat, im );
    else  
      triangle( x, y, dat, im );
  }

  /*
   *  Finally the quadrangles.
   */

  for (i=0; i < in.nquad; i++)
  {
    for (j=0; j<4; j++)
    {
      ptr = (in.etyp4==INTTYPE)? in.splist4[4*i+j] - 1
                               : in.plist4[4*i+j] - 1;
      if ( ptr < -1  ||  (ptr==-1 && j!=3)  ||  ptr > in.nv-1 )
        quit("vset2ris: Connectivity list element out of range.");
      x[j]   = in.x[ptr];
      y[j]   = in.y[ptr];
      dat[j] = in.data[ptr];
    }

    if ( ptr != -1 )
    { if ( iso == TRUE ) 
        iso_quadrangle( x, y, dat, im );
      else      
        quadrangle( x, y, dat, im );
    }
    else       /* last ptr == -1 means triangle stored in quadr conn list */
    { if ( isot == TRUE )  
        iso_triangle( x, y, dat, im );
      else  
        triangle( x, y, dat, im );
    }
  }
}

/***************************************************************
 *  triangle: Fills the raster-pixels inside the triangle      *
 *            specified by (x,y,dat)[3].  This version uses an *
 *            extension of the "intersection" approach.        *
 ***************************************************************/

#define SMALL  0.001                        /* Needed to avoid rounding problems */

void triangle( float32 x[3], float32 y[3], float32 dat[3], 
               struct Raster *im )
{ 
  int     ilo, imid, ihi, itop;             /* ordering indices of vertices */
  float32 xleft, xright, datleft, datright; /* artificially constructed "base" */

  int32   i, j, imin, imax, jmin, jmax;
  float32 lower, upper, xs, x1, x2, dat1, xx, yy, dd, dx, dy, ratio, ddatdx, tmp;

  int do_lower, do_upper;                   /* booleans */

/*
 *  step 1. Determine vertex ordering in y-direction:  
 *                   y[ilo] <= y[imid] <= y[ihi].
 *          The triple (ilo,imid,ihi) is a permutation of (0,1,2).
 */
  
  ihi = y[1]>y[0]; 
  ilo = !ihi; 

  if ( y[2]>y[ihi] ) 
  { imid = ihi;
    ihi = 2; }
  else if ( y[2]<y[ilo] )
  { imid = ilo;
    ilo = 2; }
  else
    imid = 2;

/*
 * step 2. Split the triangle into two triangles with common "base" 
 *         parallel to x-axis:  (xleft,y[imid]) -- (xright,y[imid]).
 */

  ratio = (y[imid]-y[ilo])/(y[ihi]-y[ilo]);
  xs    = x[ilo] + ratio*(x[ihi]-x[ilo]);

  if ( x[imid] < xs )
  { xleft    = x[imid];
    datleft  = dat[imid];
    xright   = xs;
    datright = dat[ilo] + ratio*(dat[ihi]-dat[ilo]); }
  else
  { xleft    = xs;
    datleft  = dat[ilo] + ratio*(dat[ihi]-dat[ilo]);
    xright   = x[imid];
    datright = dat[imid]; }

/*
 * step 3. Calculate pixel values on lines parallel to the x-axis, 
 *         but in triangle interior:
 */

  dx = (xmax-xmin) / (im->hres-1);
  dy = (ymax-ymin) / (im->vres-1);

  ddatdx = (datright-datleft) / (xright-xleft); /* data x-derivative is constant on triangle */

  lower = ( y[imid] - y[ilo] )/dy;
  upper = ( y[ihi] - y[imid] )/dy;

  do_lower = lower > SMALL/2 || lower >= upper; /* check for degeneracies of "triangle-halfs" */
  do_upper = upper > SMALL/2 || upper > lower;

  jmin = ROUNDUP(   (y[ilo]-ymin)/dy - SMALL );
  jmax = ROUNDDOWN( (y[ihi]-ymin)/dy + SMALL );       
  
  if (jmin<0)          jmin = 0;
  if (jmax>im->vres-1) jmax = im->vres-1;       /* domain clipping */
  
  for ( j=jmin; j<=jmax; j++)
  { /* 
     * Calculate the two intersections of line yy==ymin+j*dy
     * and the line through (xleft,y[imid]) --(x[itop],y[itop])
     *     and line through (xright,y[imid])--(x[itop],y[itop]) 
     * with itop = ihi or ilo (depending on yy).
     */
    yy   = ymin + j*dy;
    itop = ( yy>y[imid] && do_upper || !do_lower )? ihi: ilo;
    ratio= (yy-y[imid]) / (y[itop]-y[imid]);
  
    x1   = xleft   + (x[itop]-xleft)*ratio;
    dat1 = datleft + (dat[itop]-datleft)*ratio;
    x2   = xright  + (x[itop]-xright)*ratio; 
  
    imin = ROUNDUP(   (x1-xmin)/dx - SMALL );
    imax = ROUNDDOWN( (x2-xmin)/dx + SMALL );   /* expand slightly, otherwise some pixels may be lost */
  
    if (imin<0)          imin = 0;
    if (imax>im->hres-1) imax = im->hres-1;     /* domain clipping */
  
    for ( i=imin; i<=imax; i++)
    { xx = xmin + i*dx;
      dd = dat1 + (xx-x1)*ddatdx;
      if (dd>datamax) dd=datamax;
      if (dd<datamin) dd=datamin;               /* data clipping */
      IMAGE(i,j) = 1.5 + (dd-datamin)/(datamax-datamin)*237.9; 
    }
  }
}


/***************************************************************************
 * iso_triangle: Fills raster-pixels inside the triangle specified         *
 *               by (x,y,dat)[3].  This version uses an iso-parametric     *
 *               mapping to determine the pixel-value.                     *
 ***************************************************************************/

#define SMALL   0.001   /* Triangle is expanded slightly in isoparametric space */

void iso_triangle( float32 x[3], float32 y[3], float32 dat[3], 
                   struct Raster *im )
{ 
  int32   i, j;
  float32 a11, a12, a21, a22, det, rhs1, rhs2, xi, eta, xx, yy, dd, tmp;                

  int32   imin, imax, jmin, jmax;     /* pixel coords of rectangular hull */

  float32 dx = (xmax-xmin) / (im->hres-1);
  float32 dy = (ymax-ymin) / (im->vres-1);

  /* 
   * step 1. Determine the "rectangular hull" of the triangle
   */

  imin = ROUNDUP(   ( MIN3(x[0],x[1],x[2]) - xmin )/dx - SMALL );
  imax = ROUNDDOWN( ( MAX3(x[0],x[1],x[2]) - xmin )/dx + SMALL );
  jmin = ROUNDUP(   ( MIN3(y[0],y[1],y[2]) - ymin )/dy - SMALL );
  jmax = ROUNDDOWN( ( MAX3(y[0],y[1],y[2]) - ymin )/dy + SMALL );

  if (imin<0)          imin = 0;
  if (imax>im->hres-1) imax = im->hres-1;
  if (jmin<0)          jmin = 0;
  if (jmax>im->vres-1) jmax = im->vres-1;     /* domain clipping */
  
  /* 
   * step 2. For each pixel in the rectangular hull check if it is 
   *         inside the triangle.
   *             If it is: calculate and store pixel value
   */

  a11 = x[1] - x[0];
  a12 = x[2] - x[0];
  a21 = y[1] - y[0];
  a22 = y[2] - y[0];
  det = a11*a22 - a12*a21;

  for (i=imin; i<=imax; i++)
  { 
    xx = xmin + i*dx;
    rhs1 = xx - x[0];

    for (j=jmin; j<=jmax; j++)
    { 
      yy = ymin + j*dy;

      /* calculate iso-parametric coords (xi,eta), implicitly defined by:  
       *     xx = x[0] + xi*(x[1]-x[0]) + eta*(x[2]-x[0])
       *     yy = y[0] + xi*(y[1]-y[0]) + eta*(y[2]-y[0]) 
       * this mapping maps the triangle to a simple triangle in xi,eta space:
       *     (x[0],y[0]) ---> (0,0)
       *     (x[1],y[1]) ---> (1,0)
       *     (x[2],y[2]) ---> (0,1)
       */

      rhs2 = yy - y[0];

      xi  = (  a22*rhs1 - a12*rhs2 ) / det;
      eta = ( -a21*rhs1 + a11*rhs2 ) / det;

      if ( xi>-SMALL && eta>-SMALL && xi+eta<1.0+SMALL )   /* interior point */
      { 
        dd = dat[0] + xi*(dat[1]-dat[0]) + eta*(dat[2]-dat[0]);
        if (dd>datamax) dd=datamax;
        if (dd<datamin) dd=datamin;                        /* data clipping */
        IMAGE(i,j) = 1.5 + (dd-datamin)/(datamax-datamin)*237.9;
      }
    }
  }
}

/***************************************************************************
 * quadrangle: Fills raster-pixels inside the quadrangle specified         *
 *             by (x,y,dat)[4].  Splits the quadrangle into two triangles. *
 *    ASSUMES: THE FOUR VERTICES ARE ORDERED.                              *
 ***************************************************************************/

void quadrangle( float32 x[4], float32 y[4], float32 dat[4], 
                 struct Raster *im )
{
  float32 xt[3], yt[3], datt[3];       /* one of the two triangles */

  int id = diagonal_choice( x, y ); 
                      /* first vertex index of diagonal which must be used
                         to split quadrangle into two triangles */

  /* 
   *  first the triangle (x,y,dat)[id,id+1,id+2],  where id=0 or 1.
   */

  if ( isot == TRUE )
    iso_triangle( &x[id], &y[id], &dat[id], im );
  else
    triangle( &x[id], &y[id], &dat[id], im );

  /* 
   *  now the other triangle:  (x,y,dat)[0,2-id,3] 
   */

  xt[0] = x[0];
  xt[1] = x[2-id];
  xt[2] = x[3];

  yt[0] = y[0];
  yt[1] = y[2-id];
  yt[2] = y[3];

  datt[0] = dat[0];
  datt[1] = dat[2-id];
  datt[2] = dat[3];

  if ( isot == TRUE )
    iso_triangle( xt, yt, datt, im );
  else
    triangle( xt, yt, datt, im );
}

/***************************************************************************
 * diagonal_choice:  Determines interior diagonal in quadrangle specified  *
 *                   by (x,y)[4].  Returns 0 or 1.  Value of 0 means       *
 *                   "diagonal (x,y)[0,2]".  Value of 1 means              *
 *                   "diagonal (x,y)[1,3]".                                *
 ***************************************************************************/

#define WALK_SIGN(i,j,k) ( (x[k]-x[j]) * (y[i]-y[j]) - \
                           (y[k]-y[j]) * (x[i]-x[j]) > 0 )
/*
 *  Determines sign of directed walk around triangle with vertices (x,y)[i,j,k]:
 *    WALK_SIGN == 1 means positive (counter-clockwise),
 *    WALK_SIGN == 0 means negative (clockwise).
 */

/* if quad is convex: make an optimal choice based on area-ratio of triangles ??? */

int diagonal_choice( float32 x[4], float32 y[4] )  
{
  return ( WALK_SIGN(0,1,2) != WALK_SIGN(2,3,0) );
}

/*************************************************************************
 * convex:  Returns 1 if quadrangle specified by (x,y)[4] is convex.     * 
 *          Returns 0 if not convex.                                     *
 *************************************************************************/

int convex( float32 x[4], float32 y[4] )
{
  return ( WALK_SIGN(0,1,2) == WALK_SIGN(2,3,0) &&
           WALK_SIGN(1,2,3) == WALK_SIGN(3,0,1)  );
}

/********************************************************************************
 * simple_iso_quadrangle:                                                       *
 *                 Fills raster-pixels inside the quadrangle specified          *
 *                 by (x,y,dat)[4].  This version uses a bilinear isoparametric *
 *                 mapping, and a simple initial pixel selection.               *
 *        ASSUMES: THE FOUR VERTICES ARE ORDERED.                               *
 *           BUGS: FOR EXTERIOR PIXELS: RISK OF SINGULAR JACOBIAN.              *
 ********************************************************************************/

/* 
 * Defs needed for Newton inversion of isoparametric mapping:
 */

#define RESX(xi,eta) ( xx - ( x[0] + (xi)*(x[1]-x[0]) +\
                       (eta)*(x[3]-x[0]) + (xi)*(eta)*(x[2]+x[0]-x[1]-x[3]) ) )
#define RESY(xi,eta) ( yy - ( y[0] + (xi)*(y[1]-y[0]) +\
                       (eta)*(y[3]-y[0]) + (xi)*(eta)*(y[2]+y[0]-y[1]-y[3]) ) )
#define MAXIT   15     /* max allowed nr of Newton iterations; usually 2 or 3 needed */
#define DIVERG  1.0E3  /* divergence check in Newton iteration; 
                          should be large enough to account for exterior pixels */
#define EPS     0.001  /* convergence check in Newton iteration */

#define SMALL   0.001  /* Quadrangle is expanded slightly in isoparametric space. */



void simple_iso_quadrangle( float32 x[4], float32 y[4], float32 dat[4], 
                            struct Raster *im )
{ 
  int32   i, j, iter, noconverg;
  float32 a11, a12, a21, a22, det, rhs1, rhs2, 
          xi, eta, prev_xi, prev_eta, xx, yy, dd, tmp;                

  int32   imin, imax, jmin, jmax;   /* pixel coords of rectangular hull */

  float32 dx = (xmax-xmin) / (im->hres-1);
  float32 dy = (ymax-ymin) / (im->vres-1);

  if ( !convex(x,y) )               /* use more robust approach instead */
  { 
    printf("simple_iso_quadrangle: Quadrangle isn't convex; switching to triangularization.\n");
    quadrangle( x, y, dat, im );  
    return;
  }

  /* 
   * step 1. Determine the "rectangular hull" of the triangle
   */

  imin = ROUNDUP(   ( MIN4(x[0],x[1],x[2],x[3]) - xmin )/dx - SMALL );
  imax = ROUNDDOWN( ( MAX4(x[0],x[1],x[2],x[3]) - xmin )/dx + SMALL );
  jmin = ROUNDUP(   ( MIN4(y[0],y[1],y[2],y[3]) - ymin )/dy - SMALL );
  jmax = ROUNDDOWN( ( MAX4(y[0],y[1],y[2],y[3]) - ymin )/dy + SMALL );

  if (imin<0)          imin = 0;
  if (imax>im->hres-1) imax = im->hres-1;
  if (jmin<0)          jmin = 0;
  if (jmax>im->vres-1) jmax = im->vres-1;     /* domain clipping */
  
  /* 
   * step 2. For each pixel in the rectangular hull check if it is 
   *         inside the quadrangle. If it is: calculate and store pixel value.
   */

  for (i=imin; i<=imax; i++)
  { 
    xx = xmin + i*dx;

    for (j=jmin; j<=jmax; j++)
    { 
      yy = ymin + j*dy;

      /* Calculate iso-parametric coords (xi,eta), implicitly defined by:  
       *     xx = x[0] + xi*(x[1]-x[0]) + eta*(x[3]-x[0]) 
       *                                + xi*eta*(x[2]+x[0]-x[1]-x[3])
       *     yy = y[0] + xi*(y[1]-y[0]) + eta*(y[3]-y[0])
       *                                + xi*eta*(y[2]+y[0]-y[1]-y[3])
       * This mapping maps the quadrangle to the unit square in xi,eta space:
       *     (x[0],y[0]) ---> (0,0)
       *     (x[1],y[1]) ---> (1,0)
       *     (x[2],y[2]) ---> (1,1)
       *     (x[3],y[3]) ---> (0,1)
       * For each pixel, the "physical coords" xx,yy are known.
       * The corresponding xi,eta must be found from a Newton process, since 
       * the isoparametric mapping is nonlinear.
       */

      xi   = prev_xi  = 0.5;
      eta  = prev_eta = 0.5;   /* initial guess for Newton iteration */
               /* If nr of pixels for this quadrangle is large, 
                  a more efficient guess is to use "previous" values */
      iter = 0;

      do {
        iter++;

        /* the current residual */
        rhs1 = RESX(xi,eta);
        rhs2 = RESY(xi,eta);

        /* the Jacobian of the mapping */
        a11 = x[1] - x[0] + eta*(x[2]+x[0]-x[1]-x[3]);
        a12 = x[3] - x[0] +  xi*(x[2]+x[0]-x[1]-x[3]);
        a21 = y[1] - y[0] + eta*(y[2]+y[0]-y[1]-y[3]);
        a22 = y[3] - y[0] +  xi*(y[2]+y[0]-y[1]-y[3]);
        det = a11*a22 - a12*a21;

        xi  = xi  + (  a22*rhs1 - a12*rhs2 ) / det;
        eta = eta + ( -a21*rhs1 + a11*rhs2 ) / det;

/*      
 *       This robust convergence criterium is too complex for Macintosh MPW C:
 *
 *       noconverg = ( RESX(xi+EPS,eta)*RESX(xi-EPS,eta) > 0 ||
 *                     RESX(xi,eta+EPS)*RESX(xi,eta-EPS) > 0 ||
 *                     RESY(xi+EPS,eta)*RESY(xi-EPS,eta) > 0 ||
 *                     RESY(xi,eta+EPS)*RESY(xi,eta-EPS) > 0 );   
 */

        noconverg = ( fabs(xi-prev_xi) + fabs(eta-prev_eta) > EPS );
        prev_xi  = xi;
        prev_eta = eta;
      } 
      while( noconverg && iter<MAXIT && fabs(xi)+fabs(eta) < DIVERG );


      if (noconverg)
      { 
        printf("simple_iso_quadrangle: Newton failed; switching to triangularization.\n");
        quadrangle( x, y, dat, im );  /* use more robust approach instead */
        return;
      }

      if ( xi>-SMALL && eta>-SMALL && xi<1.0+SMALL && eta<1.0+SMALL )  /* interior point */
      { 
        dd = dat[0] + xi*(dat[1]-dat[0]) + eta*(dat[3]-dat[0]) + 
             xi*eta*(dat[2]+dat[0]-dat[1]-dat[3]);
        if (dd>datamax) dd=datamax;
        if (dd<datamin) dd=datamin;                                    /* data clipping */
        IMAGE(i,j) = 1.5 + (dd-datamin)/(datamax-datamin)*237.9;
      }
    }
  }
}
#undef DIVERG     
         /* ^ different for next implementation */


/********************************************************************************
 * iso_quadrangle: Fills raster-pixels inside the quadrangle specified          *
 *                 by (x,y,dat)[4].  This version uses a bilinear isoparametric *
 *                 mapping.                                                     *
 *        ASSUMES: THE FOUR VERTICES ARE ORDERED.                               *
 ********************************************************************************/

/* 
 * Defs needed for Newton inversion of isoparametric mapping:
 */

#define RESX(xi,eta) ( xx - ( x[0] + (xi)*(x[1]-x[0]) +\
                       (eta)*(x[3]-x[0]) + (xi)*(eta)*(x[2]+x[0]-x[1]-x[3]) ) )
#define RESY(xi,eta) ( yy - ( y[0] + (xi)*(y[1]-y[0]) +\
                       (eta)*(y[3]-y[0]) + (xi)*(eta)*(y[2]+y[0]-y[1]-y[3]) ) )
#define MAXIT   15     /* max allowed nr of Newton iterations; usually 2 or 3 needed */
#define DIVERG  3.0    /* divergence check in Newton iteration; 
                          should be slightly larger than 2.0 */
#define EPS     0.001  /* convergence check in Newton iteration */

#ifdef DEBUGGING
  int32 kount[MAXIT+1] = {0};

  void write_iter( void )
  { int i;
    printf("\n\n--- Newton iteration count ---\n");

    for (i=1; i<=MAXIT; i++)
      printf("nr times Newton needed %d iters: %ld\n", i, (long)(kount[i]) ); 
  }
#endif

#define SMALL  0.001    /* Needed to avoid rounding problems */



void iso_quadrangle( float32 x[4], float32 y[4], float32 dat[4], 
                     struct Raster *im )
{ 
  int     iter, noconverg;
  float32 a11, a12, a21, a22, det, rhs1, rhs2, 
          xi, eta, prev_xi, prev_eta;               /* isoparametric stuff */         

  int32 i, j, imin, imax, jmin, jmax;
  int   ilo, imid1, imid2, ihi;
  int   do_lower, do_mid, do_upper;                 /* booleans */

  float32 dx, dy, xx, yy, dd, upper, mid, lower, 
          x1, x2, xs1, xs2, xleft1, xleft2, xright1, xright2, ratio, tmp;

  if ( !convex(x,y) )               /* use more robust approach instead */
  { 
    printf("iso_quadrangle: Quadrangle isn't convex; switching to triangularization.\n");
    quadrangle( x, y, dat, im );  
    return;
  }

/*
 *  step 1. Determine vertex ordering in y-direction:
 *              y[ilo] <= y[imid1] <= y[imid2] <= y[ihi].
 *          The quadruple (ilo,imid1,imid2,ihi) is a permutation of (0,1,2,3).
 */

  ihi = y[1]>y[0]; 
  ilo = !ihi; 

  if ( y[2]>y[ihi] ) 
  { imid1 = ihi;
    ihi = 2; }
  else if ( y[2]>y[ilo] )
    imid1 = 2;
  else
  { imid1 = ilo;
    ilo = 2; }
  
  if ( y[3]>y[ihi] )
  { imid2 = ihi;
    ihi = 3; }
  else if ( y[3]>y[imid1] )
    imid2 = 3;
  else if ( y[3]>y[ilo] )
  { imid2 = imid1;
    imid1 = 3; }
  else
  { imid2 = imid1;
    imid1 = ilo;
    ilo = 3; }
/*
#ifdef DEBUGGING
  printf("ordering(lo->hi): %d%d%d%d\n", ilo, imid1, imid2, ihi );
#endif
*/

/*
 *  step 2. Split the CONVEX quadrangle into two triangles and a simple
 *          quadrangle, with two common "bases" parallel to the x-axis:
 *              top  =             (x[ihi],y[ihi])
 *             base2 = (xleft2,y[imid2]) -- (xright2,y[imid2])
 *             base1 = (xleft1,y[imid1]) -- (xright1,y[imid1])
 *              bot  =             (x[ilo],y[ilo])
 */

  dx = (xmax-xmin) / (im->hres-1);
  dy = (ymax-ymin) / (im->vres-1);

  /* First check for degeneracies of "quad-parts" */

  upper = ( y[ihi] - y[imid2] )/dy;
  mid   = ( y[imid2] - y[imid1] )/dy;
  lower = ( y[imid1] - y[ilo] )/dy;

  do_upper = upper > SMALL/3 || ( upper > mid && upper > lower );
  do_mid   = mid > SMALL/3   || ( mid >= upper && mid >= lower );
  do_lower = lower > SMALL/3 || ( lower > mid && lower > upper ); 

  /*
   * Two possibilities: 
   * either (x,y)[ilo]--(x,y)[ihi] is a quadrangle side, or a diagonal.
   */

  xs2 = x[ilo] + (y[imid2]-y[ilo])/(y[ihi]-y[ilo])*(x[ihi]-x[ilo]);
  xs1 = x[ilo] + (y[imid1]-y[ilo])/(y[ihi]-y[ilo])*(x[ihi]-x[ilo]);

  if ( abs(ihi-ilo) == 2 )  /* --- it's a DIAGONAL --- */
  { 
    if ( do_upper || do_mid )
    { 
      if ( x[imid2] < xs2 || x[imid1] > xs1 )
      { xleft2  = x[imid2];
        xright2 = x[imid1] + (y[imid2]-y[imid1])/(y[ihi]-y[imid1])*(x[ihi]-x[imid1]); }
      else
      { xleft2  = x[imid1] + (y[imid2]-y[imid1])/(y[ihi]-y[imid1])*(x[ihi]-x[imid1]);
        xright2 = x[imid2]; }
    }
    if ( do_lower || do_mid )
    { 
      if ( x[imid2] < xs2 || x[imid1] > xs1 )
      { xleft1  = x[ilo] + (y[imid1]-y[ilo])/(y[imid2]-y[ilo])*(x[imid2]-x[ilo]);
        xright1 = x[imid1]; }
      else
      { xleft1  = x[imid1];
        xright1 = x[ilo] + (y[imid1]-y[ilo])/(y[imid2]-y[ilo])*(x[imid2]-x[ilo]); }
    }
  }
  else                      /* --- it's a quadrangle SIDE --- */
  {
    if ( x[imid2] < xs2 ) 
    { xleft2  = x[imid2];
      xright2 = xs2; }
    else
    { xleft2  = xs2;
      xright2 = x[imid2]; }
    if ( x[imid1] < xs1 )    
    { xleft1  = x[imid1];
      xright1 = xs1; }
    else
    { xleft1  = xs1;
      xright1 = x[imid1]; }
  }
/*
#ifdef DEBUGGING
  printf("top:   (%f,%f)\n", x[ihi],y[ihi] );
  printf("base2: (%f,%f) -- (%f,%f)\n", xleft2,y[imid2], xright2,y[imid2] );
  printf("base1: (%f,%f) -- (%f,%f)\n", xleft1,y[imid1], xright1,y[imid1] );
  printf("bot:   (%f,%f)\n", x[ilo],y[ilo] );

  printf("\n do_lower, do_mid, do_upper: %d%d%d\n", do_lower, do_mid, do_upper );
#endif
*/

  /* 
   * step 3. Process each pixel in the quadrangle.
   */

  jmin = ROUNDUP(   (y[ilo]-ymin)/dy - SMALL );
  jmax = ROUNDDOWN( (y[ihi]-ymin)/dy + SMALL );       
  
  if (jmin<0)          jmin = 0;
  if (jmax>im->vres-1) jmax = im->vres-1;       /* domain clipping */
  
  for (j=jmin; j<=jmax; j++)
  { 
    yy = ymin + j*dy;
    if ( do_upper && (yy>y[imid2] || !do_mid && !do_lower) )
    { ratio = (yy-y[imid2]) / (y[ihi]-y[imid2]);
      x1 = xleft2  + (x[ihi]-xleft2)*ratio;
      x2 = xright2 + (x[ihi]-xright2)*ratio; }
    else if ( do_mid && (yy>y[imid1] || !do_lower) )
    { ratio = (yy-y[imid1]) / (y[imid2]-y[imid1]);
      x1 = xleft1  + (xleft2-xleft1)*ratio;
      x2 = xright1 + (xright2-xright1)*ratio; }
    else  
    { ratio = (yy-y[ilo]) / (y[imid1]-y[ilo]);
      x1 = x[ilo] + (xleft1-x[ilo])*ratio;
      x2 = x[ilo] + (xright1-x[ilo])*ratio; }
  
  
    imin = ROUNDUP(   (x1-xmin)/dx - SMALL );
    imax = ROUNDDOWN( (x2-xmin)/dx + SMALL );   /* expand slightly, otherwise some pixels may be lost */
  
    if (imin<0)          imin = 0;
    if (imax>im->hres-1) imax = im->hres-1;     /* domain clipping */

    for (i=imin; i<=imax; i++)
    { 
      xx = xmin + i*dx;

      /* Calculate iso-parametric coords (xi,eta), implicitly defined by:  
       *     xx = x[0] + xi*(x[1]-x[0]) + eta*(x[3]-x[0]) 
       *                                + xi*eta*(x[2]+x[0]-x[1]-x[3])
       *     yy = y[0] + xi*(y[1]-y[0]) + eta*(y[3]-y[0])
       *                                + xi*eta*(y[2]+y[0]-y[1]-y[3])
       * This mapping maps the quadrangle to the unit square in xi,eta space:
       *     (x[0],y[0]) ---> (0,0)
       *     (x[1],y[1]) ---> (1,0)
       *     (x[2],y[2]) ---> (1,1)
       *     (x[3],y[3]) ---> (0,1)
       * For each pixel, the "physical coords" xx,yy are known.
       * The corresponding xi,eta must be found from a Newton process, since 
       * the isoparametric mapping is nonlinear.
       */

      xi   = prev_xi  = 0.5;
      eta  = prev_eta = 0.5;   /* initial guess for Newton iteration */
               /* If nr of pixels for this quadrangle is large, 
                  a more efficient guess is to use "previous" values */
      iter = 0;

      do {
        iter++;

        /* the current residual */
        rhs1 = RESX(xi,eta);
        rhs2 = RESY(xi,eta);

        /* the Jacobian of the mapping */
        a11 = x[1] - x[0] + eta*(x[2]+x[0]-x[1]-x[3]);
        a12 = x[3] - x[0] +  xi*(x[2]+x[0]-x[1]-x[3]);
        a21 = y[1] - y[0] + eta*(y[2]+y[0]-y[1]-y[3]);
        a22 = y[3] - y[0] +  xi*(y[2]+y[0]-y[1]-y[3]);
        det = a11*a22 - a12*a21;

        xi  = xi  + (  a22*rhs1 - a12*rhs2 ) / det;
        eta = eta + ( -a21*rhs1 + a11*rhs2 ) / det;

/*      
 *       This robust convergence criterium is too complex for Macintosh MPW C:
 *
 *       noconverg = ( RESX(xi+EPS,eta)*RESX(xi-EPS,eta) > 0 ||
 *                     RESX(xi,eta+EPS)*RESX(xi,eta-EPS) > 0 ||
 *                     RESY(xi+EPS,eta)*RESY(xi-EPS,eta) > 0 ||
 *                     RESY(xi,eta+EPS)*RESY(xi,eta-EPS) > 0 );   
 */

        noconverg = ( fabs(xi-prev_xi) + fabs(eta-prev_eta) > EPS );
        prev_xi  = xi;
        prev_eta = eta;
      } 
      while( noconverg && iter<MAXIT && fabs(xi)+fabs(eta) < DIVERG );


      if (noconverg)
      { 
        printf("iso_quadrangle: Newton failed; switching to triangularization.\n");
        quadrangle( x, y, dat, im );  /* use more robust approach instead */
        return;
      }

      if ( xi>-SMALL && eta>-SMALL && xi<1.0+SMALL && eta<1.0+SMALL )  /* interior point */
      { 
        dd = dat[0] + xi*(dat[1]-dat[0]) + eta*(dat[3]-dat[0]) + 
             xi*eta*(dat[2]+dat[0]-dat[1]-dat[3]);
        if (dd>datamax) dd=datamax;
        if (dd<datamin) dd=datamin;             /* data clipping */
        IMAGE(i,j) = 1.5 + (dd-datamin)/(datamax-datamin)*237.9;
      }
#ifdef DEBUGGING
      else
        printf("iso_quadrangle: pixel rejected.\n");
      kount[iter]++;
#endif
    }
  }
}



