Index: /trunk/Nebulous-Server/bin/neb-voladm =================================================================== --- /trunk/Nebulous-Server/bin/neb-voladm (revision 34799) +++ /trunk/Nebulous-Server/bin/neb-voladm (revision 34800) @@ -68,5 +68,5 @@ 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]$/); } Index: /trunk/extsrc/gpcsw/gpcsrc/Make.Common =================================================================== --- /trunk/extsrc/gpcsw/gpcsrc/Make.Common (revision 34799) +++ /trunk/extsrc/gpcsw/gpcsrc/Make.Common (revision 34800) @@ -156,8 +156,8 @@ # 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 Index: /trunk/ippScripts/Build.PL =================================================================== --- /trunk/ippScripts/Build.PL (revision 34799) +++ /trunk/ippScripts/Build.PL (revision 34800) @@ -128,4 +128,5 @@ scripts/listvideocells.pl scripts/skycalibration.pl + scripts/regenerate_background.pl )], dist_abstract => 'Scripts for running the Pan-STARRS IPP', Index: /trunk/ippScripts/MANIFEST =================================================================== --- /trunk/ippScripts/MANIFEST (revision 34799) +++ /trunk/ippScripts/MANIFEST (revision 34800) @@ -46,3 +46,4 @@ scripts/skycell_jpeg.pl scripts/lap_science.pl +scripts/regenerate_background.pl t/00_distribution.t Index: /trunk/ippScripts/scripts/detrend_process_imfile.pl =================================================================== --- /trunk/ippScripts/scripts/detrend_process_imfile.pl (revision 34799) +++ /trunk/ippScripts/scripts/detrend_process_imfile.pl (revision 34800) @@ -27,4 +27,6 @@ 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,4 +126,11 @@ # 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"; Index: /trunk/ippScripts/scripts/regenerate_background.pl =================================================================== --- /trunk/ippScripts/scripts/regenerate_background.pl (revision 34800) +++ /trunk/ippScripts/scripts/regenerate_background.pl (revision 34800) @@ -0,0 +1,358 @@ +#! /usr/bin/env perl + +# Script to call the appropriate commands to regenerate the background models +# for stages subsequent to the chip stage. + +use warnings; +use strict; +use Carp; + +use IPC::Cmd 0.36 qw( can_run run ); +use File::Spec; +use File::Temp qw( tempfile ); +use PS::IPP::Config 1.01 qw( :standard ); +use Getopt::Long qw (GetOptions :config auto_help auto_version gnu_getopt ); +use Pod::Usage qw( pod2usage ); +use PS::IPP::Metadata::Config; +use PS::IPP::Metadata::List qw( parse_md_list ); + + + +my $ppConfigDump = can_run('ppConfigDump'); +my $pswarp = can_run('pswarp'); + +# Parse command line options. +my ($stage, $camera, $stage_id, $skycell_id, $dbname, $verbose, $logfile, $save_temps); +GetOptions('stage=s' => \$stage, # which analysis stage to generate bkg model + 'camera|i=s' => \$camera, # user-supplied camera name + 'stage_id=s' => \$stage_id, # id for this stage + 'skycell_id=s' => \$skycell_id, # skycell_id + 'verbose' => \$verbose, # verbose commands + 'logfile=s' => \$logfile, # destination for stdout and stderr + 'save_temps' => \$save_temps, # save temporary files. + 'dbname=s' => \$dbname # database name + ) or pod2usage ( 2 ); + +pod2usage( -msg => "Usage: $0 --camera (name) --stage (stage) --stage_id (stage_id) [--dbname (dbname)]", + -exitval => 2 ) if scalar @ARGV; + +pod2usage( -msg => "Required options: --camera (name) --stage (stage) --stage_id (stage_id)", + -exitval => 3 ) unless + defined $camera and + defined $stage and + defined $stage_id; + +pod2usage( -msg => "Required options: --camera (name) --stage (stage) --stage_id (stage_id) --skycell_id (skycell_id)", + -exitval => 3) if (($stage eq 'warp')&&(!defined($skycell_id))); + + +my $ipprc = PS::IPP::Config->new ( $camera ) or + &my_die("Unable to set up", $stage, $stage_id, $PS_EXIT_CONFIG_ERROR); # used for path/neb info + +# We can only regenerate a background for certain stages. +unless (($stage eq "warp") || ($stage eq "stack") || ($stage eq "warpstack")) { + &my_die("Invalid stage to regenerate", $stage, $stage_id, $PS_EXIT_CONFIG_ERROR); +} + +$ipprc->redirect_output($logfile) or + &my_die("Unable to redirect output to logfile", $stage, $stage_id, $PS_EXIT_UNKNOWN_ERROR) if $logfile; + +my $mdcParser = PS::IPP::Metadata::Config->new; # parser for metadata config files + + +## warp stage +if ($stage eq 'warp') { + my $warp_id = $stage_id; # Because I think I used the wrong id in places. + &my_die("--stage_id required for stage warp", $stage, $stage_id, $PS_EXIT_CONFIG_ERROR) if !$stage_id; + &my_die("--skycell_id required for stage warp", $stage, $stage_id, $PS_EXIT_CONFIG_ERROR) if !$skycell_id; + + # Configuration stuff + my $warptool = can_run('warptool') or + &my_die("Can't find warptool",$stage,$stage_id, $PS_EXIT_UNKNOWN_ERROR); + + my $astromSource; # The astrometry source + { + my $command = "$ppConfigDump -camera $camera -recipe PSWARP BACKGROUND -dump-recipe PSWARP -"; + my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = + run(command => $command, verbose => $verbose); + unless ($success) { + $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); + &my_die("Unable to perform ppConfigDump: $error_code", + $stage, $stage_id, $error_code); + } + my $metadata = $mdcParser->parse(join "", @$stdout_buf) or + &my_die("Unable to parse metadata config doc", + $stage, $stage_id, $PS_EXIT_PROG_ERROR); + $astromSource = metadataLookupStr($metadata, 'ASTROM.SOURCE'); + } + + + my $imfiles; # Array of component files + my $command = "$warptool -warped -warp_id $stage_id -skycell_id $skycell_id"; + $command .= " -dbname $dbname " if defined $dbname; + + my ( $success, $error_code, $full_buf, $stdout_buf, $stderrr_buf ) = + run(command => $command, verbose => $verbose); + unless ($success) { + $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); + &my_die("Unable to perform warptool: $error_code", $stage, $stage_id, $error_code); + } + + if (@$stdout_buf == 0) { + # No results + } + + my $in_md = $mdcParser->parse(join "", @$stdout_buf) or + &my_die("Unable to parse metadata config doc", $stage, $stage_id, $PS_EXIT_PROG_ERROR); + $imfiles = parse_md_list($in_md) or + &my_die("Unable to parse metadata", $stage, $stage_id, $PS_EXIT_PROG_ERROR); + + foreach my $imfile (@$imfiles) { + my $skycell_id = $imfile->{skycell_id}; + my $path_base = $imfile->{path_base}; + my $tess_id = $imfile->{tess_id}; + + # Parse the results to determine if there are any problems that prevent this from + # being processed. I think I might need to think about this a bit more. + my $status = 1; + $status = 0 unless defined $path_base and $path_base ne "NULL"; + my $poor_quality = $imfile->{quality} > 0; + +# my $bkg_file = $ipprc->filename("PSWARP.OUTPUT.BKGMODEL",$path_base,$skycell_id); +# my $bkg_exists = $ipprc->file_exists($bkg_file); + my $bkg_exists = $imfile->{background_model}; + if ($bkg_exists) { $status = 0; } # We do not need to remake this. + + # Construct temp files + my ($bkgList_file,$bkgList_name) = tempfile("/tmp/bkgreg.warp.bkg.${warp_id}.${skycell_id}.XXXX", + UNLINK => !$save_temps); + my ($astromList_file,$astromList_name) = tempfile("/tmp/bkgreg.warp.astrom.${warp_id}.${skycell_id}.XXXX", + UNLINK => !$save_temps); + my $chipFiles; + { + my $inputs_command = "$warptool -scmap -warp_id $warp_id "; + $inputs_command .= " -skycell_id $skycell_id "; + $inputs_command .= " -dbname $dbname " if defined $dbname; + my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = + run (command => $inputs_command, verbose => $verbose); + unless ($success) { + $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); + &my_die("Unable to perform warptool -scmap: $error_code", + $skycell_id,$warp_id,$PS_EXIT_SYS_ERROR); + } + my $md = $mdcParser->parse(join "", @$stdout_buf) or + &my_die("Unable to parse warptool -scmap: $error_code", + $skycell_id,$warp_id,$PS_EXIT_SYS_ERROR); + $chipFiles = parse_md_list($md) or + &my_die("Unable to parse warptool -scmap: $error_code", + $skycell_id,$warp_id,$PS_EXIT_SYS_ERROR); + } + my $wrote_astrom = 0; + foreach my $chipFile (@$chipFiles) { + my $bkgModel = $ipprc->filename("PSPHOT.BACKMDL", + $chipFile->{chip_path_base}, + $chipFile->{class_id}); + my $astrom = $ipprc->filename($astromSource,$chipFile->{cam_path_base}); + print $bkgList_file "$bkgModel\n"; + if (!$wrote_astrom) { + print $astromList_file "$astrom\n"; + $wrote_astrom = 1; + } + } + close ($bkgList_file); + close ($astromList_file); + + # Construct skyfile if needed + my $skyCell_name = $ipprc->filename("SKYCELL.TEMPLATE",$path_base,$skycell_id); + my $skyCell_exists = $ipprc->file_exists($skyCell_name); + if (!$skyCell_exists) { # Generate it + my $skyFile = prepare_output("SKYCELL.TEMPLATE", $path_base, $skycell_id, 1); + $ipprc->skycell_file( $tess_id, $skycell_id, $skyFile, $verbose ) or + &my_die("Unable to generate template skycell",$skycell_id,$warp_id,$PS_EXIT_SYS_ERROR); + } + my $bkgOut_log = prepare_output("LOG.EXP",$path_base . ".BKG_REG",$skycell_id,1); + + # Construct pswarp command + my $pswarp_command = "$pswarp "; + $pswarp_command .= " -list $bkgList_name -astromlist $astromList_name -bkglist $bkgList_name "; + $pswarp_command .= " ${path_base} "; + $pswarp_command .= " $skyCell_name "; + $pswarp_command .= " -F PSPHOT.PSF.SAVE PSPHOT.PSF.SKY.SAVE -F PSPHOT.OUTPUT PSPHOT.OUT.CMF.MEF "; + $pswarp_command .= " -F PSPHOT.BACKMDL PSPHOT.BACKMDL.MEF "; + $pswarp_command .= " -F SOURCE.PLOT.MOMENTS SOURCE.PLOT.SKY.MOMENTS "; + $pswarp_command .= " -F SOURCE.PLOT.PSFMODEL SOURCE.PLOT.SKY.PSFMODEL "; + $pswarp_command .= " -F SOURCE.PLOT.APRESID SOURCE.PLOT.SKY.APRESID "; + $pswarp_command .= " -recipe PSWARP WARP -Db PSF FALSE -Db BACKGROUND.MODEL TRUE "; + $pswarp_command .= " -Db SOURCES FALSE "; + $pswarp_command .= " -log $bkgOut_log -threads 1 "; + $pswarp_command .= " -dbname $dbname " if defined $dbname; + + # Execute + my ( $success, $error_code, $full_buf, $stdout_buf, $stderrr_buf ) = + run(command => $pswarp_command, verbose => $verbose); + unless ($success) { + $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); + &my_die("Unable to perform pswarp: $error_code", $stage, $stage_id, $error_code); + } + + my $update_command = "$warptool "; + $update_command .= " -updateskyfile -warp_id $warp_id -skycell_id $skycell_id "; + $update_command .= " -set_background_model 1 "; + $update_command .= " -dbname $dbname " if defined $dbname; + + ( $success, $error_code, $full_buf, $stdout_buf, $stderrr_buf ) = + run(command => $update_command, verbose => $verbose); + unless ($success) { + $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); + &my_die("Unable to perform warptool: $error_code", $stage, $stage_id, $error_code); + } + + + } +} +## stack stage +elsif ($stage eq 'stack') { + &my_die("--stage_id required for stage stack", $stage, $stage_id, $PS_EXIT_CONFIG_ERROR) if !$stage_id; + my $stack_id = $stage_id; # Same as above. Alias this so I don't make mistakes. + + # Configuration stuff + my $stacktool = can_run('stacktool') or + &my_die("Can't find stacktool",$stage,$stage_id,$PS_EXIT_UNKNOWN_ERROR); + my $ppStackMedian = can_run('ppStackMedian') or + &my_die("Can't find ppStackMedian",$stage,$stage_id,$PS_EXIT_UNKNOWN_ERROR); + + # Get the information about this run + my $st_command = "$stacktool -sumskyfile -stack_id $stack_id"; + $st_command .= " -dbname $dbname " if defined $dbname; + my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = + run(command => $st_command, verbose => $verbose); + unless($success) { + $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); + &my_die("Unable to perform stacktool -sumskyfile",$stage,$stage_id,$error_code); + } + if (@$stdout_buf == 0) { + # Nothing to do; + } + my $in_md = $mdcParser->parse(join "", @$stdout_buf) or + &my_die("Unable to parse metadata config doc", $stage, $stage_id, $PS_EXIT_PROG_ERROR); + my $stacks = parse_md_list($in_md) or + &my_die("Unable to parse metadata", $stage, $stage_id, $PS_EXIT_PROG_ERROR); + my $stack = shift(@$stacks); + my $path_base = $stack->{path_base}; + + # Get inputs + my $imfiles; # Array of component files + my $command = "$stacktool -inputskyfile -stack_id $stack_id"; + $command .= " -dbname $dbname " if defined $dbname; + ( $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 stacktool -inputskyfile",$stage,$stage_id,$error_code); + } + if (@$stdout_buf == 0) { + # Nothing to do; + } + $in_md = $mdcParser->parse(join "", @$stdout_buf) or + &my_die("Unable to parse metadata config doc", $stage, $stage_id, $PS_EXIT_PROG_ERROR); + $imfiles = parse_md_list($in_md) or + &my_die("Unable to parse metadata", $stage, $stage_id, $PS_EXIT_PROG_ERROR); + + my $num = 0; + my $expTime; + my ($inputMDC_file, $inputMDC_name) = tempfile("/tmp/bkgreg.stack.mdc.${stage_id}.XXXX", + UNLINK => !$save_temps); + foreach my $imfile (@$imfiles) { + if ($imfile->{ignored}) {next; } + unless ($imfile->{background_model}) { + &my_die("Not all inputs have valid background models", + $stage,$stage_id,$PS_EXIT_PROG_ERROR); + } + print $inputMDC_file "INPUT${num}\tMETADATA\n"; + + print $inputMDC_file "\tIMAGE\tSTR\t" . $ipprc->filename( "PSWARP.OUTPUT.BKGMODEL", + $imfile->{path_base} ) . "\n"; + print $inputMDC_file "\tSOURCES\tSTR\t" . $ipprc->filename( "PSWARP.OUTPUT.SOURCES", + $imfile->{path_base} ) . "\n"; + print $inputMDC_file "END\n\n"; + $num++; + + $expTime = $imfile->{exp_time}; # I need something here. + } + close($inputMDC_file); + + my $bkgOut_log = prepare_output("LOG.EXP",$path_base . ".BKG_REG",1); + + my $ppStack_command = "$ppStackMedian -input $inputMDC_name ${path_base}.mdl "; + $ppStack_command .= " -recipe PPSTACK STACK -recipe PPSUB STACK -recipe PSPHOT STACK "; + $ppStack_command .= " -recipe PPSTATS STACKSTATS -stack-type DEEP_STACK "; +# $ppStack_command .= " -F PSPHOT.PSF.SAVE PSPHOT.PSF.SKY.SAVE "; +# $ppStack_command .= " -F PSPHOT.OUTPUT PSPHOT.OUT.CMF.MEF -F PSPHOT.BACKMDL PSPHOT.BACKMDL.MEF "; +# $ppStack_command .= " -F SOURCE.PLOT.MOMENTS SOURCE.PLOT.SKY.MOMENTS "; +# $ppStack_command .= " -F SOURCE.PLOT.PSFMODEL SOURCE.PLOT.SKY.PSFMODEL "; +# $ppStack_command .= " -F SOURCE.PLOT.APRESID SOURCE.PLOT.SKY.APRESID "; + $ppStack_command .= " -Db PHOTOMETRY F -Db VARIANCE F -Db CONVOLVE F "; + $ppStack_command .= " -Df DEFAULT.EXPTIME $expTime "; + $ppStack_command .= " -threads 1 -log $bkgOut_log "; + $ppStack_command .= " -dbname $dbname " if defined $dbname; + + ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = + run(command => $ppStack_command, verbose => $verbose); + unless($success) { + $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); + &my_die("Unable to perform ppStackMedian",$stage,$stage_id,$error_code); + } + + my $update_command = "$stacktool "; + $update_command .= " -updatesumskyfile -stack_id $stack_id "; + $update_command .= " -set_background_model 1 -fault 0 "; + $update_command .= " -dbname $dbname " if defined $dbname; + + ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = + run(command => $update_command, verbose => $verbose); + unless ($success) { + $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); + &my_die("Unable to perform stacktool: $error_code", $stage, $stage_id, $error_code); + } + + + +} +## do both warp and stack +elsif ($stage eq 'warpstack') { + +} + +sub prepare_output +{ + my $filerule = shift; + my $outroot = shift; + my $skycell_id = shift; + my $delete = shift; + $delete = 0 if !defined $delete; + + my $error; + my $output = $ipprc->prepare_output($filerule, $outroot, $skycell_id, $delete, \$error) + or &my_die("failed to prepare output file for: $filerule", $skycell_id, $stage_id, $error); + return $output; +} + + +sub my_die { + my $msg = shift; # Warning message on die + my $stage = shift; # stage name + my $stage_id = shift; # identifier + my $exit_code = shift; # exit code + + carp($msg); + exit $exit_code; +} + + + + + + + + + Index: /trunk/ippScripts/scripts/skycell_jpeg.pl =================================================================== --- /trunk/ippScripts/scripts/skycell_jpeg.pl (revision 34799) +++ /trunk/ippScripts/scripts/skycell_jpeg.pl (revision 34800) @@ -268,4 +268,12 @@ 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,42 +289,51 @@ $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"; Index: /trunk/ippScripts/scripts/stack_skycell.pl =================================================================== --- /trunk/ippScripts/scripts/stack_skycell.pl (revision 34799) +++ /trunk/ippScripts/scripts/stack_skycell.pl (revision 34800) @@ -229,4 +229,6 @@ 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,4 +242,5 @@ 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"; Index: /trunk/ippScripts/scripts/warp_skycell.pl =================================================================== --- /trunk/ippScripts/scripts/warp_skycell.pl (revision 34799) +++ /trunk/ippScripts/scripts/warp_skycell.pl (revision 34800) @@ -128,4 +128,5 @@ # 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,4 +140,5 @@ &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,4 +170,8 @@ 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,5 +215,8 @@ 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,4 +244,9 @@ 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,5 +259,7 @@ 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,4 +285,5 @@ $command .= " -variancelist $weightName"; $command .= " -astromlist $astromName"; + $command .= " -bkglist $bkgName" if ($doBackground); $command .= " $outroot $skyFile"; $command .= " -F PSPHOT.PSF.SAVE PSPHOT.PSF.SKY.SAVE"; Index: /trunk/ippTasks/Makefile.am =================================================================== --- /trunk/ippTasks/Makefile.am (revision 34799) +++ /trunk/ippTasks/Makefile.am (revision 34800) @@ -45,5 +45,6 @@ diffphot.pro \ lap.pro \ - vp.pro + vp.pro \ + bg.regeneration.pro other_files = \ Index: /trunk/ippTasks/bg.regeneration.pro =================================================================== --- /trunk/ippTasks/bg.regeneration.pro (revision 34800) +++ /trunk/ippTasks/bg.regeneration.pro (revision 34800) @@ -0,0 +1,306 @@ +## bg.regeneration.pro : tasks for regenerating the background model for warps and stacks : -*- sh -*- + +check.globals + +if ($?POLL_LIMIT_BGREG == 0) set POLL_LIMIT_BGREG = 40 + +macro set.bgreg.poll + if ($0 != 2) + echo "USAGE:set.bgreg.poll (value)" + break; + end + + $POLL_LIMIT_BGREG = $1 +end +macro get.bgreg.poll + echo $POLL_LIMIT_BGREG +end + +## Initialize the books for the tasks to do +book init bgRegWarp +book init bgRegStack + +## Database lists +$bgRegWarp_DB = 0 +$bgRegStack_DB = 0 + +macro bgreg.warp.status + book listbook bgRegWarp +end +macro bgreg.stack.status + book listbook bgRegStack +end + +## Reset tasks +macro bgreg.reset + book init bgRegWarp + book init bgRegStack +end + +## Turn tasks on +macro bgreg.on + task bgreg.warp.load + active true + end + task bgreg.warp.run + active true + end + task bgreg.stack.load + active true + end + task bgreg.stack.run + active true + end +end +## or off +macro bgreg.off + task bgreg.warp.load + active false + end + task bgreg.warp.run + active false + end + task bgreg.stack.load + active false + end + task bgreg.stack.run + active false + end +end + +# Load task for warp +task bgreg.warp.load + host local + + periods -poll $LOADPOLL + periods -poll $LOADEXEC + periods -timeout 30 + npending 1 + + stdout NULL + stderr $LOGDIR/bgreg.warp.load + + task.exec + if ($LABEL:n == 0) break + $run = warptool -warped -background_model 0 + if ($DB:n == 0) + option DEFAULT + else + #save the DB naem for the exit tasks + option $DB:$bgRegWarp_DB + $run = $run -dbname $DB:$bgRegWarp_DB + $bgRegWarp_DB ++ + if ($bgRegWarp_DB >= $DB:n) set bgRegWarp_DB = 0 + end + add_poll_args run + add_poll_labels run + # change the limit? + $run = $run -limit $POLL_LIMIT_BGREG + command $run + end + + # success + task.exit 0 + # convert 'stdout' to book format + ipptool2book stdout bgRegWarp -key warp_id:skycell_id -uniq -setword dbname $options:0 -setword pantaskState INIT + if ($VERBOSE > 2) + book listbook bgRegWarp + end + + # delete existing entries in the appropriate pantaskStates + process_cleanup bgRegWarp + end + + # locked list + task.exit default + showcommand failure + end + + task.exit crash + showcommand crash + end + + # operation times out? + task.exit timeout + showcommand timeout + end +end + +task bgreg.warp.run + periods -poll $RUNPOLL + periods -exec $RUNEXEC + periods -timeout 60 + + task.exec + # if we can't exec, set the retry time long + periods -exec $RUNEXEC + + book npages bgRegWarp -var N + if ($N == 0) break + if ($NETWORK == 0) break + + # look for new pages + book getpage bgRegWarp 0 -var pageName -key pantaskState INIT + if ("$pageName" == "NULL") break + + book setword bgRegWarp $pageName pantaskState RUN + book getword bgRegWarp $pageName warp_id -var STAGE_ID + book getword bgRegWarp $pageName skycell_id -var SKYCELL_ID + book getword bgRegWarp $pageName camera -var CAMERA + book getword bgRegWarp $pageName dbname -var DBNAME + + set.host.for.skycell $SKYCELL_ID + + stdout $LOGDIR/bgreg.warp.log + stderr $LOGDIR/bgreg.warp.log + + $run = regenerate_background.pl --stage warp --stage_id $STAGE_ID --skycell_id $SKYCELL_ID --camera $CAMERA + + add_standard_args run + options $pageName + + if ($VERBOSE > 1) + echo command $run + end + periods -exec 0.05 + command $run + end + + # default exit status + task.exit default + process_exit bgRegWarp $options:0 $JOB_STATUS + end + + # locked list + task.exit crash + showcommand crash + echo "hostname: $JOB_HOSTNAME" + book setword bgRegWarp $options:0 pantasksState CRASH + end + + # operation timed out + task.exit timeout + showcommand timeout + book setword bgRegWarp $options:0 pantaskState TIMEOUT + end +end + +# Load task for stack +task bgreg.stack.load + host local + + periods -poll $LOADPOLL + periods -poll $LOADEXEC + periods -timeout 30 + npending 1 + + stdout NULL + stderr $LOGDIR/bgreg.stack.load + + task.exec + if ($LABEL:n == 0) break + $run = stacktool -tobkg + if ($DB:n == 0) + option DEFAULT + else + #save the DB naem for the exit tasks + option $DB:$bgRegStack_DB + $run = $run -dbname $DB:$bgRegStack_DB + $bgRegStack_DB ++ + if ($bgRegStack_DB >= $DB:n) set bgRegStack_DB = 0 + end + add_poll_args run + add_poll_labels run + # change the limit? + $run = $run -limit $POLL_LIMIT_BGREG + command $run + end + + # success + task.exit 0 + # convert 'stdout' to book format + ipptool2book stdout bgRegStack -key stack_id:skycell_id -uniq -setword dbname $options:0 -setword pantaskState INIT + if ($VERBOSE > 2) + book listbook bgRegStack + end + + # delete existing entries in the appropriate pantaskStates + process_cleanup bgRegStack + end + + # locked list + task.exit default + showcommand failure + end + + task.exit crash + showcommand crash + end + + # operation times out? + task.exit timeout + showcommand timeout + end +end + +task bgreg.stack.run + periods -poll $RUNPOLL + periods -exec $RUNEXEC + periods -timeout 60 + + task.exec + # if we can't exec, set the retry time long + periods -exec $RUNEXEC + + book npages bgRegStack -var N + if ($N == 0) break + if ($NETWORK == 0) break + + # look for new pages + book getpage bgRegStack 0 -var pageName -key pantaskState INIT + if ("$pageName" == "NULL") break + + book setword bgRegStack $pageName pantaskState RUN + book getword bgRegStack $pageName stack_id -var STAGE_ID + book getword bgRegStack $pageName skycell_id -var SKYCELL_ID + book getword bgRegStack $pageName camera -var CAMERA + book getword bgRegStack $pageName dbname -var DBNAME + set.host.for.skycell $SKYCELL_ID + + stdout $LOGDIR/bgreg.stack.log + stderr $LOGDIR/bgreg.stack.log + + $run = regenerate_background.pl --stage stack --stage_id $STAGE_ID --skycell_id $SKYCELL_ID --camera $CAMERA + + add_standard_args run + options $pageName + + if ($VERBOSE > 1) + echo command $run + end + periods -exec 0.05 + command $run + end + + # default exit status + task.exit default + process_exit bgRegStack $options:0 $JOB_STATUS + end + + # locked list + task.exit crash + showcommand crash + echo "hostname: $JOB_HOSTNAME" + book setword bgRegStack $options:0 pantasksState CRASH + end + + # operation timed out + task.exit timeout + showcommand timeout + book setword bgRegStack $options:0 pantaskState TIMEOUT + end +end + + + + Index: /trunk/ippTools/share/Makefile.am =================================================================== --- /trunk/ippTools/share/Makefile.am (revision 34799) +++ /trunk/ippTools/share/Makefile.am (revision 34800) @@ -393,4 +393,5 @@ stacktool_sassskyfile.sql \ stacktool_tosum.sql \ + stacktool_tobkg.sql \ stacktool_tosummary.sql \ stacktool_addsummary.sql \ Index: /trunk/ippTools/share/stacktool_inputskyfile.sql =================================================================== --- /trunk/ippTools/share/stacktool_inputskyfile.sql (revision 34799) +++ /trunk/ippTools/share/stacktool_inputskyfile.sql (revision 34800) @@ -4,4 +4,5 @@ rawExp.exp_id, rawExp.exp_name, + rawExp.exp_time, rawExp.object, rawExp.dateobs, Index: /trunk/ippTools/share/stacktool_tobkg.sql =================================================================== --- /trunk/ippTools/share/stacktool_tobkg.sql (revision 34800) +++ /trunk/ippTools/share/stacktool_tobkg.sql (revision 34800) @@ -0,0 +1,29 @@ +SELECT + stackRun.stack_id, + stackRun.tess_id, + stackRun.skycell_id, + stackRun.workdir, + stackRun.reduction, + stackRun.label, + stackRun.state, + stackSumSkyfile.path_base, + rawExp.camera, + IFNULL(Label.priority, 10000) AS priority +FROM stackRun +JOIN stackInputSkyfile USING(stack_id) +JOIN warpRun USING(warp_id) +JOIN warpSkyfile USING(warp_id,skycell_id) +JOIN fakeRun USING(fake_id) +JOIN camRun USING(cam_id) +JOIN chipRun USING(chip_id) +JOIN rawExp USING(exp_id) +LEFT JOIN stackSumSkyfile USING(stack_id) +LEFT JOIN Label ON Label.label = stackRun.label +WHERE + ((stackRun.state = 'full' AND stackSumSkyfile.fault = 0 AND stackSumSkyfile.quality = 0 AND stackSumSkyfile.background_model != 1)) + AND (Label.active OR Label.active IS NULL) + -- WHERE hook %s +GROUP BY stack_id +HAVING (SUM(IF(warpRun.state = 'full', 1, 0)) = COUNT(stackInputSkyfile.warp_id) AND + SUM(IF(warpSkyfile.background_model = 1, 1, 0)) = COUNT(stackInputSkyfile.warp_id)) + Index: /trunk/ippTools/src/stacktool.c =================================================================== --- /trunk/ippTools/src/stacktool.c (revision 34799) +++ /trunk/ippTools/src/stacktool.c (revision 34800) @@ -37,4 +37,5 @@ 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,4 +77,5 @@ 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,4 +954,79 @@ // 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,12 +1758,26 @@ 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); Index: /trunk/ippTools/src/stacktool.h =================================================================== --- /trunk/ippTools/src/stacktool.h (revision 34799) +++ /trunk/ippTools/src/stacktool.h (revision 34800) @@ -31,4 +31,5 @@ STACKTOOL_MODE_INPUTSKYFILE, STACKTOOL_MODE_TOSUM, + STACKTOOL_MODE_TOBKG, STACKTOOL_MODE_ADDSUMSKYFILE, STACKTOOL_MODE_SUMSKYFILE, Index: /trunk/ippTools/src/stacktoolConfig.c =================================================================== --- /trunk/ippTools/src/stacktoolConfig.c (revision 34799) +++ /trunk/ippTools/src/stacktoolConfig.c (revision 34800) @@ -150,4 +150,11 @@ 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,5 +287,6 @@ 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,4 +312,5 @@ 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); Index: /trunk/ippTools/src/warptool.c =================================================================== --- /trunk/ippTools/src/warptool.c (revision 34799) +++ /trunk/ippTools/src/warptool.c (revision 34800) @@ -2059,4 +2059,5 @@ 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,4 +2066,8 @@ 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,6 +2130,22 @@ // 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", "=="); Index: /trunk/ippTools/src/warptoolConfig.c =================================================================== --- /trunk/ippTools/src/warptoolConfig.c (revision 34799) +++ /trunk/ippTools/src/warptoolConfig.c (revision 34800) @@ -410,5 +410,5 @@ 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(); Index: /trunk/ippconfig/gpc1/pswarp.config =================================================================== --- /trunk/ippconfig/gpc1/pswarp.config (revision 34799) +++ /trunk/ippconfig/gpc1/pswarp.config (revision 34800) @@ -1,1 +1,4 @@ ASTROM.SOURCE STR PSASTRO.OUTPUT # Source file rule for astrometry, or NULL +BACKGROUND.MODEL BOOL TRUE +BKG.XGRID S32 50 +BKG.YGRID S32 50 Index: /trunk/ippconfig/recipes/filerules-mef.mdc =================================================================== --- /trunk/ippconfig/recipes/filerules-mef.mdc (revision 34799) +++ /trunk/ippconfig/recipes/filerules-mef.mdc (revision 34800) @@ -120,4 +120,5 @@ PSWARP.SKYCELL INPUT @FILES CHIP IMAGE PSWARP.ASTROM INPUT @FILES CHIP CMF +PSWARP.BKGMODEL INPUT @FILES CHIP IMAGE ## files used by ppsub @@ -138,4 +139,5 @@ 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,4 +298,5 @@ 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,4 +345,6 @@ 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,4 +375,6 @@ 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 Index: /trunk/ippconfig/recipes/filerules-simple.mdc =================================================================== --- /trunk/ippconfig/recipes/filerules-simple.mdc (revision 34799) +++ /trunk/ippconfig/recipes/filerules-simple.mdc (revision 34800) @@ -90,4 +90,5 @@ PSWARP.SKYCELL INPUT @FILES FPA IMAGE PSWARP.ASTROM INPUT @FILES FPA CMF +PSWARP.BKGMODEL INPUT @FILES CHIP IMAGE ## files used by ppsub @@ -108,4 +109,5 @@ 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,4 +253,5 @@ 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,4 +303,6 @@ 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,4 +332,6 @@ 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 Index: /trunk/ippconfig/recipes/filerules-split.mdc =================================================================== --- /trunk/ippconfig/recipes/filerules-split.mdc (revision 34799) +++ /trunk/ippconfig/recipes/filerules-split.mdc (revision 34800) @@ -108,4 +108,5 @@ PSWARP.SKYCELL INPUT @FILES FPA IMAGE PSWARP.ASTROM INPUT @FILES CHIP CMF +PSWARP.BKGMODEL INPUT @FILES CHIP IMAGE ## files used by ppsub @@ -126,4 +127,5 @@ 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,4 +280,5 @@ 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,4 +336,6 @@ 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,6 +365,8 @@ 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 Index: /trunk/ippconfig/recipes/nightly_science.config =================================================================== --- /trunk/ippconfig/recipes/nightly_science.config (revision 34799) +++ /trunk/ippconfig/recipes/nightly_science.config (revision 34800) @@ -229,4 +229,5 @@ DIFFABLE BOOL TRUE REDUCTION STR SWEETSPOT + ADDITIONAL_STACK_LABEL STR ecliptic.rp # REDUCTION STR SSTF_4SIG # WARP S16 60 @@ -240,4 +241,5 @@ DIFFABLE BOOL TRUE REDUCTION STR SWEETSPOT + ADDITIONAL_STACK_LABEL STR ecliptic.rp # REDUCTION STR SSTF_4SIG END @@ -250,4 +252,5 @@ DIFFABLE BOOL TRUE REDUCTION STR SWEETSPOT + ADDITIONAL_STACK_LABEL STR ecliptic.rp # REDUCTION STR SSTF_4SIG END @@ -260,4 +263,5 @@ DIFFABLE BOOL TRUE REDUCTION STR SWEETSPOT + ADDITIONAL_STACK_LABEL STR ecliptic.rp # REDUCTION STR SSTF_4SIG END Index: /trunk/ippconfig/recipes/ppBackground.mdc =================================================================== --- /trunk/ippconfig/recipes/ppBackground.mdc (revision 34799) +++ /trunk/ippconfig/recipes/ppBackground.mdc (revision 34800) @@ -1,4 +1,18 @@ # 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 Index: /trunk/ippconfig/recipes/ppStack.config =================================================================== --- /trunk/ippconfig/recipes/ppStack.config (revision 34799) +++ /trunk/ippconfig/recipes/ppStack.config (revision 34800) @@ -103,4 +103,6 @@ STACK.TYPE STR DEEP_STACK +DEFAULT.EXPTIME F32 NAN + # Recipe overrides for STACK STACK METADATA Index: /trunk/ippconfig/recipes/pswarp.config =================================================================== --- /trunk/ippconfig/recipes/pswarp.config (revision 34799) +++ /trunk/ippconfig/recipes/pswarp.config (revision 34800) @@ -19,4 +19,7 @@ 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 Index: /trunk/ppBackground/src/ppBackground.h =================================================================== --- /trunk/ppBackground/src/ppBackground.h (revision 34799) +++ /trunk/ppBackground/src/ppBackground.h (revision 34800) @@ -37,4 +37,13 @@ ); +/// 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( Index: /trunk/ppBackground/src/ppBackgroundRestore.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundRestore.c (revision 34799) +++ /trunk/ppBackground/src/ppBackgroundRestore.c (revision 34800) @@ -5,4 +5,32 @@ #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,5 +55,18 @@ 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,5 +90,9 @@ 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]; + } } Index: /trunk/ppSkycell/src/ppSkycellCamera.c =================================================================== --- /trunk/ppSkycell/src/ppSkycellCamera.c (revision 34799) +++ /trunk/ppSkycell/src/ppSkycellCamera.c (revision 34800) @@ -103,11 +103,27 @@ 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; } Index: /trunk/ppSkycell/src/ppSkycellLoop.c =================================================================== --- /trunk/ppSkycell/src/ppSkycellLoop.c (revision 34799) +++ /trunk/ppSkycell/src/ppSkycellLoop.c (revision 34800) @@ -15,6 +15,6 @@ ) { - 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,5 +27,5 @@ { 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,6 +38,6 @@ 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,5 +132,6 @@ 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,5 +179,5 @@ target->data.S32[i] = numProj; psTrace("ppSkycell", 3, "Image %d uses new projection\n", i); - + psArrayAdd(regionHDUs,1,hdu); numProj++; } @@ -184,10 +185,10 @@ 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,7 +199,7 @@ 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,10 +209,12 @@ 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,19 +230,24 @@ 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,15 +264,14 @@ 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,4 +294,8 @@ 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,36 +324,161 @@ 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; + + } Index: /trunk/ppStack/src/Makefile.am =================================================================== --- /trunk/ppStack/src/Makefile.am (revision 34799) +++ /trunk/ppStack/src/Makefile.am (revision 34800) @@ -1,3 +1,3 @@ -bin_PROGRAMS = ppStack +bin_PROGRAMS = ppStack ppStackMedian if HAVE_SVNVERSION @@ -25,4 +25,7 @@ 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,4 +46,35 @@ 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 \ Index: /trunk/ppStack/src/ppStack.h =================================================================== --- /trunk/ppStack/src/ppStack.h (revision 34799) +++ /trunk/ppStack/src/ppStack.h (revision 34800) @@ -38,5 +38,8 @@ 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,4 +114,7 @@ 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 Index: /trunk/ppStack/src/ppStackArguments.c =================================================================== --- /trunk/ppStack/src/ppStackArguments.c (revision 34799) +++ /trunk/ppStack/src/ppStackArguments.c (revision 34800) @@ -324,4 +324,5 @@ } + psTrace("ppStack", 1, "Done parsing arguments\n"); Index: /trunk/ppStack/src/ppStackCamera.c =================================================================== --- /trunk/ppStack/src/ppStackCamera.c (revision 34799) +++ /trunk/ppStack/src/ppStackCamera.c (revision 34800) @@ -156,5 +156,6 @@ 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,4 +216,21 @@ } + // 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,4 +478,25 @@ 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,4 +535,6 @@ } + // Define output file here. + return true; } Index: /trunk/ppStack/src/ppStackCleanup.c =================================================================== --- /trunk/ppStack/src/ppStackCleanup.c (revision 34799) +++ /trunk/ppStack/src/ppStackCleanup.c (revision 34800) @@ -22,10 +22,24 @@ 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,7 +76,6 @@ bool ppStackCleanup (pmConfig *config, ppStackOptions *options) { - psExit exitValue = ppStackExitCode(PS_EXIT_SUCCESS); // Exit code - + // Ensure everything closes if (config) { @@ -72,4 +85,6 @@ 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."); Index: /trunk/ppStack/src/ppStackCombineAlternate.c =================================================================== --- /trunk/ppStack/src/ppStackCombineAlternate.c (revision 34800) +++ /trunk/ppStack/src/ppStackCombineAlternate.c (revision 34800) @@ -0,0 +1,115 @@ +#include "ppStack.h" + +// This is the doomsday switch. +// #define TESTING // Enable test output +bool ppStackCombineMedian(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +{ + psAssert(stack, "Require stack"); + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psTimerStart("PPSTACK_BKGMED"); + + + pmReadout *outRO = options->outRO; + + psArray *inputs = psArrayAlloc(options->num); + for (int i = 0; i < options->num; i++) { + ppStackFileActivationSingle(config, PPSTACK_FILES_MEDIAN_IN, true, i); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); + pmFPAview *view = ppStackFilesIterateDown(config); + pmReadout *ro = pmFPAviewThisReadout(view,file->fpa); + inputs->data[i] = ro; + pmFPAfileClose(file,view); + } + if (!pmStackSimpleMedianCombine(outRO,inputs)) { + psFree(inputs); + return(false); + } +#if 0 + if (!ppStackWriteImage("/tmp/test_forced.median.fits", + outRO->parent->parent->parent->hdu->header, + outRO->image, + config)) { + fprintf(stderr,"Failed to write image because fail.\n"); + } +#endif + for (int i = 0; i < options->num; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); + pmFPAview *view = ppStackFilesIterateDown(config); + bool success = pmFPAfileClose(file,view); + if (!success) { + psTrace("ppStack",5,"I failed at closing a file.\n"); + } + + psFitsClose(file->fits); + file->fits = NULL; + file->header = NULL; + file->state = PM_FPA_STATE_CLOSED; + file->wrote_phu = false; + ppStackFileActivationSingle(config, PPSTACK_FILES_MEDIAN_IN, false, i); + } + psFree(inputs); + outRO->data_exists = true; + outRO->parent->data_exists = true; + outRO->parent->parent->data_exists = true; + + return(true); +} + + +bool ppStackCombineBackground(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +{ + psAssert(stack, "Require stack"); + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psTimerStart("PPSTACK_BKGMED"); + + + pmReadout *bkgRO = options->bkgRO; + + psArray *inputs = psArrayAlloc(options->num); + for (int i = 0; i < options->num; i++) { + ppStackFileActivationSingle(config, PPSTACK_FILES_BKG, true, i); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", i); + pmFPAview *view = ppStackFilesIterateDown(config); + pmReadout *ro = pmFPAviewThisReadout(view,file->fpa); + inputs->data[i] = ro; + pmFPAfileClose(file,view); + } + if (!pmStackSimpleMedianCombine(bkgRO,inputs)) { + psFree(inputs); + return(false); + } +#if 0 + if (!ppStackWriteImage("/tmp/test_forced.bkgmdl.fits", + bkgRO->parent->parent->parent->hdu->header, + bkgRO->image, + config)) { + fprintf(stderr,"Failed to write image because fail.\n"); + } +#endif + for (int i = 0; i < options->num; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", i); + pmFPAview *view = ppStackFilesIterateDown(config); + bool success = pmFPAfileClose(file,view); + if (!success) { + psTrace("ppStack",5,"I failed at closing a file.\n"); + } + + psFitsClose(file->fits); + file->fits = NULL; + file->header = NULL; + file->state = PM_FPA_STATE_CLOSED; + file->wrote_phu = false; + ppStackFileActivationSingle(config, PPSTACK_FILES_BKG, false, i); + } + psFree(inputs); + bkgRO->data_exists = true; + bkgRO->parent->data_exists = true; + bkgRO->parent->parent->data_exists = true; + + return(true); +} + Index: /trunk/ppStack/src/ppStackCombinePrepare.c =================================================================== --- /trunk/ppStack/src/ppStackCombinePrepare.c (revision 34799) +++ /trunk/ppStack/src/ppStackCombinePrepare.c (revision 34800) @@ -1,5 +1,5 @@ #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,6 +27,46 @@ 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,7 +82,19 @@ } - 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; + } + + + } Index: /trunk/ppStack/src/ppStackFiles.c =================================================================== --- /trunk/ppStack/src/ppStackFiles.c (revision 34799) +++ /trunk/ppStack/src/ppStackFiles.c (revision 34800) @@ -14,4 +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,4 +29,5 @@ "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.EXPNUM", "PPSTACK.OUTPUT.EXPWT", "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", + "PPSTACK.OUTPUT.BKGMODEL", NULL }; /// Unconvolved stack files @@ -41,4 +52,7 @@ 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); Index: /trunk/ppStack/src/ppStackLoop.c =================================================================== --- /trunk/ppStack/src/ppStackLoop.c (revision 34799) +++ /trunk/ppStack/src/ppStackLoop.c (revision 34800) @@ -63,5 +63,5 @@ // 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,5 +84,5 @@ pmCellFreeData(options->cells->data[i]); } - psFree(stack); + // psFree(stack); // Pixel rejection @@ -142,4 +142,11 @@ } + // 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,5 +178,5 @@ return false; } - psFree(stack); + // psFree(stack); psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); ppStackMemDump("cleanup"); @@ -186,5 +193,5 @@ // 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."); Index: /trunk/ppStack/src/ppStackLoop.h =================================================================== --- /trunk/ppStack/src/ppStackLoop.h (revision 34799) +++ /trunk/ppStack/src/ppStackLoop.h (revision 34800) @@ -4,4 +4,9 @@ // 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,4 +37,5 @@ 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 Index: /trunk/ppStack/src/ppStackMedian.c =================================================================== --- /trunk/ppStack/src/ppStackMedian.c (revision 34800) +++ /trunk/ppStack/src/ppStackMedian.c (revision 34800) @@ -0,0 +1,51 @@ +#include "ppStack.h" + +int main(int argc, char *argv[]) +{ + psLibInit(NULL); + psTimerStart("PPSTACK"); + psTimerStart("PPSTACK_STEPS"); + + pmErrorRegister(); + ppStackErrorRegister(); + psphotErrorRegister(); + + ppStackOptions *options = NULL; // Options for stacking + + pmConfig *config = pmConfigRead(&argc, argv, PPSTACK_RECIPE); // Configuration + if (!config) { + ppStackCleanup(config, options); + } + + ppStackVersionPrint(); + + if (!pmModelClassInit()) { + psError(PPSTACK_ERR_PROG, false, "Unable to initialise model classes."); + ppStackCleanup(config, options); + } + + if (!psphotInit()) { + psError(PPSTACK_ERR_PROG, false, "Error initialising psphot."); + ppStackCleanup(config, options); + } + + if (!ppStackArgumentsSetup(argc, argv, config)) { + ppStackCleanup(config, options); + } + + if (!ppStackCamera(config)) { + ppStackCleanup(config, options); + } + + if (!ppStackArgumentsParse(config)) { + ppStackCleanup(config, options); + } + + options = ppStackOptionsAlloc(); + if (!ppStackMedianLoop(config, options)) { + ppStackCleanup(config, options); + } + + ppStackCleanup(config, options); +} + Index: /trunk/ppStack/src/ppStackMedianLoop.c =================================================================== --- /trunk/ppStack/src/ppStackMedianLoop.c (revision 34800) +++ /trunk/ppStack/src/ppStackMedianLoop.c (revision 34800) @@ -0,0 +1,204 @@ +#include "ppStack.h" + +// static functions are defined below +/* static int stackSummary(const ppStackOptions *options, const char *place); */ + +bool ppStackMedianLoop(pmConfig *config, ppStackOptions *options) +{ + assert(config); + + psTimerStart("PPSTACK_TOTAL"); + psTimerStart("PPSTACK_STEPS"); + + // Setup + psTrace("ppStack", 1, "Setup....\n"); + if (!ppStackSetup(options, config)) { + psError(psErrorCodeLast(), false, "Unable to setup."); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 0: Setup: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("setup"); + + // Duplicate code that is in ppStackConvolve + options->cells = psArrayAlloc(options->num); // Cells for convolved images --- a handle for reading again + pmFPAview *view = ppStackFilesIterateDown(config); + for (int i = 0; i < options->num; i++) { + ppStackFileActivationSingle(config, PPSTACK_FILES_MEDIAN_IN, true, i); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest + + // pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy + if (!view) { + return false; + } + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); // Input readout + pmCell *inCell = readout->parent; // Input cell */ + options->cells->data[i] = psMemIncrRefCounter(inCell); + + } + psFree(view); + + ppStackThreadData *stack = ppStackThreadDataSetup(options, config, options->convolve); + if (!stack) { + psError(psErrorCodeLast(), false, "Unable to initialise stack threads."); + return false; + } + + // Preparation for stacking + if (!ppStackCombinePrepare("PPSTACK.OUTPUT",NULL,NULL,PPSTACK_FILES_MEDIAN_OUT,stack,options,config)) { + psError(psErrorCodeLast(), false, "Unabel to combine data."); + psFree(stack); + return false; + } + + if (!ppStackCombineMedian(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unabel to combine data."); + psFree(stack); + return false; + } +/* for (int x = 0; x < options->outRO->image->numCols; x++) { */ +/* for (int y = 0; y < options->outRO->image->numRows; y++) { */ +/* options->outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; */ +/* } */ +/* } */ + psFree(options->outRO->mask); + // Do this before performing photometry so that the cmf header gets all of the information. + options->zpInput = NULL; + options->expTimeInput = NULL; + options->airmassInput = NULL; + if (!ppStackUpdateHeader(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to update header."); + psFree(stack); + return false; + } +/* // 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"); */ +/* if (!ppStackPhotometry(options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to perform photometry."); */ +/* return false; */ +/* } */ +/* psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS")); */ +/* ppStackMemDump("photometry"); */ + + // Create JPEGS +/* if (!ppStackJPEGs(stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to make jpegs."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +/* // Assemble Stats */ +/* if (!ppStackStats(stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to assemble statistics."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ + + // Clean up + + psTrace("ppStack", 2, "Cleaning up after combination....\n"); + if (!ppStackCleanupFiles(stack, options, config, PPSTACK_FILES_MEDIAN_OUT, PPSTACK_FILES_MEDIAN_OUT, false)) { + psError(psErrorCodeLast(), false, "Unable to clean up."); + psFree(stack); + return false; + } + // psFree(stack); + psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("cleanup"); + +/* // Unconvolved stack --- it's cheap to calculate, compared to everything else! */ +/* // XXX unconvolved stack is currently using the convolved mask! oops! */ +/* if (options->convolve) { */ +/* // Start threading */ +/* ppStackThreadData *stack = ppStackThreadDataSetup(options, config, false); */ +/* if (!stack) { */ +/* psError(psErrorCodeLast(), false, "Unable to initialise stack threads."); */ +/* return false; */ +/* } */ + +/* // Prepare for combination */ +/* if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", NULL, PPSTACK_FILES_UNCONV, */ +/* stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to prepare for combination."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ + +/* // generate the unconvolved stack. NOTE: this one must be normalized since the inputs have not been */ +/* psTrace("ppStack", 2, "Stack of unconvolved images....\n"); */ +/* if (!ppStackCombineFinal(stack, options->origCovars, options, config, false, true, false)) { */ +/* psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +/* psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Unconvolved Stack: %f sec", psTimerClear("PPSTACK_STEPS")); */ +/* ppStackMemDump("unconv"); */ + +/* // Update Header */ +/* if (!ppStackUpdateHeader(stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to update header."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +/* // Clean up unconvolved stack */ +/* psTrace("ppStack", 2, "Cleaning up after unconvolved stack....\n"); */ +/* if (!ppStackCleanupFiles(stack, options, config, PPSTACK_FILES_UNCONV, PPSTACK_FILES_NONE, false)) { */ +/* psError(psErrorCodeLast(), false, "Unable to clean up."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +/* psFree(stack); */ +/* } */ + psFree(options->cells); options->cells = NULL; + + // Finish up + psTrace("ppStack", 1, "Finishing up....\n"); + if (!ppStackFinish(options, config)) { + psError(psErrorCodeLast(), false, "Unable to finish up."); + return false; + } + ppStackMemDump("finish"); + + return true; +} + +/* /// Print a summary of the inputs, and return the number of good inputs */ +/* static int stackSummary(const ppStackOptions *options, const char *place) */ +/* { */ +/* int numGood = 0; // Number of good inputs */ +/* psString summary = NULL; // Summary of images */ +/* for (int i = 0; i < options->num; i++) { */ +/* psString reason = NULL; // Reason for rejecting */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] == 0) { */ +/* psStringAppend(&reason, " Good."); */ +/* numGood++; */ +/* } else { */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_CAL) { */ +/* psStringAppend(&reason, " Calibration failed."); */ +/* } */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_PSF) { */ +/* psStringAppend(&reason, " PSF measurement failed."); */ +/* } */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_MATCH) { */ +/* psStringAppend(&reason, " PSF matching failed."); */ +/* } */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_CHI2) { */ +/* psStringAppend(&reason, " PSF matching chi^2 deviant."); */ +/* } */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_REJECT) { */ +/* psStringAppend(&reason, " Rejection exceeded threshold."); */ +/* } */ +/* } */ +/* psStringAppend(&summary, "Image %d: %s\n", i, reason); */ +/* psFree(reason); */ +/* } */ +/* psLogMsg("ppStack", PS_LOG_INFO, "Summary of images for %s:\n%s", place, summary); */ +/* psFree(summary); */ + +/* return numGood; */ +/* } */ Index: /trunk/ppStack/src/ppStackOptions.h =================================================================== --- /trunk/ppStack/src/ppStackOptions.h (revision 34799) +++ /trunk/ppStack/src/ppStackOptions.h (revision 34800) @@ -48,4 +48,7 @@ // Rejection psArray *rejected; // Rejected pixels + // Background + pmReadout *bkgRO; // Output background readout + psArray *bkgImages; // Input background images } ppStackOptions; Index: /trunk/ppStack/src/ppStackPrepare.c =================================================================== --- /trunk/ppStack/src/ppStackPrepare.c (revision 34799) +++ /trunk/ppStack/src/ppStackPrepare.c (revision 34800) @@ -132,4 +132,5 @@ psVectorInit(options->inputMask, 0); options->exposures = psVectorAlloc(options->num, PS_TYPE_F32); + psVectorInit(options->exposures, NAN); @@ -150,4 +151,8 @@ 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]; Index: /trunk/ppStack/src/ppStackReadout.c =================================================================== --- /trunk/ppStack/src/ppStackReadout.c (revision 34799) +++ /trunk/ppStack/src/ppStackReadout.c (revision 34800) @@ -194,5 +194,4 @@ - bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, const psArray *readouts, const psVector *mask, const psArray *rejected, const psVector *weightings, Index: /trunk/ppStack/src/ppStackSetup.c =================================================================== --- /trunk/ppStack/src/ppStackSetup.c (revision 34799) +++ /trunk/ppStack/src/ppStackSetup.c (revision 34800) @@ -89,5 +89,8 @@ 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,11 +100,34 @@ { // 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); Index: /trunk/ppStack/src/ppStackThread.c =================================================================== --- /trunk/ppStack/src/ppStackThread.c (revision 34799) +++ /trunk/ppStack/src/ppStackThread.c (revision 34800) @@ -40,4 +40,5 @@ psFree(stack->maskFits); psFree(stack->varianceFits); + psFree(stack->bkgFits); return; } @@ -67,5 +68,5 @@ PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL); } - + ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return psMemSetDeallocator(stack, (psFreeFunc)stackThreadDataFree); @@ -77,4 +78,5 @@ stack->maskFits = psArrayAlloc(numInputs); stack->varianceFits = psArrayAlloc(numInputs); + stack->bkgFits = psArrayAlloc(numInputs); for (int i = 0; i < numInputs; i++) { if (!cells->data[i]) { @@ -95,8 +97,12 @@ 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); + } } Index: /trunk/ppStack/src/ppStackThread.h =================================================================== --- /trunk/ppStack/src/ppStackThread.h (revision 34799) +++ /trunk/ppStack/src/ppStackThread.h (revision 34800) @@ -25,4 +25,5 @@ psArray *maskFits; // FITS file pointers for masks psArray *varianceFits; // FITS file pointers for variances + psArray *bkgFits; // FITS file pointers for background models } ppStackThreadData; Index: /trunk/ppStack/src/ppStackUpdateHeader.c =================================================================== --- /trunk/ppStack/src/ppStackUpdateHeader.c (revision 34799) +++ /trunk/ppStack/src/ppStackUpdateHeader.c (revision 34800) @@ -6,9 +6,10 @@ 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,5 +124,23 @@ 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; } Index: /trunk/ppViz/src/ppVizPSF/ppVizPSFCamera.c =================================================================== --- /trunk/ppViz/src/ppVizPSF/ppVizPSFCamera.c (revision 34799) +++ /trunk/ppViz/src/ppVizPSF/ppVizPSFCamera.c (revision 34800) @@ -19,4 +19,5 @@ files->data[0] = psStringCopy(name); if (psMetadataLookup(config->arguments, file)) { + psMetadataRemoveKey(config->arguments, file); } @@ -42,5 +43,6 @@ 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"); Index: /trunk/ppViz/src/ppVizPSF/ppVizPSFLoop.c =================================================================== --- /trunk/ppViz/src/ppVizPSF/ppVizPSFLoop.c (revision 34799) +++ /trunk/ppViz/src/ppVizPSF/ppVizPSFLoop.c (revision 34800) @@ -20,5 +20,5 @@ return NULL; } - + fprintf(stderr,"Woo!\n"); pmChip *chip; // Chip from FPA while ((chip = pmFPAviewNextChip(view, psfFile->fpa, 1))) { @@ -56,15 +56,13 @@ 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,4 +70,12 @@ 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,5 +86,5 @@ 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,4 +95,5 @@ 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,4 +125,5 @@ } if (!sources && !data->input) { + fprintf(stderr,"Here! default\n"); // Generate fake image with only a single realisation of the PSF sources = psArrayAlloc(1); @@ -151,6 +159,7 @@ 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,7 +168,9 @@ } } 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,5 +179,5 @@ } } - + psFree(sources); psFree(xOffset); Index: /trunk/psLib/src/imageops/psImageInterpolate.c =================================================================== --- /trunk/psLib/src/imageops/psImageInterpolate.c (revision 34799) +++ /trunk/psLib/src/imageops/psImageInterpolate.c (revision 34800) @@ -31,5 +31,5 @@ #include "psImageConvolve.h" -# define IS_BILIN_SEPARABLE 0 +# define IS_BILIN_SEPARABLE 1 #include "psImageInterpolate.h" @@ -189,4 +189,5 @@ switch (mode) { case PS_INTERPOLATE_FLAT: + case PS_INTERPOLATE_BILINEAR_SIMPLE: // Nothing to pre-compute break; @@ -287,6 +288,8 @@ # 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,4 +940,59 @@ +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,7 +1027,9 @@ 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,4 +1079,5 @@ switch (mode) { case PS_INTERPOLATE_FLAT: + case PS_INTERPOLATE_BILINEAR_SIMPLE: // No smearing by design return 1.0; @@ -1083,4 +1144,5 @@ 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; Index: /trunk/psLib/src/imageops/psImageInterpolate.h =================================================================== --- /trunk/psLib/src/imageops/psImageInterpolate.h (revision 34799) +++ /trunk/psLib/src/imageops/psImageInterpolate.h (revision 34800) @@ -33,4 +33,5 @@ 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; Index: /trunk/psLib/src/imageops/psImagePixelManip.c =================================================================== --- /trunk/psLib/src/imageops/psImagePixelManip.c (revision 34799) +++ /trunk/psLib/src/imageops/psImagePixelManip.c (revision 34800) @@ -216,15 +216,54 @@ } - +#define psImageOverlayLoopClean(DATATYPE,OP) { \ + for (int row=y0;rowdata.DATATYPE[row]; \ + ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-y0]; \ + for (int col=x0;coldata.DATATYPE[row]; \ + ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-y0]; \ + for (int col=x0;coldata.DATATYPE[row]; \ ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-y0]; \ - for (int col=x0;coldata.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,7 +84,10 @@ sum += inImage->data.F32[y][x]; numPix++; + + } } - + + // Values to set float imageValue; @@ -79,5 +101,11 @@ } 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; } Index: /trunk/psModules/src/camera/pmFPAfileIO.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfileIO.c (revision 34799) +++ /trunk/psModules/src/camera/pmFPAfileIO.c (revision 34800) @@ -76,5 +76,4 @@ while ((item = psMetadataGetAndIncrement (iter)) != NULL) { pmFPAfile *file = item->data.V; - switch (place) { case PM_FPA_BEFORE: Index: /trunk/psModules/src/detrend/pmDark.c =================================================================== --- /trunk/psModules/src/detrend/pmDark.c (revision 34799) +++ /trunk/psModules/src/detrend/pmDark.c (revision 34800) @@ -353,12 +353,35 @@ // 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,5 +406,5 @@ } - pmDarkVisualInit(values); + pmDarkVisualInit(values_set->data[max_orders->n - 1]); pmReadout *outReadout = output->readouts->data[0]; @@ -423,15 +446,53 @@ } - 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; Index: /trunk/psModules/src/imcombine/pmStack.c =================================================================== --- /trunk/psModules/src/imcombine/pmStack.c (revision 34799) +++ /trunk/psModules/src/imcombine/pmStack.c (revision 34800) @@ -1358,4 +1358,57 @@ } +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( Index: /trunk/psModules/src/imcombine/pmStack.h =================================================================== --- /trunk/psModules/src/imcombine/pmStack.h (revision 34799) +++ /trunk/psModules/src/imcombine/pmStack.h (revision 34800) @@ -41,4 +41,8 @@ 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 Index: /trunk/psModules/src/objects/Makefile.am =================================================================== --- /trunk/psModules/src/objects/Makefile.am (revision 34799) +++ /trunk/psModules/src/objects/Makefile.am (revision 34800) @@ -111,4 +111,5 @@ pmSourceOutputs.h \ pmSourceIO.h \ + pmSourceSatstar.h \ pmSourcePlots.h \ pmSourceVisual.h \ Index: /trunk/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- /trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 34799) +++ /trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 34800) @@ -91,5 +91,16 @@ // 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; Index: /trunk/pswarp/src/pswarp.c =================================================================== --- /trunk/pswarp/src/pswarp.c (revision 34799) +++ /trunk/pswarp/src/pswarp.c (revision 34800) @@ -46,4 +46,9 @@ goto DIE; } + if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) { + if (!pswarpDefineBackground(config)) { + goto DIE; + } + } // Open the statistics file @@ -66,4 +71,10 @@ if (!pswarpLoop(config, stats)) { goto DIE; + } + if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) { + if (!pswarpLoopBackground(config, stats)) { + fprintf(stderr,"Dying!\n"); + goto DIE; + } } Index: /trunk/pswarp/src/pswarp.h =================================================================== --- /trunk/pswarp/src/pswarp.h (revision 34799) +++ /trunk/pswarp/src/pswarp.h (revision 34800) @@ -62,4 +62,9 @@ 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,5 +86,7 @@ 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); Index: /trunk/pswarp/src/pswarpArguments.c =================================================================== --- /trunk/pswarp/src/pswarpArguments.c (revision 34799) +++ /trunk/pswarp/src/pswarpArguments.c (revision 34800) @@ -86,5 +86,6 @@ 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,4 +165,21 @@ } + 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,5 +194,9 @@ "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,4 +211,7 @@ "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"); Index: /trunk/pswarp/src/pswarpDefine.c =================================================================== --- /trunk/pswarp/src/pswarpDefine.c (revision 34799) +++ /trunk/pswarp/src/pswarpDefine.c (revision 34800) @@ -41,5 +41,5 @@ return false; } - + // open the full skycell file; no need to defer different depths. only load the header data pmFPAview *view = pmFPAviewAlloc (0); @@ -86,4 +86,6 @@ psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0"); psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0"); + + } @@ -123,9 +125,153 @@ } } - + 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; +} Index: /trunk/pswarp/src/pswarpFileNames.h =================================================================== --- /trunk/pswarp/src/pswarpFileNames.h (revision 34799) +++ /trunk/pswarp/src/pswarpFileNames.h (revision 34800) @@ -12,4 +12,5 @@ "PSWARP.MASK", "PSWARP.VARIANCE", + "PSWARP.BKGMODEL", NULL }; @@ -20,4 +21,5 @@ "PSWARP.OUTPUT.MASK", "PSWARP.OUTPUT.VARIANCE", + "PSWARP.OUTPUT.BKGMODEL", NULL }; Index: /trunk/pswarp/src/pswarpLoop.c =================================================================== --- /trunk/pswarp/src/pswarpLoop.c (revision 34799) +++ /trunk/pswarp/src/pswarpLoop.c (revision 34800) @@ -179,4 +179,6 @@ // read WCS data from the corresponding header pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa); + + if (bilevelAstrometry) { if (!pmAstromReadBilevelChip (chip, hdu->header)) { @@ -195,5 +197,5 @@ } } - + pmCell *cell; while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { @@ -227,5 +229,5 @@ pswarpTransformReadout(output, readout, config); - + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { psError(psErrorCodeLast(), false, "Unable to write files."); @@ -255,5 +257,5 @@ goto DONE; } - + pmCell *outCell = output->parent; ///< Output cell pmChip *outChip = outCell->parent; ///< Output chip @@ -450,2 +452,354 @@ 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; +} Index: /trunk/pswarp/src/pswarpParseCamera.c =================================================================== --- /trunk/pswarp/src/pswarpParseCamera.c (revision 34799) +++ /trunk/pswarp/src/pswarpParseCamera.c (revision 34800) @@ -84,4 +84,10 @@ } + 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,4 +140,20 @@ } + 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 Index: /trunk/pswarp/src/pswarpTransformReadout.c =================================================================== --- /trunk/pswarp/src/pswarpTransformReadout.c (revision 34799) +++ /trunk/pswarp/src/pswarpTransformReadout.c (revision 34800) @@ -29,4 +29,5 @@ psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.MODE"); ///< Mode for interp + int numKernels = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.NUM"); ///< Number of kernels @@ -60,5 +61,6 @@ // 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,4 +132,11 @@ 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,5 +180,4 @@ } 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,4 +191,6 @@ jacobian += args->jacobian * args->goodPixels; } + + } psFree(job); @@ -201,8 +211,10 @@ 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; } + } } Index: /trunk/pswarp/src/pswarpTransformSources.c =================================================================== --- /trunk/pswarp/src/pswarpTransformSources.c (revision 34799) +++ /trunk/pswarp/src/pswarpTransformSources.c (revision 34800) @@ -41,5 +41,5 @@ 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; Index: /trunk/pswarp/src/pswarpTransformTile.c =================================================================== --- /trunk/pswarp/src/pswarpTransformTile.c (revision 34799) +++ /trunk/pswarp/src/pswarpTransformTile.c (revision 34800) @@ -46,4 +46,8 @@ args->jacobian = NAN; + args->background_warping = false; + args->offset_x = 0.0; + args->offset_y = 0.0; + return args; } @@ -83,16 +87,32 @@ // 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;