Changeset 36075

Timestamp:

Aug 30, 2013, 4:55:55 PM (13 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/eam_branches/ipp-20130711

Files:

• . (modified) (1 prop)
• PS-IPP-PStamp/lib/PS/IPP/PStamp/Job.pm (modified) (16 diffs)
• PS-IPP-PStamp/lib/PS/IPP/PStamp/RequestFile.pm (modified) (4 diffs)
• ppSub/src/ppSubMatchPSFs.c (modified) (2 diffs)
• psModules (modified) (1 prop)
• psModules/src/imcombine/pmSubtractionMatch.c (modified) (1 diff)
• psModules/src/imcombine/pmSubtractionSimple.c (modified) (7 diffs)
• psModules/src/imcombine/pmSubtractionSimple.h (modified) (1 diff)
• psModules/src/objects/pmPSFtryFitPSF.c (modified) (1 diff)

branches/eam_branches/ipp-20130711/PS-IPP-PStamp/lib/PS/IPP/PStamp/Job.pm

@@ -34 +34 @@
 use POSIX;
 use Time::HiRes qw(gettimeofday);
+use IO::Handle;
 
 my $dvo_verbose = 0;
@@ -176 +177 @@
         $need_magic, $dateobs_begin, $dateobs_end, $filter, $data_group, $option_mask, $mjd_min, $mjd_max, 
         0, 0,   # fwhm cuts are not applied here
-        undef, undef, undef, $verbose);
+        undef, 0, # no cam run information
+        undef, undef, $verbose);
 
     return $results;
@@ -243 +245 @@
     if (isnull($filter)) {
         $filter = undef;
+    }
+
+    if ($stage eq 'stack_summary') {
+        # stack_summary jobs are so different from others that we lookup in a specialized function
+        return lookup_stack_summary($ipprc, $row, $imagedb, $tess_id, $component, $filter, $data_group, $option_mask, $mjd_min, $mjd_max, $verbose);
     }
 
@@ -322 +329 @@
         0, $dateobs_begin, $dateobs_end, $filter, $data_group, $option_mask, $mjd_min, $mjd_max,
         getValOrZero($row->{FWHM_MIN}), getValOrZero($row->{FWHM_MAX}),
-        undef, $release_name, $survey,
+        undef, 0, $release_name, $survey,
         $verbose);
 
@@ -351 +358 @@
     my $fwhm_max = shift;
     my $selectedAstrom = shift;
+    my $selected_cam_id = shift;
     my $release_name = shift;
     my $survey = shift;
@@ -602 +610 @@
         # The image selectors are such that multiple runs my have be returned for the same exposure.
         # Return only the latest one.
-        $images = filterRuns($stage, $need_magic, $images, $inverse, $verbose);
+        $images = filterRuns($stage, $choose_components, $need_magic, $images, $inverse, $verbose);
     }
 
@@ -699 +707 @@
                 $out->{astrom} = $selectedAstrom;
                 $out->{cam_path_base} = $selectedAstrom;
