Changeset 34800

Timestamp:

Dec 11, 2012, 2:04:31 PM (13 years ago)

Author:

watersc1

Message:

Merge from my branch including background restoration and stack stage projection cell binned images.

Location:

trunk

Files:

• . (modified) (1 prop)
• Nebulous-Server/bin/neb-voladm (modified) (1 diff)
• Ohana (modified) (1 prop)
• Ohana/src/libohana/src (modified) (1 prop)
• Ohana/src/opihi (modified) (1 prop)
• Ohana/src/opihi/cmd.astro (modified) (1 prop)
• Ohana/src/opihi/cmd.data (modified) (1 prop)
• extsrc/gpcsw/gpcsrc/Make.Common (modified) (1 diff)
• ippMonitor (modified) (1 prop)
• ippScripts/Build.PL (modified) (1 diff)
• ippScripts/MANIFEST (modified) (1 diff)
• ippScripts/scripts/destreak_restore_camera.pl (modified) (1 prop)
• ippScripts/scripts/detrend_process_imfile.pl (modified) (2 diffs)
• ippScripts/scripts/ipp_apply_burntool_single.pl (modified) (1 prop)
• ippScripts/scripts/magic_destreak.pl (modified) (1 prop)
• ippScripts/scripts/regenerate_background.pl (copied) (copied from branches/czw_branch/20120906/ippScripts/scripts/regenerate_background.pl )
• ippScripts/scripts/skycell_jpeg.pl (modified) (2 diffs)
• ippScripts/scripts/stack_skycell.pl (modified) (2 diffs)
• ippScripts/scripts/warp_skycell.pl (modified) (7 diffs)
• ippTasks (modified) (1 prop)
• ippTasks/Makefile.am (modified) (1 diff)
• ippTasks/bg.regeneration.pro (copied) (copied from branches/czw_branch/20120906/ippTasks/bg.regeneration.pro )
• ippToPsps (modified) (1 prop)
• ippTools/share/Makefile.am (modified) (1 diff)
• ippTools/share/camtool_find_pendingimfile.sql (modified) (1 prop)
• ippTools/share/chiptool_setimfiletoupdate.sql (modified) (1 prop)
• ippTools/share/pxadmin_create_tables.sql (modified) (1 prop)
• ippTools/share/stacktool_inputskyfile.sql (modified) (1 diff)
• ippTools/share/stacktool_tobkg.sql (copied) (copied from branches/czw_branch/20120906/ippTools/share/stacktool_tobkg.sql )
• ippTools/share/warptool_towarped.sql (modified) (1 prop)
• ippTools/src (modified) (1 prop)
• ippTools/src/magictool.c (modified) (1 prop)
• ippTools/src/stacktool.c (modified) (4 diffs)
• ippTools/src/stacktool.h (modified) (1 diff)
• ippTools/src/stacktoolConfig.c (modified) (3 diffs)
• ippTools/src/warptool.c (modified) (3 diffs)
• ippTools/src/warptoolConfig.c (modified) (1 diff)
• ippconfig (modified) (1 prop)
• ippconfig/gpc1/pswarp.config (modified) (1 diff)
• ippconfig/recipes/filerules-mef.mdc (modified) (5 diffs)
• ippconfig/recipes/filerules-simple.mdc (modified) (5 diffs)
• ippconfig/recipes/filerules-split.mdc (modified) (5 diffs)
• ippconfig/recipes/nightly_science.config (modified) (4 diffs)
• ippconfig/recipes/ppBackground.mdc (modified) (1 diff)
• ippconfig/recipes/ppStack.config (modified) (1 diff)
• ippconfig/recipes/ppSub.config (modified) (1 prop)
• ippconfig/recipes/pswarp.config (modified) (1 diff)
• ippconfig/recipes/reductionClasses.mdc (modified) (1 prop)
• ppBackground/src/ppBackground.h (modified) (1 diff)
• ppBackground/src/ppBackgroundRestore.c (modified) (3 diffs)
• ppImage/src (modified) (1 prop)
• ppSkycell/src/ppSkycellCamera.c (modified) (1 diff)
• ppSkycell/src/ppSkycellLoop.c (modified) (12 diffs)
• ppStack/src/Makefile.am (modified) (3 diffs)
• ppStack/src/ppStack.h (modified) (2 diffs)
• ppStack/src/ppStackArguments.c (modified) (1 diff)
• ppStack/src/ppStackCamera.c (modified) (4 diffs)
• ppStack/src/ppStackCleanup.c (modified) (3 diffs)
• ppStack/src/ppStackCombineAlternate.c (copied) (copied from branches/czw_branch/20120906/ppStack/src/ppStackCombineAlternate.c )
• ppStack/src/ppStackCombinePrepare.c (modified) (3 diffs)
• ppStack/src/ppStackFiles.c (modified) (3 diffs)
• ppStack/src/ppStackLoop.c (modified) (5 diffs)
• ppStack/src/ppStackLoop.h (modified) (2 diffs)
• ppStack/src/ppStackMedian.c (copied) (copied from branches/czw_branch/20120906/ppStack/src/ppStackMedian.c )
• ppStack/src/ppStackMedianLoop.c (copied) (copied from branches/czw_branch/20120906/ppStack/src/ppStackMedianLoop.c )
• ppStack/src/ppStackOptions.h (modified) (1 diff)
• ppStack/src/ppStackPrepare.c (modified) (2 diffs)
• ppStack/src/ppStackReadout.c (modified) (1 diff)
• ppStack/src/ppStackSetup.c (modified) (2 diffs)
• ppStack/src/ppStackThread.c (modified) (4 diffs)
• ppStack/src/ppStackThread.h (modified) (1 diff)
• ppStack/src/ppStackUpdateHeader.c (modified) (2 diffs)
• ppViz/src/ppVizPSF/ppVizPSFCamera.c (modified) (2 diffs)
• ppViz/src/ppVizPSF/ppVizPSFLoop.c (modified) (9 diffs)
• psLib/src/imageops/psImageInterpolate.c (modified) (7 diffs)
• psLib/src/imageops/psImageInterpolate.h (modified) (1 diff)
• psLib/src/imageops/psImagePixelManip.c (modified) (3 diffs)
• psModules (modified) (1 prop)
• psModules/src/camera/pmFPABin.c (modified) (4 diffs)
• psModules/src/camera/pmFPAfileIO.c (modified) (1 diff)
• psModules/src/detrend/pmDark.c (modified) (3 diffs)
• psModules/src/imcombine/pmStack.c (modified) (1 diff)
• psModules/src/imcombine/pmStack.h (modified) (1 diff)
• psModules/src/objects/Makefile.am (modified) (1 diff)
• psModules/src/objects/pmFootprintCullPeaks.c (modified) (1 diff)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotStackImageLoop.c (modified) (1 prop)
• pstamp/scripts (modified) (1 prop)
• psvideophot (modified) (1 prop)
• pswarp/src/pswarp.c (modified) (2 diffs)
• pswarp/src/pswarp.h (modified) (2 diffs)
• pswarp/src/pswarpArguments.c (modified) (4 diffs)
• pswarp/src/pswarpDefine.c (modified) (3 diffs)
• pswarp/src/pswarpFileNames.h (modified) (2 diffs)
• pswarp/src/pswarpLoop.c (modified) (5 diffs)
• pswarp/src/pswarpParseCamera.c (modified) (2 diffs)
• pswarp/src/pswarpTransformReadout.c (modified) (6 diffs)
• pswarp/src/pswarpTransformSources.c (modified) (1 diff)
• pswarp/src/pswarpTransformTile.c (modified) (2 diffs)

trunk/Nebulous-Server/bin/neb-voladm

@@ -68 +68 @@
         if (defined $allocate and $allocate !~ m/^[01]$/)
         or (defined $available and $available !~ m/^[01]$/)
-        or (defined $xattr and $xattr !~ m/^[012]$/);
+        or (defined $xattr and $xattr !~ m/^[0123]$/);
 }
 

trunk/extsrc/gpcsw/gpcsrc/Make.Common

@@ -156 +156 @@
 # compilation and -shared applies to the link stage.
 #
-FCDEBUG  = -g -O -fno-automatic
+FCDEBUG  = -g -O0 -fno-automatic
 FFWARN   = -Wall
 FFLAGS   = $(FCDEBUG) $(EXTRA_FFLAGS) $(FFWARN) $(CCDEFS) $(CCINCS)
-CCDEBUG  = -g -O
+CCDEBUG  = -g -O0
 CCSHARE  = -fPIC -shared
 CCWARN   = -Wall -Wstrict-prototypes # -Wshadow

trunk/ippScripts/Build.PL

@@ -128 +128 @@
         scripts/listvideocells.pl
         scripts/skycalibration.pl
+        scripts/regenerate_background.pl
     )],
     dist_abstract => 'Scripts for running the Pan-STARRS IPP',

trunk/ippScripts/MANIFEST

@@ -46 +46 @@
 scripts/skycell_jpeg.pl
 scripts/lap_science.pl
+scripts/regenerate_background.pl
 t/00_distribution.t

trunk/ippScripts/scripts/detrend_process_imfile.pl

@@ -27 +27 @@
 my $ppImage = can_run('ppImage') or (warn "Can't find ppImage" and $missing_tools = 1);
 my $ppStatsFromMetadata = can_run('ppStatsFromMetadata') or (warn "Can't find ppStatsFromMetadata" and $missing_tools = 1);
