Changeset 26321 for branches/pap/psModules/src/imcombine/pmSubtractionStamps.c

Timestamp:

Dec 2, 2009, 12:04:51 PM (16 years ago)

Author:

Paul Price

Message:

Merging Gene's work on adding a window function to the stamps. I'd been thinking about something like this for a while, so glad someone's done it. It's got visualisation turned on though.

Location:

branches/pap/psModules/src/imcombine

Files:

• . (modified) (1 prop)
• pmSubtractionStamps.c (modified) (4 diffs)

branches/pap/psModules/src/imcombine/pmSubtractionStamps.c

@@ -46 +46 @@
     psFree(list->y);
     psFree(list->flux);
+    psFree(list->window);
 }
 
@@ -220 +221 @@
     list->y = NULL;
     list->flux = NULL;
+    list->window = NULL;
     list->footprint = footprint;
     list->sysErr = sysErr;
@@ -239 +241 @@
     out->y = NULL;
     out->flux = NULL;
+    out->window = psMemIncrRefCounter(in->window);
     out->footprint = in->footprint;
 
@@ -606 +609 @@
 }
 
+
+bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
+    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
+
+    int size = kernelSize + stamps->footprint; // Size of postage stamps
+
+    psFree (stamps->window);
+    stamps->window = psKernelAlloc(-size, size, -size, size);
+    psImageInit(stamps->window->image, 0.0);
+
+    // generate normalizations for each stamp
+    psVector *norm1 = psVectorAlloc(stamps->num, PS_TYPE_F32);
+    psVector *norm2 = psVectorAlloc(stamps->num, PS_TYPE_F32);
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (!stamp) continue;
+        if (!stamp->image1) continue;
+        if (!stamp->image2) continue;
+
+        float sum1 = 0.0;
+        float sum2 = 0.0;
+        for (int y = -size; y <= size; y++) {
+            for (int x = -size; x <= size; x++) {
+                sum1 += stamp->image1->kernel[y][x];
+                sum2 += stamp->image2->kernel[y][x];
+            }
+        }
+        norm1->data.F32[i] = sum1;
+        norm2->data.F32[i] = sum2;
+    }
+
+    // storage vector for flux data
+    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psVector *flux = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
+
+    // generate the window pixels
+    for (int y = -size; y <= size; y++) {
+        for (int x = -size; x <= size; x++) {
+
+            flux->n = 0;
+            for (int i = 0; i < stamps->num; i++) {
+                pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+                if (!stamp) continue;
+                if (!stamp->image1) continue;
+                if (!stamp->image2) continue;
+                
+                psVectorAppend (flux, stamp->image1->kernel[y][x] / norm1->data.F32[i]);
+                psVectorAppend (flux, stamp->image2->kernel[y][x] / norm2->data.F32[i]);
+            }
+
+            psStatsInit (stats);
+            if (!psVectorStats (stats, flux, NULL, NULL, 0)) {
+                psAbort ("failed to generate stats");
+            }
+            stamps->window->kernel[y][x] = stats->robustMedian;
+        }
+    }
+
+    // XXX TESTING: 
+    {
+        psFits *fits = psFitsOpen ("window.fits", "w");
+        psFitsWriteImage (fits, NULL, stamps->window->image, 0, NULL);
+        psFitsClose (fits);
+    }
+
+    psFree (stats);
+    psFree (flux);
+    psFree (norm1);
+    psFree (norm2);
+    return true;
+}
 
 bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,