IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 4101


Ignore:
Timestamp:
Jun 3, 2005, 12:30:20 PM (21 years ago)
Author:
gusciora
Message:

Added bug 366 suggestions.

Location:
trunk/psLib/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psStats.c

    r4029 r4101  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.129 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-05-25 20:26:55 $
     11 *  @version $Revision: 1.130 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-06-03 22:30:20 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    16281628    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
    16291629            (stats->options & PS_STAT_ROBUST_STDEV)) {
    1630 
    1631         // Using the above (myMean, myStdev) as initial estimates, we fit a
    1632         // Gaussian to the robustHistogramVector.
    1633         psImage *covar = NULL;
    1634         /* XXX: Old, remove.
    1635                 psMinimization *min = psMinimizationAlloc(100, 0.1);
    1636                 psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
    1637                 psArray *myCoords = psArrayAlloc(robustHistogramVector->n);
    1638                 psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
    1639                 myParams->data.F32[0] = myMean;
    1640                 myParams->data.F32[1] = myStdev;
    1641                 for (i=0;i<robustHistogramVector->n;i++) {
    1642                     myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
    1643                     ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i;
    1644                     y->data.F32[i] = robustHistogramVector->data.F32[i];
    1645                 }
    1646         */
    1647 
    1648         psMinimization *min = psMinimizationAlloc(100, 0.01);
    1649         psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
    1650         psArray *myCoords = psArrayAlloc(2 * dL + 1);
     1630        // We fit a 1-D polynomial to the data.
     1631        // XXX: Since we are no longer fitting a 1-D Gaussian, we can probably
     1632        // remove some of theabove code that calculated the initial estimate
     1633        // for the mean and sigma.
     1634
    16511635        psVector *y = psVectorAlloc(2 * dL + 1, PS_TYPE_F32);
    1652         myParams->data.F32[0] = myMean;
    1653         myParams->data.F32[1] = myStdev;
    1654 
    16551636        for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
    16561637            int index = i - modeBinNum + dL;
    1657             myCoords->data[index] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
    1658             ((psVector *) (myCoords->data[index]))->data.F32[0] = PS_BIN_MIDPOINT(robustHistogram, i);
     1638            // XXX: Should this be the natural log?
    16591639            y->data.F32[index] = robustHistogramVector->data.F32[i];
    1660         }
    1661 
    1662         // XXX: Must fit a 2-D polynomial, not a Gaussian.  (bug 366)
    1663         psBool rc = psMinimizeLMChi2(min,
    1664                                      NULL,
    1665                                      myParams,
    1666                                      NULL,
    1667                                      myCoords,
    1668                                      y,
    1669                                      NULL,
    1670                                      (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);
    1671         if (rc == false) {
     1640            //            y->data.F32[index] = logf(robustHistogramVector->data.F32[i]);
     1641        }
     1642
     1643        psPolynomial1D *tmpPoly = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
     1644        // XXX: What about the NULL x argument?
     1645        tmpPoly = psVectorFitPolynomial1D(tmpPoly, NULL, y, NULL);
     1646        if (tmpPoly == NULL) {
    16721647            psLogMsg(__func__, PS_LOG_WARN,
    1673                      "WARNING: failed to minimize with psMinimizeLMChi2().\n");
    1674         }
    1675 
    1676         // XXX: Verify this with IfA
    1677         // XXX: The check on the minimization is better than the difference from myMean.
    1678         //      Do they still want this code?
     1648                     "WARNING: failed fit a 1D polynomial.\n");
     1649        }
     1650        psF32 polyFitSigma = PS_SQRT_F32(-0.5 / tmpPoly->coeff[2]);
     1651        psF32 polyFitMean = tmpPoly->coeff[1] * PS_SQR(polyFitSigma);
     1652        // psF32 polyFitNorm = exp(tmpPoly->coedd[0] + PS_SQR(polyFitMean) / (2.0 * PS_SQR(polyFitSigma)));
    16791653
    16801654        if (stats->options & PS_STAT_ROBUST_MEAN) {
    1681             if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) {
    1682                 psLogMsg(__func__, PS_LOG_WARN,
    1683                          "WARNING: the fitted Gaussian has more than 10%% error for the mean.\n");
    1684                 psLogMsg(__func__, PS_LOG_WARN,
    1685                          "WARNING: Using the calculated mean instead of Gaussian-fitted mean.");
    1686                 //                     "WARNING: Using the calculated mean instead of Gaussian-fitted mean.(calc, fit) is (%f, %f)\n",
    1687                 //                     myMean, myParams->data.F32[0]);
    1688                 stats->robustMean = myMean;
    1689             } else {
    1690                 stats->robustMean = myParams->data.F32[0];
    1691             }
     1655            stats->robustMean = polyFitMean;
    16921656        }
    16931657
    16941658        if (stats->options & PS_STAT_ROBUST_STDEV) {
    1695             if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) {
    1696                 psLogMsg(__func__, PS_LOG_WARN,
    1697                          "WARNING: the fitted Gaussian has more than 10%% error for the stdev.\n");
    1698                 psLogMsg(__func__, PS_LOG_WARN,
    1699                          "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.");
    1700                 //                     "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.(calc, fit) is (%f, %f)\n",
    1701                 //                     myStdev, myParams->data.F32[1]);
    1702                 stats->robustStdev = myStdev;
    1703             } else {
    1704                 stats->robustStdev = myParams->data.F32[1];
    1705             }
    1706         }
    1707         psFree(covar);
    1708         psFree(min);
    1709         psFree(myCoords);
     1659            stats->robustStdev = polyFitSigma;
     1660        }
    17101661        psFree(y);
    1711         psFree(myParams);
     1662        psFree(tmpPoly);
    17121663    }
    17131664

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

    r4029 r4101  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.129 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-05-25 20:26:55 $
     11 *  @version $Revision: 1.130 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-06-03 22:30:20 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    16281628    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
    16291629            (stats->options & PS_STAT_ROBUST_STDEV)) {
    1630 
    1631         // Using the above (myMean, myStdev) as initial estimates, we fit a
    1632         // Gaussian to the robustHistogramVector.
    1633         psImage *covar = NULL;
    1634         /* XXX: Old, remove.
    1635                 psMinimization *min = psMinimizationAlloc(100, 0.1);
    1636                 psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
    1637                 psArray *myCoords = psArrayAlloc(robustHistogramVector->n);
    1638                 psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
    1639                 myParams->data.F32[0] = myMean;
    1640                 myParams->data.F32[1] = myStdev;
    1641                 for (i=0;i<robustHistogramVector->n;i++) {
    1642                     myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
    1643                     ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i;
    1644                     y->data.F32[i] = robustHistogramVector->data.F32[i];
    1645                 }
    1646         */
    1647 
    1648         psMinimization *min = psMinimizationAlloc(100, 0.01);
    1649         psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
    1650         psArray *myCoords = psArrayAlloc(2 * dL + 1);
     1630        // We fit a 1-D polynomial to the data.
     1631        // XXX: Since we are no longer fitting a 1-D Gaussian, we can probably
     1632        // remove some of theabove code that calculated the initial estimate
     1633        // for the mean and sigma.
     1634
    16511635        psVector *y = psVectorAlloc(2 * dL + 1, PS_TYPE_F32);
    1652         myParams->data.F32[0] = myMean;
    1653         myParams->data.F32[1] = myStdev;
    1654 
    16551636        for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
    16561637            int index = i - modeBinNum + dL;
    1657             myCoords->data[index] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
    1658             ((psVector *) (myCoords->data[index]))->data.F32[0] = PS_BIN_MIDPOINT(robustHistogram, i);
     1638            // XXX: Should this be the natural log?
    16591639            y->data.F32[index] = robustHistogramVector->data.F32[i];
    1660         }
    1661 
    1662         // XXX: Must fit a 2-D polynomial, not a Gaussian.  (bug 366)
    1663         psBool rc = psMinimizeLMChi2(min,
    1664                                      NULL,
    1665                                      myParams,
    1666                                      NULL,
    1667                                      myCoords,
    1668                                      y,
    1669                                      NULL,
    1670                                      (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);
    1671         if (rc == false) {
     1640            //            y->data.F32[index] = logf(robustHistogramVector->data.F32[i]);
     1641        }
     1642
     1643        psPolynomial1D *tmpPoly = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
     1644        // XXX: What about the NULL x argument?
     1645        tmpPoly = psVectorFitPolynomial1D(tmpPoly, NULL, y, NULL);
     1646        if (tmpPoly == NULL) {
    16721647            psLogMsg(__func__, PS_LOG_WARN,
    1673                      "WARNING: failed to minimize with psMinimizeLMChi2().\n");
    1674         }
    1675 
    1676         // XXX: Verify this with IfA
    1677         // XXX: The check on the minimization is better than the difference from myMean.
    1678         //      Do they still want this code?
     1648                     "WARNING: failed fit a 1D polynomial.\n");
     1649        }
     1650        psF32 polyFitSigma = PS_SQRT_F32(-0.5 / tmpPoly->coeff[2]);
     1651        psF32 polyFitMean = tmpPoly->coeff[1] * PS_SQR(polyFitSigma);
     1652        // psF32 polyFitNorm = exp(tmpPoly->coedd[0] + PS_SQR(polyFitMean) / (2.0 * PS_SQR(polyFitSigma)));
    16791653
    16801654        if (stats->options & PS_STAT_ROBUST_MEAN) {
    1681             if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) {
    1682                 psLogMsg(__func__, PS_LOG_WARN,
    1683                          "WARNING: the fitted Gaussian has more than 10%% error for the mean.\n");
    1684                 psLogMsg(__func__, PS_LOG_WARN,
    1685                          "WARNING: Using the calculated mean instead of Gaussian-fitted mean.");
    1686                 //                     "WARNING: Using the calculated mean instead of Gaussian-fitted mean.(calc, fit) is (%f, %f)\n",
    1687                 //                     myMean, myParams->data.F32[0]);
    1688                 stats->robustMean = myMean;
    1689             } else {
    1690                 stats->robustMean = myParams->data.F32[0];
    1691             }
     1655            stats->robustMean = polyFitMean;
    16921656        }
    16931657
    16941658        if (stats->options & PS_STAT_ROBUST_STDEV) {
    1695             if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) {
    1696                 psLogMsg(__func__, PS_LOG_WARN,
    1697                          "WARNING: the fitted Gaussian has more than 10%% error for the stdev.\n");
    1698                 psLogMsg(__func__, PS_LOG_WARN,
    1699                          "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.");
    1700                 //                     "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.(calc, fit) is (%f, %f)\n",
    1701                 //                     myStdev, myParams->data.F32[1]);
    1702                 stats->robustStdev = myStdev;
    1703             } else {
    1704                 stats->robustStdev = myParams->data.F32[1];
    1705             }
    1706         }
    1707         psFree(covar);
    1708         psFree(min);
    1709         psFree(myCoords);
     1659            stats->robustStdev = polyFitSigma;
     1660        }
    17101661        psFree(y);
    1711         psFree(myParams);
     1662        psFree(tmpPoly);
    17121663    }
    17131664
Note: See TracChangeset for help on using the changeset viewer.