+                $out->{cam_id} = $selected_cam_id if $selected_cam_id;
                 if ($selectedAstrom =~ /\.smf$/) {
                     $out->{cam_path_base} =~ s/\.smf$//;
@@ -731 +740 @@
 
     return $output;
+}
+
+sub lookup_stack_summary {
+    my ($ipprc, $row, $imagedb, $tess_id, $component, $filter, $data_group, $option_mask, $mjd_min, $mjd_max, $verbose) = @_;
+
+    $row->{error_code} = $PSTAMP_NO_IMAGE_MATCH;
+
+    $component = '' if $component eq 'all';
+
+    my $req_type = $row->{REQ_TYPE};
+    my $results;
+    if ($req_type eq 'byid') {
+        my $sass_id = $row->{ID};
+        my $command = "$stacktool -dbname $imagedb -summary -sass_id $sass_id";
+        $results = runToolAndParse($command, $verbose);
+    } else {
+        my ($release_name, $survey, $default_tess_id) = get_release_info($row);
+        if (!$tess_id and $default_tess_id) {
+            $tess_id = $default_tess_id;
+        }
+
+        my $skycells;
+        if ($req_type eq 'bycoord') {
+            my $ra = $row->{CENTER_X};
+            my $dec = $row->{CENTER_Y};
+            $skycells = lookup_skycell_by_coords($ipprc, $tess_id, $component, $ra, $dec, $verbose);
+        } elsif ($req_type eq 'byskycell') {
+            if (!defined $tess_id or !defined $component) {
+                print STDERR "Error: TESS_ID and  COMPONENT are required for REQ_TYPE byskycell\n";
+                $row->{error_code} = $PSTAMP_INVALID_REQUEST;
+            }
+            my $skycell = {
+                tess_id    => $tess_id,
+                component => $component
+            };
+            push @$skycells, $skycell
+        } else {
+            print STDERR "Error: $req_type is not a valid REQ_TYPE for IMG_TYPE stack_summary\n";
+            $row->{error_code} = $PSTAMP_INVALID_REQUEST;
+            return undef;
+        }
+
+        if ($skycells) {
+            # We have a list of skycells. Find matching projection cells.
+            # and then matching rows in stackSummary
+
+            my %projection_cells_found;
+            foreach my $skycell (@$skycells) {
+                my $tess_id = $skycell->{tess_id};
+                my $skycell_id = $skycell->{component};
+
+                my $projection_cell = skycell_id_to_projection_cell($skycell_id);
+
+                # only need to include a projection cell once
+                next if $projection_cells_found{$projection_cell};
+
+                $projection_cells_found{$projection_cell} = 1;
+
+                # print "$tess_id $skycell_id $projection_cell\n";
+
+                my $command;
+                if ($release_name or $survey) {
+                    $command = "$releasetool -dbname $imagedb";
+                    if ($release_name) {
+                        $command .= " -release_name $release_name";
+                    } else {
+                        $command .= " -priority_order";
+                    }
+                    $command .= " -surveyName $survey";
+                } else {
+                    $command = "$stacktool -dbname $imagedb";
+                }
+                $command .= " -summary -tess_id $tess_id";
+                $command .= " -projection_cell $projection_cell";
+
+                $command .= " -filter $filter" if $filter;
+                $command .= " -data_group $data_group" if $data_group;
+
+                my $these_results = runToolAndParse($command, $verbose);
+                push @$results, @$these_results if $these_results;
+            }
+        }
+    }
+
+    if ($results) {
+        foreach my $summary (@$results) {
+            my $path_base = $summary->{path_base};
+            my $projection_cell = $summary->{projection_cell};
+
+            # we need to flesh out the hashes returned with values that the parser expects
+            $summary->{stage} = 'stack_summary';
+            $summary->{component} = $projection_cell;
+            $summary->{stage_id} = $summary->{sass_id};
+            $summary->{path_base} = "$path_base.$projection_cell";
+
+            # The stack_summary stage doesn't follow the usual ipp conventions for file rules
+            # and path_base. We leave it up to the program that runs the job to handle this.
+            # However, the parser uses the image parameter to set the job's outputBase (by removing the .fits)
+            # so we need to set it. Just use the new path_base.
+            $summary->{image} = $summary->{path_base};
+            $summary->{state} = 'full';
+            $summary->{data_state} = 'full';
+            $summary->{fault} = 0;
+            $summary->{imagedb} = $imagedb;
+            $summary->{row_index} = [];
+            push @{$summary->{row_index}}, 0;
+
+            # XXX: Are there any other required parameters? 
+            # TODO: create a function to validate the components returned
+            # and make sure that any required elements are there.
+        }
+    }
+
+    return $results;
 }
 
@@ -946 +1069 @@
                     $dateobs_begin, $dateobs_end, $filter, $data_group, $option_mask, $mjd_min, $mjd_max,
                     $fwhm_min, $fwhm_max,
-                    $chip->{astrom}, undef, undef,
+                    $chip->{astrom}, $chip->{cam_id}, undef, undef,
                     $verbose);
 
@@ -970 +1093 @@
                     $dateobs_begin, $dateobs_end, $filter, $data_group, $option_mask, $mjd_min, $mjd_max, 
                     $fwhm_min, $fwhm_max,
-                    undef, $release_name, $survey,
+                    undef, 0, $release_name, $survey,
                     $verbose);
 
@@ -1166 +1289 @@
 
     my @lines;
