Changeset 36633 for trunk/psphot/src

Timestamp:

Apr 2, 2014, 10:50:14 AM (12 years ago)

Author:

bills

Message:

Merge with ipp-20140401 (The staticsky tag)

Location:

trunk

Files:

• . (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotFullForce.SourceStats.c (modified) (1 diff)
• psphot/src/psphotFullForceReadout.c (modified) (2 diffs)
• psphot/src/psphotFullForceSummary.c (modified) (2 diffs)
• psphot/src/psphotFullForceSummaryReadout.c (modified) (3 diffs)
• psphot/src/psphotGalaxyShape.c (modified) (3 diffs)
• psphot/src/psphotStackImageLoop.c (modified) (1 prop)

trunk/psphot/src/psphotFullForce.SourceStats.c

@@ -85 +85 @@
     psAssert (markVal, "missing MARK.PSPHOT");
 
-    psArray *sources = detections->allSources;
+    // psArray *sources = detections->allSources;
+    psArray *sources = detections->newSources ? detections->newSources : detections->allSources;
 
     // generate the array of sources, define the associated pixel

trunk/psphot/src/psphotFullForceReadout.c

@@ -52 +52 @@
     psphotLoadExtSources (config, view, filerule);
 
+#ifdef notmoved
+    // XXX: moved down below
     // merge the newly selected sources into the existing list (detections->allSources)
     // NOTE: merge OLD and NEW
     psphotMergeSources (config, view, filerule);
+#endif
 
     // construct sources and measure moments and other basic stats (saved on detections->allSources)
@@ -63 +66 @@
         return psphotReadoutCleanup (config, view, filerule);
     }
+    
+    // classify sources based on moments, brightness.  if a PSF model has been loaded, the PSF
+    // clump defined for it is used not measured (detections->newSources)
+    if (!psphotRoughClass (config, view, filerule)) { // pass 1
+        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
+        return psphotReadoutCleanup (config, view, filerule);
+    }
+
+    // merge the newly selected sources into the existing list (detections->allSources)
+    // NOTE: merge OLD and NEW
+    psphotMergeSources (config, view, filerule);
+
 
     // generate a psf model for any readouts which need one

trunk/psphot/src/psphotFullForceSummary.c

@@ -78 +78 @@
         usage();
     }
+    if ((N = psArgumentGet(argc, argv, "-cff"))) {
+        if (argc <= N+1) {
+          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
+          usage();
+        }
+        psArgumentRemove(N, &argc, argv);
+
+        // We read CFF file as table not as sources. All we want are the nominal parameters input to psphotFullForce
+        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "CFF_FILE", 0, "CFF file with nominal parameters", argv[N]);
+        psArgumentRemove(N, &argc, argv);
+    } else {
+        psErrorStackPrint(stderr, "-cff must be supplied.");
+        usage();
+    }
  
     if (argc < 2) {
@@ -148 +162 @@
     }
     output->save = true;
+
+    // We read cff file as a table not as sources
+    psString cffName = psMetadataLookupStr(NULL, config->arguments, "CFF_FILE");
+    if (!cffName) {
+        psError(PSPHOT_ERR_CONFIG, true, "CFF_FILE is missing from arguments");
+        return false;
+    }
+    psString resolvedName = pmConfigConvertFilename(cffName, config, false, false);
+    if (!resolvedName) {
+        psError(PSPHOT_ERR_CONFIG, false, "failed to resolve CFF_FILE %s", cffName);
+        return false;
+    }
+    psFits *fits = psFitsOpen(resolvedName, "r");
+    if (!fits) {
+        psError(PSPHOT_ERR_CONFIG, false, "failed to open fits CFF_FILE %s", resolvedName);
+        return false;
+    }
+    psArray *inTable = psFitsReadTable(fits);
+    if (!inTable) {
+        psError(PSPHOT_ERR_CONFIG, false, "failed to read cff fits table from CFF_FILE %s", resolvedName);
+        return false;
+    }
+    psFitsClose(fits);
+
+    // Convert to an array indexed by ID
+    psArray *sortedTable = psArrayAlloc(4*inTable->n);
+    for (int i=0; i<inTable->n; i++) {
+        psMetadata *row = inTable->data[i];
+        psS32 ID = psMetadataLookupS32(&status, row, "ID");
+        if (ID >= sortedTable->n) {
+            sortedTable = psArrayRealloc(sortedTable, 2*ID);
+        }
+        if (sortedTable->data[ID]) {
+            psError(PSPHOT_ERR_CONFIG, true, "Duplicate row with ID %d", ID);
+            return false;
+        }
+        sortedTable->data[ID] = psMemIncrRefCounter(row);
+    }
+    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLE", PS_META_REPLACE, "cff table", sortedTable);
+    psFree(inTable);
+    psFree(sortedTable);
 
     return true;

