Index: branches/czw_branch/cleanup/psLib/src/fits/psFitsImage.c =================================================================== --- branches/czw_branch/cleanup/psLib/src/fits/psFitsImage.c (revision 24713) +++ branches/czw_branch/cleanup/psLib/src/fits/psFitsImage.c (revision 24951) @@ -439,4 +439,8 @@ p_psFitsReadInfo *info = p_psFitsReadInfoAlloc(fits, region, z); + if (!info) { + psError(PS_ERR_IO, false, "Unable to read FITS information"); + return NULL; + } // Size of image Index: branches/czw_branch/cleanup/psLib/src/fits/psFitsScale.c =================================================================== --- branches/czw_branch/cleanup/psLib/src/fits/psFitsScale.c (revision 24713) +++ branches/czw_branch/cleanup/psLib/src/fits/psFitsScale.c (revision 24951) @@ -46,4 +46,6 @@ psAssert(image, "impossible"); psAssert(options, "impossible"); + + psTrace("psLib.fits", 3, "Scaling image to preserve dynamic range"); double range = pow(2.0, options->bitpix) - 1.0; // Range of values for target BITPIX, reduced by the BLANK @@ -109,4 +111,6 @@ psAssert(options, "impossible"); + psTrace("psLib.fits", 3, "Scaling image by statistics"); + // Measure the mean and stdev // psImageBackground automatically excludes pixels that are non-finite, so we don't need to bother about a @@ -131,6 +135,10 @@ } + psTrace("psLib.fits", 5, "Mean: %lf Stdev: %lf", mean, stdev); + long range = 1 << options->stdevBits; // Range of values to carry standard deviation *bscale = stdev / (double) range; + + psTrace("psLib.fits", 5, "Number of bits: %ld BSCALE: %lf", range, *bscale); double imageVal; // Value on image @@ -157,4 +165,6 @@ *bzero = imageVal - *bscale * diskVal; + + psTrace("psLib.fits", 5, "Image %lf corresponds to disk %ld --> BZERO: %lf", imageVal, diskVal, *bzero); return true; Index: branches/czw_branch/cleanup/psLib/src/imageops/psImageCovariance.c =================================================================== --- branches/czw_branch/cleanup/psLib/src/imageops/psImageCovariance.c (revision 24713) +++ branches/czw_branch/cleanup/psLib/src/imageops/psImageCovariance.c (revision 24951) @@ -282,4 +282,55 @@ } +psKernel *psImageCovarianceAverageWeighted(const psArray *array, const psVector *weights) +{ + PS_ASSERT_ARRAY_NON_NULL(array, NULL); + PS_ASSERT_ARRAY_NON_EMPTY(array, NULL); + if (!weights) { + return psImageCovarianceAverage(array); + } + PS_ASSERT_VECTOR_TYPE(weights, PS_TYPE_F32, NULL); + + int xMin = INT_MAX, xMax = INT_MIN, yMin = INT_MAX, yMax = INT_MIN; // Range for covariance + double sumWeights = 0.0; // Sum of weights + for (int i = 0; i < array->n; i++) { + psKernel *covar = array->data[i]; // Covariance matrix + if (!covar) { + continue; + } + xMin = PS_MIN(xMin, covar->xMin); + xMax = PS_MAX(xMax, covar->xMax); + yMin = PS_MIN(yMin, covar->yMin); + yMax = PS_MAX(yMax, covar->yMax); + sumWeights += weights->data.F32[i]; + } + if (sumWeights == 0) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No covariance matrices supplied for summation"); + return NULL; + } + + psKernel *sum = psKernelAlloc(xMin, xMax, yMin, yMax); // Summed covariance + for (int i = 0; i < array->n; i++) { + psKernel *covar = array->data[i]; // Covariance matrix + if (!covar) { + continue; + } + for (int y = covar->yMin; y <= covar->yMax; y++) { + for (int x = covar->xMin; x <= covar->xMax; x++) { + if (!isfinite(covar->kernel[y][x])) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Non-finite covariance matrix element at %d,%d for input %d", + x, y, i); + psFree(sum); + return NULL; + } + sum->kernel[y][x] += weights->data.F32[i] * covar->kernel[y][x]; + } + } + } + psBinaryOp(sum->image, sum->image, "/", psScalarAlloc((float)sumWeights, PS_TYPE_F32)); + + return sum; +} + psKernel *psImageCovarianceTruncate(const psKernel *covar, float frac) Index: branches/czw_branch/cleanup/psLib/src/imageops/psImageCovariance.h =================================================================== --- branches/czw_branch/cleanup/psLib/src/imageops/psImageCovariance.h (revision 24713) +++ branches/czw_branch/cleanup/psLib/src/imageops/psImageCovariance.h (revision 24951) @@ -57,4 +57,10 @@ ); +/// Weighted average of multiple covariance pseudo-matrices +psKernel *psImageCovarianceAverageWeighted( + const psArray *array, ///< Array of covariance pseudo-matrices + const psVector *weights ///< Weights for each (F32) + ); + /// Truncate covariance pseudo-matrix /// Index: branches/czw_branch/cleanup/psLib/src/math/psPolynomialMD.c =================================================================== --- branches/czw_branch/cleanup/psLib/src/math/psPolynomialMD.c (revision 24713) +++ branches/czw_branch/cleanup/psLib/src/math/psPolynomialMD.c (revision 24951) @@ -321,5 +321,4 @@ } -// XXX this function should take a (psVectorMaskType markVal) argument bool psPolynomialMDClipFit(psPolynomialMD *poly, const psVector *values, const psVector *errors, const psVector *mask, psVectorMaskType maskVal, const psArray *coordsArray, Index: branches/czw_branch/cleanup/psLib/src/math/psStats.c =================================================================== --- branches/czw_branch/cleanup/psLib/src/math/psStats.c (revision 24713) +++ branches/czw_branch/cleanup/psLib/src/math/psStats.c (revision 24951) @@ -128,4 +128,17 @@ RESULT = Xt; } +# define COUNT_WARNING(LIMIT, INTERVAL, ...) { \ + static int nCalls = 1; \ + if (nCalls < LIMIT) { \ + psWarning(__VA_ARGS__); \ + } \ + if (!(nCalls % INTERVAL)) { \ + psWarning(__VA_ARGS__); \ + psWarning("(warning raised %d times)", nCalls); \ + } \ + nCalls ++; \ +} + + /*****************************************************************************/ /* TYPE DEFINITIONS */ @@ -331,5 +344,5 @@ if (count == 0) { - psLogMsg(TRACE, PS_LOG_DETAIL, "No valid data in input vector.\n"); + COUNT_WARNING(10, 100, "No valid data in input vector.\n"); stats->sampleUQ = NAN; stats->sampleLQ = NAN; @@ -394,5 +407,5 @@ // If the mean is NAN, then generate a warning and set the stdev to NAN. if (isnan(stats->sampleMean)) { - psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN.\n"); + COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN."); stats->sampleStdev = NAN; return true; @@ -438,11 +451,10 @@ // Assume that the user knows what he's doing when he masks out everything --> no error. stats->sampleStdev = NAN; - psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleStdev(): no valid psVector elements (%ld). Setting stats->sampleStdev = NAN.\n", count); + COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): no valid psVector elements (%ld). Setting stats->sampleStdev = NAN.\n", count); return true; } if (count == 1) { stats->sampleStdev = 0.0; - psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld). " - "Setting stats->sampleStdev = 0.0.\n", count); + COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld). Setting stats->sampleStdev = 0.0.\n", count); return true; } @@ -468,5 +480,5 @@ } if (isnan(stats->sampleMean)) { - psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleMoments(): sample mean is NAN.\n"); + COUNT_WARNING(10, 100, "WARNING: vectorSampleMoments(): sample mean is NAN.\n"); goto SAMPLE_MOMENTS_BAD; } @@ -475,5 +487,5 @@ } if (isnan(stats->sampleStdev) || stats->sampleStdev == 0.0) { - psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleMoments(): sample stdev is NAN or 0.\n"); + COUNT_WARNING(10, 100, "WARNING: vectorSampleMoments(): sample stdev is NAN or 0.\n"); goto SAMPLE_MOMENTS_BAD; } @@ -583,5 +595,5 @@ vectorSampleMedian(myVector, tmpMask, maskVal, stats); if (isnan(stats->sampleMedian)) { - psLogMsg(TRACE, PS_LOG_DETAIL, "Call to vectorSampleMedian returned NAN\n"); + COUNT_WARNING(10, 100, "Call to vectorSampleMedian returned NAN\n"); return true; } @@ -591,5 +603,5 @@ vectorSampleStdev(myVector, errors, tmpMask, maskVal, stats); if (isnan(stats->sampleStdev)) { - psLogMsg(TRACE, PS_LOG_DETAIL, "Call to vectorSampleStdev returned NAN\n"); + COUNT_WARNING(10, 100, "Call to vectorSampleStdev returned NAN\n"); return true; } @@ -649,5 +661,5 @@ if (isnan(stats->sampleMean) || isnan(stats->sampleStdev)) { iter = stats->clipIter; - psLogMsg(TRACE, PS_LOG_DETAIL, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n"); + COUNT_WARNING(10, 100, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n"); clippedMean = NAN; clippedStdev = NAN; @@ -747,5 +759,5 @@ if (numValid == 0 || isnan(min) || isnan(max)) { // Data range calculation failed - psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n"); + COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); goto escape; } @@ -763,5 +775,5 @@ stats->results |= PS_STAT_ROBUST_STDEV; stats->results |= PS_STAT_ROBUST_QUARTILE; - psLogMsg(TRACE, PS_LOG_DETAIL, "All data points have the same value: %f.\n", min); + COUNT_WARNING(10, 100, "All data points have the same value: %f.\n", min); psFree(mask); psFree(statsMinMax); @@ -832,5 +844,5 @@ // convert bin to bin value: this is the robust histogram median. if (isnan(stats->robustMedian)) { - psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to fit a quadratic and calculate the 50-percent position.\n"); + COUNT_WARNING(10, 100, "Failed to fit a quadratic and calculate the 50-percent position.\n"); goto escape; } @@ -854,5 +866,5 @@ if ((binLo < 0) || (binHi < 0)) { - psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n"); + COUNT_WARNING(10, 100, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n"); goto escape; } @@ -950,5 +962,5 @@ stats->results |= PS_STAT_ROBUST_STDEV; stats->results |= PS_STAT_ROBUST_QUARTILE; - psLogMsg(TRACE, PS_LOG_DETAIL, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS); + COUNT_WARNING(10, 100, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS); psFree(mask); psFree(statsMinMax); @@ -978,5 +990,5 @@ psF32 binHi25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f); if (isnan(binLo25F32) || isnan(binHi25F32)) { - psLogMsg(TRACE, PS_LOG_DETAIL, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnBin() returned a NAN.\n"); + COUNT_WARNING(10, 100, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnBin() returned a NAN.\n"); goto escape; } @@ -1263,5 +1275,5 @@ float max = statsMinMax->max; if (numValid == 0 || isnan(min) || isnan(max)) { - psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n"); + COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); @@ -1460,5 +1472,5 @@ float max = statsMinMax->max; if (numValid == 0 || isnan(min) || isnan(max)) { - psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n"); + COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); @@ -1508,5 +1520,5 @@ PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0); if (binMin == binMax) { - psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n"); + COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); return true; @@ -1722,9 +1734,5 @@ // If the mean is NAN, then generate a warning and set the stdev to NAN. - if (isnan(stats->robustMedian)) { - stats->fittedStdev = NAN; - stats->fittedStdev = NAN; - return true; - } + if (isnan(stats->robustMedian)) goto escape; float guessStdev = stats->robustStdev; // pass the guess sigma @@ -1756,8 +1764,16 @@ float max = statsMinMax->max; if (numValid == 0 || isnan(min) || isnan(max)) { - psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n"); + COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); - stats->fittedStdev = NAN; - stats->fittedStdev = NAN; + goto escape; + } + + // If all data points have the same value, then we set the appropriate members of stats and return. + if (fabs(max - min) <= FLT_EPSILON) { + COUNT_WARNING(10, 100, "All data points have the same value: %f.\n", min); + stats->fittedMean = min; + stats->fittedStdev = 0.0; + stats->results |= PS_STAT_FITTED_MEAN_V4; + stats->results |= PS_STAT_FITTED_STDEV_V4; return true; } @@ -1774,10 +1790,8 @@ psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { - psLogMsg(TRACE, PS_LOG_DETAIL, "Unable to generate histogram for fitted statistics v4.\n"); + COUNT_WARNING(10, 100, "Unable to generate histogram for fitted statistics v4.\n"); psFree(histogram); psFree(statsMinMax); - stats->fittedStdev = NAN; - stats->fittedStdev = NAN; - return true; + goto escape; } if (psTraceGetLevel("psLib.math") >= 8) { @@ -1809,9 +1823,7 @@ PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0); if (binMin == binMax) { - psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n"); + COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); - stats->fittedStdev = NAN; - stats->fittedStdev = NAN; - return true; + goto escape; } @@ -1882,15 +1894,13 @@ if (!status) { - psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to fit a gaussian to the robust histogram.\n"); + COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n"); psFree(poly); psFree(histogram); psFree(statsMinMax); - stats->fittedStdev = NAN; - stats->fittedStdev = NAN; - return true; + goto escape; } if (poly->coeff[2] >= 0.0) { - psLogMsg(TRACE, PS_LOG_MINUTIA, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); + COUNT_WARNING(10, 100, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); psFree(poly); psFree(histogram); @@ -1906,8 +1916,6 @@ } - psLogMsg(TRACE, PS_LOG_DETAIL, "fit did not converge\n"); - stats->fittedStdev = NAN; - stats->fittedStdev = NAN; - return true; + COUNT_WARNING(10, 100, "fit did not converge\n"); + goto escape; } @@ -1975,11 +1983,9 @@ if (!status) { - psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to fit a gaussian to the robust histogram.\n"); + COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n"); psFree(poly); psFree(histogram); psFree(statsMinMax); - stats->fittedStdev = NAN; - stats->fittedStdev = NAN; - return true; + goto escape; } @@ -2026,4 +2032,12 @@ psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev); + stats->results |= PS_STAT_FITTED_MEAN_V4; + stats->results |= PS_STAT_FITTED_STDEV_V4; + + return true; + +escape: + stats->fittedMean = NAN; + stats->fittedStdev = NAN; stats->results |= PS_STAT_FITTED_MEAN_V4; stats->results |= PS_STAT_FITTED_STDEV_V4; Index: branches/czw_branch/cleanup/psLib/src/mathtypes/psVector.c =================================================================== --- branches/czw_branch/cleanup/psLib/src/mathtypes/psVector.c (revision 24713) +++ branches/czw_branch/cleanup/psLib/src/mathtypes/psVector.c (revision 24951) @@ -122,5 +122,4 @@ } - if (vector->nalloc == nalloc) { // No need to realloc to same size @@ -131,5 +130,5 @@ elementSize = PSELEMTYPE_SIZEOF(elemType); - long nstart = vector->n; + long nallocOld = vector->nalloc; if (nalloc < vector->n) { vector->n = nalloc; @@ -139,8 +138,8 @@ P_PSVECTOR_SET_NALLOC(vector,nalloc); - // fill newly allocated range with zeros: - if (nstart < nalloc) { - long nNew = nalloc - nstart; - memset (&vector->data.U8[nstart*elementSize], 0, nNew*elementSize); + // fill newly allocated range with zeros: + if (nallocOld < nalloc) { + long nNew = nalloc - nallocOld; + memset (&vector->data.U8[nallocOld*elementSize], 0, nNew*elementSize); }