IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 17, 2006, 5:05:47 PM (20 years ago)
Author:
Paul Price
Message:

Updating polynomial fitting to work properly with masks.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r7133 r7135  
    1010 *  @author EAM, IfA
    1111 *
    12  *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2006-05-18 01:21:16 $
     12 *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2006-05-18 03:05:47 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    379379        for (int i = 0; i < numTerms; i++) {
    380380            if (termMask[i]) {
     381                matrix[i][i] = 1.0;
    381382                continue;
    382383            }
    … …  
    560561        for (int i = 0; i < nTerm; i++) {
    561562            if (termMask[i]) {
     563                matrix[i][i] = 1.0;
    562564                continue;
    563565            }
    … …  
    575577    }
    576578    psFree(xSums);
    577 
    578     #if 0
    579     // set masked elements in A,B appropriately
    580     // PAP: Is this necessary, given we initialise to zero, and we have already masked on the terms?
    581     for (int i = 0; i < nTerm; i++) {
    582         if (!termMask[i])
    583             continue;
    584         vector[i] = 0;
    585         matrix[i][i] = 1.0;
    586         for (int j = i + 1; j < nTerm; j++) {
    587             matrix[i][j] = 0.0;
    588             matrix[j][i] = 0.0;
    589         }
    590     }
    591     #endif
    592579
    593580    if (psTraceGetLevel(__func__) >= 4) {
    … …  
    10271014        // Iterating over the matrix
    10281015        for (int i = 0; i < nTerm; i++) {
    1029             int l = i % nXterm;         // x index
    1030             int m = i / nXterm;         // y index
     1016            int l = i / nYterm;         // x index
     1017            int m = i % nYterm;         // y index
    10311018            if (termMask[l][m]) {
     1019                matrix[i][i] = 1.0;
    10321020                continue;
    10331021            }
    … …  
    10351023            matrix[i][i] += sums[2*l][2*m] * wt; // The diagonal entry
    10361024            for (int j = i + 1; j < nTerm; j++) { // Doing the upper diagonal only: we will use symmetry
    1037                 int p = j % nXterm;     // x index
    1038                 int q = j / nXterm;     // y index
     1025                int p = j / nYterm;     // x index
     1026                int q = j % nYterm;     // y index
    10391027                if (termMask[p][q]) {
    10401028                    continue;
    … …  
    10541042    } else {
    10551043        // select the appropriate solution entries
    1056         for (int n = 0; n < nXterm; n++) {
    1057             for (int m = 0; m < nYterm; m++) {
    1058                 myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
    1059                 myPoly->coeffErr[n][m] = sqrt(A->data.F64[n+m*nXterm][n+m*nXterm]);
    1060             }
     1044        for (int i = 0; i < nTerm; i++) {
     1045            int l = i / nYterm;         // x index
     1046            int m = i % nYterm;         // y index
     1047            myPoly->coeff[l][m] = B->data.F64[i];
     1048            myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]);
    10611049        }
    10621050    }
    … …  
    14511439    }
    14521440    psU8 ***termMask = myPoly->mask;    // Mask for polynomial terms
    1453     int nXYterm = nXterm * nYterm;      // Multiplication of the numbers, to calculate the index
     1441    int nYZterm = nYterm * nZterm;      // Multiplication of the numbers, to calculate the index
    14541442
    14551443    // Build the B and A data structs.
    … …  
    14711459
    14721460        for (int i = 0; i < nTerm; i++) {
    1473             int iz = i / nXYterm; // z index
    1474             int iy = (i % nXYterm) / nXterm; // y index
    1475             int ix = (i % nXYterm) % nXterm; // x index
     1461            int ix = i / nYZterm; // x index
     1462            int iy = (i % nYZterm) / nZterm; // y index
     1463            int iz = (i % nYZterm) % nZterm; // z index
    14761464            if (termMask[ix][iy][iz]) {
     1465                matrix[i][i] = 1.0;
    14771466                continue;
    14781467            }
    … …  
    14811470            matrix[i][i] += Sums[2*ix][2*iy][2*iz] * wt;
    14821471            for (int j = i + 1; j < nTerm; j++) {
    1483                 int jz = j / (nXterm * nYterm); // z index
    1484                 int jy = (j % (nXterm * nYterm)) / nXterm; // y index
    1485                 int jx = (j % (nXterm * nYterm)) % nXterm; // x index
     1472                int jx = j / (nYZterm); // x index
     1473                int jy = (j % nYZterm) / nZterm; // y index
     1474                int jz = (j % nYZterm) % nZterm; // z index
    14861475                if (termMask[jx][jy][jz]) {
    14871476                    continue;
    … …  
    15151504        }
    15161505        // select the appropriate solution entries
    1517         for (psS32 ix = 0; ix < nXterm; ix++) {
    1518             for (psS32 iy = 0; iy < nYterm; iy++) {
    1519                 for (psS32 iz = 0; iz < nZterm; iz++) {
    1520                     psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
    1521                     myPoly->coeff[ix][iy][iz] = B->data.F64[nx];
    1522                     myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
    1523                 }
    1524             }
     1506        for (int i = 0; i < nTerm; i++) {
     1507            int ix = i / nYZterm; // x index
     1508            int iy = (i % nYZterm) / nZterm; // y index
     1509            int iz = (i % nYZterm) % nZterm; // z index
     1510            myPoly->coeff[ix][iy][iz] = B->data.F64[i];
     1511            myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[i][i]);
    15251512        }
    15261513    } else {
    … …  
    19331920    PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F64, NULL);
    19341921    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
    1935     if (mask != NULL) {
     1922    if (mask) {
    19361923        PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL);
    19371924        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
    … …  
    19761963    }
    19771964    psU8 ****termMask = myPoly->mask;    // Mask for polynomial terms
    1978     int nXYZterm = nXterm * nYterm * nZterm; // Multiplication of the numbers, for calculating the index
    1979     int nXYterm = nXterm * nYterm;      // Multiplication of the numbers, for calculating the index
     1965    int nYZTterm = nYterm * nZterm * nTterm; // Multiplication of the numbers, for calculating the index
     1966    int nZTterm = nZterm * nTterm;      // Multiplication of the numbers, for calculating the index
    19801967
    19811968    // Build the B and A data structs.
    … …  
    19971984
    19981985        for (int i = 0; i < nTerm; i++) {
    1999             int it = i / (nXYZterm); // t index
    2000             int iz = (i % (nXYZterm)) / (nXYterm); // z index
    2001             int iy = ((i % (nXYZterm)) % (nXYterm)) / nXterm; // y index
    2002             int ix = ((i % (nXYZterm)) % (nXYterm)) % nXterm; // x index
     1986            int ix = i / (nYZTterm); // x index
     1987            int iy = (i % (nYZTterm)) / (nZTterm); // y index
     1988            int iz = ((i % (nYZTterm)) % (nZTterm)) / nTterm; // z index
     1989            int it = ((i % (nYZTterm)) % (nZTterm)) % nTterm; // t index
    20031990            if (termMask[ix][iy][iz][it]) {
     1991                matrix[i][i] = 1.0;
    20041992                continue;
    20051993            }
    … …  
    20081996            matrix[i][i] += Sums[2*ix][2*iy][2*iz][2*it] * wt;
    20091997            for (int j = i + 1; j < nTerm; j++) {
    2010                 int jt = j / nXYZterm; // t index
    2011                 int jz = (j % nXYZterm) / nXYterm; // z index
    2012                 int jy = ((j % nXYZterm) % nXYterm) / nXterm; // y index
    2013                 int jx = ((j % nXYZterm) % nXYterm) % nXterm; // x index
     1998                int jx = j / nYZTterm; // x index
     1999                int jy = (j % nYZTterm) / nZTterm; // y index
     2000                int jz = ((j % nYZTterm) % nZTterm) / nTterm; // z index
     2001                int jt = ((j % nYZTterm) % nZTterm) % nTterm; // t index
    20142002                if (termMask[jx][jy][jz][jt]) {
    20152003                    continue;
    … …  
    20472035
    20482036        // select the appropriate solution entries
    2049         for (psS32 ix = 0; ix < nXterm; ix++) {
    2050             for (psS32 iy = 0; iy < nYterm; iy++) {
    2051                 for (psS32 iz = 0; iz < nZterm; iz++) {
    2052                     for (psS32 it = 0; it < nTterm; it++) {
    2053                         psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
    2054                         myPoly->coeff[ix][iy][iz][it] = B->data.F64[nx];
    2055                         myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
    2056                     }
    2057                 }
    2058             }
     2037        for (int i = 0; i < nTerm; i++) {
     2038            int ix = i / nYZTterm; // x index
     2039            int iy = (i % nYZTterm) / nZTterm; // y index
     2040            int iz = ((i % nYZTterm) % nZTterm) / nTterm; // z index
     2041            int it = ((i % nYZTterm) % nZTterm) % nTterm; // t index
     2042            myPoly->coeff[ix][iy][iz][it] = B->data.F64[i];
     2043            myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[i][i]);
    20592044        }
    20602045    } else {
    … …  
    21362121    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
    21372122    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
    2138     if (mask != NULL) {
     2123    if (mask) {
    21392124        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
    21402125        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
    … …  
    23072292    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
    23082293    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
    2309     if (mask != NULL) {
     2294    if (mask) {
    23102295        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
    23112296        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
Note: See TracChangeset for help on using the changeset viewer.