trunk/psphot/src/psphotFullForceSummaryReadout.c

@@ -14 +14 @@
     psArray *zeroPt;
     psArray *exptime;
+    psArray *cffTable;
 } galaxyShapeOptions;
 
@@ -69 +70 @@
             if (!setOptions(&options, readout, recipe, true)) {
                 psError (PS_ERR_UNKNOWN, false, "problem determining galaxy shape options.");
+                return false;
+            }
+            options.cffTable = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLE");
+            if (!options.cffTable) {
+                psError (PS_ERR_UNKNOWN, true, "Cannot find cff table in arguments.");
                 return false;
             }
@@ -309 +315 @@
                 // copy the best fit params to the model
                 psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type);
-                axes.major *= options->rMajor->data.F32[min_j];
-                axes.minor *= options->rMinor->data.F32[min_j];
+                if (outSrc->seq >= options->cffTable->n) {
+                    psError(PS_ERR_PROGRAMMING, true, "object sequence number %d is larger than cff table length %ld",
+                        outSrc->seq, options->cffTable->n);
+                    return NULL;
+                }
+                psMetadata *row = options->cffTable->data[outSrc->seq];
+                if (!row) {
+                    psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d", outSrc->seq);
+                    psFree(outSrc);
+                    return NULL;
+                }
+                bool mdok;
+                psF32 rMajor = psMetadataLookupF32(&mdok, row, "R_MAJOR");
+                if (!mdok) {
+                    psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d", outSrc->seq);
+                    psFree(outSrc);
+                    return NULL;
+                }
+                psF32 rMinor = psMetadataLookupF32(&mdok, row, "R_MINOR");
+                if (!mdok) {
+                    psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d", outSrc->seq);
+                    psFree(outSrc);
+                    return NULL;
+                }
+                axes.major = rMajor * options->rMajor->data.F32[min_j];
+                axes.minor = rMinor * options->rMinor->data.F32[min_j];
                 pmPSF_AxesToModel (outSrc->modelEXT->params->data.F32, axes, outSrc->modelEXT->type);
                 outSrc->modelEXT->chisq = minChisq; 

trunk/psphot/src/psphotGalaxyShape.c

@@ -205 +205 @@
         if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
         if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+        if (!source->modelPSF) continue;
 
         // psphotSetRadiusMomentsExact sets the radius based on Mrf
@@ -274 +275 @@
 
             int i = source->galaxyFits->chisq->n - 1;
-            float thisChisq = source->galaxyFits->chisq->data.F32[i];
-            if (isfinite(thisChisq) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
-                chisqBest = thisChisq;
-                fRmajorBest = fRmajor;
-                fRminorBest = fRminor;
+            float flux = source->galaxyFits->Flux->data.F32[i];
+            if (isfinite(flux)) {
+                float thisChisq = source->galaxyFits->chisq->data.F32[i];
+                if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
+                    chisqBest = thisChisq;
+                    fRmajorBest = fRmajor;
+                    fRminorBest = fRminor;
+                }
             }
+            // reset I0 to avoid potential problems on the next iteration
+            PAR[PM_PAR_I0] = 1.0;
         }
     }
 
+#define SAVE_BEST_MODEL
 #ifdef SAVE_BEST_MODEL
     // Save model with smallest chisq
@@ -301 +308 @@
         psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, false);
 
-        // Replace modelEXT with this model
-        // XXX: only do this if the model is good
-        psFree (source->modelEXT);
-
-        source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
-        source->type = PM_SOURCE_TYPE_EXTENDED;
-        source->mode |= PM_SOURCE_MODE_EXTMODEL;
-        source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
-
-        // cache the model flux
-        pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
+        // Replace modelEXT with this model, if the model is good
+        if (isfinite(PAR[PM_PAR_I0])) {
+            psFree (source->modelEXT);
+
+            source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
+            source->type = PM_SOURCE_TYPE_EXTENDED;
+            source->mode |= PM_SOURCE_MODE_EXTMODEL;
+            source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
+
+            // cache the model flux
+            pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
+        }
     }