Changeset 36680 for branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmSubtractionSimple.c

Timestamp:

Apr 21, 2014, 5:42:34 AM (12 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ps2-tc3-20130727

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/imcombine/pmSubtractionSimple.c (modified) (11 diffs)

branches/eam_branches/ps2-tc3-20130727/psModules/src/imcombine/pmSubtractionSimple.c

@@ -34 +34 @@
 #include "pmSourceSatstar.h"
 
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 
@@ -63 +64 @@
 }  
 
+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 +107 @@
                               psImageMaskType maskVal,
                               psImageMaskType maskBad,
-                              psImageMaskType maskPoor
+                              psImageMaskType maskPoor,
+                              float deconvolveThreshold
                               ) {
   //
@@ -115 +133 @@
   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 +210 @@
   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 +229 @@
   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 +273 @@
   // 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 +294 @@
       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 +317 @@
       conv1->covariance = psMemIncrRefCounter(ro1->covariance);
     }
-    pmSubtractionBorder(imageC2,varC2,maskC2,maskBox,maskBlank);
-    pmSubtractionMaskApply(imageC2,varC2,maskC2,PM_SUBTRACTION_MODE_2);
   }    
 
@@ -294 +334 @@
     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 +357 @@
     //    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]);
-    
   }