+my $nebrepair = can_run('neb-repair') or (warn "Can't find neb-repair" and $missing_tools = 1);
+
 if ($missing_tools) {
     warn("Can't find required tools.");
@@ -124 +126 @@
 # Run ppImage
 unless ($no_op) {
+    my $repair_cmd = "$nebrepair $input_uri";
+    my ($repair_success, $repair_error_code, $repair_full_buf, $repair_stdout_buf, $repair_stderr_buf ) = run(command => $repair_cmd, verbose => $verbose);
+    unless ($repair_success) {
+        &my_die("Unable to attempt repair: $input_uri $repair_error_code", $det_id, $exp_id, $class_id $PS_EXIT_SYS_ERROR);
+    }
+
+
     my $command = "$ppImage -file $input_uri $outroot";
     $command .= " -recipe PPIMAGE $ppimage_recipe";

trunk/ippScripts/scripts/skycell_jpeg.pl

@@ -268 +268 @@
 
     my %tangents = ();
+
+    my %products = ('image' => "PPSTACK.UNCONV",
+                    'mask'  => "PPSTACK.UNCONV.MASK",
+                    'variance' => "PPSTACK.UNCONV.VARIANCE",
+                    'exp'   => "PPSTACK.UNCONV.EXP",
+                    'num'   => "PPSTACK.UNCONV.EXPNUM",
+                    'bkg'   => "PPSTACK.OUTPUT.BKGMODEL"
+        );
     
     foreach my $imfile (@$imfiles) {
@@ -281 +289 @@
 
         $projection_cell =~ s/^(.*)\..*$/$1/;
-        
+
         unless (exists($tangents{$projection_cell})) {
             # Make a temp file and fill, but be sure to save 
-            ($tempFile, $tempName) = tempfile("/tmp/skycell.$projection_cell.XXXX",
-                                                 UNLINK => !$save_temps);
-            $tangents{$projection_cell}{FILE} = $tempFile;
-            $tangents{$projection_cell}{NAME} = $tempName;
-            if ($masks) {
-                my ($maskFile, $maskName) = tempfile("/tmp/skycell.$projection_cell.masks.XXXX",
-                                                     UNLINK => !$save_temps);
-                $tangents{$projection_cell}{MFILE} = $maskFile;
-                $tangents{$projection_cell}{MNAME} = $maskName;
-            }           
-        }
-        print "$skycell_id $projection_cell\n"; 
-        my $file = $ipprc->filename("PPSTACK.OUTPUT", $path_base, $skycell_id);
-        print "$file $state $quality\n";
-        my $f_fh = $tangents{$projection_cell}{FILE};
-        print $f_fh "$file\n";
-        if ($masks) {
-            my $mask = $ipprc->filename("PPSTACK.OUTPUT.MASK", $path_base, $skycell_id);
-            print "$mask\n";
-            my $m_fh = $tangents{$projection_cell}{MFILE};
-            print $m_fh "$mask\n";
+            foreach my $key (keys %products) {
+                ($tempFile, $tempName) = tempfile("/tmp/skycell.$projection_cell.$key.XXXX",
+                                                  UNLINK => !$save_temps);
+                $tangents{$projection_cell}{$key}{FILE} = $tempFile;
+                $tangents{$projection_cell}{$key}{NAME} = $tempName;
+            }
+        }
+        foreach my $key (keys %products) {
+            print "$skycell_id $projection_cell\n";     
+            my $file = $ipprc->filename($products{$key}, $path_base, $skycell_id);
+            print "$file $state $quality\n";
+            my $f_fh = $tangents{$projection_cell}{$key}{FILE};
+            print $f_fh "$file\n";
         }
     }
     foreach my $projection_cell (keys %tangents) {
-        $command = "$ppSkycell -images $tangents{$projection_cell}{NAME}";
-        if ($masks) {
-            $command .= " -masks $tangents{$projection_cell}{MNAME} ";
-        }
-        $command .= " ${outroot}.${projection_cell} ";
-        print "$command\n";
-        ( $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 ppSkycell: $error_code", $stage_id, $error_code);
-        }
+        ## Loop over results here.
+        # Images
+        # Masks
+        # Variances
+        # Nexptime
+        # Nexp
+        # Backgrounds
+        
+        foreach my $key (keys %products) {
+            $command = "$ppSkycell -images $tangents{$projection_cell}{$key}{NAME}";
+            $command .= " ${outroot}.${projection_cell}.${key} ";
+            if ($key eq 'bkg') {
+                $command .= " -Di BIN1 1 -Di BIN2 1 ";
+            }
+            elsif ($key eq 'image') {
+                $command .= " -masks $tangents{$projection_cell}{mask}{NAME} ";
+            }
+            elsif ($key eq 'mask') { 
+                next; # This should be made with the images.
+            }
+            print "$command\n";
+            ( $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 ppSkycell: $error_code", $stage_id, $error_code);
+            }
+        }
+
         # Update database:
         $command = "$stacktool -addsummary -sass_id $stage_id -projection_cell $projection_cell -path_base $outroot";

trunk/ippScripts/scripts/stack_skycell.pl

@@ -229 +229 @@
     my $sources = $ipprc->filename("PSWARP.OUTPUT.SOURCES", $file->{path_base}); # Sources name
 
+    my $bkgmodel = $ipprc->filename("PSWARP.OUTPUT.BKGMODEL", $file->{path_base});
+
     &my_die("Image $image does not exist", $stack_id, $PS_EXIT_SYS_ERROR) unless $ipprc->file_exists( $image );
     &my_die("Mask $mask does not exist", $stack_id, $PS_EXIT_SYS_ERROR) unless $ipprc->file_exists( $mask );
@@ -240 +242 @@
     print $listFile "\tPSF\tSTR\t" . $psf . "\n" if $convolve;
     print $listFile "\tSOURCES\tSTR\t" . $sources . "\n";
+    print $listFile "\tBKGMODEL\tSTR\t" . $bkgmodel . "\n" if $ipprc->file_exists( $bkgmodel );
 
     print $listFile "END\n\n";

trunk/ippScripts/scripts/warp_skycell.pl

@@ -128 +128 @@
 # Where do we get the astrometry source from?
 my $astromSource;               # The astrometry source
+my $doBackground;               # Do we want to make background models?
 {
     my $command = "$ppConfigDump -camera $camera -recipe PSWARP $recipe_pswarp -dump-recipe PSWARP -";
@@ -139 +140 @@
         &my_die("Unable to parse metadata config doc", $warp_id, $skycell_id, $tess_dir, $PS_EXIT_PROG_ERROR);
     $astromSource = metadataLookupStr($metadata, 'ASTROM.SOURCE');
+    $doBackground = metadataLookupBool($metadata, 'BACKGROUND.MODEL');    
 }
 
@@ -168 +170 @@
 if ($do_stats) {
     $outputStats = prepare_output ("SKYCELL.STATS", $outroot, $skycell_id, 1) if $do_stats;
+}
+my $outputBKGs;
+if ($doBackground) {
+    $outputBKGs = prepare_output ("PSWARP.OUTPUT.BKGMODEL", $outroot, $skycell_id, 1);
 }
 my $configuration;
@@ -209 +215 @@
 my ($weightFile, $weightName) = tempfile( "$tempOutRoot.weight.list.XXXX", UNLINK => !$save_temps);
 my ($astromFile, $astromName) = tempfile( "$tempOutRoot.astrom.list.XXXX", UNLINK => !$save_temps);
-
+my ($bkgFile, $bkgName);
+if ($doBackground) {
+    ($bkgFile, $bkgName) = tempfile( "$tempOutRoot.bkg.list.XXXX", UNLINK => !$save_temps);
+}
 my $wrote_astrom = 0;
 foreach my $imfile (@$imfiles) {
@@ -235 +244 @@
     print $maskFile   "$mask\n";
     print $weightFile "$weight\n";
+    my $bkg;
+    if ($doBackground) {
+        $bkg    = $ipprc->filename("PSPHOT.BACKMDL", $imfile->{chip_path_base}, $imfile->{class_id});
+        print $bkgFile "$bkg\n";
+    }
 
     if (!$wrote_astrom) {
@@ -245 +259 @@
 close $weightFile;
 close $astromFile;
-
+if ($doBackground) {
+    close($bkgFile);
+}
 # We need the recipe to determine if we care whether the PSF is generated or not
 my $recipe;
@@ -269 +285 @@
     $command .= " -variancelist $weightName";
     $command .= " -astromlist $astromName";
+    $command .= " -bkglist $bkgName" if ($doBackground);
     $command .= " $outroot $skyFile";
     $command .= " -F PSPHOT.PSF.SAVE PSPHOT.PSF.SKY.SAVE";

trunk/ippTasks/Makefile.am

@@ -45 +45 @@
         diffphot.pro \
         lap.pro \
-        vp.pro
+        vp.pro \
+        bg.regeneration.pro
 
 other_files = \

trunk/ippTools/share/Makefile.am

@@ -393 +393 @@
         stacktool_sassskyfile.sql \
         stacktool_tosum.sql \
+        stacktool_tobkg.sql \
         stacktool_tosummary.sql \
         stacktool_addsummary.sql \

trunk/ippTools/share/stacktool_inputskyfile.sql

@@ -4 +4 @@
     rawExp.exp_id,
     rawExp.exp_name,
+    rawExp.exp_time,
     rawExp.object,
     rawExp.dateobs,

trunk/ippTools/src/stacktool.c

@@ -37 +37 @@
 static bool inputskyfileMode(pxConfig *config);
 static bool tosumMode(pxConfig *config);
+static bool tobkgMode(pxConfig *config);
 static bool addsumskyfileMode(pxConfig *config);
 static bool sumskyfileMode(pxConfig *config);
@@ -76 +77 @@
         MODECASE(STACKTOOL_MODE_INPUTSKYFILE,          inputskyfileMode);
         MODECASE(STACKTOOL_MODE_TOSUM,                 tosumMode);
+        MODECASE(STACKTOOL_MODE_TOBKG,                 tobkgMode);
         MODECASE(STACKTOOL_MODE_ADDSUMSKYFILE,         addsumskyfileMode);
         MODECASE(STACKTOOL_MODE_SUMSKYFILE,            sumskyfileMode);
@@ -952 +954 @@
         // negative simple so the default is true
         if (!ippdbPrintMetadatas(stdout, output, "stackSumSkyfile", !simple)) {
+            psError(PS_ERR_UNKNOWN, false, "failed to print array");
+            psFree(output);
+            return false;
+        }
+    }
+
+    psFree(output);
+
+    return true;
+}
+
+static bool tobkgMode(pxConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(config, false);
+
+    psMetadata *where = psMetadataAlloc();
+    PXOPT_COPY_S64(config->args, where, "-stack_id", "stack_id", "==");
+    pxAddLabelSearchArgs (config, where, "-label", "stackRun.label", "==");
+
+    PXOPT_LOOKUP_U64(limit, config->args, "-limit", false, false);
+    PXOPT_LOOKUP_BOOL(simple, config->args, "-simple", false);
+
+    psString query = pxDataGet("stacktool_tobkg.sql");
+    if (!query) {
+        psError(PXTOOLS_ERR_SYS, false, "failed to retreive SQL statement");
+        return false;
+    }
+
+    psString whereClause = psStringCopy(""); // WHERE conditions to add
+    if (psListLength(where->list)) {
+        psString new = psDBGenerateWhereConditionSQL(where, NULL);
+        psStringAppend(&whereClause, "\nAND %s", new);
+        psFree(new);
+    }
+    psFree(where);
+
+    psStringAppend(&query, "\nORDER by priority DESC, stack_id");
+
+    // treat limit == 0 as "no limit"
+    if (limit) {
+        psString limitString = psDBGenerateLimitSQL(limit);
+        psStringAppend(&query, " %s", limitString);
+        psFree(limitString);
+    }
+
+    if (!p_psDBRunQueryF(config->dbh, query, whereClause)) {
+        psError(PS_ERR_UNKNOWN, false, "database error");
+        psFree(query);
+        return false;
+    }
+    psFree(query);
+
+    psArray *output = p_psDBFetchResult(config->dbh);
+    if (!output) {
+        psErrorCode err = psErrorCodeLast();
+        switch (err) {
+            case PS_ERR_DB_CLIENT:
+                psError(PXTOOLS_ERR_SYS, false, "database error");
+            case PS_ERR_DB_SERVER:
+                psError(PXTOOLS_ERR_PROG, false, "database error");
+            default:
+                psError(PXTOOLS_ERR_PROG, false, "unknown error");
+        }
+
+        return false;
+    }
+    if (!psArrayLength(output)) {
+        psTrace("stacktool", PS_LOG_INFO, "no rows found");
+        psFree(output);
+        return true;
+    }
+
+    if (psArrayLength(output)) {
+        // negative simple so the default is true
+        if (!ippdbPrintMetadatas(stdout, output, "stackBkgSkyfile", !simple)) {
             psError(PS_ERR_UNKNOWN, false, "failed to print array");
             psFree(output);
@@ -1681 +1758 @@
     PXOPT_LOOKUP_S16(fault, config->args, "-fault", true, false);
     PXOPT_LOOKUP_S16(quality, config->args, "-set_quality", false, false);
+    PXOPT_LOOKUP_S16(background_model, config->args, "-set_background_model", false, false);
 
     psMetadata *where = psMetadataAlloc();
     PXOPT_COPY_S64(config->args, where, "-stack_id",   "stack_id",   "==");
 
-    if (!pxSetFaultCode(config->dbh, "stackSumSkyfile", where, fault, quality)) {
+    if (background_model) {
+      psMetadata *values = psMetadataAlloc();
+      PXOPT_COPY_S16(config->args, values, "-set_background_model", "background_model", "==");
+      long rows = psDBUpdateRows(config->dbh,"stackSumSkyfile", where, values);
+      psFree(values);
+      if (!rows) {
+        // This maybe should rollback and error if rows != 1
+        psError(PS_ERR_UNKNOWN, true, "no rows changed");
+        return false;
+      }
+    }
+    else {    
+      if (!pxSetFaultCode(config->dbh, "stackSumSkyfile", where, fault, quality)) {
         psError(PS_ERR_UNKNOWN, false, "failed to set set fault flag");
         psFree (where);
         return false;
+      }
     }
     psFree (where);

trunk/ippTools/src/stacktool.h

@@ -31 +31 @@
     STACKTOOL_MODE_INPUTSKYFILE,
     STACKTOOL_MODE_TOSUM,
+    STACKTOOL_MODE_TOBKG,
     STACKTOOL_MODE_ADDSUMSKYFILE,
     STACKTOOL_MODE_SUMSKYFILE,

trunk/ippTools/src/stacktoolConfig.c

@@ -150 +150 @@
     psMetadataAddBool(tosumArgs, PS_LIST_TAIL, "-simple",  0,            "use the simple output format", false);
 
+    // -tobkg
+    psMetadata *tobkgArgs = psMetadataAlloc();
+    psMetadataAddS64(tobkgArgs, PS_LIST_TAIL, "-stack_id", 0,            "search by stack ID", 0);
+    psMetadataAddStr(tobkgArgs, PS_LIST_TAIL, "-label", PS_META_DUPLICATE_OK, "search by label", NULL);
+    psMetadataAddU64(tobkgArgs, PS_LIST_TAIL, "-limit",  0,            "limit result set to N items", 0);
+    psMetadataAddBool(tobkgArgs, PS_LIST_TAIL, "-simple",  0,            "use the simple output format", false);
+    
     // -addsumskyfile
     psMetadata *addsumskyfileArgs = psMetadataAlloc();
@@ -280 +287 @@
     psMetadataAddS16(updatesumskyfileArgs, PS_LIST_TAIL, "-fault", 0,            "set fault code (required)", 0);
     psMetadataAddS16(updatesumskyfileArgs, PS_LIST_TAIL, "-set_quality", 0,            "set quality", 0);
-
+    psMetadataAddS16(updatesumskyfileArgs, PS_LIST_TAIL, "-set_background_model", 0,   "set background model", 0);
+    
     // -exportrun
     psMetadata *exportrunArgs = psMetadataAlloc();
@@ -304 +312 @@
     PXOPT_ADD_MODE("-inputskyfile",    "", STACKTOOL_MODE_INPUTSKYFILE,    inputskyfileArgs);
     PXOPT_ADD_MODE("-tosum",           "", STACKTOOL_MODE_TOSUM,          tosumArgs);
+    PXOPT_ADD_MODE("-tobkg",           "", STACKTOOL_MODE_TOBKG,          tobkgArgs);
     PXOPT_ADD_MODE("-addsumskyfile",   "", STACKTOOL_MODE_ADDSUMSKYFILE,   addsumskyfileArgs);
     PXOPT_ADD_MODE("-sumskyfile",      "list results of stackRun", STACKTOOL_MODE_SUMSKYFILE,      sumskyfileArgs);

trunk/ippTools/src/warptool.c

@@ -2059 +2059 @@
         PXOPT_LOOKUP_S64(magicked, config->args, "-set_magicked", false, false);
         PXOPT_LOOKUP_S16(quality, config->args, "-set_quality", false, false);
+
         if (magicked) {
             psStringAppend(&set_magicked_skyfile, "\n , warpSkyfile.magicked = %" PRId64, magicked);
@@ -2065 +2066 @@
         if (quality) {
           psStringAppend(&set_magicked_skyfile, "\n , warpSkyfile.quality = %"PRId16, quality);
+        }
+        PXOPT_LOOKUP_S16(background_model, config->args, "-set_background_model", false, false);
+        if (background_model) {
+          psStringAppend(&set_magicked_skyfile, "\n , warpSkyfile.background_model = %"PRId16, background_model);
         }
     } else if (!strcmp(data_state, "cleaned") || !strcmp(data_state, "purged")) {
@@ -2125 +2130 @@
     // warp_id, skycell_id, fault are required
     PXOPT_LOOKUP_STR(state, config->args, "-set_state", false, false);
-
-    if (!state) {
+    PXOPT_LOOKUP_S16(background_model, config->args, "-set_background_model", false, false);
+    if (background_model) {
+      // CZW 2012-12-06: I'm unclear why we don't use this form for all updates?
+      psMetadata *where = psMetadataAlloc();
+      psMetadata *values = psMetadataAlloc();
+      PXOPT_COPY_S64(config->args, where, "-warp_id", "warp_id", "==");
+      PXOPT_COPY_STR(config->args, where, "-skycell_id", "skycell_id", "==");
+      PXOPT_COPY_S16(config->args, values, "-set_background_model", "background_model", "==");
+      long rows = psDBUpdateRows(config->dbh,"warpSkyfile", where, values);
+      psFree(values);
+      psFree(where);
+      if (!rows) {
+        // This maybe should rollback and error if rows != 1
+        psError(PS_ERR_UNKNOWN, true, "no rows changed");
+        return false;
+      }
+    }
+    else if (!state) {
       psMetadata *where = psMetadataAlloc();
       PXOPT_COPY_S64(config->args, where, "-warp_id", "warp_id", "==");

trunk/ippTools/src/warptoolConfig.c

@@ -410 +410 @@
     psMetadataAddS16(updateskyfileArgs, PS_LIST_TAIL, "-set_quality",  0,"new quality value", 0);
     psMetadataAddStr(updateskyfileArgs, PS_LIST_TAIL, "-set_state", 0,   "set state", 0);
-
+    psMetadataAddS16(updateskyfileArgs, PS_LIST_TAIL, "-set_background_model", 0, "set the background_model value", 0);
     // -exportrun
     psMetadata *exportrunArgs = psMetadataAlloc();

trunk/ippconfig/gpc1/pswarp.config

@@ -1 +1 @@
 ASTROM.SOURCE           STR     PSASTRO.OUTPUT  # Source file rule for astrometry, or NULL
+BACKGROUND.MODEL        BOOL    TRUE
+BKG.XGRID               S32     50
+BKG.YGRID               S32     50

trunk/ippconfig/recipes/filerules-mef.mdc

@@ -120 +120 @@
 PSWARP.SKYCELL          INPUT    @FILES        CHIP       IMAGE
 PSWARP.ASTROM           INPUT    @FILES        CHIP       CMF
+PSWARP.BKGMODEL           INPUT    @FILES        CHIP       IMAGE
 
 ## files used by ppsub
@@ -138 +139 @@
 PPSTACK.INPUT.PSF       INPUT    @FILES        CHIP       PSF
 PPSTACK.INPUT.SOURCES   INPUT    @FILES        FPA        CMF
+PPSTACK.INPUT.BKGMODEL    INPUT    @FILES        FPA        IMAGE
 
 ## files used by ppstamp
@@ -296 +298 @@
 PSWARP.BIN2             OUTPUT {OUTPUT}.b2.fits                  IMAGE     COMP_IMG   FPA        TRUE      NONE
 PSWARP.CONFIG           OUTPUT {OUTPUT}.pswarp.mdc               TEXT      NONE       FPA        TRUE      NONE
+PSWARP.OUTPUT.BKGMODEL  OUTPUT {OUTPUT}.mdl.fits                      IMAGE           NONE       FPA        TRUE      NONE
                                                                                      
 SKYCELL.STATS           OUTPUT {OUTPUT}.stats                    STATS     NONE       FPA        TRUE      NONE
@@ -342 +345 @@
 PPSTACK.OUTPUT.JPEG1    OUTPUT {OUTPUT}.b1.jpg                   JPEG      NONE       FPA        TRUE      NONE
 PPSTACK.OUTPUT.JPEG2    OUTPUT {OUTPUT}.b2.jpg                   JPEG      NONE       FPA        TRUE      NONE
+PPSTACK.OUTPUT.BKGMODEL OUTPUT {OUTPUT}.mdl.fits                 IMAGE     NONE       FPA        TRUE      NONE
+PPSTACK.OUTPUT.BKGREST  OUTPUT {OUTPUT}.bkgrest.fits             IMAGE            NONE       FPA        TRUE      NONE
 PPSTACK.CONFIG          OUTPUT {OUTPUT}.ppStack.mdc              TEXT      NONE       FPA        TRUE      NONE
                                                                                      
@@ -370 +375 @@
 PPSKYCELL.BIN1          OUTPUT {OUTPUT}.{FILE.INDEX}.b1.fits     IMAGE     NONE       FPA        TRUE      NONE
 PPSKYCELL.BIN2          OUTPUT {OUTPUT}.{FILE.INDEX}.b2.fits     IMAGE     NONE       FPA        TRUE      NONE
+PPSKYCELL.BIN1.MASK          OUTPUT {OUTPUT}.{FILE.INDEX}.b1.mk.fits     MASK     COMP_MASK       FPA        TRUE      NONE
+PPSKYCELL.BIN2.MASK          OUTPUT {OUTPUT}.{FILE.INDEX}.b2.mk.fits     MASK     COMP_MASK       FPA       TRUE      NONE
 
 LOG.IMFILE              OUTPUT {OUTPUT}.{CHIP.NAME}.log          TEXT      NONE       CHIP       TRUE      NONE

trunk/ippconfig/recipes/filerules-simple.mdc

@@ -90 +90 @@
 PSWARP.SKYCELL            INPUT    @FILES        FPA        IMAGE     
 PSWARP.ASTROM             INPUT    @FILES        FPA        CMF
+PSWARP.BKGMODEL           INPUT    @FILES        CHIP       IMAGE
 
 ## files used by ppsub
@@ -108 +109 @@
 PPSTACK.INPUT.PSF         INPUT    @FILES        CHIP       PSF
 PPSTACK.INPUT.SOURCES     INPUT    @FILES        FPA        CMF
+PPSTACK.INPUT.BKGMODEL    INPUT    @FILES        FPA        IMAGE
 
 ## files used by ppstamp
@@ -251 +253 @@
 PSWARP.BIN2                  OUTPUT {OUTPUT}.b2.fits              IMAGE           NONE       FPA        TRUE      NONE
 PSWARP.CONFIG                OUTPUT {OUTPUT}.pswarp.mdc           TEXT            NONE       FPA        TRUE      NONE
+PSWARP.OUTPUT.BKGMODEL       OUTPUT {OUTPUT}.mdl.fits             IMAGE           NONE   FPA        TRUE      NONE
                                                      
 SKYCELL.STATS                OUTPUT {OUTPUT}.stats                STATS           NONE       FPA        TRUE      NONE
@@ -300 +303 @@
 PPSTACK.OUTPUT.JPEG1         OUTPUT {OUTPUT}.b1.jpg               JPEG            NONE       FPA        TRUE      NONE
 PPSTACK.OUTPUT.JPEG2         OUTPUT {OUTPUT}.b2.jpg               JPEG            NONE       FPA        TRUE      NONE
+PPSTACK.OUTPUT.BKGMODEL      OUTPUT {OUTPUT}.mdl.fits             IMAGE           NONE       FPA        TRUE      NONE
+PPSTACK.OUTPUT.BKGREST       OUTPUT {OUTPUT}.bkgrest.fits         IMAGE           NONE       FPA        TRUE      NONE
 PPSTACK.CONFIG               OUTPUT {OUTPUT}.ppStack.mdc          TEXT            NONE       FPA        TRUE      NONE
                                              
@@ -327 +332 @@
 PPSKYCELL.BIN1               OUTPUT {OUTPUT}.{FILE.INDEX}.b1.fits     IMAGE     NONE       FPA        TRUE      NONE
 PPSKYCELL.BIN2               OUTPUT {OUTPUT}.{FILE.INDEX}.b2.fits     IMAGE     NONE       FPA        TRUE      NONE
+PPSKYCELL.BIN1.MASK          OUTPUT {OUTPUT}.b1.mk.fits     MASK     COMP_MASK       FPA        TRUE      NONE
+PPSKYCELL.BIN2.MASK          OUTPUT {OUTPUT}.b2.mk.fits     MASK     COMP_MASK       FPA       TRUE      NONE
 
 LOG.IMFILE                   OUTPUT {OUTPUT}.imfile.log           TEXT            NONE       FPA        TRUE      NONE

trunk/ippconfig/recipes/filerules-split.mdc

@@ -108 +108 @@
 PSWARP.SKYCELL            INPUT    @FILES        FPA        IMAGE
 PSWARP.ASTROM             INPUT    @FILES        CHIP       CMF
+PSWARP.BKGMODEL           INPUT    @FILES        CHIP       IMAGE
                           
 ## files used by ppsub    
@@ -126 +127 @@
 PPSTACK.INPUT.PSF         INPUT    @FILES        CHIP       PSF
 PPSTACK.INPUT.SOURCES     INPUT    @FILES        FPA        CMF
+PPSTACK.INPUT.BKGMODEL    INPUT    @FILES        FPA        IMAGE
                           
 ## files used by ppstamp  
@@ -278 +280 @@
 PSWARP.CONFIG                OUTPUT {OUTPUT}.pswarp.mdc               TEXT            NONE       FPA        TRUE      NONE
 PSWARP.OUTPUT.UNCOMPRESSED   OUTPUT {OUTPUT}.fits                     IMAGE           NONE       FPA        TRUE      NONE
+PSWARP.OUTPUT.BKGMODEL       OUTPUT {OUTPUT}.mdl.fits                 IMAGE           NONE       FPA        TRUE      NONE
                                                                                                    
 SKYCELL.STATS                OUTPUT {OUTPUT}.stats                    STATS           NONE       FPA        TRUE      NONE
@@ -333 +336 @@
 PPSTACK.OUTPUT.JPEG1         OUTPUT {OUTPUT}.b1.jpg                   JPEG            NONE       FPA        TRUE      NONE
 PPSTACK.OUTPUT.JPEG2         OUTPUT {OUTPUT}.b2.jpg                   JPEG            NONE       FPA        TRUE      NONE
+PPSTACK.OUTPUT.BKGMODEL      OUTPUT {OUTPUT}.mdl.fits                 IMAGE           NONE       FPA        TRUE      NONE
+PPSTACK.OUTPUT.BKGREST       OUTPUT {OUTPUT}.bkgrest.fits             IMAGE           NONE       FPA        TRUE      NONE
 PPSTACK.CONFIG               OUTPUT {OUTPUT}.ppStack.mdc              TEXT            NONE       FPA        TRUE      NONE
 
@@ -360 +365 @@
 PPSKYCELL.JPEG1              OUTPUT {OUTPUT}.0.b1.jpeg     JPEG            NONE       CHIP       TRUE      NONE
 PPSKYCELL.JPEG2              OUTPUT {OUTPUT}.0.b2.jpeg     JPEG            NONE       CHIP       TRUE      NONE
-PPSKYCELL.BIN1          OUTPUT {OUTPUT}.0.b1.fits     IMAGE     COMP_IMG       FPA        TRUE      NONE
-PPSKYCELL.BIN2          OUTPUT {OUTPUT}.0.b2.fits     IMAGE     COMP_IMG       FPA       TRUE      NONE
+PPSKYCELL.BIN1               OUTPUT {OUTPUT}.b1.fits     IMAGE     COMP_IMG       FPA        TRUE      NONE
+PPSKYCELL.BIN2               OUTPUT {OUTPUT}.b2.fits     IMAGE     COMP_IMG       FPA       TRUE      NONE
+PPSKYCELL.BIN1.MASK          OUTPUT {OUTPUT}.b1.mk.fits     MASK     COMP_MASK       FPA        TRUE      NONE
+PPSKYCELL.BIN2.MASK          OUTPUT {OUTPUT}.b2.mk.fits     MASK     COMP_MASK       FPA       TRUE      NONE
 
 LOG.IMFILE                   OUTPUT {OUTPUT}.{CHIP.NAME}.log          TEXT            NONE       CHIP       TRUE      NONE

trunk/ippconfig/recipes/nightly_science.config

@@ -229 +229 @@
   DIFFABLE  BOOL TRUE
   REDUCTION STR SWEETSPOT
+  ADDITIONAL_STACK_LABEL STR ecliptic.rp
 #  REDUCTION STR SSTF_4SIG
 #  WARP      S16 60
@@ -240 +241 @@
   DIFFABLE    BOOL TRUE
   REDUCTION STR SWEETSPOT
+  ADDITIONAL_STACK_LABEL STR ecliptic.rp
 #  REDUCTION STR SSTF_4SIG
 END
@@ -250 +252 @@
   DIFFABLE    BOOL TRUE
   REDUCTION STR SWEETSPOT
+  ADDITIONAL_STACK_LABEL STR ecliptic.rp
 #  REDUCTION STR SSTF_4SIG
 END
@@ -260 +263 @@
   DIFFABLE    BOOL TRUE
   REDUCTION STR SWEETSPOT
+  ADDITIONAL_STACK_LABEL STR ecliptic.rp
 #  REDUCTION STR SSTF_4SIG
 END

trunk/ippconfig/recipes/ppBackground.mdc

@@ -1 +1 @@
 # Recipe options for ppBackground
 
+BINNING_RECIPE  STR     PSPHOT
+BINNING_XNAME   STR     BACKGROUND.XBIN
+BINNING_YNAME   STR     BACKGROUND.YBIN
+
+DEFAULT    METADATA
+END
+
+
+STACK       METADATA
+  BINNING_RECIPE        STR     PSWARP
+  BINNING_XNAME STR     BKG.XGRID       
+  BINNING_YNAME STR     BKG.YGRID
+END
+            
 
 BACKGROUND       METADATA

trunk/ippconfig/recipes/ppStack.config

@@ -103 +103 @@
 STACK.TYPE      STR     DEEP_STACK
 
+DEFAULT.EXPTIME  F32    NAN
+
 # Recipe overrides for STACK
 STACK   METADATA

trunk/ippconfig/recipes/pswarp.config

@@ -19 +19 @@
 MASKSTAT.ADVISORY   U32 0x080
 
+BACKGROUND.MODEL    BOOL   FALSE             # Construct a warped version of the chip stage background
+BKG.XGRID              S32    4               # These need to be tuned
+BKG.YGRID              S32    4
 
 # Default recipe for warping

trunk/ppBackground/src/ppBackground.h

@@ -37 +37 @@
     );
 
+/// Determine the binning from the recipe if available.
+psImageBinning *ppBackgroundBinningByRecipe(const psImage *image, // Image for which to generate a bg model
+                                            const pmConfig *config, // Configuration
+                                            psString recipe_name,
+                                            psString Xbin_name,
+                                            psString Ybin_name
+                                            );
+
+
 /// Restore the background to an image
 bool ppBackgroundRestore(

trunk/ppBackground/src/ppBackgroundRestore.c

@@ -5 +5 @@
 
 #include "ppBackground.h"
+
+psImageBinning *ppBackgroundBinningByRecipe(const psImage *image, // Image for which to generate a bg model
+                                            const pmConfig *config, // Configuration
+                                            psString recipe_name,
+                                            psString Xbin_name,
+                                            psString Ybin_name
+                                            )
+{
+  bool status = true;
+
+  psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, recipe_name);
+  assert (recipe);
+
+  // I have the fine image size, I know the binning factor, determine the ruff image size
+  psImageBinning *binning = psImageBinningAlloc();
+  binning->nXfine = image->numCols;
+  binning->nYfine = image->numRows;
+  binning->nXbin  = psMetadataLookupS32(&status, recipe, Xbin_name);
+  binning->nYbin  = psMetadataLookupS32(&status, recipe, Ybin_name);
+
+  psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER);
+  psImageBinningSetSkip(binning, image);
+
+  return binning;
+}
+
+  
+
 
 bool ppBackgroundRestore(pmChip *chip, const pmChip *background, const pmChip *pattern,
@@ -27 +55 @@
     if (background) {
         pmReadout *bgRO = pmFPAviewThisReadout(view, background->parent); // Readout with background
-        psImageBinning *binning = psphotBackgroundBinning(image, config);
+
+        psImageBinning *binning;
+        psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPBACKGROUND_RECIPE); // Recipe
+        if (!recipe) {
+          binning = psphotBackgroundBinning(image, config);
+        }
+        else {
+          binning = ppBackgroundBinningByRecipe(image,config,
+                                                psMetadataLookupStr(NULL,recipe,"BINNING_RECIPE"),
+                                                psMetadataLookupStr(NULL,recipe,"BINNING_XNAME"),
+                                                psMetadataLookupStr(NULL,recipe,"BINNING_YNAME"));
+        }
+                                                
+        fprintf(stderr,"%d %d %d %d\n",binning->nXfine,binning->nYfine,binning->nXbin,binning->nYbin);
         if (!binning) {
             psError(psErrorCodeLast(), false, "Unable to find background binning");
@@ -49 +90 @@
         for (int y = 0; y < numRows; y++) {
             for (int x = 0; x < numCols; x++) {
+              if ((y % 250 == 0)&&(x % 250 == 0)) {
+                printf("%d %d %g\n",x,y,bgImage->data.F32[y][x]);
+              }
                 image->data.F32[y][x] += bgImage->data.F32[y][x];
+                
             }
         }

trunk/ppSkycell/src/ppSkycellCamera.c

@@ -103 +103 @@
         return false;
     }
-    if (!pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.BIN1")) {
+    pmFPAfile *bin1 = pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.BIN1");
+    if (!bin1) {
         psError(psErrorCodeLast(), false, "Unable to define output.");
         return false;
     }
-    if (!pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.BIN2")) {
+    pmFPAfile *bin2 = pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.BIN2");
+    if (!bin2) {
         psError(psErrorCodeLast(), false, "Unable to define output.");
         return false;
+    }
+    if (data->masksName) {
+      pmFPAfile *mask1 = pmFPAfileDefineOutput(data->config, bin1->fpa , "PPSKYCELL.BIN1.MASK");
+      if (!mask1) {
+        psError(psErrorCodeLast(), false, "Unable to define output.");
+        return false;
+      }
+      mask1->save = true;
+      pmFPAfile *mask2 = pmFPAfileDefineOutput(data->config, bin2->fpa , "PPSKYCELL.BIN2.MASK");
+      if (!mask2) {
+        psError(psErrorCodeLast(), false, "Unable to define output.");
+        return false;
+      }
+      mask2->save = true;
     }
 

trunk/ppSkycell/src/ppSkycellLoop.c

@@ -15 +15 @@
                          )
 {
-    base->x0 = PS_MIN(base->x0, compare->x0);
-    base->x1 = PS_MAX(base->x1, compare->x1);
+    base->x0 = PS_MAX(base->x0, compare->x0);
+    base->x1 = PS_MIN(base->x1, compare->x1);
     base->y0 = PS_MIN(base->y0, compare->y0);
     base->y1 = PS_MAX(base->y1, compare->y1);
@@ -27 +27 @@
 {
     psPlane *fromCoords = psPlaneAlloc(), *toCoords = psPlaneAlloc(); // Coordinates for transforms
-    psRegion *region = psRegionAlloc(INFINITY, -INFINITY, INFINITY, -INFINITY); // Region for skycell
+    psRegion *region = psRegionAlloc(-INFINITY, INFINITY, INFINITY, -INFINITY); // Region for skycell
 
 // Limit the region using the nominated point (X,Y)
@@ -38 +38 @@
     toCoords->x += wcs->crpix1; \
     toCoords->y += wcs->crpix2; \
-    region->x0 = PS_MIN(region->x0, toCoords->x); \
-    region->x1 = PS_MAX(region->x1, toCoords->x); \
+    region->x0 = PS_MAX(region->x0, toCoords->x); \
+    region->x1 = PS_MIN(region->x1, toCoords->x); \
     region->y0 = PS_MIN(region->y0, toCoords->y); \
     region->y1 = PS_MAX(region->y1, toCoords->y);
@@ -132 +132 @@
     psVector *target = psVectorAlloc(data->numInputs, PS_TYPE_S32); // Target for each input
     psArray *imageRegions = psArrayAlloc(data->numInputs); // Region for image
-
+    psArray *regionHDUs = psArrayAlloc(data->numInputs);
+    // Determine which projection cells we have to deal with.
     for (int i = 0; i < data->numInputs; i++) {
         pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", i); // File to examine
@@ -178 +179 @@
             target->data.S32[i] = numProj;
             psTrace("ppSkycell", 3, "Image %d uses new projection\n", i);
-
+            psArrayAdd(regionHDUs,1,hdu);
             numProj++;
         }
@@ -184 +185 @@
 
     pmFPAfileActivate(data->config->files, false, NULL);
-
+    // Loop over projections
     for (int i = 0; i < numProj; i++) {
         psRegion *projRegion = projRegions->data[i]; // Region for skycell projection
         psTrace("ppSkycell", 2, "Projection %d: [%.0f:%.0f,%.0f:%.0f]\n",
                 i, projRegion->x0, projRegion->x1, projRegion->y0, projRegion->y1);
-        int xSize = projRegion->x1 - projRegion->x0 + 1; // Size of unbinned image
+        int xSize = -projRegion->x1 + projRegion->x0 + 1; // Size of unbinned image
         int ySize = projRegion->y1 - projRegion->y0 + 1; // Size of unbinned image
         // Size of binned image 1
@@ -198 +199 @@
         psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); // Binned image
         psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); // Binned image
-        psImageInit(image1, 0);
-        psImageInit(image2, 0);
-
+        psImageInit(image1,NAN);
+        psImageInit(image2,NAN);
+        
         psImage *mask1 = NULL, *mask2 = NULL; // Binned masks
         if (data->masksName) {
@@ -208 +209 @@
             psImageInit(mask2, 0xFF);
         }
-
+        pmHDU *projhdu = NULL;
+        int modify_wcs1 = 1;
+        int modify_wcs2 = 1;
+        // Loop over inputs to this projection.
         for (int j = 0; j < data->numInputs; j++) {
             if (target->data.S32[j] != i) {
                 continue;
             }
-
             pmFPAfileActivateSingle(data->config->files, true, "PPSKYCELL.IMAGE", j);
             if (data->masksName) {
@@ -227 +230 @@
 
             pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", j);
+            if (!projhdu) {
+              projhdu = file->fpa->hdu;
+              //              psMetadataCopy(projhdu->header,file->fpa->hdu->header);
+            }
+#if 0
+            else {
+              if ((psMetadataLookupF32(NULL,file->fpa->hdu->header,"CRPIX1") >=
+                   psMetadataLookupF32(NULL,projhdu->header,"CRPIX1"))&&
+                  (psMetadataLookupF32(NULL,file->fpa->hdu->header,"CRPIX2") >=
+                   psMetadataLookupF32(NULL,projhdu->header,"CRPIX1"))) {
+                projhdu = file->fpa->hdu;
+                psMetadataCopy(projhdu->header,file->fpa->hdu->header);
+              }
+            }
+#endif    
+            
             pmReadout *inRO = pmFPAviewThisReadout(view, file->fpa); // Readout with input
             psFree(view);
 
-            // Flip images; no idea why this has to be done, but apparently it does
-            {
-                psImage *rot = psImageRotate(NULL, inRO->image, M_PI, NAN, PS_INTERPOLATE_BILINEAR);
-                psFree(inRO->image);
-                inRO->image = rot;
-            }
-            if (inRO->mask) {
-                psImage *rot = psImageRotate(NULL, inRO->mask, M_PI, 0, PS_INTERPOLATE_FLAT);
-                psFree(inRO->mask);
-                inRO->mask = rot;
-            }
-
+            //      data->maskVal = 0xffff;
             pmReadout *bin1RO = pmReadoutAlloc(NULL), *bin2RO = pmReadoutAlloc(NULL); // Binned readouts
             if (!pmReadoutRebin(bin1RO, inRO, data->maskVal, data->bin1, data->bin1)) {
@@ -256 +264 @@
             psRegion *imageRegion = imageRegions->data[j]; // Region for image
             // Offsets for image on skycell
-            int xOffset1 = (imageRegion->x0 - projRegion->x0) / (float)data->bin1;
-            int yOffset1 = (imageRegion->y0 - projRegion->y0) / (float)data->bin1;
+            int xOffset1 = (-imageRegion->x0 + projRegion->x0) / (float)data->bin1;
+            int yOffset1 = (-imageRegion->y1 + projRegion->y1) / (float)data->bin1;
             int xOffset2 = xOffset1 / (float)data->bin2, yOffset2 = yOffset1 / (float)data->bin2;
-
             // XXX Completely neglecting rotations
             // The skycells are divided up neatly with them all having the same orientation
-            psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "=");
-            psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "=");
+            psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "E");
+            psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "E");
             if (data->masksName) {
-                psImageOverlaySection(mask1, bin1RO->mask, xOffset1, yOffset1, "=");
-                psImageOverlaySection(mask2, bin2RO->mask, xOffset2, yOffset2, "=");
+                psImageOverlaySection(mask1, bin1RO->mask, xOffset1, yOffset1, "M");
+                psImageOverlaySection(mask2, bin2RO->mask, xOffset2, yOffset2, "M");
             }
 
@@ -287 +294 @@
           pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN1");
           pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN2");
+          if (data->masksName) {
+            pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN1.MASK");
+            pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN2.MASK");
+          }
         }
         
@@ -313 +324 @@
 
         if (data->doFits) {
+
+          pmFPAfile *fits1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1", 0);
+          // Copy header from projection root hdu
+          fits1->fpa->hdu = pmHDUAlloc(NULL);
+          fits1->fpa->hdu->header = psMetadataAlloc();
+          psMetadataCopy(fits1->fpa->hdu->header,projhdu->header);
+
+          // Change wcs here
+#define WCS_DEBUG 0
+          if (modify_wcs1) {
+            pmAstromWCS *WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header);
+            double cd1f = 1.0 * data->bin1;
+            double cd2f = 1.0 * data->bin1;
+#if WCS_DEBUG
+            fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2,
+                    cd1f,cd2f,WCS->cdelt1,WCS->cdelt2,
+                    WCS->crpix1,WCS->crpix2
+                    );
+#endif
+            WCS->cdelt1 *= cd1f;
+            WCS->cdelt2 *= cd2f;
+            WCS->crpix1 = WCS->crpix1 / cd1f;
+            WCS->crpix2 = WCS->crpix2 / cd2f;
+#if WCS_DEBUG
+            fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g) %d %d %d %d\n",data->bin1,data->bin2,
+                    cd1f,cd2f,WCS->cdelt1,WCS->cdelt2,
+                    WCS->crpix1,WCS->crpix2,
+                    WCS->trans->x->nX,              WCS->trans->x->nY,
+                    WCS->trans->y->nX,              WCS->trans->y->nY
+                    );
+#endif
+            for (int q = 0; q <= WCS->trans->x->nX; q++) {
+              for (int r = 0; r <= WCS->trans->x->nY; r++) {
+                WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
+#if WCS_DEBUG
+                fprintf(stderr,"PC: %d %d %g %g\n",
+                        q,r,pow(cd1f,q) * pow(cd2f,r),
+                        WCS->trans->x->coeff[q][r]);
+#endif
+              }
+            }
+            for (int q = 0; q <= WCS->trans->y->nX; q++) {
+              for (int r = 0; r <= WCS->trans->y->nY; r++) {
+                WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
+#if WCS_DEBUG
+                fprintf(stderr,"PC: %d %d %g %g\n",
+                         q,r,pow(cd1f,q) * pow(cd2f,r),
+                         WCS->trans->y->coeff[q][r]);
+#endif
+              }
+            }
+            pmAstromWCStoHeader (fits1->fpa->hdu->header,WCS);
+            WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header);
+#if WCS_DEBUG
+            fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2,
+                    cd1f,cd2f,WCS->cdelt1,WCS->cdelt2,
+                    WCS->crpix1,WCS->crpix2
+                    );
+#endif
+            modify_wcs1 = 0;
+          }
+
+          
+          pmChip *Fchip1 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN1");
+          psMetadataAddS32(Fchip1->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1);
+          psMetadataAddS32(Fchip1->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1);
+
           pmCell *Fcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1"); // Rebinned cell 1
-          pmCell *Fcell2 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN2"); // Rebinned cell 2
-
-          // This is a hack to get a functioning header created so the fits images can be written out.
+/*        // This is a hack to get a functioning header created so the fits images can be written out. */
           psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1);
           psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1);
           psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1);
+
+          pmReadout *Fro1 = pmReadoutAlloc(Fcell1);
+          Fro1->image = image1;
+          Fro1->mask = mask1;
+          Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true;
+          
+          fits1->save = true;
+          fits1->fileIndex = i;
+
+          // Do the mask if we need to
+          if (data->masksName) {
+/*          pmFPAfile *Mask1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1.MASK", 0); */
+/*          //      Mask1->fpa->hdu = pmHDUAlloc(NULL); */
+/*          //      Mask1->fpa->hdu->header = psMetadataAlloc(); */
+/*          //      psMetadataCopy(Mask1->fpa->hdu->header,fits1->fpa->hdu->header); */
+/*          pmChip *Mchip1 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN1.MASK"); */
+/*          psMetadataAddS32(Mchip1->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); */
+/*          psMetadataAddS32(Mchip1->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); */
+            
+/*          pmCell *Mcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1.MASK"); // Rebinned cell 1 */
+/*          /\*           // This is a hack to get a functioning header created so the fits images can be written out. *\/ */
+/*          psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1); */
+/*          psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); */
+/*          psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); */
+            
+/*          pmReadout *Mro1 = pmReadoutAlloc(Mcell1); */
+/*          Mro1->image = image1; */
+/*          Mro1->mask = mask1; */
+/*          Mro1->data_exists = Mcell1->data_exists = Mcell1->parent->data_exists = true; */
+            
+/*          Mask1->save = true; */
+/*          Mask1->fileIndex = i; */
+          }
+            
+            
+          
+          // Repeat with second binned image
+          pmFPAfile *fits2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN2", 0);
+          fits2->fpa->hdu = pmHDUAlloc(NULL);
+          fits2->fpa->hdu->header = psMetadataAlloc();
+          psMetadataCopy(fits2->fpa->hdu->header,projhdu->header);
+
+          if (modify_wcs2) {
+            pmAstromWCS *WCS = pmAstromWCSfromHeader(fits2->fpa->hdu->header);
+            double cd1f = 1.0 * data->bin2 * data->bin1;
+            double cd2f = 1.0 * data->bin2 * data->bin1;
+            
+            WCS->cdelt1 *= cd1f;
+            WCS->cdelt2 *= cd2f;
+            WCS->crpix1 = WCS->crpix1 / cd1f;
+            WCS->crpix2 = WCS->crpix2 / cd2f;
+            
+            for (int q = 0; q < WCS->trans->x->nX; q++) {
+              for (int r = 0; r < WCS->trans->x->nY; r++) {
+                WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
+              }
+            }
+            for (int q = 0; q < WCS->trans->y->nX; q++) {
+              for (int r = 0; r < WCS->trans->y->nY; r++) {
+                WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
+              }
+            }
+            pmAstromWCStoHeader (fits2->fpa->hdu->header,WCS);
+            modify_wcs2 = 0;
+          }
+          
+          pmChip *Fchip2 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN2");
+          psMetadataAddS32(Fchip2->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1);
+          psMetadataAddS32(Fchip2->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1);
+
+          
+          pmCell *Fcell2 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN2"); // Rebinned cell 2
           psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1);
           psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1);
           psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1);
