IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 11755


Ignore:
Timestamp:
Feb 12, 2007, 2:19:51 PM (19 years ago)
Author:
Paul Price
Message:

Can really only code a single rejection iteration in psClip, because the user needs to be able to re-run his analysis to produce a new list of values to clip at each iteration. Moved the min-max rejection into a separate function.

Location:
trunk/psLib/src/math
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r11619 r11755  
    1515#include "psClip.h"
    1616
    17 // No-op; purpose of function is to identify the type
     17// No-op; purpose of function is merely to identify the type
    1818static void clipParamsFree(psClipParams *params)
    1919{
    … …  
    4444
    4545
    46 long psClip(psClipParams *params, psVector *values, psVector *mask, const psVector *errors)
     46long psClipMinMax(const psClipParams *params, const psVector *values, psVector *mask)
     47{
     48    PS_ASSERT_VECTOR_NON_NULL(values, -1);
     49    PS_ASSERT_VECTOR_NON_NULL(mask, -1);
     50    PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_MASK, -1);
     51    PS_ASSERT_PTR(params, -1);
     52    PS_ASSERT_VECTORS_SIZE_EQUAL(values, mask, -1);
     53    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(params->fracLow, 0.0, -1);
     54    PS_ASSERT_FLOAT_LESS_THAN(params->fracLow, 1.0, -1);
     55    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(params->fracHigh, 0.0, -1);
     56    PS_ASSERT_FLOAT_LESS_THAN(params->fracHigh, 1.0, -1);
     57    PS_ASSERT_INT_NONZERO(params->clipped, -1);
     58
     59    if (params->fracLow == 0.0 && params->fracHigh == 0.0) {
     60        // No min-max rejection desired
     61        return 0;
     62    }
     63
     64    psMaskType masked = params->masked; // Indicates masked values
     65    psMaskType clipped = params->clipped; // Indicates clipped values
     66    masked |= clipped;                  // Make sure we're also masking clipped values
     67    psMaskType *maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference mask
     68    long totalMasked = 0;               // Total number of pixels masked
     69
     70    // Apply fracLow,fracHigh if there are enough values
     71
     72    // Run through to get number of operational values
     73    long numValid = 0;                  // Number of valid values
     74    for (long i = 0; i < mask->n; i++) {
     75        if (!(maskData[i] & masked)) {
     76            numValid++;
     77        }
     78    }
     79    psTrace("psLib.math", 2, "%ld valid values.\n", numValid);
     80
     81    // XXX: Not sure if sorting provides the fastest implementation.  It might be quicker to do a linear pass
     82    // through the data, pulling out the highest M and lowest N values,
     83
     84    float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of values to keep
     85    if (numValid * keepFrac >= params->numKeep) {
     86        psTrace("psLib.math", 1, "Applying min/max clipping.\n");
     87        psVector *index = psVectorSortIndex(NULL, values);
     88        int numLow = numValid * params->fracLow; // Number of low values to clip
     89        int numHigh = numValid * params->fracHigh; // Number of high values to clip
     90        // Low values
     91        psS32 *indexData = index->data.S32; // Dereference index
     92        long numMasked = 0;             // Number masked
     93        for (long i = 0; i < index->n && numMasked < numLow; i++) {
     94            // Don't count the ones that are already masked
     95            if (! (maskData[indexData[i]] & masked)) {
     96                maskData[indexData[i]] |= clipped;
     97                numMasked++;
     98            }
     99        }
     100        totalMasked += numMasked;
     101        numMasked = 0;
     102        // High values
     103        for (long i = values->n - 1;  i >= 0 && numMasked < numHigh; i--) {
     104            // Don't count the ones that are already masked
     105            if (! (maskData[indexData[i]] & masked)) {
     106                maskData[indexData[i]] |= clipped;
     107                numMasked++;
     108            }
     109        }
     110        totalMasked += numMasked;
     111        psFree(index);
     112    }
     113
     114    return totalMasked;
     115}
     116
     117
     118long psClipReject(psClipParams *params, const psVector *values, psVector *mask, const psVector *errors)
    47119{
    48120    PS_ASSERT_VECTOR_NON_NULL(values, -1);
    … …  
    57129        PS_ASSERT_VECTOR_TYPE_EQUAL(values, errors, -1);
    58130    }
    59     if (params->iter > 0 && (!isfinite(params->rej) || params->rej <= 0)) {
     131    PS_ASSERT_INT_NONZERO(params->meanStat, -1);
     132    PS_ASSERT_INT_NONZERO(params->stdevStat, -1);
     133    if (!isfinite(params->rej) || params->rej <= 0) {
    60134        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Rejection limit (%f) is not valid.\n", params->rej);
    61135        return -1;
    … …  
    65139    psMaskType clipped = params->clipped; // Indicates clipped values
    66140    masked |= clipped;                  // Make sure we're also masking clipped values
    67     psMaskType *maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference mask
    68     long totalMasked = 0;               // Total number of pixels masked
    69141
    70     // Immediate min/max rejection: apply fracLow,fracHigh if there are enough values
    71     if (params->fracLow != 0.0 || params->fracHigh != 0.0) {
    72         // Run through to get number of operational values
    73         long numValid = 0;                  // Number of valid values
    74         for (long i = 0; i < mask->n; i++) {
    75             if (!(maskData[i] & masked)) {
    76                 numValid++;
    77             }
    78         }
    79         psTrace("psLib.math", 2, "%ld valid values.\n", numValid);
    80 
    81         float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of values to keep
    82         if (numValid * keepFrac >= params->numKeep) {
    83             psTrace("psLib.math", 1, "Applying min/max clipping.\n");
    84             psVector *index = psVectorSortIndex(NULL, values);
    85             int numLow = numValid * params->fracLow; // Number of low values to clip
    86             int numHigh = numValid * params->fracHigh; // Number of high values to clip
    87             // Low values
    88             psS32 *indexData = index->data.S32; // Dereference index
    89             long numMasked = 0;             // Number masked
    90             for (long i = 0; i < index->n && numMasked < numLow; i++) {
    91                 // Don't count the ones that are already masked
    92                 if (! (maskData[indexData[i]] & masked)) {
    93                     maskData[indexData[i]] |= clipped;
    94                     numMasked++;
    95                 }
    96             }
    97             totalMasked += numMasked;
    98             numMasked = 0;
    99             // High values
    100             for (long i = values->n - 1;  i >= 0 && numMasked < numHigh; i--) {
    101                 // Don't count the ones that are already masked
    102                 if (! (maskData[indexData[i]] & masked)) {
    103                     maskData[indexData[i]] |= clipped;
    104                     numMasked++;
    105                 }
    106             }
    107             totalMasked += numMasked;
    108             psFree(index);
    109 
    110             // Turn off min/max rejection so that future calls will not trigger it
    111             params->fracHigh = 0.0;
    112             params->fracLow = 0.0;
    113         }
    114     }
    115 
    116     // Iterate, with rejection
     142    // Get statistics
    117143    psStats *stats = psStatsAlloc(params->meanStat | params->stdevStat);
    118144    if (!psVectorStats(stats, values, errors, mask, masked)) {
    … …  
    124150    params->stdev = psStatsGetValue(stats, params->stdevStat);
    125151
    126 
    127152    #define REJECT_CASE(TYPE) \
    128153case PS_TYPE_##TYPE: { \
    129154        ps##TYPE *valuesData = values->data.TYPE; /* Dereference for speed */ \
    130         for (int i = 0; i < params->iter; i++) { \
    131             if (errors) { \
    132                 ps##TYPE *errorsData = errors->data.TYPE; \
    133                 for (int j = 0; j < values->n; j++) { \
    134                     if (!(maskData[j] & masked) && \
    135                             fabs(valuesData[j] - params->mean) > params->rej * errorsData[j]) { \
    136                         maskData[j] |= clipped; \
    137                         totalMasked++; \
    138                     } \
    139                 } \
    140             } else { \
    141                 ps##TYPE limit = params->rej * params->stdev; /* Limit on deviation */ \
    142                 for (int j = 0; j < values->n; j++) { \
    143                     if (!(maskData[j] & masked) && fabs(valuesData[j] - params->mean) > limit) { \
    144                         maskData[j] |= clipped; \
    145                         totalMasked++; \
    146                     } \
     155        psMaskType *maskData = mask->data.PS_TYPE_MASK_DATA; /* Dereference mask for speed */ \
     156        if (errors) { \
     157            ps##TYPE *errorsData = errors->data.TYPE; \
     158            for (int j = 0; j < values->n; j++) { \
     159                if (!(maskData[j] & masked) && \
     160                        fabs(valuesData[j] - params->mean) > params->rej * errorsData[j]) { \
     161                    maskData[j] |= clipped; \
     162                    totalMasked++; \
    147163                } \
    148164            } \
    149             \
    150             if (!psVectorStats(stats, values, errors, mask, masked)) { \
    151                 psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on values.\n"); \
    152                 psFree(stats); \
    153                 return -1; \
     165        } else { \
     166            ps##TYPE limit = params->rej * params->stdev; /* Limit on deviation */ \
     167            for (int j = 0; j < values->n; j++) { \
     168                if (!(maskData[j] & masked) && fabs(valuesData[j] - params->mean) > limit) { \
     169                    maskData[j] |= clipped; \
     170                    totalMasked++; \
     171                } \
    154172            } \
    155             params->mean = psStatsGetValue(stats, params->meanStat); \
    156             params->stdev = psStatsGetValue(stats, params->stdevStat); \
    157173        } \
    158174    }
    159175
     176    long totalMasked = 0;               // Total number of pixels masked
    160177    switch (values->type.type) {
    161178        REJECT_CASE(S8);

  • trunk/psLib/src/math/psClip.h

    r11248 r11755  
    22 * @brief vector clipping functions
    33 *
    4  * $Revision: 1.2 $ $Name: not supported by cvs2svn $
    5  * $Date: 2007-01-23 22:47:23 $
    6  * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     4 * @author Paul Price, IfA.
     5 *
     6 * $Revision: 1.3 $ $Name: not supported by cvs2svn $
     7 * $Date: 2007-02-13 00:19:51 $
     8 * Copyright 2007 Institute for Astronomy, University of Hawaii
    79 */
    810
    … …  
    2123    float fracLow;                      ///< Fraction of low values to clip
    2224    int numKeep;                        ///< Minimum number of values to keep from clipping
    23     int iter;                           ///< Number of rejection iterations
     25    int iter;                           ///< Number of rejection iterations; unused by psClip functions
    2426    float rej;                          ///< Rejection limit (standard deviations)
    2527    psMaskType masked;                  ///< Mask value for entries already masked
    2628    psMaskType clipped;                 ///< Mask value to give to clipped entries
    27     float mean;                         ///< Resultant mean
    28     float stdev;                        ///< Resultant stdev
     29    double mean;                        ///< Resultant mean
     30    double stdev;                       ///< Resultant stdev
    2931}
    3032psClipParams;
    3133
    3234
    33 long psClip(psClipParams *params,       ///< Clip parameters
    34             psVector *values,           ///< Values to inspect and clip
    35             psVector *mask,             ///< Mask for values
    36             const psVector *errors      ///< Errors for values
    37            );
     35
     36/// Apply min-max clipping to a list of values
     37///
     38/// The specified fraction of high and low values are identified as clipped in the mask.  Errors are not used
     39/// in this step.
     40long psClipMinMax(const psClipParams *params, ///< Clip parameters
     41                  const psVector *values, ///< Values to inspect and clip
     42                  psVector *mask        ///< Mask for values
     43    );
     44
     45
     46
     47/// Apply a rejection iteration to a list of values
     48///
     49/// The specified rejection limit is applied to the values and errors; discrepant values are identified as
     50/// clipped in the mask.  This function only applies a single rejection iteration.
     51long psClipReject(psClipParams *params, ///< Clip parameters
     52                  const psVector *values, ///< Values to inspect and clip
     53                  psVector *mask,       ///< Mask for values
     54                  const psVector *errors ///< Errors for values, or NULL
     55    );
    3856
    3957/// @}
Note: See TracChangeset for help on using the changeset viewer.