IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

Merging in pap_branch_080320 --- modernised version of ppMerge.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppMerge/src/ppMergeMask.c

    r15937 r17227  
    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     }
     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    psVector *values = psVectorAlloc(sample, PS_TYPE_F32); // Pixel values for statistics
     47    pmFPAview *inView = pmFPAviewAlloc(0); // View for input
     48    for (int i = 0; i < numFiles; i++) {
     49        pmFPAfileActivate(config->files, false, NULL);
     50        psArray *files = ppMergeFileActivateSingle(config, PPMERGE_FILES_INPUT, true, i); // Input files
     51        pmFPAfile *input = files->data[0]; // Input file
     52        psFree(files);
     53        pmFPA *inFPA = input->fpa;  // Input FPA
     54        *inView = *view;
     55
     56        int valueIndex = 0;             // Index for vector of pixel values
     57
     58        pmCell *inCell;                 // Input cell
     59        while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) {
     60            pmHDU *hdu = pmHDUFromCell(inCell); // HDU for cell
     61            if (!hdu || hdu->blankPHU) {
     62                // No data here
     63                continue;
     64            }
     65            psTrace("ppMerge", 1, "Getting suspect pixels for file %d chip %d cell %d",
     66                    i, inView->chip, inView->cell);
     67
     68            if (!pmFPAfileIOChecks(config, inView, PM_FPA_BEFORE)) {
     69                psFree(inView);
     70                goto MERGE_MASK_ERROR;
     71            }
     72
     73            if (!ppMergeFileOpenInput(config, view, i)) {
     74                psError(PS_ERR_UNKNOWN, false, "Unable to open file %d", i);
     75                psFree(inView);
     76                goto MERGE_MASK_ERROR;
     77            }
     78
     79            if (inCell->readouts->n > 1) {
     80                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     81                        "File %d chip %d cell %d contains more than one readout (%ld)",
     82                        i, inView->chip, inView->cell, inCell->readouts->n);
     83                psFree(inView);
     84                goto MERGE_MASK_ERROR;
     85            }
     86
     87            pmReadout *readout;
     88            if (inCell->readouts && inCell->readouts->n == 1) {
     89                readout = psMemIncrRefCounter(inCell->readouts->data[0]); // Input readout
     90            } else {
     91                readout = pmReadoutAlloc(inCell);
     92            }
     93
     94            if (!ppMergeFileReadInput(config, readout, i, 0)) {
     95                psError(PS_ERR_UNKNOWN, false, "Unable to read readout %d", i);
     96                psFree(inView);
     97                psFree(readout);
     98                goto MERGE_MASK_ERROR;
     99            }
     100
     101            pmCell *outCell = pmFPAfileThisCell(config->files, inView, "PPMERGE.OUTPUT.MASK"); // Output cell
     102            pmReadout *outRO = NULL;    // Output readout
     103            if (outCell->readouts && outCell->readouts->n == 1) {
     104                outRO = psMemIncrRefCounter(outCell->readouts->data[0]);
     105            } else {
     106                outRO = pmReadoutAlloc(outCell);
     107            }
     108            psImage *outMask = outRO->mask;    // Output mask image (for iterative generation of mask)
     109
     110            psImage *image = readout->image, *mask = readout->mask; // Image and mask
     111            int numCols = readout->image->numCols, numRows = readout->image->numRows; // Image size
     112            int numPix = numCols * numRows; // Number of pixels
     113            int numCells = chipStats ? readout->parent->parent->cells->n : 1; // Number of cells
     114            int num = PS_MIN(numPix, sample / numCells); // Number of values to add
     115            if (!chipStats) {
     116                valueIndex = 0;
     117            }
     118            for (int i = 0; i < num; i++) {
     119                int pixel = numPix * psRandomUniform(rng);
     120                int x = pixel % numCols;
     121                int y = pixel / numCols;
     122                if ((mask && (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal)) ||
     123                    !isfinite(image->data.F32[y][x]) ||
     124                    (outMask && (outMask->data.PS_TYPE_MASK_DATA[y][x] & maskVal))) {
     125                    continue;
     126                }
     127
     128                values->data.F32[valueIndex++] = image->data.F32[y][x];
     129            }
     130
     131            if (!chipStats) {
     132                values->n = valueIndex;
     133                if (!psVectorStats(statistics, values, NULL, NULL, 0)) {
     134                    psError(PS_ERR_UNKNOWN, false, "Unable to do statistics on readout.");
     135                    psFree(inView);
     136                    psFree(readout);
     137                    goto MERGE_MASK_ERROR;
     138                }
     139
     140                if (!pmMaskFlagSuspectPixels(outRO, readout, psStatsGetValue(statistics, meanStat),
     141                                             psStatsGetValue(statistics, stdevStat), maskSuspect, maskVal)) {
     142                    psError(PS_ERR_UNKNOWN, false, "Unable to find suspect values in file %d", i);
     143                    psFree(inView);
     144                    psFree(readout);
     145                    goto MERGE_MASK_ERROR;
     146                }
     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            psFree(outRO);
     157        }
     158
     159        // Additional run through cells if we want chip-level statistics
     160        if (chipStats && valueIndex > 0) {
     161            values->n = valueIndex;
     162            if (!psVectorStats(statistics, values, NULL, NULL, 0)) {
     163                psError(PS_ERR_UNKNOWN, false, "Unable to do statistics on chip.");
     164                goto MERGE_MASK_ERROR;
     165            }
     166            inView->cell = -1;
     167            while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) {
     168                pmHDU *hdu = pmHDUFromCell(inCell); // HDU for cell
     169                if (!hdu || hdu->blankPHU) {
     170                    // No data here
     171                    continue;
     172                }
     173                pmReadout *readout = inCell->readouts->data[0]; // Readout of interest
     174
     175                inView->readout = 0;
     176                pmReadout *outRO = pmFPAfileThisReadout(config->files, inView, "PPMERGE.OUTPUT.MASK");
     177
     178                if (!pmMaskFlagSuspectPixels(outRO, readout, psStatsGetValue(statistics, meanStat),
     179                                             psStatsGetValue(statistics, stdevStat), maskSuspect, maskVal)) {
     180                    psError(PS_ERR_UNKNOWN, false, "Unable to find suspect values in file %d", i);
     181                    goto MERGE_MASK_ERROR;
     182                }
     183                pmCellFreeData(inCell);
     184
     185                inView->readout = -1;
     186                if (!pmFPAfileIOChecks(config, inView, PM_FPA_AFTER)) {
     187                    psFree(inView);
     188                    goto MERGE_MASK_ERROR;
     189                }
     190            }
     191        }
     192    }
     193    psFree(inView);
     194    psFree(statistics); statistics = NULL;
     195    psFree(values); values = NULL;
     196
     197
     198    // Another run through the chip to threshold on the suspects
     199    pmFPAfileActivate(config->files, false, NULL);
     200    if (writeOut) {
     201        ppMergeFileActivate(config, PPMERGE_FILES_OUTPUT, true);
     202    }
     203    pmFPA *outFPA = outChip->parent;    // Output FPA
     204    pmCell *outCell;                    // Output cell
     205    pmFPAview *outView = pmFPAviewAlloc(0); // View into output FPA
     206    *outView = *view;
     207    while ((outCell = pmFPAviewNextCell(outView, outFPA, 1))) {
     208        pmHDU *hdu = pmHDUFromCell(outCell); // HDU for cell
     209        if (!hdu || hdu->blankPHU) {
     210            // No data here
     211            continue;
     212        }
     213
     214        psTrace("ppMerge", 1, "Getting bad pixels for chip %d cell %d", outView->chip, outView->cell);
     215
     216        assert(outCell->readouts && outCell->readouts->n == 1);
     217        pmReadout *outRO = outCell->readouts->data[0]; // Output readout
     218        if (!pmMaskIdentifyBadPixels(outRO, maskVal, maskBad, maskMode)) {
     219            psError(PS_ERR_UNKNOWN, false, "Unable to mask bad pixels");
     220            goto MERGE_MASK_ERROR;
     221        }
     222
     223        // Supplementary outputs
     224        {
     225            // The counts image is fairly useless, but it preserves the model
     226            pmCell *countsCell = pmFPAfileThisCell(config->files, outView, "PPMERGE.OUTPUT.COUNT");
     227            pmReadout *countsRO = pmReadoutAlloc(countsCell); // Readout with count of inputs per pixel
     228            countsRO->image = psImageAlloc(outRO->mask->numCols, outRO->mask->numRows, PS_TYPE_F32);
     229            psImageInit(countsRO->image, numFiles);
     230            countsRO->data_exists = countsCell->data_exists = countsCell->parent->data_exists = true;
     231            psFree(countsRO);
     232
     233            pmCell *sigmaCell = pmFPAfileThisCell(config->files, outView, "PPMERGE.OUTPUT.SIGMA");
     234            pmReadout *sigmaRO = pmReadoutAlloc(sigmaCell); // Readout with suspect image
     235            psImage *suspect = psMetadataLookupPtr(NULL, outRO->analysis, PM_MASK_ANALYSIS_SUSPECT);
     236            sigmaRO->image = psImageCopy(sigmaRO->image, suspect, PS_TYPE_F32);
     237            psMetadataRemoveKey(outRO->analysis, PM_MASK_ANALYSIS_SUSPECT);
     238            sigmaRO->data_exists = sigmaCell->data_exists = sigmaCell->parent->data_exists = true;
     239            psFree(sigmaRO);
     240        }
     241
     242        if (maskGrowVal > 0) {
     243            psImage *grown = psImageGrowMask(NULL, outRO->mask, maskVal, maskGrow, maskGrowVal); // Grown mask
     244            psFree(outRO->mask);
     245            outRO->mask = grown;
     246        }
     247
     248        if (writeOut) {
     249            if (!pmFPAfileIOChecks(config, outView, PM_FPA_BEFORE)) {
     250                psFree(outView);
     251                goto MERGE_MASK_ERROR;
     252            }
     253
     254            outRO->data_exists = outCell->data_exists = outChip->data_exists = true;
     255
     256            // Average concepts
     257            psList *cells = psListAlloc(NULL); // List of cells, for concept averaging
     258            for (int i = 0; i < numFiles; i++) {
     259                pmFPAfile *inFile = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file
     260                pmCell *inCell = pmFPAviewThisCell(outView, inFile->fpa); // Input cell
     261                psListAdd(cells, PS_LIST_TAIL, inCell);
     262            }
     263            if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, true)) {
     264                psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts.");
     265                psFree(cells);
     266                psFree(outRO);
     267                psFree(outView);
     268                return false;
     269            }
     270            psFree(cells);
     271
     272            // Statistics on the merged cell using a fake image
     273            outRO->image = psImageAlloc(outRO->mask->numCols, outRO->mask->numRows, PS_TYPE_F32);
     274            psImageInit(outRO->image, 1.0);
     275            if (!ppStatsFPA(stats, outRO->parent->parent->parent, outView,
     276                            maskVal | pmConfigMask("BLANK", config), config)) {
     277                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to generate stats for image.");
     278                psFree(outRO);
     279                psFree(outView);
     280                return false;
     281            }
     282            psFree(outRO->image);
     283            outRO->image = NULL;
     284
     285            // Write
     286            if (!pmFPAfileIOChecks(config, outView, PM_FPA_AFTER)) {
     287                psFree(outView);
     288                goto MERGE_MASK_ERROR;
     289            }
     290        }
     291    }
     292    psFree(outView);
     293
     294    ppMergeFileActivate(config, PPMERGE_FILES_ALL, true);
    50295
    51296    return true;
     297
     298
     299MERGE_MASK_ERROR:
     300    psFree(statistics);
     301    psFree(values);
     302    return false;
    52303}
    53304
    54 bool ppMergeMaskReadoutStats(const pmReadout *readout,
    55                              const pmReadout *output,
    56                              ppMergeOptions *options, // Options
    57                              psRandom *rng)
     305
     306
     307
     308
     309bool ppMergeMask(pmConfig *config)
    58310{
    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 
     311    assert(config);
     312
     313    bool mdok;                          // Status of MD lookup
     314    int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
     315    bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks?
     316    bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); // Do we have weights?
     317    int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Number of rejection iterations
     318
     319    PS_ASSERT_INT_POSITIVE(iter, false);
     320
     321    pmFPAview *view = pmFPAviewAlloc(0); // View to component of interest
     322    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
     323
     324    psMetadata *stats = NULL;           // Statistics for output
     325    if (psMetadataLookup(config->arguments, "STATS.NAME")) {
     326        stats = psMetadataAlloc();
     327        psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "STATS.DATA", 0, "Statistics output", stats);
     328    }
     329
     330    psArray *inputs = ppMergeFileDataLevel(config, "PPMERGE.INPUT", PM_FPA_LEVEL_READOUT); // Input images
     331    psFree(inputs);
     332    if (haveMasks) {
     333        psArray *masks = ppMergeFileDataLevel(config, "PPMERGE.INPUT.MASK", PM_FPA_LEVEL_READOUT);
     334        psFree(masks);
     335    }
     336    if (haveWeights) {
     337        psArray *weights = ppMergeFileDataLevel(config, "PPMERGE.INPUT.WEIGHT", PM_FPA_LEVEL_READOUT);
     338        psFree(weights);
     339    }
     340
     341    if (!ppMergeFileActivate(config, PPMERGE_FILES_INPUT, true)) {
     342        psError(PS_ERR_UNKNOWN, false, "Unable to activate files.");
     343        goto PPMERGE_MASK_ERROR;
     344    }
     345    if (!ppMergeFileActivate(config, PPMERGE_FILES_OUTPUT, true)) {
     346        psError(PS_ERR_UNKNOWN, false, "Unable to activate files.");
     347        goto PPMERGE_MASK_ERROR;
     348    }
     349
     350    psString outName = ppMergeOutputFile(config); // Name of output file
     351    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, outName); // Output file
     352    psFree(outName);
     353    assert(output && output->fpa);
     354    pmFPA *outFPA = output->fpa;        // Output FPA
     355
     356    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     357        goto PPMERGE_MASK_ERROR;
     358    }
     359    pmChip *outChip;                    // Chip of interest
     360    while ((outChip = pmFPAviewNextChip(view, outFPA, 1))) {
     361        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     362            goto PPMERGE_MASK_ERROR;
     363        }
     364
     365        for (int i = 0; i < iter; i++) {
     366            if (!mergeMask(config, view, (i == iter - 1), rng, stats)) {
     367                psError(PS_ERR_UNKNOWN, false, "Unable to merge chip %d", view->chip);
     368                goto PPMERGE_MASK_ERROR;
     369            }
     370        }
     371
     372        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     373            goto PPMERGE_MASK_ERROR;
     374        }
     375    }
     376
     377    psList *fpaList = psListAlloc(NULL);// List of FPAs for concept averaging
     378    for (int i = 0; i < numFiles; i++) {
     379        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file
     380        psListAdd(fpaList, PS_LIST_TAIL, file->fpa);
     381    }
     382    if (!pmConceptsAverageFPAs(outFPA, fpaList)) {
     383        psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts.");
     384        psFree(fpaList);
     385        goto PPMERGE_MASK_ERROR;
     386    }
     387    psFree(fpaList);
     388
     389    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     390        goto PPMERGE_MASK_ERROR;
     391    }
     392
     393    psFree(stats);
     394    psFree(view);
    102395    psFree(rng);
     396
    103397    return true;
     398
     399PPMERGE_MASK_ERROR:
     400    psFree(stats);
     401    psFree(view);
     402    psFree(rng);
     403    return false;
    104404}
    105405
    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.