/* Mex version of Roger Lukas's ZGRID routine, translated into C using f2c, minimally modified (including switch to double precision), and wrapped in a mex gateway. The modifications were done in part to eliminate the need to link with any f2c libraries. Some of the diagnostic output was deleted to make this easier. This mex function is intended to be called by an m-file, not directly by the user. It omits most input parameter error checking, so be careful what you feed it! Eric Firing Sat 97/04/26 */ #include "mex.h" #include "matrix.h" #include #include #include /* Defines that substitute mex versions of C library functions. */ #define calloc mxCalloc #define printf mexPrintf #define free mxFree /* Defines for f2c types, converting everything to double precision */ #define integer int #define doublereal double #define float double #define real double /* Defines for f2c functions. */ #define dabs(x) fabs(x) #define dmax(x,y) ((x)>(y)?(x):(y)) /* Global declarations. */ /* common block variables were turned into structures of pointers */ /* Input data */ typedef struct { double *xp, *yp, *zp; int n; } data_1_TYPE; /* Blanking polygon, which zgrid uses also as the output array. */ typedef struct { double *z; int nrow; int nr1; /* nrow + 1, for convenience in indexing */ int ncol; } grid_TYPE; int zgrid_(integer *nx, integer *ny, real *x1, real *y1, real *dx, real *dy, integer *n, real *del, real *cay, integer *nrng, integer *idbug, grid_TYPE grid, data_1_TYPE data_1); /* Mex inputs: z, xp, yp, zp, x1, y1, dx, dy, del, cay, nrng, idbug 0 1 2 3 4 5 6 7 8 9 10 11 */ /*----------------------------------------------------*/ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray* prhs[]) { int nrng, idbug; double x1, y1, dx, dy; int nr, nc, nn2, nn3, ngrid; double del, cay; int ret; int i; double *zb; double NaN; data_1_TYPE data_1; grid_TYPE grid; if (nlhs != 1) mexErrMsgTxt("Need 1 output: gridded z."); if (nrhs != 12) mexErrMsgTxt("Need 12 inputs to zgrid mex."); NaN = mxGetNaN(); zb = mxGetPr(prhs[0]); grid.nrow = mxGetM(prhs[0]); grid.ncol = mxGetN(prhs[0]); ngrid = grid.nrow*grid.ncol; grid.z = mxCalloc(ngrid, sizeof(double)); for (i = 0; i < ngrid; i++) { if (mxIsNaN(zb[i])) grid.z[i] = 1E35; else grid.z[i] = zb[i]; } data_1.xp = mxGetPr(prhs[1]); data_1.yp = mxGetPr(prhs[2]); data_1.zp = mxGetPr(prhs[3]); data_1.n = mxGetM(prhs[1]) * mxGetN(prhs[1]); nn2 = mxGetM(prhs[2]) * mxGetN(prhs[2]); nn3 = mxGetM(prhs[3]) * mxGetN(prhs[3]); if (nn2 != nn3 || nn2 != data_1.n) mexErrMsgTxt("Lengths of xp, yp, zp don't match."); x1 = mxGetScalar(prhs[4]); y1 = mxGetScalar(prhs[5]); dx = mxGetScalar(prhs[6]); dy = mxGetScalar(prhs[7]); del = mxGetScalar(prhs[8]); cay = mxGetScalar(prhs[9]); nrng = (int) mxGetScalar(prhs[10]); idbug = (int) mxGetScalar(prhs[11]); plhs[0] = mxCreateDoubleMatrix(grid.nrow, grid.ncol, mxREAL); ret = zgrid_(&grid.nrow, &grid.ncol, &x1, &y1, &dx, &dy, &data_1.n, &del, &cay, &nrng, &idbug, grid, data_1); if (ret != 0) mexErrMsgTxt("zgrid failed"); for (i = 0; i < ngrid; i++) { if (grid.z[i] >= 9e34 ) grid.z[i] = NaN; } mxSetPr(plhs[0], grid.z); } /* mexFunction */ int zgrid_(integer *nx, integer *ny, real *x1, real *y1, real *dx, real *dy, integer *n, real *del, real *cay, integer *nrng, integer *idbug, grid_TYPE grid, data_1_TYPE data_1) { int *jmnew; /* Initialized data */ static real big = 9e34; static real eps = 0.001; static integer itmax = 100; static integer nopt = 30; static integer lp = 16; /* System generated locals */ integer i__1, i__2, i__3; real r__1; double sqrt(doublereal); /* Local variables */ static real zadd; static integer iter, nnew; static real zijn, zimm, zmin, zjmm, zmax, zipp, zjpp, zsum, zul200; static integer i, j, k, idivg; static real zbase, pcent, relax; static integer imnew; static real dzmax, dzrms, rtrms, caysd4; static integer im, jm; static real z00, dz; static integer nblank; static real cayscl, corect, zrange, epslon; static integer ngrpts; static real dzrmsp, zulovf, den, abz; static integer nij, npg; static real zij, zim, zjm, wgt; static integer nxm, nym; static real zpk, zip, zjp, zul, del2, del4; /* ********************************************************************** */ /* PGM = ZGRID */ /* ROGER B. LUKAS */ /* DEPT. OF OCEANOGRAPHY */ /* UNIV. OF HAWAII */ /* 12/2/80 */ /* ********************************************************************** */ /* ZGRID interpolates irregularly distributed data to a */ /* rectangular grid. The interpolation is controlled by */ /* variation of the grid mesh and the parameters CAY and DEL. */ /* The method used is to first assign each data point to */ /* the closest grid point. If more than one data point is */ /* assigned to a grid point, the data are averaged at that */ /* grid point. Presently, ZGRID can handle up to 200 data */ /* points at one grid point. Any data points after 200 have */ /* been assigned to a grid point are ignored for that grid */ /* point. A grid point can be "blanked out" by */ /* assigning a value of 1.0E35 or greater to it prior to */ /* entry to ZGRID. */ /* All unblanked, unfilled grid points are interpolated */ /* from the filled grid points as the solution to the */ /* partial differential equation */ /* (D2Z/DX2) + (D2Z/DY2) - K ( (D4Z/DX4) + (D4Z/DY4) ) = 0. */ /* The solution to this equation is found by using the */ /* method of successive over-relaxation (See Numerical */ /* Methods for Partial Differential Equations by William */ /* Ames). */ /* An initial estimate of the solution is made by assigning */ /* each unfilled grid point the value of the closest */ /* filled grid point, as long as there is no more than */ /* NRNG grid points between them (measured along the axes */ /* directions only, no diagonal distances). If a grid */ /* point is still unfilled this process, it is set */ /* to 1.0E35 and blanked out. */ /* initially, Gauss-Seidel iteration is carried out */ /* until enough information about the rate of convergence */ /* to a solution can be obtained. Then an attempt is */ /* made to estimate the optimum relaxation parameter */ /* to the convergence. The process of refining */ /* the estimate the optimum relaxation parameter */ /* is continued every NOPT iterations. */ /* If the solution does not converge within ITMAX iterations, */ /* the user is warned, and the most recent grid is */ /* returned. This can happen because either */ /* 1) there are large areas of the grid */ /* without any data, */ /* 2) the regions without data are not */ /* "simply" connected, or */ /* 3) or the boundary conditions are inconsistent. */ /* There are two types of boundaries: those where the */ /* points outside the boundary are undefined (the edges */ /* of the grid and interior blanked regions), and those */ /* where the points on and outside the boundary are filled */ /* with data. For the first type, the boundary conditions */ /* are that DZ/DN = D2Z/DN2 = 0. For the second type, */ /* the boundary conditions are Z = ZB and D2Z/DN2 = 0. */ /* here D/DN is the partial derivative normal to the */ /* boundary. */ /* Subroutine arguments are: */ /* Z The grid to which the data is to interpolated. */ /* NX The number of grid points along the abcissa. */ /* NY The number of grid points along the ordinate. */ /* X1 The abcissa value of Z(1,1) in data units */ /* Y1 The ordinate value of Z(1,1) in data units */ /* DX,DY The length of the sides of the grid cells */ /* IN DATA UNITS */ /* N The number of data points to be gridded. */ /* CAY The parameter controlling how much the */ /* fourth order part of the differential equation */ /* is to be used in the interpolation. */ /* DEL The parameter controlling the aspect ratio of */ /* the interpolation. DEL is defined as DX/DY */ /* strictly, but has been made a parameter to */ /* allow the user to exercise more control over */ /* the interpolation. */ /* NRNG The maximum distance in grid points which will be */ /* used to search for valid data as a first estimate to */ /* the solution. */ /* IDBUG Controls the level of debugging information */ /* that is printed out. 0 gives none, 1 records */ /* entry and exit from the subroutine, and 2 records */ /* the iteration process and lists the final grid */ /* values. */ /* JMNEW should be dimensioned the same as the row dimension */ /* of Z. */ /* *********************************************************************** */ grid.nr1 = grid.nrow + 1; if (*idbug > 0) { printf("nx %d, ny %d, n %d, cay %f, nrng %d\n", *nx, *ny, *n, *cay, *nrng); } /* CHECK FOR EXISTENCE OF DATA */ /* IF NO DATA, ISSUE ERROR MESSAGE AND TERMINATE */ L10: if (*n == 0) { printf("No data\n"); return(-1); } /* COMPUTE LOCAL CONSTANTS FROM GLOBAL ARGUMENTS */ nxm = *nx - 1; nym = *ny - 1; del2 = *del * *del; del4 = del2 * del2; cayscl = *cay / (*dx * *dx); caysd4 = cayscl * del4; /* ZERO GRID POINTS THAT ARE NOT BLANKED OUT */ ngrpts = *nx * *ny; nblank = 0; i__1 = *ny; for (j = 1; j <= i__1; ++j) { i__2 = *nx; for (i = 1; i <= i__2; ++i) { if (grid.z[i + j * grid.nrow - grid.nr1] < big) { goto L39; } ++nblank; goto L40; L39: grid.z[i + j * grid.nrow - grid.nr1] = 0.0; L40: ; } /* L41: */ } /* CHECK TO SEE IF THE WHOLE GRID IS BLANKED OUT */ if (nblank == ngrpts) { printf("Everything is blanked out.\n"); return(-1); } /* GET ZBASE WHICH WILL MAKE ALL ZP VALUES POSITIVE BY AT LEAST */ /* .25(ZMAX-ZMIN) AND FILL IN GRID WITH ZEROS. */ /* *********************************************************************** */ /* FIND MAXIMUM AND MINIMUM FOR Z */ zmin = big; zmax = -(doublereal)big; i__1 = *n; for (k = 1; k <= i__1; ++k) { zpk = data_1.zp[k - 1]; if (zpk > zmax) { zmax = zpk; } if (zpk < zmin) { zmin = zpk; } /* L20: */ } /* DETERMINE RANGE OF Z VALUES */ /* SET UP METHOD TO SET RANGES */ zrange = zmax - zmin; /* ABEND IF ZRANGE IS NOT POSITIVE */ if (zrange <= 0.0) { printf("Range of input values = 0; check input.\n"); return(-1); } zbase = -(doublereal)zmin + zrange * 0.25; zul = zbase + zmax; zul200 = zul * 200.0; zulovf = zul200 * 200.0; /* AFFIX EACH POINT ZP TO NEAREST GRID PT. TAKE AVG IF MORE THAN */ /* ONE NEAR PT. ADD ZBASE PLUS 10*ZRANGE AND MAKE NEGATIVE. */ /* INITIALLY SET EACH UNSET GRID PT TO VALUE OF NEAREST KNOWN PT */ /* ZUL200 IS USED TO KEEP TRACK OF # OF POINTS ASSIGNED TO EACH */ /* GRID POINT. ZULOVF CHECKS TO MAKE SURE AVERAGE WILL BE CORRECT */ /* BY REJECTING ANY NEW POINTS. PRESENTLY SET UP SO THAT NO MORE */ /* THAN 200 POINTS CAN BE ASSIGNED TO ONE GRID POINT. */ /* *********************************************************************** */ zadd = zbase + zul200; i__1 = *n; for (k = 1; k <= i__1; ++k) { i = (data_1.xp[k - 1] - *x1) / *dx + 1.5; if (i <= 0 || i > *nx) { goto L110; } j = (data_1.yp[k - 1] - *y1) / *dy + 1.5; if (j <= 0 || j > *ny) { goto L110; } if (grid.z[i + j * grid.nrow - grid.nr1] >= zulovf) { goto L110; } grid.z[i + j * grid.nrow - grid.nr1] = grid.z[i + j * grid.nrow - grid.nr1] + data_1.zp[ k - 1] + zadd; L110: ; } /* AVERAGE DATA ASSIGNED TO EACH GRID POINT. */ /* ALSO, SEARCH GRID FOR GRID POINTS WITH NO DATA. */ /* NPG AT END OF SEARCH IS THE NUMBER OF GRID POINTS */ /* WITHOUT DATA. */ npg = 0; zadd = zrange * 10.0 - zul200; i__1 = *ny; for (j = 1; j <= i__1; ++j) { i__2 = *nx; for (i = 1; i <= i__2; ++i) { zij = grid.z[i + j * grid.nrow - grid.nr1]; /* SKIP BLANKED POINTS */ if (zij >= big) { goto L150; } nij = zij / zul200; if (nij <= 0) { goto L145; } grid.z[i + j * grid.nrow - grid.nr1] = -(doublereal)(zij / nij + zadd); goto L150; /* NO DATA AT THIS GRID POINT, FLAG IT. */ L145: grid.z[i + j * grid.nrow - grid.nr1] = -1e35f; ++npg; L150: ; } /* L151: */ } /* IF EVERY GRID POINT HAS AT LEAST ONE DATA POINT */ /* ASSIGNED TO IT, THEN THERE IS NO NEED TO DO ANY */ /* INTERPOLATION, SO RETURN. */ if (npg == 0) { printf("Every grid point has a data point; no interpolation needed.\n"); goto L2120; } /* IF 5% OR LESS OF THE UNBLANKED GRID POINTS HAVE VALID DATA, THEN */ /* TERMINATE EXECUTION SINCE THE INTERPOLATION WILL TAKE A LONG TIME */ /* TO CONVERGE TO A SOLUTION, IF AT ALL. */ L155: pcent = (real) (npg / (ngrpts - nblank)); if (pcent >= 0.98) { printf("Less than 2%% of the grid has data; quitting.\n"); return(-1); } /* FOR EACH GRID POINT WITHOUT DATA, ASSIGN IT THE */ /* VALUE OF THE NEAREST GRID POINT AS A FIRST */ /* APPROXIMATION TO THE SOLUTION, BUT DON'T EXTRAPOLATE MORE */ /* THAN NRNG GRID CELLS. ANY GRID POINTS STILL UNASSIGNED */ /* AFTER EXTRAPOLATION WILL BE BLANKED OUT. */ /* IMNEW AND JMNEW KEEP TRACK OF EXTRAPOLATED POINTS */ /* jmnew is indexed by i, the row-dimension index, so its length is nx, or grid.nrow. */ if (*nrng <= 0) { goto L200; /* Skip this step entirely. */ } jmnew = calloc(grid.nrow,sizeof(int)); if (jmnew == NULL) { printf("Failed to allocate memory for jmnew.\n"); return (-1); } i__1 = *nrng; for (iter = 1; iter <= i__1; ++iter) { nnew = 0; i__2 = *ny; for (j = 1; j <= i__2; ++j) { /* start outer index loop: column index */ imnew = 0; i__3 = *nx; for (i = 1; i <= i__3; ++i) { /* start inner loop: row index */ if (grid.z[i + j * grid.nrow - grid.nr1] > -(doublereal)big) { goto L192; } /* UNFILLED GRID POINT, GET VALUE FROM NEAREST GR ID POINT */ /* WITH DATA */ if (i == 1 || imnew > 0) { goto L162; } zijn = (r__1 = grid.z[i - 1 + j * grid.nrow - grid.nr1], dabs(r__1)); if (zijn < big) { goto L195; } L162: if (j == 1 || jmnew[i - 1] > 0) { goto L172; } zijn = (r__1 = grid.z[i + (j - 1) * grid.nrow - grid.nr1], dabs(r__1)); if (zijn < big) { goto L195; } L172: if (i == *nx) { goto L182; } zijn = (r__1 = grid.z[i + 1 + j * grid.nrow - grid.nr1], dabs(r__1)); if (zijn < big) { goto L195; } L182: if (j == *ny) { goto L192; } zijn = (r__1 = grid.z[i + (j + 1) * grid.nrow - grid.nr1], dabs(r__1)); if (zijn < big) { goto L195; } L192: jmnew[i - 1] = 0; imnew = 0; goto L197; L195: jmnew[i - 1] = 1; imnew = 1; grid.z[i + j * grid.nrow - grid.nr1] = zijn; ++nnew; L197: ; } /* end of loop through rows: i */ /* L198: */ } /* end of column loop: j */ /* IF NO NEW GRID POINTS HAVE BEEN FILLED, FINISHED */ if (nnew <= 0) { goto L199; } } /* end of iteration loop: iter */ L199: free(jmnew); L200: /* BLANK ANY GRID POINTS STILL UNFILLED */ i__1 = *ny; for (j = 1; j <= i__1; ++j) { i__2 = *nx; for (i = 1; i <= i__2; ++i) { abz = (r__1 = grid.z[i + j * grid.nrow - grid.nr1], dabs(r__1)); if (abz >= big) { grid.z[i + j * grid.nrow - grid.nr1] = abz; } /* L202: */ } /* L203: */ } /* IMPROVE THE NON-DATA POINTS BY APPLYING POINT OVER-RELAXATION */ /* USING THE LAPLACE-SPLINE EQUATION (CARRE'S METHOD IS USED) */ /* EVERY NOPT ITERATIONS, ESTIMATE THE OPTIMUM RELAXATION */ /* PARAMETER. */ /* *********************************************************************** */ epslon = eps * zrange; dzrmsp = zrange * .1f; relax = 1.0; idivg = 0; i__1 = itmax; for (iter = 1; iter <= i__1; ++iter) { dzrms = 0.0; dzmax = 0.0; i__2 = *ny; for (j = 1; j <= i__2; ++j) { i__3 = *nx; for (i = 1; i <= i__3; ++i) { z00 = grid.z[i + j * grid.nrow - grid.nr1]; /* IF BLANKED OR FILLED, SKIP THIS GRID POINT */ if (z00 >= big || z00 < 0.0) { goto L2000; } wgt = 0.0; zsum = 0.0; im = 0; if (i == 1) { goto L570; } zim = (r__1 = grid.z[i - 1 + j * grid.nrow - grid.nr1], dabs(r__1)); if (zim >= big) { goto L570; } im = 1; wgt += 1.0; zsum += zim; if (i == 2) { goto L570; } zimm = (r__1 = grid.z[i - 2 + j * grid.nrow - grid.nr1], dabs(r__1)); if (zimm >= big) { goto L570; } wgt += cayscl; zsum -= cayscl * (zimm - zim * 2.0); L570: if (i == *nx) { goto L700; } zip = (r__1 = grid.z[i + 1 + j * grid.nrow - grid.nr1], dabs(r__1)); if (zip >= big) { goto L700; } wgt += 1.0; zsum += zip; if (im == 0) { goto L620; } wgt += cayscl * 4.f; zsum += cayscl * 2.0 * (zim + zip); L620: if (i == nxm) { goto L700; } zipp = (r__1 = grid.z[i + 2 + j * grid.nrow - grid.nr1], dabs(r__1)); if (zipp >= big) { goto L700; } wgt += cayscl; zsum -= cayscl * (zipp - zip * 2.0); L700: jm = 0; if (j == 1) { goto L1570; } zjm = (r__1 = grid.z[i + (j - 1) * grid.nrow - grid.nr1], dabs(r__1)); if (zjm >= big) { goto L1570; } jm = 1; wgt += del2; zsum += zjm * del2; if (j == 2) { goto L1570; } zjmm = (r__1 = grid.z[i + (j - 2) * grid.nrow - grid.nr1], dabs(r__1)); if (zjmm >= big) { goto L1570; } wgt += caysd4; zsum -= caysd4 * (zjmm - zjm * 2.0); L1570: if (j == *ny) { goto L1700; } zjp = (r__1 = grid.z[i + (j + 1) * grid.nrow - grid.nr1], dabs(r__1)); if (zjp >= big) { goto L1700; } wgt += del2; zsum += zjp * del2; if (jm == 0) { goto L1620; } wgt += caysd4 * 4.f; zsum += caysd4 * 2.0 * (zjm + zjp); L1620: if (j == nym) { goto L1700; } zjpp = (r__1 = grid.z[i + (j + 2) * grid.nrow - grid.nr1], dabs(r__1)); if (zjpp >= big) { goto L1700; } wgt += caysd4; zsum -= caysd4 * (zjpp - zjp * 2.0); L1700: dz = (zsum / wgt - z00) * relax; dzrms += dz * dz; /* Computing MAX */ r__1 = dabs(dz); dzmax = dmax(r__1,dzmax); grid.z[i + j * grid.nrow - grid.nr1] = z00 + dz; L2000: ; } /* L2001: */ } dzrms = sqrt(dzrms / npg); rtrms = dzrms / dzrmsp; dzrmsp = dzrms; if (rtrms >= 1.0) { ++idivg; } if (rtrms < 1.0) { idivg = 0; } if (idivg >= 10) { printf("Fails to converge: idivg >= 10\n"); return(-1); } /* TEST STOPPING CRITERION */ if (dzmax < epslon) { goto L2120; } if (iter != nopt) { goto L2100; } /* RECOMPUTE THE RELAXATION PARAMETER IN HOPES */ /* OF SPEEDING UP THE CONVERGENCE */ idivg = 0; den = sqrt(1.0 - rtrms * rtrms) + 1.0; relax = 2.0 / den; L2100: ; } /* SOLUTION HASN'T CONVERGED IN ITMAX ITERATIONS, SO */ /* WRITE OUT APPROPRIATE WARNING MESSAGE AND CONTINUE */ /* EXECUTION. */ L2120: /* ALL FINISHED WITH INTERPOLATION. CORRECT ALL GRID POINTS */ /* FOR THE DISPLACEMENT MADE IN GRIDDING THE DATA. */ /* *********************************************************************** */ corect = zbase + zrange * 10.0; i__1 = *ny; for (j = 1; j <= i__1; ++j) { i__2 = *nx; for (i = 1; i <= i__2; ++i) { if (grid.z[i + j * grid.nrow - grid.nr1] >= big) { goto L2500; } grid.z[i + j * grid.nrow - grid.nr1] = (r__1 = grid.z[i + j * grid.nrow - grid.nr1], dabs(r__1)) - corect; L2500: ; } /* L2501: */ } return 0; } /* zgrid_ */ /*============================================================*/ /* Below is a simple test routine for the main zgrid_ function. */ #if 0 void main() { int i,j, ret; int nx, ny, n, nrng, idbug; double x1, y1, dx, dy, del, cay; x1 = 0.0; y1 = 0.0; dx = 0.5; dy = 2; del = 0.1; cay = 1; n = 100; nx = 20; ny = 5; nrng = 1; idbug = 1; data_1.xp = malloc(100*sizeof(double)); data_1.yp = malloc(100*sizeof(double)); data_1.zp = malloc(100*sizeof(double)); grid.nrow = nx; grid.z = malloc(grid.nrow*ny*sizeof(double)); if (grid.z == NULL) { printf("Malloc failed\n"); return; } for (i = 0; i < nx; i++) for (j = 0; j < ny; j++) grid.z[i + grid.nrow*j] = 0.0; /* data_1 is also being set up fortran-fashion */ for (i = 0; i<10; i++) { for (j = 0; j < 10; j++) { data_1.xp[i+10*j] = i; data_1.yp[i+10*j] = j; data_1.zp[i+10*j] = sin(i) * (j+1); } } ret = zgrid_(&nx, &ny, &x1, &y1, &dx, &dy, &n, &del, &cay, &nrng, &idbug); printf("Return: %d\n", ret); for (i = 0; i < nx; i++) for (j = 0; j < ny; j++) printf("%d %d %f\n", i, j, grid.z[i + grid.nrow*j]); } #endif