Index: trunk/ippScripts/scripts/camera_exp.pl =================================================================== --- trunk/ippScripts/scripts/camera_exp.pl (revision 35074) +++ trunk/ippScripts/scripts/camera_exp.pl (revision 35081) @@ -41,5 +41,5 @@ my ( $exp_tag, $cam_id, $camera, $outroot, $dbname, $reduction, $dvodb, $verbose, $no_update, - $no_op, $redirect, $save_temps, $run_state, $skip_binned, $skip_masks); + $no_op, $redirect, $save_temps, $run_state, $skip_binned, $skip_masks, $bkg_only); GetOptions( 'exp_tag=s' => \$exp_tag, # Exposure identifier @@ -53,4 +53,5 @@ 'skip-binned' => \$skip_binned, # override recipe - don't create binned images 'skip-refmask' => \$skip_masks, # override recipe - don't create refmask + 'bkg-only' => \$bkg_only, # override recipe - only do background continuity 'verbose' => \$verbose, # Print to stdout 'no-update' => \$no_update, # Update the database? @@ -111,4 +112,5 @@ &my_die("Unrecognised PSASTRO recipe", $cam_id, $PS_EXIT_CONFIG_ERROR) unless defined $recipe_psastro; +my $bkg_recipe = 'PPIMAGE_BKGCONT'; # Add to reduction? my $mdcParser = PS::IPP::Metadata::Config->new; # Parser for metadata config files @@ -215,6 +217,10 @@ my ($list4File, $list4Name) = tempfile( "/tmp/$exp_tag.cm.$cam_id.b4.list.XXXX", UNLINK => !$save_temps ); # For astrometry +### Create temp file for background models +my ($list5File, $list5Name) = tempfile( "/tmp/$exp_tag.cm.$cam_id.b5.list.XXXX", UNLINK => !$save_temps ); # For background models + # XXX we perform astrometry iff photometry output exists my @outMasks; # Names of output masks +my @bkg_models; # Names of output background models foreach my $file (@$files) { next if $file->{quality} != 0; @@ -226,11 +232,13 @@ my $chipObjects = $ipprc->filename("PSPHOT.OUTPUT", $file->{path_base}, $class_id); my $chipMask = $ipprc->filename("PPIMAGE.CHIP.MASK", $file->{path_base}, $class_id); - + print $list1File ($ipprc->filename("PPIMAGE.BIN1", $file->{path_base}, $class_id) . "\n"); print $list2File ($ipprc->filename("PPIMAGE.BIN2", $file->{path_base}, $class_id) . "\n"); print $list3File ($chipObjects . "\n"); print $list4File ($chipMask . "\n"); + print $list5File ($ipprc->filename("PSPHOT.BACKMDL", $file->{path_base}, $class_id) . "\n"); push @outMasks, prepare_output("PSASTRO.OUTPUT.MASK", $outroot, $class_id, 1) if $produceMasks; + push @bkg_models, prepare_output("PPIMAGE.BACKMDL", $outroot, $class_id, 1); } close $list1File; @@ -238,5 +246,5 @@ close $list3File; close $list4File; - +close $list5File; # Output products $ipprc->outroot_prepare($outroot); @@ -291,5 +299,5 @@ } - { + if (!$bkg_only) { # run psastro on the chipObjects, producing fpaObjects my $command; @@ -353,4 +361,22 @@ } } + } + # Construct FPA continuity corrected background images + { + my $command; + $command = "$ppImage"; + $command .= " -list $list5Name $outroot"; + $command .= " -recipe PPIMAGE $bkg_recipe"; + $command .= " -dbname $dbname" if defined $dbname; + + my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = + run(command => $command, verbose => $verbose); + unless ($success) { + $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); + &my_die("Unable to perform ppImage to fix background: $error_code", $cam_id, $error_code); + } + foreach my $bkgModel (@bkg_models) { + check_output($bkgModel, $replicateOutputs); + } } } Index: trunk/ippTools/share/stacktool_tobkg.sql =================================================================== --- trunk/ippTools/share/stacktool_tobkg.sql (revision 35074) +++ trunk/ippTools/share/stacktool_tobkg.sql (revision 35081) @@ -25,5 +25,6 @@ -- WHERE hook %s GROUP BY stack_id -HAVING (SUM(IF(warpRun.state = 'full', 1, 0)) = COUNT(stackInputSkyfile.warp_id) AND +HAVING ((SUM(IF(warpRun.state = 'cleaned', 1, 0)) = COUNT(stackInputSkyfile.warp_id) OR + SUM(IF(warpRun.state = 'full', 1, 0)) = COUNT(stackInputSkyfile.warp_id)) AND SUM(IF(warpSkyfile.background_model = 1, 1, 0)) = COUNT(stackInputSkyfile.warp_id)) Index: trunk/ippconfig/recipes/filerules-mef.mdc =================================================================== --- trunk/ippconfig/recipes/filerules-mef.mdc (revision 35074) +++ trunk/ippconfig/recipes/filerules-mef.mdc (revision 35081) @@ -214,5 +214,5 @@ PPIMAGE.PATTERN OUTPUT {OUTPUT}.ptn PATTERN NONE CHIP TRUE MEF PPIMAGE.STATS OUTPUT {OUTPUT}.stats STATS NONE FPA TRUE MEF - +PPIMAGE.BACKMDL OUTPUT {OUTPUT}.{CHIP.NAME}.mdl.fits IMAGE NONE CHIP TRUE NONE ## note: these use the same output naming convention since they are used for different ppMerge runs Index: trunk/ippconfig/recipes/filerules-simple.mdc =================================================================== --- trunk/ippconfig/recipes/filerules-simple.mdc (revision 35074) +++ trunk/ippconfig/recipes/filerules-simple.mdc (revision 35081) @@ -178,5 +178,5 @@ PPIMAGE.PATTERN OUTPUT {OUTPUT}.ptn PATTERN NONE FPA TRUE NONE PPIMAGE.STATS OUTPUT {OUTPUT}.stats STATS NONE FPA TRUE NONE - +PPIMAGE.BACKMDL OUTPUT {OUTPUT}.{CHIP.NAME}.mdl.fits IMAGE NONE CHIP TRUE NONE ## note: these use the same output naming convention since they are used for different ppMerge runs Index: trunk/ippconfig/recipes/filerules-split.mdc =================================================================== --- trunk/ippconfig/recipes/filerules-split.mdc (revision 35074) +++ trunk/ippconfig/recipes/filerules-split.mdc (revision 35081) @@ -198,4 +198,5 @@ PPIMAGE.PATTERN OUTPUT {OUTPUT}.{CHIP.NAME}.ptn PATTERN NONE CHIP TRUE NONE PPIMAGE.STATS OUTPUT {OUTPUT}.{CHIP.NAME}.stats STATS NONE CHIP TRUE NONE +PPIMAGE.BACKMDL OUTPUT {OUTPUT}.{CHIP.NAME}.mdl.fits IMAGE NONE CHIP TRUE NONE ## note: these use the same output naming convention since they are used for different ppMerge runs Index: trunk/ippconfig/recipes/ppImage.config =================================================================== --- trunk/ippconfig/recipes/ppImage.config (revision 35074) +++ trunk/ippconfig/recipes/ppImage.config (revision 35081) @@ -42,5 +42,5 @@ MASK.STATS BOOL FALSE # calculate mask statistics. GAIN.OVERRIDE BOOL TRUE # Override a non-finite gain? - +BACKGROUND.CONTINUITY BOOL FALSE # Apply continuity correction to background models for an FPA. TILTYSTREAK.BY.CLASS METADATA # apply the 'tiltystreak' tool END @@ -1147,4 +1147,38 @@ MASK.SATURATED BOOL FALSE # DO NOT Mask the saturated pixels MASK.LOW BOOL FALSE # DO NOT Mask the saturated pixels +END + +# Calculate new background models with continuity correction applied +PPIMAGE_BKGCONT METADATA + OVERSCAN BOOL FALSE # Overscan subtraction + NONLIN BOOL FALSE # Non-linearity correction; not implemented + BIAS BOOL FALSE # Bias subtraction + DARK BOOL FALSE # Dark subtraction + SHUTTER BOOL FALSE # Shutter correction + FLAT BOOL FALSE # Flat-field normalisation + MASK BOOL FALSE # Mask bad pixels + MASK.BUILD BOOL FALSE # Build internal mask? + FRINGE BOOL FALSE # Fringe subtraction + PHOTOM BOOL FALSE # Source identification and photometry + ASTROM.CHIP BOOL FALSE # Astrometry per chip? + ASTROM.MOSAIC BOOL FALSE # Astrometry for mosaic? + BASE.FITS BOOL FALSE # Save base image? + BASE.MASK.FITS BOOL FALSE # Save base detrended image? + BASE.VARIANCE.FITS BOOL FALSE # Save base detrended image? + CHIP.FITS BOOL FALSE # Save chip-mosaic-ed image? + CHIP.MASK.FITS BOOL FALSE # Save chip-mosaic-ed image? + CHIP.VARIANCE.FITS BOOL FALSE # Save chip-mosaic-ed image? + BIN1.FITS BOOL FALSE # Save 1st binned chip image? + BIN2.FITS BOOL FALSE # Save 2nd binned chip image? + BIN1.JPEG BOOL FALSE # Save 1st binned jpeg? + BIN2.JPEG BOOL FALSE # Save 2nd binned jpeg? + BIN1.XBIN S32 1 # Image is already binned + BIN1.YBIN S32 1 # Image is already binned + BIN2.XBIN S32 1 # Image is already binned + BIN2.YBIN S32 1 # Image is already binned + MASK.SATURATED BOOL FALSE # DO NOT Mask the saturated pixels + MASK.LOW BOOL FALSE # DO NOT Mask the saturated pixels + BACKGROUND.CONTINUITY BOOL TRUE + PATTERN.CONTINUITY.WIDTH S32 2 END Index: trunk/ppImage/src/Makefile.am =================================================================== --- trunk/ppImage/src/Makefile.am (revision 35074) +++ trunk/ppImage/src/Makefile.am (revision 35081) @@ -44,4 +44,5 @@ ppImageRebinReadout.c \ ppImageMosaic.c \ + ppImageMosaicBackground.c \ ppImageMaskStats.c \ ppImagePhotom.c \ Index: trunk/ppImage/src/ppImage.h =================================================================== --- trunk/ppImage/src/ppImage.h (revision 35074) +++ trunk/ppImage/src/ppImage.h (revision 35081) @@ -41,4 +41,5 @@ bool doPatternCell; // Cell pattern correction bool doPatternContinuity; // Cell continuity correction + bool doBackgroundContinuity; // Do mosaic continuity correction bool doFringe; // Fringe subtraction bool doPhotom; // Source identification and photometry @@ -50,6 +51,5 @@ bool checkNoise; // measure cell-level variance bool applyParity; // Apply Cell parities - - bool doMaskStats; // Calculate mask statistics + bool doMaskStats; // Calculate mask statistics bool doCrosstalkMeasure; // measure crosstalk signal @@ -183,4 +183,7 @@ bool ppImageDetrendPatternApply(pmConfig *config, pmChip *chip, const pmFPAview *inputView, const ppImageOptions *options); +// Do background continuity step +bool ppImageMosaicBackground(pmConfig *config, const ppImageOptions *options); + // Record which detrend file was used for the detrending bool ppImageDetrendRecord( Index: trunk/ppImage/src/ppImageLoop.c =================================================================== --- trunk/ppImage/src/ppImageLoop.c (revision 35074) +++ trunk/ppImage/src/ppImageLoop.c (revision 35081) @@ -41,5 +41,5 @@ ESCAPE("load failure for FPA"); } - + pmChip *chip; // Chip from FPA while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) { @@ -233,5 +233,8 @@ ESCAPE("Unable to bin chip (level 2)."); } - + if (options->doBackgroundContinuity) { + pmFPAfile *out = psMetadataLookupPtr(NULL, config->files, "PPIMAGE.BACKMDL"); + pmFPAfileCopyView(out->fpa,out->src,view); + } // Close cells (XXX shouldn't pmFPAfileClose iterate down as needed?) view->cell = -1; @@ -251,4 +254,9 @@ } + // Do background model continuity updates + if (!ppImageMosaicBackground(config, options )) { + ESCAPE("failure in Background Mosaic"); + } + // generate the full-scale FPA mosaic if (!ppImageMosaicFPA(config, options, "PPIMAGE.OUTPUT.FPA1", "PPIMAGE.BIN1")) { Index: trunk/ppImage/src/ppImageMosaicBackground.c =================================================================== --- trunk/ppImage/src/ppImageMosaicBackground.c (revision 35081) +++ trunk/ppImage/src/ppImageMosaicBackground.c (revision 35081) @@ -0,0 +1,77 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "ppImage.h" + + +bool ppImageMosaicBackground(pmConfig *config, const ppImageOptions *options) { + assert(config); + assert(options); + + if (options->doBackgroundContinuity) { + bool status; + pmFPAfile *in = psMetadataLookupPtr(&status, config->files, "PPIMAGE.INPUT"); + if (!status) { + psErrorStackPrint(stderr, "Can't find required I/O file!\n"); + exit(EXIT_FAILURE); + } + pmFPAfile *out = psMetadataLookupPtr(&status, config->files, "PPIMAGE.BACKMDL"); + if (!status) { + psErrorStackPrint(stderr, "Can't find required I/O file!\n"); + exit(EXIT_FAILURE); + } + //pmReadout *model = pmFPAGenerateReadout(config, view, psphotGetFilerule("PSPHOT.BACKMDL"), inFPA, binning, index); + + pmFPAview *view = pmFPAviewAlloc(0); + // pmFPAAddSourceFromView(out->fpa, view, out->format); + // psFree(view); + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + // ESCAPE("load failure for Chip"); + } + + pmChip *chip; + + if (!pmPatternContinuityBackground(in,out,options->patternCellBG,options->patternCellMean, + options->maskValue, options->darkMask, options->patternContinuityEdgeWidth)) { + // Free things? + psFree(view); + return(false); + } + + // Write out output models + psFree(view); + view = pmFPAviewAlloc(0); + while ((chip = pmFPAviewNextChip(view, out->fpa, 1)) != NULL) { + pmCell *cell; + while ((cell = pmFPAviewNextCell(view, out->fpa, 1)) != NULL) { + pmReadout *readout; + while ((readout = pmFPAviewNextReadout(view, out->fpa, 1)) != NULL) { + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { + psError(PS_ERR_UNKNOWN, false, "I/O failure in ppImageMosaicBackground"); + psFree(view); + return(false); + } + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(PS_ERR_UNKNOWN, false, "I/O failure in ppImageMosaicBackground"); + psFree(view); + return(false); + } + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(PS_ERR_UNKNOWN, false, "I/O failure in ppImageMosaicBackground"); + psFree(view); + return(false); + } + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(PS_ERR_UNKNOWN, false, "I/O failure in ppImageMosaicBackground"); + psFree(view); + return(false); + } + psFree(view); + } + + return(true); +} Index: trunk/ppImage/src/ppImageOptions.c =================================================================== --- trunk/ppImage/src/ppImageOptions.c (revision 35074) +++ trunk/ppImage/src/ppImageOptions.c (revision 35081) @@ -35,4 +35,5 @@ options->doPatternCell = false; // Cell pattern correction options->doPatternContinuity = false; // Cell continuity correction + options->doBackgroundContinuity = false; // Chip level background continuity correction options->doFringe = false; // Fringe subtraction options->doPhotom = false; // Source identification and photometry @@ -429,4 +430,7 @@ } + // Option to enable the background continuity + options->doBackgroundContinuity = psMetadataLookupBool(NULL, recipe, "BACKGROUND.CONTINUITY"); + // Remnance options Index: trunk/ppImage/src/ppImageParseCamera.c =================================================================== --- trunk/ppImage/src/ppImageParseCamera.c (revision 35074) +++ trunk/ppImage/src/ppImageParseCamera.c (revision 35081) @@ -303,4 +303,5 @@ } + // chipImage -> psphotInput (pmFPAfileDefineFromFile) : fpa is constructed // psphotInput -> psphotOutput (pmFPAfileDefineOutputFromFile) : fpa is equated @@ -387,4 +388,25 @@ // the input data is the same as the outImage data : force the free levels to match input->freeLevel = PS_MIN(outImage->freeLevel, input->freeLevel); + + // define the continuity corrected background model files + if (options->doBackgroundContinuity) { + pmFPAfile *bkgMosaicModel = pmFPAfileDefineFromFPA(config,input->fpa, 1, 1, "PPIMAGE.BACKMDL"); + if (!bkgMosaicModel) { + psError(PS_ERR_IO, false, _("Unable to generate new file from PPIMAGE.BACKMDL")); + psFree(options); + return NULL; + } + if (bkgMosaicModel->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPIMAGE.BACKMDL is not of type IMAGE"); + psFree(options); + return NULL; + + } + bkgMosaicModel->save = options->doBackgroundContinuity; + // bkgMosaicModel->freeLevel = PS_MIN(bkgMosaicModel->freeLevel, PM_FPA_LEVEL_FPA); + bkgMosaicModel->freeLevel = PM_FPA_LEVEL_FPA; + bkgMosaicModel->dataLevel = bkgMosaicModel->dataLevel; + bkgMosaicModel->fileLevel = PS_MIN(bkgMosaicModel->fileLevel, bkgMosaicModel->dataLevel); + } // define the binned target files (which may just be carriers for some camera configurations) @@ -497,5 +519,5 @@ chipVariance->save = false; } - + if (psTraceGetLevel("ppImage.config") > 0) { // Get a look inside all the files. Index: trunk/ppStack/src/ppStackCamera.c =================================================================== --- trunk/ppStack/src/ppStackCamera.c (revision 35074) +++ trunk/ppStack/src/ppStackCamera.c (revision 35081) @@ -47,5 +47,6 @@ } bool convolve = psMetadataLookupBool(NULL, recipe, "CONVOLVE"); // Convolve images before stack? - + bool doBackground = psMetadataLookupBool(NULL, recipe, "BACKGROUND.MODEL"); + psArray *runImages = pmFPAfileDefineMultipleFromRun(&status, NULL, config, "PPSTACK.INPUT"); // Input images from previous run @@ -217,16 +218,25 @@ // Grab bkgmodel information here - if (!bkgmodel || strlen(bkgmodel) == 0) { - // We have no background models. + if (doBackground) { + if (!bkgmodel || strlen(bkgmodel) == 0) { + // We have no background models. + } + else { + pmFPAfile *inputBKG = defineFile(config,NULL,"PPSTACK.INPUT.BKGMODEL",bkgmodel, + PM_FPA_FILE_IMAGE); + if (!inputBKG) { + // We failed to generate an pmFPAfile, so disable the background construction and continue. + doBackground = false; + psMetadataAddBool(recipe, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, + "Do photometry on stacked image?", false); + +#if (0) + psError(psErrorCodeLast(), false, + "Unable to define file from bkgmodel %d (%s)",i,bkgmodel); + return(false); +#endif + } + }// End bkgmodel } - else { - pmFPAfile *inputBKG = defineFile(config,NULL,"PPSTACK.INPUT.BKGMODEL",bkgmodel, - PM_FPA_FILE_IMAGE); - if (!inputBKG) { - psError(psErrorCodeLast(), false, - "Unable to define file from bkgmodel %d (%s)",i,bkgmodel); - return(false); - } - }// End bkgmodel Index: trunk/ppStack/src/ppStackCombineAlternate.c =================================================================== --- trunk/ppStack/src/ppStackCombineAlternate.c (revision 35074) +++ trunk/ppStack/src/ppStackCombineAlternate.c (revision 35081) @@ -198,5 +198,5 @@ offset = (Sy * Sxx - Sx * Sxy) / D; scale = (S * Sxy - Sx * Sy) / D; - fprintf(stderr,"Scales: %d %g %g\n",i,offset,scale); + fprintf(stderr,"Scales: %d %g %g %g %g %g %g %g %g\n",i,offset,scale,D,Sx,Sy,Sxx,Syy,Sxy); // Apply scaling factors for (int y = minInputRows; y < maxInputRows; y++) { Index: trunk/psModules/src/detrend/pmPattern.c =================================================================== --- trunk/psModules/src/detrend/pmPattern.c (revision 35074) +++ trunk/psModules/src/detrend/pmPattern.c (revision 35081) @@ -938,3 +938,366 @@ } - +// This should reuse code more efficiently +bool pmPatternContinuityBackground(pmFPAfile *in, pmFPAfile *out, psStatsOptions bgStat, psStatsOptions cellStat, + psImageMaskType maskVal, psImageMaskType maskBad, int edgeWidth) +{ + PS_ASSERT_PTR_NON_NULL(in, false); + PS_ASSERT_PTR_NON_NULL(out, false); + + int numChips = out->fpa->chips->n; // Number of cells + + psVector *meanMask = psVectorAlloc(numChips, PS_TYPE_VECTOR_MASK); // Mask for means + psVectorInit(meanMask, 0); + + // Mask bits + enum { + PM_PATTERN_IGNORE = 0x01, // Ignore this cell + PM_PATTERN_TWEAK = 0x02, // Tweak this cell + PM_PATTERN_ERROR = 0x04, // Error in calculating background + PM_PATTERN_ALL = 0xFF, // All causes + }; + + // Measure mean of each cell edge, and use that to determine the cell offsets. + psStatsOptions stat = cellStat; // Define which statistic to use. + + psStats *bgStats = psStatsAlloc(stat); // Statistics on background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + + psRegion region = {0,0,0,0}; + + /* These images hold the edge data for the OTA structure. */ + psImage *A = psImageAlloc(8,8,PS_TYPE_F64); // Top edge + psImage *B = psImageAlloc(8,8,PS_TYPE_F64); // Bottom edge + psImage *C = psImageAlloc(8,8,PS_TYPE_F64); // Right edge + psImage *D = psImageAlloc(8,8,PS_TYPE_F64); // Left edge + psImageInit(A,0.0); + psImageInit(B,0.0); + psImageInit(C,0.0); + psImageInit(D,0.0); + + // Corners don't exist. + A->data.F64[0][0] = NAN; + B->data.F64[0][0] = NAN; + C->data.F64[0][0] = NAN; + D->data.F64[0][0] = NAN; + A->data.F64[7][0] = NAN; + B->data.F64[7][0] = NAN; + C->data.F64[7][0] = NAN; + D->data.F64[7][0] = NAN; + A->data.F64[0][7] = NAN; + B->data.F64[0][7] = NAN; + C->data.F64[0][7] = NAN; + D->data.F64[0][7] = NAN; + A->data.F64[7][7] = NAN; + B->data.F64[7][7] = NAN; + C->data.F64[7][7] = NAN; + D->data.F64[7][7] = NAN; + + for (int i = 0; i < numChips; i++) { + pmChip *chip = out->fpa->chips->data[i]; + pmCell *cell = chip->cells->data[0]; // Cell of interest + pmReadout *ro = cell->readouts->data[0]; // Readout of interest + + psStatsInit(bgStats); + + // Convert cell iterator i into an xy coordinate on the grid of cells + // This is wrong for chips + const char *chipName = psMetadataLookupStr(NULL,chip->concepts, "CHIP.NAME"); + int xParity = psMetadataLookupS16(NULL,chip->concepts, "CHIP.XPARITY"); + int yParity = psMetadataLookupS16(NULL,chip->concepts, "CHIP.YPARITY"); + int x = chipName[2] - '0'; + int y = chipName[3] - '0'; + + + for (int j = 0; j < 4; j++) { + if (j == 0) { // Region B + region = psRegionSet(0,ro->image->numCols, + 0,edgeWidth); + } + else if (j == 1) { // Region A + region = psRegionSet(0,ro->image->numCols, + ro->image->numRows - edgeWidth,ro->image->numRows); + } + else if (j == 2) { // Region D + region = psRegionSet(0,edgeWidth, + 0,ro->image->numRows); + } + else if (j == 3) { // Region C + region = psRegionSet(ro->image->numCols - edgeWidth,ro->image->numCols, + 0,ro->image->numRows); + } + psImage *subset = psImageSubset(ro->image,region); + + if (!psImageBackground(bgStats, NULL, subset, NULL, maskVal, rng)) { + psWarning("Unable to measure background for cell %d on edge %d\n", i, j); + psErrorClear(); + meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PM_PATTERN_ERROR; + if (j == 0) { B->data.F64[y][x] = NAN; } + else if (j == 1) { A->data.F64[y][x] = NAN; } + else if (j == 2) { C->data.F64[y][x] = NAN; } + else if (j == 3) { D->data.F64[y][x] = NAN; } + psFree(subset); + continue; // Move on to next edge, as only part of this cell may be a problem + } + + // If the returned value is zero, assume something is wrong. Do I still need this? + if (psStatsGetValue(bgStats,stat) < 1e-6) { + if (j == 0) { B->data.F64[y][x] = NAN; } + else if (j == 1) { A->data.F64[y][x] = NAN; } + else if (j == 2) { C->data.F64[y][x] = NAN; } + else if (j == 3) { D->data.F64[y][x] = NAN; } + } + // If we have an error for this cell/edge, make sure we mask the value + if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_ERROR) { + if (j == 0) { B->data.F64[y][x] = NAN; } + else if (j == 1) { A->data.F64[y][x] = NAN; } + else if (j == 2) { C->data.F64[y][x] = NAN; } + else if (j == 3) { D->data.F64[y][x] = NAN; } + } + else { // Set the value to match what we got from the edge box. + if (xParity == -1) { + if (j == 2) { C->data.F64[y][x] = psStatsGetValue(bgStats,stat); } + else if (j == 3) { D->data.F64[y][x] = psStatsGetValue(bgStats,stat); } + } + if (xParity == 1) { + if (j == 3) { C->data.F64[y][x] = psStatsGetValue(bgStats,stat); } + else if (j == 2) { D->data.F64[y][x] = psStatsGetValue(bgStats,stat); } + } + if (yParity == 1) { + if (j == 0) { B->data.F64[y][x] = psStatsGetValue(bgStats,stat); } + else if (j == 1) { A->data.F64[y][x] = psStatsGetValue(bgStats,stat); } + } + if (yParity == -1) { + if (j == 1) { B->data.F64[y][x] = psStatsGetValue(bgStats,stat); } + else if (j == 0) { A->data.F64[y][x] = psStatsGetValue(bgStats,stat); } + } + + } +#if (0) + for (int u = 0; u < subset->numCols; u++) { + for (int v = 0; v < subset->numRows; v++) { + psTrace("psModules.detrend.cont",10,"BOX: %d %d (%d %d) (%d %d) %f %d", + i,j,x,y,u,v,subset->data.F32[v][u],0); + } + } +#endif + psFree(subset); + + } + psTrace("psModules.detrend.cont",5, "OTA: %d (%d %d) (%d %d) A: %f B: %f C: %f D: %f", + i,x,y,xParity,yParity, + A->data.F64[y][x],B->data.F64[y][x],C->data.F64[y][x],D->data.F64[y][x]); + } + psFree(bgStats); + psFree(rng); + + // We've now allocated all the edge values, so we can now minimize the offsets. + // This involves solving the equation A x = b, where + // A is the (64x64 for GPC1) matrix containing the edges that match for each cell + // x is the solution vector + // b is the combination of offsets across each cell boundary for each cell. + // Below "XX" is used as the matrix A, and "solution" is used as both b and x + // (due to the way psMatrixLUSolve operates). + psVector *solution = psVectorAlloc(64,PS_TYPE_F64); + psImage *XX = psImageAlloc(64,64,PS_TYPE_F64); + psVectorInit(solution,0.0); + psImageInit(XX,0.0); + + + for (int i = 0; i < numChips; i++) { + // Accumulate all the possible edge differences we can for this cell. + // As we do so, make a note of the correlations by incrementing the element of the matrix. + pmChip *chip = out->fpa->chips->data[i]; + const char *chipName = psMetadataLookupStr(NULL,chip->concepts, "CHIP.NAME"); + int x = chipName[2] - '0'; + int y = chipName[3] - '0'; + + int j; + int k = 8 * x + y; + double critical_value = 0.0; + + if (x + 1 < 8) { // We have a neighbor adjacent in the +x direction + j = 8 * (x + 1) + y; // Determine that neighbor's index + + if (fabs(C->data.F64[y][x]) > fabs(D->data.F64[y][x+1])) { + critical_value = 2.0 * fabs(D->data.F64[y][x+1]); + } + else { + critical_value = 2.0 * fabs(C->data.F64[y][x]); + } + if (critical_value < 25) { critical_value = 25; } + psTrace("psModules.detrend.cont",5,"CmD %d %d %d %d %g %g %g", // diagnostic + i,x,y,j, + C->data.F64[y][x], + D->data.F64[y][x+1], + critical_value + ); + if (// If there are no errors with the neighbor, + (isfinite(C->data.F64[y][x]))&&(isfinite(D->data.F64[y][x+1]))&& // and all edges have valid values, + (fabs(C->data.F64[y][x] - D->data.F64[y][x+1]) < critical_value) // and there are no large discontinuities, + ) { + solution->data.F64[k] += C->data.F64[y][x] - D->data.F64[y][x+1]; // Take the difference + XX->data.F64[k][k] += 1; // increment our relation with ourself + XX->data.F64[k][j] += -1; // decrement our relation with the neighbor + } + } + if (x - 1 > -1) { // etc. + j = 8 * (x - 1) + y; + if (fabs(C->data.F64[y][x-1]) > fabs(D->data.F64[y][x])) { + critical_value = 2.0 * fabs(D->data.F64[y][x]); + } + else { + critical_value = 2.0 * fabs(C->data.F64[y][x-1]); + } + if (critical_value < 25) { critical_value = 25; } + psTrace("psModules.detrend.cont",5,"DmC %d %d %d %d %g %g %g", + i,x,y,j, + D->data.F64[y][x], + C->data.F64[y][x-1], + critical_value + ); + + if ( + (isfinite(D->data.F64[y][x]))&&(isfinite(C->data.F64[y][x-1]))&& + (fabs(D->data.F64[y][x] - C->data.F64[y][x-1]) < critical_value) + ) { + solution->data.F64[k] += D->data.F64[y][x] - C->data.F64[y][x-1]; + XX->data.F64[k][k] += 1; + XX->data.F64[k][j] += -1; + } + } + if (y + 1 < 8) { + j = 8 * x + (y + 1); + psTrace("psModules.detrend.cont",5,"AmB %d %d %d %d %g %g", + i,x,y,j, + A->data.F64[y][x], + B->data.F64[y+1][x] + ); + if (fabs(A->data.F64[y][x]) > fabs(B->data.F64[y+1][x])) { + critical_value = 2.0 * fabs(B->data.F64[y+1][x]); + } + else { + critical_value = 2.0 * fabs(A->data.F64[y][x]); + } + if (critical_value < 25) { critical_value = 25; } + if ( + (isfinite(A->data.F64[y][x]))&&(isfinite(B->data.F64[y+1][x]))&& + (fabs(A->data.F64[y][x] - B->data.F64[y+1][x]) < critical_value) + ) { + solution->data.F64[k] += A->data.F64[y][x] - B->data.F64[y+1][x]; + XX->data.F64[k][k] += 1; + XX->data.F64[k][j] += -1; + } + } + if (y - 1 > -1) { + j = 8 * x + (y - 1); + psTrace("psModules.detrend.cont",5,"BmA %d %d %d %d %g %g", + i,x,y,j, + B->data.F64[y][x], + A->data.F64[y-1][x] + ); + if (fabs(A->data.F64[y-1][x]) > fabs(B->data.F64[y][x])) { + critical_value = 2.0 * fabs(B->data.F64[y][x]); + } + else { + critical_value = 2.0 * fabs(A->data.F64[y-1][x]); + } + if (critical_value < 25) { critical_value = 25; } + if ( + (isfinite(B->data.F64[y][x]))&&(isfinite(A->data.F64[y-1][x]))&& + (fabs(B->data.F64[y][x] - A->data.F64[y-1][x]) < critical_value) + ) { + solution->data.F64[k] += B->data.F64[y][x] - A->data.F64[y-1][x]; + XX->data.F64[k][k] += 1; + XX->data.F64[k][j] += -1; + } + } + } + double max_XX = 0; + double solution_V = 0; + int i_peak = -1; + for (int i = 0; i < numChips + 4; i++) { // If any cells have no value of themself, set the matrix to 1.0. + if (XX->data.F64[i][i] == 0.0) { + XX->data.F64[i][i] = 1.0; + } + if (XX->data.F64[i][i] > max_XX) { + max_XX = XX->data.F64[i][i]; + solution_V = solution->data.F64[i]; + i_peak = i; + } + } + psTrace("psModules.detrend.cont",5,"fixed point: %d %g\n", + i_peak,solution_V); + +#if (1) + for (int i = 0; i < numChips + 4; i++) { // print matrix A + psTrace("psModules.detrend.cont",5,"A: %3d % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f", + i, + XX->data.F64[i][0], XX->data.F64[i][1], XX->data.F64[i][2], XX->data.F64[i][3], + XX->data.F64[i][4], XX->data.F64[i][5], XX->data.F64[i][6], XX->data.F64[i][7], + XX->data.F64[i][8], XX->data.F64[i][9], XX->data.F64[i][10], XX->data.F64[i][11], + XX->data.F64[i][12], XX->data.F64[i][13], XX->data.F64[i][14], XX->data.F64[i][15], + XX->data.F64[i][16], XX->data.F64[i][17], XX->data.F64[i][18], XX->data.F64[i][19], + XX->data.F64[i][20], XX->data.F64[i][21], XX->data.F64[i][22], XX->data.F64[i][23], + XX->data.F64[i][24], XX->data.F64[i][25], XX->data.F64[i][26], XX->data.F64[i][27], + XX->data.F64[i][28], XX->data.F64[i][29], XX->data.F64[i][30], XX->data.F64[i][31], + XX->data.F64[i][32], XX->data.F64[i][33], XX->data.F64[i][34], XX->data.F64[i][35], + XX->data.F64[i][36], XX->data.F64[i][37], XX->data.F64[i][38], XX->data.F64[i][39], + XX->data.F64[i][40], XX->data.F64[i][41], XX->data.F64[i][42], XX->data.F64[i][43], + XX->data.F64[i][44], XX->data.F64[i][45], XX->data.F64[i][46], XX->data.F64[i][47], + XX->data.F64[i][48], XX->data.F64[i][49], XX->data.F64[i][50], XX->data.F64[i][51], + XX->data.F64[i][52], XX->data.F64[i][53], XX->data.F64[i][54], XX->data.F64[i][55], + XX->data.F64[i][56], XX->data.F64[i][57], XX->data.F64[i][58], XX->data.F64[i][59], + XX->data.F64[i][60], XX->data.F64[i][61], XX->data.F64[i][62], XX->data.F64[i][63] + ); + } + + for (int i = 0; i < numChips + 4; i++) { // print vector b + psTrace("psModules.detrend.cont",5,"b: %d %f", + i, + solution->data.F64[i] + ); + } +#endif + + // Solve the Ax=b equation + // psMatrixLUSolve(XX,solution); + psMatrixGJSolve(XX,solution); +#if (1) + for (int i = 0; i < numChips + 4; i++) { // print vector b + psTrace("psModules.detrend.cont",5,"x: %d %f", + i, + solution->data.F64[i] + ); + } +#endif + + // Cleanup + psFree(XX); + psFree(A); + psFree(B); + psFree(C); + psFree(D); + + // Correct cells based on the offsets calculated, and store the result in the analysis metadata. + for (int i = 0; i < numChips; i++) { + pmChip *chip = out->fpa->chips->data[i]; + pmCell *cell = chip->cells->data[0]; // Cell of interest + pmReadout *ro = cell->readouts->data[0]; // Readout of interest + const char *chipName = psMetadataLookupStr(NULL,chip->concepts, "CHIP.NAME"); + int x = chipName[2] - '0'; + int y = chipName[3] - '0'; + int k = 8 * x + y; + float correction = solution->data.F64[k]; + psLogMsg("psModules.detrend", PS_LOG_DETAIL, "Correcting background of chip %s by %f", + chipName, correction); + psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(correction, PS_TYPE_F32)); + psMetadataAddF32(ro->analysis, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE, + "Pattern chip correction solution", correction); + } + + psFree(solution); + psFree(meanMask); + + return true; +} Index: trunk/psModules/src/detrend/pmPattern.h =================================================================== --- trunk/psModules/src/detrend/pmPattern.h (revision 35074) +++ trunk/psModules/src/detrend/pmPattern.h (revision 35081) @@ -14,4 +14,10 @@ #include #include +#include +#include +#include +#include +#include +#include /// @addtogroup detrend Detrend Creation and Application @@ -65,4 +71,12 @@ ); +bool pmPatternContinuityBackground(pmFPAfile *in, + pmFPAfile *out, + psStatsOptions bgStat, + psStatsOptions cellStat, + psImageMaskType maskVal, + psImageMaskType maskBad, + int edgeWidth); + /// Apply previously measured cell pattern correction bool pmPatternContinuityApply(pmReadout *ro, ///< Readout to correct Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 35074) +++ trunk/psModules/src/imcombine/pmStack.c (revision 35081) @@ -1385,4 +1385,5 @@ psVector *pixelData = psVectorAlloc(input->n,PS_TYPE_F32); + psVector *pixelMask = psVectorAlloc(input->n,PS_TYPE_VECTOR_MASK); psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); @@ -1392,10 +1393,25 @@ pmReadout *ro = stack->data[i]; psImage *image = ro->image; + pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; pixelData->data.F32[i] = image->data.F32[y][x]; + if (isfinite(image->data.F32[y][x])&& + (fabs(image->data.F32[y][x]) < 1e5)) { + pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; + } + else { + pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1; + } +#if (0) + if ((x == 59)&&(y > 40)&&(y < 50)) { + fprintf(stderr,"%d %d %d %d %g\n", + x,y,i,pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i],pixelData->data.F32[i]); + } } - if (!psVectorStats(stats,pixelData,NULL,NULL,0)) { +#endif + if (!psVectorStats(stats,pixelData,NULL,pixelMask,1)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate median"); psFree(stats); psFree(pixelData); + psFree(pixelMask); psFree(stack); return(false); @@ -1403,4 +1419,11 @@ combined->image->data.F32[y][x] = stats->robustMedian; combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; +#if (0) + if ((x == 59)&&(y > 40)&&(y < 50)) { + fprintf(stderr,"%d %d %d %d %g\n", + x,y,-1,combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x], + combined->image->data.F32[y][x]); + } +#endif } } @@ -1408,4 +1431,5 @@ psFree(stats); psFree(pixelData); + psFree(pixelMask); psFree(stack); return (true);