Changeset 29599 for trunk/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Oct 28, 2010, 5:51:55 PM (16 years ago)

Author:

eugene

Message:

compile-time option to use the windowed fft

File:

• trunk/psModules/src/imcombine/pmSubtraction.c (modified) (8 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -27 +27 @@
 #include "pmSubtraction.h"
 
+# define FFT_WINDOW 0
 //#define TESTING
 
@@ -96 +97 @@
     return;
 }
-
 
 // Generate an image of the solved kernel
@@ -187 +187 @@
     }
 
+    /* we have three possible implementations for handling flux conservation and photometric scaling.
+       
+       1) the original implementation of Alard-Lupton subtracted the 0th kernel gaussian from
+          the rest of the kernels.  the non-zero integral of the 0th kernel allows for some
+          flux re-scaling if needed.
+
+       2) the original ppSub implementation subtracted a delta function from the kernels.  this
+          seems to make the solution sensitive to noise terms.
+
+       3) accept the measured flux normalization and make all kernels have zero integral, but
+          use a gaussian to zero the flux.  
+
+          *** the operation below is only needed if we use option (2) ***
+    */
+
     if (normalise) {
         // Put in the normalisation component
         kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels));
     }
-
     return kernel;
 }
@@ -243 +257 @@
     psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Subimage mask
 
+# if (FFT_WINDOW)
+    psImage *convolved = psImageConvolveFFTwithWindow(NULL, subImage, subMask, maskVal, kernel); // Convolution
+# else
     psImage *convolved = psImageConvolveFFT(NULL, subImage, subMask, maskVal, kernel); // Convolution
+# endif
 
     psFree(subImage);
@@ -287 +305 @@
 
     // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once
+# if (FFT_WINDOW)
+    psImage *convVariance = psImageConvolveFFTwithWindow(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance
+    psImage *convKE = subKE ? psImageConvolveFFTwithWindow(NULL, subKE, subMask, maskVal, kernel) : NULL; // Conv KE
+# else
     psImage *convVariance = psImageConvolveFFT(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance
     psImage *convKE = subKE ? psImageConvolveFFT(NULL, subKE, subMask, maskVal, kernel) : NULL; // Conv KE
+# endif
 
     psFree(subVariance);
@@ -539 +562 @@
     PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
 
+# if (FFT_WINDOW)
+    psImage *conv = psImageConvolveFFTwithWindow(NULL, image->image, NULL, 0, kernel); // Convolved image
+# else
     psImage *conv = psImageConvolveFFT(NULL, image->image, NULL, 0, kernel); // Convolved image
+# endif
 
     // note: do not attempt to renormalize kernels here: cannot have different stars with
@@ -978 +1005 @@
 }
 
+// generate an image of the convolution kernel realized at the given coordinate
+// if 'wantDual' is set, solution2 is supplied
 psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)
 {
@@ -1285 +1314 @@
 
             psRegion *subRegion = psRegionAlloc(i, xSubMax, j, ySubMax); // Bounds of subtraction
-            if (threaded) {
+            // for a TEST, do not run threaded for testing
+            // if (false && threaded) {
+            if (threaded) {
                 psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CONVOLVE");
                 psArray *args = job->args;