#define VERSION "2.1" /* don't touch this */ #define VDATE "December 13, 1992" /* --------------------------------------------------------------------------- Purpose: --------------------------------------------------------------------------- vs2ris produces 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 | | color raster (single HDF Vset file ---> | vs2ris | ------> (HDF RIS8 file) or three HDF files) | | ----------- Presently, only float32 data and coordinates are supported. Connectivity lists must have the type int16 or int32. --------------------------------------------------------------------------- Source Availability: --------------------------------------------------------------------------- This software was developed and made available by: Fred Walsteijn Institute for Marine and Atmospheric Research Utrecht University 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- : --------------------------------------------------------------------------- VSsizeof() is used to distinguish "int16" from "int32" connectivity lists. The total vdata element size returned in arg vsize by VSinquire() includes the size of other fields if present. Therefore the conn list element size cannot be obtained from VSinquire() in general. VSsizeof() is also used to verify that only float32 input is present. --------------------------------------------------------------------------- Implementation Notes -3- : --------------------------------------------------------------------------- Here is a list of preprocessor constants which modify the behaviour of vs2ris. Use them in the C compile command/Makefile if necessary: -DNO_VSET --------- If NO_VSET is defined the Vset functions are not accessed: then only curvilinear data can be used as input. Useful on systems where the Vset library functions are not implemented but the regular HDF functions are. -DDEBUG ----------- Write small amount of info to stdout: command line args, read_vset() stuff, Newton iteration count. You don't need this option. -DDEBUGGING ----------- Writes a huge amount of uninteresting info to stdout: details of the pixel generating function(s) for *each* element. Don't enable it. -DMAC, -DMPW ------------ On the Mac, vs2ris does some handy things: sets the output file type/creator to that of NCSA Image; rotates the cursor. --------------------------------------------------------------------------- Speed comparisons: --------------------------------------------------------------------------- Benchmark performed on a MacIIsi with Apple's NuBus/coprocessor card and built-in colors disabled. Test cases with rectangular input data are included to compare "vs2ris" with the standard NCSA HDF tool "fp2hdf". Both tools were 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. 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. Use something better than VSsizeof() to determine the integral type of the connectivity lists in read_vset(). 3. Include support for float64 scales, data, ...; and other ints too ? Maybe "void *" vars instead of "float32 *", "int16 *", "int32 *" etc; cast to correct type when actually used in calculations ? --------------------------------------------------------------------------- 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. If you encounter a bug, report it to: Fred Walsteijn walsteyn@fys.ruu.nl --------------------------------------------------------------------------- History: --------------------------------------------------------------------------- Jun 9, 1991 -- first beta version Sep 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 Apr 20, 1992 -- every strcpy() replaced by pointer assignment: fixes problems with long pathnames/args. FNAME_MAX removed; STR_MAX replaced by "vg.h" constants. Mac cursor now rotates. NO_VSET introduced. Sep 24, 1992 -- vs2ris now interfaces to NCSA's HDF 3.2r1 library. Dec 13, 1992 -- vs2ris now interfaces to NCSA's HDF 3.2r2 library; read_vset(): uses Vinitialize(), Vfinish(), DFACC_READ; checks float size with VSsizeof(); read_crv(): checks float size with DFSDgetNT(). (by: Fred Walsteijn, walsteyn@fys.ruu.nl) */ /* ---------------------------- includes --------------------------- */ #include #include #include #include #include #ifdef NO_VSET #include "hdf.h" #else #include "vg.h" #endif #include "dfsd.h" #if ( defined(MAC) && defined(MPW) ) #include /* InitCursorCtl(), RotateCursor() */ #endif /* ---------------------------- 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: float32 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 uint8 /* image pixel type; "df.h": uint8 == unsigned char */ #define PIXELPTR VOIDP /* pointer to image pixel type; needed in DFR8putimage() call */ #define PAL uint8 /* palette entry type */ #define PALPTR VOIDP /* pointer to palette entry type; needed in DFPgetpal() call */ #define BACKGR_COL 0 /* background color in image */ #define CMAX 10 /* maximum number of connectivity lists to look for */ #define MIN_RES 3 /* minimum number of pixels in x & y direction */ #undef FALSE #undef TRUE #undef CRV #undef VSET #undef INT16TYPE #undef INT32TYPE enum boolean { FALSE, TRUE }; enum filetype { CRV, VSET }; /* input file(s) type: CRV == 3*HDF, VSET == 1 HDF Vset */ enum elmtype { INT16TYPE, INT32TYPE }; /* integral type of HDF Vset connectivity lists */ /* * 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 */ *y, /* vertical coords of vertices */ *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 == INT16TYPE then triangles in splist3; otherwise plist3 */ etyp4; /* etyp4 == INT16TYPE then quadrangles in splist4; otherwise plist4 */ int32 *plist3, *plist4; /* pointers to int32 connectivity lists */ int16 *splist3, *splist4; /* pointers to int16 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, /* input curvi-scalar HDF file name || input HDF Vset file name */ *gridx, *gridy, /* input curvi-coordinate HDF file names, CRV only */ *outfile, /* output HDF RIS8 file name */ *palfile; /* palette HDF 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, VSET only */ int32 ncon; /* number of connectivity lists specified, VSET only */ char *groupname, *pxname, *pyname, *connectname[CMAX], *scalarname; /* vgroup name and vdata field 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.2r2 prototypes // function result --------------------------------------------------------------------------- In "hdf.h": #define SUCCEED 0 #define FAIL (-1) In "vg.h": typedef int32 HFILEID; --------------------------------------------------------------------------- hproto.h <- hdf.h <- vg.h --------------------------------------------------------------------------- int32 Hopen (char *path, intn access, int16 ndds); // FAIL, file_id intn Hclose (int32 file_id); // FAIL, SUCCEED intn DFPgetpal (char *filename, VOIDP palette); // FAIL, SUCCEED int DFR8setpalette (uint8 *pal); //always SUCCEED int DFR8putimage (char *filename, VOIDP image, int32 xdim, int32 ydim, uint16 compress); // FAIL, SUCCEED --------------------------------------------------------------------------- dfsd.h --------------------------------------------------------------------------- int DFSDgetdims (char *filename, intn *prank, int32 sizes[], intn maxrank); // FAIL, SUCCEED int DFSDgetrange (void *pmax, void *pmin); // FAIL || -1, 0 int DFSDgetdata (char *filename, intn rank, int32 maxsizes[], void *data); // FAIL, SUCCEED int DFSDgetNT(int32 *pnumbertype); // FAIL, 0 --------------------------------------------------------------------------- *VSET* vproto.h <- vg.h --------------------------------------------------------------------------- void Vinitialize(HFILEID f); void Vfinish(HFILEID f); VGROUP *Vattach (HFILEID f, int32 vgid, char *accesstype); // NULL, ptr void Vdetach (VGROUP *vg); void Vgetname (VGROUP *vg, char *vgname); int32 Visvs (VGROUP *vg, int32 id); // 0, 1 (=="id refers to vdata") int32 Vgetid (HFILEID f, int32 vgid); // FAIL, >=0 (==id of next vgroup) int32 Vgetnext (VGROUP *vg, int32 id); // FAIL, >=0 (==id of next entry) VDATA * VSattach (HFILEID f, int32 vsid, char *accesstype); // NULL, ptr void VSdetach (VDATA *vs); int32 VSread (VDATA *vs, BYTE buf[], int32 nelt, int32 interlace); // FAIL, >=0 (==nr elems read) int32 VSfexist (VDATA *vs, char *fields); // FAIL, 1 int32 VSinquire (VDATA *vs,int32 *nelt, int32 *interlace, char *fields, int32 *eltsize, char *vsname); // FAIL, SUCCEED int32 VSsetfields (VDATA *vs, char *fields); // FAIL, SUCCEED int32 VSsizeof (VDATA *vs, char *fields); // FAIL, >0 (==fieldssize) --------------------------------------------------------------------------- */ /********************************************************************************* * 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 */ #if ( defined(MAC) && defined(MPW) ) InitCursorCtl(NULL); RotateCursor(32); #endif 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 DEBUG write_opt( opt ); #endif setup_image( &im ); /* initialize 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: Internal error: unknown 'in.ftyp' value."); write_image( opt.outfile, im ); #ifdef DEBUG 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 = INT16TYPE; in->dim0 = in->dim1 = in->ntri = in->nquad = in->nv = 0; im->hres = im->vres = 0; im->image = NULL; opt->infile = opt->gridx = opt->gridy = opt->outfile = opt->palfile = ""; } /**************************************************************************** * 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 float32 data\n"); printf( "or 'curvilinear' HDF float32 data.\n"); printf( "By: Fred Walsteijn, Utrecht University, The Netherlands (walsteyn@fys.ruu.nl)\n\n"); printf( "SYNTAX:\n" ); printf( " %s -vs \n", tool ); printf( " -o -res [OPTIONS ...] [VS-OPTIONS ...]\n\n" ); printf( "OR:\n" ); printf( " %s -crv \n", tool ); printf( " -o -res [OPTIONS ...]\n\n" ); printf( "where:\n" ); printf( " -vs \n" ); printf( " Input HDF Vset file containing 2-D vertex grid & data.\n" ); printf( " -crv \n" ); printf( " Three HDF input files specifying data on a curvilinear 2-D grid.\n" ); printf( " Each file contains a single SDS describing respectively the data values,\n" ); printf( " x coordinates, and y coordinates of the gridpoints.\n" ); printf( " -o \n" ); printf( " HDF RIS8 output file containing a raster image.\n" ); printf( " -res \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 \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's. Default: min/max of coordinates in \n" ); printf( " input file(s).\n" ); printf( " -range \n" ); printf( " The palette is scaled on this data range. Allows clipping.\n" ); printf( " Default: min/max of input file data.\n" ); printf( " -p \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 \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 \n" ); printf( " Name of vdata field containing x coords. Default: 'px'.\n" ); printf( " -py \n" ); printf( " Name of vdata field containing y coords. Default: 'py'.\n" ); printf( " { -connect } ...\n" ); printf( " Name(s) of vdata field(s) containing a triangle and/or quadrangle\n" ); printf( " connectivity list of type int16 or int32. 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 \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 name and vdata fieldnames similar (not identical) to PolyView defaults: */ opt->groupname = ""; /* "the 1st vgroup" */ opt->pxname = "px"; opt->pyname = "py"; opt->scalarname = "scalar"; opt->connectname[0] = "plist"; opt->connectname[1] = "plist3"; 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; iinfile = argv[++i]; opt->gridx = argv[++i]; opt->gridy = argv[++i]; } else if ( !strcmp(argv[i],"-vs") && iinfile = argv[++i]; } else if ( !strcmp(argv[i],"-o") && ioutfile = argv[++i]; } else if ( !strcmp(argv[i],"-iso") ) iso = TRUE; else if ( !strcmp(argv[i],"-isot") ) isot = TRUE; else if ( !strcmp(argv[i],"-domain") && idomain = 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") && iminmax = 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") && ihres = hres; im->vres = vres; } else if ( !strcmp(argv[i],"-p") && ipal = TRUE; opt->palfile = argv[++i]; } else if (!strcmp(argv[i],"-group") && igroupname = argv[++i]; else if (!strcmp(argv[i],"-px") && ipxname = argv[++i]; else if (!strcmp(argv[i],"-py") && ipyname = argv[++i]; else if (!strcmp(argv[i],"-scalar") && iscalarname = argv[++i]; else if (!strcmp(argv[i],"-connect") && incon < CMAX ) { opt->connect = TRUE; 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 DEBUG /************************************************************ * 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==INT16TYPE)?"INT16TYPE":((arg==INT32TYPE)?"INT32TYPE":"*TRASH*") ) void write_opt( struct Options opt ) { int32 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; iinfile, &rank, hdfdims, 2) ) quit("read_crv: Unable to get dimensions from data file."); if (DFSDgetNT(&numtype)==FAIL || numtype!=DFNT_FLOAT32) quit("read_crv: Sorry, I can only read float32 scalar-data."); 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 && !DFSDgetrange(&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 (DFSDgetNT(&numtype)==FAIL || numtype!=DFNT_FLOAT32) quit("read_crv: Sorry, I can only read float32 gridx-data."); if ( DFSDgetdims(opt->gridy, &rank, hdfdims, 2) || rank != 2 || hdfdims[0] != in->dim0 || hdfdims[1] != in->dim1 ) quit("read_crv: Incompatible gridy file."); if (DFSDgetNT(&numtype)==FAIL || numtype!=DFNT_FLOAT32) quit("read_crv: Sorry, I can only read float32 gridy-data."); 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. * ****************************************************/ #ifdef NO_VSET void read_vset( struct Input *in, struct Options *opt ) { quit("read_vset: Sorry, your system does not support Vset files."); } #else /* !NO_VSET */ void read_vset( struct Input *in, struct Options *opt ) { VGROUP *vg; VDATA *vs; int32 f; char vgname[VGNAMELENMAX+1], vsname[VSNAMELENMAX+1], fields[VSFIELDMAX*(FIELDNAMELENMAX+1)]; int32 vgid, vsid; enum boolean found, pxfound, pyfound, scalarfound, connectfound[CMAX], confound, oneconfound, con3found, con4found; int32 i, csize, vsize, interlace, nelements, npx, npy, nscalar; if ( (f = Hopen( opt->infile, DFACC_READ, 0 )) == FAIL) quit("read_vset: Cannot open input file."); Vinitialize(f); /* * 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)) != FAIL) { 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 0 ) { #ifdef DEBUG 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) == SUCCEED && (csize = VSsizeof(vs, opt->pxname)) > 0) { if ( csize != sizeof(float32) ) quit("read_vset: Sorry, I can only read float32 px-data."); if ((in->x = (float32 *) malloc((size_t) (nelements*sizeof(float32)))) == NULL) quit("read_vset: Cannot allocate mem for px."); VSread(vs, (BYTE *)(in->x), nelements, NO_INTERLACE ); pxfound = TRUE; npx = nelements; #ifdef DEBUG printf(" px read\n" ); #endif } if ( pyfound == FALSE && VSfexist(vs, opt->pyname) == 1 && VSsetfields(vs, opt->pyname) == SUCCEED && (csize = VSsizeof(vs, opt->pyname)) > 0) { if ( csize != sizeof(float32) ) quit("read_vset: Sorry, I can only read float32 py-data."); if ((in->y = (float32 *) malloc((size_t) (nelements*sizeof(float32)))) == NULL) quit("read_vset: Cannot allocate mem for py."); VSread(vs, (BYTE *)(in->y), nelements, NO_INTERLACE ); pyfound = TRUE; npy = nelements; #ifdef DEBUG printf(" py read\n" ); #endif } if ( scalarfound == FALSE && VSfexist(vs, opt->scalarname) == 1 && VSsetfields(vs, opt->scalarname) == SUCCEED && (csize = VSsizeof(vs, opt->scalarname)) > 0) { if ( csize != sizeof(float32) ) quit("read_vset: Sorry, I can only read float32 scalar-data."); if ((in->data = (float32 *) malloc((size_t) (nelements*sizeof(float32)))) == NULL) quit("read_vset: Cannot allocate mem for scalar."); VSread(vs, (BYTE *)(in->data), nelements, NO_INTERLACE ); scalarfound = TRUE; nscalar = nelements; #ifdef DEBUG 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]) == SUCCEED && (csize = VSsizeof(vs, opt->connectname[i])) > 0) { #ifdef DEBUG printf(" attempting conn list named %s\n", opt->connectname[i] ); printf(" size found %d\n", csize ); #endif if ( con3found == FALSE && csize == 3*sizeof(int16) ) /* int16 triangle conn list */ { if ((in->splist3 = (int16 *) malloc((size_t) (nelements*csize))) == NULL) quit("read_vset: Cannot allocate mem for connectivity list."); VSread(vs, (BYTE *)(in->splist3), nelements, NO_INTERLACE ); in->etyp3 = INT16TYPE; in->ntri = nelements; con3found = TRUE; #ifdef DEBUG printf(" it's a tri-INT16TYPE\n" ); #endif } else if ( con4found == FALSE && csize == 4*sizeof(int16) ) /* int16 quadrangle conn list */ { if ((in->splist4 = (int16 *) malloc((size_t) (nelements*csize))) == NULL) quit("read_vset: Cannot allocate mem for connectivity list."); VSread(vs, (BYTE *)(in->splist4), nelements, NO_INTERLACE ); in->etyp4 = INT16TYPE; in->nquad = nelements; con4found = TRUE; #ifdef DEBUG printf(" it's a quad-INT16TYPE\n" ); #endif } else if ( con3found == FALSE && csize == 3*sizeof(int32) ) /* int32 triangle conn list */ { if ((in->plist3 = (int32 *) malloc((size_t) (nelements*csize))) == NULL) quit("read_vset: Cannot allocate mem for connectivity list."); VSread(vs, (BYTE *)(in->plist3), nelements, NO_INTERLACE ); in->etyp3 = INT32TYPE; in->ntri = nelements; con3found = TRUE; #ifdef DEBUG printf(" it's a tri-INT32TYPE\n" ); #endif } else if ( con4found == FALSE && csize == 4*sizeof(int32) ) /* int32 quadrangle conn list */ { if ((in->plist4 = (int32 *) malloc((size_t) (nelements*csize))) == NULL) quit("read_vset: Cannot allocate mem for connectivity list."); VSread(vs, (BYTE *)(in->plist4), nelements, NO_INTERLACE ); in->etyp4 = INT32TYPE; in->nquad = nelements; con4found = TRUE; #ifdef DEBUG printf(" it's a quad-INT32TYPE\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); Vfinish(f); Hclose(f); if ( pxfound == FALSE ) quit("read_vset: Cannot find px field."); if ( pyfound == FALSE ) quit("read_vset: Cannot find py field."); if ( scalarfound == FALSE ) quit("read_vset: Cannot find scalar field."); if ( oneconfound == FALSE ) quit("read_vset: Cannot find any connect field."); 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."); } #endif /* !NO_VSET */ /***************************************************************** * 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, (PIXELPTR) 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.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 ); #if ( defined(MAC) && defined(MPW) ) if ( i % cursor3 == 0 ) RotateCursor(32); #endif } /* * Finally the quadrangles. */ for (i=0; i < in.nquad; i++) { for (j=0; j<4; j++) { ptr = (in.etyp4==INT16TYPE)? 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 ); } #if ( defined(MAC) && defined(MPW) ) if ( i % cursor4 == 0 ) RotateCursor(32); #endif } } /*************************************************************** * 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]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 (ddhres-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 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 triangulation.\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-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 (ddy[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-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