-
-          psMetadataAddS32(Fcell1->parent->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1);
-          psMetadataAddS32(Fcell1->parent->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1);
-          psMetadataAddS32(Fcell2->parent->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1);
-          psMetadataAddS32(Fcell2->parent->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1);
-
-          pmReadout *Fro1 = pmReadoutAlloc(Fcell1), *Fro2 = pmReadoutAlloc(Fcell2); // Binned readouts
-          
-          Fro1->image = image1;
+          
+          pmReadout *Fro2 = pmReadoutAlloc(Fcell2); 
           Fro2->image = image2;
-          Fro1->mask = mask1;
           Fro2->mask = mask2;
-          
-          Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true;
           Fro2->data_exists = Fcell2->data_exists = Fcell2->parent->data_exists = true;
-          
-          pmFPAfile *fits1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1", 0);
-          fits1->save = true;
-          fits1->fileIndex = i;
-          pmFPAfile *fits2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN2", 0);
+
           fits2->save = true;
           fits2->fileIndex = i;
+
+          
         }
 

trunk/ppStack/src/Makefile.am

@@ -1 @@
-bin_PROGRAMS = ppStack
+bin_PROGRAMS = ppStack ppStackMedian
 
 if HAVE_SVNVERSION
@@ -25 +25 @@
 ppStack_LDFLAGS = $(PSLIB_LIBS)   $(PSMODULE_LIBS)   $(PSPHOT_LIBS)   $(PPSTATS_LIBS)   $(PPSTACK_LIBS)
 
