Changeset 27517 for branches/eam_branches/20100225/ppStack/src/ppStackSources.c

Timestamp:

Mar 29, 2010, 3:55:49 PM (16 years ago)

Author:

eugene

Message:

update merges from trunk

Location:

branches/eam_branches/20100225

Files:

• . (modified) (1 prop)
• ppStack/src/ppStackSources.c (modified) (13 diffs)

branches/eam_branches/20100225/ppStack/src/ppStackSources.c

@@ -8 +8 @@
 
 //#define TESTING                         // Enable debugging output
+
+//#define ASTROMETRY                    // Correct astrometry?
 
 #ifdef TESTING
@@ -62 +64 @@
 
     if (!options->matchZPs && !options->photometry) {
-        int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
-        options->norm = psVectorAlloc(num, PS_TYPE_F32);
-        psVectorInit (options->norm, 0.0);
-
-        // XXX do I need to set this?
-        // options->sumExposure = sumExpTime;
-
+        options->norm = psVectorAlloc(options->num, PS_TYPE_F32);
+        psVectorInit(options->norm, 0.0);
         return true;
     }
@@ -90 +87 @@
             }
             source->psfMag += 1.0;
+#ifdef ASTROMETRY
+            if (source->modelPSF) {
+                source->modelPSF->params->data.F32[PM_PAR_XPOS] += 1.0;
+                source->modelPSF->params->data.F32[PM_PAR_YPOS] += 1.0;
+            }
+            if (source->peak) {
+                source->peak->xf += 1.0;
+                source->peak->yf += 1.0;
+            }
+#endif
         }
     }
@@ -119 +126 @@
         return false;
     }
-
-    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+    psMetadata *zpTargetMenu = psMetadataLookupMetadata(NULL, recipe, "ZP.TARGET"); // Target zero point terms
+    if (!zpTargetMenu) {
+        psError(PPSTACK_ERR_CONFIG, false, "Unable to find ZP.TARGET in recipe.");
+        return false;
+    }
+
+    int num = options->num;             // Number of inputs
     psAssert(num == sourceLists->n, "Wrong number of source lists: %ld\n", sourceLists->n);
 
-    psVector *zp = psVectorAlloc(num, PS_TYPE_F32); // Zero points for each image
+    psVector *zp = psVectorAlloc(num, PS_TYPE_F32); // Relative zero points for each image
+    psVector *zpExp = psVectorAlloc(num, PS_TYPE_F32); // Measured zero points for each image (maybe)
+    int zpExpNum = 0;                                  // Number of measured zero points
     const char *filter = NULL;          // Filter name
     float airmassTerm = NAN;            // Airmass term
-    float sumExpTime = 0.0;             // Sum of the exposure time
+    float zpTarget = NAN;               // Target zero point
+    int numGoodImages = 0;              // Number of good images
     for (int i = 0; i < num; i++) {
+        psArray *sources = sourceLists->data[i]; // Source list
+        if (!sources || sources->n == 0) {
+            psLogMsg("ppStack", PS_LOG_WARN, "Image %d has no sources for transparency measurement.", i);
+            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_CAL;
+            zp->data.F32[i] = NAN;
+            continue;
+        }
+        numGoodImages++;
+
         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
-        pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
 
 #if defined(TESTING) && 0
@@ -146 +169 @@
 #endif
 
-
-        float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time
+        float exptime = options->exposures->data.F32[i]; // Exposure time
         float airmass = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.AIRMASS"); // Airmass
         const char *expFilter = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.FILTER"); // Filter name
+        zpExp->data.F32[i] = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.ZP"); // Exposure zero point
+        psLogMsg("ppStack", PS_LOG_INFO,
+                 "Image %d: %.2f sec exposure in %s at airmass %.2f with zero point %.2f",
+                 i, exptime, expFilter, airmass, zpExp->data.F32[i]);
         if (!isfinite(exptime) || exptime == 0 || !isfinite(airmass) || airmass == 0 ||
             !expFilter || strlen(expFilter) == 0) {
@@ -158 +184 @@
             return false;
         }
+        if (isfinite(zpExp->data.F32[i])) {
+            zpExp->data.F32[i] += 2.5 * log10(exptime);
+            zpExpNum++;
+        }
 
         if (!filter) {
             filter = expFilter;
-            airmassTerm = psMetadataLookupF32(NULL, airmassZP, filter);
-            if (!isfinite(airmassTerm)) {
+            bool mdok;
+            airmassTerm = psMetadataLookupF32(&mdok, airmassZP, filter);
+            if (!mdok || !isfinite(airmassTerm)) {
                 psError(PPSTACK_ERR_CONFIG, false,
                         "Unable to find airmass term (ZP.AIRMASS) for filter %s", filter);
@@ -168 +199 @@
                 return false;
             }
+            zpTarget = psMetadataLookupF32(&mdok, zpTargetMenu, filter);
+            if (!mdok || !isfinite(zpTarget)) {
+                psError(PPSTACK_ERR_CONFIG, false,
+                        "Unable to find target zero point (ZP.TARGET) for filter %s", filter);
+                psFree(zp);
+                return false;
+            }
         } else if (strcmp(filter, expFilter) != 0) {
             psError(PPSTACK_ERR_CONFIG, false, "Filters don't match: %s vs %s", filter, expFilter);
@@ -174 +212 @@
         }
 
-        zp->data.F32[i] = airmassTerm * airmass - 2.5 * log10(exptime);
-        sumExpTime += exptime;
-    }
-
-    options->sumExposure = sumExpTime;
+        zp->data.F32[i] = airmassTerm * airmass + 2.5 * log10(exptime);
+    }
+
+    if (numGoodImages == 0) {
+        psLogMsg("ppStack", PS_LOG_WARN, "No images with sources to measure transparency.");
+        options->quality = PPSTACK_ERR_REJECTED;
+        psFree(zp);
+        psFree(zpExp);
+        return true;
+    }
+    if (numGoodImages == 1) {
+        psArray *sources = NULL;        // Sources
+        for (int i = 0; i < num && !sources; i++) {
+            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+                continue;
+            }
+            sources = sourceLists->data[i];
+        }
+        options->sources = psMemIncrRefCounter(sources);
+        psLogMsg("ppStack", PS_LOG_WARN, "Single image with sources --- no need to match transparency.");
+        psFree(zp);
+        psFree(zpExp);
+        return true;
+    }
+
+    if (zpExpNum == numGoodImages) {
+        for (int i = 0; i < num; i++) {
+            zp->data.F32[i] = zpExp->data.F32[i];
+        }
+    }
 
     psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches
