IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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.

File:
1 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);
Note: See TracChangeset for help on using the changeset viewer.