+ppStackMedian_CFLAGS    = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) $(PPSTATS_CFLAGS) $(PPSTACK_CFLAGS)
+ppStackMedian_LDFLAGS = $(PSLIB_LIBS)   $(PSMODULE_LIBS)   $(PSPHOT_LIBS)   $(PPSTATS_LIBS)   $(PPSTACK_LIBS)
+
 ppStack_SOURCES =               \
         ppStack.c               \
@@ -43 +46 @@
         ppStackCombinePrepare.c \
         ppStackCombineInitial.c \
+        ppStackCombineAlternate.c       \
+        ppStackReject.c         \
+        ppStackCombineFinal.c   \
+        ppStackCleanup.c        \
+        ppStackPhotometry.c     \
+        ppStackFinish.c         \
+        ppStackTarget.c         \
+        ppStackUpdateHeader.c   \
+        ppStackJPEGs.c          \
+        ppStackStats.c          \
+        ppStackErrorCodes.c
+
+ppStackMedian_SOURCES =         \
+        ppStackMedian.c         \
+        ppStackArguments.c      \
+        ppStackCamera.c         \
+        ppStackFiles.c          \
+        ppStackMedianLoop.c             \
+        ppStackPSF.c            \
+        ppStackReadout.c        \
+        ppStackVersion.c        \
+        ppStackMatch.c          \
+        ppStackSources.c        \
+        ppStackThread.c         \
+        ppStackOptions.c        \
+        ppStackSetup.c          \
+        ppStackPrepare.c        \
+        ppStackConvolve.c       \
+        ppStackCombinePrepare.c \
+        ppStackCombineInitial.c \
+        ppStackCombineAlternate.c       \
         ppStackReject.c         \
         ppStackCombineFinal.c   \

trunk/ppStack/src/ppStack.h

@@ -38 +38 @@
     PPSTACK_FILES_STACK,                // Stack files
     PPSTACK_FILES_UNCONV,               // Unconvolved stack files
-    PPSTACK_FILES_PHOT                  // Files for photometry
+    PPSTACK_FILES_PHOT,                 // Files for photometry
+    PPSTACK_FILES_BKG,                  // Files for bkg
+    PPSTACK_FILES_MEDIAN_IN,                // Files for median only stacks.
+    PPSTACK_FILES_MEDIAN_OUT                // Files for median only stacks.
 } ppStackFileList;
 
@@ -111 +114 @@
 bool ppStackReadoutFinalThread(psThreadJob *job // Job to process
     );
+// Perform median stacking for background
+bool ppStackCombineBackground(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config);
+bool ppStackCombineMedian(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config);
 
 // Return software version

trunk/ppStack/src/ppStackArguments.c

@@ -324 +324 @@
     }
 
+    
     psTrace("ppStack", 1, "Done parsing arguments\n");
 

trunk/ppStack/src/ppStackCamera.c

@@ -156 +156 @@
             psString psf = psMetadataLookupStr(&mdok, input, "PSF"); // Name of PSF
             psString sources = psMetadataLookupStr(&mdok, input, "SOURCES"); // Name of sources
-
+            psString bkgmodel = psMetadataLookupStr(&mdok, input, "BKGMODEL"); // Name of warped background model
+            
             pmFPAfile *imageFile = defineFile(config, NULL, "PPSTACK.INPUT",
                                               image, PM_FPA_FILE_IMAGE); // File for image
@@ -215 +216 @@
             }
 
+            // Grab bkgmodel information here
+            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) {
+                psError(psErrorCodeLast(), false,
+                        "Unable to define file from bkgmodel %d (%s)",i,bkgmodel);
+                return(false);
+              }
+            }// End bkgmodel
+
+
+            
+            
             i++;
         }
@@ -460 +478 @@
     jpeg2->save = true;
 
+    // Output background
+    pmFPAfile *outBkg = pmFPAfileDefineOutput(config,NULL,"PPSTACK.OUTPUT.BKGMODEL");
+    if (!outBkg) {
+      psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.BKGMODEL"));
+      return(false);
+    }
+    if (outBkg->type != PM_FPA_FILE_IMAGE) {
+      psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.BKGMODEL is not of type IMAGE");
+      return(false);
+    }
+    outBkg->save = true;
+    pmFPAfile *outBkgRest = pmFPAfileDefineOutput(config,NULL,"PPSTACK.OUTPUT.BKGREST");
+    if (!outBkgRest) {
+      psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.BKGREST"));
+      return(false);
+    }
+    if (outBkgRest->type != PM_FPA_FILE_IMAGE) {
+      psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.BKGREST is not of type IMAGE");
+      return(false);
+    }
+    
     // For photometry, we operate on the chip-mosaicked image
     // we create a copy of the mosaicked image for psphot so we can write out a clean image
@@ -496 +535 @@
     }
 
+    // Define output file here.
+    
     return true;
 }

trunk/ppStack/src/ppStackCleanup.c

@@ -22 +22 @@
     options->outRO = NULL;
 
-    options->expRO->data_exists = false;
-    options->expRO->parent->data_exists = false;
-    options->expRO->parent->parent->data_exists = false;
-    psFree(options->expRO);
-    options->expRO = NULL;
+    if (options->expRO) {
+      options->expRO->data_exists = false;
+      if (options->expRO->parent) {
+        options->expRO->parent->data_exists = false;
+        options->expRO->parent->parent->data_exists = false;
+      }
+      psFree(options->expRO);
+      options->expRO = NULL;
+    }
 
