IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 30, 2009, 5:20:29 PM (17 years ago)
Author:
watersc1
Message:

Finished up my edits to the detrend cleanup, and some changes to my
copy of burntool and the pslib astrometry. Detrend cleanup has not
been tested yet. That's up next.

Location:
branches/czw_branch/cleanup
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/cleanup

  • branches/czw_branch/cleanup/psLib/src/fits/psFitsImage.c

    r22729 r24951  
    439439
    440440    p_psFitsReadInfo *info = p_psFitsReadInfoAlloc(fits, region, z);
     441    if (!info) {
     442        psError(PS_ERR_IO, false, "Unable to read FITS information");
     443        return NULL;
     444    }
    441445
    442446    // Size of image

  • branches/czw_branch/cleanup/psLib/src/fits/psFitsScale.c

    r23259 r24951  
    4646    psAssert(image, "impossible");
    4747    psAssert(options, "impossible");
     48
     49    psTrace("psLib.fits", 3, "Scaling image to preserve dynamic range");
    4850
    4951    double range = pow(2.0, options->bitpix) - 1.0; // Range of values for target BITPIX, reduced by the BLANK
    … …  
    109111    psAssert(options, "impossible");
    110112
     113    psTrace("psLib.fits", 3, "Scaling image by statistics");
     114
    111115    // Measure the mean and stdev
    112116    // psImageBackground automatically excludes pixels that are non-finite, so we don't need to bother about a
    … …  
    131135    }
    132136
     137    psTrace("psLib.fits", 5, "Mean: %lf Stdev: %lf", mean, stdev);
     138
    133139    long range = 1 << options->stdevBits;  // Range of values to carry standard deviation
    134140    *bscale = stdev / (double) range;
     141
     142    psTrace("psLib.fits", 5, "Number of bits: %ld BSCALE: %lf", range, *bscale);
    135143
    136144    double imageVal;                    // Value on image
    … …  
    157165
    158166    *bzero = imageVal - *bscale * diskVal;
     167
     168    psTrace("psLib.fits", 5, "Image %lf corresponds to disk %ld --> BZERO: %lf", imageVal, diskVal, *bzero);
    159169
    160170    return true;

  • branches/czw_branch/cleanup/psLib/src/imageops/psImageCovariance.c

    r23882 r24951  
    282282}
    283283
     284psKernel *psImageCovarianceAverageWeighted(const psArray *array, const psVector *weights)
     285{
     286    PS_ASSERT_ARRAY_NON_NULL(array, NULL);
     287    PS_ASSERT_ARRAY_NON_EMPTY(array, NULL);
     288    if (!weights) {
     289        return psImageCovarianceAverage(array);
     290    }
     291    PS_ASSERT_VECTOR_TYPE(weights, PS_TYPE_F32, NULL);
     292
     293    int xMin = INT_MAX, xMax = INT_MIN, yMin = INT_MAX, yMax = INT_MIN; // Range for covariance
     294    double sumWeights = 0.0;            // Sum of weights
     295    for (int i = 0; i < array->n; i++) {
     296        psKernel *covar = array->data[i]; // Covariance matrix
     297        if (!covar) {
     298            continue;
     299        }
     300        xMin = PS_MIN(xMin, covar->xMin);
     301        xMax = PS_MAX(xMax, covar->xMax);
     302        yMin = PS_MIN(yMin, covar->yMin);
     303        yMax = PS_MAX(yMax, covar->yMax);
     304        sumWeights += weights->data.F32[i];
     305    }
     306    if (sumWeights == 0) {
     307        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No covariance matrices supplied for summation");
     308        return NULL;
     309    }
     310
     311    psKernel *sum = psKernelAlloc(xMin, xMax, yMin, yMax); // Summed covariance
     312    for (int i = 0; i < array->n; i++) {
     313        psKernel *covar = array->data[i]; // Covariance matrix
     314        if (!covar) {
     315            continue;
     316        }
     317        for (int y = covar->yMin; y <= covar->yMax; y++) {
     318            for (int x = covar->xMin; x <= covar->xMax; x++) {
     319                if (!isfinite(covar->kernel[y][x])) {
     320                    psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     321                            "Non-finite covariance matrix element at %d,%d for input %d",
     322                            x, y, i);
     323                    psFree(sum);
     324                    return NULL;
     325                }
     326                sum->kernel[y][x] += weights->data.F32[i] * covar->kernel[y][x];
     327            }
     328        }
     329    }
     330    psBinaryOp(sum->image, sum->image, "/", psScalarAlloc((float)sumWeights, PS_TYPE_F32));
     331
     332    return sum;
     333}
     334
    284335
    285336psKernel *psImageCovarianceTruncate(const psKernel *covar, float frac)

  • branches/czw_branch/cleanup/psLib/src/imageops/psImageCovariance.h

    r23882 r24951  
    5757    );
    5858
     59/// Weighted average of multiple covariance pseudo-matrices
     60psKernel *psImageCovarianceAverageWeighted(
     61    const psArray *array,               ///< Array of covariance pseudo-matrices
     62    const psVector *weights             ///< Weights for each (F32)
     63    );
     64
    5965/// Truncate covariance pseudo-matrix
    6066///

  • branches/czw_branch/cleanup/psLib/src/math/psPolynomialMD.c

    r24090 r24951  
    321321}
    322322
    323 // XXX this function should take a (psVectorMaskType markVal) argument
    324323bool psPolynomialMDClipFit(psPolynomialMD *poly, const psVector *values, const psVector *errors,
    325324                           const psVector *mask, psVectorMaskType maskVal, const psArray *coordsArray,

  • branches/czw_branch/cleanup/psLib/src/math/psStats.c

    r24092 r24951  
    128128        RESULT = Xt; }
    129129
     130# define COUNT_WARNING(LIMIT, INTERVAL, ...) { \
     131        static int nCalls = 1; \
     132        if (nCalls < LIMIT) { \
     133            psWarning(__VA_ARGS__); \
     134        } \
     135        if (!(nCalls % INTERVAL)) { \
     136            psWarning(__VA_ARGS__); \
     137            psWarning("(warning raised %d times)", nCalls); \
     138        } \
     139        nCalls ++; \
     140}
     141 
     142
    130143/*****************************************************************************/
    131144/* TYPE DEFINITIONS                                                          */
    … …  
    331344
    332345    if (count == 0) {
    333         psLogMsg(TRACE, PS_LOG_DETAIL, "No valid data in input vector.\n");
     346        COUNT_WARNING(10, 100, "No valid data in input vector.\n");
    334347        stats->sampleUQ = NAN;
    335348        stats->sampleLQ = NAN;
    … …  
    394407    // If the mean is NAN, then generate a warning and set the stdev to NAN.
    395408    if (isnan(stats->sampleMean)) {
    396         psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN.\n");
     409        COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN.");
    397410        stats->sampleStdev = NAN;
    398411        return true;
    … …  
    438451        // Assume that the user knows what he's doing when he masks out everything --> no error.
    439452        stats->sampleStdev = NAN;
    440         psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleStdev(): no valid psVector elements (%ld). Setting stats->sampleStdev = NAN.\n", count);
     453        COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): no valid psVector elements (%ld). Setting stats->sampleStdev = NAN.\n", count);
    441454        return true;
    442455    }
    443456    if (count == 1) {
    444457        stats->sampleStdev = 0.0;
    445         psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld).  "
    446                 "Setting stats->sampleStdev = 0.0.\n", count);
     458        COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld). Setting stats->sampleStdev = 0.0.\n", count);
    447459        return true;
    448460    }
    … …  
    468480    }
    469481    if (isnan(stats->sampleMean)) {
    470         psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleMoments(): sample mean is NAN.\n");
     482        COUNT_WARNING(10, 100, "WARNING: vectorSampleMoments(): sample mean is NAN.\n");
    471483        goto SAMPLE_MOMENTS_BAD;
    472484    }
    … …  
    475487    }
    476488    if (isnan(stats->sampleStdev) || stats->sampleStdev == 0.0) {
    477         psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleMoments(): sample stdev is NAN or 0.\n");
     489        COUNT_WARNING(10, 100, "WARNING: vectorSampleMoments(): sample stdev is NAN or 0.\n");
    478490        goto SAMPLE_MOMENTS_BAD;
    479491    }
    … …  
    583595    vectorSampleMedian(myVector, tmpMask, maskVal, stats);
    584596    if (isnan(stats->sampleMedian)) {
    585         psLogMsg(TRACE, PS_LOG_DETAIL, "Call to vectorSampleMedian returned NAN\n");
     597        COUNT_WARNING(10, 100, "Call to vectorSampleMedian returned NAN\n");
    586598        return true;
    587599    }
    … …  
    591603    vectorSampleStdev(myVector, errors, tmpMask, maskVal, stats);
    592604    if (isnan(stats->sampleStdev)) {
    593         psLogMsg(TRACE, PS_LOG_DETAIL, "Call to vectorSampleStdev returned NAN\n");
     605        COUNT_WARNING(10, 100, "Call to vectorSampleStdev returned NAN\n");
    594606        return true;
    595607    }
    … …  
    649661        if (isnan(stats->sampleMean) || isnan(stats->sampleStdev)) {
    650662            iter = stats->clipIter;
    651             psLogMsg(TRACE, PS_LOG_DETAIL, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n");
     663            COUNT_WARNING(10, 100, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n");
    652664            clippedMean = NAN;
    653665            clippedStdev = NAN;
    … …  
    747759        if (numValid == 0 || isnan(min) || isnan(max)) {
    748760            // Data range calculation failed
    749             psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n");
     761            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    750762            goto escape;
    751763        }
    … …  
    763775            stats->results |= PS_STAT_ROBUST_STDEV;
    764776            stats->results |= PS_STAT_ROBUST_QUARTILE;
    765             psLogMsg(TRACE, PS_LOG_DETAIL, "All data points have the same value: %f.\n", min);
     777            COUNT_WARNING(10, 100, "All data points have the same value: %f.\n", min);
    766778            psFree(mask);
    767779            psFree(statsMinMax);
    … …  
    832844        // convert bin to bin value: this is the robust histogram median.
    833845        if (isnan(stats->robustMedian)) {
    834             psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to fit a quadratic and calculate the 50-percent position.\n");
     846            COUNT_WARNING(10, 100, "Failed to fit a quadratic and calculate the 50-percent position.\n");
    835847            goto escape;
    836848        }
    … …  
    854866
    855867        if ((binLo < 0) || (binHi < 0)) {
    856             psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n");
     868            COUNT_WARNING(10, 100, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n");
    857869            goto escape;
    858870        }
    … …  
    950962                stats->results |= PS_STAT_ROBUST_STDEV;
    951963                stats->results |= PS_STAT_ROBUST_QUARTILE;
    952                 psLogMsg(TRACE, PS_LOG_DETAIL, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS);
     964                COUNT_WARNING(10, 100, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS);
    953965                psFree(mask);
    954966                psFree(statsMinMax);
    … …  
    978990    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f);
    979991    if (isnan(binLo25F32) || isnan(binHi25F32)) {
    980         psLogMsg(TRACE, PS_LOG_DETAIL, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnBin() returned a NAN.\n");
     992        COUNT_WARNING(10, 100, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnBin() returned a NAN.\n");
    981993        goto escape;
    982994    }
    … …  
    12631275        float max = statsMinMax->max;
    12641276        if (numValid == 0 || isnan(min) || isnan(max)) {
    1265             psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n");
     1277            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    12661278            psFree(statsMinMax);
    12671279            psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    … …  
    14601472        float max = statsMinMax->max;
    14611473        if (numValid == 0 || isnan(min) || isnan(max)) {
    1462             psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n");
     1474            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    14631475            psFree(statsMinMax);
    14641476            psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    … …  
    15081520        PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0);
    15091521        if (binMin == binMax) {
    1510             psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n");
     1522            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    15111523            psFree(statsMinMax);
    15121524            return true;
    … …  
    17221734
    17231735    // If the mean is NAN, then generate a warning and set the stdev to NAN.
    1724     if (isnan(stats->robustMedian)) {
    1725         stats->fittedStdev = NAN;
    1726         stats->fittedStdev = NAN;
    1727         return true;
    1728     }
     1736    if (isnan(stats->robustMedian)) goto escape;
    17291737
    17301738    float guessStdev = stats->robustStdev;  // pass the guess sigma
    … …  
    17561764        float max = statsMinMax->max;
    17571765        if (numValid == 0 || isnan(min) || isnan(max)) {
    1758             psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n");
     1766            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    17591767            psFree(statsMinMax);
    1760             stats->fittedStdev = NAN;
    1761             stats->fittedStdev = NAN;
     1768            goto escape;
     1769        }
     1770
     1771        // If all data points have the same value, then we set the appropriate members of stats and return.
     1772        if (fabs(max - min) <= FLT_EPSILON) {
     1773            COUNT_WARNING(10, 100, "All data points have the same value: %f.\n", min);
     1774            stats->fittedMean = min;
     1775            stats->fittedStdev = 0.0;
     1776            stats->results |= PS_STAT_FITTED_MEAN_V4;
     1777            stats->results |= PS_STAT_FITTED_STDEV_V4;
    17621778            return true;
    17631779        }
    … …  
    17741790        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
    17751791        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
    1776             psLogMsg(TRACE, PS_LOG_DETAIL, "Unable to generate histogram for fitted statistics v4.\n");
     1792            COUNT_WARNING(10, 100, "Unable to generate histogram for fitted statistics v4.\n");
    17771793            psFree(histogram);
    17781794            psFree(statsMinMax);
    1779             stats->fittedStdev = NAN;
    1780             stats->fittedStdev = NAN;
    1781             return true;
     1795            goto escape;
    17821796        }
    17831797        if (psTraceGetLevel("psLib.math") >= 8) {
    … …  
    18091823        PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0);
    18101824        if (binMin == binMax) {
    1811             psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n");
     1825            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    18121826            psFree(statsMinMax);
    1813             stats->fittedStdev = NAN;
    1814             stats->fittedStdev = NAN;
    1815             return true;
     1827            goto escape;
    18161828        }
    18171829
    … …  
    18821894
    18831895            if (!status) {
    1884                 psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to fit a gaussian to the robust histogram.\n");
     1896                COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
    18851897                psFree(poly);
    18861898                psFree(histogram);
    18871899                psFree(statsMinMax);
    1888                 stats->fittedStdev = NAN;
    1889                 stats->fittedStdev = NAN;
    1890                 return true;
     1900                goto escape;
    18911901            }
    18921902
    18931903            if (poly->coeff[2] >= 0.0) {
    1894                 psLogMsg(TRACE, PS_LOG_MINUTIA, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
     1904                COUNT_WARNING(10, 100, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    18951905                psFree(poly);
    18961906                psFree(histogram);
    … …  
    19061916                }
    19071917
    1908                 psLogMsg(TRACE, PS_LOG_DETAIL, "fit did not converge\n");
    1909                 stats->fittedStdev = NAN;
    1910                 stats->fittedStdev = NAN;
    1911                 return true;
     1918                COUNT_WARNING(10, 100, "fit did not converge\n");
     1919                goto escape;
    19121920            }
    19131921
    … …  
    19751983
    19761984            if (!status) {
    1977                 psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to fit a gaussian to the robust histogram.\n");
     1985                COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
    19781986                psFree(poly);
    19791987                psFree(histogram);
    19801988                psFree(statsMinMax);
    1981                 stats->fittedStdev = NAN;
    1982                 stats->fittedStdev = NAN;
    1983                 return true;
     1989                goto escape;
    19841990            }
    19851991
    … …  
    20262032    psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev);
    20272033
     2034    stats->results |= PS_STAT_FITTED_MEAN_V4;
     2035    stats->results |= PS_STAT_FITTED_STDEV_V4;
     2036
     2037    return true;
     2038
     2039escape:
     2040    stats->fittedMean = NAN;
     2041    stats->fittedStdev = NAN;
    20282042    stats->results |= PS_STAT_FITTED_MEAN_V4;
    20292043    stats->results |= PS_STAT_FITTED_STDEV_V4;

  • branches/czw_branch/cleanup/psLib/src/mathtypes/psVector.c

    r24619 r24951  
    122122    }
    123123
    124    
    125124    if (vector->nalloc == nalloc) {     
    126125        // No need to realloc to same size
    … …  
    131130    elementSize = PSELEMTYPE_SIZEOF(elemType);
    132131
    133     long nstart = vector->n;
     132    long nallocOld = vector->nalloc;
    134133    if (nalloc < vector->n) {
    135134        vector->n = nalloc;
    … …  
    139138    P_PSVECTOR_SET_NALLOC(vector,nalloc);
    140139
    141     // fill newly allocated range with zeros:
    142     if (nstart < nalloc) {
    143         long nNew = nalloc - nstart;
    144         memset (&vector->data.U8[nstart*elementSize], 0, nNew*elementSize);
     140    // fill newly allocated range with zeros: 
     141    if (nallocOld < nalloc) {
     142        long nNew = nalloc - nallocOld;
     143        memset (&vector->data.U8[nallocOld*elementSize], 0, nNew*elementSize);
    145144    }
    146145
Note: See TracChangeset for help on using the changeset viewer.