IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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]);
Note: See TracChangeset for help on using the changeset viewer.