+    if (options->bkgRO) {
+      options->bkgRO->data_exists = false;
+      if (options->bkgRO->parent) {
+        options->bkgRO->parent->data_exists = false;
+        options->bkgRO->parent->parent->data_exists = false;
+      }
+      psFree(options->bkgRO);
+      options->bkgRO = NULL;
+    }
+    
     for (int i = 0; i < options->num; i++) {
         pmCellFreeData(options->cells->data[i]);
@@ -62 +76 @@
 
 bool ppStackCleanup (pmConfig *config, ppStackOptions *options) {
-
     psExit exitValue = ppStackExitCode(PS_EXIT_SUCCESS); // Exit code
-
+    
     // Ensure everything closes
     if (config) {
@@ -72 +85 @@
         ppStackFileActivation(config, PPSTACK_FILES_UNCONV, true);
         ppStackFileActivation(config, PPSTACK_FILES_PHOT, true);
+        ppStackFileActivation(config, PPSTACK_FILES_MEDIAN_IN, true);
+        ppStackFileActivation(config, PPSTACK_FILES_MEDIAN_OUT, true);
         if (!ppStackFilesIterateUp(config)) {
             psError(psErrorCodeLast(), false, "Unable to close files.");

trunk/ppStack/src/ppStackCombinePrepare.c

@@ -1 +1 @@
 #include "ppStack.h"
 
-bool ppStackCombinePrepare(const char *outName, const char *expName,
+bool ppStackCombinePrepare(const char *outName, const char *expName, const char *bkgName,
                            ppStackFileList files, ppStackThreadData *stack,
                            ppStackOptions *options, pmConfig *config)
@@ -27 +27 @@
     options->outRO = pmReadoutAlloc(cell); // Output readout
 
-    pmCell *expCell = pmFPAfileThisCell(config->files, view, expName); // Exposure cell
-    options->expRO = pmReadoutAlloc(expCell); // Output readout
+    if (expName) {
+      pmCell *expCell = pmFPAfileThisCell(config->files, view, expName); // Exposure cell
+      options->expRO = pmReadoutAlloc(expCell); //Output readout
+    }
+/*     else { */
+/*       options->expRO = NULL; */
+/*     } */
+    pmCell *bkgCell;
+    int bkg_r0,bkg_c0;
+    int bkg_nC,bkg_nR;
+    if (bkgName) {
+      ppStackFileActivationSingle(config, PPSTACK_FILES_BKG, true, 0);
+      pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", 0);
+      pmFPAview *view = ppStackFilesIterateDown(config);
+      pmReadout *ro = pmFPAviewThisReadout(view,file->fpa);
+
+      bkg_r0 = ro->image->row0;
+      bkg_c0 = ro->image->col0;
+      bkg_nC = ro->image->numCols;
+      bkg_nR = ro->image->numRows;
+      bkgCell = pmFPAfileThisCell(config->files, view, bkgName); // Bkg cell
+      
+      options->bkgRO = pmReadoutAlloc(bkgCell); // BKG readout
+      //      if (!pmHDUGenerateForFPA(options->bkgRO->parent->parent->parent)) {
+      options->bkgRO->parent->parent->parent->hdu = pmHDUAlloc(NULL);
+      if (!options->bkgRO->parent->parent->parent->hdu) {
+        fprintf(stderr,"failed to generate a HDU for this thing.\n");
+      }
+      options->bkgRO->parent->parent->parent->hdu->header = psMetadataCopy(options->bkgRO->parent->parent->parent->hdu->header,
+                                                                           ro->parent->parent->parent->hdu->header);
+
+      options->bkgRO->parent->concepts = psMetadataCopy(options->bkgRO->parent->concepts,
+                                                        ro->parent->concepts);
+      options->bkgRO->parent->parent->concepts = psMetadataCopy(options->bkgRO->parent->parent->concepts,
+                                                                ro->parent->parent->concepts);
+      options->bkgRO->parent->parent->parent->concepts = psMetadataCopy(options->bkgRO->parent->parent->parent->concepts,
+                                                                        ro->parent->parent->parent->concepts);
+
+    }
+    else {
+      options->bkgRO = NULL;
+    }
 
     psFree(view);
@@ -42 +82 @@
     }
 
-    if (!pmReadoutStackDefineOutput(options->expRO, col0, row0, numCols, numRows, true, true, 0)) {
+    if (expName) {
+      if (!pmReadoutStackDefineOutput(options->expRO, col0, row0, numCols, numRows, true, true, 0)) {
         psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output.");
         return false;
+      }
+    }
+
+    if (bkgName) {
+      if (!pmReadoutStackDefineOutput(options->bkgRO, bkg_c0, bkg_r0, bkg_nC, bkg_nR, false, false, 0)) {
+        psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output.");
+        return false;
+      }
+      
+
+      
     }
 

trunk/ppStack/src/ppStackFiles.c

@@ -14 +14 @@
 /// Files required for the convolution
 static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL };
+
+/// Files required for the background
+static char *filesBkg[] = { "PPSTACK.INPUT.BKGMODEL", NULL };
+
+/// Files required for median only stacking
+static char *filesMedianIn[] =  { "PPSTACK.INPUT",                                                           
+                              NULL };
+/// Files required for median only stacking
+static char *filesMedianOut[] =  { "PPSTACK.OUTPUT",                                                         
+                              NULL };
 
 /// Regular (convolved) stack files
@@ -19 +29 @@
                               "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.EXPNUM", "PPSTACK.OUTPUT.EXPWT",
                               "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2",
+                              "PPSTACK.OUTPUT.BKGMODEL",
                               NULL };
 /// Unconvolved stack files
@@ -41 +52 @@
       case PPSTACK_FILES_UNCONV:   return filesUnconv;
       case PPSTACK_FILES_PHOT:     return filesPhot;
+    case PPSTACK_FILES_BKG:        return filesBkg;
+    case PPSTACK_FILES_MEDIAN_IN:        return filesMedianIn;
+    case PPSTACK_FILES_MEDIAN_OUT:        return filesMedianOut;
       default:
         psAbort("Unrecognised file list: %x", list);

trunk/ppStack/src/ppStackLoop.c

@@ -63 +63 @@
 
     // Prepare for combination
-    if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", PPSTACK_FILES_STACK, stack, options, config)) {
+    if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.BKGMODEL", PPSTACK_FILES_STACK, stack, options, config)) {
         psError(psErrorCodeLast(), false, "Unable to prepare for combination.");
         psFree(stack);
@@ -84 +84 @@
         pmCellFreeData(options->cells->data[i]);
     }
-    psFree(stack);
+    //    psFree(stack);
 
     // Pixel rejection
@@ -142 +142 @@
     }
 
+    // Generate median background stack here.
+    if (!ppStackCombineBackground(stack, options, config)) {
+      psError(psErrorCodeLast(), false, "Unable to generate median of background images.");
+      psFree(stack);
+      return false;
+    }
+    ppStackFileActivation(config,PPSTACK_FILES_BKG ,false);    
     // Photometry
     psTrace("ppStack", 1, "Photometering stacked image....\n");
@@ -171 +178 @@
         return false;
     }
-    psFree(stack);
+    //    psFree(stack);
     psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));
     ppStackMemDump("cleanup");
@@ -186 +193 @@
 
         // Prepare for combination
-        if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", PPSTACK_FILES_UNCONV,
+        if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", NULL, PPSTACK_FILES_UNCONV,
                                    stack, options, config)) {
             psError(psErrorCodeLast(), false, "Unable to prepare for combination.");

trunk/ppStack/src/ppStackLoop.h

@@ -4 +4 @@
 // Loop over the inputs, doing the combination
 bool ppStackLoop(
+    pmConfig *config,                    // Configuration
+    ppStackOptions *options             // Options for stacking
+    );
+// Median only loop.
+bool ppStackMedianLoop(
     pmConfig *config,                    // Configuration
     ppStackOptions *options             // Options for stacking
@@ -32 +37 @@
     const char *outName,                // Name of output file
     const char *expName,                // Name of exposure file
+    const char *bkgName,                // Name of background file
     ppStackFileList files,              // Files of interest
     ppStackThreadData *stack,           // Stack

trunk/ppStack/src/ppStackOptions.h

@@ -48 +48 @@
     // Rejection
     psArray *rejected;                  // Rejected pixels
+    // Background
+    pmReadout *bkgRO;                   // Output background readout
+    psArray   *bkgImages;               // Input background images
 } ppStackOptions;
 

trunk/ppStack/src/ppStackPrepare.c

@@ -132 +132 @@
     psVectorInit(options->inputMask, 0);
     options->exposures = psVectorAlloc(options->num, PS_TYPE_F32);
+
     psVectorInit(options->exposures, NAN);
 
@@ -150 +151 @@
 
         options->exposures->data.F32[i] = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE");
+        if ((options->exposures->data.F32[i] == 0)||
+            (!(isfinite(options->exposures->data.F32[i])))){
+          options->exposures->data.F32[i] = psMetadataLookupF32(NULL,recipe,"DEFAULT.EXPTIME");
+        }
         options->sumExposure += options->exposures->data.F32[i];
 

trunk/ppStack/src/ppStackReadout.c

@@ -194 +194 @@
 
 
-
 bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, const psArray *readouts,
                          const psVector *mask, const psArray *rejected, const psVector *weightings,

trunk/ppStack/src/ppStackSetup.c

@@ -89 +89 @@
     options->origMasks = psArrayAlloc(num);
     options->origVariances = psArrayAlloc(num);
+    options->bkgImages = psArrayAlloc(num);
     pmFPAview *view = pmFPAviewAlloc(0);
+    int nullMasks = 0;
+    int nullVariances = 0;
     for (int i = 0; i < num; i++) {
         {
@@ -97 +100 @@
         {
             // We want the convolved mask, since that defines the area that has been tested for outliers
+          if (options->convolve) {
             options->origMasks->data[i] = psMemIncrRefCounter(options->convMasks->data[i]);
+          }
+          else {
+            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.MASK", i);
+            options->origMasks->data[i] = pmFPAfileName(file, view, config);
+          }
+          if (!(options->origMasks->data[i])) {
+            nullMasks++;
+          }
         }
         {
             pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.VARIANCE", i);
             options->origVariances->data[i] = pmFPAfileName(file, view, config);
+            if (!(options->origVariances->data[i])) {
+              nullVariances++;
+            }
         }
+/*      { */
+/*        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", i); */
+/*        options->bkgImages->data[i] = pmFPAfileName(file, view, config); */
+/*      } */
     }
+    if (nullMasks == num) {
+      psFree(options->origMasks);
+    }
+    if (nullVariances == num) {
+      psFree(options->origVariances);
+    }
+    
     psFree(view);
 

trunk/ppStack/src/ppStackThread.c

@@ -40 +40 @@
     psFree(stack->maskFits);
     psFree(stack->varianceFits);
+    psFree(stack->bkgFits);
     return;
 }
@@ -67 +68 @@
         PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL);
     }
-
+    
     ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return
     psMemSetDeallocator(stack, (psFreeFunc)stackThreadDataFree);
@@ -77 +78 @@
     stack->maskFits   = psArrayAlloc(numInputs);
     stack->varianceFits = psArrayAlloc(numInputs);
+    stack->bkgFits    = psArrayAlloc(numInputs);
     for (int i = 0; i < numInputs; i++) {
         if (!cells->data[i]) {
@@ -95 +97 @@
             psFree(resolved); \
         }
-
+                
         IMAGE_OPEN(imageNames, stack->imageFits, i);
-        IMAGE_OPEN(maskNames, stack->maskFits, i);
-        IMAGE_OPEN(varianceNames, stack->varianceFits, i);
+        if (maskNames) {
+          IMAGE_OPEN(maskNames, stack->maskFits, i);
+        }
+        if (varianceNames) {
+          IMAGE_OPEN(varianceNames, stack->varianceFits, i);
+        }
     }
 

trunk/ppStack/src/ppStackThread.h

@@ -25 +25 @@
     psArray *maskFits;                  // FITS file pointers for masks
     psArray *varianceFits;              // FITS file pointers for variances
+    psArray *bkgFits;                   // FITS file pointers for background models
 } ppStackThreadData;
 

trunk/ppStack/src/ppStackUpdateHeader.c

@@ -6 +6 @@
 
     pmReadout *outRO = options->outRO;                                      // Output readout
+    pmReadout *expRO = options->expRO;
 
     // Propagate WCS
     bool wcsDone = false;           // Have we done the WCS?
     for (int i = 0; i < options->num && !wcsDone; i++) {
-        if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+      if ((options->inputMask)&&(options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i])) {
             continue;
         }
@@ -123 +124 @@
         snprintf (field, 64, "AIR_%04d", i);
         psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image airmass", value);
-    }   
+    }
+
+    // Copy information into expRO, because it should be there too.
+    if ((expRO)&&(expRO->parent)) {
+      pmHDU *expROhdu = pmHDUFromCell(expRO->parent);
+      if (!expROhdu->header) {
+        expROhdu->header = psMetadataAlloc();
+      }
+      expRO->parent->parent->parent->hdu->header = psMetadataCopy(expRO->parent->parent->parent->hdu->header,
+                                                                  outRO->parent->parent->parent->hdu->header);
+      
+      expRO->parent->concepts = psMetadataCopy(expRO->parent->concepts,
+                                               outRO->parent->concepts);
+      expRO->parent->parent->concepts = psMetadataCopy(expRO->parent->parent->concepts,
+                                                       outRO->parent->concepts);
+      expRO->parent->parent->parent->concepts = psMetadataCopy(expRO->parent->parent->parent->concepts,
+                                                               outRO->parent->parent->parent->concepts);
+    }
+    
     return true;
 }

trunk/ppViz/src/ppVizPSF/ppVizPSFCamera.c

@@ -19 +19 @@
     files->data[0] = psStringCopy(name);
     if (psMetadataLookup(config->arguments, file)) {
+
         psMetadataRemoveKey(config->arguments, file);
     }
@@ -42 +43 @@
         fileArguments("SOURCES", data->sourcesName, "Input sources", data->config);
         pmFPAfile *srcs = pmFPAfileBindFromArgs(&status, psf, data->config,
-                                                "PSPHOT.INPUT.CMF", "SOURCES"); // File
+                                                "PSPHOT.INPUT.CMF", "SOURCES"); // File
+        fprintf(stderr,"%ld %d\n",(long) srcs, status);
         if (!status || !srcs) {
             psError(PS_ERR_IO, false, "Failed to build file from PSPHOT.INPUT.CMF");

trunk/ppViz/src/ppVizPSF/ppVizPSFLoop.c

@@ -20 +20 @@
         return NULL;
     }
-
+    fprintf(stderr,"Woo!\n");
     pmChip *chip;                       // Chip from FPA
     while ((chip = pmFPAviewNextChip(view, psfFile->fpa, 1))) {
@@ -56 +56 @@
                 psWarning("More than one readout present for chip %d, cell %d", view->chip, view->cell);
             }
-
+            fprintf(stderr,"Woo!\n");
             pmReadout *readout;         // Readout from cell
             while ((readout = pmFPAviewNextReadout(view, psfFile->fpa, 1))) {
+              fprintf(stderr,"Woo?\n");
                 if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
                     psError(PS_ERR_UNKNOWN, false, "Error loading data from files.");
                     return false;
                 }
-                if (!readout->data_exists) {
-                    continue;
-                }
-
+                fprintf(stderr,"Woo2?\n");
                 pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF");              // PSF
                 assert(psf);
@@ -72 +70 @@
                 bool mdok;              // Status of MD lookup
                 psArray *sources = psMetadataLookupPtr(&mdok, readout->analysis, "PSPHOT.SOURCES"); // Sources
+
+                if (sources && !readout->data_exists) {  // This fails if -sources not specified
+                    continue;
+                }
+                readout->data_exists = true;
+                fprintf(stderr,"Woo! %ld\n", sources ? sources->n : -1);
+
+                fprintf(stderr,"%d\n",mdok);
                 int numCols = 0, numRows = 0;              // Size of image
                 psVector *xOffset = NULL, *yOffset = NULL; // Offset from source to true position
@@ -80 +86 @@
                     psLogMsg("ppVizPSF", PS_LOG_INFO, "Generating %dx%d image", numCols, numRows);
                 }
