IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 20, 2008, 5:08:17 PM (18 years ago)
Author:
Paul Price
Message:

Modernising ppMerge to use the pmFPAfile system of reading/writing, and to use a similar configuration system. The hope is that this will make things easier to maintain.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/pap_branch_080320/ppMerge/src/ppMergeMask.c

    r15937 r17083  
    55#include <stdio.h>
    66#include <string.h>
    7 #include <unistd.h>
    87#include <assert.h>
     8
    99#include <pslib.h>
    1010#include <psmodules.h>
    … …  
    1212
    1313#include "ppMerge.h"
    14 #include "ppMergeData.h"
    15 #include "ppMergeMask.h"
    16 
    17 
    18 // Generate a mask
    19 bool ppMergeMask(ppMergeData *data,  // Data
    20                  ppMergeOptions *options, // Options
    21                  pmConfig *config    // Configuration
    22     )
     14
     15static bool mergeMask(pmConfig *config, // Configuration
     16                      const pmFPAview *view, // View to chip
     17                      bool writeOut,     // Write output?
     18                      psRandom *rng,    // Random number generator
     19                      psMetadata *stats // Statistics output
     20                      )
    2321{
    24     for (int i = 0; i < 2; i++) {
    25         if (!ppMergeMaskSuspect (data, options, config)) {
    26             psError(PS_ERR_UNKNOWN, true, "failed on mask suspect.\n");
    27             return false;
    28         }
    29 
    30         if (!ppMergeMaskBad (data, options, config)) {
    31             psError(PS_ERR_UNKNOWN, true, "failed on mask bad.\n");
    32             return false;
    33         }
    34     }
    35 
    36     if (!ppMergeMaskAverageConcepts (data, options, config)) {
    37         psError(PS_ERR_UNKNOWN, true, "failed on average concepts.\n");
    38         return false;
    39     }
    40 
    41     if (!ppMergeMaskGrow (data, options, config)) {
    42         psError(PS_ERR_UNKNOWN, true, "failed on mask grow.\n");
    43         return false;
    44     }
    45 
    46     if (!ppMergeMaskWrite (data, options, config)) {
    47         psError(PS_ERR_UNKNOWN, true, "failed on mask write.\n");
    48         return false;
    49     }
    50 
     22    assert(config);
     23    assert(view);
     24    assert(view->chip != -1 && view->cell == -1 && view->readout == -1);
     25
     26    bool mdok;                          // Status of MD lookup
     27    int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of input files
     28    psStatsOptions meanStat = psMetadataLookupS32(NULL, config->arguments, "MEAN"); // Statistic for mean
     29    psStatsOptions stdevStat = psMetadataLookupS32(NULL, config->arguments, "STDEV"); // Statistic for stdev
     30    psMaskType maskVal = psMetadataLookupU8(NULL, config->arguments, "MASKVAL"); // Value to mask
     31    int sample = psMetadataLookupS32(NULL, config->arguments, "SAMPLE"); // Size of sample for statistics
     32    bool chipStats = psMetadataLookupBool(&mdok, config->arguments, "MASK.CHIPSTATS"); // Statistics on chip?
     33    float maskSuspect = psMetadataLookupF32(NULL, config->arguments, "MASK.SUSPECT"); // Threshold for suspect pixels
     34    float maskBad = psMetadataLookupF32(NULL, config->arguments, "MASK.BAD"); // Threshold for bad pixels
     35    pmMaskIdentifyMode maskMode = psMetadataLookupS32(NULL, config->arguments, "MASK.MODE"); // Mode for identifying bad pixels
     36    int maskGrow = psMetadataLookupS32(NULL, config->arguments, "MASK.GROW"); // Radius to grow mask
     37    psMaskType maskGrowVal = psMetadataLookupU8(NULL, config->arguments, "MASK.GROWVAL"); // Value for grown mask
     38
     39    psStats *statistics = psStatsAlloc(meanStat | stdevStat); // Statistics for background
     40
     41    psString outName = ppMergeOutputFile(config); // Name of output file
     42    pmChip *outChip = pmFPAfileThisChip(config->files, view, outName); // Output chip
     43    psFree(outName);
     44
     45    // For each input file, get the statistics, which can be calculated at the chip or cell levels
     46    psArray *suspects = psArrayAlloc(outChip->cells->n); // Images with suspected pixels
     47    psVector *values = psVectorAlloc(sample, PS_TYPE_F32); // Pixel values for statistics
     48    pmFPAview *inView = pmFPAviewAlloc(0); // View for input
     49    for (int i = 0; i < numFiles; i++) {
     50        pmFPAfileActivate(config->files, false, NULL);
     51        psArray *files = ppMergeFileActivateSingle(config, PPMERGE_FILES_INPUT, true, i); // Input files
     52        pmFPAfile *input = files->data[0]; // Input file
     53        psFree(files);
     54        pmFPA *inFPA = input->fpa;  // Input FPA
     55        *inView = *view;
     56
     57        int valueIndex = 0;             // Index for vector of pixel values
     58
     59        pmCell *inCell;                 // Input cell
     60        while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) {
     61            pmHDU *hdu = pmHDUFromCell(inCell); // HDU for cell
     62            if (!hdu || hdu->blankPHU) {
     63                // No data here
     64                continue;
     65            }
     66            psTrace("ppMerge", 1, "Getting suspect pixels for file %d chip %d cell %d",
     67                    i, inView->chip, inView->cell);
     68
     69            if (!pmFPAfileIOChecks(config, inView, PM_FPA_BEFORE)) {
     70                psFree(inView);
     71                goto MERGE_MASK_ERROR;
     72            }
     73
     74            if (!ppMergeFileOpenInput(config, view, i)) {
     75                psError(PS_ERR_UNKNOWN, false, "Unable to open file %d", i);
     76                psFree(inView);
     77                goto MERGE_MASK_ERROR;
     78            }
     79
     80            if (inCell->readouts->n > 1) {
     81                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     82                        "File %d chip %d cell %d contains more than one readout (%ld)",
     83                        i, inView->chip, inView->cell, inCell->readouts->n);
     84                psFree(inView);
     85                goto MERGE_MASK_ERROR;
     86            }
     87
     88            pmReadout *readout;
     89            if (inCell->readouts && inCell->readouts->n == 1) {
     90                readout = psMemIncrRefCounter(inCell->readouts->data[0]); // Input readout
     91            } else {
     92                readout = pmReadoutAlloc(inCell);
     93            }
     94
     95            if (!ppMergeFileReadInput(config, readout, i, 0)) {
     96                psError(PS_ERR_UNKNOWN, false, "Unable to read readout %d", i);
     97                psFree(inView);
     98                psFree(readout);
     99                goto MERGE_MASK_ERROR;
     100            }
     101
     102            psImage *outMask = NULL;    // Output mask image (for iterative generation of mask)
     103            pmCell *outCell = pmFPAfileThisCell(config->files, inView, "PPMERGE.OUTPUT.MASK"); // Output cell
     104            if (outCell->readouts && outCell->readouts->n == 1) {
     105                pmReadout *outRO = outCell->readouts->data[0]; // Output readout
     106                if (outRO) {
     107                    outMask = outRO->mask;
     108                }
     109            }
     110
     111            psImage *image = readout->image, *mask = readout->mask; // Image and mask
     112            int numCols = readout->image->numCols, numRows = readout->image->numRows; // Image size
     113            int numPix = numCols * numRows; // Number of pixels
     114            int numCells = chipStats ? readout->parent->parent->cells->n : 1; // Number of cells
     115            int num = PS_MIN(numPix, sample / numCells); // Number of values to add
     116            if (!chipStats) {
     117                valueIndex = 0;
     118            }
     119            for (int i = 0; i < num; i++) {
     120                int pixel = numPix * psRandomUniform(rng);
     121                int x = pixel % numCols;
     122                int y = pixel / numCols;
     123                if ((mask && (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal)) ||
     124                    !isfinite(image->data.F32[y][x]) ||
     125                    (outMask && (outMask->data.PS_TYPE_MASK_DATA[y][x] & maskVal))) {
     126                    continue;
     127                }
     128
     129                values->data.F32[valueIndex++] = image->data.F32[y][x];
     130            }
     131
     132            if (!chipStats) {
     133                values->n = valueIndex;
     134                if (!psVectorStats(statistics, values, NULL, NULL, 0)) {
     135                    psError(PS_ERR_UNKNOWN, false, "Unable to do statistics on readout.");
     136                    psFree(inView);
     137                    psFree(readout);
     138                    goto MERGE_MASK_ERROR;
     139                }
     140
     141                psMetadataAddF32(readout->analysis, PS_LIST_TAIL, PM_MASK_ANALYSIS_MEAN, PS_META_REPLACE,
     142                                 "Mean value of readout", psStatsGetValue(statistics, meanStat));
     143                psMetadataAddF32(readout->analysis, PS_LIST_TAIL, PM_MASK_ANALYSIS_STDEV, PS_META_REPLACE,
     144                                 "Stdev value of readout", psStatsGetValue(statistics, stdevStat));
     145                suspects->data[inView->cell] = pmMaskFlagSuspectPixels(suspects->data[inView->cell],
     146                                                                       readout, maskSuspect, maskVal);
     147                pmCellFreeData(inCell);
     148
     149                if (!pmFPAfileIOChecks(config, inView, PM_FPA_AFTER)) {
     150                    psFree(inView);
     151                    psFree(readout);
     152                    goto MERGE_MASK_ERROR;
     153                }
     154            }
     155            psFree(readout);
     156        }
     157
     158        // Additional run through cells if we want chip-level statistics
     159        if (chipStats && valueIndex > 0) {
     160            values->n = valueIndex;
     161            if (!psVectorStats(statistics, values, NULL, NULL, 0)) {
     162                psError(PS_ERR_UNKNOWN, false, "Unable to do statistics on chip.");
     163                goto MERGE_MASK_ERROR;
     164            }
     165            inView->cell = -1;
     166            while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) {
     167                pmHDU *hdu = pmHDUFromCell(inCell); // HDU for cell
     168                if (!hdu || hdu->blankPHU) {
     169                    // No data here
     170                    continue;
     171                }
     172                pmReadout *readout = inCell->readouts->data[0]; // Readout of interest
     173
     174                psMetadataAddF32(readout->analysis, PS_LIST_TAIL, PM_MASK_ANALYSIS_MEAN, PS_META_REPLACE,
     175                                 "Mean value of chip", psStatsGetValue(statistics, meanStat));
     176                psMetadataAddF32(readout->analysis, PS_LIST_TAIL, PM_MASK_ANALYSIS_STDEV, PS_META_REPLACE,
     177                                 "Stdev value of chip", psStatsGetValue(statistics, stdevStat));
     178                suspects->data[inView->cell] = pmMaskFlagSuspectPixels(suspects->data[inView->cell],
     179                                                                       readout, maskSuspect, maskVal);
     180                pmCellFreeData(inCell);
     181
     182                if (!pmFPAfileIOChecks(config, inView, PM_FPA_AFTER)) {
     183                    psFree(inView);
     184                    goto MERGE_MASK_ERROR;
     185                }
     186            }
     187        }
     188    }
     189    psFree(inView);
     190    psFree(statistics); statistics = NULL;
     191    psFree(values); values = NULL;
     192
     193
     194    // Another run through the chip to threshold on the suspects
     195    pmFPAfileActivate(config->files, false, NULL);
     196    if (writeOut) {
     197        ppMergeFileActivate(config, PPMERGE_FILES_OUTPUT, true);
     198    }
     199    pmFPA *outFPA = outChip->parent;    // Output FPA
     200    pmCell *outCell;                    // Output cell
     201    pmFPAview *outView = pmFPAviewAlloc(0); // View into output FPA
     202    *outView = *view;
     203    while ((outCell = pmFPAviewNextCell(outView, outFPA, 1))) {
     204        pmHDU *hdu = pmHDUFromCell(outCell); // HDU for cell
     205        if (!hdu || hdu->blankPHU) {
     206            // No data here
     207            continue;
     208        }
     209
     210        psTrace("ppMerge", 1, "Getting bad pixels for chip %d cell %d", outView->chip, outView->cell);
     211
     212        pmReadout *outRO = NULL;    // Output readout
     213        if (outCell->readouts && outCell->readouts->n == 1) {
     214            outRO = psMemIncrRefCounter(outCell->readouts->data[0]);
     215        } else {
     216            outRO = pmReadoutAlloc(outCell);
     217        }
     218
     219        psImage *mask = pmMaskIdentifyBadPixels(suspects->data[outView->cell], maskVal, numFiles, maskBad,
     220                                                maskMode);
     221        psFree(suspects->data[outView->cell]);
     222        suspects->data[outView->cell] = NULL;
     223        if (maskGrowVal > 0) {
     224            outRO->mask = psImageGrowMask(outRO->mask, mask, maskVal, maskGrow, maskGrowVal);
     225            psFree(mask);
     226        } else {
     227            psFree(outRO->mask);
     228            outRO->mask = mask;
     229        }
     230
     231        outRO->data_exists = outCell->data_exists = outChip->data_exists = true;
     232
     233        if (writeOut) {
     234            if (!pmFPAfileIOChecks(config, outView, PM_FPA_BEFORE)) {
     235                psFree(outView);
     236                goto MERGE_MASK_ERROR;
     237            }
     238
     239            // Average concepts
     240            psList *cells = psListAlloc(NULL); // List of cells, for concept averaging
     241            for (int i = 0; i < numFiles; i++) {
     242                pmFPAfile *inFile = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file
     243                pmCell *inCell = pmFPAviewThisCell(outView, inFile->fpa); // Input cell
     244                psListAdd(cells, PS_LIST_TAIL, inCell);
     245            }
     246            if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, true)) {
     247                psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts.");
     248                psFree(cells);
     249                psFree(outRO);
     250                psFree(outView);
     251                return false;
     252            }
     253            psFree(cells);
     254
     255            // Statistics on the merged cell using a fake image
     256            outRO->image = psImageAlloc(outRO->mask->numCols, outRO->mask->numRows, PS_TYPE_F32);
     257            psImageInit(outRO->image, 1.0);
     258            if (!ppStatsFPA(stats, outRO->parent->parent->parent, outView,
     259                            maskVal | pmConfigMask("BLANK", config), config)) {
     260                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to generate stats for image.");
     261                psFree(outRO);
     262                psFree(outView);
     263                return false;
     264            }
     265            psFree(outRO->image);
     266            outRO->image = NULL;
     267
     268            // Write
     269            if (!pmFPAfileIOChecks(config, outView, PM_FPA_AFTER)) {
     270                psFree(outView);
     271                goto MERGE_MASK_ERROR;
     272            }
     273        }
     274        psFree(outRO);
     275    }
     276    psFree(outView);
     277
     278    ppMergeFileActivate(config, PPMERGE_FILES_ALL, true);
     279
     280    psFree(suspects);
    51281    return true;
     282
     283
     284MERGE_MASK_ERROR:
     285    psFree(suspects);
     286    psFree(statistics);
     287    psFree(values);
     288    return false;
    52289}
    53290
    54 bool ppMergeMaskReadoutStats(const pmReadout *readout,
    55                              const pmReadout *output,
    56                              ppMergeOptions *options, // Options
    57                              psRandom *rng)
     291
     292
     293
     294
     295bool ppMergeMask(pmConfig *config)
    58296{
    59     PS_ASSERT_PTR_NON_NULL(readout, false);
    60     PS_ASSERT_IMAGE_NON_NULL(readout->image, false);
    61     PS_ASSERT_IMAGE_NON_EMPTY(readout->image, false);
    62     PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false);
    63     if (readout->mask) {
    64         PS_ASSERT_IMAGE_NON_EMPTY(readout->mask, false);
    65         PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, readout->mask, false);
    66         PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_MASK, false);
    67     }
    68 
    69     psImage *mask = NULL;
    70     psImage *image = readout->image;    // Image of interest
    71 
    72     if (output) {
    73         mask = output->mask;      // Corresponding mask
    74     }
    75 
    76     if (rng) {
    77         psMemIncrRefCounter(rng);
    78     } else {
    79         rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
    80     }
    81 
    82     // XXX note that this now will accept any of several stats options
    83     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    84     stats->nSubsample = options->sample;
    85 
    86     if (!psImageBackground(stats, NULL, image, mask, options->combine->maskVal, rng)) {
    87         psError(PS_ERR_UNKNOWN, false, "Failure to measure image statistics.\n");
    88         psFree(stats);
    89         psFree(rng);
    90         return false;
    91     }
    92     if (!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) {
    93         psError(PS_ERR_UNKNOWN, false, "invalide image statistics (nan).\n");
    94         psFree(stats);
    95         psFree(rng);
    96         return false;
    97     }
    98 
    99     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.MEDIAN", PS_META_REPLACE, "image stats", stats->robustMedian);
    100     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.STDEV",  PS_META_REPLACE, "image stats", stats->robustStdev);
    101 
     297    assert(config);
     298
     299    bool mdok;                          // Status of MD lookup
     300    int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
     301    bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks?
     302    bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); // Do we have weights?
     303    int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Number of rejection iterations
     304
     305    PS_ASSERT_INT_POSITIVE(iter, false);
     306
     307    pmFPAview *view = pmFPAviewAlloc(0); // View to component of interest
     308    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
     309
     310    psMetadata *stats = NULL;           // Statistics for output
     311    if (psMetadataLookup(config->arguments, "STATS.NAME")) {
     312        stats = psMetadataAlloc();
     313        psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "STATS.DATA", 0, "Statistics output", stats);
     314    }
     315
     316    psArray *inputs = ppMergeFileDataLevel(config, "PPMERGE.INPUT", PM_FPA_LEVEL_READOUT); // Input images
     317    psFree(inputs);
     318    if (haveMasks) {
     319        psArray *masks = ppMergeFileDataLevel(config, "PPMERGE.INPUT.MASK", PM_FPA_LEVEL_READOUT);
     320        psFree(masks);
     321    }
     322    if (haveWeights) {
     323        psArray *weights = ppMergeFileDataLevel(config, "PPMERGE.INPUT.WEIGHT", PM_FPA_LEVEL_READOUT);
     324        psFree(weights);
     325    }
     326
     327    if (!ppMergeFileActivate(config, PPMERGE_FILES_INPUT, true)) {
     328        psError(PS_ERR_UNKNOWN, false, "Unable to activate files.");
     329        goto PPMERGE_MASK_ERROR;
     330    }
     331    if (!ppMergeFileActivate(config, PPMERGE_FILES_OUTPUT, false)) {
     332        psError(PS_ERR_UNKNOWN, false, "Unable to activate files.");
     333        goto PPMERGE_MASK_ERROR;
     334    }
     335
     336    psString outName = ppMergeOutputFile(config); // Name of output file
     337    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, outName); // Output file
     338    psFree(outName);
     339    assert(output && output->fpa);
     340    pmFPA *outFPA = output->fpa;        // Output FPA
     341
     342    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     343        goto PPMERGE_MASK_ERROR;
     344    }
     345    pmChip *outChip;                    // Chip of interest
     346    while ((outChip = pmFPAviewNextChip(view, outFPA, 1))) {
     347        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     348            goto PPMERGE_MASK_ERROR;
     349        }
     350
     351        for (int i = 0; i < iter; i++) {
     352            if (!mergeMask(config, view, (i == iter - 1), rng, stats)) {
     353                psError(PS_ERR_UNKNOWN, false, "Unable to merge chip %d", view->chip);
     354                goto PPMERGE_MASK_ERROR;
     355            }
     356        }
     357
     358        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     359            goto PPMERGE_MASK_ERROR;
     360        }
     361    }
     362
     363    psList *fpaList = psListAlloc(NULL);// List of FPAs for concept averaging
     364    for (int i = 0; i < numFiles; i++) {
     365        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file
     366        psListAdd(fpaList, PS_LIST_TAIL, file->fpa);
     367    }
     368    if (!pmConceptsAverageFPAs(outFPA, fpaList)) {
     369        psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts.");
     370        psFree(fpaList);
     371        goto PPMERGE_MASK_ERROR;
     372    }
     373    psFree(fpaList);
     374
     375    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     376        goto PPMERGE_MASK_ERROR;
     377    }
     378
     379    psFree(stats);
     380    psFree(view);
    102381    psFree(rng);
     382
    103383    return true;
     384
     385PPMERGE_MASK_ERROR:
     386    psFree(stats);
     387    psFree(view);
     388    psFree(rng);
     389    return false;
    104390}
    105391
    106 bool ppMergeMaskChipStats (const pmChip *chip,
    107                            const pmChip *output,
    108                            ppMergeOptions *options,
    109                            psRandom *rng)
    110 {
    111     PS_ASSERT_PTR_NON_NULL(chip, false);
    112 
    113     // XXX note that this now will accept any of several stats options
    114     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    115 
    116     if (rng) {
    117         psMemIncrRefCounter(rng);
    118     } else {
    119         rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
    120     }
    121 
    122     // accumulate a vector of data values using options->sample per readout
    123     psVector *values = psVectorAllocEmpty(options->sample, PS_TYPE_F32); // Vector containing subsample
    124 
    125     for (int nCell = 0; nCell < chip->cells->n; nCell++) {
    126         pmCell *cell = chip->cells->data[nCell];
    127         if (cell->readouts->n == 0) continue;
    128 
    129         for (int nReadout = 0; nReadout < cell->readouts->n; nReadout++) {
    130             pmReadout *readout = cell->readouts->data[nReadout];
    131             if (!readout->image) continue;
    132 
    133             psTrace("ppMerge", 4, "Measure statistics for cell %d, readout %d\n", nCell, nReadout);
    134 
    135             psImage *image = readout->image;    // Image of interest
    136 
    137             pmCell *cellOutput = output->cells->data[nCell];
    138             pmReadout *readoutOutput = NULL;
    139             psImage *mask = NULL;
    140             if (cellOutput->readouts->n > 0) {
    141                 readoutOutput = cellOutput->readouts->data[0];
    142                 mask = readoutOutput->mask;      // Corresponding mask
    143             }
    144 
    145             // Size of image
    146             long nx = image->numCols;
    147             long ny = image->numRows;
    148             const int Npixels = nx*ny;  // Total number of pixels
    149             const int Nsubset = PS_MIN(options->sample, Npixels); // Number of pixels in subset
    150            
    151             values = psVectorRealloc (values, values->n + Nsubset); // make sure we have enough space
    152 
    153             // select a subset of the image pixels to measure the stats
    154             for (long i = 0; i < Nsubset; i++) {
    155                 double frnd = psRandomUniform(rng);
    156                 int pixel = Npixels * frnd;
    157                 int ix = pixel % nx;
    158                 int iy = pixel / nx;
    159 
    160                 if (!isfinite(image->data.F32[iy][ix])) continue;
    161                 if (mask && (mask->data.U8[iy][ix] & options->combine->maskVal)) continue;
    162 
    163                 float value = image->data.F32[iy][ix];
    164                 values->data.F32[values->n] = value;
    165                 values->n ++;
    166             }
    167         }
    168     }
    169 
    170     // no valid data, skip the chip
    171     if (!values->n) {
    172         psFree(values);
    173         psFree(stats);
    174         psFree(rng);
    175         return true;
    176     }
    177 
    178     // calculate the statistics
    179     if (!psVectorStats (stats, values, NULL, NULL, 0)) {
    180         psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for chip");
    181         psFree(values);
    182         psFree(stats);
    183         psFree(rng);
    184         return false;
    185     }
    186     if (!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) {
    187         psError(PS_ERR_UNKNOWN, false, "invalid image statistics (nan).\n");
    188         psFree(values);
    189         psFree(stats);
    190         psFree(rng);
    191         return false;
    192     }
    193 
    194     // supply the stats to the readout analysis metadata
    195     for (int nCell = 0; nCell < chip->cells->n; nCell++) {
    196         pmCell *cell = chip->cells->data[nCell];
    197         if (cell->readouts->n == 0) continue;
    198 
    199         for (int nReadout = 0; nReadout < cell->readouts->n; nReadout++) {
    200             pmReadout *readout = cell->readouts->data[nReadout];
    201             if (!readout->image) continue;
    202 
    203             psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.MEDIAN", PS_META_REPLACE, "image stats", stats->robustMedian);
    204             psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.STDEV",  PS_META_REPLACE, "image stats", stats->robustStdev);
    205         }
    206     }
    207 
    208     psLogMsg ("ppMerge", PS_LOG_INFO, "statistics for chip: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
    209 
    210     psFree(values);
    211     psFree(stats);
    212     psFree(rng);
    213     return true;
    214 }
Note: See TracChangeset for help on using the changeset viewer.