@@ -212 +275 @@
             options->norm = psVectorAlloc(num, PS_TYPE_F32);
             for (int i = 0; i < num; i++) {
-                if (!isfinite(trans->data.F32[i])) {
+                if (inputMask->data.U8[i] || !isfinite(trans->data.F32[i])) {
                     continue;
                 }
                 psArray *sources = sourceLists->data[i]; // Sources of interest
-                float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];
-                options->norm->data.F32[i] = magCorr;
-                psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n",
-                         i, magCorr);
+                float magCorr = zp->data.F32[i] + trans->data.F32[i] - 2.5*log10(options->sumExposure);
+                if (zpExpNum == numGoodImages) {
+                    // Using measured zero points, so attempt to set target zero point
+                    magCorr -= zpTarget;
+                }
+                options->norm->data.F32[i] = -magCorr;
+                psLogMsg("ppStack", PS_LOG_INFO,
+                         "Applying scale correction to image %d: %f mag (%f)\n",
+                         i, magCorr, trans->data.F32[i]);
 
                 for (int j = 0; j < sources->n; j++) {
@@ -226 +294 @@
                         continue;
                     }
-                    source->psfMag += magCorr;
-                }
-            }
-        }
+                    source->psfMag -= magCorr;
+                }
+            }
+        }
+
+        if (zpExpNum == numGoodImages) {
+            // Producing image with target zero point
+            options->zp = zpTarget;
+        } else {
+            options->zp = NAN;
+        }
+
 
 #ifdef TESTING
@@ -239 +315 @@
         // Double check: all transparencies should be zero
         {
-            psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
+            psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches
             if (!matches) {
                 psError(PPSTACK_ERR_DATA, false, "Unable to match sources");
@@ -245 +321 @@
                 return false;
             }
-            psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
-                                                   transThresh, starRej, starSys);
+            psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
+                                                   iter2, starRej2, starSys2, starLimit,
+                                                   transIter, transRej, transThresh); // Transparencies
             for (int i = 0; i < num; i++) {
                 fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
             }
             psFree(trans);
-        }
-#endif
-    }
-
-
-#if 0
+            psFree(matches);
+        }
+#endif
+    }
+
+    psFree(zp);
+    psFree(zpExp);
+
+#ifdef ASTROMETRY
     // Position offsets
     {
@@ -266 +346 @@
         }
         for (int i = 0; i < num; i++) {
-            psVector *offset = offsets->data[i];
-            fprintf(stderr, "Offset %d: %f,%f\n", i, offset->data.F32[0], offset->data.F32[1]);
+            if (options->inputMask->data.U8[i]) {
+                continue;
+            }
+            psArray *sources = sourceLists->data[i]; // Sources of interest
+            psVector *offset = offsets->data[i];                      // Offsets for image
+            float dx = offset->data.F32[0], dy = offset->data.F32[1]; // Offsets to apply
+            if (!isfinite(dx) || !isfinite(dy)) {
+                continue;
+            }
+            psLogMsg("ppStack", PS_LOG_INFO, "Applying astrometric correction to image %d: %f,%f\n",
+                     i, dx, dy);
+            for (int j = 0; j < sources->n; j++) {
+                pmSource *source = sources->data[j]; // Source of interest
+                if (!source) {
+                    continue;
+                }
+                if (source->modelPSF) {
+                    source->modelPSF->params->data.F32[PM_PAR_XPOS] -= dx;
+                    source->modelPSF->params->data.F32[PM_PAR_YPOS] -= dy;
+                }
+                if (source->peak) {
+                    source->peak->xf -= dx;
+                    source->peak->yf -= dy;
+                }
+            }
         }
         psFree(offsets);
     }
 #endif
+
+#if (defined TESTING && defined ASTROMETRY)
+        // Double check: all offsets should be zero
+        {
+            psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches
+            if (!matches) {
+                psError(PPSTACK_ERR_DATA, false, "Unable to match sources");
+                psFree(zp);
+                return false;
+            }
+            psArray *offsets = pmSourceMatchRelastro(matches, num, tol, iter1, starRej1,
+                                                     iter2, starRej2, starLimit); // Shifts for each image
+            for (int i = 0; i < num; i++) {
+                psVector *offset = offsets->data[i]; // Offsets for image
+                fprintf(stderr, "Offset of image %d: %f,%f\n", i, offset->data.F32[0], offset->data.F32[1]);
+            }
+            psFree(offsets);
+            psFree(matches);
+        }
+#endif
+
 
     if (options->photometry) {