-                if (sources) {
+                if (sources && !sources->n) {
                     psMemIncrRefCounter(sources);
                     psLogMsg("ppVizPSF", PS_LOG_INFO, "Using %ld input sources from CMF file", sources->n);
@@ -89 +95 @@
 
                 if (data->fakeNum > 0 && isfinite(data->fakeMag)) {
+                  fprintf(stderr,"Here! fakes\n");
                     long numOld = 0; // Old number of sources
                     long numNew = -1; // New number of sources
@@ -118 +125 @@
                 }
                 if (!sources && !data->input) {
+                  fprintf(stderr,"Here! default\n");
                     // Generate fake image with only a single realisation of the PSF
                     sources = psArrayAlloc(1);
@@ -151 +159 @@
                 psFree(psf->residuals);
                 psf->residuals = NULL;
-
+                fprintf(stderr,"still here");
                 if (sources) {
+                  fprintf(stderr,"Here! sources?\n");
                     if (!pmReadoutFakeFromSources(readout, numCols, numRows, sources, 0, xOffset, yOffset,
                                                   psf, data->minFlux, 0, false, true)) {
@@ -159 +168 @@
                     }
                 } else if (data->input) {
+
                     psVector *x = data->input->data[0]; // x coordinates
                     psVector *y = data->input->data[1]; // y coordinates
                     psVector *mag = data->input->data[2]; // Magnitudes
+                    fprintf(stderr,"Here! input? %d %d %g %d %ld\n",numCols,numRows,data->minFlux,1,x->n);
                     if (!pmReadoutFakeFromVectors(readout, numCols, numRows, x, y, mag, xOffset, yOffset,
                                                   psf, data->minFlux, 0, false, true)) {
@@ -168 +179 @@
                     }
                 }
-
+                
                 psFree(sources);
                 psFree(xOffset);

trunk/psLib/src/imageops/psImageInterpolate.c

@@ -31 +31 @@
 #include "psImageConvolve.h"
 
-# define IS_BILIN_SEPARABLE 0
+# define IS_BILIN_SEPARABLE 1
 
 #include "psImageInterpolate.h"
@@ -189 +189 @@
     switch (mode) {
       case PS_INTERPOLATE_FLAT:
+      case PS_INTERPOLATE_BILINEAR_SIMPLE:
         // Nothing to pre-compute
         break;
@@ -287 +288 @@
 # endif
       case PS_INTERPOLATE_BIQUADRATIC:
+      case PS_INTERPOLATE_BILINEAR_SIMPLE:
         // 2D kernel, would cost too much memory to pre-calculate
         break;
+
       default:
         psAbort("Unsupported interpolation mode: %x", mode);
@@ -937 +940 @@
 
 
+psImageInterpolateStatus interpolateJustWork(double *imageValue, double *varianceValue,
+                                             psImageMaskType *maskValue, float x, float y,
+                                             const psImageInterpolation *interp) {
+  float u1,u2;
+  int xl,xh;
+  int yl,yh;
+  const psImage *image = interp->image; // Image to interpolate
+
+  xl = floor(x);
+  if (xl < 0) { xl = 0; }
+  if (xl >= image->numCols) {xl = image->numCols - 1; }
+  xh = xl + 1;
+  if (xh >= image->numCols) {xh = image->numCols - 1; }
+
+  yl = floor(y);
+  if (yl < 0) { yl = 0; }
+  if (yl >= image->numRows) {yl = image->numRows - 1; }
+  yh = yl + 1;
+  if (yh >= image->numRows) {yh = image->numRows - 1; }
+
+  if (imageValue && image) {
+    if (image->data.F32[yl][xl] == image->data.F32[yl][xh]) {
+      u1 = image->data.F32[yl][xh];
+    }
+    else {
+      u1 = (xh - x) * image->data.F32[yl][xl] + (x - xl) * image->data.F32[yl][xh];
+    }
+    if (image->data.F32[yh][xl] == image->data.F32[yh][xh]) {
+      u2 = image->data.F32[yh][xh];
+    }
+    else {
+      u2 = (xh - x) * image->data.F32[yh][xl] + (x - xl) * image->data.F32[yh][xh];
+    }
+    if (u1 == u2) {
+      *imageValue = u1;
+    }
+    else {
+      *imageValue = (yh - y) * u1 + (y - yl) * u2;
+    }
+    if ((floor(x) >= image->numCols)||
+        (floor(y) >= image->numRows)) {
+      *imageValue = NAN;
+    }
+  }
+#if (0)
+  if (*imageValue == 0.0) {
+    fprintf(stderr,"IJK: Zero!: %g %g [%d %d %d %d] %g %g (%g %g %g %g)\n",
+            x,y,xl,xh,yl,yh,u1,u2,
+            image->data.F32[yl][xl],image->data.F32[yl][xh],image->data.F32[yh][xl],image->data.F32[yh][xh]
+            );
+  }
+#endif
+  return PS_INTERPOLATE_STATUS_GOOD;
+}
+  
 
 psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue,
@@ -969 +1027 @@
       case PS_INTERPOLATE_BIQUADRATIC:
         return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp);
+      case PS_INTERPOLATE_BILINEAR_SIMPLE:
+        return interpolateJustWork(imageValue, varianceValue, maskValue, x, y, interp);
       case PS_INTERPOLATE_BILINEAR:
 # if (!IS_BILIN_SEPARABLE)
-        return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp);
+        return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp);
 # endif
       case PS_INTERPOLATE_GAUSS:
@@ -1019 +1079 @@
     switch (mode) {
       case PS_INTERPOLATE_FLAT:
+      case PS_INTERPOLATE_BILINEAR_SIMPLE:
         // No smearing by design
         return 1.0;
@@ -1083 +1144 @@
     if (!strcasecmp(name, "FLAT"))     return PS_INTERPOLATE_FLAT;
     if (!strcasecmp(name, "BILINEAR")) return PS_INTERPOLATE_BILINEAR;
+    if (!strcasecmp(name, "SIMPLEBILINEAR")) return PS_INTERPOLATE_BILINEAR_SIMPLE;
     if (!strcasecmp(name, "BIQUADRATIC")) return PS_INTERPOLATE_BIQUADRATIC;
     if (!strcasecmp(name, "GAUSS"))    return PS_INTERPOLATE_GAUSS;

trunk/psLib/src/imageops/psImageInterpolate.h

@@ -33 +33 @@
     PS_INTERPOLATE_LANCZOS3,            ///< Sinc interpolation with 6x6 pixel kernel
     PS_INTERPOLATE_LANCZOS4,            ///< Sinc interpolation with 8x8 pixel kernel
+    PS_INTERPOLATE_BILINEAR_SIMPLE,     ///< Simple manual bilinear interpolation
 } psImageInterpolateMode;
 

trunk/psLib/src/imageops/psImagePixelManip.c

@@ -216 +216 @@
     }
 
-
+#define psImageOverlayLoopClean(DATATYPE,OP) {   \
+      for (int row=y0;row<imageRowLimit;row++) {            \
+        ps##DATATYPE* imageRow = image->data.DATATYPE[row];        \
+        ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-y0]; \
+        for (int col=x0;col<imageColLimit;col++) {                 \
+          if (!isfinite(imageRow[col])) {                          \
+            imageRow[col] OP overlayRow[col-x0];                   \
+          }                                                        \
+          else if (!isfinite(overlayRow[col-x0])) {                \
+            imageRow[col] = imageRow[col];                         \
+          }                                                        \
+          else {                                                   \
+            imageRow[col] = overlayRow[col-x0];                    \
+          }                                                             \
+        }                                                               \
+      }                                                                 \
+      pixelsOverlaid += (imageRowLimit - y0) * (imageColLimit - x0);    \
+    }
+
+#define psImageOverlayLoopMask(DATATYPE) {                  \
+      for (int row=y0;row<imageRowLimit;row++) {            \
+        ps##DATATYPE* imageRow = image->data.DATATYPE[row];        \
+        ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-y0]; \
+        for (int col=x0;col<imageColLimit;col++) {                 \
+          if (!isfinite(imageRow[col])) {                          \
+            imageRow[col] = overlayRow[col-x0];            \
+          }                                                        \
+          else if (!isfinite(overlayRow[col-x0])) {                \
+            imageRow[col] = imageRow[col];                         \
+          }                                                        \
+          else {                                                   \
+            imageRow[col] &= overlayRow[col-x0];                   \
+          }                                                             \
+        }                                                               \
+      }                                                                 \
+      pixelsOverlaid += (imageRowLimit - y0) * (imageColLimit - x0);    \
+    }
+        
+    
+    
     #define psImageOverlayLoop(DATATYPE,OP) { \
         for (int row=y0;row<imageRowLimit;row++) { \
             ps##DATATYPE* imageRow = image->data.DATATYPE[row]; \
             ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-y0]; \
-            for (int col=x0;col<imageColLimit;col++) { \
-                imageRow[col] OP overlayRow[col-x0]; \
-            } \
-        } \
+            for (int col=x0;col<imageColLimit;col++) {                 \
+                imageRow[col] OP overlayRow[col-x0];                   \
+              }                                                        \
+        }                                                              \
         pixelsOverlaid += (imageRowLimit - y0) * (imageColLimit - x0); \
-    }
+      }
 
     #define psImageOverlayLoopDivide(DATATYPE,BADVALUE) { \
@@ -273 +312 @@
         psImageOverlayLoop(DATATYPE,*=); \
         break; \
+    case 'E': \
+      psImageOverlayLoopClean(DATATYPE,=); \
+      break;                               \
     case '/': \
         psImageOverlayLoopDivide(DATATYPE,BADVALUE); \
         break; \
     case '=': \
-        psImageOverlaySetLoop(DATATYPE); \
+        psImageOverlaySetLoop(DATATYPE);                \
         break; \
     default: \
@@ -286 +328 @@
     } \
     break;
+    #define psImageOverlayCaseInt(DATATYPE,BADVALUE) \
+case PS_TYPE_##DATATYPE: \
+    switch (*op) { \
+    case '+': \
+        psImageOverlayLoop(DATATYPE,+=); \
+        break; \
+    case '-': \
+        psImageOverlayLoop(DATATYPE,-=); \
+        break; \
+    case '*': \
+        psImageOverlayLoop(DATATYPE,*=); \
+        break; \
+    case 'E': \
+      psImageOverlayLoopClean(DATATYPE,=); \
+      break;                               \
+    case 'M':                              \
+      psImageOverlayLoopMask(DATATYPE);    \
+      break;                               \
+    case '/': \
+        psImageOverlayLoopDivide(DATATYPE,BADVALUE); \
+        break; \
+    case '=': \
+        psImageOverlaySetLoop(DATATYPE);                \
+        break; \
+    default: \
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
+                _("Specified operation, '%s', is not supported."), \
+                op); \
+        return pixelsOverlaid; \
+    } \
+    break;
 
     switch (type) {
-        psImageOverlayCase(U8, 0);
-        psImageOverlayCase(U16,0);
-        psImageOverlayCase(U32,0);       // Not a requirement
-        psImageOverlayCase(U64,0);       // Not a requirement
-        psImageOverlayCase(S8, 0);
-        psImageOverlayCase(S16,0);
-        psImageOverlayCase(S32,0);       // Not a requirement
-        psImageOverlayCase(S64,0);       // Not a requirement
+        psImageOverlayCaseInt(U8, 0);
+        psImageOverlayCaseInt(U16,0);
+        psImageOverlayCaseInt(U32,0);       // Not a requirement
+        psImageOverlayCaseInt(U64,0);       // Not a requirement
+        psImageOverlayCaseInt(S8, 0);
+        psImageOverlayCaseInt(S16,0);
+        psImageOverlayCaseInt(S32,0);       // Not a requirement
+        psImageOverlayCaseInt(S64,0);       // Not a requirement
         psImageOverlayCase(F32,NAN);
         psImageOverlayCase(F64,NAN);

trunk/psModules/src/camera/pmFPABin.c

@@ -43 +43 @@
     }
 
+    int Nbits = (int) (ceil(log(maskVal)/log(2)) + 1);
+    int *bitcounter = malloc(sizeof(int) * Nbits);
+    int pxlcount;
+
     int xLast = numColsIn - 1, yLast = numRowsIn - 1; // Last index
     int yStart = psImageBinningGetFineY(binning, 0); // Starting input y for binning
@@ -55 +59 @@
             float sum = 0.0;            // Sum of pixels
             int numPix = 0;             // Number of pixels
+
+            for (int j = 0; j < Nbits; j++) { // Reset bit counter
+              bitcounter[j] = 0;
+            }
+            pxlcount = 0;
+            
             for (int y = yStart; y < yStop; y++) {
                 for (int x = xStart; x < xStop; x++) {
+                  if (inMask && (inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] != 0)) {
+                      for (int j = 0; j < Nbits; j++) {
+                        psImageMaskType M = (psImageMaskType) pow(2,j);
+                        if (inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & M) {
+                          bitcounter[j]++;
+                        }
+                      }
+                    }
+                  
                     if (inMask && (inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) {
                         continue;
@@ -65 +84 @@
                     sum += inImage->data.F32[y][x];
                     numPix++;
+
+
                 }
             }
-
+            
+            
             // Values to set
             float imageValue;
@@ -79 +101 @@
             }
             outImage->data.F32[yOut][xOut] = imageValue;
-            outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = maskValue;
+            //            outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = maskValue;
+            outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = 0;
+            for (int j = 0; j < Nbits; j++) {
+              if (bitcounter[j] > 0.5 * pxlcount) {
+                outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] |= (int) pow(2,j);
+              }
+            }
             xStart = xStop;
         }

trunk/psModules/src/camera/pmFPAfileIO.c

@@ -76 +76 @@
     while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
         pmFPAfile *file = item->data.V;
-
         switch (place) {
           case PM_FPA_BEFORE:

trunk/psModules/src/detrend/pmDark.c

@@ -353 +353 @@
 
     // retrieve the required parameter vectors
-    psArray *values  = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");
-    psAssert(values, "values not supplied");
+    psArray *in_values  = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");
+    psAssert(in_values, "values not supplied");
     psVector *roMask = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK");
     psAssert(roMask, "roMask not supplied");
-    psVector *orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS");
-    psAssert(orders, "orders not supplied");
-
-    psPolynomialMD *poly = psPolynomialMDAlloc(orders); // Polynomial for fitting
+    psVector *max_orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS");
+    psAssert(max_orders, "orders not supplied");
+
+    psArray *values_set = psArrayAlloc(max_orders->n);
+    psArray *poly_set = psArrayAlloc(max_orders->n);
+    psVector *logL = psVectorAlloc(max_orders->n,PS_TYPE_F64);
+
+    for (int i = 0; i < max_orders->n; i++) {
+      psVector *orders = psVectorAlloc(i+1,PS_TYPE_U8);
+      for (int j = 0; j < orders->n; j++) {
+        orders->data.U8[j] = max_orders->data.U8[j];
+      }
+      poly_set->data[i] =  psPolynomialMDAlloc(orders); // Polynomial for fitting
+      
+      psArray *values = psArrayAlloc(in_values->n);
+      
+      for (int j = 0; j < values->n; j++) {
+        psVector *these_values = psVectorAlloc(i+1,PS_TYPE_F32);
+        psVector *input_values = in_values->data[j];
+
+        for (int k = 0; k < orders->n; k++) {
+          these_values->data.F32[k] = input_values->data.F32[k];
+        }
+        values->data[j] = these_values;
+      }
+      values_set->data[i] = values;
+    }
 
     // retrieve the norm vector, if supplied
@@ -383 +406 @@
     }
 
-    pmDarkVisualInit(values);
+    pmDarkVisualInit(values_set->data[max_orders->n - 1]);
 
     pmReadout *outReadout = output->readouts->data[0];
@@ -423 +446 @@
             }
 
