Changeset 35941

Timestamp:

Aug 13, 2013, 4:01:17 PM (13 years ago)

Author:

watersc1

Message:

Merge from trunk of stacking changes. 1) Option to set target PSF to be marginally larger than the largest accepted input PSF. 2) Fix inconsistencies in the masking for simple convolution.

Location:

tags/ipp-20130712

Files:

• ippconfig/gpc1/ppStack.config (modified) (1 diff, 1 prop)
• ippconfig/recipes/ppStack.config (modified) (1 diff, 1 prop)
• ppStack/src/ppStackPrepare.c (modified) (2 diffs)
• psModules/src/imcombine/pmSubtractionMatch.c (modified) (1 diff)
• psModules/src/imcombine/pmSubtractionSimple.c (modified) (10 diffs, 1 prop)
• psModules/src/imcombine/pmSubtractionSimple.h (modified) (1 diff)

tags/ipp-20130712/ippconfig/gpc1/ppStack.config

@@ -63 +63 @@
     PSF.INPUT.THRESH        F32   10.0
     PSF.INPUT.ASYMMETRY     F32   0.2
+    PSF.TARGET.AS.MAX       BOOL  TRUE
    MATCH.REJ        F32    4.0
    SAFE          BOOL   F

tags/ipp-20130712/ippconfig/recipes/ppStack.config

@@ -95 +95 @@
 PSF.INPUT.THRESH        F32  NAN        # Set minimum limit below which we do not exclude an input (defaults to 0.0)
 PSF.INPUT.ASYMMETRY     F32  NAN        # Set difference in mixture model populations to consider equal.
+PSF.TARGET.AS.MAX       BOOL F          # Set the target PSF FWHM as the maximum of accepted input FWHM values.
+PSF.TARGET.AS.MAX.EPSILON F32 0.1       # Amount to set the target PSF FWHM larger than the maximum input. Target = eps + max(input)
+
 
 TEMP.IMAGE      STR     conv.im.fits    # Suffix for temporary convolved images

tags/ipp-20130712/ppStack/src/ppStackPrepare.c

@@ -257 +257 @@
 
     bool mdok = false;
-    bool  simpleClip = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.SIMPLE");
+    bool  simpleClip = psMetadataLookupBool(&mdok, recipe, "PSF.INPUT.CLIP.SIMPLE");
     float maxFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.MAX"); // max allowed input fwhm
     float clipFWHMnSig = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.NSIGMA"); // sigma clipping of inputs
     float threshFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.THRESH"); // FWHM size that we refuse to clip below
     float asymmetryFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.ASYMMETRY"); // max bimodal asymmetry
+
     
     psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message
@@ -477 +478 @@
                          "Target PSF for stack", options->psf);
         options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * numCols, 0.5 * numRows); // FWHM for target
-        psLogMsg("ppStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing);
-
+
+
+        bool psfTargetAsMax = psMetadataLookupBool(&mdok, recipe, "PSF.TARGET.AS.MAX");
+        if (psfTargetAsMax) { // Should we use the largest input as the target?
+          float psfTargetEpsilon = psMetadataLookupF32(&mdok,recipe,"PSF.TARGET.AS.MAX.EPSILON");
+          if (!mdok) { psfTargetEpsilon = 0.0; }
+          options->targetSeeing = 0.0;
+          for (int i = 0; i < num; i++) {
+            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
+            options->targetSeeing = PS_MAX(options->targetSeeing,options->inputSeeing->data.F32[i]);
+          }
+          psLogMsg("ppStack", PS_LOG_INFO, "Using MAX accepted input FWHM as target (max: %f epsilon: %f target %f)\n",
+                   options->targetSeeing,psfTargetEpsilon,options->targetSeeing + psfTargetEpsilon);
+          options->targetSeeing = options->targetSeeing + psfTargetEpsilon;
+        }
+        psLogMsg("ppStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing);    
+        
         pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.TARGET.PSF"); // Output chip
         psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN,

tags/ipp-20130712/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;
       }

tags/ipp-20130712/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
@@ -90 +106 @@
                               psImageMaskType maskVal,
                               psImageMaskType maskBad,
-                              psImageMaskType maskPoor
+                              psImageMaskType maskPoor,
+                              float deconvolveThreshold
                               ) {
   //
@@ -115 +132 @@
   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;
@@ -190 +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;
   }
   
@@ -205 +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?
-
+  int maskPoorVal = 16384; // Another value that should be found elsewhere.
   //
   // Make a fake pmSubtractionKernels element so we can add it appropriately.
@@ -249 +272 @@
   // 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);
-    conv1->covariance = psImageCovarianceCalculate(kernel,ro1->covariance);
+    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);
@@ -260 +293 @@
       conv2->covariance = psMemIncrRefCounter(ro2->covariance);
     }
-    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);
-    conv2->covariance = psImageCovarianceCalculate(kernel,ro2->covariance);
+    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);
@@ -275 +316 @@
       conv1->covariance = psMemIncrRefCounter(ro1->covariance);
     }
-    pmSubtractionBorder(imageC2,varC2,maskC2,maskBox,maskBlank);
-    pmSubtractionMaskApply(imageC2,varC2,maskC2,PM_SUBTRACTION_MODE_2);
   }    
 
@@ -294 +333 @@
     float flux1,flux2;
 
-    if (convolution_direction == 1) {
+    if (conv1) {
       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) {
+    }
+    else {
+      simple_do_boxphot(&nPix1,&flux1,source,image1,mask1,maskBad,photRadius);
+    }
+
+    if (conv2) {
       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);
-      }
-    }
+    }
+    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;
@@ -321 +356 @@
     //    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]);
-    
   }
 

tags/ipp-20130712/psModules/src/imcombine/pmSubtractionSimple.h

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