Index: branches/eam_branches/ipp-20130711/PS-IPP-PStamp/lib/PS/IPP/PStamp/Job.pm =================================================================== --- branches/eam_branches/ipp-20130711/PS-IPP-PStamp/lib/PS/IPP/PStamp/Job.pm (revision 36074) +++ branches/eam_branches/ipp-20130711/PS-IPP-PStamp/lib/PS/IPP/PStamp/Job.pm (revision 36075) @@ -34,4 +34,5 @@ use POSIX; use Time::HiRes qw(gettimeofday); +use IO::Handle; my $dvo_verbose = 0; @@ -176,5 +177,6 @@ $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,4 +245,9 @@ 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,5 +329,5 @@ 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,4 +358,5 @@ my $fwhm_max = shift; my $selectedAstrom = shift; + my $selected_cam_id = shift; my $release_name = shift; my $survey = shift; @@ -602,5 +610,5 @@ # 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,4 +707,5 @@ $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,4 +740,118 @@ 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,5 +1069,5 @@ $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,5 +1093,5 @@ $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,5 +1289,11 @@ 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,4 +1431,5 @@ $image->{astrom} = $astromFile; $image->{cam_path_base} = $camRoot; + $image->{cam_id} = $camRun->{cam_id}; return 1; @@ -1649,4 +1779,5 @@ sub filterRuns { my $stage = shift; + my $drop_duplicate_runs = shift; my $need_magic = shift; my $inputs = shift; @@ -1694,5 +1825,5 @@ } - $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,5 +1833,5 @@ } - 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,4 +2004,195 @@ 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; +} + Index: branches/eam_branches/ipp-20130711/PS-IPP-PStamp/lib/PS/IPP/PStamp/RequestFile.pm =================================================================== --- branches/eam_branches/ipp-20130711/PS-IPP-PStamp/lib/PS/IPP/PStamp/RequestFile.pm (revision 36074) +++ branches/eam_branches/ipp-20130711/PS-IPP-PStamp/lib/PS/IPP/PStamp/RequestFile.pm (revision 36075) @@ -29,4 +29,6 @@ $PSTAMP_SELECT_BACKMDL $PSTAMP_SELECT_JPEG + $PSTAMP_SELECT_EXP + $PSTAMP_SELECT_NUM $PSTAMP_SELECT_UNCOMPRESSED $PSTAMP_SELECT_INVERSE @@ -71,6 +73,6 @@ 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,12 +99,13 @@ 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,5 +146,5 @@ PSTAMP_NO_OVERLAP PSTAMP_NOT_AUTHORIZED -PSTAMP_NOT_AUTHORIZED +PSTAMP_NO_VALID_PIXELS ); Index: branches/eam_branches/ipp-20130711/ppSub/src/ppSubMatchPSFs.c =================================================================== --- branches/eam_branches/ipp-20130711/ppSub/src/ppSubMatchPSFs.c (revision 36074) +++ branches/eam_branches/ipp-20130711/ppSub/src/ppSubMatchPSFs.c (revision 36075) @@ -57,5 +57,5 @@ if (!pmFPACopy(photFile->fpa, ro->parent->parent->parent)) { psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry"); - return false; + return NAN; } @@ -76,5 +76,5 @@ psWarning("Unable to determine PSF."); psFree(view); - return true; + return NAN; } psFree(view); Index: branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionMatch.c (revision 36074) +++ branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionMatch.c (revision 36075) @@ -540,5 +540,5 @@ // 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; } Index: branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.c =================================================================== --- branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.c (revision 36074) +++ branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.c (revision 36075) @@ -63,4 +63,39 @@ } +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,29 +106,32 @@ 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,5 +139,5 @@ if (conv1) { - conv1->covariance = psMemIncrRefCounter(ro1->covariance); + // conv1->covariance = psMemIncrRefCounter(ro1->covariance); if (!conv1->image) { conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); @@ -118,5 +156,5 @@ } if (conv2) { - conv2->covariance = psMemIncrRefCounter(ro2->covariance); + // if (!conv2->image) { conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); @@ -171,13 +209,17 @@ 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,114 +228,9 @@ 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,4 +266,126 @@ 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 Index: branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.h =================================================================== --- branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.h (revision 36074) +++ branches/eam_branches/ipp-20130711/psModules/src/imcombine/pmSubtractionSimple.h (revision 36075) @@ -16,5 +16,6 @@ psImageMaskType maskVal, psImageMaskType maskBad, - psImageMaskType maskPoor + psImageMaskType maskPoor, + float deconvolveThreshold ); Index: branches/eam_branches/ipp-20130711/psModules/src/objects/pmPSFtryFitPSF.c =================================================================== --- branches/eam_branches/ipp-20130711/psModules/src/objects/pmPSFtryFitPSF.c (revision 36074) +++ branches/eam_branches/ipp-20130711/psModules/src/objects/pmPSFtryFitPSF.c (revision 36075) @@ -75,5 +75,5 @@ 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; }