IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 15254


Ignore:
Timestamp:
Oct 9, 2007, 9:27:04 AM (19 years ago)
Author:
eugene
Message:

changed mask to coeffMask, define PS_POLY_MASK_FIT,SET; distinguish between coeffs masked for fitting (FIT) and masked for evaluation (SET)

Location:
trunk
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/astro/psCoord.c

    r15050 r15254  
    1010 *  @author GLG, MHPCC
    1111 *
    12  *  @version $Revision: 1.140 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2007-09-28 00:37:27 $
     12 *  @version $Revision: 1.141 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2007-10-09 19:25:44 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    765765                for (psS32 j = 0 ; j < orderY+1 ; j++) {
    766766                    myPT->x->coeff[i][j] = 0.0;
    767                     myPT->x->mask[i][j] = 0;
    768767                    myPT->y->coeff[i][j] = 0.0;
    769                     myPT->y->mask[i][j] = 0;
     768                    myPT->x->coeffMask[i][j] = PS_POLY_MASK_NONE;
     769                    myPT->y->coeffMask[i][j] = PS_POLY_MASK_NONE;
    770770                }
    771771            }
    … …  
    807807        for (psS32 t2y = 0 ; t2y < (trans2->x->nY + 1) ; t2y++) {
    808808            psTrace("psLib.astro", 6, "X: -------------------- (t2x, t2y) (%d, %d) --------------------\n", t2x, t2y);
    809             if (trans2->x->mask[t2x][t2y] == 0) {
    810                 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);
    811                 psPolynomial2D *newPoly = multiplyDPoly2D(trans1XPolys[t2x], trans1YPolys[t2y]);
    812 
    813                 if (psTraceGetLevel("psLib.astro") >= 6) {
    814                     PS_POLY_PRINT_2D(newPoly);
    815                 }
    816 
    817                 // Set the appropriate coeffs in myPT->x
    818                 for (psS32 i = 0 ; i < (1 + newPoly->nX) ; i++) {
    819                     for (psS32 j = 0 ; j < (1 + newPoly->nY) ; j++) {
    820                         myPT->x->coeff[i][j]+= newPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y];
    821                     }
    822                 }
    823 
    824                 if (psTraceGetLevel("psLib.astro") >= 6) {
    825                     PS_POLY_PRINT_2D(myPT->x);
    826                 }
    827                 psFree(newPoly);
    828             }
     809            if (trans2->x->coeffMask[t2x][t2y] & PS_POLY_MASK_SET) {
     810                continue;
     811            }
     812
     813            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);
     814            psPolynomial2D *newPoly = multiplyDPoly2D(trans1XPolys[t2x], trans1YPolys[t2y]);
     815
     816            if (psTraceGetLevel("psLib.astro") >= 6) {
     817                PS_POLY_PRINT_2D(newPoly);
     818            }
     819
     820            // Set the appropriate coeffs in myPT->x
     821            for (psS32 i = 0 ; i < (1 + newPoly->nX) ; i++) {
     822                for (psS32 j = 0 ; j < (1 + newPoly->nY) ; j++) {
     823                    myPT->x->coeff[i][j]+= newPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y];
     824                }
     825            }
     826
     827            if (psTraceGetLevel("psLib.astro") >= 6) {
     828                PS_POLY_PRINT_2D(myPT->x);
     829            }
     830            psFree(newPoly);
    829831        }
    830832    }
    … …  
    841843        for (psS32 t2y = 0 ; t2y < (trans2->y->nY + 1) ; t2y++) {
    842844            psTrace("psLib.astro", 5, "Y: -------------------- (t2x, t2y) (%d, %d) --------------------\n", t2x, t2y);
    843             if (trans2->y->mask[t2x][t2y] == 0) {
    844                 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);
    845                 psPolynomial2D *newPoly = multiplyDPoly2D(trans1XPolys[t2x], trans1YPolys[t2y]);
    846 
    847                 if (psTraceGetLevel("psLib.astro") >= 6) {
    848                     PS_POLY_PRINT_2D(newPoly);
    849                 }
    850 
    851                 // Set the appropriate coeffs in myPT->x
    852                 for (psS32 i = 0 ; i < (1 + newPoly->nX) ; i++) {
    853                     for (psS32 j = 0 ; j < (1 + newPoly->nY) ; j++) {
    854                         myPT->y->coeff[i][j]+= newPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y];
    855                     }
    856                 }
    857                 if (psTraceGetLevel("psLib.astro") >= 6) {
    858                     PS_POLY_PRINT_2D(myPT->x);
    859                 }
    860                 psFree(newPoly);
    861             }
     845            if (trans2->y->coeffMask[t2x][t2y] & PS_POLY_MASK_SET) {
     846                continue;
     847            }
     848            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);
     849            psPolynomial2D *newPoly = multiplyDPoly2D(trans1XPolys[t2x], trans1YPolys[t2y]);
     850
     851            if (psTraceGetLevel("psLib.astro") >= 6) {
     852                PS_POLY_PRINT_2D(newPoly);
     853            }
     854
     855            // Set the appropriate coeffs in myPT->x
     856            for (psS32 i = 0 ; i < (1 + newPoly->nX) ; i++) {
     857                for (psS32 j = 0 ; j < (1 + newPoly->nY) ; j++) {
     858                    myPT->y->coeff[i][j]+= newPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y];
     859                }
     860            }
     861            if (psTraceGetLevel("psLib.astro") >= 6) {
     862                PS_POLY_PRINT_2D(myPT->x);
     863            }
     864            psFree(newPoly);
    862865        }
    863866    }

  • trunk/psLib/src/imageops/psImagePixelInterpolate.c

    r14924 r15254  
    1111 *  @author Eugene Magnier, IfA
    1212 *
    13  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2007-09-20 23:54:25 $
     13 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2007-10-09 19:25:44 $
    1515 *
    1616 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    153153    // allocate a 2D polynomial to fit a quadratic to the valid neighbor pixels.
    154154    psPolynomial2D *poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
    155     poly->mask[2][2] = 1;
    156     poly->mask[2][1] = 1;
    157     poly->mask[1][2] = 1;
     155    poly->coeffMask[2][2] = PS_POLY_MASK_SET;
     156    poly->coeffMask[2][1] = PS_POLY_MASK_SET;
     157    poly->coeffMask[1][2] = PS_POLY_MASK_SET;
    158158
    159159    for (int iy = 0; iy < state->numRows; iy++) {
    … …  
    275275    // allocate a 2D polynomial to fit a quadratic to the valid neighbor pixels.
    276276    psPolynomial2D *poly2o = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
    277     poly2o->mask[2][2] = 1;
    278     poly2o->mask[2][1] = 1;
    279     poly2o->mask[1][2] = 1;
     277    poly2o->coeffMask[2][2] = PS_POLY_MASK_SET;
     278    poly2o->coeffMask[2][1] = PS_POLY_MASK_SET;
     279    poly2o->coeffMask[1][2] = PS_POLY_MASK_SET;
    280280
    281281    // allocate a 2D polynomial to fit a plane to the valid neighbor pixels.
    282282    psPolynomial2D *poly1o = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
    283     poly2o->mask[1][1] = 1;
     283    poly2o->coeffMask[1][1] = PS_POLY_MASK_SET;
    284284
    285285    for (int iy = 0; iy < state->numRows; iy++) {

  • trunk/psLib/src/math/psMinimizePolyFit.c

    r13396 r15254  
    1010 *  @author EAM, IfA
    1111 *
    12  *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2007-05-16 16:16:48 $
     12 *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2007-10-09 19:25:45 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    340340    psImage *sums = psImageAlloc(numData, numTerms, PS_TYPE_F64);
    341341    for (int i = 0; i < numTerms; i++) {
    342         if (myPoly->mask[i]) {
     342        if (myPoly->coeffMask[i] & PS_POLY_MASK_BOTH) {
    343343            continue;
    344344        }
    … …  
    358358        dataMask = mask->data.U8;
    359359    }
    360     psU8 *termMask = myPoly->mask;      // Mask for polynomial terms
     360    psU8 *coeffMask = myPoly->coeffMask;      // Mask for polynomial terms
    361361    psF64 *yData = y->data.F64;         // Coordinate data
    362362    psF64 *yErrData = NULL;             // Errors in the coordinate
    … …  
    380380
    381381        for (int i = 0; i < numTerms; i++) {
    382             if (termMask[i]) {
     382            if (coeffMask[i] & PS_POLY_MASK_BOTH) {
    383383                matrix[i][i] = 1.0;
    384384                continue;
    … …  
    387387            matrix[i][i] += sumsData[i][k] * sumsData[i][k] * wt; // The diagonal entry
    388388            for (int j = i + 1; j < numTerms; j++) { // The upper diagonal only: we will use symmetry
    389                 if (termMask[j]) {
     389                if (coeffMask[j] & PS_POLY_MASK_BOTH) {
    390390                    continue;
    391391                }
    … …  
    533533        dataMask = mask->data.U8;
    534534    }
    535     psU8 *termMask = myPoly->mask;      // Dereferenced version of mask for polynomial terms
     535    psU8 *coeffMask = myPoly->coeffMask;      // Dereferenced version of mask for polynomial terms
    536536    psF64 *ordinates = NULL;            // Dereferenced version of ordinate data
    537537    if (x) {
    … …  
    567567
    568568        for (int i = 0; i < nTerm; i++) {
    569             if (termMask[i]) {
     569            if (coeffMask[i] & PS_POLY_MASK_SET) {
    570570                matrix[i][i] = 1.0;
    571571                continue;
    … …  
    574574            matrix[i][i] += sums[2 * i] * wt; // The diagonal entry
    575575            for (int j = i + 1; j < nTerm; j++) { // The upper diagonal only: we will use symmetry
    576                 if (termMask[j]) {
     576                if (coeffMask[j] & PS_POLY_MASK_SET) {
    577577                    continue;
    578578                }
    … …  
    584584    }
    585585    psFree(xSums);
     586
     587    // elements which are masked for fitting need to be subtracted from the vector
     588    for (int i = 0; i < nTerm; i++) {
     589        if (coeffMask[i] & PS_POLY_MASK_BOTH) {
     590            continue;
     591        }
     592        for (int j = 0; j < nTerm; j++) { // The upper diagonal only: we will use symmetry
     593            if (coeffMask[j] & PS_POLY_MASK_SET) {
     594                continue;
     595            }
     596            if (!(coeffMask[j] & PS_POLY_MASK_FIT)) {
     597                continue;
     598            }
     599            vector[i] -= matrix[i][j]*myPoly->coeff[j];
     600        }
     601    }
     602   
     603    // set the un-fitted and un-set elements to 0 or 1 for pivots
     604    for (int i = 0; i < nTerm; i++) {
     605        if (coeffMask[i] & PS_POLY_MASK_BOTH) {
     606            for (int j = 0; j < nTerm; j++) { // The upper diagonal only: we will use symmetry
     607                matrix[i][j] = 0.0;
     608                matrix[j][i] = 0.0;
     609            }
     610            matrix[i][i] = 1.0;
     611            continue;
     612        }
     613    }
    586614
    587615    if (psTraceGetLevel("psLib.math") >= 4) {
    … …  
    610638            // polynomial coefficients.  this is only true for the 1D case
    611639            for (psS32 k = 0; k < nTerm; k++) {
    612                 myPoly->coeff[k] = B->data.F64[k];
    613                 myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
    614             }
    615         }
    616 
     640                if (coeffMask[k] & PS_POLY_MASK_FIT) continue;
     641                myPoly->coeff[k] = B->data.F64[k];
     642                myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
     643            }
     644        }
    617645    } else {
    618646        // LUD version of the fit
    … …  
    633661            } else {
    634662                for (psS32 k = 0; k < nTerm; k++) {
     663                    if (coeffMask[k] & PS_POLY_MASK_FIT) continue;
    635664                    myPoly->coeff[k] = coeffs->data.F64[k];
    636665                    // XXX LUD does not give inverse of A
    … …  
    9901019        psFree(A);
    9911020        psFree(B);
    992         psTrace("psLib.math", 4, "---- %s() End ----\n", __func__);
     1021        psTrace("psLib.math", 6, "---- %s() End ----\n", __func__);
    9931022        return false;
    9941023    }
    … …  
    9971026    psF64 **matrix = A->data.F64;       // Dereference the least-squares matrix
    9981027    psF64 *vector = B->data.F64;        // Dereference the least-squares vector
    999     psU8 **termMask = myPoly->mask;     // Dereference mask for polynomial terms
     1028    psU8 **coeffMask = myPoly->coeffMask;     // Dereference mask for polynomial terms
    10001029    psU8 *dataMask = NULL;              // Dereference mask for data
    10011030    if (mask) {
    … …  
    10311060            int l = i / nYterm;         // x index
    10321061            int m = i % nYterm;         // y index
    1033             if (termMask[l][m]) {
     1062            if (coeffMask[l][m] & PS_POLY_MASK_SET) {
    10341063                matrix[i][i] = 1.0;
    10351064                continue;
    … …  
    10401069                int p = j / nYterm;     // x index
    10411070                int q = j % nYterm;     // y index
    1042                 if (termMask[p][q]) {
     1071                if (coeffMask[p][q] & PS_POLY_MASK_SET) {
    10431072                    continue;
    10441073                }
    … …  
    10501079    }
    10511080    psFree(xySums);
     1081
     1082    // elements which are masked for fitting need to be subtracted from the vector
     1083    for (int i = 0; i < nTerm; i++) {
     1084        int ix = i / nYterm;         // x index
     1085        int iy = i % nYterm;         // y index
     1086        if (coeffMask[ix][iy] & PS_POLY_MASK_BOTH) {
     1087            continue;
     1088        }
     1089        for (int j = 0; j < nTerm; j++) { // The upper diagonal only: we will use symmetry
     1090            int jx = j / nYterm;         // x index
     1091            int jy = j % nYterm;         // y index
     1092            if (coeffMask[jx][jy] & PS_POLY_MASK_SET) {
     1093                continue;
     1094            }
     1095            if (!(coeffMask[jx][jy] & PS_POLY_MASK_FIT)) {
     1096                continue;
     1097            }
     1098            vector[i] -= matrix[i][j]*myPoly->coeff[jx][jy];
     1099        }
     1100    }
     1101   
     1102    // set the un-fitted and un-set elements to 0 or 1 for pivots
     1103    for (int i = 0; i < nTerm; i++) {
     1104        int ix = i / nYterm;         // x index
     1105        int iy = i % nYterm;         // y index
     1106        if (coeffMask[ix][iy] & PS_POLY_MASK_BOTH) {
     1107            for (int j = 0; j < nTerm; j++) { // The upper diagonal only: we will use symmetry
     1108                matrix[i][j] = 0.0;
     1109                matrix[j][i] = 0.0;
     1110            }
     1111            matrix[i][i] = 1.0;
     1112            continue;
     1113        }
     1114    }
     1115
     1116    if (psTraceGetLevel("psLib.math") >= 4) {
     1117        printf("Least-squares vector:\n");
     1118        for (int i = 0; i < nTerm; i++) {
     1119            printf("%f ", B->data.F64[i]);
     1120        }
     1121        printf("\n");
     1122        printf("Least-squares matrix:\n");
     1123        for (int i = 0; i < nTerm; i++) {
     1124            for (int j = 0; j < nTerm; j++) {
     1125                printf("%f ", A->data.F64[i][j]);
     1126            }
     1127            printf("\n");
     1128        }
     1129    }
    10521130
    10531131    if (!psMatrixGJSolve(A, B)) {
    … …  
    10621140        int l = i / nYterm;         // x index
    10631141        int m = i % nYterm;         // y index
    1064         myPoly->coeff[l][m] = B->data.F64[i];
    1065         myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]);
     1142
     1143        // retain the incoming values if masked on the fit
     1144        if (coeffMask[l][m] & PS_POLY_MASK_FIT) continue;
     1145        myPoly->coeff[l][m] = B->data.F64[i];
     1146        myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]);
    10661147    }
    10671148    psFree(A);
    10681149    psFree(B);
    10691150
    1070     psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
     1151    psTrace("psLib.math", 6, "---- %s() end ----\n", __func__);
    10711152    return true;
    10721153}
    … …  
    13871468        dataMask = mask->data.U8;
    13881469    }
    1389     psU8 ***termMask = myPoly->mask;    // Mask for polynomial terms
     1470    psU8 ***coeffMask = myPoly->coeffMask;    // Mask for polynomial terms
    13901471    int nYZterm = nYterm * nZterm;      // Multiplication of the numbers, to calculate the index
    13911472
    … …  
    14111492            int iy = (i % nYZterm) / nZterm; // y index
    14121493            int iz = (i % nYZterm) % nZterm; // z index
    1413             if (termMask[ix][iy][iz]) {
     1494            if (coeffMask[ix][iy][iz] & PS_POLY_MASK_BOTH) {
    14141495                matrix[i][i] = 1.0;
    14151496                continue;
    … …  
    14221503                int jy = (j % nYZterm) / nZterm; // y index
    14231504                int jz = (j % nYZterm) % nZterm; // z index
    1424                 if (termMask[jx][jy][jz]) {
     1505                if (coeffMask[jx][jy][jz] & PS_POLY_MASK_BOTH) {
    14251506                    continue;
    14261507                }
    … …  
    14551536            int iy = (i % nYZterm) / nZterm; // y index
    14561537            int iz = (i % nYZterm) % nZterm; // z index
     1538            if (coeffMask[ix][iy][iz] & PS_POLY_MASK_FIT) continue;
    14571539            myPoly->coeff[ix][iy][iz] = B->data.F64[i];
    14581540            myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[i][i]);
    … …  
    14801562                        for (psS32 iz = 0; iz < nZterm; iz++) {
    14811563                            psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
     1564                            if (coeffMask[ix][iy][iz] & PS_POLY_MASK_FIT) continue;
    14821565                            myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx];
    14831566                            // XXX myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
    … …  
    18391922        dataMask = mask->data.U8;
    18401923    }
    1841     psU8 ****termMask = myPoly->mask;    // Mask for polynomial terms
     1924    psU8 ****coeffMask = myPoly->coeffMask;    // Mask for polynomial terms
    18421925    int nYZTterm = nYterm * nZterm * nTterm; // Multiplication of the numbers, for calculating the index
    18431926    int nZTterm = nZterm * nTterm;      // Multiplication of the numbers, for calculating the index
    … …  
    18651948            int iz = ((i % (nYZTterm)) % (nZTterm)) / nTterm; // z index
    18661949            int it = ((i % (nYZTterm)) % (nZTterm)) % nTterm; // t index
    1867             if (termMask[ix][iy][iz][it]) {
     1950            if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_BOTH) {
    18681951                matrix[i][i] = 1.0;
    18691952                continue;
    … …  
    18771960                int jz = ((j % nYZTterm) % nZTterm) / nTterm; // z index
    18781961                int jt = ((j % nYZTterm) % nZTterm) % nTterm; // t index
    1879                 if (termMask[jx][jy][jz][jt]) {
     1962                if (coeffMask[jx][jy][jz][jt] & PS_POLY_MASK_BOTH) {
    18801963                    continue;
    18811964                }
    … …  
    19182001            int iz = ((i % nYZTterm) % nZTterm) / nTterm; // z index
    19192002            int it = ((i % nYZTterm) % nZTterm) % nTterm; // t index
     2003            if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_FIT) continue;
    19202004            myPoly->coeff[ix][iy][iz][it] = B->data.F64[i];
    19212005            myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[i][i]);
    … …  
    19442028                            for (psS32 it = 0; it < nTterm; it++) {
    19452029                                psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
     2030                                if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_FIT) continue;
    19462031                                myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx];
    19472032                                // myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);

  • trunk/psLib/src/math/psPolynomialMetadata.c

    r11669 r15254  
    1212 *  @author Ross Harman, MHPCC
    1313 *
    14  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    15  *  @date $Date: 2007-02-06 21:55:28 $
     14 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     15 *  @date $Date: 2007-10-09 19:25:45 $
    1616 *
    1717 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    6262            // an undefined component implies the component was masked
    6363            // this is symmetrical with the 1DtoMD function
    64             poly->mask[nx] = 1;
    6564            poly->coeff[nx] = 0;
    6665            poly->coeffErr[nx] = 0;
     66            poly->coeffMask[nx] = PS_POLY_MASK_SET;
    6767        } else {
    68             poly->mask[nx] = 0;
     68            poly->coeffMask[nx] = PS_POLY_MASK_NONE;
    6969            nElements ++;
    7070        }
    … …  
    123123    // place polynomial entries on folder
    124124    for (int nx = 0; nx < poly->nX + 1; nx++) {
    125         if (poly->mask[nx] == 0) {
     125        if (!(poly->coeffMask[nx] & PS_POLY_MASK_SET)) {
    126126            sprintf(namespace, "VAL_X%02d", nx);
    127127            sprintf(namespace_err, "ERR_X%02d", nx);
    … …  
    174174                // an undefined component implies the component was masked
    175175                // this is symmetrical with the 2DtoMD function
    176                 poly->mask[nx][ny] = 1;
    177176                poly->coeff[nx][ny] = 0;
    178177                poly->coeffErr[nx][ny] = 0;
     178                poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
    179179            } else {
    180                 poly->mask[nx][ny] = 0;
     180                poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
    181181                nElements ++;
    182182            }
    … …  
    203203    PS_ASSERT_PTR_NON_NULL(md, false);
    204204    PS_ASSERT_PTR_NON_NULL(poly, false);
    205     //XXX:  Current implementation only supports ordinary polynomials.
     205
     206    // XXX Current implementation only supports ordinary polynomials.
    206207    if (poly->type != PS_POLYNOMIAL_ORD)
    207208        return false;
    208 
    209     // XXX I'm puzzled by this test.  a polynomial of 0 order with a value of 0 is a
    210     // perfectly valid polynomial and can be written out.  a polynomial with all elements
    211     // masked still carries information.
    212     //Make sure polynomial isn't 0, completely empty
    213     //if (poly->nX == 0 && poly->nY == 0 && poly->coeff[0][0] == 0)
    214     //return false;
    215209
    216210    int Nbyte;
    … …  
    245239    for (int nx = 0; nx < poly->nX + 1; nx++) {
    246240        for (int ny = 0; ny < poly->nY + 1; ny++) {
    247             if (poly->mask[nx][ny] == 0) {
     241            if (!(poly->coeffMask[nx][ny] & PS_POLY_MASK_SET)) {
    248242                sprintf(namespace, "VAL_X%02d_Y%02d", nx, ny);
    249243                sprintf(namespace_err, "ERR_X%02d_Y%02d", nx, ny);
    … …  
    304298                    // an undefined component implies the component was masked
    305299                    // this is symmetrical with the 3DtoMD function
    306                     poly->mask[nx][ny][nz] = 1;
    307300                    poly->coeff[nx][ny][nz] = 0;
    308301                    poly->coeffErr[nx][ny][nz] = 0;
     302                    poly->coeffMask[nx][ny][nz] = PS_POLY_MASK_SET;
    309303                } else {
    310                     poly->mask[nx][ny][nz] = 0;
     304                    poly->coeffMask[nx][ny][nz] = PS_POLY_MASK_NONE;
    311305                    nElements ++;
    312306                }
    … …  
    370364        for (int ny = 0; ny < poly->nY + 1; ny++) {
    371365            for (int nz = 0; nz < poly->nZ + 1; nz++) {
    372                 if (poly->mask[nx][ny][nz] == 0) {
     366                if (!(poly->coeffMask[nx][ny][nz] & PS_POLY_MASK_SET)) {
    373367                    sprintf(namespace, "VAL_X%02d_Y%02d_Z%02d", nx, ny, nz);
    374368                    sprintf(namespace_err, "ERR_X%02d_Y%02d_Z%02d", nx, ny, nz);
    … …  
    437431                        // an undefined component implies the component was masked
    438432                        // this is symmetrical with the 4DtoMD function
    439                         poly->mask[nx][ny][nz][nt] = 1;
    440433                        poly->coeff[nx][ny][nz][nt] = 0;
    441434                        poly->coeffErr[nx][ny][nz][nt] = 0;
     435                        poly->coeffMask[nx][ny][nz][nt] = PS_POLY_MASK_SET;
    442436                    } else {
    443                         poly->mask[nx][ny][nz][nt] = 0;
     437                        poly->coeffMask[nx][ny][nz][nt] = PS_POLY_MASK_NONE;
    444438                        nElements ++;
    445439                    }
    … …  
    506500            for (int nz = 0; nz < poly->nZ + 1; nz++) {
    507501                for (int nt = 0; nt < poly->nT + 1; nt++) {
    508                     if (poly->mask[nx][ny][nz][nt] == 0) {
     502                    if (!(poly->coeffMask[nx][ny][nz][nt] & PS_POLY_MASK_SET)) {
    509503                        sprintf(namespace, "VAL_X%02d_Y%02d_Z%02d_T%02d", nx, ny, nz, nt);
    510504                        sprintf(namespace_err, "ERR_X%02d_Y%02d_Z%02d_T%02d", nx, ny, nz, nt);

  • trunk/psLib/src/math/psPolynomialUtils.c

    r15049 r15254  
    159159
    160160    psPolynomial2D *poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
    161     poly->mask[2][2] = 1;
    162     poly->mask[1][2] = 1;
    163     poly->mask[2][1] = 1;
     161    poly->coeffMask[2][2] = PS_POLY_MASK_SET;
     162    poly->coeffMask[1][2] = PS_POLY_MASK_SET;
     163    poly->coeffMask[2][1] = PS_POLY_MASK_SET;
    164164
    165165    poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ;
    … …  
    212212    for (int i = 0; i < nXout + 1; i++) {
    213213        for (int j = 0; j < nYout + 1; j++) {
    214             out->mask[i][j] = poly->mask[i+1][j];
     214            out->coeffMask[i][j] = poly->coeffMask[i+1][j];
    215215            out->coeff[i][j] = poly->coeff[i+1][j] * (i+1);
    216216        }
    … …  
    240240    for (int i = 0; i < nXout + 1; i++) {
    241241        for (int j = 0; j < nYout + 1; j++) {
    242             out->mask[i][j] = poly->mask[i][j+1];
     242            out->coeffMask[i][j] = poly->coeffMask[i][j+1];
    243243            out->coeff[i][j] = poly->coeff[i][j+1] * (j+1);
    244244        }

  • trunk/psModules/src/astrom/pmAstrometryDistortion.c

    r13653 r15254  
    77*  @author EAM, IfA
    88*
    9 *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2007-06-06 00:43:19 $
     9*  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2007-10-09 19:27:04 $
    1111*
    1212*  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    6262
    6363    psPolynomial2D *local = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
    64     local->mask[1][1] = 1;
     64    local->coeffMask[1][1] = PS_POLY_MASK_SET;
    6565
    6666    // measure gradient for fractional chip regions
    … …  
    211211    for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
    212212        for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
    213             if (fpa->toTPA->x->mask[i][j]) {
    214                 if ((i > 0) && (i <= fpa->toTPA->x->nX)) {
    215                     localX->mask[i-1][j] = 1;
    216                 }
    217                 if ((j > 0) && (j <= fpa->toTPA->x->nY)) {
    218                     localY->mask[i][j-1] = 1;
    219                 }
    220             }
     213            if ((i > 0) && (i <= fpa->toTPA->x->nX)) {
     214                localX->coeffMask[i-1][j] = fpa->toTPA->x->coeffMask[i][j];
     215            }
     216            if ((j > 0) && (j <= fpa->toTPA->x->nY)) {
     217                localY->coeffMask[i][j-1] = fpa->toTPA->x->coeffMask[i][j];
     218            }
    221219        }
    222220    }
    … …  
    240238    fpa->toTPA->x->coeff[0][0] = 0;
    241239    for (int i = 1; i <= fpa->toTPA->x->nX; i++) {
    242         if (!fpa->toTPA->x->mask[i][0]) {
    243             fpa->toTPA->x->coeff[i][0] = localX->coeff[i-1][0] / i;
    244         }
     240        if (fpa->toTPA->x->coeffMask[i][0] & PS_POLY_MASK_SET) {
     241            continue;
     242        }
     243        fpa->toTPA->x->coeff[i][0] = localX->coeff[i-1][0] / i;
    245244    }
    246245    for (int j = 1; j <= fpa->toTPA->x->nY; j++) {
    247         if (!fpa->toTPA->x->mask[0][j]) {
    248             fpa->toTPA->x->coeff[0][j] = localY->coeff[0][j-1] / j;
    249         }
     246        if (fpa->toTPA->x->coeffMask[0][j] & PS_POLY_MASK_SET) {
     247            continue;
     248        }
     249        fpa->toTPA->x->coeff[0][j] = localY->coeff[0][j-1] / j;
    250250    }
    251251    for (int i = 1; i <= fpa->toTPA->x->nX; i++) {
    252252        for (int j = 1; j <= fpa->toTPA->x->nY; j++) {
    253             if (!fpa->toTPA->x->mask[i][j]) {
    254                 fpa->toTPA->x->coeff[i][j] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
    255             }
     253            if (fpa->toTPA->x->coeffMask[i][j] & PS_POLY_MASK_SET) {
     254                continue;
     255            }
     256            fpa->toTPA->x->coeff[i][j] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
    256257        }
    257258    }
    … …  
    268269    for (int i = 0; i < fpa->toTPA->y->nX; i++) {
    269270        for (int j = 0; j < fpa->toTPA->y->nY; j++) {
    270             if (fpa->toTPA->y->mask[i][j]) {
    271                 if ((i > 0) && (i <= fpa->toTPA->y->nX)) {
    272                     localX->mask[i-1][j] = 1;
    273                 }
    274                 if ((j > 0) && (j <= fpa->toTPA->y->nY)) {
    275                     localY->mask[i][j-1] = 1;
    276                 }
    277             }
     271            if ((i > 0) && (i <= fpa->toTPA->y->nX)) {
     272                localX->coeffMask[i-1][j] = fpa->toTPA->y->coeffMask[i][j];
     273            }
     274            if ((j > 0) && (j <= fpa->toTPA->y->nY)) {
     275                localY->coeffMask[i][j-1] = fpa->toTPA->y->coeffMask[i][j];
     276            }
    278277        }
    279278    }
    … …  
    286285    fpa->toTPA->y->coeff[0][0] = 0;
    287286    for (int i = 1; i <= fpa->toTPA->y->nX; i++) {
    288         if (!fpa->toTPA->y->mask[i][0]) {
    289             fpa->toTPA->y->coeff[i][0] = localX->coeff[i-1][0] / i;
    290         }
     287        if (fpa->toTPA->y->coeffMask[i][0] & PS_POLY_MASK_SET) {
     288            continue;
     289        }
     290        fpa->toTPA->y->coeff[i][0] = localX->coeff[i-1][0] / i;
    291291    }
    292292    for (int j = 1; j <= fpa->toTPA->y->nY; j++) {
    293         if (!fpa->toTPA->y->mask[0][j]) {
    294             fpa->toTPA->y->coeff[0][j] = localY->coeff[0][j-1] / j;
    295         }
     293        if (fpa->toTPA->y->coeffMask[0][j] & PS_POLY_MASK_SET) {
     294            continue;
     295        }
     296        fpa->toTPA->y->coeff[0][j] = localY->coeff[0][j-1] / j;
    296297    }
    297298    for (int i = 1; i <= fpa->toTPA->y->nX; i++) {
    298299        for (int j = 1; j <= fpa->toTPA->y->nY; j++) {
    299             if (!fpa->toTPA->y->mask[i][j]) {
    300                 fpa->toTPA->y->coeff[i][j] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
     300            if (fpa->toTPA->y->coeffMask[i][j] & PS_POLY_MASK_SET) {
     301                continue;
    301302            }
     303            fpa->toTPA->y->coeff[i][j] = 0.5*(localX->coeff[i-1][j] / i + localY->coeff[i][j-1] / j);
    302304        }
    303305    }

  • trunk/psModules/src/astrom/pmAstrometryUtils.c

    r12696 r15254  
    77 *  @author EAM, IfA
    88 *
    9  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-03-30 21:12:56 $
     9 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-10-09 19:27:04 $
    1111 *
    1212 *  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    105105        psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx);
    106106        for (int j = 0; j <= input->x->nY; j++) {
    107             output->x->mask[i][j] = input->x->mask[i][j];
    108             output->y->mask[i][j] = input->y->mask[i][j];
    109             output->x->coeff[i][j] = (output->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1);
    110             output->y->coeff[i][j] = (output->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1);
     107            output->x->coeffMask[i][j] = input->x->coeffMask[i][j];
     108            output->y->coeffMask[i][j] = input->y->coeffMask[i][j];
     109            output->x->coeff[i][j] = (output->x->coeffMask[i][j] & PS_POLY_MASK_SET) ? 0 : psPolynomial2DEval (xPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1);
     110            output->y->coeff[i][j] = (output->y->coeffMask[i][j] & PS_POLY_MASK_SET) ? 0 : psPolynomial2DEval (yPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1);
    111111
    112112            // take the next derivative wrt y, catch output (is NULL on last pass)
    … …  
    146146        for (int j = 0; j <= order; j++) {
    147147            if (i + j > order) {
    148                 transform->x->mask [i][j] = 1;
    149                 transform->y->mask [i][j] = 1;
     148                transform->x->coeffMask [i][j] = PS_POLY_MASK_SET;
     149                transform->y->coeffMask [i][j] = PS_POLY_MASK_SET;
    150150            }
    151151        }
    … …  
    153153    transform->x->coeff[1][0] = 1;
    154154    transform->y->coeff[0][1] = 1;
    155     transform->x->mask[1][0] = 0;
    156     transform->y->mask[0][1] = 0;
     155    transform->x->coeffMask[1][0] = PS_POLY_MASK_NONE;
     156    transform->y->coeffMask[0][1] = PS_POLY_MASK_NONE;
    157157
    158158    return transform;
    … …  
    197197            if (i + j > order) {
    198198                // high-order cross terms must be masked (eg, x^3 y^2)
    199                 status &= transform->x->mask[i][j];
     199                status &= (transform->x->coeffMask[i][j] & PS_POLY_MASK_SET);
    200200            } else {
    201                 status &= !transform->x->mask[i][j];
     201                status &= !(transform->x->coeffMask[i][j] & PS_POLY_MASK_SET);
    202202                if ((i == 1) && (i + j == 1)) {
    203203                    // linear, diagonal terms must be non-zero
    … …  
    216216            if (i + j > order) {
    217217                // high-order cross terms must be masked (eg, x^3 y^2)
    218                 status &= transform->y->mask[i][j];
     218                status &= (transform->y->coeffMask[i][j] & PS_POLY_MASK_SET);
    219219            } else {
    220                 status &= !transform->y->mask[i][j];
     220                status &= !(transform->y->coeffMask[i][j] & PS_POLY_MASK_SET);
    221221                if ((j == 1) && (i + j == 1)) {
    222222                    // linear, diagonal terms must be 1.0
    … …  
    249249        for (int j = 0; j <= order; j++) {
    250250            if (i + j > order) {
    251                 distort->x->mask [i][j][0][0] = 1;
    252                 distort->y->mask [i][j][0][0] = 1;
     251                distort->x->coeffMask [i][j][0][0] = PS_POLY_MASK_SET;
     252                distort->y->coeffMask [i][j][0][0] = PS_POLY_MASK_SET;
    253253            }
    254254        }
    … …  
    306306            if (i + j > order) {
    307307                // high-order cross terms must be masked (eg, x^3 y^2)
    308                 status &= distort->x->mask[i][j][0][0];
     308                status &= (distort->x->coeffMask[i][j][0][0] & PS_POLY_MASK_SET);
    309309            } else {
    310                 status &= !distort->x->mask[i][j][0][0];
     310                status &= !(distort->x->coeffMask[i][j][0][0] & PS_POLY_MASK_SET);
    311311                if ((i == 1) && (i + j == 1)) {
    312312                    // linear, diagonal terms must be 1.0
    … …  
    325325            if (i + j > order) {
    326326                // high-order cross terms must be masked (eg, x^3 y^2)
    327                 status &= distort->y->mask[i][j][0][0];
     327                status &= (distort->y->coeffMask[i][j][0][0] & PS_POLY_MASK_SET);
    328328            } else {
    329                 status &= !distort->y->mask[i][j][0][0];
     329                status &= !(distort->y->coeffMask[i][j][0][0] & PS_POLY_MASK_SET);
    330330                if ((j == 1) && (i + j == 1)) {
    331331                    // linear, diagonal terms must be 1.0

  • trunk/psModules/src/astrom/pmAstrometryWCS.c

    r12838 r15254  
    77 *  @author EAM, IfA
    88 *
    9  *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-04-17 19:39:04 $
     9 *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2007-10-09 19:27:04 $
    1111 *
    1212 *  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    298298                        continue;
    299299                    if (i + j > fitOrder) {
    300                         wcs->trans->x->mask[i][j] = 1;
    301                         wcs->trans->y->mask[i][j] = 1;
     300                        wcs->trans->x->coeffMask[i][j] = PS_POLY_MASK_SET;
     301                        wcs->trans->y->coeffMask[i][j] = PS_POLY_MASK_SET;
    302302                        continue;
    303303                    }

  • trunk/psModules/src/detrend/pmShutterCorrection.c

    r14935 r15254  
    230230
    231231    // mask out the terms we will not fit
    232     line->mask[0][0] = 1;
    233     line->mask[1][1] = 1;
     232    line->coeffMask[0][0] = PS_POLY_MASK_SET;
     233    line->coeffMask[1][1] = PS_POLY_MASK_SET;
    234234    line->coeff[0][0] = 0;
    235235    line->coeff[1][1] = 0;

  • trunk/psModules/src/objects/pmPSF_IO.c

    r15230 r15254  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2007-10-06 01:02:49 $
     8 *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2007-10-09 19:26:25 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    358358                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
    359359                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
    360                         psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);
     360                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->coeffMask[ix][iy]);
    361361
    362362                        psArrayAdd (psfTable, 100, row);
    … …  
    710710            assert (xPow <= poly->nX);
    711711            assert (yPow <= poly->nY);
    712             poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
    713             poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
    714             poly->mask[xPow][yPow]    = psMetadataLookupU8  (&status, row, "MASK");
     712            poly->coeff[xPow][yPow]     = psMetadataLookupF32 (&status, row, "VALUE");
     713            poly->coeffErr[xPow][yPow]  = psMetadataLookupF32 (&status, row, "ERROR");
     714            poly->coeffMask[xPow][yPow] = psMetadataLookupU8  (&status, row, "MASK");
    715715        }
    716716    }

  • trunk/psModules/src/objects/pmPSFtry.c

    r15016 r15254  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.48 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2007-09-25 22:05:05 $
     7 *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2007-10-09 19:26:25 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    280280    // linear clipped fit of ApResid to r2rflux
    281281    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
    282     poly->mask[1] = 1; // fit only a constant offset (no SKYBIAS)
     282    poly->coeffMask[1] = PS_POLY_MASK_SET; // fit only a constant offset (no SKYBIAS)
    283283
    284284    // XXX replace this with a psVectorStats call?  since we are not fitting the trend

  • trunk/psModules/src/objects/pmTrend2D.c

    r15000 r15254  
    33 *  @author EAM, IfA
    44 *
    5  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2007-09-24 21:27:58 $
     5 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2007-10-09 19:26:25 $
    77 *
    88 *  Copyright 2004 Institute for Astronomy, University of Hawaii
    … …  
    4848            for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
    4949                if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
    50                     trend->poly->mask[nx][ny] = 1;
     50                    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
    5151                } else {
    52                     trend->poly->mask[nx][ny] = 0;
     52                    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
    5353                }
    5454            }
    … …  
    9696            for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
    9797                if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
    98                     trend->poly->mask[nx][ny] = 1;
     98                    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
    9999                } else {
    100                     trend->poly->mask[nx][ny] = 0;
     100                    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
    101101                }
    102102            }
    … …  
    137137            for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
    138138                if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
    139                     trend->poly->mask[nx][ny] = 1;
     139                    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
    140140                } else {
    141                     trend->poly->mask[nx][ny] = 0;
     141                    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
    142142                }
    143143            }
Note: See TracChangeset for help on using the changeset viewer.