IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 18811


Ignore:
Timestamp:
Jul 30, 2008, 6:01:40 PM (18 years ago)
Author:
eugene
Message:

thread-friendly implementation of pmReadoutCombine (init with pmReadoutCombinePrepare)

Location:
branches/eam_branch_20080719/psModules/src/imcombine
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branch_20080719/psModules/src/imcombine/pmReadoutCombine.c

    r17228 r18811  
    4242}
    4343
     44// check the input parameters and set up the output images
     45bool pmReadoutCombinePrepare(pmReadout *output, const psArray *inputs, const pmCombineParams *params)
     46{
     47    // Check inputs
     48    PS_ASSERT_PTR_NON_NULL(output, false);
     49    PS_ASSERT_ARRAY_NON_NULL(inputs, false);
     50    PS_ASSERT_PTR_NON_NULL(params, false);
     51    PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracLow, 0.0, 1.0, false);
     52    PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracHigh, 0.0, 1.0, false);
     53
     54    // valid combintion statistic?
     55    bool valid = false;
     56    valid |= (params->combine == PS_STAT_SAMPLE_MEAN);
     57    valid |= (params->combine == PS_STAT_SAMPLE_MEDIAN);
     58    valid |= (params->combine == PS_STAT_ROBUST_MEDIAN);
     59    valid |= (params->combine == PS_STAT_FITTED_MEAN);
     60    valid |= (params->combine == PS_STAT_CLIPPED_MEAN);
     61    if (!valid) {
     62        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Combination method is not SAMPLE_MEAN, SAMPLE_MEDIAN, "
     63                "ROBUST_MEDIAN, FITTED_MEAN or CLIPPED_MEAN.\n");
     64        return false;
     65    }
     66
     67    // weights exist if weights desired?
     68    for (int i = 0; i < inputs->n; i++) {
     69        pmReadout *readout = inputs->data[i]; // Readout of interest
     70        if (params->weights && !readout->weight) {
     71            psError(PS_ERR_UNEXPECTED_NULL, true,
     72                    "Rejection based on weights requested, but no weights supplied for image %d.\n", i);
     73            return false;
     74        }
     75    }
     76
     77    pmHDU *hdu = pmHDUFromReadout(output); // Output HDU
     78    if (!hdu) {
     79        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find HDU for readout.\n");
     80        return false;
     81    }
     82
     83    //  set the output header metadata
     84    psString comment = NULL;        // Comment to add to header
     85    psStringAppend(&comment, "Combining using statistic: %x", params->combine);
     86    if (!hdu->header) {
     87        hdu->header = psMetadataAlloc();
     88    }
     89    psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
     90    psFree(comment);
     91
     92    // note the clipping parameters, if used
     93    if (params->combine == PS_STAT_CLIPPED_MEAN) {
     94        psString comment = NULL;    // Comment to add to header
     95        psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", params->iter, params->rej);
     96        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
     97        psFree(comment);
     98    }
     99
     100    // note the use of weights
     101    if (params->weights) {
     102        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
     103                         "Using input weights to combine images", "");
     104    }
     105
     106    // note the rejection fraction
     107    float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of pixels to keep
     108    if (keepFrac != 1.0) {
     109        psString comment = NULL;        // Comment to add to header
     110        psStringAppend(&comment, "Min/max rejection: %f high, %f low, keep %d",
     111                       params->fracHigh, params->fracLow, params->nKeep);
     112        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
     113        psFree(comment);
     114    }
     115
     116    // note the mask value actually used
     117    psMaskType maskVal = params->maskVal; // The mask value
     118    if (maskVal) {
     119        psString comment = NULL;        // Comment to add to header
     120        psStringAppend(&comment, "Mask for combination: %x", maskVal);
     121        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
     122        psFree(comment);
     123    }
     124
     125    // determine the output image size based on the input images
     126    int row0, col0, numCols, numRows;
     127    if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, inputs)) {
     128        psError(PS_ERR_UNKNOWN, false, "problem setting output readout size.");
     129        return false;
     130    }
     131
     132    // generate the required output images based on the specified sizes
     133    pmReadoutStackDefineOutput(output, col0, row0, numCols, numRows, true, params->weights, params->blank);
     134    psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0);
     135
     136    // these calls allocate and save the requested images on the output analysis metadata
     137    psImage *counts = pmReadoutSetAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_COUNT, numCols, numRows, PS_TYPE_U16, 0);
     138    if (!counts) {
     139        return false;
     140    }
     141    psImage *sigma = pmReadoutSetAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_SIGMA, numCols, numRows, PS_TYPE_F32, NAN);
     142    if (!sigma) {
     143        return false;
     144    }
     145
     146    // Update the "concepts"
     147    psList *inputCells = psListAlloc(NULL); // List of cells
     148    for (long i = 0; i < inputs->n; i++) {
     149        pmReadout *readout = inputs->data[i]; // Readout of interest
     150        psListAdd(inputCells, PS_LIST_TAIL, readout->parent);
     151    }
     152    bool success = pmConceptsAverageCells(output->parent, inputCells, NULL, NULL, true);
     153    psFree(inputCells);
     154
     155    // set these even though the values are not yet set
     156    output->data_exists = true;
     157    output->parent->data_exists = true;
     158    output->parent->parent->data_exists = true;
     159
     160    return success;
     161}
    44162
    45163// XXX: Maybe add support for S16 and S32 types.  Currently, only F32 supported.
    … …  
    61179    PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracLow, 0.0, 1.0, false);
    62180    PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracHigh, 0.0, 1.0, false);
    63     if (params->combine != PS_STAT_SAMPLE_MEAN && params->combine != PS_STAT_SAMPLE_MEDIAN &&
    64             params->combine != PS_STAT_ROBUST_MEDIAN && params->combine != PS_STAT_FITTED_MEAN &&
    65             params->combine != PS_STAT_CLIPPED_MEAN) {
    66         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Combination method is not SAMPLE_MEAN, SAMPLE_MEDIAN, "
    67                 "ROBUST_MEDIAN, FITTED_MEAN or CLIPPED_MEAN.\n");
    68         return false;
    69     }
    70     for (int i = 0; i < inputs->n; i++) {
    71         pmReadout *readout = inputs->data[i]; // Readout of interest
    72         if (params->weights && !readout->weight) {
    73             psError(PS_ERR_UNEXPECTED_NULL, true,
    74                     "Rejection based on weights requested, but no weights supplied for image %d.\n", i);
    75             return false;
    76         }
    77     }
    78 
    79     bool first = !output->image;        // First pass through?
    80181
    81182    pmHDU *hdu = pmHDUFromReadout(output); // Output HDU
    … …  
    85186    }
    86187
    87     if (first) {
    88         psString comment = NULL;        // Comment to add to header
    89         psStringAppend(&comment, "Combining using statistic: %x", params->combine);
    90         if (!hdu->header) {
    91             hdu->header = psMetadataAlloc();
    92         }
    93         psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
    94         psFree(comment);
    95     }
     188    pthread_t id = pthread_self();
     189    char name[64];
     190    sprintf (name, "%x", (unsigned int) id);
     191    psTimerStart (name);
    96192
    97193    psStatsOptions combineStdev = 0; // Statistics option for weights
    … …  
    118214        stats->clipSigma = params->rej;
    119215        stats->clipIter = params->iter;
    120 
    121         if (first) {
    122             psString comment = NULL;    // Comment to add to header
    123             psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma",
    124                            params->iter, params->rej);
    125             psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
    126             psFree(comment);
    127         }
    128     }
    129     if (params->weights && first) {
    130         psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
    131                          "Using input weights to combine images", "");
    132     }
     216    }
     217
     218    psImage *counts = pmReadoutGetAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_COUNT);
     219    if (!counts) {
     220        return false;
     221    }
     222    psImage *sigma = pmReadoutGetAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_SIGMA);
     223    if (!sigma) {
     224        return false;
     225    }
     226
     227    stats->options |= combineStdev;
    133228
    134229    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
    135230    int xSize, ySize;                   // Size of the output image
    136     if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
    137                                 inputs)) {
     231    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, inputs)) {
    138232        psError(PS_ERR_UNKNOWN, false, "No valid input readouts.");
    139233        return false;
    140234    }
    141 
    142     pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true, params->weights,
    143                         params->blank);
    144     psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0);
    145 
    146     psImage *counts = pmReadoutAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_COUNT, xSize, ySize,
    147                                              PS_TYPE_U16, 0);
    148     if (!counts) {
    149         return false;
    150     }
    151     psImage *sigma = pmReadoutAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_SIGMA, xSize, ySize,
    152                                             PS_TYPE_F32, NAN);
    153     if (!sigma) {
    154         psFree(counts);
    155         return false;
    156     }
    157 
    158     stats->options |= combineStdev;
    159235
    160236    // We loop through each pixel in the output image.  We loop through each input readout.  We determine if
    … …  
    180256
    181257    float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of pixels to keep
    182     if (keepFrac != 1.0 && first) {
    183         psString comment = NULL;        // Comment to add to header
    184         psStringAppend(&comment, "Min/max rejection: %f high, %f low, keep %d",
    185                        params->fracHigh, params->fracLow, params->nKeep);
    186         psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
    187         psFree(comment);
    188     }
    189 
    190258    psMaskType maskVal = params->maskVal; // The mask value
    191     if (maskVal && first) {
    192         psString comment = NULL;        // Comment to add to header
    193         psStringAppend(&comment, "Mask for combination: %x", maskVal);
    194         psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
    195         psFree(comment);
    196     }
    197259
    198260    #ifndef PS_NO_TRACE
    … …  
    226288        int yOut = i - output->row0; // y position on output readout
    227289        #ifdef SHOW_BUSY
    228 
     290       
    229291        if (psTraceGetLevel("psModules.imcombine") > 9) {
    230292            printf("Processing row %d\r", i);
    … …  
    242304                int xIn = j - readout->col0; // x position on input readout
    243305                psImage *image = readout->image; // The readout image
    244 
    245                 #if 0 // This should have been taken care of already:
    246                 // Check bounds
    247                 if (xIn < 0 || xIn >= image->numCols || yIn < 0 || yIn >= image->numRows) {
    248                     continue;
    249                 }
    250                 #endif
    251306
    252307                pixelsData[r] = image->data.F32[yIn][xIn];
    … …  
    348403    }
    349404    #endif
     405
    350406    psFree(index);
    351407    psFree(pixels);
    … …  
    355411    psFree(stats);
    356412    psFree(invScale);
    357     psFree(counts);
    358     psFree(sigma);
    359 
    360     // Update the "concepts"
    361     psList *inputCells = psListAlloc(NULL); // List of cells
    362     for (long i = 0; i < inputs->n; i++) {
    363         pmReadout *readout = inputs->data[i]; // Readout of interest
    364         psListAdd(inputCells, PS_LIST_TAIL, readout->parent);
    365     }
    366     bool success = pmConceptsAverageCells(output->parent, inputCells, NULL, NULL, true);
    367     psFree(inputCells);
    368 
    369     output->data_exists = true;
    370     output->parent->data_exists = true;
    371     output->parent->parent->data_exists = true;
    372 
    373     return success;
     413
     414    // fprintf (stderr, "done with combine %x : %f sec\n", (unsigned int) id, psTimerMark (name));
     415    return true;
    374416}
    375417

  • branches/eam_branch_20080719/psModules/src/imcombine/pmReadoutCombine.h

    r17228 r18811  
    55 * @author Paul Price, IfA
    66 *
    7  * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
    8  * @date $Date: 2008-03-29 03:10:17 $
     7 * @version $Revision: 1.13.20.1 $ $Name: not supported by cvs2svn $
     8 * @date $Date: 2008-07-31 04:01:40 $
    99 * Copyright 2004-2006 Institute for Astronomy, University of Hawaii
    1010 */
    … …  
    3636                                     );
    3737
     38// check the input parameters and set up the output images
     39bool pmReadoutCombinePrepare(pmReadout *output, const psArray *inputs, const pmCombineParams *params);
     40
    3841/// Combine multiple readouts, applying zero and scale, with optional minmax clipping
    3942bool pmReadoutCombine(pmReadout *output,///< Output readout; altered and returned
Note: See TracChangeset for help on using the changeset viewer.