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

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)
• psModules (modified) (1 prop)
• psModules/src/imcombine/pmSubtractionSimple.c (modified) (7 diffs)

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