Index: trunk/psLib/src/astro/psCoord.c =================================================================== --- trunk/psLib/src/astro/psCoord.c (revision 15253) +++ trunk/psLib/src/astro/psCoord.c (revision 15254) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.140 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-28 00:37:27 $ + * @version $Revision: 1.141 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-09 19:25:44 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -765,7 +765,7 @@ for (psS32 j = 0 ; j < orderY+1 ; j++) { myPT->x->coeff[i][j] = 0.0; - myPT->x->mask[i][j] = 0; myPT->y->coeff[i][j] = 0.0; - myPT->y->mask[i][j] = 0; + myPT->x->coeffMask[i][j] = PS_POLY_MASK_NONE; + myPT->y->coeffMask[i][j] = PS_POLY_MASK_NONE; } } @@ -807,24 +807,26 @@ for (psS32 t2y = 0 ; t2y < (trans2->x->nY + 1) ; t2y++) { psTrace("psLib.astro", 6, "X: -------------------- (t2x, t2y) (%d, %d) --------------------\n", t2x, t2y); - if (trans2->x->mask[t2x][t2y] == 0) { - psTrace("psLib.astro", 6, "In this iteration, we raise trans1->x to the %d power and trans1->y to the %d-power.\n", t2x, t2y); - psPolynomial2D *newPoly = multiplyDPoly2D(trans1XPolys[t2x], trans1YPolys[t2y]); - - if (psTraceGetLevel("psLib.astro") >= 6) { - PS_POLY_PRINT_2D(newPoly); - } - - // Set the appropriate coeffs in myPT->x - for (psS32 i = 0 ; i < (1 + newPoly->nX) ; i++) { - for (psS32 j = 0 ; j < (1 + newPoly->nY) ; j++) { - myPT->x->coeff[i][j]+= newPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y]; - } - } - - if (psTraceGetLevel("psLib.astro") >= 6) { - PS_POLY_PRINT_2D(myPT->x); - } - psFree(newPoly); - } + if (trans2->x->coeffMask[t2x][t2y] & PS_POLY_MASK_SET) { + continue; + } + + psTrace("psLib.astro", 6, "In this iteration, we raise trans1->x to the %d power and trans1->y to the %d-power.\n", t2x, t2y); + psPolynomial2D *newPoly = multiplyDPoly2D(trans1XPolys[t2x], trans1YPolys[t2y]); + + if (psTraceGetLevel("psLib.astro") >= 6) { + PS_POLY_PRINT_2D(newPoly); + } + + // Set the appropriate coeffs in myPT->x + for (psS32 i = 0 ; i < (1 + newPoly->nX) ; i++) { + for (psS32 j = 0 ; j < (1 + newPoly->nY) ; j++) { + myPT->x->coeff[i][j]+= newPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y]; + } + } + + if (psTraceGetLevel("psLib.astro") >= 6) { + PS_POLY_PRINT_2D(myPT->x); + } + psFree(newPoly); } } @@ -841,23 +843,24 @@ for (psS32 t2y = 0 ; t2y < (trans2->y->nY + 1) ; t2y++) { psTrace("psLib.astro", 5, "Y: -------------------- (t2x, t2y) (%d, %d) --------------------\n", t2x, t2y); - if (trans2->y->mask[t2x][t2y] == 0) { - psTrace("psLib.astro", 5, "In this iteration, we raise trans1->x to the %d power and trans1->y to the %d-power.\n", t2x, t2y); - psPolynomial2D *newPoly = multiplyDPoly2D(trans1XPolys[t2x], trans1YPolys[t2y]); - - if (psTraceGetLevel("psLib.astro") >= 6) { - PS_POLY_PRINT_2D(newPoly); - } - - // Set the appropriate coeffs in myPT->x - for (psS32 i = 0 ; i < (1 + newPoly->nX) ; i++) { - for (psS32 j = 0 ; j < (1 + newPoly->nY) ; j++) { - myPT->y->coeff[i][j]+= newPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y]; - } - } - if (psTraceGetLevel("psLib.astro") >= 6) { - PS_POLY_PRINT_2D(myPT->x); - } - psFree(newPoly); - } + if (trans2->y->coeffMask[t2x][t2y] & PS_POLY_MASK_SET) { + continue; + } + psTrace("psLib.astro", 5, "In this iteration, we raise trans1->x to the %d power and trans1->y to the %d-power.\n", t2x, t2y); + psPolynomial2D *newPoly = multiplyDPoly2D(trans1XPolys[t2x], trans1YPolys[t2y]); + + if (psTraceGetLevel("psLib.astro") >= 6) { + PS_POLY_PRINT_2D(newPoly); + } + + // Set the appropriate coeffs in myPT->x + for (psS32 i = 0 ; i < (1 + newPoly->nX) ; i++) { + for (psS32 j = 0 ; j < (1 + newPoly->nY) ; j++) { + myPT->y->coeff[i][j]+= newPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y]; + } + } + if (psTraceGetLevel("psLib.astro") >= 6) { + PS_POLY_PRINT_2D(myPT->x); + } + psFree(newPoly); } } Index: trunk/psLib/src/imageops/psImagePixelInterpolate.c =================================================================== --- trunk/psLib/src/imageops/psImagePixelInterpolate.c (revision 15253) +++ trunk/psLib/src/imageops/psImagePixelInterpolate.c (revision 15254) @@ -11,6 +11,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-20 23:54:25 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-09 19:25:44 $ * * Copyright 2007 Institute for Astronomy, University of Hawaii @@ -153,7 +153,7 @@ // allocate a 2D polynomial to fit a quadratic to the valid neighbor pixels. psPolynomial2D *poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); - poly->mask[2][2] = 1; - poly->mask[2][1] = 1; - poly->mask[1][2] = 1; + poly->coeffMask[2][2] = PS_POLY_MASK_SET; + poly->coeffMask[2][1] = PS_POLY_MASK_SET; + poly->coeffMask[1][2] = PS_POLY_MASK_SET; for (int iy = 0; iy < state->numRows; iy++) { @@ -275,11 +275,11 @@ // allocate a 2D polynomial to fit a quadratic to the valid neighbor pixels. psPolynomial2D *poly2o = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); - poly2o->mask[2][2] = 1; - poly2o->mask[2][1] = 1; - poly2o->mask[1][2] = 1; + poly2o->coeffMask[2][2] = PS_POLY_MASK_SET; + poly2o->coeffMask[2][1] = PS_POLY_MASK_SET; + poly2o->coeffMask[1][2] = PS_POLY_MASK_SET; // allocate a 2D polynomial to fit a plane to the valid neighbor pixels. psPolynomial2D *poly1o = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1); - poly2o->mask[1][1] = 1; + poly2o->coeffMask[1][1] = PS_POLY_MASK_SET; for (int iy = 0; iy < state->numRows; iy++) { Index: trunk/psLib/src/math/psMinimizePolyFit.c =================================================================== --- trunk/psLib/src/math/psMinimizePolyFit.c (revision 15253) +++ trunk/psLib/src/math/psMinimizePolyFit.c (revision 15254) @@ -10,6 +10,6 @@ * @author EAM, IfA * - * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-05-16 16:16:48 $ + * @version $Revision: 1.32 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-09 19:25:45 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -340,5 +340,5 @@ psImage *sums = psImageAlloc(numData, numTerms, PS_TYPE_F64); for (int i = 0; i < numTerms; i++) { - if (myPoly->mask[i]) { + if (myPoly->coeffMask[i] & PS_POLY_MASK_BOTH) { continue; } @@ -358,5 +358,5 @@ dataMask = mask->data.U8; } - psU8 *termMask = myPoly->mask; // Mask for polynomial terms + psU8 *coeffMask = myPoly->coeffMask; // Mask for polynomial terms psF64 *yData = y->data.F64; // Coordinate data psF64 *yErrData = NULL; // Errors in the coordinate @@ -380,5 +380,5 @@ for (int i = 0; i < numTerms; i++) { - if (termMask[i]) { + if (coeffMask[i] & PS_POLY_MASK_BOTH) { matrix[i][i] = 1.0; continue; @@ -387,5 +387,5 @@ matrix[i][i] += sumsData[i][k] * sumsData[i][k] * wt; // The diagonal entry for (int j = i + 1; j < numTerms; j++) { // The upper diagonal only: we will use symmetry - if (termMask[j]) { + if (coeffMask[j] & PS_POLY_MASK_BOTH) { continue; } @@ -533,5 +533,5 @@ dataMask = mask->data.U8; } - psU8 *termMask = myPoly->mask; // Dereferenced version of mask for polynomial terms + psU8 *coeffMask = myPoly->coeffMask; // Dereferenced version of mask for polynomial terms psF64 *ordinates = NULL; // Dereferenced version of ordinate data if (x) { @@ -567,5 +567,5 @@ for (int i = 0; i < nTerm; i++) { - if (termMask[i]) { + if (coeffMask[i] & PS_POLY_MASK_SET) { matrix[i][i] = 1.0; continue; @@ -574,5 +574,5 @@ matrix[i][i] += sums[2 * i] * wt; // The diagonal entry for (int j = i + 1; j < nTerm; j++) { // The upper diagonal only: we will use symmetry - if (termMask[j]) { + if (coeffMask[j] & PS_POLY_MASK_SET) { continue; } @@ -584,4 +584,32 @@ } psFree(xSums); + + // elements which are masked for fitting need to be subtracted from the vector + for (int i = 0; i < nTerm; i++) { + if (coeffMask[i] & PS_POLY_MASK_BOTH) { + continue; + } + for (int j = 0; j < nTerm; j++) { // The upper diagonal only: we will use symmetry + if (coeffMask[j] & PS_POLY_MASK_SET) { + continue; + } + if (!(coeffMask[j] & PS_POLY_MASK_FIT)) { + continue; + } + vector[i] -= matrix[i][j]*myPoly->coeff[j]; + } + } + + // set the un-fitted and un-set elements to 0 or 1 for pivots + for (int i = 0; i < nTerm; i++) { + if (coeffMask[i] & PS_POLY_MASK_BOTH) { + for (int j = 0; j < nTerm; j++) { // The upper diagonal only: we will use symmetry + matrix[i][j] = 0.0; + matrix[j][i] = 0.0; + } + matrix[i][i] = 1.0; + continue; + } + } if (psTraceGetLevel("psLib.math") >= 4) { @@ -610,9 +638,9 @@ // polynomial coefficients. this is only true for the 1D case for (psS32 k = 0; k < nTerm; k++) { - myPoly->coeff[k] = B->data.F64[k]; - myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]); - } - } - + if (coeffMask[k] & PS_POLY_MASK_FIT) continue; + myPoly->coeff[k] = B->data.F64[k]; + myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]); + } + } } else { // LUD version of the fit @@ -633,4 +661,5 @@ } else { for (psS32 k = 0; k < nTerm; k++) { + if (coeffMask[k] & PS_POLY_MASK_FIT) continue; myPoly->coeff[k] = coeffs->data.F64[k]; // XXX LUD does not give inverse of A @@ -990,5 +1019,5 @@ psFree(A); psFree(B); - psTrace("psLib.math", 4, "---- %s() End ----\n", __func__); + psTrace("psLib.math", 6, "---- %s() End ----\n", __func__); return false; } @@ -997,5 +1026,5 @@ psF64 **matrix = A->data.F64; // Dereference the least-squares matrix psF64 *vector = B->data.F64; // Dereference the least-squares vector - psU8 **termMask = myPoly->mask; // Dereference mask for polynomial terms + psU8 **coeffMask = myPoly->coeffMask; // Dereference mask for polynomial terms psU8 *dataMask = NULL; // Dereference mask for data if (mask) { @@ -1031,5 +1060,5 @@ int l = i / nYterm; // x index int m = i % nYterm; // y index - if (termMask[l][m]) { + if (coeffMask[l][m] & PS_POLY_MASK_SET) { matrix[i][i] = 1.0; continue; @@ -1040,5 +1069,5 @@ int p = j / nYterm; // x index int q = j % nYterm; // y index - if (termMask[p][q]) { + if (coeffMask[p][q] & PS_POLY_MASK_SET) { continue; } @@ -1050,4 +1079,53 @@ } psFree(xySums); + + // elements which are masked for fitting need to be subtracted from the vector + for (int i = 0; i < nTerm; i++) { + int ix = i / nYterm; // x index + int iy = i % nYterm; // y index + if (coeffMask[ix][iy] & PS_POLY_MASK_BOTH) { + continue; + } + for (int j = 0; j < nTerm; j++) { // The upper diagonal only: we will use symmetry + int jx = j / nYterm; // x index + int jy = j % nYterm; // y index + if (coeffMask[jx][jy] & PS_POLY_MASK_SET) { + continue; + } + if (!(coeffMask[jx][jy] & PS_POLY_MASK_FIT)) { + continue; + } + vector[i] -= matrix[i][j]*myPoly->coeff[jx][jy]; + } + } + + // set the un-fitted and un-set elements to 0 or 1 for pivots + for (int i = 0; i < nTerm; i++) { + int ix = i / nYterm; // x index + int iy = i % nYterm; // y index + if (coeffMask[ix][iy] & PS_POLY_MASK_BOTH) { + for (int j = 0; j < nTerm; j++) { // The upper diagonal only: we will use symmetry + matrix[i][j] = 0.0; + matrix[j][i] = 0.0; + } + matrix[i][i] = 1.0; + continue; + } + } + + if (psTraceGetLevel("psLib.math") >= 4) { + printf("Least-squares vector:\n"); + for (int i = 0; i < nTerm; i++) { + printf("%f ", B->data.F64[i]); + } + printf("\n"); + printf("Least-squares matrix:\n"); + for (int i = 0; i < nTerm; i++) { + for (int j = 0; j < nTerm; j++) { + printf("%f ", A->data.F64[i][j]); + } + printf("\n"); + } + } if (!psMatrixGJSolve(A, B)) { @@ -1062,11 +1140,14 @@ int l = i / nYterm; // x index int m = i % nYterm; // y index - myPoly->coeff[l][m] = B->data.F64[i]; - myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]); + + // retain the incoming values if masked on the fit + if (coeffMask[l][m] & PS_POLY_MASK_FIT) continue; + myPoly->coeff[l][m] = B->data.F64[i]; + myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]); } psFree(A); psFree(B); - psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); + psTrace("psLib.math", 6, "---- %s() end ----\n", __func__); return true; } @@ -1387,5 +1468,5 @@ dataMask = mask->data.U8; } - psU8 ***termMask = myPoly->mask; // Mask for polynomial terms + psU8 ***coeffMask = myPoly->coeffMask; // Mask for polynomial terms int nYZterm = nYterm * nZterm; // Multiplication of the numbers, to calculate the index @@ -1411,5 +1492,5 @@ int iy = (i % nYZterm) / nZterm; // y index int iz = (i % nYZterm) % nZterm; // z index - if (termMask[ix][iy][iz]) { + if (coeffMask[ix][iy][iz] & PS_POLY_MASK_BOTH) { matrix[i][i] = 1.0; continue; @@ -1422,5 +1503,5 @@ int jy = (j % nYZterm) / nZterm; // y index int jz = (j % nYZterm) % nZterm; // z index - if (termMask[jx][jy][jz]) { + if (coeffMask[jx][jy][jz] & PS_POLY_MASK_BOTH) { continue; } @@ -1455,4 +1536,5 @@ int iy = (i % nYZterm) / nZterm; // y index int iz = (i % nYZterm) % nZterm; // z index + if (coeffMask[ix][iy][iz] & PS_POLY_MASK_FIT) continue; myPoly->coeff[ix][iy][iz] = B->data.F64[i]; myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[i][i]); @@ -1480,4 +1562,5 @@ for (psS32 iz = 0; iz < nZterm; iz++) { psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm; + if (coeffMask[ix][iy][iz] & PS_POLY_MASK_FIT) continue; myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx]; // XXX myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]); @@ -1839,5 +1922,5 @@ dataMask = mask->data.U8; } - psU8 ****termMask = myPoly->mask; // Mask for polynomial terms + psU8 ****coeffMask = myPoly->coeffMask; // Mask for polynomial terms int nYZTterm = nYterm * nZterm * nTterm; // Multiplication of the numbers, for calculating the index int nZTterm = nZterm * nTterm; // Multiplication of the numbers, for calculating the index @@ -1865,5 +1948,5 @@ int iz = ((i % (nYZTterm)) % (nZTterm)) / nTterm; // z index int it = ((i % (nYZTterm)) % (nZTterm)) % nTterm; // t index - if (termMask[ix][iy][iz][it]) { + if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_BOTH) { matrix[i][i] = 1.0; continue; @@ -1877,5 +1960,5 @@ int jz = ((j % nYZTterm) % nZTterm) / nTterm; // z index int jt = ((j % nYZTterm) % nZTterm) % nTterm; // t index - if (termMask[jx][jy][jz][jt]) { + if (coeffMask[jx][jy][jz][jt] & PS_POLY_MASK_BOTH) { continue; } @@ -1918,4 +2001,5 @@ int iz = ((i % nYZTterm) % nZTterm) / nTterm; // z index int it = ((i % nYZTterm) % nZTterm) % nTterm; // t index + if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_FIT) continue; myPoly->coeff[ix][iy][iz][it] = B->data.F64[i]; myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[i][i]); @@ -1944,4 +2028,5 @@ for (psS32 it = 0; it < nTterm; it++) { psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm; + if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_FIT) continue; myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx]; // myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]); Index: trunk/psLib/src/math/psPolynomialMetadata.c =================================================================== --- trunk/psLib/src/math/psPolynomialMetadata.c (revision 15253) +++ trunk/psLib/src/math/psPolynomialMetadata.c (revision 15254) @@ -12,6 +12,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-02-06 21:55:28 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-10-09 19:25:45 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -62,9 +62,9 @@ // an undefined component implies the component was masked // this is symmetrical with the 1DtoMD function - poly->mask[nx] = 1; poly->coeff[nx] = 0; poly->coeffErr[nx] = 0; + poly->coeffMask[nx] = PS_POLY_MASK_SET; } else { - poly->mask[nx] = 0; + poly->coeffMask[nx] = PS_POLY_MASK_NONE; nElements ++; } @@ -123,5 +123,5 @@ // place polynomial entries on folder for (int nx = 0; nx < poly->nX + 1; nx++) { - if (poly->mask[nx] == 0) { + if (!(poly->coeffMask[nx] & PS_POLY_MASK_SET)) { sprintf(namespace, "VAL_X%02d", nx); sprintf(namespace_err, "ERR_X%02d", nx); @@ -174,9 +174,9 @@ // an undefined component implies the component was masked // this is symmetrical with the 2DtoMD function - poly->mask[nx][ny] = 1; poly->coeff[nx][ny] = 0; poly->coeffErr[nx][ny] = 0; + poly->coeffMask[nx][ny] = PS_POLY_MASK_SET; } else { - poly->mask[nx][ny] = 0; + poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE; nElements ++; } @@ -203,14 +203,8 @@ PS_ASSERT_PTR_NON_NULL(md, false); PS_ASSERT_PTR_NON_NULL(poly, false); - //XXX: Current implementation only supports ordinary polynomials. + + // XXX Current implementation only supports ordinary polynomials. if (poly->type != PS_POLYNOMIAL_ORD) return false; - - // XXX I'm puzzled by this test. a polynomial of 0 order with a value of 0 is a - // perfectly valid polynomial and can be written out. a polynomial with all elements - // masked still carries information. - //Make sure polynomial isn't 0, completely empty - //if (poly->nX == 0 && poly->nY == 0 && poly->coeff[0][0] == 0) - //return false; int Nbyte; @@ -245,5 +239,5 @@ for (int nx = 0; nx < poly->nX + 1; nx++) { for (int ny = 0; ny < poly->nY + 1; ny++) { - if (poly->mask[nx][ny] == 0) { + if (!(poly->coeffMask[nx][ny] & PS_POLY_MASK_SET)) { sprintf(namespace, "VAL_X%02d_Y%02d", nx, ny); sprintf(namespace_err, "ERR_X%02d_Y%02d", nx, ny); @@ -304,9 +298,9 @@ // an undefined component implies the component was masked // this is symmetrical with the 3DtoMD function - poly->mask[nx][ny][nz] = 1; poly->coeff[nx][ny][nz] = 0; poly->coeffErr[nx][ny][nz] = 0; + poly->coeffMask[nx][ny][nz] = PS_POLY_MASK_SET; } else { - poly->mask[nx][ny][nz] = 0; + poly->coeffMask[nx][ny][nz] = PS_POLY_MASK_NONE; nElements ++; } @@ -370,5 +364,5 @@ for (int ny = 0; ny < poly->nY + 1; ny++) { for (int nz = 0; nz < poly->nZ + 1; nz++) { - if (poly->mask[nx][ny][nz] == 0) { + if (!(poly->coeffMask[nx][ny][nz] & PS_POLY_MASK_SET)) { sprintf(namespace, "VAL_X%02d_Y%02d_Z%02d", nx, ny, nz); sprintf(namespace_err, "ERR_X%02d_Y%02d_Z%02d", nx, ny, nz); @@ -437,9 +431,9 @@ // an undefined component implies the component was masked // this is symmetrical with the 4DtoMD function - poly->mask[nx][ny][nz][nt] = 1; poly->coeff[nx][ny][nz][nt] = 0; poly->coeffErr[nx][ny][nz][nt] = 0; + poly->coeffMask[nx][ny][nz][nt] = PS_POLY_MASK_SET; } else { - poly->mask[nx][ny][nz][nt] = 0; + poly->coeffMask[nx][ny][nz][nt] = PS_POLY_MASK_NONE; nElements ++; } @@ -506,5 +500,5 @@ for (int nz = 0; nz < poly->nZ + 1; nz++) { for (int nt = 0; nt < poly->nT + 1; nt++) { - if (poly->mask[nx][ny][nz][nt] == 0) { + if (!(poly->coeffMask[nx][ny][nz][nt] & PS_POLY_MASK_SET)) { sprintf(namespace, "VAL_X%02d_Y%02d_Z%02d_T%02d", nx, ny, nz, nt); sprintf(namespace_err, "ERR_X%02d_Y%02d_Z%02d_T%02d", nx, ny, nz, nt); Index: trunk/psLib/src/math/psPolynomialUtils.c =================================================================== --- trunk/psLib/src/math/psPolynomialUtils.c (revision 15253) +++ trunk/psLib/src/math/psPolynomialUtils.c (revision 15254) @@ -159,7 +159,7 @@ psPolynomial2D *poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2); - poly->mask[2][2] = 1; - poly->mask[1][2] = 1; - poly->mask[2][1] = 1; + poly->coeffMask[2][2] = PS_POLY_MASK_SET; + poly->coeffMask[1][2] = PS_POLY_MASK_SET; + poly->coeffMask[2][1] = PS_POLY_MASK_SET; poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ; @@ -212,5 +212,5 @@ for (int i = 0; i < nXout + 1; i++) { for (int j = 0; j < nYout + 1; j++) { - out->mask[i][j] = poly->mask[i+1][j]; + out->coeffMask[i][j] = poly->coeffMask[i+1][j]; out->coeff[i][j] = poly->coeff[i+1][j] * (i+1); } @@ -240,5 +240,5 @@ for (int i = 0; i < nXout + 1; i++) { for (int j = 0; j < nYout + 1; j++) { - out->mask[i][j] = poly->mask[i][j+1]; + out->coeffMask[i][j] = poly->coeffMask[i][j+1]; out->coeff[i][j] = poly->coeff[i][j+1] * (j+1); }