-            if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, 0xff, values, iter, rej)) {
+            int k_best = 0;
+            for (int k = 0; k < max_orders->n; k++) {
+              psPolynomialMD *poly = poly_set->data[k];
+              psArray *values = values_set->data[k];
+              
+              if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, 0xff, values, iter, rej)) {
                 psErrorClear();         // Nothing we can do about it
                 psVectorInit(poly->coeff, NAN);
-            }
-
-            pmDarkVisualPixelFit(pixels, mask);
-            pmDarkVisualPixelModel(poly, values);
-
-            for (int k = 0; k < poly->coeff->n; k++) {
+              }
+
+              pmDarkVisualPixelFit(pixels, mask);
+              pmDarkVisualPixelModel(poly, values);
+
+              // Insert math here to choose optimum model.
+              logL->data.F64[k] = 0.0;
+              psPolynomialMD *polySig = poly_set->data[0];
+              for (int m = 0; m < poly->deviations->n; m++) {
+                logL->data.F64[k] += pow(poly->deviations->data.F32[m] / polySig->stdevFit,2);
+/*              if ((xOut == 20) && (yOut == 256)) { */
+/*                psTrace("psModules.detrend",3,"pmDarkCombine DEV: %d %d: input %d models: Norders: %d logL: %g value: %g\n", */
+/*                        xOut,yOut,m,k,logL->data.F64[k],poly->deviations->data.F32[m]); */
+/*              } */
+              }
+              if (k > 0) {
+                if ( ( logL->data.F64[k - 1] - logL->data.F64[k] ) > 1) { // Hard coded criterion for a ~5% limit with one degree of freedom
+                  k_best = k;
+                }
+              }
+              if ((xOut <= 600) && (yOut <= 600)) {
+                psTrace("psModules.detrend",3,"pmDarkCombine: %d %d: models: Norders: %d logL: %g BestOrders: %d\n",
+                        xOut,yOut,k,logL->data.F64[k],k_best);
+              }
+            }
+            if (k_best > 1) {
+              k_best = 1;
+            }
+/*          k_best = 1; */
+            // Select the polynomial that seems best.
+            psPolynomialMD *poly = poly_set->data[k_best];
+              
+            //            for (int k = 0; k < poly->coeff->n; k++) {
+            for (int k = 0; k < max_orders->n + 1; k++) { // There is one more coefficient than is stored here.
                 pmReadout *ro = output->readouts->data[k]; // Readout of interest
-                ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k];
+                if (k < poly->coeff->n) {
+                  ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k];
+                }
+                else {
+                  ro->image->data.F32[yOut][xOut] = 0.0;
+                }
             }
             counts->data.U16[yOut][xOut] = poly->numFit;

trunk/psModules/src/imcombine/pmStack.c

@@ -1358 +1358 @@
 }
 
+bool pmStackSimpleMedianCombine(
+                                pmReadout *combined,
+                                psArray *input) {
+  int num = input->n;
+  //  int numCols, numRows;
+  int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
+  int xSize, ySize;                   // Size of the output image
+
+  psArray *stack = psArrayAlloc(num); // Stack of readouts  
+  for (int i = 0; i < num; i++) {
+    //    pmStackData *data = input->data[i]; // Stack data for this input
+    pmReadout *ro = input->data[i]; // data->readout;  // Readout of interest
+    if (!ro) {
+      continue;
+    }
+    stack->data[i] = psMemIncrRefCounter(ro);
+  }    
+
+  if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
+                              stack)) {
+    psError(psErrorCodeLast(), false, "Input stack is not valid.");
+    psFree(stack);
+    return false;
+  }
+
+  psVector *pixelData = psVectorAlloc(input->n,PS_TYPE_F32);
+  psStats  *stats     = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+
+  for (int y = minInputRows; y < maxInputRows; y++) {
+    for (int x = minInputCols; x < maxInputCols; x++) {
+      for (int i = 0; i < input->n; i++) {
+        pmReadout *ro  = stack->data[i];
+        psImage *image = ro->image;
+        pixelData->data.F32[i] = image->data.F32[y][x];
+      }
+      if (!psVectorStats(stats,pixelData,NULL,NULL,0)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate median");
+        psFree(stats);
+        psFree(pixelData);
+        psFree(stack);
+        return(false);
+      }
+      combined->image->data.F32[y][x] = stats->robustMedian;
+      combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
+    }
+  }
+  
+  psFree(stats);
+  psFree(pixelData);
+  psFree(stack);
+  return (true);
+}
+
 /// Stack input images
 bool pmStackCombine(

trunk/psModules/src/imcombine/pmStack.h

@@ -41 +41 @@
                               float addVariance ///< Additional variance when rejecting
     );
+/// Stack input images simply
+bool pmStackSimpleMedianCombine(pmReadout *combined, ///< Combined readout (output)
+                                psArray *input       ///< Input array of pmStackData
+                                );
 
 /// Stack input images

trunk/psModules/src/objects/Makefile.am

@@ -111 +111 @@
         pmSourceOutputs.h \
         pmSourceIO.h \
+        pmSourceSatstar.h \ 
         pmSourcePlots.h \
         pmSourceVisual.h \

trunk/psModules/src/objects/pmFootprintCullPeaks.c

@@ -91 +91 @@
 
         // max flux is above threshold for brightest peak
-        pmPeak *maxPeak = fp->peaks->data[0];
+      pmPeak *maxPeak = NULL;
+      for (int i = 0; i < fp->peaks->n; i++) {
+        pmPeak *testPeak = fp->peaks->data[i];
+        float this_peak = useSmoothedImage ? testPeak->smoothFlux : testPeak->rawFlux;
+        
+        if (isfinite(this_peak)) {
+          maxPeak = fp->peaks->data[i];
+          break;
+        }
+      }
+      psAssert(maxPeak,"maxPeak was not set in these peaks");
+      //      = fp->peaks->data[0];
         float maxFlux = useSmoothedImage ? maxPeak->smoothFlux : maxPeak->rawFlux;
 

trunk/pswarp/src/pswarp.c

@@ -46 +46 @@
         goto DIE;
     }
+    if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) {
+      if (!pswarpDefineBackground(config)) {
+        goto DIE;
+      }
+    }
 
     // Open the statistics file
@@ -66 +71 @@
     if (!pswarpLoop(config, stats)) {
         goto DIE;
+    }
+    if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) {
+      if (!pswarpLoopBackground(config, stats)) {
+        fprintf(stderr,"Dying!\n");
+        goto DIE;
+      }
     }
 

trunk/pswarp/src/pswarp.h

@@ -62 +62 @@
     psImage *region;
 
+    /** values which are needed to control the background model warping. */
+    bool background_warping;
+    double offset_x;
+    double offset_y;
+  
     /** input values for this tile */
     int gridX;
@@ -81 +86 @@
 bool pswarpParseCamera (pmConfig *config);
 bool pswarpDefine (pmConfig *config);
+bool pswarpDefineBackground (pmConfig *config);
 bool pswarpLoop (pmConfig *config, psMetadata *stats);
+bool pswarpLoopBackground (pmConfig *config, psMetadata *stats);
 psExit pswarpExitCode(psExit exitValue);
 bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config);

trunk/pswarp/src/pswarpArguments.c

@@ -86 +86 @@
     pmConfigFileSetsMD (config->arguments, &argc, argv, "MASK", "-mask", "-masklist");
     pmConfigFileSetsMD (config->arguments, &argc, argv, "VARIANCE", "-variance", "-variancelist");
-
+    pmConfigFileSetsMD (config->arguments, &argc, argv, "BACKGROUND", "-background", "-bkglist");
+    
     if (argc != 3) {
         usage();
@@ -164 +165 @@
     }
 
+    bool doBKG = psMetadataLookupBool(&status,recipe, "BACKGROUND.MODEL"); ///< Generate the warped background model?
+    if (!status) {
+      doBKG = false;
+      psWarning("BACKGROUND.MODEL is not set in the %s recipe -- defaulting to FALSE.", PSWARP_RECIPE);
+    }
+    int bkgXgrid = psMetadataLookupS32(&status,recipe, "BKG.XGRID"); ///< Xsize of background model
+    if (!status) {
+      bkgXgrid = 10;
+      psWarning("BKG.XGRID is not set in the %s recipe -- defaultint to %d.",PSWARP_RECIPE,bkgXgrid);
+    }
+    int bkgYgrid = psMetadataLookupS32(&status,recipe, "BKG.YGRID"); ///< Xsize of background model
+    if (!status) {
+      bkgYgrid = 10;
+      psWarning("BKG.YGRID is not set in the %s recipe -- defaultint to %d.",PSWARP_RECIPE,bkgYgrid);
+    }
+
+    
     // Set recipe values in the recipe (since we've possibly altered some)
     psMetadataAddS32(recipe, PS_LIST_TAIL, "GRID.NX", PS_META_REPLACE,
@@ -176 +194 @@
                      "Fraction of bad flux for a pixel to be marked as poor", poorFrac);
     psMetadataAddBool(recipe, PS_LIST_TAIL, "PSF", PS_META_REPLACE, "Generate a PSF Model?", PSF);
-
+    psMetadataAddBool(recipe, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, "Generate the warped background model?", doBKG);
+    psMetadataAddS32(recipe, PS_LIST_TAIL, "BKG.XGRID", PS_META_REPLACE, "Xsize of background model", bkgXgrid);
+    psMetadataAddS32(recipe, PS_LIST_TAIL, "BKG.YGRID", PS_META_REPLACE, "Ysize of background model", bkgYgrid);
+    
+    
     // Set recipe values in the arguments
     psMetadataAddS32(config->arguments, PS_LIST_TAIL, "GRID.NX", 0,
@@ -189 +211 @@
                      "Fraction of bad flux for a pixel to be marked as poor", poorFrac);
     psMetadataAddBool(config->arguments, PS_LIST_TAIL, "PSF", PS_META_REPLACE, "Generate a PSF Model?", PSF);
+    psMetadataAddBool(config->arguments, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, "Generate the warped background model?", doBKG);
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BKG.XGRID", PS_META_REPLACE, "Xsize of background model", bkgXgrid);
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BKG.YGRID", PS_META_REPLACE, "Ysize of background model", bkgYgrid);
 
     psTrace("pswarp", 1, "Done with pswarpArguments...\n");

trunk/pswarp/src/pswarpDefine.c

@@ -41 +41 @@
         return false;
     }
-
+    
     // open the full skycell file; no need to defer different depths. only load the header data
     pmFPAview *view = pmFPAviewAlloc (0);
@@ -86 +86 @@
         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0");
         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0");
+
+        
     }
 
@@ -123 +125 @@
         }
     }
-
+    
     view->chip = view->cell = view->readout = -1;
     pmFPAAddSourceFromView(output->fpa, view, output->format);
+
 
     psFree (view);
     return true;
 }
+
+bool pswarpDefineBackground (pmConfig *config) {
+
+    // load the PSWARP recipe
+    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
+    if (!recipe) {
+        psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n");
+        return false;
+    }
+
+    // select the input data sources
+    pmFPAfile *skycell = psMetadataLookupPtr (NULL, config->files, "PSWARP.SKYCELL");
+    if (!skycell) {
+        psError(PSWARP_ERR_CONFIG, false, "Can't find skycell data!\n");
+        return false;
+    }
+    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PSWARP.OUTPUT.BKGMODEL");
+    if (!output) {
+        psError(PSWARP_ERR_CONFIG, false, "Can't find output data!\n");
+        return false;
+    }
+    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.INPUT");
+    if (!input) {
+        psError(PSWARP_ERR_CONFIG, false, "Can't find input data!\n");
+        return false;
+    }
+    
+    // open the full skycell file; no need to defer different depths. only load the header data
+    pmFPAview *view = pmFPAviewAlloc (0);
+    pmFPAfileOpen (skycell, view, config);
+
+    // Read header and create target
+    {
+        if (!pmFPAReadHeaderSet(skycell->fpa, skycell->fits, config)) {
+            psError(psErrorCodeLast(), false, "Unable to read headers for skycell.");
+            psFree(view);
+            return false;
+        }
+        view->chip = 0;
+        view->cell = 0;
+        view->readout = 0;
+        pmCell *source = pmFPAfileThisCell(config->files, view, "PSWARP.SKYCELL"); ///< Source cell
+        pmHDU *hdu = pmHDUFromCell(source); ///< HDU for source
+        if (!hdu || !hdu->header) {
+            psError(PM_ERR_PROG, false, "Unable to find header for sky cell.");
+            psFree(view);
+            return false;
+        }
+        int numCols = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); ///< Number of columns
+        int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); ///< Number of rows
+        if ((numCols == 0) || (numRows == 0)) {
+            psError(PSWARP_ERR_DATA, false, "skycell has invalid dimensions %d x %d", numCols, numRows);
+            psFree(view);
+            return false;
+        }
+
+        pmCell *target = pmFPAviewThisCell(view, output->fpa); ///< Target cell
+        pmReadout *readout = pmReadoutAlloc(target); ///< Target readout
+        readout->image = psImageAlloc(numCols / output->xBin, numRows / output->yBin, PS_TYPE_F32);
+        psImageInit(readout->image, NAN);
+        psFree(readout);                // Drop reference
+
+        bool status = false;
+        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XBIN");
+        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YBIN");
+/*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.XBIN",PS_META_REPLACE,"",output->xBin); */
+/*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.YBIN",PS_META_REPLACE,"",output->yBin); */
+        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XSIZE");
+        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YSIZE");
+/*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.XSIZE",PS_META_REPLACE,"",numCols / output->xBin); */
+/*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.YSIZE",PS_META_REPLACE,"",numRows / output->yBin); */
+        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XPARITY");
+        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YPARITY");
+        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0");
+        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0");
+
+        
+    }
+
+    // XXX this is not a sufficient test
+    view->chip = 0;
+    view->cell = 0;
+    view->readout = -1;
+    pmHDU *phu = pmFPAviewThisPHU(view, skycell->fpa); ///< Skycell PHU
+    pmHDU *hdu = pmFPAviewThisHDU(view, skycell->fpa); ///< Skycell header
+
+    pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
+
+    double cd1f = 1.0 * output->xBin;
+    double cd2f = 1.0 * output->yBin;
+    
+    WCS->crpix1 = WCS->crpix1 / cd1f;
+    WCS->crpix2 = WCS->crpix2 / cd2f;
+    
+    WCS->cdelt1 *= cd1f;
+    WCS->cdelt2 *= cd2f;
+
+    WCS->trans->x->coeff[1][0] *= cd1f;
+    WCS->trans->x->coeff[0][1] *= cd2f;
+    WCS->trans->y->coeff[1][0] *= cd1f;
+    WCS->trans->y->coeff[0][1] *= cd2f;
+    
+    
+    pmAstromWCStoHeader (hdu->header,WCS);
+
+    bool bilevelAstrometry = false;
+    if (phu) {
+        char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
+        if (ctype) {
+            bilevelAstrometry = !strcmp(&ctype[4], "-DIS");
+        }
+    }
+
+    // We read from the skycell into the output.  i.e., the output receives the desired astrometry.
+    pmChip *outputChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output
+    if (bilevelAstrometry) {
+        if (!pmAstromReadBilevelMosaic(output->fpa, phu->header)) {
+            psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for skycell.");
+            psFree(view);
+            return false;
+        }
+        if (!pmAstromReadBilevelChip(outputChip, hdu->header)) {
+            psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for skycell.");
+            psFree(view);
+            return false;
+        }
+    } else {
+        // we use a default FPA pixel scale of 1.0
+      if (!pmAstromReadWCS(output->fpa, outputChip, hdu->header, 1.0)) {
+            psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for skycell.");
+            psFree(view);
+            return false;
+        }
+    }
+
+    
+    view->chip = view->cell = view->readout = -1;
+    pmFPAAddSourceFromView(output->fpa, view, output->format);
+
+
+    psFree (view);
+    return true;
+}

trunk/pswarp/src/pswarpFileNames.h

@@ -12 +12 @@
   "PSWARP.MASK",
   "PSWARP.VARIANCE",
+  "PSWARP.BKGMODEL",
   NULL
 };
@@ -20 +21 @@
   "PSWARP.OUTPUT.MASK",
   "PSWARP.OUTPUT.VARIANCE",
+  "PSWARP.OUTPUT.BKGMODEL",
   NULL
 };

trunk/pswarp/src/pswarpLoop.c

@@ -179 +179 @@
         // read WCS data from the corresponding header
         pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
+
+        
         if (bilevelAstrometry) {
             if (!pmAstromReadBilevelChip (chip, hdu->header)) {
@@ -195 +197 @@
             }
         }
-
+        
         pmCell *cell;
         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
@@ -227 +229 @@
 
                 pswarpTransformReadout(output, readout, config);
-
+                
                 if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
                     psError(psErrorCodeLast(), false, "Unable to write files.");
@@ -255 +257 @@
         goto DONE;
     }
-
+    
     pmCell *outCell = output->parent;   ///< Output cell
     pmChip *outChip = outCell->parent;  ///< Output chip
@@ -450 +452 @@
     return true;
 }
