Index: /trunk/ppMerge/src/ppMergeMaskByImageStats.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskByImageStats.c (revision 15804) +++ /trunk/ppMerge/src/ppMergeMaskByImageStats.c (revision 15804) @@ -0,0 +1,87 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include +#include + +#include "ppMerge.h" +#include "ppMergeData.h" +#include "ppMergeMask.h" + + +// Generate a mask image consisting of pixels which are outliers from the image mean +bool ppMergeMaskByImageStats(ppMergeData *data, // Data + ppMergeOptions *options, // Options + pmConfig *config // Configuration + ) +{ + psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names + assert(filenames); // It should be here --- it's put here in ppMergeConfig + + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator + pmFPA *fpaOut = data->out; // Output FPA + + // Iterate over each file + for (int i = 0; i < filenames->n; i++) { + if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { + continue; + } + psFits *fits = data->files->data[i]; // FITS file handle + if (!fits) { + continue; + } + psTrace("ppMerge", 3, "File %d: %s\n", i, (const char*)filenames->data[i]); + + pmFPA *fpaIn = data->in->data[i]; // Input FPA + pmFPAview *view = pmFPAviewAlloc(0); // View of FPA, for iteration + pmChip *chipIn; // Input chip of interest + while ((chipIn = pmFPAviewNextChip(view, fpaIn, 1))) { + pmCell *cellIn; // Input cell of interest + while ((cellIn = pmFPAviewNextCell(view, fpaIn, 1))) { + if (!pmCellRead(cellIn, fits, config->database)) { + continue; + } + if (cellIn->readouts->n == 0) { + continue; + } + + pmCell *cellOut = pmFPAviewThisCell(view, fpaOut); // Output cell + // Suspect pixels image + psImage *suspect = psMetadataLookupPtr(NULL, cellOut->analysis, "MASK.SUSPECT"); + bool first = suspect ? false : true; + + pmReadout *roIn; // Input readout of interest + while ((roIn = pmFPAviewNextReadout(view, fpaIn, 1))) { + if (!roIn->image) { + continue; + } + psTrace("ppMerge", 4, "Flagging suspect pixels in chip %d, cell %d, ro %d\n", + view->chip, view->cell, view->readout); + float frac = options->sample / (float)(roIn->image->numCols * roIn->image->numRows); + suspect = pmMaskFlagSuspectPixels(suspect, roIn, options->maskSuspect, + options->combine->maskVal, frac, rng); + } + + if (first) { + psMetadataAddImage(cellOut->analysis, PS_LIST_TAIL, "MASK.SUSPECT", 0, + "Suspect pixels", suspect); + psFree(suspect); + } + + pmCellFreeData(cellIn); + } + pmChipFreeData(chipIn); + } + pmFPAFreeData(fpaIn); + psFree(view); + } + psFree(rng); + + return true; +} Index: /trunk/ppMerge/src/ppMergeMaskByPixelStats.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskByPixelStats.c (revision 15804) +++ /trunk/ppMerge/src/ppMergeMaskByPixelStats.c (revision 15804) @@ -0,0 +1,177 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include +#include + +#include "ppMerge.h" +#include "ppMergeData.h" +#include "ppMergeCombine.h" +#include "ppMergeVersion.h" + +// XXX for the moment, this function and ppMergeMaskByImageStats both load (and free) the input +// data + +// Combine the inputs +bool ppMergeMaskByPixelStats(ppMergeData *data, // Data + ppMergeOptions *options, // Options + pmConfig *config // Configuration + ) +{ + psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names + assert(filenames); // It should be here --- it's put here in ppMergeConfig + + // Iterate over the FPA + pmFPA *fpa = data->out; // Output FPA + pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration + int cellNum = -1; // Cell number in the whole FPA + if (data->out->hdu) { + pmFPAUpdateNames(data->out, NULL, NULL); + } + + pmChip *chip; // Chip of interest + pmHDU *lastHDU = NULL; // Last HDU to be updated + + while ((chip = pmFPAviewNextChip(view, fpa, 1))) { + if (chip->hdu) { + // Data will exist soon + pmFPAUpdateNames(data->out, chip, NULL); + chip->data_exists = true; + } + pmCell *cell; // Cell of interest + while ((cell = pmFPAviewNextCell(view, fpa, 1))) { + cellNum++; + if (cell->hdu) { + // Data will exist soon + pmFPAUpdateNames(data->out, chip, cell); + chip->data_exists = cell->data_exists = true; + } + pmReadout *readout = pmReadoutAlloc(cell); // Output readout of interest + psArray *stack = psArrayAlloc(filenames->n); // Stack of readouts to combine + + // Read bit by bit + int numRead; // Number of inputs read + int numScan = 0; + + // Put version metadata into header + pmHDU *hdu = pmHDUFromCell(cell); + if (hdu && hdu != lastHDU) { + if (!hdu->header) { + hdu->header = psMetadataAlloc(); + } + ppMergeVersionMetadata(hdu->header); + lastHDU = hdu; + } + + while (1) { + numRead = 0; + for (int i = 0; i < filenames->n; i++) { + if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { + continue; + } + psFits *fits = data->files->data[i]; // FITS file handle + if (!fits) { + continue; + } + + if (!stack->data[i]) { + pmFPA *fpaIn = data->in->data[i]; // Input FPA + pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip + pmCell *cellIn = chipIn->cells->data[view->cell]; // Input cell + stack->data[i] = pmReadoutAlloc(cellIn); // Input readout + } + + // Only reading and writing the first readout in each cell (plane 0) + bool readOK; + if (pmReadoutReadNext(&readOK, stack->data[i], fits, 0, options->rows)) { + if (!readOK) { + psError(PS_ERR_IO, false, "Failed to read concepts for cell.\n"); + psErrorStackPrint(stderr, "trouble reading data!\n"); + exit (1); + } + numRead++; + } else { + psTrace("ppMerge", 3, "Unable to read from file %d for chip %d, " + "cell %d, scan %d\n", i, view->chip, view->cell, numScan); + } + } + + psTrace("ppMerge", 5, "Chip %d, cell %d, scan %d\n", view->chip, view->cell, numScan); + if (numRead == 0) { + break; + } + ppMergeMaskReadoutByPixelStats(readout, stack, options->combine); + numScan++; + } + + // Get list of cells for concepts averaging + // XXX only do this operation once between Image and Pixel Stats + psList *inCells = psListAlloc(NULL); // List of cells + for (int i = 0; i < filenames->n; i++) { + if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { + continue; + } + pmCell *cellIn = pmFPAviewThisCell(view, data->in->data[i]); // Input cell + psListAdd(inCells, PS_LIST_TAIL, cellIn); + } + if (!pmConceptsAverageCells(cell, inCells, NULL, NULL, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts."); + psFree(inCells); + return false; + } + psFree(inCells); + psFree(stack); + psFree(readout); // Drop reference + + // Blow away the cell data + for (int i = 0; i < filenames->n; i++) { + pmFPA *fpaIn = data->in->data[i]; // Input FPA + if (!fpaIn) { continue; } // was not a valid input file + pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip + pmCell *cellIn = chipIn->cells->data[view->cell]; // Input cell + pmCellFreeData(cellIn); + } + } + + // Blow away the chip data + for (int i = 0; i < filenames->n; i++) { + pmFPA *fpaIn = data->in->data[i]; // Input FPA + if (!fpaIn) { continue; } // was not a valid input file + pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip + pmChipFreeData(chipIn); + } + } + + // Blow away the FPA data + for (int i = 0; i < filenames->n; i++) { + pmFPA *fpaIn = data->in->data[i]; // Input FPA + if (!fpaIn) { continue; } // was not a valid input file + pmFPAFreeData(fpaIn); + } + + psFree(view); + + // Get list of FPAs for concepts averaging + psList *inFPAs = psListAlloc(NULL); // List of FPAs + for (int i = 0; i < filenames->n; i++) { + if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { + continue; + } + pmFPA *fpaIn = data->in->data[i]; // Input FPA + psListAdd(inFPAs, PS_LIST_TAIL, fpaIn); + } + if (!pmConceptsAverageFPAs(fpa, inFPAs)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts."); + psFree(inFPAs); + return false; + } + psFree(inFPAs); + + return true; +} Index: /trunk/ppMerge/src/ppMergeMaskOutput.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskOutput.c (revision 15804) +++ /trunk/ppMerge/src/ppMergeMaskOutput.c (revision 15804) @@ -0,0 +1,157 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include +#include + +#include "ppMerge.h" +#include "ppMergeData.h" +#include "ppMergeMask.h" + + +// write out the mask: assumes it has been saved on the cell->analysis as MASK.UNION +bool ppMergeMaskOutput(ppMergeData *data, // Data + ppMergeOptions *options, // Options + pmConfig *config // Configuration + ) +{ + psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names + assert(filenames); // It should be here --- it's put here in ppMergeConfig + + pmFPA *fpaOut = data->out; // Output FPA + + pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration + if (fpaOut->hdu) { + pmFPAUpdateNames(fpaOut, NULL, NULL); + } + pmFPAWriteMask(fpaOut, data->outFile, config->database, true, false); // Write header only + pmChip *chipOut; // Output chip of interest + while ((chipOut = pmFPAviewNextChip(view, fpaOut, 1))) { + if (chipOut->hdu) { + chipOut->data_exists = true; + pmFPAUpdateNames(fpaOut, chipOut, NULL); + } + pmChipWriteMask(chipOut, data->outFile, config->database, true, false); // Write header only + pmCell *cellOut; // Output cell of interest + while ((cellOut = pmFPAviewNextCell(view, fpaOut, 1))) { + if (cellOut->hdu) { + chipOut->data_exists = cellOut->data_exists = true; + pmFPAUpdateNames(fpaOut, chipOut, cellOut); + } + pmCellWriteMask(cellOut, data->outFile, config->database, true); // Write header only + + psImage *maskUnion = psMetadataLookupPtr(NULL, cellOut->analysis, "MASK.UNION"); + if (! maskUnion) { + continue; + } + + pmReadout *roOut = pmReadoutAlloc(cellOut); // Output readout + roOut->mask = maskUnion; + roOut->data_exists = cellOut->data_exists = chipOut->data_exists = true; + + // XXX skip this? Get list of cells for concepts averaging + { + psList *inCells = psListAlloc(NULL); // List of cells + for (int i = 0; i < filenames->n; i++) { + if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { + continue; + } + pmCell *cellIn = pmFPAviewThisCell(view, data->in->data[i]); // Input cell + psListAdd(inCells, PS_LIST_TAIL, cellIn); + } + if (!pmConceptsAverageCells(cellOut, inCells, NULL, NULL, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts."); + psFree(inCells); + return false; + } + psFree(inCells); + } + + { + // Add MD5 information for cell + pmHDU *hdu = pmHDUFromCell(cell); // HDU that owns the cell + const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); + + psString headerName = NULL; // Header name for MD5 + psStringAppend(&headerName, "MD5_%s_%s", chipName, cellName); + + psVector *md5 = psImageMD5(readout->image); // md5 hash + psString md5string = psMD5toString(md5); // String + psFree(md5); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, + "Image MD5", md5string); + psFree(md5string); + psFree(headerName); + } + + // Statistics on the merged cell + if (data->statsFile) { + if (!data->stats) { + data->stats = psMetadataAlloc(); + } + + // Build a fake image to do statistics + roOut->image = psImageAlloc(roOut->mask->numCols, roOut->mask->numRows, PS_TYPE_F32); + psImageInit(roOut->image, 1.0); + if (!ppStatsFPA(data->stats, data->out, view, + options->combine->maskVal | pmConfigMask("BLANK", config), + config)) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to generate stats for image.\n"); + return false; + } + psFree(roOut->image); + roOut->image = NULL; + } + + psFree(roOut); // Drop reference + + if (cellOut->hdu && !cellOut->hdu->blankPHU) { + psTrace("ppMerge", 5, "Writing out cell HDU.\n"); + pmCellWriteMask(cellOut, data->outFile, config->database, false); + pmCellFreeData(cellOut); + } + } + + if (chipOut->hdu && !chipOut->hdu->blankPHU) { + psTrace("ppMerge", 5, "Writing out chip HDU.\n"); + pmChipWriteMask(chipOut, data->outFile, config->database, false, false); + pmChipFreeData(chipOut); + } + } + + // Get list of FPAs for concepts averaging + { + psList *inFPAs = psListAlloc(NULL); // List of FPAs + for (int i = 0; i < filenames->n; i++) { + if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { + continue; + } + pmFPA *fpaIn = data->in->data[i]; // Input FPA + psListAdd(inFPAs, PS_LIST_TAIL, fpaIn); + } + + if (!pmConceptsAverageFPAs(fpaOut, inFPAs)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts."); + psFree(inFPAs); + return false; + } + psFree(inFPAs); + } + + if (fpaOut->hdu && !fpaOut->hdu->blankPHU) { + psTrace("ppMerge", 5, "Writing out FPA HDU.\n"); + pmFPAWriteMask(fpaOut, data->outFile, config->database, false, false); + } + pmFPAFreeData(fpaOut); + + psFree(view); + + return true; +} Index: /trunk/ppMerge/src/ppMergeMaskReadoutByPixelStats.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskReadoutByPixelStats.c (revision 15804) +++ /trunk/ppMerge/src/ppMergeMaskReadoutByPixelStats.c (revision 15804) @@ -0,0 +1,327 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include + +// XXX: Maybe add support for S16 and S32 types. Currently, only F32 supported. +bool ppMergeMaskReadoutByPixelStats(pmReadout *output, const psArray *inputs, const pmCombineParams *params) +{ + // Check inputs + PS_ASSERT_PTR_NON_NULL(output, false); + PS_ASSERT_ARRAY_NON_NULL(inputs, false); + PS_ASSERT_PTR_NON_NULL(params, false); + PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracLow, 0.0, 1.0, false); + PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracHigh, 0.0, 1.0, false); + if (params->combine != PS_STAT_SAMPLE_MEAN && params->combine != PS_STAT_SAMPLE_MEDIAN && + params->combine != PS_STAT_ROBUST_MEDIAN && params->combine != PS_STAT_FITTED_MEAN && + params->combine != PS_STAT_CLIPPED_MEAN) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Combination method is not SAMPLE_MEAN, SAMPLE_MEDIAN, " + "ROBUST_MEDIAN, FITTED_MEAN or CLIPPED_MEAN.\n"); + return false; + } + + // XXX is it possible / desired to use the weight in this analysis? + for (int i = 0; i < inputs->n; i++) { + pmReadout *readout = inputs->data[i]; // Readout of interest + if (params->weights && !readout->weight) { + psError(PS_ERR_UNEXPECTED_NULL, true, + "Rejection based on weights requested, but no weights supplied for image %d.\n", i); + return false; + } + } + + bool first = !output->image; // First pass through? + + pmHDU *hdu = pmHDUFromReadout(output); // Output HDU + if (!hdu) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find HDU for readout.\n"); + return false; + } + + if (first) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Combining using statistic: %x", params->combine); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + + psStats *stats = psStatsAlloc(params->combine); // The statistics to use in the combination + if (params->combine == PS_STAT_CLIPPED_MEAN) { + stats->clipSigma = params->rej; + stats->clipIter = params->iter; + + if (first) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", + params->iter, params->rej); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + } + + int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine + int xSize, ySize; // Size of the output image + if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, + inputs)) { + psError(PS_ERR_UNKNOWN, false, "No valid input readouts."); + return false; + } + + pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true); + psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0); + + psStatsOptions combineStdev = 0; // Statistics option for weights + if (params->weights) { + + if (!output->weight) { + output->weight = psImageAlloc(xSize, ySize, PS_TYPE_F32); + } + if (output->weight->numCols < xSize || output->weight->numRows < ySize) { + psImage *newWeight = psImageAlloc(xSize, ySize, PS_TYPE_F32); + psImageInit(newWeight, 0.0); + psImageOverlaySection(newWeight, output->weight, output->col0, output->row0, "="); + psFree(output->weight); + output->weight = newWeight; + } + + if (first) { + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, + "Using input weights to combine images", ""); + } + + // Get the correct statistics option for weights + switch (params->combine) { + case PS_STAT_SAMPLE_MEAN: + case PS_STAT_SAMPLE_MEDIAN: + combineStdev = PS_STAT_SAMPLE_STDEV; + break; + case PS_STAT_ROBUST_MEDIAN: + combineStdev = PS_STAT_ROBUST_STDEV; + break; + case PS_STAT_FITTED_MEAN: + combineStdev = PS_STAT_FITTED_STDEV; + break; + case PS_STAT_CLIPPED_MEAN: + combineStdev = PS_STAT_CLIPPED_STDEV; + break; + default: + psAbort("Should never get here --- checked params->combine before.\n"); + } + stats->options |= combineStdev; + } + + // We loop through each pixel in the output image. We loop through each input readout. We determine if + // that output pixel is contained in the image from that readout. If so, we save it in psVector pixels. + // If not, we set a mask for that element in pixels. Then, we mask off pixels not between fracLow and + // fracHigh. Then we call the vector stats routine on those pixels/mask. Then we set the output pixel + // value to the result of the stats call. + + psVector *pixels = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of pixels + psF32 *pixelsData = pixels->data.F32; // Dereference pixels + + psVector *mask = psVectorAlloc(inputs->n, PS_TYPE_U8); // Mask for stack + psU8 *maskData = mask->data.U8; // Dereference mask + + psVector *weights = NULL; // Stack of weights + psVector *errors = NULL; // Stack of errors (sqrt of variance/weights), for psVectorStats + psF32 *weightsData = NULL; // Dereference weights + if (params->weights) { + weights = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of weights + weightsData = weights->data.F32; + } + psVector *index = NULL; // The indices to sort the pixels + + float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of pixels to keep + if (keepFrac != 1.0 && first) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Min/max rejection: %f high, %f low, keep %d", + params->fracHigh, params->fracLow, params->nKeep); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + + psMaskType maskVal = params->maskVal; // The mask value + if (maskVal && first) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Mask for combination: %x", maskVal); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + + #ifndef PS_NO_TRACE + psTrace("psModules.imcombine", 3, "Iterating output: %d --> %d, %d --> %d\n", + minInputCols - output->col0, maxInputCols - output->col0, + minInputRows - output->row0, maxInputRows - output->row0); + if (psTraceGetLevel("psModules.imcombine") >= 3) { + for (int r = 0; r < inputs->n; r++) { + pmReadout *readout = inputs->data[r]; // Input readout + psTrace("psModules.imcombine", 3, "Iterating input %d: %d --> %d, %d --> %d\n", r, + minInputCols - readout->col0, maxInputCols - readout->col0, + minInputRows - readout->row0, maxInputRows - readout->row0); + } + } + #endif + + // Dereference output products + psF32 **outputImage = output->image->data.F32; // Output image + psU8 **outputMask = output->mask->data.U8; // Output mask + psF32 **outputWeight = NULL; // Output weight map + if (output->weight) { + outputWeight = output->weight->data.F32; + } + + psVector *invScale = NULL; // Inverse scale; pre-calculated for efficiency + if (scale) { + invScale = (psVector*)psBinaryOp(NULL, psScalarAlloc(1.0, PS_TYPE_F32), "/", (const psPtr)scale); + } + + for (int i = minInputRows; i < maxInputRows; i++) { + int yOut = i - output->row0; // y position on output readout + #ifdef SHOW_BUSY + + if (psTraceGetLevel("psModules.imcombine") > 9) { + printf("Processing row %d\r", i); + fflush(stdout); + } + #endif + for (int j = minInputCols; j < maxInputCols; j++) { + int xOut = j - output->col0; // x position on output readout + + int numValid = 0; // Number of valid pixels in the stack + memset(maskData, 0, mask->n * sizeof(psU8)); // Reset the mask + for (int r = 0; r < inputs->n; r++) { + pmReadout *readout = inputs->data[r]; // Input readout + int yIn = i - readout->row0; // y position on input readout + int xIn = j - readout->col0; // x position on input readout + psImage *image = readout->image; // The readout image + + #if 0 // This should have been taken care of already: + // Check bounds + if (xIn < 0 || xIn >= image->numCols || yIn < 0 || yIn >= image->numRows) { + continue; + } + #endif + + pixelsData[r] = image->data.F32[yIn][xIn]; + if (!isfinite(pixelsData[r])) { + maskData[r] = 1; + continue; + } + + // Check mask + psImage *roMask = readout->mask; // The mask image + if (roMask && roMask->data.U8[yIn][xIn] & maskVal) { + maskData[r] = 1; + continue; + } + + if (params->weights) { + weightsData[r] = readout->weight->data.F32[yIn][xIn]; + } + + if (zero) { + pixelsData[r] -= zero->data.F32[r]; + } + if (scale) { + pixelsData[r] *= invScale->data.F32[r]; + if (params->weights) { + weightsData[r] *= invScale->data.F32[r] * invScale->data.F32[r]; + } + } + + numValid++; + } + + if (numValid == 0) { + outputMask[yOut][xOut] = params->blank; + outputImage[yOut][xOut] = NAN; + continue; + } + + // Apply fracLow,fracHigh if there are enough pixels + if (numValid * keepFrac >= params->nKeep && keepFrac != 1.0) { + index = psVectorSortIndex(index, pixels); + int numLow = numValid * params->fracLow; // Number of low pixels to clip + int numHigh = numValid * params->fracHigh; // Number of high pixels to clip + // Low pixels + psS32 *indexData = index->data.S32; // Dereference index + for (int k = 0, numMasked = 0; numMasked < numLow && k < index->n; k++) { + // Don't count the ones that are already masked + if (!maskData[indexData[k]]) { + maskData[indexData[k]] = 1; + numMasked++; + } + } + // High pixels + for (int k = pixels->n - 1, numMasked = 0; numMasked < numHigh && k >= 0; k--) { + // Don't count the ones that are already masked + if (! maskData[indexData[k]]) { + maskData[indexData[k]] = 1; + numMasked++; + } + } + } + + // XXXXX this step probably is very expensive : convert errors to variance everywhere? + if (params->weights) { + errors = (psVector*)psUnaryOp(errors, weights, "sqrt"); + } + + // Combination + if (!psVectorStats(stats, pixels, errors, mask, 1)) { + // Can't do much about it, but it's not worth worrying about + psErrorClear(); + outputImage[yOut][xOut] = NAN; + outputMask[yOut][xOut] = params->blank; + if (params->weights) { + outputWeight[yOut][xOut] = NAN; + } + } else { + outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine); + if (!isfinite(outputImage[yOut][xOut])) { + outputMask[yOut][xOut] = params->blank; + } + if (params->weights) { + float stdev = psStatsGetValue(stats, combineStdev); + outputWeight[yOut][xOut] = PS_SQR(stdev); // Variance + // XXXX this is not the correct formal error. + // also, the weighted mean is not obviously the correct thing here + } + } + } + } + #ifdef SHOW_BUSY + if (psTraceGetLevel("psModules.imcombine") > 9) { + printf("\n"); + } + #endif + psFree(index); + psFree(pixels); + psFree(mask); + psFree(weights); + psFree(errors); + psFree(stats); + psFree(invScale); + + // Update the "concepts" + psList *inputCells = psListAlloc(NULL); // List of cells + for (long i = 0; i < inputs->n; i++) { + pmReadout *readout = inputs->data[i]; // Readout of interest + psListAdd(inputCells, PS_LIST_TAIL, readout->parent); + } + bool success = pmConceptsAverageCells(output->parent, inputCells, NULL, NULL, true); + psFree(inputCells); + + output->data_exists = true; + output->parent->data_exists = true; + output->parent->parent->data_exists = true; + + return success; +} +