IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 35081


Ignore:
Timestamp:
Feb 1, 2013, 5:01:08 PM (13 years ago)
Author:
watersc1
Message:

Merging changes that implement FPA level background continuity.

Location:
trunk
Files:
1 added
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/ippScripts/scripts/camera_exp.pl

    r33576 r35081  
    4141
    4242my ( $exp_tag, $cam_id, $camera, $outroot, $dbname, $reduction, $dvodb, $verbose, $no_update,
    43      $no_op, $redirect, $save_temps, $run_state, $skip_binned, $skip_masks);
     43     $no_op, $redirect, $save_temps, $run_state, $skip_binned, $skip_masks, $bkg_only);
    4444GetOptions(
    4545    'exp_tag=s'         => \$exp_tag, # Exposure identifier
    … …  
    5353    'skip-binned'       => \$skip_binned, # override recipe - don't create binned images
    5454    'skip-refmask'      => \$skip_masks, # override recipe - don't create refmask
     55    'bkg-only'          => \$bkg_only,  # override recipe - only do background continuity
    5556    'verbose'           => \$verbose,   # Print to stdout
    5657    'no-update'         => \$no_update, # Update the database?
    … …  
    111112&my_die("Unrecognised PSASTRO recipe", $cam_id, $PS_EXIT_CONFIG_ERROR) unless defined $recipe_psastro;
    112113
     114my $bkg_recipe     = 'PPIMAGE_BKGCONT'; # Add to reduction?
    113115my $mdcParser = PS::IPP::Metadata::Config->new; # Parser for metadata config files
    114116
    … …  
    215217my ($list4File, $list4Name) = tempfile( "/tmp/$exp_tag.cm.$cam_id.b4.list.XXXX", UNLINK => !$save_temps ); # For astrometry
    216218
     219### Create temp file for background models
     220my ($list5File, $list5Name) = tempfile( "/tmp/$exp_tag.cm.$cam_id.b5.list.XXXX", UNLINK => !$save_temps ); # For background models
     221
    217222# XXX we perform astrometry iff photometry output exists
    218223my @outMasks;                   # Names of output masks
     224my @bkg_models;                 # Names of output background models
    219225foreach my $file (@$files) {
    220226    next if $file->{quality} != 0;
    … …  
    226232    my $chipObjects = $ipprc->filename("PSPHOT.OUTPUT", $file->{path_base}, $class_id);
    227233    my $chipMask   = $ipprc->filename("PPIMAGE.CHIP.MASK", $file->{path_base}, $class_id);
    228 
     234   
    229235    print $list1File ($ipprc->filename("PPIMAGE.BIN1", $file->{path_base}, $class_id) . "\n");
    230236    print $list2File ($ipprc->filename("PPIMAGE.BIN2", $file->{path_base}, $class_id) . "\n");
    231237    print $list3File ($chipObjects . "\n");
    232238    print $list4File ($chipMask . "\n");
     239    print $list5File ($ipprc->filename("PSPHOT.BACKMDL", $file->{path_base}, $class_id) . "\n");
    233240
    234241    push @outMasks, prepare_output("PSASTRO.OUTPUT.MASK", $outroot, $class_id, 1) if $produceMasks;
     242    push @bkg_models, prepare_output("PPIMAGE.BACKMDL", $outroot, $class_id, 1);
    235243}
    236244close $list1File;
    … …  
    238246close $list3File;
    239247close $list4File;
    240 
     248close $list5File;
    241249# Output products
    242250$ipprc->outroot_prepare($outroot);
    … …  
    291299    }
    292300
    293     {
     301    if (!$bkg_only) {
    294302        # run psastro on the chipObjects, producing fpaObjects
    295303        my $command;
    … …  
    353361            }
    354362        }
     363    }
     364    # Construct FPA continuity corrected background images
     365    {
     366        my $command;
     367        $command = "$ppImage";
     368        $command .= " -list $list5Name $outroot";
     369        $command .= " -recipe PPIMAGE $bkg_recipe";
     370        $command .= " -dbname $dbname" if defined $dbname;
     371
     372        my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
     373            run(command => $command, verbose => $verbose);
     374        unless ($success) {
     375            $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
     376            &my_die("Unable to perform ppImage to fix background: $error_code", $cam_id, $error_code);
     377        }
     378        foreach my $bkgModel (@bkg_models) {
     379            check_output($bkgModel, $replicateOutputs);
     380        }
    355381    }
    356382}

  • trunk/ippTools/share/stacktool_tobkg.sql

    r34800 r35081  
    2525    -- WHERE hook %s
    2626GROUP BY stack_id
    27 HAVING (SUM(IF(warpRun.state = 'full', 1, 0)) = COUNT(stackInputSkyfile.warp_id) AND
     27HAVING ((SUM(IF(warpRun.state = 'cleaned', 1, 0)) = COUNT(stackInputSkyfile.warp_id) OR
     28         SUM(IF(warpRun.state = 'full', 1, 0)) = COUNT(stackInputSkyfile.warp_id)) AND
    2829        SUM(IF(warpSkyfile.background_model = 1, 1, 0)) = COUNT(stackInputSkyfile.warp_id))
    2930

  • trunk/ippconfig/recipes/filerules-mef.mdc

    r34800 r35081  
    214214PPIMAGE.PATTERN         OUTPUT {OUTPUT}.ptn                      PATTERN   NONE       CHIP       TRUE      MEF
    215215PPIMAGE.STATS           OUTPUT {OUTPUT}.stats                    STATS     NONE       FPA        TRUE      MEF
    216 
     216PPIMAGE.BACKMDL              OUTPUT {OUTPUT}.{CHIP.NAME}.mdl.fits     IMAGE           NONE       CHIP       TRUE      NONE
    217217
    218218## note: these use the same output naming convention since they are used for different ppMerge runs

  • trunk/ippconfig/recipes/filerules-simple.mdc

    r34800 r35081  
    178178PPIMAGE.PATTERN              OUTPUT {OUTPUT}.ptn                  PATTERN         NONE       FPA        TRUE      NONE
    179179PPIMAGE.STATS                OUTPUT {OUTPUT}.stats                STATS           NONE       FPA        TRUE      NONE
    180 
     180PPIMAGE.BACKMDL              OUTPUT {OUTPUT}.{CHIP.NAME}.mdl.fits     IMAGE           NONE       CHIP       TRUE      NONE
    181181
    182182## note: these use the  same output naming convention since they are used for different ppMerge runs

  • trunk/ippconfig/recipes/filerules-split.mdc

    r34800 r35081  
    198198PPIMAGE.PATTERN              OUTPUT {OUTPUT}.{CHIP.NAME}.ptn          PATTERN         NONE       CHIP       TRUE      NONE
    199199PPIMAGE.STATS                OUTPUT {OUTPUT}.{CHIP.NAME}.stats        STATS           NONE       CHIP       TRUE      NONE
     200PPIMAGE.BACKMDL              OUTPUT {OUTPUT}.{CHIP.NAME}.mdl.fits     IMAGE           NONE       CHIP       TRUE      NONE
    200201
    201202## note: these use the  same output naming convention since they are used for different ppMerge runs

  • trunk/ippconfig/recipes/ppImage.config

    r34063 r35081  
    4242MASK.STATS         BOOL    FALSE           # calculate mask statistics.
    4343GAIN.OVERRIDE           BOOL    TRUE            # Override a non-finite gain?
    44 
     44BACKGROUND.CONTINUITY BOOL FALSE           # Apply continuity correction to background models for an FPA.
    4545TILTYSTREAK.BY.CLASS METADATA              # apply the 'tiltystreak' tool
    4646END
    … …  
    11471147  MASK.SATURATED   BOOL    FALSE           # DO NOT Mask the saturated pixels
    11481148  MASK.LOW         BOOL    FALSE           # DO NOT Mask the saturated pixels
     1149END
     1150
     1151# Calculate new background models with continuity correction applied
     1152PPIMAGE_BKGCONT    METADATA
     1153  OVERSCAN         BOOL    FALSE           # Overscan subtraction
     1154  NONLIN           BOOL    FALSE           # Non-linearity correction; not implemented
     1155  BIAS             BOOL    FALSE           # Bias subtraction
     1156  DARK             BOOL    FALSE           # Dark subtraction
     1157  SHUTTER          BOOL    FALSE           # Shutter correction
     1158  FLAT             BOOL    FALSE           # Flat-field normalisation
     1159  MASK             BOOL    FALSE           # Mask bad pixels
     1160  MASK.BUILD       BOOL    FALSE           # Build internal mask?
     1161  FRINGE           BOOL    FALSE           # Fringe subtraction
     1162  PHOTOM           BOOL    FALSE           # Source identification and photometry
     1163  ASTROM.CHIP      BOOL    FALSE           # Astrometry per chip?
     1164  ASTROM.MOSAIC    BOOL    FALSE           # Astrometry for mosaic?
     1165  BASE.FITS        BOOL    FALSE           # Save base image?
     1166  BASE.MASK.FITS   BOOL    FALSE           # Save base detrended image?
     1167  BASE.VARIANCE.FITS BOOL    FALSE           # Save base detrended image?
     1168  CHIP.FITS        BOOL    FALSE           # Save chip-mosaic-ed image?
     1169  CHIP.MASK.FITS   BOOL    FALSE           # Save chip-mosaic-ed image?
     1170  CHIP.VARIANCE.FITS BOOL    FALSE           # Save chip-mosaic-ed image?
     1171  BIN1.FITS        BOOL    FALSE           # Save 1st binned chip image?
     1172  BIN2.FITS        BOOL    FALSE           # Save 2nd binned chip image?
     1173  BIN1.JPEG        BOOL    FALSE           # Save 1st binned jpeg?
     1174  BIN2.JPEG        BOOL    FALSE           # Save 2nd binned jpeg?
     1175  BIN1.XBIN        S32     1               # Image is already binned
     1176  BIN1.YBIN        S32     1               # Image is already binned
     1177  BIN2.XBIN        S32     1               # Image is already binned
     1178  BIN2.YBIN        S32     1               # Image is already binned
     1179  MASK.SATURATED   BOOL    FALSE           # DO NOT Mask the saturated pixels
     1180  MASK.LOW         BOOL    FALSE           # DO NOT Mask the saturated pixels
     1181  BACKGROUND.CONTINUITY BOOL TRUE
     1182  PATTERN.CONTINUITY.WIDTH S32 2
    11491183END
    11501184

  • trunk/ppImage/src/Makefile.am

    r33590 r35081  
    4444        ppImageRebinReadout.c \
    4545        ppImageMosaic.c \
     46        ppImageMosaicBackground.c \
    4647        ppImageMaskStats.c \
    4748        ppImagePhotom.c \

  • trunk/ppImage/src/ppImage.h

    r33848 r35081  
    4141    bool doPatternCell;                 // Cell pattern correction
    4242    bool doPatternContinuity;           // Cell continuity correction
     43    bool doBackgroundContinuity;        // Do mosaic continuity correction
    4344    bool doFringe;                      // Fringe subtraction
    4445    bool doPhotom;                      // Source identification and photometry
    … …  
    5051    bool checkNoise;                    // measure cell-level variance
    5152    bool applyParity;                   // Apply Cell parities
    52 
    53   bool doMaskStats;                      // Calculate mask statistics
     53    bool doMaskStats;                   // Calculate mask statistics
    5454 
    5555    bool doCrosstalkMeasure;            // measure crosstalk signal
    … …  
    183183bool ppImageDetrendPatternApply(pmConfig *config, pmChip *chip, const pmFPAview *inputView, const ppImageOptions *options);
    184184
     185// Do background continuity step
     186bool ppImageMosaicBackground(pmConfig *config, const ppImageOptions *options);
     187
    185188// Record which detrend file was used for the detrending
    186189bool ppImageDetrendRecord(

  • trunk/ppImage/src/ppImageLoop.c

    r33590 r35081  
    4141        ESCAPE("load failure for FPA");
    4242    }
    43 
     43   
    4444    pmChip *chip;                       // Chip from FPA
    4545    while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) {
    … …  
    233233            ESCAPE("Unable to bin chip (level 2).");
    234234        }
    235 
     235        if (options->doBackgroundContinuity) {
     236          pmFPAfile *out = psMetadataLookupPtr(NULL, config->files, "PPIMAGE.BACKMDL");
     237          pmFPAfileCopyView(out->fpa,out->src,view);
     238        }
    236239        // Close cells (XXX shouldn't pmFPAfileClose iterate down as needed?)
    237240        view->cell = -1;
    … …  
    251254    }
    252255
     256    // Do background model continuity updates
     257    if (!ppImageMosaicBackground(config, options )) {
     258      ESCAPE("failure in Background Mosaic");
     259    }
     260   
    253261    // generate the full-scale FPA mosaic
    254262    if (!ppImageMosaicFPA(config, options, "PPIMAGE.OUTPUT.FPA1", "PPIMAGE.BIN1")) {

  • trunk/ppImage/src/ppImageOptions.c

    r33848 r35081  
    3535    options->doPatternCell   = false;   // Cell pattern correction
    3636    options->doPatternContinuity = false; // Cell continuity correction
     37    options->doBackgroundContinuity = false; // Chip level background continuity correction
    3738    options->doFringe        = false;   // Fringe subtraction
    3839    options->doPhotom        = false;   // Source identification and photometry
    … …  
    429430    }
    430431
     432    // Option to enable the background continuity
     433    options->doBackgroundContinuity = psMetadataLookupBool(NULL, recipe, "BACKGROUND.CONTINUITY");
     434
    431435
    432436    // Remnance options

  • trunk/ppImage/src/ppImageParseCamera.c

    r33848 r35081  
    303303    }
    304304
     305
    305306    // chipImage    -> psphotInput  (pmFPAfileDefineFromFile)       : fpa is constructed
    306307    // psphotInput  -> psphotOutput (pmFPAfileDefineOutputFromFile) : fpa is equated
    … …  
    387388    // the input data is the same as the outImage data : force the free levels to match
    388389    input->freeLevel = PS_MIN(outImage->freeLevel, input->freeLevel);
     390
     391    // define the continuity corrected background model files
     392    if (options->doBackgroundContinuity) {
     393      pmFPAfile *bkgMosaicModel = pmFPAfileDefineFromFPA(config,input->fpa, 1, 1,  "PPIMAGE.BACKMDL");
     394      if (!bkgMosaicModel) {
     395        psError(PS_ERR_IO, false, _("Unable to generate new file from PPIMAGE.BACKMDL"));
     396        psFree(options);
     397        return NULL;
     398      }
     399      if (bkgMosaicModel->type != PM_FPA_FILE_IMAGE) {
     400        psError(PS_ERR_IO, true, "PPIMAGE.BACKMDL is not of type IMAGE");
     401        psFree(options);
     402        return NULL;
     403
     404      }
     405      bkgMosaicModel->save = options->doBackgroundContinuity;
     406      //      bkgMosaicModel->freeLevel = PS_MIN(bkgMosaicModel->freeLevel, PM_FPA_LEVEL_FPA);
     407      bkgMosaicModel->freeLevel = PM_FPA_LEVEL_FPA;
     408      bkgMosaicModel->dataLevel = bkgMosaicModel->dataLevel;
     409      bkgMosaicModel->fileLevel = PS_MIN(bkgMosaicModel->fileLevel, bkgMosaicModel->dataLevel);
     410    }
    389411
    390412    // define the binned target files (which may just be carriers for some camera configurations)
    … …  
    497519        chipVariance->save = false;
    498520    }
    499 
     521   
    500522    if (psTraceGetLevel("ppImage.config") > 0) {
    501523        // Get a look inside all the files.

  • trunk/ppStack/src/ppStackCamera.c

    r34800 r35081  
    4747    }
    4848    bool convolve = psMetadataLookupBool(NULL, recipe, "CONVOLVE"); // Convolve images before stack?
    49 
     49    bool doBackground = psMetadataLookupBool(NULL, recipe, "BACKGROUND.MODEL");
     50   
    5051    psArray *runImages = pmFPAfileDefineMultipleFromRun(&status, NULL, config,
    5152                                                        "PPSTACK.INPUT"); // Input images from previous run
    … …  
    217218
    218219            // Grab bkgmodel information here
    219             if (!bkgmodel ||  strlen(bkgmodel) == 0) {
    220               // We have no background models.
     220            if (doBackground) {
     221              if (!bkgmodel ||  strlen(bkgmodel) == 0) {
     222                // We have no background models.
     223              }
     224              else {
     225                pmFPAfile *inputBKG = defineFile(config,NULL,"PPSTACK.INPUT.BKGMODEL",bkgmodel,
     226                                                 PM_FPA_FILE_IMAGE);
     227                if (!inputBKG) {
     228                  // We failed to generate an pmFPAfile, so disable the background construction and continue.
     229                  doBackground = false;
     230                  psMetadataAddBool(recipe, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE,
     231                                    "Do photometry on stacked image?", false);
     232
     233#if (0)
     234                  psError(psErrorCodeLast(), false,
     235                          "Unable to define file from bkgmodel %d (%s)",i,bkgmodel);
     236                  return(false);
     237#endif
     238                }
     239              }// End bkgmodel
    221240            }
    222             else {
    223               pmFPAfile *inputBKG = defineFile(config,NULL,"PPSTACK.INPUT.BKGMODEL",bkgmodel,
    224                                                PM_FPA_FILE_IMAGE);
    225               if (!inputBKG) {
    226                 psError(psErrorCodeLast(), false,
    227                         "Unable to define file from bkgmodel %d (%s)",i,bkgmodel);
    228                 return(false);
    229               }
    230             }// End bkgmodel
    231241
    232242

  • trunk/ppStack/src/ppStackCombineAlternate.c

    r34848 r35081  
    198198      offset = (Sy * Sxx - Sx * Sxy) / D;
    199199      scale  = (S * Sxy - Sx * Sy) / D;
    200       fprintf(stderr,"Scales: %d %g %g\n",i,offset,scale);
     200      fprintf(stderr,"Scales: %d %g %g %g %g %g %g %g %g\n",i,offset,scale,D,Sx,Sy,Sxx,Syy,Sxy);
    201201      // Apply scaling factors
    202202      for (int y = minInputRows; y < maxInputRows; y++) {

  • trunk/psModules/src/detrend/pmPattern.c

    r33340 r35081  
    938938}
    939939
    940 
     940// This should reuse code more efficiently
     941bool pmPatternContinuityBackground(pmFPAfile *in, pmFPAfile *out,  psStatsOptions bgStat, psStatsOptions cellStat,
     942                                   psImageMaskType maskVal, psImageMaskType maskBad, int edgeWidth)
     943{
     944    PS_ASSERT_PTR_NON_NULL(in, false);
     945    PS_ASSERT_PTR_NON_NULL(out, false);
     946
     947    int numChips = out->fpa->chips->n;            // Number of cells
     948
     949    psVector *meanMask = psVectorAlloc(numChips, PS_TYPE_VECTOR_MASK); // Mask for means
     950    psVectorInit(meanMask, 0);
     951
     952    // Mask bits
     953    enum {
     954        PM_PATTERN_IGNORE = 0x01,       // Ignore this cell
     955        PM_PATTERN_TWEAK  = 0x02,       // Tweak this cell
     956        PM_PATTERN_ERROR  = 0x04,       // Error in calculating background
     957        PM_PATTERN_ALL    = 0xFF,       // All causes
     958    };
     959
     960    // Measure mean of each cell edge, and use that to determine the cell offsets.
     961    psStatsOptions stat = cellStat;     // Define which statistic to use.
     962   
     963    psStats *bgStats = psStatsAlloc(stat); // Statistics on background
     964    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     965
     966    psRegion region = {0,0,0,0};
     967
     968    /* These images hold the edge data for the OTA structure.  */
     969    psImage *A = psImageAlloc(8,8,PS_TYPE_F64); // Top edge
     970    psImage *B = psImageAlloc(8,8,PS_TYPE_F64); // Bottom edge
     971    psImage *C = psImageAlloc(8,8,PS_TYPE_F64); // Right edge
     972    psImage *D = psImageAlloc(8,8,PS_TYPE_F64); // Left edge
     973    psImageInit(A,0.0);
     974    psImageInit(B,0.0);
     975    psImageInit(C,0.0);
     976    psImageInit(D,0.0);
     977
     978    // Corners don't exist.
     979    A->data.F64[0][0] = NAN;
     980    B->data.F64[0][0] = NAN;
     981    C->data.F64[0][0] = NAN;
     982    D->data.F64[0][0] = NAN;
     983    A->data.F64[7][0] = NAN;
     984    B->data.F64[7][0] = NAN;
     985    C->data.F64[7][0] = NAN;
     986    D->data.F64[7][0] = NAN;
     987    A->data.F64[0][7] = NAN;
     988    B->data.F64[0][7] = NAN;
     989    C->data.F64[0][7] = NAN;
     990    D->data.F64[0][7] = NAN;
     991    A->data.F64[7][7] = NAN;
     992    B->data.F64[7][7] = NAN;
     993    C->data.F64[7][7] = NAN;
     994    D->data.F64[7][7] = NAN;
     995   
     996    for (int i = 0; i < numChips; i++) {
     997      pmChip *chip = out->fpa->chips->data[i];
     998      pmCell *cell = chip->cells->data[0]; // Cell of interest
     999      pmReadout *ro = cell->readouts->data[0]; // Readout of interest
     1000
     1001      psStatsInit(bgStats);
     1002
     1003      // Convert cell iterator i into an xy coordinate on the grid of cells
     1004      // This is wrong for chips
     1005      const char *chipName = psMetadataLookupStr(NULL,chip->concepts, "CHIP.NAME");
     1006      int xParity          = psMetadataLookupS16(NULL,chip->concepts, "CHIP.XPARITY");
     1007      int yParity          = psMetadataLookupS16(NULL,chip->concepts, "CHIP.YPARITY");
     1008      int x = chipName[2] - '0';
     1009      int y = chipName[3] - '0';
     1010
     1011     
     1012      for (int j = 0; j < 4; j++) {
     1013        if (j == 0) {  // Region B
     1014          region = psRegionSet(0,ro->image->numCols,
     1015                               0,edgeWidth);
     1016        }
     1017        else if (j == 1) { // Region A
     1018          region = psRegionSet(0,ro->image->numCols,
     1019                               ro->image->numRows - edgeWidth,ro->image->numRows);
     1020        }
     1021        else if (j == 2) { // Region D
     1022          region = psRegionSet(0,edgeWidth,
     1023                               0,ro->image->numRows);
     1024        }
     1025        else if (j == 3) { // Region C
     1026          region = psRegionSet(ro->image->numCols - edgeWidth,ro->image->numCols,
     1027                               0,ro->image->numRows);
     1028        }
     1029        psImage *subset  = psImageSubset(ro->image,region);
     1030
     1031        if (!psImageBackground(bgStats, NULL, subset, NULL, maskVal, rng)) {
     1032          psWarning("Unable to measure background for cell %d on edge %d\n", i, j);
     1033          psErrorClear();
     1034          meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PM_PATTERN_ERROR;
     1035          if (j == 0)      { B->data.F64[y][x] = NAN; }
     1036          else if (j == 1) { A->data.F64[y][x] = NAN; }
     1037          else if (j == 2) { C->data.F64[y][x] = NAN; }
     1038          else if (j == 3) { D->data.F64[y][x] = NAN; }
     1039          psFree(subset);
     1040          continue; // Move on to next edge, as only part of this cell may be a problem
     1041        }
     1042       
     1043        // If the returned value is zero, assume something is wrong.  Do I still need this?
     1044        if (psStatsGetValue(bgStats,stat) < 1e-6) {
     1045          if (j == 0)      { B->data.F64[y][x] = NAN; }
     1046          else if (j == 1) { A->data.F64[y][x] = NAN; }
     1047          else if (j == 2) { C->data.F64[y][x] = NAN; }
     1048          else if (j == 3) { D->data.F64[y][x] = NAN; }
     1049        }
     1050        // If we have an error for this cell/edge, make sure we mask the value
     1051        if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_ERROR) {
     1052          if (j == 0)      { B->data.F64[y][x] = NAN; }
     1053          else if (j == 1) { A->data.F64[y][x] = NAN; }
     1054          else if (j == 2) { C->data.F64[y][x] = NAN; }
     1055          else if (j == 3) { D->data.F64[y][x] = NAN; }
     1056        }
     1057        else { // Set the value to match what we got from the edge box.
     1058          if (xParity == -1) {
     1059            if (j == 2)      { C->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1060            else if (j == 3) { D->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1061          }
     1062          if (xParity == 1) {
     1063            if (j == 3)      { C->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1064            else if (j == 2) { D->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1065          }
     1066          if (yParity == 1) {
     1067            if (j == 0)      { B->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1068            else if (j == 1) { A->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1069          }
     1070          if (yParity == -1) {
     1071            if (j == 1)      { B->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1072            else if (j == 0) { A->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1073          }
     1074         
     1075        }
     1076#if (0)
     1077        for (int u = 0; u < subset->numCols; u++) {
     1078          for (int v = 0; v < subset->numRows; v++) {
     1079            psTrace("psModules.detrend.cont",10,"BOX: %d %d (%d %d) (%d %d) %f %d",
     1080                    i,j,x,y,u,v,subset->data.F32[v][u],0);
     1081          }
     1082        }
     1083#endif 
     1084        psFree(subset);
     1085       
     1086      }
     1087      psTrace("psModules.detrend.cont",5, "OTA: %d (%d %d) (%d %d) A: %f B: %f C: %f D: %f",
     1088              i,x,y,xParity,yParity,
     1089              A->data.F64[y][x],B->data.F64[y][x],C->data.F64[y][x],D->data.F64[y][x]);         
     1090    }
     1091    psFree(bgStats);
     1092    psFree(rng);
     1093   
     1094    // We've now allocated all the edge values, so we can now minimize the offsets.
     1095    // This involves solving the equation A x = b, where
     1096    // A is the (64x64 for GPC1) matrix containing the edges that match for each cell
     1097    // x is the solution vector
     1098    // b is the combination of offsets across each cell boundary for each cell.
     1099    // Below "XX" is used as the matrix A, and "solution" is used as both b and x
     1100    //   (due to the way psMatrixLUSolve operates).
     1101    psVector *solution = psVectorAlloc(64,PS_TYPE_F64);
     1102    psImage  *XX       = psImageAlloc(64,64,PS_TYPE_F64);
     1103    psVectorInit(solution,0.0);
     1104    psImageInit(XX,0.0);
     1105
     1106
     1107    for (int i = 0; i < numChips; i++) {
     1108      // Accumulate all the possible edge differences we can for this cell.
     1109      // As we do so, make a note of the correlations by incrementing the element of the matrix.
     1110      pmChip *chip = out->fpa->chips->data[i];
     1111      const char *chipName = psMetadataLookupStr(NULL,chip->concepts, "CHIP.NAME");
     1112      int x = chipName[2] - '0';
     1113      int y = chipName[3] - '0';
     1114
     1115      int j;
     1116      int k = 8 * x + y;
     1117      double critical_value = 0.0;
     1118
     1119      if (x + 1 < 8) {  // We have a neighbor adjacent in the +x direction
     1120        j = 8 * (x + 1) + y; // Determine that neighbor's index
     1121       
     1122        if (fabs(C->data.F64[y][x]) > fabs(D->data.F64[y][x+1])) {
     1123          critical_value = 2.0 * fabs(D->data.F64[y][x+1]);
     1124        }
     1125        else {
     1126          critical_value = 2.0 * fabs(C->data.F64[y][x]);
     1127        }
     1128        if (critical_value < 25) { critical_value = 25; }
     1129        psTrace("psModules.detrend.cont",5,"CmD %d %d %d %d %g %g %g", // diagnostic
     1130                i,x,y,j,
     1131                C->data.F64[y][x],
     1132                D->data.F64[y][x+1],
     1133                critical_value
     1134                );
     1135        if (// If there are no errors with the neighbor,
     1136            (isfinite(C->data.F64[y][x]))&&(isfinite(D->data.F64[y][x+1]))&&     // and all edges have valid values,
     1137            (fabs(C->data.F64[y][x] - D->data.F64[y][x+1]) < critical_value)     // and there are no large discontinuities,
     1138            ) {   
     1139          solution->data.F64[k] += C->data.F64[y][x] - D->data.F64[y][x+1];     // Take the difference
     1140          XX->data.F64[k][k] += 1;                                              // increment our relation with ourself
     1141          XX->data.F64[k][j] += -1;                                             // decrement our relation with the neighbor
     1142        }
     1143      }
     1144      if (x - 1 > -1) { // etc.
     1145        j = 8 * (x - 1) + y;
     1146        if (fabs(C->data.F64[y][x-1]) > fabs(D->data.F64[y][x])) {
     1147          critical_value = 2.0 * fabs(D->data.F64[y][x]);
     1148        }
     1149        else {
     1150          critical_value = 2.0 * fabs(C->data.F64[y][x-1]);
     1151        }
     1152        if (critical_value < 25) { critical_value = 25; }
     1153        psTrace("psModules.detrend.cont",5,"DmC %d %d %d %d %g %g %g",
     1154                i,x,y,j,
     1155                D->data.F64[y][x],
     1156                C->data.F64[y][x-1],
     1157                critical_value
     1158                );
     1159
     1160        if (
     1161            (isfinite(D->data.F64[y][x]))&&(isfinite(C->data.F64[y][x-1]))&&
     1162            (fabs(D->data.F64[y][x] - C->data.F64[y][x-1]) < critical_value)
     1163            ) {
     1164          solution->data.F64[k] += D->data.F64[y][x] - C->data.F64[y][x-1];
     1165          XX->data.F64[k][k] += 1;
     1166          XX->data.F64[k][j] += -1;
     1167        }
     1168      }
     1169      if (y + 1 < 8) {
     1170        j = 8 * x + (y + 1);
     1171        psTrace("psModules.detrend.cont",5,"AmB %d %d %d %d %g %g",
     1172                i,x,y,j,
     1173                A->data.F64[y][x],
     1174                B->data.F64[y+1][x]
     1175                );
     1176        if (fabs(A->data.F64[y][x]) > fabs(B->data.F64[y+1][x])) {
     1177          critical_value = 2.0 * fabs(B->data.F64[y+1][x]);
     1178        }
     1179        else {
     1180          critical_value = 2.0 * fabs(A->data.F64[y][x]);
     1181        }
     1182        if (critical_value < 25) { critical_value = 25; }
     1183        if (
     1184            (isfinite(A->data.F64[y][x]))&&(isfinite(B->data.F64[y+1][x]))&&
     1185            (fabs(A->data.F64[y][x] - B->data.F64[y+1][x]) < critical_value)
     1186            ) {
     1187          solution->data.F64[k] += A->data.F64[y][x] - B->data.F64[y+1][x];
     1188          XX->data.F64[k][k] += 1;
     1189          XX->data.F64[k][j] += -1;
     1190        }
     1191      }
     1192      if (y - 1 > -1) {
     1193        j = 8 * x +  (y - 1);
     1194        psTrace("psModules.detrend.cont",5,"BmA %d %d %d %d %g %g",
     1195                i,x,y,j,
     1196                B->data.F64[y][x],
     1197                A->data.F64[y-1][x]
     1198                );
     1199        if (fabs(A->data.F64[y-1][x]) > fabs(B->data.F64[y][x])) {
     1200          critical_value = 2.0 * fabs(B->data.F64[y][x]);
     1201        }
     1202        else {
     1203          critical_value = 2.0 * fabs(A->data.F64[y-1][x]);
     1204        }
     1205        if (critical_value < 25) { critical_value = 25; }
     1206        if (
     1207            (isfinite(B->data.F64[y][x]))&&(isfinite(A->data.F64[y-1][x]))&&
     1208            (fabs(B->data.F64[y][x] - A->data.F64[y-1][x]) < critical_value)
     1209            ) {
     1210          solution->data.F64[k] += B->data.F64[y][x] - A->data.F64[y-1][x];
     1211          XX->data.F64[k][k] += 1;
     1212          XX->data.F64[k][j] += -1;
     1213        }
     1214      }
     1215    }
     1216    double max_XX = 0;
     1217    double solution_V = 0;
     1218    int i_peak = -1;
     1219    for (int i = 0; i < numChips + 4; i++) { // If any cells have no value of themself, set the matrix to 1.0.
     1220      if (XX->data.F64[i][i] == 0.0) {
     1221        XX->data.F64[i][i] = 1.0;
     1222      }
     1223      if (XX->data.F64[i][i] > max_XX) {
     1224        max_XX = XX->data.F64[i][i];
     1225        solution_V = solution->data.F64[i];
     1226        i_peak = i;
     1227      }
     1228    }
     1229    psTrace("psModules.detrend.cont",5,"fixed point: %d %g\n",
     1230            i_peak,solution_V);
     1231   
     1232#if (1)
     1233    for (int i = 0; i < numChips + 4; i++) { // print matrix A
     1234      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",
     1235              i,
     1236              XX->data.F64[i][0],             XX->data.F64[i][1],             XX->data.F64[i][2],             XX->data.F64[i][3],
     1237              XX->data.F64[i][4],             XX->data.F64[i][5],             XX->data.F64[i][6],             XX->data.F64[i][7],
     1238              XX->data.F64[i][8],             XX->data.F64[i][9],             XX->data.F64[i][10],            XX->data.F64[i][11],
     1239              XX->data.F64[i][12],            XX->data.F64[i][13],            XX->data.F64[i][14],            XX->data.F64[i][15],
     1240              XX->data.F64[i][16],            XX->data.F64[i][17],            XX->data.F64[i][18],            XX->data.F64[i][19],
     1241              XX->data.F64[i][20],            XX->data.F64[i][21],            XX->data.F64[i][22],            XX->data.F64[i][23],
     1242              XX->data.F64[i][24],            XX->data.F64[i][25],            XX->data.F64[i][26],            XX->data.F64[i][27],
     1243              XX->data.F64[i][28],            XX->data.F64[i][29],            XX->data.F64[i][30],            XX->data.F64[i][31],
     1244              XX->data.F64[i][32],            XX->data.F64[i][33],            XX->data.F64[i][34],            XX->data.F64[i][35],
     1245              XX->data.F64[i][36],            XX->data.F64[i][37],            XX->data.F64[i][38],            XX->data.F64[i][39],
     1246              XX->data.F64[i][40],            XX->data.F64[i][41],            XX->data.F64[i][42],            XX->data.F64[i][43],
     1247              XX->data.F64[i][44],            XX->data.F64[i][45],            XX->data.F64[i][46],            XX->data.F64[i][47],
     1248              XX->data.F64[i][48],            XX->data.F64[i][49],            XX->data.F64[i][50],            XX->data.F64[i][51],
     1249              XX->data.F64[i][52],            XX->data.F64[i][53],            XX->data.F64[i][54],            XX->data.F64[i][55],
     1250              XX->data.F64[i][56],            XX->data.F64[i][57],            XX->data.F64[i][58],            XX->data.F64[i][59],
     1251              XX->data.F64[i][60],            XX->data.F64[i][61],            XX->data.F64[i][62],            XX->data.F64[i][63]
     1252              );
     1253    }
     1254
     1255    for (int i = 0; i < numChips + 4; i++) { // print vector b
     1256      psTrace("psModules.detrend.cont",5,"b: %d %f",
     1257              i,
     1258              solution->data.F64[i]
     1259              );
     1260    }
     1261#endif   
     1262   
     1263    // Solve the Ax=b equation
     1264    //    psMatrixLUSolve(XX,solution);
     1265    psMatrixGJSolve(XX,solution);
     1266#if (1)
     1267    for (int i = 0; i < numChips + 4; i++) { // print vector b
     1268      psTrace("psModules.detrend.cont",5,"x: %d %f",
     1269              i,
     1270              solution->data.F64[i]
     1271              );
     1272    }
     1273#endif
     1274   
     1275    // Cleanup
     1276    psFree(XX);
     1277    psFree(A);
     1278    psFree(B);
     1279    psFree(C);
     1280    psFree(D);
     1281
     1282    // Correct cells based on the offsets calculated, and store the result in the analysis metadata.
     1283    for (int i = 0; i < numChips; i++) {
     1284        pmChip *chip = out->fpa->chips->data[i];
     1285        pmCell *cell = chip->cells->data[0]; // Cell of interest
     1286        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
     1287        const char *chipName = psMetadataLookupStr(NULL,chip->concepts, "CHIP.NAME");
     1288        int x = chipName[2] - '0';
     1289        int y = chipName[3] - '0';
     1290        int k = 8 * x + y;
     1291        float correction = solution->data.F64[k];
     1292        psLogMsg("psModules.detrend", PS_LOG_DETAIL, "Correcting background of chip %s by %f",
     1293                 chipName, correction);
     1294        psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(correction, PS_TYPE_F32));
     1295        psMetadataAddF32(ro->analysis, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE,
     1296                         "Pattern chip correction solution", correction);
     1297    }
     1298
     1299    psFree(solution);
     1300    psFree(meanMask);
     1301
     1302    return true;
     1303}

  • trunk/psModules/src/detrend/pmPattern.h

    r33243 r35081  
    1414#include <pslib.h>
    1515#include <pmFPA.h>
     16#include <pmFPAview.h>
     17#include <pmFPAfile.h>
     18#include <pmFPAfileFitsIO.h>
     19#include <pmFPAHeader.h>
     20#include <pmHDU.h>
     21#include <pmHDUUtils.h>
    1622
    1723/// @addtogroup detrend Detrend Creation and Application
    … …  
    6571    );
    6672
     73bool pmPatternContinuityBackground(pmFPAfile *in,
     74                                   pmFPAfile *out,
     75                                   psStatsOptions bgStat,
     76                                   psStatsOptions cellStat,
     77                                   psImageMaskType maskVal,
     78                                   psImageMaskType maskBad,
     79                                   int edgeWidth);
     80
    6781/// Apply previously measured cell pattern correction
    6882bool pmPatternContinuityApply(pmReadout *ro,          ///< Readout to correct

  • trunk/psModules/src/imcombine/pmStack.c

    r34842 r35081  
    13851385
    13861386  psVector *pixelData = psVectorAlloc(input->n,PS_TYPE_F32);
     1387  psVector *pixelMask = psVectorAlloc(input->n,PS_TYPE_VECTOR_MASK);
    13871388  psStats  *stats     = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    13881389
    … …  
    13921393        pmReadout *ro  = stack->data[i];
    13931394        psImage *image = ro->image;
     1395        pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
    13941396        pixelData->data.F32[i] = image->data.F32[y][x];
     1397        if (isfinite(image->data.F32[y][x])&&
     1398            (fabs(image->data.F32[y][x]) < 1e5)) {
     1399          pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
     1400        }
     1401        else {
     1402          pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1;
     1403        }
     1404#if (0)
     1405      if ((x == 59)&&(y > 40)&&(y < 50)) {
     1406          fprintf(stderr,"%d %d %d %d %g\n",
     1407                  x,y,i,pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i],pixelData->data.F32[i]);
     1408        }
    13951409      }
    1396       if (!psVectorStats(stats,pixelData,NULL,NULL,0)) {
     1410#endif
     1411      if (!psVectorStats(stats,pixelData,NULL,pixelMask,1)) {
    13971412        psError(PS_ERR_UNKNOWN, false, "Unable to calculate median");
    13981413        psFree(stats);
    13991414        psFree(pixelData);
     1415        psFree(pixelMask);
    14001416        psFree(stack);
    14011417        return(false);
    … …  
    14031419      combined->image->data.F32[y][x] = stats->robustMedian;
    14041420      combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
     1421#if (0)
     1422      if ((x == 59)&&(y > 40)&&(y < 50)) {
     1423        fprintf(stderr,"%d %d %d %d %g\n",
     1424                x,y,-1,combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x],
     1425                combined->image->data.F32[y][x]);
     1426      }
     1427#endif
    14051428    }
    14061429  }
    … …  
    14081431  psFree(stats);
    14091432  psFree(pixelData);
     1433  psFree(pixelMask);
    14101434  psFree(stack);
    14111435  return (true);
Note: See TracChangeset for help on using the changeset viewer.