+
+// Loop over the inputs, warp them to the output skycell and then write out the output.
+bool pswarpLoopBackground(pmConfig *config, psMetadata *stats)
+{
+    bool status;
+    bool mdok;                          // Status of MD lookup
+    const char *skyCamera = psMetadataLookupStr(NULL, config->arguments,
+                                                "SKYCELL.CAMERA");  // Name of camera for skycell
+    pmConfigCamerasCull(config, skyCamera);
+    pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
+
+    // load the recipe
+    psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE);
+    if (!recipe) {
+        psError(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE);
+        return false;
+    }
+
+    if (!pswarpSetMaskBits(config)) {
+        psError(psErrorCodeLast(), false, "failed to set mask bits");
+        return NULL;
+    }
+
+    // select the input data sources
+    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.BKGMODEL");
+    if (!input) {
+        psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n");
+        return false;
+    }
+
+    // use the external astrometry source if supplied
+    pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM");
+    if (!astrom) {
+        astrom = input;
+    }
+
+    if (astrom->camera != input->camera) {
+        psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match.");
+        return false;
+    }
+
+    // select the output readout
+    pmFPAview *view = pmFPAviewAlloc(0);
+    view->chip = 0;
+    view->cell = 0;
+    view->readout = 0;
+    pmReadout *output = pmFPAfileThisReadout(config->files, view, "PSWARP.OUTPUT.BKGMODEL");
+    if (!output) {
+        psError(PSWARP_ERR_CONFIG, true, "Can't find output background data!\n");
+        return false;
+    }
+    psFree (view);
+    // Turn all skycell files on to generate them, and then turn them off for the loop over the input images
+    // the input, which is in a different format.
+    {
+        pswarpFileActivation(config, detectorFiles, false);
+        pswarpFileActivation(config, photFiles, false);
+        pswarpFileActivation(config, independentFiles, false);
+        pswarpFileActivation(config, skycellFiles, true);
+        if (!pswarpIOChecksBefore(config)) {
+            psError(psErrorCodeLast(), false, "Unable to read files.");
+            goto DONE;
+        }
+        pswarpFileActivation(config, skycellFiles, false);
+    }
+    // Read the input astrometry
+    // XXX rather than use the activations here, this should just explicitly loop over the desired filerule
+    {
+
+      pmFPAfileActivate(config->files, true, "PSWARP.ASTROM");
+
+        pmChip *chip;
+        pmFPAview *view = pmFPAviewAlloc(0);
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+            psError(psErrorCodeLast(), false, "Unable to read files.");
+            goto DONE;
+        }
+
+        while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
+#if 0
+          // This needs to be removed because otherwise it throws an error of duplicate PSPHOT.DETECTIONS.
+            if (!chip->process || !chip->file_exists) { continue; }
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                psError(psErrorCodeLast(), false, "Unable to read files.");
+                goto DONE;
+            }
+#endif
+            pmCell *cell;
+            while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
+              psTrace ("pswarp", 4, "ACell %d: %x %x %d\n", view->cell, cell->file_exists, cell->process,psErrorCodeLast());
+                if (!cell->process || !cell->file_exists) { continue; }
+                if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) ||
+                    !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
+                    psError(psErrorCodeLast(), false, "Unable to read files.");
+                    goto DONE;
+                }
+            }
+            if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
+                psError(psErrorCodeLast(), false, "Unable to write files.");
+                goto DONE;
+            }
+        }
+        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
+            psError(psErrorCodeLast(), false, "Unable to write files.");
+            goto DONE;
+        }
+        psFree(view);
+        pswarpFileActivation(config, detectorFiles, true);
+        pmFPAfileActivate(config->files, false, "PSWARP.ASTROM");
+    }
+
+    // Don't care about the skycell anymore --- we've read it, and that's all we need to do.
+    pmFPAfileActivate(config->files, false, "PSWARP.SKYCELL");
+    view = pmFPAviewAlloc(0);
+
+    // find the FPA phu
+    bool bilevelAstrometry = false;
+    pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa);
+
+    //    pmAstromWCS *WCSF = pmAstromWCSfromHeader(phu->header);
+    
+    if (phu) {
+        char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
+        if (ctype) {
+            bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
+        }
+    }
+    if (bilevelAstrometry) {
+        if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) {
+            psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA.");
+            psFree(view);
+            psFree(stats);
+            goto DONE;
+        }
+    }
+
+    psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging
+
+    // files associated with the science image
+    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+        psError(psErrorCodeLast(), false, "Unable to read files.");
+        goto DONE;
+    }
+    pmChip *chip;
+    while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
+        psTrace ("pswarp", 4, "DChip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) { continue; }
+        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+            psError(psErrorCodeLast(), false, "Unable to read files.");
+            goto DONE;
+        }
+
+        // Pull information from the header of the background files so we can use it to set values.
+        pmHDU *hdu = pmFPAviewThisHDU(view,input->fpa);
+        psMetadata *header = hdu->header;
+        
+        int IMAXIS1 = psMetadataLookupS32(NULL,header,"IMNAXIS1");
+        int IMAXIS2 = psMetadataLookupS32(NULL,header,"IMNAXIS2");
+        int NAXIS1 = psMetadataLookupS32(NULL,header,"NAXIS1");
+        int NAXIS2 = psMetadataLookupS32(NULL,header,"NAXIS2");
+        char *CCDSUM = psMetadataLookupStr(NULL,header,"CCDSUM");
+        int CCDSUM1 = atoi(strtok(CCDSUM," "));
+        int CCDSUM2 = atoi(strtok(NULL," "));
+        
+        psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_XOFFSET", PS_META_REPLACE,
+                         "xoffset for background model data", (NAXIS1 * CCDSUM1 - IMAXIS1) / (2.0 * CCDSUM1));
+        psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_YOFFSET", PS_META_REPLACE,
+                         "yoffset for background model data", (NAXIS2 * CCDSUM2 - IMAXIS2) / (2.0 * CCDSUM2));
+        psTrace("pswarp",5,"%d %d %d %d %d %d %g %g %d %d",
+                psMetadataLookupS32(NULL,header,"IMNAXIS1"),
+                psMetadataLookupS32(NULL,header,"IMNAXIS2"),
+                psMetadataLookupS32(NULL,header,"NAXIS1"),
+                psMetadataLookupS32(NULL,header,"NAXIS2"),
+                CCDSUM1,CCDSUM2,
+                psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_XOFFSET"),
+                psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_YOFFSET"),
+                psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"),
+                psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID"));
+        
+        
+        // read WCS data from the corresponding header
+        hdu = pmFPAviewThisHDU (view, astrom->fpa);
+
+        pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
+
+        double cd1f = (1.0 * CCDSUM1);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"));
+        double cd2f = (1.0 * CCDSUM2);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID"));
+
+        WCS->cdelt1 *= cd1f;
+        WCS->cdelt2 *= cd2f;
+        WCS->crpix1 = WCS->crpix1 / cd1f;
+        WCS->crpix2 = WCS->crpix2 / cd2f;
+
+        // WCS->trans->x->nX/nY
+        for (int q = 0; q <= WCS->trans->x->nX; q++) {
+          for (int r = 0; r <= WCS->trans->x->nY; r++) {
+            WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
+          }
+        }
+        for (int q = 0; q <= WCS->trans->y->nX; q++) {
+          for (int r = 0; r <= WCS->trans->y->nY; r++) {
+            WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
+          }
+        }
+        pmAstromWCStoHeader (hdu->header,WCS);
+        
+        if (bilevelAstrometry) {
+            if (!pmAstromReadBilevelChip (chip, hdu->header)) {
+                psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
+                psFree(view);
+                psFree(stats);
+                goto DONE;
+            }
+        } else {
+            // we use a default FPA pixel scale of 1.0
+            if (!pmAstromReadWCS (astrom->fpa, chip, hdu->header, 1.0)) {
+                psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
+                psFree(view);
+                psFree(stats);
+                goto DONE;
+            }
+        }
+        // Modify structure here.
+
+        
+        pmCell *cell;
+        while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
+            psTrace ("pswarp", 4, "DCell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            if (!cell->process || !cell->file_exists) { continue; }
+            if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+                psError(psErrorCodeLast(), false, "Unable to read files.");
+                goto DONE;
+            }
+
+            psListAdd(cells, PS_LIST_TAIL, cell);
+
+            // process each of the readouts
+            pmReadout *readout;
+            while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                    psError(psErrorCodeLast(), false, "Unable to read files.");
+                    goto DONE;
+                }
+                if (!readout->data_exists) {
+                    continue;
+                }
+
+                for (int x = 0; x < readout->image->numCols; x++) {
+                  for (int y = 0; y < readout->image->numRows; y++) {
+                    readout->image->data.F32[y][x] = readout->image->data.F32[y][x] * (cd1f * cd2f) /
+                      (psMetadataLookupS32(&mdok,config->arguments,"BKG.XGRID") *
+                       psMetadataLookupS32(&mdok,config->arguments,"BKG.YGRID"));
+                  }
+                }
+                
+                psMetadataAddS32(config->arguments,PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 8);
+
+                psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", true);
+                pswarpTransformReadout(output, readout, config);
+                psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", false);
+
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                    psError(psErrorCodeLast(), false, "Unable to write files.");
+                    goto DONE;
+                }
+            }
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                psError(psErrorCodeLast(), false, "Unable to write files.");
+                goto DONE;
+            }
+        }
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+            psError(psErrorCodeLast(), false, "Unable to write files.");
+            goto DONE;
+        }
+    }
+    if (!output->data_exists) {
+        psWarning("No overlap between input and skycell.");
+        psphotFilesActivate(config, false);
+        psFree(cells);
+        psFree(view);
+        goto DONE;
+    }
+    pmCell *outCell = output->parent;   ///< Output cell
+    pmChip *outChip = outCell->parent;  ///< Output chip
+    pmFPA *outFPA = outChip->parent;    ///< Output FP
+
+/*     if (!pswarpPixelsLit(output, stats, config)) { */
+/*         psError(psErrorCodeLast(), false, "Unable to calculate pixel regions."); */
+/*         psFree(cells); */
+/*         psFree(view); */
+/*         goto DONE; */
+/*     } */
+    psRegion *trimsec = psMetadataLookupPtr(NULL, outCell->concepts, "CELL.TRIMSEC"); ///< Trim section
+    trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
+
+    if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) {
+        psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
+        psFree(stats);
+        psFree(view);
+        goto DONE;
+    }
+
+    pmHDU *hdu = outFPA->hdu;           ///< HDU for the output warped image
+
+    // Copy header from target
+    {
+        pmFPAview *skyView = pmFPAviewAlloc(0); ///< View into skycell
+        skyView->chip = skyView->cell = 0;
+        pmCell *cell = pmFPAfileThisCell(config->files, skyView, "PSWARP.SKYCELL"); // Skycell cell
+        psFree(skyView);
+        pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU
+        if (!skyHDU) {
+            psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU.");
+            psFree(view);
+            goto DONE;
+        }
+        hdu->header = psMetadataCopy(hdu->header, skyHDU->header);
+    }
+    pswarpVersionHeader(hdu->header);
+    
+    if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
+        psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
+        psFree(stats);
+        goto DONE;
+    }
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        psError(psErrorCodeLast(), false, "Unable to write files.");
+        goto DONE;
+    }
+    // Done with the detector side of things
+    pswarpFileActivation(config, detectorFiles, false);
+    pswarpFileActivation(config, independentFiles, false);
+
+
+    // Add MD5 information for readout
+    const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME");
+    const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME");
+    psString headerName = NULL; ///< Header name for MD5
+    psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
+    psVector *md5 = psImageMD5(output->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);
+    psFree(view);
+ DONE:
+
+    return true;
+}

trunk/pswarp/src/pswarpParseCamera.c

@@ -84 +84 @@
     }
 
+    pmFPAfile *inBackground = defineInputFile(config, NULL, "PSWARP.BKGMODEL", "BACKGROUND",
+                                              PM_FPA_FILE_IMAGE);
+    if (!inBackground) {
+      psLogMsg("pswarp", PS_LOG_INFO, "No background models supplied");
+    }
+    
     // The input skycell is a required argument: it defines the output image
     // XXX we may need a different skycell structure here
@@ -134 +140 @@
     }
 
+    if (inBackground && psMetadataLookupBool(&mdok, recipe, "BACKGROUND.MODEL")) {
+      pmFPAfile *outBackground = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT.BKGMODEL");
+/*       pmFPAfile *outBackground = pmFPAfileDefineFromFPA(config,output->fpa, */
+/*                                                      psMetadataLookupS32(&mdok,recipe,"BKG.XGRID"), */
+/*                                                      psMetadataLookupS32(&mdok,recipe,"BKG.YGRID"), */
+/*                                                      "PSWARP.OUTPUT.BKGMODEL"); */
+      outBackground->xBin = psMetadataLookupS32(&mdok, recipe, "BKG.XGRID");
+      outBackground->yBin = psMetadataLookupS32(&mdok, recipe, "BKG.YGRID");
+      
+      if (!outBackground) {
+        psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.BKGMODEL");
+        return false;
+      }
+      outBackground->save = true;
+    }
+    
     if (psMetadataLookupBool(&mdok, recipe, "PSF")) {
         // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by

trunk/pswarp/src/pswarpTransformReadout.c

@@ -29 +29 @@
     psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments,
                                                                    "INTERPOLATION.MODE"); ///< Mode for interp
+
     int numKernels = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.NUM"); ///< Number of kernels
 
@@ -60 +61 @@
     // output coordinates to input coordinates
     pswarpMapGrid *grid = pswarpMapGridFromImage(input, output, nGridX, nGridY);
-
+    //    if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) {
+    
     // XXX optionally modify the grid based on this result and force the maxError < XXX
     double maxError = pswarpMapGridMaxError(grid); // Maximum (positional) error from using grid
@@ -130 +132 @@
             args->goodPixels = 0;
 
+            if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) {
+              args->background_warping = true;
+              args->offset_x = psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_XOFFSET");
+              args->offset_y = psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_YOFFSET");
+            }
+
+            
             // allocate a job
             psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE");
@@ -171 +180 @@
         } else {
             pswarpTransformTileArgs *args = job->args->data[0];
-            // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n);
             goodPixels += args->goodPixels;
             xMin = PS_MIN(args->xMin, xMin);
@@ -183 +191 @@
                 jacobian += args->jacobian * args->goodPixels;
             }
+
+            
         }
         psFree(job);
@@ -201 +211 @@
 
     if (goodPixels > 0 && psMetadataLookupBool(&mdok, recipe, "SOURCES")) {
+      if (!psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) {
         if (!pswarpTransformSources(output, input, config)) {
-            psError(psErrorCodeLast(), false, "Unable to interpolate image.");
-            return false;
+          psError(psErrorCodeLast(), false, "Unable to interpolate image.");
+          return false;
         }
+      }
     }
 

trunk/pswarp/src/pswarpTransformSources.c

@@ -41 +41 @@
     if (!outDetections) {
         outDetections = pmDetectionsAlloc();
-        psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Warped sources", outDetections);
+        psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY , "Warped sources", outDetections);
     }
     psArray *outSources = outDetections->allSources;

trunk/pswarp/src/pswarpTransformTile.c

@@ -46 +46 @@
     args->jacobian = NAN;
 
+    args->background_warping = false;
+    args->offset_x = 0.0;
+    args->offset_y = 0.0;
+    
     return args;
 }
@@ -83 +87 @@
     // Iterate over the output image pixels (parent frame)
     long goodPixels = 0;                ///< Number of input pixels landing on the output image
+
     for (int y = yMin; y < yMax; y++) {
         for (int x = xMin; x < xMax; x++) {
-
             // Only transform those pixels requested
             if (region && region->data.PS_TYPE_IMAGE_MASK_DATA[y][x]) {
                 continue;
             }
-
             // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
             // both are in the parent frames of the input and output images.
             double xIn, yIn;            // Input pixel coordinates
             pswarpMapApply(&xIn, &yIn, map, x + 0.5, y + 0.5);
+
+            // This needs to use a more reliable method to do this offset and limiting
+            //      if (args->interp->mode == 8) {
+            if (args->background_warping) {
+              //              double xOffset = 177.0 / 400.0; // (modelsize * modelbinning - xsize) / 2.0
+              //              double yOffset = 166.0 / 400.0; // (modelsize * modelbinning - ysize) / 2.0
+              xIn += args->offset_x;
+              yIn += args->offset_y;
+
+              if ((xIn > inNumCols - args->offset_x)||
+                  (yIn > inNumRows - args->offset_y)||
+                  (xIn < args->offset_x)||
+                  (yIn < args->offset_y)) {
+                continue;
+              }
+            }
+            
             if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) {
                 continue;