Index: /trunk/psLib/src/math/psClip.c =================================================================== --- /trunk/psLib/src/math/psClip.c (revision 11754) +++ /trunk/psLib/src/math/psClip.c (revision 11755) @@ -15,5 +15,5 @@ #include "psClip.h" -// No-op; purpose of function is to identify the type +// No-op; purpose of function is merely to identify the type static void clipParamsFree(psClipParams *params) { @@ -44,5 +44,77 @@ -long psClip(psClipParams *params, psVector *values, psVector *mask, const psVector *errors) +long psClipMinMax(const psClipParams *params, const psVector *values, psVector *mask) +{ + PS_ASSERT_VECTOR_NON_NULL(values, -1); + PS_ASSERT_VECTOR_NON_NULL(mask, -1); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_MASK, -1); + PS_ASSERT_PTR(params, -1); + PS_ASSERT_VECTORS_SIZE_EQUAL(values, mask, -1); + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(params->fracLow, 0.0, -1); + PS_ASSERT_FLOAT_LESS_THAN(params->fracLow, 1.0, -1); + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(params->fracHigh, 0.0, -1); + PS_ASSERT_FLOAT_LESS_THAN(params->fracHigh, 1.0, -1); + PS_ASSERT_INT_NONZERO(params->clipped, -1); + + if (params->fracLow == 0.0 && params->fracHigh == 0.0) { + // No min-max rejection desired + return 0; + } + + psMaskType masked = params->masked; // Indicates masked values + psMaskType clipped = params->clipped; // Indicates clipped values + masked |= clipped; // Make sure we're also masking clipped values + psMaskType *maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference mask + long totalMasked = 0; // Total number of pixels masked + + // Apply fracLow,fracHigh if there are enough values + + // Run through to get number of operational values + long numValid = 0; // Number of valid values + for (long i = 0; i < mask->n; i++) { + if (!(maskData[i] & masked)) { + numValid++; + } + } + psTrace("psLib.math", 2, "%ld valid values.\n", numValid); + + // XXX: Not sure if sorting provides the fastest implementation. It might be quicker to do a linear pass + // through the data, pulling out the highest M and lowest N values, + + float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of values to keep + if (numValid * keepFrac >= params->numKeep) { + psTrace("psLib.math", 1, "Applying min/max clipping.\n"); + psVector *index = psVectorSortIndex(NULL, values); + int numLow = numValid * params->fracLow; // Number of low values to clip + int numHigh = numValid * params->fracHigh; // Number of high values to clip + // Low values + psS32 *indexData = index->data.S32; // Dereference index + long numMasked = 0; // Number masked + for (long i = 0; i < index->n && numMasked < numLow; i++) { + // Don't count the ones that are already masked + if (! (maskData[indexData[i]] & masked)) { + maskData[indexData[i]] |= clipped; + numMasked++; + } + } + totalMasked += numMasked; + numMasked = 0; + // High values + for (long i = values->n - 1; i >= 0 && numMasked < numHigh; i--) { + // Don't count the ones that are already masked + if (! (maskData[indexData[i]] & masked)) { + maskData[indexData[i]] |= clipped; + numMasked++; + } + } + totalMasked += numMasked; + psFree(index); + } + + return totalMasked; +} + + +long psClipReject(psClipParams *params, const psVector *values, psVector *mask, const psVector *errors) { PS_ASSERT_VECTOR_NON_NULL(values, -1); @@ -57,5 +129,7 @@ PS_ASSERT_VECTOR_TYPE_EQUAL(values, errors, -1); } - if (params->iter > 0 && (!isfinite(params->rej) || params->rej <= 0)) { + PS_ASSERT_INT_NONZERO(params->meanStat, -1); + PS_ASSERT_INT_NONZERO(params->stdevStat, -1); + if (!isfinite(params->rej) || params->rej <= 0) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Rejection limit (%f) is not valid.\n", params->rej); return -1; @@ -65,54 +139,6 @@ psMaskType clipped = params->clipped; // Indicates clipped values masked |= clipped; // Make sure we're also masking clipped values - psMaskType *maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference mask - long totalMasked = 0; // Total number of pixels masked - // Immediate min/max rejection: apply fracLow,fracHigh if there are enough values - if (params->fracLow != 0.0 || params->fracHigh != 0.0) { - // Run through to get number of operational values - long numValid = 0; // Number of valid values - for (long i = 0; i < mask->n; i++) { - if (!(maskData[i] & masked)) { - numValid++; - } - } - psTrace("psLib.math", 2, "%ld valid values.\n", numValid); - - float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of values to keep - if (numValid * keepFrac >= params->numKeep) { - psTrace("psLib.math", 1, "Applying min/max clipping.\n"); - psVector *index = psVectorSortIndex(NULL, values); - int numLow = numValid * params->fracLow; // Number of low values to clip - int numHigh = numValid * params->fracHigh; // Number of high values to clip - // Low values - psS32 *indexData = index->data.S32; // Dereference index - long numMasked = 0; // Number masked - for (long i = 0; i < index->n && numMasked < numLow; i++) { - // Don't count the ones that are already masked - if (! (maskData[indexData[i]] & masked)) { - maskData[indexData[i]] |= clipped; - numMasked++; - } - } - totalMasked += numMasked; - numMasked = 0; - // High values - for (long i = values->n - 1; i >= 0 && numMasked < numHigh; i--) { - // Don't count the ones that are already masked - if (! (maskData[indexData[i]] & masked)) { - maskData[indexData[i]] |= clipped; - numMasked++; - } - } - totalMasked += numMasked; - psFree(index); - - // Turn off min/max rejection so that future calls will not trigger it - params->fracHigh = 0.0; - params->fracLow = 0.0; - } - } - - // Iterate, with rejection + // Get statistics psStats *stats = psStatsAlloc(params->meanStat | params->stdevStat); if (!psVectorStats(stats, values, errors, mask, masked)) { @@ -124,38 +150,29 @@ params->stdev = psStatsGetValue(stats, params->stdevStat); - #define REJECT_CASE(TYPE) \ case PS_TYPE_##TYPE: { \ ps##TYPE *valuesData = values->data.TYPE; /* Dereference for speed */ \ - for (int i = 0; i < params->iter; i++) { \ - if (errors) { \ - ps##TYPE *errorsData = errors->data.TYPE; \ - for (int j = 0; j < values->n; j++) { \ - if (!(maskData[j] & masked) && \ - fabs(valuesData[j] - params->mean) > params->rej * errorsData[j]) { \ - maskData[j] |= clipped; \ - totalMasked++; \ - } \ - } \ - } else { \ - ps##TYPE limit = params->rej * params->stdev; /* Limit on deviation */ \ - for (int j = 0; j < values->n; j++) { \ - if (!(maskData[j] & masked) && fabs(valuesData[j] - params->mean) > limit) { \ - maskData[j] |= clipped; \ - totalMasked++; \ - } \ + psMaskType *maskData = mask->data.PS_TYPE_MASK_DATA; /* Dereference mask for speed */ \ + if (errors) { \ + ps##TYPE *errorsData = errors->data.TYPE; \ + for (int j = 0; j < values->n; j++) { \ + if (!(maskData[j] & masked) && \ + fabs(valuesData[j] - params->mean) > params->rej * errorsData[j]) { \ + maskData[j] |= clipped; \ + totalMasked++; \ } \ } \ - \ - if (!psVectorStats(stats, values, errors, mask, masked)) { \ - psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on values.\n"); \ - psFree(stats); \ - return -1; \ + } else { \ + ps##TYPE limit = params->rej * params->stdev; /* Limit on deviation */ \ + for (int j = 0; j < values->n; j++) { \ + if (!(maskData[j] & masked) && fabs(valuesData[j] - params->mean) > limit) { \ + maskData[j] |= clipped; \ + totalMasked++; \ + } \ } \ - params->mean = psStatsGetValue(stats, params->meanStat); \ - params->stdev = psStatsGetValue(stats, params->stdevStat); \ } \ } + long totalMasked = 0; // Total number of pixels masked switch (values->type.type) { REJECT_CASE(S8); Index: /trunk/psLib/src/math/psClip.h =================================================================== --- /trunk/psLib/src/math/psClip.h (revision 11754) +++ /trunk/psLib/src/math/psClip.h (revision 11755) @@ -2,7 +2,9 @@ * @brief vector clipping functions * - * $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * $Date: 2007-01-23 22:47:23 $ - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + * @author Paul Price, IfA. + * + * $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * $Date: 2007-02-13 00:19:51 $ + * Copyright 2007 Institute for Astronomy, University of Hawaii */ @@ -21,19 +23,35 @@ float fracLow; ///< Fraction of low values to clip int numKeep; ///< Minimum number of values to keep from clipping - int iter; ///< Number of rejection iterations + int iter; ///< Number of rejection iterations; unused by psClip functions float rej; ///< Rejection limit (standard deviations) psMaskType masked; ///< Mask value for entries already masked psMaskType clipped; ///< Mask value to give to clipped entries - float mean; ///< Resultant mean - float stdev; ///< Resultant stdev + double mean; ///< Resultant mean + double stdev; ///< Resultant stdev } psClipParams; -long psClip(psClipParams *params, ///< Clip parameters - psVector *values, ///< Values to inspect and clip - psVector *mask, ///< Mask for values - const psVector *errors ///< Errors for values - ); + +/// Apply min-max clipping to a list of values +/// +/// The specified fraction of high and low values are identified as clipped in the mask. Errors are not used +/// in this step. +long psClipMinMax(const psClipParams *params, ///< Clip parameters + const psVector *values, ///< Values to inspect and clip + psVector *mask ///< Mask for values + ); + + + +/// Apply a rejection iteration to a list of values +/// +/// The specified rejection limit is applied to the values and errors; discrepant values are identified as +/// clipped in the mask. This function only applies a single rejection iteration. +long psClipReject(psClipParams *params, ///< Clip parameters + const psVector *values, ///< Values to inspect and clip + psVector *mask, ///< Mask for values + const psVector *errors ///< Errors for values, or NULL + ); /// @}