IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 7256


Ignore:
Timestamp:
May 31, 2006, 4:26:32 PM (20 years ago)
Author:
Paul Price
Message:

pmReadoutCombine is working, but it's kind of slow.

Location:
trunk/psModules/src/imcombine
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmReadoutCombine.c

    r7254 r7256  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2006-05-31 22:51:15 $
     7 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2006-06-01 02:26:32 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2525//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    2626
    27 #if 1
    2827// Return the statistic of interest
    29 static double getStat(const psStats *stats, // Statistics structure
    30                       psStatsOptions option // Statistics option
    31                      )
     28static inline double getStat(const psStats *stats, // Statistics structure
     29                             psStatsOptions option // Statistics option
     30                            )
    3231{
    3332    switch (option) {
    … …  
    5958    return NAN;
    6059}
    61 #endif
    62 
    63 // Mask for combination --- used only locally
    64 typedef enum {
    65     PM_READOUT_COMBINE_CLEAR      = 0x00, // No problem
    66     PM_READOUT_COMBINE_NO_IMAGE   = 0x01, // No image available
    67     PM_READOUT_COMBINE_MASKED     = 0x02, // Pixel is masked
    68     PM_READOUT_COMBINE_BAD        = 0x03, // Pixel is bad (NO_IMAGE | MASKED)
    69     PM_READOUT_COMBINE_CLIPPED    = 0x04  // Pixel has been clipped
    70 } pmReadoutCombineMask;
    71 
    7260
    7361//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    9684 *****************************************************************************/
    9785bool pmReadoutCombine(pmReadout *output,
    98                       const psArray *inputs,
     86                      psArray *inputs,
    9987                      const psVector *zero,
    10088                      const psVector *scale,
    … …  
    144132    psVector *mask   = psVectorAlloc(inputs->n, PS_TYPE_U8); // Mask for stack
    145133    mask->n = inputs->n;
    146     psVectorInit(mask, 0);
     134
    147135    bool haveWeights = false;           // Do we have weight images?
    148136    bool valid = false;                 // Do we have a single valid input?
    … …  
    150138        pmReadout *readout = inputs->data[i]; // Readout of interest
    151139        if (!readout || !readout->image) {
    152             mask->data.U8[i] = PM_READOUT_COMBINE_NO_IMAGE;
    153             continue;
     140            psError(PS_ERR_UNEXPECTED_NULL, true, "Input readout %d is NULL or has NULL image.\n", i);
     141            return false;
    154142        }
    155143
    … …  
    179167        psError(PS_ERR_IO, true, "No valid inputs.\n");
    180168        return NULL;
     169    }
     170
     171    // Apply scale and zero
     172    if (scale || zero) {
     173        for (int i = 0; i < inputs->n; i++) {
     174            pmReadout *readout = inputs->data[i]; // The particular readout
     175            if (zero) {
     176                psBinaryOp(readout->image, readout->image, "-",
     177                           psScalarAlloc(zero->data.F32[i], PS_TYPE_F32));
     178            }
     179            if (scale) {
     180                float scaleFactor = 1.0 / scale->data.F32[i]; // Scaling
     181                psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(scaleFactor, PS_TYPE_F32));
     182                if (haveWeights) {
     183                    psBinaryOp(readout->weight, readout->weight, "*",
     184                               psScalarAlloc(scaleFactor * scaleFactor, PS_TYPE_F32));
     185                }
     186            }
     187        }
    181188    }
    182189
    … …  
    277284        }
    278285        for (int j = minInputCols; j < maxInputCols; j++) {
    279 
    280286            int numValid = 0;           // Number of valid pixels in the stack
     287            memset(mask->data.U8, 0, mask->n * sizeof(psU8)); // Reset the mask
    281288            for (int r = 0; r < inputs->n; r++) {
    282289                pmReadout *readout = inputs->data[r]; // Input readout
    … …  
    284291                int xIn = j - readout->col0; // x position on input readout
    285292
    286                 // Check existence and bounds
    287                 if (mask->data.U8[r] & PM_READOUT_COMBINE_NO_IMAGE ||
    288                         xIn < 0 || xIn >= readout->image->numCols ||
    289                         yIn < 0 || yIn >= readout->image->numRows) {
     293                psImage *image = readout->image; // The readout image
     294
     295                // Check bounds
     296                if (xIn < 0 || xIn >= image->numCols || yIn < 0 || yIn >= image->numRows) {
    290297                    continue;
    291298                }
     299
     300                pixels->data.F32[r] = image->data.F32[yIn][xIn];
    292301
    293302                // Check mask
    294303                if (readout->mask && readout->mask->data.U8[yIn][xIn] & maskVal) {
    295                     mask->data.U8[r] &= PM_READOUT_COMBINE_MASKED;
     304                    mask->data.U8[r] = 1;
    296305                    continue;
    297306                }
    298307
    299                 pixels->data.F32[r] = readout->image->data.F32[yIn][xIn];
    300 
    301                 // Apply zero and scale
    302                 if (zero) {
    303                     pixels->data.F32[r] -= zero->data.F32[r];
    304                 }
    305                 if (scale) {
    306                     pixels->data.F32[r] /= scale->data.F32[r];
    307                 }
    308 
    309                 if (haveWeights) {
    310                     weights->data.F32[r] = readout->weight->data.F32[yIn][xIn];
    311                     if (scale) {
    312                         weights->data.F32[r] /= scale->data.F32[r] * scale->data.F32[r];
    313                     }
    314                 }
    315308                numValid++;
    316309            }
    … …  
    331324                int numHigh = numValid * params->fracHigh; // Number of high pixels to clip
    332325                // Low pixels
    333                 for (int k = 0, numMasked = 0; numMasked < numLow; k++) {
     326                for (int k = 0, numMasked = 0; numMasked < numLow && k < index->n; k++) {
    334327                    // Don't count the ones that are already masked
    335                     if (! mask->data.U8[index->data.S32[k]] & PM_READOUT_COMBINE_BAD) {
    336                         mask->data.U8[index->data.S32[k]] |= PM_READOUT_COMBINE_CLIPPED;
     328                    if (!mask->data.U8[index->data.S32[k]]) {
     329                        mask->data.U8[index->data.S32[k]] = 1;
    337330                        numMasked++;
    338331                    }
    339332                }
    340333                // High pixels
    341                 for (int k = pixels->n, numMasked = 0; numMasked < numHigh; k++) {
     334                for (int k = pixels->n - 1, numMasked = 0; numMasked < numHigh && k < index->n; k++) {
    342335                    // Don't count the ones that are already masked
    343                     if (! mask->data.U8[index->data.S32[k]] & PM_READOUT_COMBINE_BAD) {
    344                         mask->data.U8[index->data.S32[k]] |= PM_READOUT_COMBINE_CLIPPED;
     336                    if (! mask->data.U8[index->data.S32[k]]) {
     337                        mask->data.U8[index->data.S32[k]] = 1;
    345338                        numMasked++;
    346339                    }
    … …  
    349342
    350343            // Combination
    351             psVectorStats(stats, pixels, weights, mask, PM_READOUT_COMBINE_BAD | PM_READOUT_COMBINE_CLIPPED);
     344            psVectorStats(stats, pixels, weights, mask, 1);
    352345            output->image->data.F32[yOut][xOut] = getStat(stats, params->combine);
    353346            if (haveWeights) {
    354                 output->weight->data.F32[yOut][xOut] = PS_SQR(getStat(stats, combineStdev)); // Variance
    355             }
    356 
    357             // Clear the clipping
    358             psBinaryOp(mask, mask, "&", psScalarAlloc(~PM_READOUT_COMBINE_CLIPPED, PS_TYPE_U8));
    359 
     347                float stdev = getStat(stats, combineStdev);
     348                output->weight->data.F32[yOut][xOut] = PS_SQR(stdev); // Variance
     349            }
    360350        }
    361351    }

  • trunk/psModules/src/imcombine/pmReadoutCombine.h

    r7059 r7256  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2006-05-04 02:38:20 $
     7 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2006-06-01 02:26:32 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4242// Combine multiple readouts, applying zero and scale, with optional minmax clipping
    4343bool pmReadoutCombine(pmReadout *output,// Output readout; altered and returned
    44                       const psArray *inputs, // Array of input readouts (F32 image and weight, U8 mask)
     44                      psArray *inputs, // Array of input readouts (F32 image and weight, U8 mask)
    4545                      const psVector *zero, // Zero corrections to subtract from input, or NULL
    4646                      const psVector *scale, // Scale corrections to divide into input, or NULL
Note: See TracChangeset for help on using the changeset viewer.