-    {
+    my $looked_up_fast;
+    if ($requested_tess_id) {
+        # try using the last method first
+        $looked_up_fast = lookup_skycells_fast($ipprc, \@lines, $requested_tess_id, $ra, $dec, $verbose);
+    } 
+    
+    if (!$looked_up_fast) {
         my $command = "$whichimage $ra $dec";
         $command .= " --tess_id $requested_tess_id" if $requested_tess_id;
@@ -1302 +1431 @@
     $image->{astrom} = $astromFile;
     $image->{cam_path_base} = $camRoot;
+    $image->{cam_id} = $camRun->{cam_id};
 
     return 1;
@@ -1649 +1779 @@
 sub filterRuns {
     my $stage      = shift;
+    my $drop_duplicate_runs = shift;
     my $need_magic = shift;
     my $inputs     = shift;
@@ -1694 +1825 @@
         }
 
-        $printed = 1 if (($exp_id == $last_exp_id) and ($run_id == $last_run_id));
+        # $printed = 1 if (($exp_id == $last_exp_id) and ($run_id == $last_run_id));
 
         # skip if we need magicked run and this one has never been magicked
@@ -1702 +1833 @@
         }
 
-        if (($exp_id == $last_exp_id) and ($run_id != $last_run_id)) {
+        if (($exp_id == $last_exp_id) and ($drop_duplicate_runs || ($run_id != $last_run_id))) {
             print "Skipping duplicate ${stage}Run $run_id for $exp_id\n" if !$printed;
             next;
@@ -1873 +2004 @@
     return $diff_mode;
 }
+
+sub setErrorCodesForRows {
+    my $rowList = shift;
+    my $error_code = shift;
+
+    foreach my $row (@$rowList) {
+        $row->{error_code} = $error_code;
+    }
+
+}
+
+# Convert from skycell_id to projection_cell
+# This code is not particularly elegant but I think that it matches the SQL used
+# to create stackAssocations does. 
+
+sub skycell_id_to_projection_cell {
+    my $skycell_id = shift;
+
+    # split skycell_id into '.' separated components
+    my @strs = split '\.', $skycell_id;
+
+    my $num = scalar @strs;
+
+    # projection_cell is 'skycell_id' plus all .num values except the last one
+    my $projection_cell = $strs[0];
+    for (my $i=1; $i < $num - 1; $i++) {
+        $projection_cell .= ".$strs[$i]";
+    }
+
+    return $projection_cell;
+}
+
+# lookup_skycells_fast() 
+# Lookup using a "skycell server" process which is implemented by running dvoImagesAtCoords in "server" mode.
+# Repeated calls to this function will reuse an existing process. This avoids having to loading the
+# tessellation over and over again for each skycell lookup.
+# When this parse process exists, the pipes are closed which causes the skycell server to exit.
+
+my $scs_tess_id;    # tess_id for existing scs (skycell server) (if any)
+my $scs_out;        # reference filehandle for sending data to the scs
+my $scs_in;         # filehandle for response from scs
+
+sub lookup_skycells_fast {
+    my ($ipprc, $results, $tess_id, $ra, $dec, $verbose) = @_;
+
+    # see if we need to start the scs (skycell server)
+
+    if (!$scs_tess_id or ($scs_tess_id ne $tess_id)) {
+        if ($scs_tess_id) {
+            # tess_id has changed, close down existing scs
+            close $scs_out;
+            $scs_out = undef;
+            close $scs_in;
+            $scs_in = undef;
+            $scs_tess_id = undef;
+        }
+
+        # convert tess_id to a directory name
+        my $tess_dir = $ipprc->tessellation_catdir( $tess_id );
+        unless ($tess_dir) {
+            print STDERR "Unrecognized tess_id: $tess_id\n";
+            return 0;
+        }
+
+        # convert tess_dir to an absolute directory name
+        my $tess_dir_resolved = $ipprc->convert_filename_absolute( $tess_dir );
+        unless ($tess_dir_resolved and -d $tess_dir_resolved) {
+            print STDERR "failed to resolve tessellation directory for $tess_id\n";
+            return 0;
+        }
+
+        # start an scs
+
+        # create pipes for communicating with the server
+        pipe $scs_in, SERVER_WRITER;
+        pipe SERVER_READER, $scs_out;
+
+        # set our output to be unbuffered
+        $scs_out->autoflush(1);
+
+        # FORK
+        my $scs_pid = fork();
+
+        if (!$scs_pid) {
+            unless (defined $scs_pid) {
+                # This is still in parent process. No joy.
+                print STDERR "fork of skycell server failed: $!";
+                return 0;
+            }
+
+            # This code is running in the child process - the skycell server
+            # Close file handles for the parent's ends of the pipes
+            close $scs_in;
+            close $scs_out;
+            # close stdio filehandles. we are about to redirect them
+            # (STDERR stays open pointing to the parser log)
+            close STDIN;
+            close STDOUT;
+
+            # redirect our stdio file handles to the proper end of the pipes
+            open (STDIN,  "<&SERVER_READER")   or die "\nSCS: failed to redirect stdin";
+            open (STDOUT, ">>&SERVER_WRITER")  or die "\nSCS: failed to redirect stdout";
+
+            # All set. Now exec dvoImagesAtCoords
+
+            my $command = "$dvoImagesAtCoords -D CATDIR $tess_dir_resolved -coords -";
+
+            print STDERR "SCS: execing $command\n";
+
+            unless(exec $command) {
+                # note parent will notice that we are gone due to the pipes getting broken
+                # when this child dies.
+                die "SCS: failed to exec $command";
+            }
+        }
+
+        # This code is running in the parent process (the parser).
+        # We have succesfully launched the scs process.
+        # Close file handles for the scs' end of the pipes.
+        close SERVER_WRITER;
+        close SERVER_READER;
+
+        # remember the tess_id 
+        $scs_tess_id = $tess_id;
+    }
+
+    # send coordinates to the skycell server and wait for the results 
+    # (which come nearly instantly which is the point of doing this)
+
+    # If pipes to the skycell server die, don't abort. 
+    # Our reads and writes detect errors and act sensibly.
+    local $SIG{PIPE} = 'IGNORE';
+
+    # write the coordinates to the pipe
+    my $write_ok = print $scs_out "1 $ra $dec\n";
+    if (!$write_ok) {
+        # server process probably didn't start properly or died. Fail.
+        # Caller will attempt to use the conventional method.
+        print STDERR "write to skycell server failed. Error is $!\n";
+        return 0;
+    }
+
+    my $very_verbose = 0;
+    print STDERR "  sent coordinates to skycell server\n" if $very_verbose;
+
+    my $received_done = 0;
+
+    # Wait for response. If pipe breaks the read returns with no data.
+    while (my $line = <$scs_in>) {
+        chomp $line;
+        if ($line eq 'DONE') {
+            # No more results for these coordinates.
+            $received_done = 1;
+            print STDERR "    received DONE\n" if $very_verbose;
+            last;
+        }
+        print STDERR "    received $line\n" if $very_verbose;
+
+        # The caller expects the output to be in the format used by whichimage
+        # which omits the point number and includes the tess_id.
+        my ($ptnum, $raout, $decout, $skycell_id) = split " ", $line;
+
+        my $out = "$raout $decout $tess_id $skycell_id";
+
+        # print the result to the log
+        print "$out\n" if $verbose;
+
+        push @$results, $out;
+    }
+
+    print STDERR "Out of wait for responses loop\n" if $very_verbose;
+
+    if (!$received_done) {
+        # something has gone wrong.
+        print STDERR "Read loop exited without receiving DONE.\n";
+
+        die "BAILING out" if $very_verbose;
+
+        # forget about the existing server instance
+        $scs_tess_id = undef;
+        close $scs_out;
+        close $scs_in;
+
+        # drop any results received so far
+        @$results = ();
+        return 0;
+    }
+
+    return 1;
+}
+
         
 

branches/eam_branches/ipp-20130711/PS-IPP-PStamp/lib/PS/IPP/PStamp/RequestFile.pm

@@ -29 +29 @@
                     $PSTAMP_SELECT_BACKMDL
                     $PSTAMP_SELECT_JPEG
+                    $PSTAMP_SELECT_EXP
+                    $PSTAMP_SELECT_NUM
                     $PSTAMP_SELECT_UNCOMPRESSED
                     $PSTAMP_SELECT_INVERSE
@@ -71 +73 @@
 our $PSTAMP_SELECT_BACKMDL   = 32;
 our $PSTAMP_SELECT_JPEG      = 64;
-# unused 128
-# unused 256
+our $PSTAMP_SELECT_EXP       = 128;
+our $PSTAMP_SELECT_NUM       = 256;
 our $PSTAMP_SELECT_UNCOMPRESSED = 512;
 our $PSTAMP_SELECT_INVERSE      = 1024;
@@ -97 +99 @@
 our $PSTAMP_INVALID_REQUEST  = 21;
 our $PSTAMP_UNKNOWN_PROJECT  = 22;
-our $PSTAMP_UNKNOWN_PRODUCT  = 22;  #this error code was mis-named it is left for compatabiliyt
+our $PSTAMP_UNKNOWN_PRODUCT  = 22;  #this error code was a typo it is left for compatabiliy
 our $PSTAMP_NO_IMAGE_MATCH   = 23;
 our $PSTAMP_NOT_DESTREAKED   = 24;
 our $PSTAMP_NOT_AVAILABLE    = 25;
-our $PSTAMP_GONE             = 26;  # this value is used in ippTools
+our $PSTAMP_GONE             = 26;  # this value is also used in ippTools
 our $PSTAMP_NO_JOBS_QUEUED   = 27;
 our $PSTAMP_NO_OVERLAP       = 28;
 our $PSTAMP_NOT_AUTHORIZED   = 29;
+our $PSTAMP_NO_VALID_PIXELS  = 30;
 
 
@@ -143 +146 @@
 PSTAMP_NO_OVERLAP
 PSTAMP_NOT_AUTHORIZED
-PSTAMP_NOT_AUTHORIZED
+PSTAMP_NO_VALID_PIXELS
 );
 

branches/eam_branches/ipp-20130711/ppSub/src/ppSubMatchPSFs.c

@@ -57 +57 @@
     if (!pmFPACopy(photFile->fpa, ro->parent->parent->parent)) {
         psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry");
-        return false;
+        return NAN;
     }
 
@@ -76 +76 @@
         psWarning("Unable to determine PSF.");
         psFree(view);
-        return true;
+        return NAN;
     }
     psFree(view);

branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionMatch.c

@@ -540 +540 @@
     // Bail here if we're doing the simple matching
     if (type == PM_SUBTRACTION_KERNEL_SIMPLE) {
-      if (!pmSubtractionSimpleMatch(conv1,conv2,ro1,ro2,sources,size,maskVal,maskBad,maskPoor)) {
+      if (!pmSubtractionSimpleMatch(conv1,conv2,ro1,ro2,sources,size,maskVal,maskBad,maskPoor,optThreshold)) {
         return false;
       }

branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.c

@@ -63 +63 @@
 }  
 
+bool simple_apply_mask(psImage *image, psImage *weight, psImage *mask,
+                       psImageMaskType maskVal) {
+  for (int y = 0; y < mask->numRows; y++) {
+    for (int x = 0; x < mask->numCols; x++) {
+      if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) {
+        image->data.F32[y][x] = NAN;
+        if (weight) {
+          weight->data.F32[y][x] = NAN;
+        }
+      }
+    }
+  }
+  return(true);
+}
+                       
+
+// Copied from pmSubtraction
+static void solvedKernelPreCalc(psKernel *kernel, // Kernel, updated
+                                const pmSubtractionKernels *kernels, // Kernel basis functions
+                                float value,                         // Normalisation value for basis function
+                                int index                  // Index of basis function of interest
+                                )
+{
+  int size = kernels->size;           // Kernel half-size
+  pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated values
+  for (int v = -size; v <= size; v++) {
+    for (int u = -size; u <= size; u++) {
+      kernel->kernel[v][u] +=  value * preCalc->kernel->kernel[v][u];
+    }
+  }
+
+  return;
+}
+//End copy
+
 bool pmSubtractionSimpleMatch(pmReadout *conv1,
                               pmReadout *conv2,
@@ -71 +106 @@
                               psImageMaskType maskVal,
                               psImageMaskType maskBad,
-                              psImageMaskType maskPoor
+                              psImageMaskType maskPoor,
+                              float deconvolveThreshold
                               ) {
   //
   // We've already validated the input values at this level
   float sig2fwhm = 2.0 * sqrt(2.0 * log(2.0));
-  float fwhm1,fwhm2;
-  float sigma1,sigma2,sigmaKern;
+  float fwhm1 = 0,fwhm2 = 0;
+  float sigma1 = 0,sigma2 = 0,sigmaKern = 0;
   float chisq = 1.0;
   int convolution_direction = 0;
-  psImage *image1;
-  psImage *mask1;
-  psImage *var1;
-
-  psImage *image2;
-  psImage *mask2;
-  psImage *var2;
-
-  psImage *imageC1;
-  psImage *maskC1;
-  psImage *varC1;
-
-  psImage *imageC2;
-  psImage *maskC2;
-  psImage *varC2;
-
+  psImage *image1 = NULL;
+  psImage *mask1 = NULL;
+  psImage *var1 = NULL;
+
+  psImage *image2 = NULL;
+  psImage *mask2 = NULL;
+  psImage *var2 = NULL;
+
+  psImage *imageC1 = NULL;
+  psImage *maskC1 = NULL;
+  psImage *varC1 = NULL;
+
+  psImage *imageC2 = NULL;
+  psImage *maskC2 = NULL;
+  psImage *varC2 = NULL;
+
+  psImage *maskTemp = NULL;
+  
   // Allocate images, as this is usually done by subtractionMatchAlloc after this function is called.  
   int numCols = ro1->image->numCols;
@@ -101 +139 @@
 
   if (conv1) {
-    conv1->covariance = psMemIncrRefCounter(ro1->covariance);
+    //    conv1->covariance = psMemIncrRefCounter(ro1->covariance);
     if (!conv1->image) {
       conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
@@ -118 +156 @@
   }
   if (conv2) {
-    conv2->covariance = psMemIncrRefCounter(ro2->covariance);
+    //    
     if (!conv2->image) {
       conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
@@ -171 +209 @@
   if (!conv1) {
     if (convolution_direction == 1) {
-      chisq = 100;
-    }
-    convolution_direction = 2;
+      if (sigma1 - sigma2 > deconvolveThreshold) {
+        chisq = 100;
+      }
+    }
+    //    convolution_direction = 2;
   }
   if (!conv2) {
     if (convolution_direction == 2) {
-      chisq = 100;
-    }
-    convolution_direction = 1;
+      if (sigma2 - sigma1 > deconvolveThreshold) {
+        chisq = 100;
+      }
+    }
+    //    convolution_direction = 1;
   }
   
@@ -186 +228 @@
   int maskBox = (int) ceil(sigmaKern * 1.1774); // diameter is 1/2 FWHM
   int maskBlank = 8;  // I should be able to get this from a reference, right?
-
-  //
-  // Construct required kernel.  No longer needed as we can direct convolve
-  //  psVector *kernelVec = pmSubtractionKernelSIMPLE(sigmaKern,0,size); // This is normalized to unity.
-  //  psFree(kernelVec);
-
-  //
-  // Do convolutions
-  if (convolution_direction == 1) {
-    psImageSmoothMask_Threaded(imageC1,image1,mask1,maskVal,sigmaKern,6,1e-6);
-    psImageSmoothMask_Threaded(varC1,var1,mask1,maskVal,sigmaKern * M_SQRT1_2,6,1e-6);
-    maskC1 = psImageConvolveMask(maskC1,mask1,maskVal,maskBad,
-                                 -maskBox,maskBox,-maskBox,maskBox);
-    if (conv2) {
-      imageC2 = psImageCopy(imageC2,image2,PS_TYPE_F32);
-      varC2   = psImageCopy(varC2,var2,PS_TYPE_F32);
-      maskC2  = psImageCopy(maskC2,mask2,PS_TYPE_IMAGE_MASK);
-    }
-    pmSubtractionBorder(imageC1,varC1,maskC1,maskBox,maskBlank);
-    pmSubtractionMaskApply(imageC1,varC1,maskC1,PM_SUBTRACTION_MODE_1);
-  }
-  else if (convolution_direction == 2) {
-    psImageSmoothMask_Threaded(imageC2,image2,mask2,maskVal,sigmaKern,6,1e-6);
-    psImageSmoothMask_Threaded(varC2,var2,mask2,maskVal,sigmaKern * M_SQRT1_2,6,1e-6);
-    maskC2 = psImageConvolveMask(maskC2,mask2,maskVal,maskBad,
-                                 -maskBox,maskBox,-maskBox,maskBox);
-    if (conv1) {
-      imageC1 = psImageCopy(imageC1,image1,PS_TYPE_F32);
-      varC1   = psImageCopy(varC1,var1,PS_TYPE_F32);
-      maskC1  = psImageCopy(maskC1,mask1,PS_TYPE_IMAGE_MASK);
-    }
-    pmSubtractionBorder(imageC2,varC2,maskC2,maskBox,maskBlank);
-    pmSubtractionMaskApply(imageC2,varC2,maskC2,PM_SUBTRACTION_MODE_2);
-  }    
-
-  //
-  // Do normalization
-  float normalization = 1.0;
-
-  // Scan source list, do box photometry on peaks, and then solve the linear relation.
-  int photRadius = (int) floor(PS_MAX(sigma1,sigma2) * 2.0 * sqrt(2.0 * log(2.0))); // Go out a FWHM diameter from the center.
-  psVector *logFluxDifferences = psVectorAlloc(sources->n,PS_TYPE_F32);
-  psVector *fitMask = psVectorAlloc(sources->n,PS_TYPE_VECTOR_MASK);
-  for (int i = 0; i < sources->n; i++) {
-    pmSource *source = sources->data[i];
-    int nPix1,nPix2;
-    float flux1,flux2;
-
-    if (convolution_direction == 1) {
-      simple_do_boxphot(&nPix1,&flux1,source,imageC1,maskC1,maskBad,photRadius);
-      if (conv2) {
-        simple_do_boxphot(&nPix2,&flux2,source,imageC2,maskC2,maskBad,photRadius);
-      }
-      else {
-        simple_do_boxphot(&nPix2,&flux2,source,image2,mask2,maskBad,photRadius);
-      }
-    }
-    else if (convolution_direction == 2) {
-      simple_do_boxphot(&nPix2,&flux2,source,imageC2,maskC2,maskBad,photRadius);
-      if (conv1) {
-        simple_do_boxphot(&nPix1,&flux1,source,imageC1,maskC1,maskBad,photRadius);
-      }
-      else {
-        simple_do_boxphot(&nPix1,&flux1,source,image1,mask1,maskBad,photRadius);
-      }
-    }
-    logFluxDifferences->data.F32[i] = flux2 - flux1;
-    fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
-    if ((PS_MIN(nPix1,nPix2) <= 0.75 * PS_MAX(nPix1,nPix2))||
-        (!isfinite(flux1))||(!isfinite(flux2))) {
-      fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
-    }
-
-    //    fprintf(stderr,"SOURCES: %d %g %g %g -> %d %d %g %g %d %g\n",i,source->peak->xf,source->peak->yf,source->psfMag,
-    //      nPix1,nPix2,flux1,flux2,fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i],logFluxDifferences->data.F32[i]);
-    
-  }
-
-  // Given the differences in log-flux space, the normalization factor is just the exponential of the median difference
-  psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-  if (!psVectorStats(stats,logFluxDifferences,NULL,fitMask,0xff)) {
-    // This should complain.
-    normalization = 1.0;
-  }
-
-  normalization = pow(10,stats->robustMedian);
-  // fprintf(stderr,"NORM: %g+/-%g\n",stats->robustMedian,stats->robustStdev);
-  
-  psFree(stats);
-  psFree(logFluxDifferences);
-  psFree(fitMask);
-
-  // Apply normalization
-  if (normalization != 1.0) {
-    if ((conv1)&&((convolution_direction == 1))) {
-      psBinaryOp(imageC1,imageC1,"*",psScalarAlloc((float) normalization, PS_TYPE_F32));
-      psBinaryOp(varC1,varC1,"*",psScalarAlloc((float) PS_SQR(normalization), PS_TYPE_F32));
-    }
-    else if ((conv2)&&(convolution_direction == 2)) {
-      normalization = 1.0 / normalization; // Because we fit one way, but are using it in the other.
-      psBinaryOp(imageC2,imageC2,"*",psScalarAlloc((float) normalization, PS_TYPE_F32));
-      psBinaryOp(varC2,varC2,"*",psScalarAlloc((float) PS_SQR(normalization), PS_TYPE_F32));
-    }
-  }
-  
-
+  int maskPoorVal = 16384; // Another value that should be found elsewhere.
   //
   // Make a fake pmSubtractionKernels element so we can add it appropriately.
-  // I call it fake because we've successfully done everything at this point
-  // without having to define these things.
+  // Defining everything here is a bit clunky, but it's necessary to get the covariance
+  // correct.
   psVector *fwhms = psVectorAlloc(1,PS_TYPE_F32);
   fwhms->data.F32[0] = sigmaKern * sig2fwhm;
@@ -329 +266 @@
   kernels->numStamps = sources->n;
   
+  psKernel *kernel = psKernelAlloc(-size,size,-size,size);
+  solvedKernelPreCalc(kernel,kernels,1.0,0);
+  
+  //
+  // Do convolutions
+  if (convolution_direction == 1) {
+    if (conv1) {
+      psImageSmoothMask_Threaded(imageC1,image1,mask1,maskVal,sigmaKern,6,1e-6);
+      psImageSmoothMask_Threaded(varC1,var1,mask1,maskVal,sigmaKern * M_SQRT1_2,6,1e-6);
+
+      maskTemp = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+      maskTemp = psImageConvolveMask(maskTemp,mask1,maskVal,maskBad,      // Mask bad values
+                                   -maskBox,maskBox,-maskBox,maskBox);
+      maskC1 = psImageConvolveMask(maskC1,maskTemp,maskPoorVal,maskPoor,  // Mask poor values
+                                   -maskBox,maskBox,-maskBox,maskBox);
+      psFree(maskTemp);
+
+      conv1->covariance = psImageCovarianceCalculate(kernel,ro1->covariance);
+      pmSubtractionBorder(imageC1,varC1,maskC1,maskBox,maskBlank);
+      simple_apply_mask(imageC1,varC1,maskC1,maskBad);
+    }
+    if (conv2) {
+      imageC2 = psImageCopy(imageC2,image2,PS_TYPE_F32);
+      varC2   = psImageCopy(varC2,var2,PS_TYPE_F32);
+      maskC2  = psImageCopy(maskC2,mask2,PS_TYPE_IMAGE_MASK);
+      conv2->covariance = psMemIncrRefCounter(ro2->covariance);
+    }
+  }
+  else if (convolution_direction == 2) {
+    if (conv2) {
+      psImageSmoothMask_Threaded(imageC2,image2,mask2,maskVal,sigmaKern,6,1e-6);
+      psImageSmoothMask_Threaded(varC2,var2,mask2,maskVal,sigmaKern * M_SQRT1_2,6,1e-6);
+
+      maskTemp = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+      maskTemp = psImageConvolveMask(maskTemp,mask2,maskVal,maskBad,      // Mask bad values
+                                   -maskBox,maskBox,-maskBox,maskBox);
+      maskC2 = psImageConvolveMask(maskC2,maskTemp,maskPoorVal,maskPoor,  // Mask poor values
+                                   -maskBox,maskBox,-maskBox,maskBox);
+      psFree(maskTemp);
+
+      conv2->covariance = psImageCovarianceCalculate(kernel,ro2->covariance);
+      pmSubtractionBorder(imageC2,varC2,maskC2,maskBox,maskBlank);
+      simple_apply_mask(imageC2,varC2,maskC2,maskBad);
+    }
+    if (conv1) {
+      imageC1 = psImageCopy(imageC1,image1,PS_TYPE_F32);
+      varC1   = psImageCopy(varC1,var1,PS_TYPE_F32);
+      maskC1  = psImageCopy(maskC1,mask1,PS_TYPE_IMAGE_MASK);
+      conv1->covariance = psMemIncrRefCounter(ro1->covariance);
+    }
+  }    
+
+  psFree(kernel); // No longer needed after doing covariance calculation
+
+  //
+  // Do normalization
+  float normalization = 1.0;
+
+  // Scan source list, do box photometry on peaks, and then solve the linear relation.
+  int photRadius = (int) floor(PS_MAX(sigma1,sigma2) * 2.0 * sqrt(2.0 * log(2.0))); // Go out a FWHM diameter from the center.
+  psVector *logFluxDifferences = psVectorAlloc(sources->n,PS_TYPE_F32);
+  psVector *fitMask = psVectorAlloc(sources->n,PS_TYPE_VECTOR_MASK);
+  for (int i = 0; i < sources->n; i++) {
+    pmSource *source = sources->data[i];
+    int nPix1,nPix2;
+    float flux1,flux2;
+
+    if (conv1) {
+      simple_do_boxphot(&nPix1,&flux1,source,imageC1,maskC1,maskBad,photRadius);
+    }
+    else {
+      simple_do_boxphot(&nPix1,&flux1,source,image1,mask1,maskBad,photRadius);
+    }
+
+    if (conv2) {
+      simple_do_boxphot(&nPix2,&flux2,source,imageC2,maskC2,maskBad,photRadius);
+    }
+    else {
+      simple_do_boxphot(&nPix2,&flux2,source,image2,mask2,maskBad,photRadius);
+    }
+
+    logFluxDifferences->data.F32[i] = flux2 - flux1;
+    fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
+    if ((PS_MIN(nPix1,nPix2) <= 0.75 * PS_MAX(nPix1,nPix2))||
+        (!isfinite(flux1))||(!isfinite(flux2))) {
+      fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
+    }
+
+    //    fprintf(stderr,"SOURCES: %d %g %g %g -> %d %d %g %g %d %g\n",i,source->peak->xf,source->peak->yf,source->psfMag,
+    //      nPix1,nPix2,flux1,flux2,fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i],logFluxDifferences->data.F32[i]);
+  }
+
+  // Given the differences in log-flux space, the normalization factor is just the exponential of the median difference
+  psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+  if (!psVectorStats(stats,logFluxDifferences,NULL,fitMask,0xff)) {
+    // This should complain.
+    normalization = 1.0;
+  }
+
+  normalization = pow(10,stats->robustMedian);
+  // fprintf(stderr,"NORM: %g+/-%g\n",stats->robustMedian,stats->robustStdev);
+  
+  psFree(stats);
+  psFree(logFluxDifferences);
+  psFree(fitMask);
+
+  // Apply normalization
+  if (normalization != 1.0) {
+    if ((conv1)&&((convolution_direction == 1))) {
+      psBinaryOp(imageC1,imageC1,"*",psScalarAlloc((float) normalization, PS_TYPE_F32));
+      psBinaryOp(varC1,varC1,"*",psScalarAlloc((float) PS_SQR(normalization), PS_TYPE_F32));
+    }
+    else if ((conv2)&&(convolution_direction == 2)) {
+      normalization = 1.0 / normalization; // Because we fit one way, but are using it in the other.
+      psBinaryOp(imageC2,imageC2,"*",psScalarAlloc((float) normalization, PS_TYPE_F32));
+      psBinaryOp(varC2,varC2,"*",psScalarAlloc((float) PS_SQR(normalization), PS_TYPE_F32));
+    }
+  }
+  
+
+  //
+  
   //
   // Actually add it to the headers

branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.h

@@ -16 +16 @@
                               psImageMaskType maskVal,
                               psImageMaskType maskBad,
-                              psImageMaskType maskPoor
+                              psImageMaskType maskPoor,
+                              float deconvolveThreshold
                               );
 

branches/eam_branches/ipp-20130711/psModules/src/objects/pmPSFtryFitPSF.c

@@ -75 +75 @@
             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y);
-            return false;
+            continue;
         }