IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 7107


Ignore:
Timestamp:
May 10, 2006, 3:24:57 AM (20 years ago)
Author:
Paul Price
Message:

Updating polynomial fitting for 3 and 4 dimensions

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r7104 r7107  
    1010 *  @author EAM, IfA
    1111 *
    12  *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2006-05-10 11:38:55 $
     12 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2006-05-10 13:24:57 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    14151415    }
    14161416
    1417     psImage    *A = NULL;
    1418     psVector   *B = NULL;
    1419     psF64 ***Sums = NULL;
    1420     psF64 wt;
    1421     psS32 nTerm;
    1422 
    1423     psS32 nXterm = 1 + myPoly->nX;
    1424     psS32 nYterm = 1 + myPoly->nY;
    1425     psS32 nZterm = 1 + myPoly->nZ;
    1426     nTerm = nXterm * nYterm * nZterm;
    1427 
    1428     A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
    1429     B = psVectorAlloc(nTerm, PS_TYPE_F64);
     1417
     1418    int nXterm = 1 + myPoly->nX;        // Number of x terms
     1419    int nYterm = 1 + myPoly->nY;        // Number of y terms
     1420    int nZterm = 1 + myPoly->nZ;        // Number of z terms
     1421    int nTerm = nXterm * nYterm * nZterm; // Total number of terms
     1422    int nData = x->n;                   // Number of data points
     1423    psImage    *A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); // Least-squares matrix
     1424    psVector   *B = psVectorAlloc(nTerm, PS_TYPE_F64); // Least-squares vector
    14301425    B->n = B->nalloc;
     1426
    14311427    // Initialize data structures.
    14321428    if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) {
    … …  
    14391435    }
    14401436
    1441     // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ...
     1437    // Dereference points for speed in the loop
     1438    psF64 **matrix = A->data.F64;       // Least-squares matrix
     1439    psF64 *vector = B->data.F64;        // Least-squares vector
     1440    psF64 *xData = x->data.F64;         // x
     1441    psF64 *yData = y->data.F64;         // y
     1442    psF64 *zData = z->data.F64;         // z
     1443    psF64 *fData = f->data.F64;         // f
     1444    psF64 *fErrData = NULL;             // Error in f
     1445    if (fErr) {
     1446        fErrData = fErr->data.F64;
     1447    }
     1448    psU8 *dataMask = NULL;              // Mask for data
     1449    if (mask) {
     1450        dataMask = mask->data.U8;
     1451    }
     1452    psU8 ***termMask = myPoly->mask;    // Mask for polynomial terms
     1453    int nXYterm = nXterm * nYterm;      // Multiplication of the numbers, to calculate the index
    14421454
    14431455    // Build the B and A data structs.
    1444     for (psS32 k  = 0; k < x->n; k++) {
    1445         if ((mask != NULL) && (mask->data.U8[k] & maskValue)) {
     1456    psF64 ***Sums = NULL;         // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ...
     1457    for (int k = 0; k < nData; k++) {
     1458        if (dataMask && dataMask[k] & maskValue) {
    14461459            continue;
    14471460        }
    14481461
    1449         Sums = BuildSums3D(Sums, x->data.F64[k], y->data.F64[k], z->data.F64[k], nXterm, nYterm, nZterm);
    1450 
     1462        Sums = BuildSums3D(Sums, xData[k], yData[k], zData[k], nXterm, nYterm, nZterm);
     1463
     1464        double wt;
    14511465        if (fErr == NULL) {
    14521466            wt = 1.0;
    14531467        } else {
    14541468            // this filters fErr == 0 values
    1455             wt = (fErr->data.F64[k] == 0.0) ? 0.0 : 1.0 / PS_SQR(fErr->data.F64[k]);
    1456         }
    1457 
    1458         // we could skip half of the array and assign at the end
    1459         // we must handle masked orders
    1460         for (psS32 ix = 0; ix < nXterm; ix++) {
    1461             for (psS32 iy = 0; iy < nYterm; iy++) {
    1462                 for (psS32 iz = 0; iz < nZterm; iz++) {
    1463                     if (myPoly->mask[ix][iy][iz]) {
    1464                         continue;
    1465                     } else {
    1466                         psS32 nx = ix + iy*nXterm + iz*nXterm*nYterm;
    1467                         B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz] * wt;
    1468                     }
    1469                 }
    1470             }
    1471         }
    1472 
    1473         for (psS32 ix = 0; ix < nXterm; ix++) {
    1474             for (psS32 iy = 0; iy < nYterm; iy++) {
    1475                 for (psS32 iz = 0; iz < nZterm; iz++) {
    1476                     if (myPoly->mask[ix][iy][iz])
    1477                         continue;
    1478                     psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
    1479                     for (psS32 jx = 0; jx < nXterm; jx++) {
    1480                         for (psS32 jy = 0; jy < nYterm; jy++) {
    1481                             for (psS32 jz = 0; jz < nZterm; jz++) {
    1482                                 if (myPoly->mask[jx][jy][jz])
    1483                                     continue;
    1484                                 psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm;
    1485                                 A->data.F64[nx][ny] += Sums[ix+jx][iy+jy][iz+jz] * wt;
    1486                             }
    1487                         }
    1488                     }
    1489                 }
    1490             }
    1491         }
    1492     }
    1493 
    1494     for (psS32 ix = 0; ix < nXterm; ix++) {
    1495         for (psS32 iy = 0; iy < nYterm; iy++) {
    1496             for (psS32 iz = 0; iz < nZterm; iz++) {
    1497                 if (!myPoly->mask[ix][iy][iz])
     1469            wt = (fErr->data.F64[k] == 0.0) ? 0.0 : 1.0 / PS_SQR(fErrData[k]);
     1470        }
     1471
     1472        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
     1476            if (termMask[ix][iy][iz]) {
     1477                continue;
     1478            }
     1479
     1480            vector[i] += fData[k] * Sums[ix][iy][iz] * wt;
     1481            matrix[i][i] += Sums[2*ix][2*iy][2*iz] * wt;
     1482            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
     1486                if (termMask[jx][jy][jz]) {
    14981487                    continue;
    1499                 psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
    1500                 B->data.F64[nx] = 0;
    1501                 for (psS32 jx = 0; jx < nXterm; jx++) {
    1502                     for (psS32 jy = 0; jy < nYterm; jy++) {
    1503                         for (psS32 jz = 0; jz < nZterm; jz++) {
    1504                             psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm;
    1505                             A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
    1506                         }
    1507                     }
    1508                 }
    1509             }
    1510         }
    1511     }
     1488                }
     1489                double value = Sums[ix+jx][iy+jy][iz+jz] * wt;
     1490                matrix[i][j] += value;
     1491                matrix[j][i] += value;
     1492            }
     1493        }
     1494    }
     1495
     1496    // Free the sums
     1497    for (psS32 ix = 0; ix < 2*nXterm; ix++) {
     1498        for (psS32 iy = 0; iy < 2*nYterm; iy++) {
     1499            psFree(Sums[ix][iy]);
     1500        }
     1501        psFree(Sums[ix]);
     1502    }
     1503    psFree(Sums);
     1504
    15121505
    15131506    // XXX: rel10_ifa used psMatrixGJSolve().  However, this failed tests.  So, I'm using psMatrixLUD().
    … …  
    15181511            psFree(A);
    15191512            psFree(B);
    1520 
    1521             for (psS32 ix = 0; ix < 2*nXterm; ix++) {
    1522                 for (psS32 iy = 0; iy < 2*nYterm; iy++) {
    1523                     psFree(Sums[ix][iy]);
    1524                 }
    1525                 psFree(Sums[ix]);
    1526             }
    1527             psFree(Sums);
    15281513            psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n");
    15291514            return(NULL);
    … …  
    15781563    psFree(A);
    15791564    psFree(B);
    1580 
    1581     for (psS32 ix = 0; ix < 2*nXterm; ix++) {
    1582         for (psS32 iy = 0; iy < 2*nYterm; iy++) {
    1583             psFree(Sums[ix][iy]);
    1584         }
    1585         psFree(Sums[ix]);
    1586     }
    1587     psFree(Sums);
    15881565
    15891566    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    … …  
    19611938    }
    19621939
    1963     // I think this is 1 dimension down
    1964     psImage     *A = NULL;
    1965     psVector    *B = NULL;
    1966     psF64 ****Sums = NULL;
    1967     psF64 wt;
    1968     psS32 nTerm;
    1969 
    1970     psS32 nXterm = 1 + myPoly->nX;
    1971     psS32 nYterm = 1 + myPoly->nY;
    1972     psS32 nZterm = 1 + myPoly->nZ;
    1973     psS32 nTterm = 1 + myPoly->nZ;
    1974     nTerm = nXterm * nYterm * nZterm * nTterm;
    1975 
    1976     A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
    1977     B = psVectorAlloc(nTerm, PS_TYPE_F64);
     1940
     1941    int nXterm = 1 + myPoly->nX;        // Number of x terms
     1942    int nYterm = 1 + myPoly->nY;        // Number of y terms
     1943    int nZterm = 1 + myPoly->nZ;        // Number of z terms
     1944    int nTterm = 1 + myPoly->nT;        // Number of t terms
     1945    int nTerm = nXterm * nYterm * nZterm * nTterm; // Total number of terms
     1946    int nData = x->n;                   // Number of data points
     1947    psImage    *A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); // Least-squares matrix
     1948    psVector   *B = psVectorAlloc(nTerm, PS_TYPE_F64); // Least-squares vector
    19781949    B->n = B->nalloc;
     1950
    19791951    // Initialize data structures.
    19801952    if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) {
    … …  
    19871959    }
    19881960
    1989     // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ...
     1961    // Dereference points for speed in the loop
     1962    psF64 **matrix = A->data.F64;       // Least-squares matrix
     1963    psF64 *vector = B->data.F64;        // Least-squares vector
     1964    psF64 *xData = x->data.F64;         // x
     1965    psF64 *yData = y->data.F64;         // y
     1966    psF64 *zData = z->data.F64;         // z
     1967    psF64 *tData = t->data.F64;         // t
     1968    psF64 *fData = f->data.F64;         // f
     1969    psF64 *fErrData = NULL;             // Error in f
     1970    if (fErr) {
     1971        fErrData = fErr->data.F64;
     1972    }
     1973    psU8 *dataMask = NULL;              // Mask for data
     1974    if (mask) {
     1975        dataMask = mask->data.U8;
     1976    }
     1977    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
    19901980
    19911981    // Build the B and A data structs.
    1992     for (psS32 k  = 0; k < x->n; k++) {
    1993         if ((mask != NULL) && (mask->data.U8[k] & maskValue)) {
     1982    psF64 ****Sums = NULL;        // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ...
     1983    for (int k = 0; k < nData; k++) {
     1984        if (dataMask && dataMask[k] & maskValue) {
    19941985            continue;
    19951986        }
    19961987
    1997         Sums = BuildSums4D(Sums, x->data.F64[k], y->data.F64[k], z->data.F64[k], t->data.F64[k], nXterm, nYterm, nZterm, nTterm);
    1998 
     1988        Sums = BuildSums4D(Sums, xData[k], yData[k], zData[k], tData[k], nXterm, nYterm, nZterm, nTterm);
     1989
     1990        double wt;
    19991991        if (fErr == NULL) {
    20001992            wt = 1.0;
    20011993        } else {
    20021994            // this filters fErr == 0 values
    2003             wt = (fErr->data.F64[k] == 0.0) ? 0.0 : 1.0 / PS_SQR(fErr->data.F64[k]);
    2004         }
    2005 
    2006         // we could skip half of the array and assign at the end
    2007         // we must handle masked orders
    2008         for (psS32 ix = 0; ix < nXterm; ix++) {
    2009             for (psS32 iy = 0; iy < nYterm; iy++) {
    2010                 for (psS32 iz = 0; iz < nZterm; iz++) {
    2011                     for (psS32 it = 0; it < nTterm; it++) {
    2012                         if (myPoly->mask[ix][iy][iz][it])
    2013                             continue;
    2014                         psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
    2015                         B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz][it] * wt;
    2016                     }
    2017                 }
    2018             }
    2019         }
    2020 
    2021         for (psS32 ix = 0; ix < nXterm; ix++) {
    2022             for (psS32 iy = 0; iy < nYterm; iy++) {
    2023                 for (psS32 iz = 0; iz < nZterm; iz++) {
    2024                     for (psS32 it = 0; it < nTterm; it++) {
    2025                         if (myPoly->mask[ix][iy][iz][it])
    2026                             continue;
    2027                         psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
    2028                         for (psS32 jx = 0; jx < nXterm; jx++) {
    2029                             for (psS32 jy = 0; jy < nYterm; jy++) {
    2030                                 for (psS32 jz = 0; jz < nZterm; jz++) {
    2031                                     for (psS32 jt = 0; jt < nTterm; jt++) {
    2032                                         if (myPoly->mask[jx][jy][jz][jt])
    2033                                             continue;
    2034                                         psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
    2035                                         A->data.F64[nx][ny]+= Sums[ix+jx][iy+jy][iz+jz][it+jt] * wt;
    2036                                     }
    2037                                 }
    2038                             }
    2039                         }
    2040                     }
    2041                 }
    2042             }
    2043         }
    2044     }
    2045 
    2046     for (psS32 ix = 0; ix < nXterm; ix++) {
    2047         for (psS32 iy = 0; iy < nYterm; iy++) {
    2048             for (psS32 iz = 0; iz < nZterm; iz++) {
    2049                 for (psS32 it = 0; it < nTterm; it++) {
    2050                     if (!myPoly->mask[ix][iy][iz][it])
    2051                         continue;
    2052                     psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
    2053                     B->data.F64[nx] = 0;
    2054                     for (psS32 jx = 0; jx < nXterm; jx++) {
    2055                         for (psS32 jy = 0; jy < nYterm; jy++) {
    2056                             for (psS32 jz = 0; jz < nZterm; jz++) {
    2057                                 for (psS32 jt = 0; jt < nTterm; jt++) {
    2058                                     psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
    2059                                     A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
    2060                                 }
    2061                             }
    2062                         }
    2063                     }
    2064                 }
    2065             }
    2066         }
    2067     }
     1995            wt = (fErr->data.F64[k] == 0.0) ? 0.0 : 1.0 / PS_SQR(fErrData[k]);
     1996        }
     1997
     1998        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
     2003            if (termMask[ix][iy][iz][it]) {
     2004                continue;
     2005            }
     2006
     2007            vector[i] += fData[k] * Sums[ix][iy][iz][it] * wt;
     2008            matrix[i][i] += Sums[2*ix][2*iy][2*iz][2*it] * wt;
     2009            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
     2014                if (termMask[jx][jy][jz][jt]) {
     2015                    continue;
     2016                }
     2017                double value = Sums[ix+jx][iy+jy][iz+jz][it+jt] * wt;
     2018                matrix[i][j] += value;
     2019                matrix[j][i] += value;
     2020            }
     2021        }
     2022    }
     2023
     2024    // Free the sums
     2025    for (int ix = 0; ix < 2*nXterm; ix++) {
     2026        for (int iy = 0; iy < 2*nYterm; iy++) {
     2027            for (int iz = 0; iz < 2*nZterm; iz++) {
     2028                psFree(Sums[ix][iy][iz]);
     2029            }
     2030            psFree(Sums[ix][iy]);
     2031        }
     2032        psFree(Sums[ix]);
     2033    }
     2034    psFree(Sums);
    20682035
    20692036
    … …  
    20722039        // does the solution in place
    20732040        // The GaussJordan version was overflowing, so I'm using LUD.
    2074         if (psMatrixGJSolve(A, B)) {
     2041        if (!psMatrixGJSolve(A, B)) {
    20752042            psFree(A);
    20762043            psFree(B);
    2077             for (psS32 ix = 0; ix < 2*nXterm; ix++) {
    2078                 for (psS32 iy = 0; iy < 2*nYterm; iy++) {
    2079                     for (psS32 iz = 0; iz < 2*nZterm; iz++) {
    2080                         psFree(Sums[ix][iy][iz]);
    2081                     }
    2082                     psFree(Sums[ix][iy]);
    2083                 }
    2084                 psFree(Sums[ix]);
    2085             }
    2086             psFree(Sums);
    20872044            psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n");
    20882045            return(NULL);
    … …  
    21452102    psFree(B);
    21462103
    2147     for (psS32 ix = 0; ix < 2*nXterm; ix++) {
    2148         for (psS32 iy = 0; iy < 2*nYterm; iy++) {
    2149             for (psS32 iz = 0; iz < 2*nZterm; iz++) {
    2150                 psFree(Sums[ix][iy][iz]);
    2151             }
    2152             psFree(Sums[ix][iy]);
    2153         }
    2154         psFree(Sums[ix]);
    2155     }
    2156     psFree(Sums);
    2157 
    21582104    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    21592105    return (myPoly);
Note: See TracChangeset for help on using the changeset viewer.