Index: /branches/pap_branch_080320/ippScripts/scripts/detrend_norm_calc.pl
===================================================================
--- /branches/pap_branch_080320/ippScripts/scripts/detrend_norm_calc.pl	(revision 17139)
+++ /branches/pap_branch_080320/ippScripts/scripts/detrend_norm_calc.pl	(revision 17139)
@@ -0,0 +1,233 @@
+#!/usr/bin/env perl
+
+use Carp;
+use warnings;
+use strict;
+
+## report the program and machine
+use Sys::Hostname;
+my $host = hostname();
+print "\n\n";
+print "Starting script $0 on $host\n\n";
+
+use vars qw( $VERSION );
+$VERSION = '0.01';
+
+use IPC::Cmd 0.36 qw( can_run );
+use IPC::Run 0.36 qw( run );
+use PS::IPP::Metadata::Config;
+use PS::IPP::Metadata::List qw( parse_md_list );
+
+use PS::IPP::Config qw($PS_EXIT_SUCCESS
+		       $PS_EXIT_UNKNOWN_ERROR
+		       $PS_EXIT_SYS_ERROR
+		       $PS_EXIT_CONFIG_ERROR
+		       $PS_EXIT_PROG_ERROR
+		       $PS_EXIT_DATA_ERROR
+		       $PS_EXIT_TIMEOUT_ERROR
+		       );
+
+use Getopt::Long qw( GetOptions :config auto_help auto_version gnu_getopt );
+use Pod::Usage qw( pod2usage );
+
+# Parse command-line arguments
+my ($det_id, $iter, $detType, $outroot, $dbname, $verbose, $no_update, $no_op);
+GetOptions(
+	'det_id|d=s'	=> \$det_id,    # Detrend id			     
+	'iteration|i=s'	=> \$iter,	# Iteration			     
+	'det_type|t=s'  => \$detType,	# Detrend type			     
+        'outroot|w=s'   => \$outroot,   # output file base name
+ 	'dbname|d=s'    => \$dbname,	# Database name			     
+        'verbose'       => \$verbose,   # Print to stdout
+        'no-update'     => \$no_update,	# Don't update the database?	     
+	'no-op'         => \$no_op,	# Don't do operations                
+	) or pod2usage( 2 );
+
+pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV;
+pod2usage( -msg => "Required options --det_id --iteration --det_type --outroot",
+	   -exitval => 3,
+	   ) unless 
+    defined $det_id  and 
+    defined $iter    and 
+    defined $detType and
+    defined $outroot;
+
+use constant STATISTIC => 'bg'; # Background statistic to use from the database
+
+# Define which detrend types we normalise
+use constant NORMALIZE => {
+    'bias'     => 0,
+    'dark'     => 0,
+    'shutter'  => 0,
+    'flat'     => 1,
+    'domeflat' => 1,
+    'skyflat'  => 1,
+    'fringe'   => 0,
+    'mask'     => 0,
+    'darkmask' => 0,
+    'flatmask' => 0,
+    }; 
+
+# Look for programs we need
+my $missing_tools;
+my $dettool = can_run('dettool') or (warn "Can't find dettool" and $missing_tools = 1);
+my $ppNormCalc = can_run('ppNormCalc') or (warn "Can't find ppNormCalc" and $missing_tools = 1);
+if ($missing_tools) { 
+    warn("Can't find required tools.");
+    exit($PS_EXIT_CONFIG_ERROR); 
+}
+
+&my_die("Unrecognised detrend type: $detType", $det_id, $iter, $PS_EXIT_PROG_ERROR) unless exists NORMALIZE()->{lc($detType)};
+
+my $mdcParser = PS::IPP::Metadata::Config->new; # Parser for metadata config files
+
+# Get the list of inputs
+my @files;			# The input files
+{
+    my $command = "$dettool -processedimfile -det_id $det_id"; # Command to run
+    $command .= " -dbname $dbname" if defined $dbname;
+    my @command = split /\s+/, $command;
+    my ( $stdin, $stdout, $stderr ); # Buffers for running program
+    print "Running [$command]...\n" if $verbose;
+    if (not run(\@command, \$stdin, \$stdout, \$stderr)) {
+	&my_die("Unable to perform dettool -processedimfile on detrend $det_id/$iter: $?",
+		$det_id, $iter, $PS_EXIT_SYS_ERROR);
+    }
+    print $stdout . "\n" if $verbose;
+    
+    # Because of the length, need to split into individual metadatas --- it parses SO much quicker!
+    my @whole = split /\n/, $stdout;
+    my @single = ();
+    while ( scalar @whole > 0 ) {
+	my $value = shift @whole;
+	push @single, $value;
+	if ($value =~ /^\s*END\s*$/) {
+	    push @single, "\n";
+
+	    my $list = parse_md_list( $mdcParser->parse( join( "\n", @single ) ) );
+	    &my_die("Unable to parse output from dettool", $det_id, $iter, $PS_EXIT_PROG_ERROR) unless
+		defined $list;
+	    push @files, @$list;
+	    @single = ();
+	}
+    }
+}
+
+my $norms;			# MDC with normalisations
+if (NORMALIZE()->{lc($detType)} and not $no_op) {
+
+    my %matrix; # Matrix of statistics as a function of exposures and classes
+    foreach my $file (@files) {
+	my $exp_id = $file->{'exp_id'}; # Exposure ID
+	my $class_id = $file->{'class_id'}; # Class ID
+	my $stat = $file->{STATISTIC()}; # Statistic of interest
+	
+	# Create matrix elements
+	$matrix{$exp_id} = {} if not defined $matrix{$exp_id};
+	$matrix{$exp_id}->{$class_id} = $stat;
+    }
+    
+    # Generate the input for ppNormCalc
+    my $normData;			# Normalisation data
+    foreach my $exp_id (keys %matrix) {
+	$normData .= "$exp_id\tMETADATA\n";
+	foreach my $class_id (keys %{$matrix{$exp_id}}) {
+	    $normData .= "\t" . $class_id . "\tF32\t" . $matrix{$exp_id}->{$class_id} . "\n";
+	}
+	$normData .= "END\n\n";
+    }
+
+    # Run ppNormCalc
+    {
+	my ( $stdout, $stderr ); # Buffers for running program
+	my @command = split /\s+/, $ppNormCalc;
+	print "Running [$ppNormCalc]...\n" if $verbose;
+	if (not run(\@command, \$normData, \$stdout, \$stderr)) {
+	    &my_die("Unable to perform ppNormCalc: $?", $det_id, $iter, $PS_EXIT_SYS_ERROR);
+	}
+	print $stdout . "\n" if $verbose;
+	
+	# Parse the output
+	$norms = $mdcParser->parse($stdout);
+	&my_die("Unable to parse metadata config doc", $det_id, $iter, $PS_EXIT_PROG_ERROR) unless $norms;
+    }
+
+} else {
+    # It's something that doesn't need normalisation --- just push in a normalisation of 1
+    my %classes;		# List of unique classes
+    foreach my $file (@files) {
+	my $class_id = $file->{'class_id'}; # Class Id
+	$classes{$class_id} = 1;
+    }
+    
+    foreach my $class_id (keys %classes) {
+	my %mdValue;	# Metadata value for this class id
+	$mdValue{name} = $class_id;
+	$mdValue{value} = 1.0;
+	push @$norms, \%mdValue;
+    }
+}
+
+my $commandBase = "$dettool -addnormalizedstat";
+$commandBase .= " -det_id $det_id";
+$commandBase .= " -iteration $iter";
+$commandBase .= " -dbname $dbname" if defined $dbname;
+
+# Process output normalisations
+unless ($no_update) {
+    foreach my $normItem (@$norms) {
+
+	my $className = $normItem->{name}; # Name of component
+	my $normalisation = $normItem->{value}; # Normalisation for component
+
+	if ($normalisation == 0.0 or lc($normalisation) eq 'nan') {
+	    warn("Class $className has bad normalisation: $normalisation");
+	    exit($PS_EXIT_SYS_ERROR);
+	}
+
+	my $command = $commandBase;
+	$command .= " -class_id $className";
+	$command .= " -norm $normalisation";
+
+	my @command = split /\s+/, $command;
+
+	my ( $stdin, $stdout, $stderr ); # Buffers for running program
+	print "Running [$command]...\n" if $verbose;
+	if (not run \@command, \$stdin, \$stdout, \$stderr) {
+	    warn("Unable to perform dettool -addnormstat for $className: $?");
+	    exit($PS_EXIT_SYS_ERROR);
+	}
+	print $stdout . "\n" if $verbose;
+    }
+} else {
+    print "skipping command: $commandBase\n";
+}
+
+sub my_die
+{
+    my $msg = shift;		# Warning message on die
+    my $det_id = shift;		# Detrend identifier
+    my $iter = shift;		# Iteration
+    my $exit_code = shift;	# Exit code to add
+
+    carp($msg);
+    if (defined $det_id and defined $iter and not $no_update) {
+	my $command = "$dettool -addnormalizedstat";
+	$command .= " -det_id $det_id";
+	$command .= " -iteration $iter";
+	$command .= " -code $exit_code";
+	$command .= " -dbname $dbname" if defined $dbname;
+        system ($command);
+    }
+    exit $exit_code;
+}
+
+
+END {
+    my $status = $?;
+    system("sync") == 0
+        or die "failed to execute sync: $!" ;
+    $? = $status;
+}
+
+__END__
Index: /branches/pap_branch_080320/ippScripts/scripts/ipp_detrend_combine.pl
===================================================================
--- /branches/pap_branch_080320/ippScripts/scripts/ipp_detrend_combine.pl	(revision 17139)
+++ /branches/pap_branch_080320/ippScripts/scripts/ipp_detrend_combine.pl	(revision 17139)
@@ -0,0 +1,228 @@
+#!/usr/bin/env perl
+
+use warnings;
+use strict;
+
+## report the program and machine
+use Sys::Hostname;
+my $host = hostname();
+print "\n\n";
+print "Starting script $0 on $host\n\n";
+
+use vars qw( $VERSION );
+$VERSION = '0.01';
+
+use Data::Dumper;
+use IPC::Cmd 0.36 qw( can_run run );
+use PS::IPP::Metadata::Config;
+use PS::IPP::Metadata::List qw( parse_md_list );
+
+use PS::IPP::Config qw($PS_EXIT_SUCCESS
+		       $PS_EXIT_UNKNOWN_ERROR
+		       $PS_EXIT_SYS_ERROR
+		       $PS_EXIT_CONFIG_ERROR
+		       $PS_EXIT_PROG_ERROR
+		       $PS_EXIT_DATA_ERROR
+		       $PS_EXIT_TIMEOUT_ERROR
+		       caturi
+		       );
+my $ipprc = PS::IPP::Config->new(); # IPP configuration
+
+use PS::IPP::Metadata::Stats;
+use Getopt::Long qw( GetOptions :config auto_help auto_version gnu_getopt );
+use Pod::Usage qw( pod2usage );
+
+my $RECIPE_PPSTATS = 'CHIPSTATS'; # Recipe to use with ppStats
+
+# Parse command-line arguments
+my ($det_type, $filelevel, $inst, $telescope, $filter,
+    $det_id1, $iter1, $det_id2, $iter2, $operation, $mask,
+    $workdir, $dbname, $no_update);
+GetOptions(
+	   'det_type=s'    => \$det_type, # Detrend type for new detrend
+	   'filelevel=s'   => \$filelevel, # File level for new detrend
+	   'inst=s'        => \$inst, # Instrument for new detrend
+	   'telescope=s'   => \$telescope, # Telescope for new detrend
+	   'filter=s'      => \$filter,	# Filter name for new detrend
+	   'det_id1=s'	   => \$det_id1, # Detrend id for detrend 1
+	   'iteration1=s'  => \$iter1, # Iteration for detrend 1
+	   'det_id2=s'	   => \$det_id2, # Detrend id for detrend 2
+	   'iteration2=s'  => \$iter2, # Iteration for detrend 2
+	   'operation=s'   => \$operation, # Operation to perform on files
+	   'mask'          => \$mask, # Operation is on a mask
+	   'workdir=s'     => \$workdir, # Working directory for output files
+	   'dbname=s'      => \$dbname,	# Database name
+	   'no-update'     => \$no_update, # Don't update the database
+	   ) or pod2usage( 2 );
+
+pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV;
+pod2usage( -msg => "Required options --det_type --filelevel --inst --telescope --det_id1 --iteration1 --det_id2 --iteration2 --workdir",
+	   -exitval => 3,
+	   )
+    unless defined $det_type
+    and defined $filelevel
+    and defined $inst
+    and defined $telescope
+    and defined $det_id1
+    and defined $iter1
+    and defined $det_id2
+    and defined $iter2
+    and defined $operation
+    and defined $workdir;
+
+$ipprc->define_camera($inst);
+
+my $STATS = 
+   [   
+       #          PPSTATS KEYWORD         STATISTIC          CHIPTOOL FLAG
+       { name => "ROBUST_MEDIAN",  type => "mean",  flag => "-bg",             dtype => "float" },
+       { name => "ROBUST_MEDIAN",  type => "stdev", flag => "-bg_mean_stdev",  dtype => "float" },
+       { name => "ROBUST_STDEV",   type => "rms",   flag => "-bg_stdev",       dtype => "float" },
+   ];
+
+# Look for programs we need
+my $missing_tools;
+my $detselect = can_run('detselect') or (warn "Can't find detselect" and $missing_tools = 1);
+my $dettool = can_run('dettool') or (warn "Can't find dettool" and $missing_tools = 1);
+my $ppArith = can_run('ppArith') or (warn "Can't find ppArith" and $missing_tools = 1);
+if ($missing_tools) {
+    warn("Can't find required tools.");
+    exit($PS_EXIT_CONFIG_ERROR);
+}
+
+my $mdcParser = PS::IPP::Metadata::Config->new; # Parser for metadata config files
+
+# Get the list of inputs
+my $files1 = filelist($det_id1, $iter1); # Hash of input files for detrend 1
+my $files2 = filelist($det_id2, $iter2); # Hash of input files for detrend 2
+die("File lists for detrends have differing lengths") unless scalar keys %$files1 == scalar keys %$files2;
+
+my ($det_id, $iter);	      # Detrend identifier for the new detrend
+unless ($no_update) {
+    my $command = "$dettool -register_detrend -det_type $det_type -filelevel $filelevel -workdir $workdir " .
+	"-inst $inst -telescope $telescope"; # Command to run
+    $command .= " -filter $filter" if defined $filter;
+    $command .= " -dbname $dbname" if defined $dbname;
+    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
+	run(command => $command, verbose => 1);
+    unless ($success) {
+	$error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
+	die("Unable to run dettool -register_detrend: $error_code");
+    }
+
+    my $metadata = $mdcParser->parse(join "", @$stdout_buf) or die("Unable to parse metadata config doc\n");
+    my $md = parse_md_list($metadata) or die("Unable to parse metadata list\n");
+
+    $det_id = $$md[0]->{det_id};
+    $iter = $$md[0]->{iteration};
+
+    die("Unable to get det_id and iteration for new detrend.\n") unless defined $det_id and defined $iter;
+} else {
+    $det_id = 'DUMMY_DET_ID';
+    $iter = 'DUMMY_ITER';
+}
+
+my $outRoot = caturi($workdir, "$inst.$det_id.$iter"); # Output root name
+my $filerule = (defined $mask ? "PPARITH.OUTPUT.MASK" : "PPARITH.OUTPUT.IMAGE"); # File rule for ppArith
+
+foreach my $class_id ( keys %$files1 ) {
+    my $md1 = $$files1{$class_id};
+    my $md2 = $$files2{$class_id};
+    die("Class_id=$class_id not defined for det_id=$det_id2") unless defined $md2;
+
+    my $uri1 = $$md1[0]->{uri};
+    my $uri2 = $$md2[0]->{uri};
+
+    die("Unable to find input file $uri1\n") unless $ipprc->file_exists($uri1);
+    die("Unable to find input file $uri2\n") unless $ipprc->file_exists($uri2);
+
+    my $outName = $ipprc->filename($filerule, $outRoot, $class_id);
+    my $outStats = $outRoot . '.stats';
+
+    my $command = "$ppArith -file1 $uri1 -op \'$operation\' -file2 $uri2 $outRoot"; # Command to run
+    $command .= " -stats $outStats -recipe PPSTATS $RECIPE_PPSTATS";
+    $command .= ' -mask' if defined $mask;
+    $command .= " -dbname $dbname" if defined $dbname;
+    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
+	run(command => $command, verbose => 1);
+    unless ($success) {
+	$error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
+	die("Unable to run ppArith: $error_code");
+    }
+
+    die("Unable to find ppArith product: $outName\n") unless $ipprc->file_exists($outName);
+    die("Unable to find ppArith product: $outStats\n") unless $ipprc->file_exists($outStats);
+
+    # Get the statistics on the processed image
+    my $stats = PS::IPP::Metadata::Stats->new($STATS); # Stats parser
+    {
+	my $statsFile;		# File handle
+	open $statsFile, $ipprc->file_resolve($outStats) or die("Can't open stats file $outStats: $!");
+	my @contents = <$statsFile>; # Contents of file
+	close $statsFile;
+	
+	my $metadata = $mdcParser->parse(join "", @contents) or die("Unable to parse metadata config doc");
+
+	unless ($stats->parse($metadata)) {
+	    &my_die("Failure extracting metadata from the statistics output file.\n");
+	}
+    }
+
+    # Register the imfile
+    unless ($no_update) {
+	my $command = "$dettool -register_detrend_imfile -det_id $det_id "; # Command to run
+	$command .= " -class_id $class_id -uri $outName -path_base $outRoot";
+	$command .= $stats->cmdflags();
+	$command .= " -dbname $dbname" if defined $dbname;
+	my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
+	    run(command => $command, verbose => 1);
+	unless ($success) {
+	    $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
+	    die("Unable to run dettool -register_detrend_imfile: $error_code");
+	}
+    }
+}
+
+
+### Pau.
+
+
+# Get a list of files for the given detrend
+sub filelist
+{
+    my $det_id = shift;		# Detrend identifier
+    my $iter = shift;		# Iteration
+
+    my $command = "$detselect -select -det_id $det_id -iteration $iter"; # Command to run
+    $command .= " -dbname $dbname" if defined $dbname;
+    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
+	run(command => $command, verbose => 1);
+    unless ($success) {
+	$error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
+	die("Unable to run detselect: $error_code");
+    }
+
+    # Because of the length, need to split into individual metadatas --- it parses SO much quicker!
+    my %files;
+
+    my $md = $mdcParser->parse( join( "", @$stdout_buf ) ); # Parsed metadata
+    my $list = parse_md_list( $md );
+
+    foreach my $item ( @$list ) {
+	my $class_id = $item->{class_id};
+	die("Multiple definitions of class_id=$class_id found for det_id=$det_id, iteration=$iter\n") if
+	    defined $files{$class_id};
+	$files{$class_id} = parse_md_list($md);
+    }
+
+    return \%files;
+}
+
+END {
+    my $status = $?;
+    system("sync") == 0
+        or die "failed to execute sync: $!" ;
+    $? = $status;
+}
+
+__END__
Index: /branches/pap_branch_080320/ippconfig/gpc1/filerules.mdc
===================================================================
--- /branches/pap_branch_080320/ippconfig/gpc1/filerules.mdc	(revision 17139)
+++ /branches/pap_branch_080320/ippconfig/gpc1/filerules.mdc	(revision 17139)
@@ -0,0 +1,138 @@
+PSASTRO.INPUT      STR PSASTRO.INPUT.CMF
+PSASTRO.OUTPUT     STR PSASTRO.OUT.CMF.SPL
+PSASTRO.OUTPUT.MEF STR PSASTRO.OUT.CMF.MEF
+PSPHOT.OUTPUT      STR PSPHOT.OUT.CMF.SPL
+
+### input file definitions
+TYPE                    INPUT    FILENAME.RULE                 DATA.LEVEL FILE.TYPE
+PPIMAGE.INPUT           INPUT    none.fits                     CHIP       IMAGE
+INPUT.MASK              INPUT    none.fits                     CHIP       MASK
+INPUT.WEIGHT            INPUT    none.fits                     CHIP       WEIGHT
+INPUT.PSF               INPUT    none.fits                     CHIP       PSF
+INPUT.SRC               INPUT    none.fits                     CHIP       CMF
+PPIMAGE.BIAS            INPUT    @DETDB                        CHIP       IMAGE
+PPIMAGE.DARK            INPUT    @DETDB                        CHIP       DARK
+PPIMAGE.SHUTTER         INPUT    @DETDB                        CHIP       IMAGE
+PPIMAGE.FLAT            INPUT    @DETDB                        CHIP       IMAGE
+PPIMAGE.MASK            INPUT    @DETDB                        CHIP       MASK
+
+## files used by psphot
+PSPHOT.LOAD             INPUT    @FILES                        CHIP       IMAGE
+PSPHOT.INPUT            INPUT    @FILES                        CHIP       IMAGE
+PSPHOT.MASK             INPUT    @FILES                        CHIP       MASK     
+PSPHOT.WEIGHT           INPUT    @FILES                        CHIP       WEIGHT     
+PSPHOT.PSF.LOAD         INPUT    @FILES                        CHIP       PSF       
+PSPHOT.INPUT.CMF        INPUT    @FILES                        CHIP       CMF       
+
+PSWARP.INPUT            INPUT    none.fits                     CHIP       IMAGE
+PSWARP.WEIGHT           INPUT    none.fits                     CHIP       WEIGHT
+PSWARP.MASK             INPUT    none.fits                     CHIP       MASK
+PSWARP.SKYCELL          INPUT    none.fits                     FPA        IMAGE
+PSWARP.ASTROM           INPUT    none.fits                     CHIP       CMF
+
+PSASTRO.WCS             INPUT    none.fits                     CHIP       CMF
+PSASTRO.MODEL           INPUT    @DETDB                        FPA        ASTROM
+
+PSASTRO.INPUT.CMP       INPUT    none.fits                     CHIP       CMP
+PSASTRO.INPUT.CMF       INPUT    none.fits                     CHIP       CMF
+
+PPSUB.INPUT             INPUT    none.fits                     FPA        IMAGE
+PPSUB.INPUT.MASK        INPUT    none.fits                     FPA        MASK
+PPSUB.INPUT.WEIGHT      INPUT    none.fits                     FPA        WEIGHT
+PPSUB.REF               INPUT    none.fits                     FPA        IMAGE
+PPSUB.REF.MASK          INPUT    none.fits                     FPA        MASK
+PPSUB.REF.WEIGHT        INPUT    none.fits                     FPA        WEIGHT
+
+PPSTACK.INPUT           INPUT    @FILES                        FPA        IMAGE
+PPSTACK.INPUT.MASK      INPUT    @FILES                        FPA        MASK
+PPSTACK.INPUT.WEIGHT    INPUT    @FILES                        FPA        WEIGHT
+PPSTACK.SOURCES         INPUT    @FILES                        FPA        CMF
+
+PPSTAMP.INPUT           INPUT    none.fits                     CHIP       IMAGE
+
+PPARITH.INPUT.IMAGE     INPUT    none.fits                     CHIP       IMAGE
+PPARITH.INPUT.MASK      INPUT    none.fits                     CHIP       MASK
+
+### output file definitions
+TYPE                   OUTPUT FILENAME.RULE                     FILE.TYPE FITS.TYPE DATA.LEVEL FILE.SAVE FILE.FORMAT
+PPIMAGE.OUTPUT         OUTPUT {OUTPUT}.{CHIP.NAME}.fits         IMAGE     COMP_POS   CHIP      TRUE      NONE
+PPIMAGE.OUTPUT.MASK    OUTPUT {OUTPUT}.{CHIP.NAME}.mk.fits      MASK      COMP_MASK  CHIP      TRUE      NONE
+PPIMAGE.OUTPUT.WEIGHT  OUTPUT {OUTPUT}.{CHIP.NAME}.wt.fits      WEIGHT    COMP_POS   CHIP      TRUE      NONE
+PPIMAGE.OUTPUT.DETMASK OUTPUT {OUTPUT}.{CHIP.NAME}.fits         IMAGE     COMP_MASK  CHIP      TRUE      NONE
+
+PPIMAGE.CHIP.RAW       OUTPUT {OUTPUT}.{CHIP.NAME}.ch.fits      IMAGE     COMP_POS   CHIP      TRUE      NONE
+# PPIMAGE.CHIP.MK      OUTPUT {OUTPUT}.{CHIP.NAME}.ch.fits      IMAGE     COMP_MASK  CHIP      TRUE      NONE
+PPIMAGE.CHIP           OUTPUT {OUTPUT}.{CHIP.NAME}.ch.fits      IMAGE     COMP_POS   CHIP      TRUE      NONE
+PPIMAGE.CHIP.MASK      OUTPUT {OUTPUT}.{CHIP.NAME}.ch.mk.fits   MASK      COMP_MASK  CHIP      TRUE      NONE
+PPIMAGE.CHIP.WEIGHT    OUTPUT {OUTPUT}.{CHIP.NAME}.ch.wt.fits   WEIGHT    COMP_POS   CHIP      TRUE      NONE
+
+PPIMAGE.OUTPUT.FPA1    OUTPUT {OUTPUT}.b1.fits                  IMAGE     NONE       FPA       TRUE      NONE
+PPIMAGE.OUTPUT.FPA2    OUTPUT {OUTPUT}.b2.fits                  IMAGE     NONE       FPA       TRUE      NONE
+
+PPIMAGE.JPEG1          OUTPUT {OUTPUT}.b1.jpg                   JPEG      NONE       FPA       TRUE      NONE
+PPIMAGE.JPEG2          OUTPUT {OUTPUT}.b2.jpg                   JPEG      NONE       FPA       TRUE      NONE
+PPIMAGE.BIN1           OUTPUT {OUTPUT}.{CHIP.NAME}.b1.fits      IMAGE     NONE       CHIP      TRUE      NONE
+PPIMAGE.BIN2           OUTPUT {OUTPUT}.{CHIP.NAME}.b2.fits      IMAGE     NONE       CHIP      TRUE      NONE
+
+PPIMAGE.STATS          OUTPUT {OUTPUT}.{CHIP.NAME}.stats        STATS     NONE       CHIP      TRUE      NONE
+
+PPMERGE.OUTPUT         OUTPUT {OUTPUT}.{CHIP.NAME}.fits         IMAGE     NONE       CHIP      TRUE      NONE
+
+DVOCORR.OUTPUT         OUTPUT {OUTPUT}.{CHIP.NAME}.fc.fits      IMAGE     NONE       CHIP      TRUE      NONE
+DVOFLAT.OUTPUT         OUTPUT {OUTPUT}.{CHIP.NAME}.co.fits      IMAGE     NONE       CHIP      TRUE      NONE
+
+PSPHOT.RESID           OUTPUT {OUTPUT}.{CHIP.NAME}.res.fits     IMAGE     COMP_SUB   CHIP      TRUE      NONE
+PSPHOT.BACKGND         OUTPUT {OUTPUT}.{CHIP.NAME}.bck.fits     IMAGE     COMP_POS   CHIP      TRUE      NONE
+PSPHOT.BACKSUB         OUTPUT {OUTPUT}.{CHIP.NAME}.sub.fits     IMAGE     COMP_SUB   CHIP      TRUE      NONE
+PSPHOT.BACKMDL         OUTPUT {OUTPUT}.{CHIP.NAME}.mdl.fits     IMAGE     COMP_POS   CHIP      TRUE      NONE
+
+PSPHOT.OUTPUT.RAW      OUTPUT {OUTPUT}.{CHIP.NAME}              RAW       NONE       CHIP      TRUE      NONE
+PSPHOT.OUTPUT.SX       OUTPUT {OUTPUT}.{CHIP.NAME}.sx           SX        NONE       CHIP      TRUE      NONE
+PSPHOT.OUTPUT.OBJ      OUTPUT {OUTPUT}.{CHIP.NAME}.obj          OBJ       NONE       CHIP      TRUE      NONE
+PSPHOT.OUTPUT.CMP      OUTPUT {OUTPUT}.{CHIP.NAME}.cmp          CMP       NONE       CHIP      TRUE      NONE
+PSPHOT.OUT.CMF.SPL     OUTPUT {OUTPUT}.{CHIP.NAME}.cmf          CMF       NONE       CHIP      TRUE      NONE
+PSPHOT.OUT.CMF.MEF     OUTPUT {OUTPUT}.cmf                      CMF       NONE       FPA       TRUE      MEF
+
+PSPHOT.PSF.SAVE        OUTPUT {OUTPUT}.psf                      PSF       NONE       CHIP      TRUE      MEF
+
+SOURCE.PLOT.MOMENTS    OUTPUT {OUTPUT}.{CHIP.NAME}.mnt.png      KAPA      NONE       CHIP      TRUE      NONE
+SOURCE.PLOT.PSFMODEL   OUTPUT {OUTPUT}.{CHIP.NAME}.psf.png      KAPA      NONE       CHIP      TRUE      NONE
+SOURCE.PLOT.APRESID    OUTPUT {OUTPUT}.{CHIP.NAME}.dap.png      KAPA      NONE       CHIP      TRUE      NONE
+
+PSASTRO.OUTPUT.CMP     OUTPUT {OUTPUT}.{CHIP.NAME}.smp          CMP       NONE       CHIP      TRUE      NONE
+PSASTRO.OUT.CMF.SPL    OUTPUT {OUTPUT}.{CHIP.NAME}.smf          CMF       NONE       CHIP      TRUE      NONE
+PSASTRO.OUT.CMF.MEF    OUTPUT {OUTPUT}.smf                      CMF       NONE       FPA       TRUE      MEF
+PSASTRO.OUT.MODEL      OUTPUT {OUTPUT}.asm               ASTROM.MODEL     NONE       FPA       TRUE      NONE
+PSASTRO.OUT.REFSTARS   OUTPUT {OUTPUT}.aref.fits         ASTROM.REFSTARS  NONE       FPA       TRUE      NONE
+
+PSWARP.OUTPUT          OUTPUT {OUTPUT}.fits                     IMAGE     COMP_POS   FPA       TRUE      NONE
+PSWARP.OUTPUT.MASK     OUTPUT {OUTPUT}.mask.fits                MASK      COMP_MASK  FPA       TRUE      NONE
+PSWARP.OUTPUT.WEIGHT   OUTPUT {OUTPUT}.wt.fits                  WEIGHT    COMP_POS   FPA       TRUE      NONE
+PSWARP.OUTPUT.SOURCES  OUTPUT {OUTPUT}.cmf                      CMF       NONE       FPA       TRUE      NONE
+PSWARP.BIN1            OUTPUT {OUTPUT}.b1.fits                  IMAGE     NONE       FPA       TRUE      NONE
+PSWARP.BIN2            OUTPUT {OUTPUT}.b2.fits                  IMAGE     NONE       FPA       TRUE      NONE
+
+SKYCELL.STATS          OUTPUT {OUTPUT}.stats                    STATS     NONE       FPA       TRUE      NONE
+SKYCELL.TEMPLATE       OUTPUT {OUTPUT}.skycell                  SKYCELL   NONE       FPA       TRUE      NONE
+
+PPSUB.OUTPUT           OUTPUT {OUTPUT}.fits                     IMAGE     COMP_SUB   FPA       TRUE      NONE
+PPSUB.OUTPUT.MASK      OUTPUT {OUTPUT}.mask.fits                MASK      COMP_MASK  FPA       TRUE      NONE
+PPSUB.OUTPUT.WEIGHT    OUTPUT {OUTPUT}.wt.fits                  WEIGHT    COMP_POS   FPA       TRUE      NONE
+
+PPSTACK.OUTPUT         OUTPUT {OUTPUT}.fits                     IMAGE     COMP_POS   FPA       TRUE      NONE
+PPSTACK.OUTPUT.MASK    OUTPUT {OUTPUT}.mask.fits                MASK      COMP_MASK  FPA       TRUE      NONE
+PPSTACK.OUTPUT.WEIGHT  OUTPUT {OUTPUT}.wt.fits                  WEIGHT    COMP_POS   FPA       TRUE      NONE
+
+PPSTAMP.OUTPUT         OUTPUT {OUTPUT}.fits                     IMAGE     NONE       FPA       TRUE      NONE
+PPSTAMP.CHIP           OUTPUT {OUTPUT}.ch.fits                  IMAGE     NONE      CHIP       FALSE     MEF
+
+PPSIM.OUTPUT           OUTPUT {OUTPUT}.fits                     IMAGE     NONE       FPA       TRUE      MEF
+
+PPARITH.OUTPUT.IMAGE   OUTPUT {OUTPUT}.fits                     IMAGE     COMP_POS  CHIP       TRUE      NONE
+PPARITH.OUTPUT.MASK    OUTPUT {OUTPUT}.fits                     MASK      COMP_MASK CHIP       TRUE      NONE
+
+LOG.IMFILE             OUTPUT {OUTPUT}.{CHIP.NAME}.log          TEXT      NONE       CHIP      TRUE      NONE
+LOG.EXP                OUTPUT {OUTPUT}.log                      TEXT      NONE       FPA       TRUE      NONE
+
+TRACE.IMFILE           OUTPUT {OUTPUT}.{CHIP.NAME}.trace        TEXT      NONE       CHIP      TRUE      NONE
+TRACE.EXP              OUTPUT {OUTPUT}.trace                    TEXT      NONE       FPA       TRUE      NONE
Index: /branches/pap_branch_080320/ippconfig/megacam/ppMerge.config
===================================================================
--- /branches/pap_branch_080320/ippconfig/megacam/ppMerge.config	(revision 17139)
+++ /branches/pap_branch_080320/ippconfig/megacam/ppMerge.config	(revision 17139)
@@ -0,0 +1,40 @@
+
+# Bias combination --- don't want min/max rejection
+PPMERGE_BIAS	METADATA
+   REJ		F32	3.0		# Rejection threshold (sigma)
+   ITER		S32	2		# Number of rejection iterations
+   FRACHIGH	F32	0.0		# Fraction of high pixels to reject immediately
+   FRACLOW	F32	0.0		# Fraction of low pixels to reject immediately
+   WEIGHTS	BOOL	FALSE		# Use image weights?
+   COMBINE	STR	CLIPPED		# Statistic to use for combination: 
+END
+
+# Dark combination --- don't want min/max rejection
+# More aggressive clipping than bias, so as to remove CRs
+PPMERGE_DARK	METADATA
+  REJ		F32	2.0		# Rejection threshold (sigma)
+  ITER		S32	4		# Number of rejection iterations
+  FRACHIGH	F32	0.0		# Fraction of high pixels to reject immediately
+  FRACLOW	F32	0.0		# Fraction of low pixels to reject immediately
+  WEIGHTS	BOOL	FALSE		# Use image weights?
+  COMBINE	STR	CLIPPED		# Statistic to use for combination: 
+END
+
+# Flat combination --- use min/max rejection
+PPMERGE_FLAT	METADATA
+	REJ		F32	3.0		# Rejection threshold (sigma)
+	ITER		S32	1		# Number of rejection iterations
+	FRACHIGH	F32	0.3		# Fraction of high pixels to reject immediately
+	FRACLOW		F32	0.1		# Fraction of low pixels to reject immediately
+	NKEEP		S32	5		# Minimum number of pixels in stack to keep
+	WEIGHTS		BOOL	TRUE		# Use image weights?
+	COMBINE		STR	MEAN		# Statistic to use for combination: 
+END
+
+
+# Fringe combination --- already included in default, above
+PPMERGE_FRINGE	METADATA
+	FRACHIGH	F32	0.1		# Fraction of high pixels to reject immediately
+	WEIGHTS		BOOL	TRUE		# Use image weights?
+END
+
Index: /branches/pap_branch_080320/ippconfig/simmosaic/camera.config
===================================================================
--- /branches/pap_branch_080320/ippconfig/simmosaic/camera.config	(revision 17139)
+++ /branches/pap_branch_080320/ippconfig/simmosaic/camera.config	(revision 17139)
@@ -0,0 +1,376 @@
+# Camera configuration file for simulated mosaic
+
+# File formats that we know about
+FORMATS         METADATA
+	TOGETHER STR	simmosaic/format_together.config
+	SPLIT	STR	simmosaic/format_split.config
+END
+ 
+# Description of camera --- all the chips and the cells that comprise them
+FPA     METADATA
+        Chip00		STR	Cell00 Cell01 Cell10 Cell11
+        Chip01		STR	Cell00 Cell01 Cell10 Cell11
+        Chip10		STR	Cell00 Cell01 Cell10 Cell11
+        Chip11		STR	Cell00 Cell01 Cell10 Cell11
+END
+
+# move this elsewhere?  we need a lookup table to go from filter ID to abstract name
+FILTER.ID       METADATA
+	NONE	STR	NONE
+	B	STR	B
+	V	STR	V
+	R	STR	R
+	I	STR	I
+	g	STR	g
+	r	STR	r
+	i	STR	i
+	z	STR	z
+END
+
+# Table of strings to use for the class identifier (see ippdb) according to the file level.
+CLASSID		METADATA
+	FPA	STR	fpa
+	CHIP	STR	{CHIP.NAME}
+	CELL	STR	{CHIP.NAME}.{CELL.NAME}
+	fpa	STR	fpa
+	chip	STR	{CHIP.NAME}
+	cell	STR	{CHIP.NAME}.{CELL.NAME}
+END
+
+DVO.CAMERADIR	STR	simmosaic		# Camera directory for DVO
+
+# convert supplied FPA.OBSTYPE values to abstract exptype names
+OBSTYPE.TABLE METADATA
+  bias 	   STR BIAS
+  zero 	   STR BIAS
+  dark 	   STR DARK
+  flat 	   STR SKYFLAT
+  skyflat  STR SKYFLAT
+  domeflat STR DOMEFLAT
+  object   STR OBJECT
+  science  STR OBJECT
+END
+
+# Recipe options
+RECIPES		METADATA
+        PPIMAGE         STR     simmosaic/ppImage.config	# Default: all (normal) options on
+	PSPHOT		STR	simmosaic/psphot.config		# psphot details
+	PSASTRO		STR	simmosaic/psastro.config	# psastro details
+	REJECTIONS	STR	simmosaic/rejections.config	# Rejection limits
+END
+
+# Reduction classes
+REDUCTION	METADATA
+	# Detrend processing
+	DETREND		METADATA
+		BIAS_PROCESS	STR	PPIMAGE_O
+		BIAS_RESID	STR	PPIMAGE_B
+		BIAS_VERIFY	STR	PPIMAGE_OB
+		BIAS_STACK	STR	PPMERGE_BIAS
+		DARK_PROCESS	STR	PPIMAGE_OB
+		DARK_RESID	STR	PPIMAGE_D
+		DARK_VERIFY	STR	PPIMAGE_OBD
+		DARK_STACK	STR	PPMERGE_DARK
+		SHUTTER_PROCESS	STR	PPIMAGE_OBD
+		SHUTTER_RESID	STR	PPIMAGE_S
+		SHUTTER_VERIFY	STR	PPIMAGE_OBDS
+		SHUTTER_STACK	STR	PPMERGE_SHUTTER
+		FLAT_PROCESS	STR	PPIMAGE_OBDS
+		FLAT_RESID	STR	PPIMAGE_F
+		FLAT_VERIFY	STR	PPIMAGE_OBDSF
+		FLAT_STACK	STR	PPMERGE_FLAT
+		FRINGE_PROCESS	STR	PPIMAGE_OBDSF
+		FRINGE_RESID	STR	PPIMAGE_R
+		FRINGE_VERIFY	STR	PPIMAGE_OBDSFR
+		FRINGE_STACK	STR	PPMERGE_FRINGE
+
+		# Generation of pixel masks from darks and flats
+		DARKMASK_PROCESS	STR	PPIMAGE_OBD
+		DARKMASK_RESID		STR	PPIMAGE_N
+		DARKMASK_VERIFY		STR	PPIMAGE_OBD
+		DARKMASK_STACK		STR	PPMERGE_DARKMASK
+		FLATMASK_PROCESS	STR	PPIMAGE_OBDSF
+		FLATMASK_RESID		STR	PPIMAGE_N
+		FLATMASK_VERIFY		STR	PPIMAGE_OBDSF
+		FLATMASK_STACK		STR	PPMERGE_FLATMASK
+		JPEG_BIN1_IMAGE_DARKMASK STR	PPIMAGE_J1_IMAGE_M
+		JPEG_BIN2_IMAGE_DARKMASK STR	PPIMAGE_J2_IMAGE_M
+		JPEG_BIN1_IMAGE_FLATMASK STR	PPIMAGE_J1_IMAGE_M
+		JPEG_BIN2_IMAGE_FLATMASK STR	PPIMAGE_J2_IMAGE_M
+		JPEG_BIN1_RESID_DARKMASK STR	PPIMAGE_J1_RESID_M
+		JPEG_BIN2_RESID_DARKMASK STR	PPIMAGE_J2_RESID_M
+		JPEG_BIN1_RESID_FLATMASK STR	PPIMAGE_J1_RESID_M
+		JPEG_BIN2_RESID_FLATMASK STR	PPIMAGE_J2_RESID_M
+
+ 		JPEG_BIN1_IMAGE_BIAS     STR  PPIMAGE_J1_IMAGE_B
+		JPEG_BIN1_IMAGE_DARK     STR  PPIMAGE_J1_IMAGE_B
+		JPEG_BIN1_IMAGE_SHUTTER  STR  PPIMAGE_J1_IMAGE_F
+		JPEG_BIN1_IMAGE_FLAT     STR  PPIMAGE_J1_IMAGE_F
+		JPEG_BIN1_IMAGE_DOMEFLAT STR  PPIMAGE_J1_IMAGE_F
+		JPEG_BIN1_IMAGE_SKYFLAT  STR  PPIMAGE_J1_IMAGE_F
+		JPEG_BIN1_IMAGE_FRINGE   STR  PPIMAGE_J1_IMAGE_R
+ 		JPEG_BIN2_IMAGE_BIAS     STR  PPIMAGE_J2_IMAGE_B
+		JPEG_BIN2_IMAGE_DARK     STR  PPIMAGE_J2_IMAGE_B
+		JPEG_BIN2_IMAGE_SHUTTER  STR  PPIMAGE_J2_IMAGE_F
+		JPEG_BIN2_IMAGE_FLAT     STR  PPIMAGE_J2_IMAGE_F
+		JPEG_BIN2_IMAGE_DOMEFLAT STR  PPIMAGE_J2_IMAGE_F
+		JPEG_BIN2_IMAGE_SKYFLAT  STR  PPIMAGE_J2_IMAGE_F
+		JPEG_BIN2_IMAGE_FRINGE   STR  PPIMAGE_J2_IMAGE_R
+
+ 		JPEG_BIN1_RESID_BIAS     STR  PPIMAGE_J1_RESID_B
+		JPEG_BIN1_RESID_DARK     STR  PPIMAGE_J1_RESID_B
+		JPEG_BIN1_RESID_SHUTTER  STR  PPIMAGE_J1_RESID_F
+		JPEG_BIN1_RESID_FLAT     STR  PPIMAGE_J1_RESID_F
+		JPEG_BIN1_RESID_DOMEFLAT STR  PPIMAGE_J1_RESID_F
+		JPEG_BIN1_RESID_SKYFLAT  STR  PPIMAGE_J1_RESID_F
+		JPEG_BIN1_RESID_FRINGE   STR  PPIMAGE_J1_RESID_R
+ 		JPEG_BIN2_RESID_BIAS     STR  PPIMAGE_J2_RESID_B
+		JPEG_BIN2_RESID_DARK     STR  PPIMAGE_J2_RESID_B
+		JPEG_BIN2_RESID_SHUTTER  STR  PPIMAGE_J2_RESID_F
+		JPEG_BIN2_RESID_FLAT     STR  PPIMAGE_J2_RESID_F
+		JPEG_BIN2_RESID_DOMEFLAT STR  PPIMAGE_J2_RESID_F
+		JPEG_BIN2_RESID_SKYFLAT  STR  PPIMAGE_J2_RESID_F
+		JPEG_BIN2_RESID_FRINGE   STR  PPIMAGE_J2_RESID_R
+	END
+	# Processing raw data
+	DEFAULT		METADATA
+		CHIP		STR	PPIMAGE_OBDSFRA
+ 		JPEG_BIN1       STR     PPIMAGE_J1
+ 		JPEG_BIN2       STR     PPIMAGE_J2
+	END
+END
+
+FITS    METADATA
+# BITPIX is the bits per pixel for writing the output data
+# COMP = NONE|RICE|GZIP|HCOMPRESS|PLIO is the compression algorithm
+# TILE.[XYZ] are the tile sizes.  0 means entire the dimension, so (0,1,1) forms tiles from rows
+# NOISE [0..16] is the number of "noise bits" to preserve when quantising floating point data; 16 for no loss
+# HSCALE is the scale factor for lossy compression with HCOMPRESS; 0 or 1 for none; 2*RMS --> 10x compression
+# HSMOOTH is the smoothing to apply to HCOMPRESSed data when decompressing; 0 for none
+
+# BITPIX(S32) is the bits per pixel for writing the output data
+# COMPRESSION(STR) = NONE|RICE|GZIP|HCOMPRESS|PLIO is the compression algorithm
+# TILE.[XYZ](S32) are the tile sizes.  0 means entire the dimension, so (0,1,1) forms tiles from rows
+# NOISE(S32) [0..16] is the number of "noise bits" to preserve when quantising floating point data
+# HSCALE(S32) is the scale factor for lossy compression with HCOMPRESS; 0 or 1 for none; 2*RMS --> 10x compression
+# HSMOOTH(S32) is the smoothing to apply to HCOMPRESSed data when decompressing; 0 for none
+# SCALING(STR) = NONE|RANGE|STDEV_POSITIVE|STDEV_NEGATIVE|STDEV_BOTH|MANUAL is the scaling scheme
+# BSCALE(F32) is the manual scaling to apply (when SCALING = MANUAL)
+# BZERO(F32) is the manual zero-point to apply (when SCALING = MANUAL)
+# STDEV.BITS(S32) is the number of bits to map to a standard deviation (when SCALING = STDEV_*)
+# STDEV.NUM(F32) is the number of standard deviations to the edge (when SCALING = STDEV_NEGATIVE|STDEV_POSITIVE)
+# FLOAT(STR) is the name of a custom floating-point type
+
+	DET_IMAGE	METADATA
+		BITPIX		S32	-32
+	END
+	DET_MASK	METADATA
+		BITPIX		S32	8
+	END
+	DET_WEIGHT	METADATA
+		BITPIX		S32	-32
+	END
+
+	SKY_IMAGE	METADATA
+		BITPIX		S32	-32
+	END
+	SKY_MASK	METADATA
+		BITPIX		S32	8
+	END
+	SKY_WEIGHT	METADATA
+		BITPIX		S32	-32
+	END
+
+	COMPRESSED_POSITIVE	METADATA
+		BITPIX		S32	16
+		SCALING		STR	STDEV_POSITIVE
+		STDEV.BITS	S32	4
+		STDEV.NUM	F32	10
+		COMPRESSION	STR	RICE
+		TILE.X		S32	0
+		TILE.Y		S32	1
+		TILE.Z		S32	1
+		NOISE		S32	8
+	END
+	COMPRESSED_MASK		METADATA
+		COMPRESSION	STR	PLIO
+		TILE.X		S32	0
+		TILE.Y		S32	1
+		TILE.Z		S32	1
+		NOISE		S32	8
+	END
+	COMPRESSED_SUBTRACTION	METADATA
+		BITPIX		S32	16
+		SCALING		STR	STDEV_BOTH
+		STDEV.BITS	S32	4
+		STDEV.NUM	F32	5
+		COMPRESSION	STR	RICE
+		TILE.X		S32	0
+		TILE.Y		S32	1
+		TILE.Z		S32	1
+		NOISE		S32	8
+	END
+
+END
+
+FILERULES METADATA
+   ### Redirections
+   PSASTRO.INPUT      STR PSASTRO.INPUT.CMF
+   PSASTRO.OUTPUT     STR PSASTRO.OUT.CMF.SPL
+   PSASTRO.OUTPUT.MEF STR PSASTRO.OUT.CMF.MEF
+   PSPHOT.OUTPUT      STR PSPHOT.OUT.CMF.SPL
+
+   ### input file definitions
+   ### use @DETDB entries to get the detrend images from the database
+   ### replace @DETDB with @FILES if you want to require it from the 
+   ### command line, or with an explicit name to require a specific file
+   TYPE               INPUT FILENAME.RULE DATA.LEVEL FILE.TYPE 
+
+   ## files used by ppImage
+   PPIMAGE.INPUT      INPUT @FILES        CHIP       IMAGE     
+   PPIMAGE.MASK       INPUT @DETDB        CELL       IMAGE     
+   PPIMAGE.BIAS       INPUT @DETDB        CELL       IMAGE     
+   PPIMAGE.DARK       INPUT @DETDB        CELL       DARK     
+   PPIMAGE.FLAT       INPUT @DETDB        CELL       IMAGE     
+   PPIMAGE.FRINGE     INPUT @DETDB        CHIP       FRINGE     
+   PPIMAGE.SHUTTER    INPUT @DETDB        CELL       IMAGE     
+
+   ## files used to build and apply the flat-field correction images
+   DVOCORR.INPUT      INPUT @FILES        CHIP       IMAGE     
+   DVOCORR.REFHEAD    INPUT @FILES        CHIP       HEADER     
+   DVOFLAT.INPUT      INPUT @FILES        CHIP       IMAGE 	
+   DVOFLAT.CORR       INPUT @DETDB        CHIP       IMAGE 	
+
+   ## files used by psphot 
+   PSPHOT.LOAD        INPUT @FILES        CHIP       IMAGE     
+   PSPHOT.INPUT       INPUT @FILES        CHIP       IMAGE     
+   PSPHOT.MASK        INPUT @FILES        CHIP       MASK     
+   PSPHOT.WEIGHT      INPUT @FILES        CHIP       WEIGHT     
+   PSPHOT.PSF.LOAD    INPUT @FILES        CHIP	     PSF       
+   PSPHOT.INPUT.CMF   INPUT @FILES        CHIP	     CMF       
+
+   ## files used by psastro 
+   PSASTRO.INPUT.CMF  INPUT @FILES        CHIP       CMF       
+
+   ## files used by pswarp
+   PSWARP.INPUT       INPUT @FILES        CHIP       IMAGE     
+   PSWARP.SKYCELL     INPUT @FILES        FPA        IMAGE     
+   PSWARP.WEIGHT      INPUT @FILES        CHIP       WEIGHT
+   PSWARP.MASK        INPUT @FILES        CHIP       MASK
+   PSWARP.ASTROM      INPUT @FILES        CHIP       CMF       
+
+   PPSUB.INPUT        INPUT    none.fits  FPA	     IMAGE
+   PPSUB.INPUT.MASK   INPUT    none.fits  FPA	     MASK
+   PPSUB.INPUT.WEIGHT INPUT    none.fits  FPA	     WEIGHT
+   PPSUB.REF          INPUT    none.fits  FPA	     IMAGE
+   PPSUB.REF.MASK     INPUT    none.fits  FPA	     MASK
+   PPSUB.REF.WEIGHT   INPUT    none.fits  FPA	     WEIGHT
+
+   PPSTACK.INPUT      INPUT    none.fits  FPA	     IMAGE
+   PPSTACK.INPUT.MASK INPUT    none.fits  FPA	     MASK
+   PPSTACK.INPUT.WEIGHT INPUT  none.fits  FPA	     WEIGHT
+
+   PPSTAMP.INPUT      INPUT @FILES        CHIP       IMAGE
+
+   PPARITH.INPUT.IMAGE INPUT @FILES        CHIP       IMAGE
+   PPARITH.INPUT.MASK  INPUT @FILES        CHIP       MASK
+
+
+   ### output file definitions
+   TYPE                OUTPUT   FILENAME.RULE      	       FILE.TYPE FITS.TYPE DATA.LEVEL FILE.SAVE FILE.FORMAT
+   PPIMAGE.OUTPUT      OUTPUT   {OUTPUT}.{CHIP.NAME}.fits      IMAGE     NONE      CHIP       TRUE      SPLIT
+   PPIMAGE.OUTPUT.MASK OUTPUT   {OUTPUT}.{CHIP.NAME}.mask.fits MASK      NONE      CHIP       TRUE      SPLIT
+   PPIMAGE.OUTPUT.WEIGHT OUTPUT {OUTPUT}.{CHIP.NAME}.wt.fits   WEIGHT    NONE      CHIP       TRUE      SPLIT
+   PPIMAGE.CHIP        OUTPUT 	{OUTPUT}.{CHIP.NAME}.ch.fits   IMAGE     NONE      CHIP       TRUE      SPLIT
+   PPIMAGE.CHIP.MASK   OUTPUT 	{OUTPUT}.{CHIP.NAME}.ch.mask.fits MASK   NONE      CHIP       TRUE      SPLIT
+   PPIMAGE.CHIP.WEIGHT OUTPUT 	{OUTPUT}.{CHIP.NAME}.ch.wt.fits WEIGHT   NONE      CHIP       TRUE      SPLIT
+   PPIMAGE.OUTPUT.FPA1 OUTPUT 	{OUTPUT}.fpa1.fits             IMAGE     NONE      FPA        TRUE      NONE
+   PPIMAGE.OUTPUT.FPA2 OUTPUT 	{OUTPUT}.fpa2.fits             IMAGE     NONE      FPA        TRUE      NONE
+   PPIMAGE.STATS       OUTPUT   {OUTPUT}.stats                 STATS     NONE      FPA        TRUE      NONE
+
+   PPIMAGE.JPEG1       OUTPUT   {OUTPUT}.b1.jpg    	       JPEG      NONE      FPA        TRUE      NONE
+   PPIMAGE.JPEG2       OUTPUT   {OUTPUT}.b2.jpg    	       JPEG      NONE      FPA        TRUE      NONE
+   PPIMAGE.BIN1        OUTPUT   {OUTPUT}.{CHIP.NAME}.b1.fits   IMAGE     NONE      CHIP       TRUE      SPLIT
+   PPIMAGE.BIN2        OUTPUT   {OUTPUT}.{CHIP.NAME}.b2.fits   IMAGE     NONE      CHIP       TRUE      SPLIT
+
+   PPMERGE.OUTPUT      OUTPUT   {OUTPUT}.{CHIP.NAME}.fits      IMAGE     NONE      CHIP       TRUE      NONE
+
+   DVOCORR.OUTPUT      OUTPUT 	{OUTPUT}.{CHIP.NAME}.fc.fits   IMAGE     NONE      CHIP       TRUE      NONE
+   DVOFLAT.OUTPUT      OUTPUT 	{OUTPUT}.{CHIP.NAME}.co.fits   IMAGE     NONE      CHIP       TRUE      NONE
+
+   PSPHOT.RESID        OUTPUT   {OUTPUT}.{CHIP.NAME}.res.fits  IMAGE     NONE      CHIP       TRUE      NONE
+   PSPHOT.BACKGND      OUTPUT   {OUTPUT}.{CHIP.NAME}.bck.fits  IMAGE     NONE      CHIP       TRUE      NONE
+   PSPHOT.BACKSUB      OUTPUT   {OUTPUT}.{CHIP.NAME}.sub.fits  IMAGE     NONE      CHIP       TRUE      NONE
+   PSPHOT.BACKMDL      OUTPUT   {OUTPUT}.{CHIP.NAME}.mdl.fits  IMAGE     NONE      CHIP       TRUE      NONE
+   PSPHOT.BACKMDL.STDEV OUTPUT   {OUTPUT}.{CHIP.NAME}.mdd.fits IMAGE     NONE      CHIP       TRUE      NONE
+
+   PSPHOT.OUTPUT.RAW   OUTPUT   {OUTPUT}.{CHIP.NAME}           RAW       NONE      CHIP       TRUE      NONE
+   PSPHOT.OUTPUT.SX    OUTPUT   {OUTPUT}.{CHIP.NAME}.sx        SX        NONE      CHIP       TRUE      NONE
+   PSPHOT.OUTPUT.OBJ   OUTPUT   {OUTPUT}.{CHIP.NAME}.obj       OBJ       NONE      CHIP       TRUE      NONE
+   PSPHOT.OUT.CMF.SPL  OUTPUT   {OUTPUT}.{CHIP.NAME}.cmf       CMF       NONE      CHIP       TRUE      NONE
+   PSPHOT.OUT.CMF.MEF  OUTPUT   {OUTPUT}.cmf                   CMF       NONE      FPA        TRUE      NONE
+
+   PSPHOT.PSF.SAVE     OUTPUT   {OUTPUT}.{CHIP.NAME}.psf       PSF       NONE      CHIP       TRUE      NONE
+
+   SOURCE.PLOT.MOMENTS  OUTPUT  {OUTPUT}.{CHIP.NAME}.mnt.png   KAPA      NONE      CHIP       TRUE      NONE
+   SOURCE.PLOT.PSFMODEL OUTPUT  {OUTPUT}.{CHIP.NAME}.psf.png   KAPA      NONE      CHIP       TRUE      NONE
+   SOURCE.PLOT.APRESID  OUTPUT  {OUTPUT}.{CHIP.NAME}.dap.png   KAPA      NONE      CHIP       TRUE      NONE
+
+   PSASTRO.OUTPUT.CMP   OUTPUT   {OUTPUT}.{CHIP.NAME}.smp      CMP       NONE      CHIP       TRUE      SPLIT
+   PSASTRO.OUT.CMF.SPL  OUTPUT   {OUTPUT}.{CHIP.NAME}.smf      CMF       NONE      CHIP       TRUE      SPLIT
+   PSASTRO.OUT.CMF.MEF  OUTPUT   {OUTPUT}.smf		       CMF       NONE      FPA        TRUE      TOGETHER
+
+   PSWARP.OUTPUT       OUTPUT   {OUTPUT}.fits      	       IMAGE     NONE      FPA        TRUE      NONE
+   PSWARP.OUTPUT.MASK  OUTPUT   {OUTPUT}.mask.fits             MASK      NONE      FPA        TRUE      NONE
+   PSWARP.OUTPUT.WEIGHT OUTPUT  {OUTPUT}.wt.fits               WEIGHT    NONE      FPA        TRUE      NONE
+   PSWARP.BIN1         OUTPUT   {OUTPUT}.b1.fits   	       IMAGE     NONE      FPA        TRUE      NONE
+   PSWARP.BIN2         OUTPUT   {OUTPUT}.b2.fits   	       IMAGE     NONE      FPA        TRUE      NONE
+
+   SKYCELL.STATS       OUTPUT   {OUTPUT}.stats                 STATS     NONE      FPA        TRUE      NONE
+   SKYCELL.TEMPLATE    OUTPUT   {OUTPUT}.skycell               SKYCELL   NONE      FPA        TRUE      NONE
+
+   PPSUB.OUTPUT        OUTPUT   {OUTPUT}.fits                  IMAGE     NONE      FPA        TRUE      NONE
+   PPSUB.OUTPUT.MASK   OUTPUT   {OUTPUT}.mask.fits             MASK      NONE      FPA        TRUE      NONE
+   PPSUB.OUTPUT.WEIGHT OUTPUT   {OUTPUT}.wt.fits               WEIGHT    NONE      FPA        TRUE      NONE
+
+   PPSTACK.OUTPUT      OUTPUT   {OUTPUT}.fits                  IMAGE     NONE      FPA        TRUE      NONE
+   PPSTACK.OUTPUT.MASK OUTPUT   {OUTPUT}.mask.fits             MASK      NONE      FPA        TRUE      NONE
+   PPSTACK.OUTPUT.WEIGHT OUTPUT {OUTPUT}.weight.fits           WEIGHT    NONE      FPA        TRUE      NONE
+
+   PPSTAMP.OUTPUT        OUTPUT {OUTPUT}.fits                  IMAGE     NONE      FPA        TRUE      NONE
+   PPSTAMP.CHIP.MEF      OUTPUT {OUTPUT}.ch.fits               IMAGE     NONE      CHIP       FALSE     MEF
+
+   PPSIM.OUTPUT        OUTPUT   {OUTPUT}.{CHIP.NAME}.fits      IMAGE     NONE      CHIP       TRUE      SPLIT
+   PPSIM.SOURCES       OUTPUT   {OUTPUT}.cmf                   CMF       NONE      FPA        TRUE      NONE
+
+   PPARITH.OUTPUT.IMAGE  OUTPUT {OUTPUT}.fits                  IMAGE     NONE      CHIP       TRUE      NONE
+   PPARITH.OUTPUT.MASK   OUTPUT {OUTPUT}.fits                  MASK     NONE      CHIP       TRUE      NONE
+
+   LOG.IMFILE            OUTPUT {OUTPUT}.{CHIP.NAME}.log       TEXT      NONE      CHIP       TRUE      NONE
+   LOG.EXP               OUTPUT {OUTPUT}.log                   TEXT      NONE      FPA        TRUE      NONE
+END
+
+# FPA file defines properties of a possible input|output object
+# user can set the filename (I|O), filename rules (O), or abstract source (@FILES, @DETDB) (I) 
+# user can set the extension name, if used
+# user can set the file type (IMAGE, JPEG, RAW, SX, OBJ, CMP, CMF) : but these are not variable in most cases!
+# user can set the file depth: only valid for output files
+# user can set the data depth: must be >= file depth
+# user can set the file format: only valid for newly created FPAs
+# user can set the colormap, scaling method, scaling range (JPEG only)
+# user can set the extension name for the data and header segments (CMF only)
+
+
+EXTNAME.RULES	METADATA
+	CMF.HEAD	STR	{CHIP.NAME}.hdr
+	CMF.DATA	STR	{CHIP.NAME}.psf
+	PSF.HEAD	STR	{CHIP.NAME}.hdr
+	PSF.TABLE	STR	{CHIP.NAME}.psf_model
+	PSF.RESID	STR	{CHIP.NAME}.psf_resid
+END
+
+BLANK.HEADERS	METADATA
+	FPA.TIME	STR	MJD-OBS
+	FPA.EXPOSURE	STR	EXPTIME
+	FPA.AIRMASS	STR	AIRMASS
+END
Index: /branches/pap_branch_080320/psModules/src/imcombine/pmReadoutCombine.c
===================================================================
--- /branches/pap_branch_080320/psModules/src/imcombine/pmReadoutCombine.c	(revision 17139)
+++ /branches/pap_branch_080320/psModules/src/imcombine/pmReadoutCombine.c	(revision 17139)
@@ -0,0 +1,356 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+#include <pslib.h>
+
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmHDUUtils.h"
+#include "pmFPAMaskWeight.h"
+#include "pmConceptsAverage.h"
+#include "pmReadoutStack.h"
+
+#include "pmReadoutCombine.h"
+
+//#define SHOW_BUSY 1                   // Show that the function is busy
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Public functions
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+// Allocator for pmCombineParams
+pmCombineParams *pmCombineParamsAlloc(psStatsOptions combine)
+{
+    pmCombineParams *params = psAlloc(sizeof(pmCombineParams));
+
+    params->combine = combine;
+    params->maskVal = 0;
+    params->blank = 0;
+    params->nKeep = 0;
+    params->fracHigh = 0.0;
+    params->fracHigh = 0.0;
+    params->iter = 1;
+    params->rej = INFINITY;
+    params->weights = false;
+
+    return params;
+}
+
+
+// XXX: Maybe add support for S16 and S32 types.  Currently, only F32 supported.
+bool pmReadoutCombine(pmReadout *output, const psArray *inputs, const psVector *zero, const psVector *scale,
+                      const pmCombineParams *params)
+{
+    // Check inputs
+    PS_ASSERT_PTR_NON_NULL(output, false);
+    PS_ASSERT_ARRAY_NON_NULL(inputs, false);
+    PS_ASSERT_PTR_NON_NULL(params, false);
+    if (zero) {
+        PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, false);
+        PS_ASSERT_VECTOR_SIZE(zero, inputs->n, false);
+    }
+    if (scale) {
+        PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, false);
+        PS_ASSERT_VECTOR_SIZE(scale, inputs->n, false);
+    }
+    PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracLow, 0.0, 1.0, false);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracHigh, 0.0, 1.0, false);
+    if (params->combine != PS_STAT_SAMPLE_MEAN && params->combine != PS_STAT_SAMPLE_MEDIAN &&
+            params->combine != PS_STAT_ROBUST_MEDIAN && params->combine != PS_STAT_FITTED_MEAN &&
+            params->combine != PS_STAT_CLIPPED_MEAN) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Combination method is not SAMPLE_MEAN, SAMPLE_MEDIAN, "
+                "ROBUST_MEDIAN, FITTED_MEAN or CLIPPED_MEAN.\n");
+        return false;
+    }
+    for (int i = 0; i < inputs->n; i++) {
+        pmReadout *readout = inputs->data[i]; // Readout of interest
+        if (params->weights && !readout->weight) {
+            psError(PS_ERR_UNEXPECTED_NULL, true,
+                    "Rejection based on weights requested, but no weights supplied for image %d.\n", i);
+            return false;
+        }
+    }
+
+    bool first = !output->image;        // First pass through?
+
+    pmHDU *hdu = pmHDUFromReadout(output); // Output HDU
+    if (!hdu) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find HDU for readout.\n");
+        return false;
+    }
+
+    if (first) {
+        psString comment = NULL;        // Comment to add to header
+        psStringAppend(&comment, "Combining using statistic: %x", params->combine);
+        if (!hdu->header) {
+            hdu->header = psMetadataAlloc();
+        }
+        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
+        psFree(comment);
+    }
+
+    psStats *stats = psStatsAlloc(params->combine); // The statistics to use in the combination
+    if (params->combine == PS_STAT_CLIPPED_MEAN) {
+        stats->clipSigma = params->rej;
+        stats->clipIter = params->iter;
+
+        if (first) {
+            psString comment = NULL;    // Comment to add to header
+            psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma",
+                           params->iter, params->rej);
+            psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
+            psFree(comment);
+        }
+    }
+
+    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
+    int xSize, ySize;                   // Size of the output image
+    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
+                                inputs)) {
+        psError(PS_ERR_UNKNOWN, false, "No valid input readouts.");
+        return false;
+    }
+
+    pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true, params->weights,
+                        params->blank);
+    psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0);
+
+    psStatsOptions combineStdev = 0; // Statistics option for weights
+    if (params->weights) {
+        if (first) {
+            psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
+                             "Using input weights to combine images", "");
+        }
+
+        // Get the correct statistics option for weights
+        switch (params->combine) {
+        case PS_STAT_SAMPLE_MEAN:
+        case PS_STAT_SAMPLE_MEDIAN:
+            combineStdev = PS_STAT_SAMPLE_STDEV;
+            break;
+        case PS_STAT_ROBUST_MEDIAN:
+            combineStdev = PS_STAT_ROBUST_STDEV;
+            break;
+        case PS_STAT_FITTED_MEAN:
+            combineStdev = PS_STAT_FITTED_STDEV;
+            break;
+        case PS_STAT_CLIPPED_MEAN:
+            combineStdev = PS_STAT_CLIPPED_STDEV;
+            break;
+        default:
+            psAbort("Should never get here --- checked params->combine before.\n");
+        }
+        stats->options |= combineStdev;
+    }
+
+    // We loop through each pixel in the output image.  We loop through each input readout.  We determine if
+    // that output pixel is contained in the image from that readout.  If so, we save it in psVector pixels.
+    // If not, we set a mask for that element in pixels.  Then, we mask off pixels not between fracLow and
+    // fracHigh.  Then we call the vector stats routine on those pixels/mask.  Then we set the output pixel
+    // value to the result of the stats call.
+
+    psVector *pixels = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of pixels
+    psF32 *pixelsData = pixels->data.F32; // Dereference pixels
+
+    psVector *mask   = psVectorAlloc(inputs->n, PS_TYPE_U8); // Mask for stack
+    psU8 *maskData = mask->data.U8;     // Dereference mask
+
+    psVector *weights = NULL;           // Stack of weights
+    psVector *errors = NULL;            // Stack of errors (sqrt of variance/weights), for psVectorStats
+    psF32 *weightsData = NULL;          // Dereference weights
+    if (params->weights) {
+        weights = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of weights
+        weightsData = weights->data.F32;
+    }
+    psVector *index = NULL;             // The indices to sort the pixels
+
+    float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of pixels to keep
+    if (keepFrac != 1.0 && first) {
+        psString comment = NULL;        // Comment to add to header
+        psStringAppend(&comment, "Min/max rejection: %f high, %f low, keep %d",
+                       params->fracHigh, params->fracLow, params->nKeep);
+        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
+        psFree(comment);
+    }
+
+    psMaskType maskVal = params->maskVal; // The mask value
+    if (maskVal && first) {
+        psString comment = NULL;        // Comment to add to header
+        psStringAppend(&comment, "Mask for combination: %x", maskVal);
+        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
+        psFree(comment);
+    }
+
+    #ifndef PS_NO_TRACE
+    psTrace("psModules.imcombine", 3, "Iterating output: %d --> %d, %d --> %d\n",
+            minInputCols - output->col0, maxInputCols - output->col0,
+            minInputRows - output->row0, maxInputRows - output->row0);
+    if (psTraceGetLevel("psModules.imcombine") >= 3) {
+        for (int r = 0; r < inputs->n; r++) {
+            pmReadout *readout = inputs->data[r]; // Input readout
+            psTrace("psModules.imcombine", 3, "Iterating input %d: %d --> %d, %d --> %d\n", r,
+                    minInputCols - readout->col0, maxInputCols - readout->col0,
+                    minInputRows - readout->row0, maxInputRows - readout->row0);
+        }
+    }
+    #endif
+
+    // Dereference output products
+    psF32 **outputImage  = output->image->data.F32; // Output image
+    psU8  **outputMask   = output->mask->data.U8; // Output mask
+    psF32 **outputWeight = NULL; // Output weight map
+    if (output->weight) {
+        outputWeight = output->weight->data.F32;
+    }
+
+    psVector *invScale = NULL;          // Inverse scale; pre-calculated for efficiency
+    if (scale) {
+        invScale = (psVector*)psBinaryOp(NULL, psScalarAlloc(1.0, PS_TYPE_F32), "/", (const psPtr)scale);
+    }
+
+    for (int i = minInputRows; i < maxInputRows; i++) {
+        int yOut = i - output->row0; // y position on output readout
+        #ifdef SHOW_BUSY
+
+        if (psTraceGetLevel("psModules.imcombine") > 9) {
+            printf("Processing row %d\r", i);
+            fflush(stdout);
+        }
+        #endif
+        for (int j = minInputCols; j < maxInputCols; j++) {
+            int xOut = j - output->col0; // x position on output readout
+
+            int numValid = 0;           // Number of valid pixels in the stack
+            memset(maskData, 0, mask->n * sizeof(psU8)); // Reset the mask
+            for (int r = 0; r < inputs->n; r++) {
+                pmReadout *readout = inputs->data[r]; // Input readout
+                int yIn = i - readout->row0; // y position on input readout
+                int xIn = j - readout->col0; // x position on input readout
+                psImage *image = readout->image; // The readout image
+
+                #if 0 // This should have been taken care of already:
+                // Check bounds
+                if (xIn < 0 || xIn >= image->numCols || yIn < 0 || yIn >= image->numRows) {
+                    continue;
+                }
+                #endif
+
+                pixelsData[r] = image->data.F32[yIn][xIn];
+                if (!isfinite(pixelsData[r])) {
+                    maskData[r] = 1;
+                    continue;
+                }
+
+                // Check mask
+                psImage *roMask = readout->mask; // The mask image
+                if (roMask && roMask->data.U8[yIn][xIn] & maskVal) {
+                    maskData[r] = 1;
+                    continue;
+                }
+
+                if (params->weights) {
+                    weightsData[r] = readout->weight->data.F32[yIn][xIn];
+                }
+
+                if (zero) {
+                    pixelsData[r] -= zero->data.F32[r];
+                }
+                if (scale) {
+                    pixelsData[r] *= invScale->data.F32[r];
+                    if (params->weights) {
+                        weightsData[r] *= invScale->data.F32[r] * invScale->data.F32[r];
+                    }
+                }
+
+                numValid++;
+            }
+
+            if (numValid == 0) {
+                outputMask[yOut][xOut] = params->blank;
+                outputImage[yOut][xOut] = NAN;
+                continue;
+            }
+
+            // Apply fracLow,fracHigh if there are enough pixels
+            if (numValid * keepFrac >= params->nKeep && keepFrac != 1.0) {
+                index = psVectorSortIndex(index, pixels);
+                int numLow = numValid * params->fracLow; // Number of low pixels to clip
+                int numHigh = numValid * params->fracHigh; // Number of high pixels to clip
+                // Low pixels
+                psS32 *indexData = index->data.S32; // Dereference index
+                for (int k = 0, numMasked = 0; numMasked < numLow && k < index->n; k++) {
+                    // Don't count the ones that are already masked
+                    if (!maskData[indexData[k]]) {
+                        maskData[indexData[k]] = 1;
+                        numMasked++;
+                    }
+                }
+                // High pixels
+                for (int k = pixels->n - 1, numMasked = 0; numMasked < numHigh && k >= 0; k--) {
+                    // Don't count the ones that are already masked
+                    if (! maskData[indexData[k]]) {
+                        maskData[indexData[k]] = 1;
+                        numMasked++;
+                    }
+                }
+            }
+
+            // XXXXX this step probably is very expensive : convert errors to variance everywhere?
+            if (params->weights) {
+                errors = (psVector*)psUnaryOp(errors, weights, "sqrt");
+            }
+
+            // Combination
+            if (!psVectorStats(stats, pixels, errors, mask, 1)) {
+                // Can't do much about it, but it's not worth worrying about
+                psErrorClear();
+                outputImage[yOut][xOut] = NAN;
+                outputMask[yOut][xOut] = params->blank;
+                if (params->weights) {
+                    outputWeight[yOut][xOut] = NAN;
+                }
+            } else {
+                outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine);
+                outputMask[yOut][xOut] = isfinite(outputImage[yOut][xOut]) ? 0 : params->blank;
+                if (params->weights) {
+                    float stdev = psStatsGetValue(stats, combineStdev);
+                    outputWeight[yOut][xOut] = PS_SQR(stdev); // Variance
+                    // XXXX this is not the correct formal error.
+                    // also, the weighted mean is not obviously the correct thing here
+                }
+            }
+        }
+    }
+    #ifdef SHOW_BUSY
+    if (psTraceGetLevel("psModules.imcombine") > 9) {
+        printf("\n");
+    }
+    #endif
+    psFree(index);
+    psFree(pixels);
+    psFree(mask);
+    psFree(weights);
+    psFree(errors);
+    psFree(stats);
+    psFree(invScale);
+
+    // Update the "concepts"
+    psList *inputCells = psListAlloc(NULL); // List of cells
+    for (long i = 0; i < inputs->n; i++) {
+        pmReadout *readout = inputs->data[i]; // Readout of interest
+        psListAdd(inputCells, PS_LIST_TAIL, readout->parent);
+    }
+    bool success = pmConceptsAverageCells(output->parent, inputCells, NULL, NULL, true);
+    psFree(inputCells);
+
+    output->data_exists = true;
+    output->parent->data_exists = true;
+    output->parent->parent->data_exists = true;
+
+    return success;
+}
+