Changeset 42830

Timestamp:

May 8, 2025, 4:53:16 PM (12 months ago)

Author:

eugene

Message:

add exp Weights for pmStackCombineByPercentile (UNIONS stacks)

File:

• trunk/psModules/src/imcombine/pmStack.c (modified) (4 diffs)

trunk/psModules/src/imcombine/pmStack.c

@@ -1527 +1527 @@
             pixelVariances->data.F32[A] = pixelVariances->data.F32[B];  \
             pixelVariances->data.F32[B] = tempVar;                      \
+            psF32 tempWgt = expWeight->data.F32[A];             \
+            expWeight->data.F32[A] = expWeight->data.F32[B];    \
+            expWeight->data.F32[B] = tempWgt;                                   \
             if (expTime) {                                              \
                 psF32 tempExp = expTime->data.F32[A];                   \
@@ -1589 +1592 @@
     psVector *pixelVariances = psVectorAlloc(stackData->n, PS_TYPE_F32);
 
+    // storage vector for exposures weights (for pixels with few inputs)
+    psVector *expWeight = psVectorAlloc(stackData->n, PS_TYPE_F32);
+
     // if we are asking for the exptime maps, generate a storage vector expTime
     psVector *expTime        = expmaps && expmaps->image ? psVectorAlloc(stackData->n, PS_TYPE_F32) : NULL;
@@ -1640 +1646 @@
                 // accumulate exposure times if required
                 if (expTime)  { expTime->data.F32[nGood] = data->exp; }
+                expWeight->data.F32[nGood] = data->exp;
                 nGood ++;
             }
@@ -1729 +1736 @@
             }
 
+            // XXX should I use the expWeight regardless of the number of inputs??
+
             float sum = 0.0;
             float varSum = 0.0;
+            float wtSum = 0.0;
             for (int n = Ns; n < Ne; n++) {
-                sum += pixelData->data.F32[n];
-                varSum += pixelVariances->data.F32[n];
-            }
-            float mean = sum / (float) Npt;
-            float varValue = varSum / (float) (nGoodClip*nGoodClip);
+                psF32 wt = expWeight->data.F32[n];
+                psF32 wt2 = wt*wt;
+                sum += pixelData->data.F32[n]*wt;
+                varSum += pixelVariances->data.F32[n]*wt2;
+                wtSum += wt;
+            }
+            // float mean = sum / (float) Npt;
+            // float varValue = varSum / (float) (nGoodClip*nGoodClip);
+            float mean = sum / wtSum;
+            float varValue = varSum / PS_SQR(wtSum);
 
             // alternative: calculate the stdev of the pixel values