Changeset 42507

Timestamp:

Aug 16, 2023, 10:31:04 AM (3 years ago)

Author:

eugene

Message:

add test code to examine the PSF residual inputs

Location:

branches/eam_branches/ipp-20230313/psphot/src

Files:

• psphotArguments.c (modified) (2 diffs)
• psphotFullForceArguments.c (modified) (2 diffs)
• psphotFullForceReadout.c (modified) (1 diff)
• psphotMakeResiduals.c (modified) (12 diffs)
• psphotOutput.c (modified) (1 diff)

branches/eam_branches/ipp-20230313/psphot/src/psphotArguments.c

@@ -97 +97 @@
 }
 
+void psphotSaveConfig (pmConfig *config);
+
 pmConfig *psphotArguments(int argc, char **argv) {
 
@@ -244 +246 @@
     psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "", argv[1]);
 
+    // save a copy for deep retrieval
+    psphotSaveConfig (config);
+
     psTrace("psphot", 1, "Done with psphotArguments...\n");
     return (config);

branches/eam_branches/ipp-20230313/psphot/src/psphotFullForceArguments.c

@@ -3 +3 @@
 static void usage(const char *program, psMetadata *arg, pmConfig *config, int exitCode);
 static void writeHelpInfo(const char* program, pmConfig* config, FILE* ofile);
+void psphotSaveConfig (pmConfig *config);
 
 pmConfig *psphotFullForceArguments(int argc, char **argv) {
@@ -111 +112 @@
     // output position is fixed
     psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "", argv[1]);
+
+    psphotSaveConfig (config);
 
     psTrace("psphot", 1, "Done with psphotFullForceArguments...\n");

branches/eam_branches/ipp-20230313/psphot/src/psphotFullForceReadout.c

@@ -106 +106 @@
     psphotKronFlux (config, view, filerule);
 
+    // Option to do the non-linear fitting for the brighter sources
+    if (1) {
+      // identify CRs and extended sources (only unmeasured sources are measured)
+      psphotSourceSize (config, view, filerule, true); // pass 1 (detections->allSources)
+
+      // non-linear PSF and EXT fit to brighter sources
+      // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
+      // XXX: can leave faulted job in done queue
+      psphotBlendFit (config, view, filerule); // pass 1 (detections->allSources)
+
+      // replace all sources
+      psphotReplaceAllSources (config, view, filerule, false); // pass 1 (detections->allSources)
+
+      // linear fit to include all sources (subtract again)
+      // NOTE : apply to ALL sources (extended + psf)
+      // do NOT skip negative-flux sources (last arg = false)
+      psphotFitSourcesLinear (config, view, filerule, true, false); // pass 2 (detections->allSources)
+    }
+
     psphotChipParams (config, view, filerule);
 

branches/eam_branches/ipp-20230313/psphot/src/psphotMakeResiduals.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+pmConfig *psphotGetConfig ();
+bool _psphotSaveCube (char *basename, char *extname, psArray *cube);
 
 # define RESIDUAL_SOFTENING 0.005 
@@ -93 +95 @@
     psVector *yC = psVectorAllocEmpty (100, PS_TYPE_F32);
 
+# define TESTRESID 1
+# if (TESTRESID)
+    psArray *sourceRawCube = psArrayAllocEmpty (1);
+    psArray *sourceSigCube = psArrayAllocEmpty (1);
+    psArray *sourceVarCube = psArrayAllocEmpty (1);
+    psArray *sourceMskCube = psArrayAllocEmpty (1);
+# endif
+
     // build (DATA - MODEL) [an image] for each psf star
     psArray *input = psArrayAllocEmpty (100);
@@ -108 +118 @@
         psImage *mask     = psImageCopy (NULL, source->maskView ? source->maskView : source->maskObj,  PS_TYPE_IMAGE_MASK);
         psImage *variance = psImageCopy (NULL, source->variance ? source->variance : source->pixels,   PS_TYPE_F32);
+
+# if (TESTRESID)
+        psImage *tmpImage = psImageCopy(NULL, image, PS_TYPE_F32);
+        psArrayAdd (sourceRawCube, 1, tmpImage);
+        psFree (tmpImage);
+# endif
+
         pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal);
 
@@ -114 +131 @@
         psBinaryOp (image, image, "/", psScalarAlloc(Io, PS_TYPE_F32));
         psBinaryOp (variance, variance, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
+
+# if (TESTRESID)
+        psArrayAdd (sourceSigCube, 1, image);
+        psArrayAdd (sourceVarCube, 1, variance);
+        psArrayAdd (sourceMskCube, 1, mask);
+# endif
 
         // we interpolate the image and variance - include the mask or not?
@@ -141 +164 @@
     pmResiduals *resid = pmResidualsAlloc (xSize, ySize, xBin, yBin);
     psImageInit (resid->mask, 0);
+
+
+# if (TESTRESID)
+    psArray *inputSigCube = psArrayAllocEmpty (1);
+    psArray *inputVarCube = psArrayAllocEmpty (1);
+    psArray *inputMskCube = psArrayAllocEmpty (1);
+
+    // below, we are generating the residual images with interpolation
+    // to save these in an output cube, I need to generate place-holder images here
+    for (int i = 0; i < input->n; i++) {
+      psImage *inputSig = psImageAlloc (resid->Ro->numCols, resid->Ro->numRows, PS_TYPE_F32);
+      psImage *inputVar = psImageAlloc (resid->Ro->numCols, resid->Ro->numRows, PS_TYPE_F32);
+      psImage *inputMsk = psImageAlloc (resid->Ro->numCols, resid->Ro->numRows, PS_TYPE_IMAGE_MASK);
+
+      psImageInit (inputSig, NAN);
+      psImageInit (inputVar, NAN);
+      psImageInit (inputMsk, 0);
+
+      psArrayAdd (inputSigCube, 1, inputSig);
+      psArrayAdd (inputVarCube, 1, inputVar);
+      psArrayAdd (inputMskCube, 1, inputMsk);
+
+      // XXX free these here, right?
+      psFree (inputSig);
+      psFree (inputVar);
+      psFree (inputMsk);
+    }
+# endif
 
     // x(resid) = (x(image) - Xo)*xBin + xCenter
@@ -194 +245 @@
                     nGoodPixel ++;
                 }
+
+# if (TESTRESID)
+                psImage *inputSig = inputSigCube->data[i]; inputSig->data.F32[oy][ox] = flux;
+                psImage *inputVar = inputVarCube->data[i]; inputVar->data.F32[oy][ox] = flux;
+                psImage *inputMsk = inputMskCube->data[i]; inputMsk->data.PS_TYPE_IMAGE_MASK_DATA[oy][ox] = fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i];
+# endif
             }
+
+            if (nGoodPixel < 0.05*input->n) {
+              fprintf (stderr, "warning: %d, %d : residual pixel with few input pixels: %d vs %d\n", ox, oy, nGoodPixel, (int) input->n);
+            }
 
             // skip pixels with insufficient data
@@ -220 +281 @@
             }
 
+            if (fabs(fluxClip->robustMedian) > 1.5) {
+              fprintf (stderr, "warning: %d, %d : residual pixel has funny initial median value: %f\n", ox, oy, fluxClip->robustMedian);
+            }
+
             // mark input pixels which are more than N sigma from the median
             int nKeep = 0;
@@ -233 +298 @@
                 if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++;
             }
+
+            if (nKeep < 5) {
+              fprintf (stderr, "warning: %d, %d : residual pixel with few good pixels: %d vs %d\n", ox, oy, nKeep, (int) input->n);
+            }
 
             if (SPATIAL_ORDER == 0) {
@@ -242 +311 @@
                 }
 
+                // XXX this test should go at the top of the loop
                 float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
                 if (radius > radiusMax) {
@@ -304 +374 @@
                 float dRo = sqrt(A->data.F32[0][0]);
 
+                if (fabs(B->data.F64[0]) > 1.5) {
+                  fprintf (stderr, "warning: %d, %d : residual pixel has funny fit value: %f\n", ox, oy, B->data.F64[0]);
+                }
+
                 if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) {
                   resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
@@ -313 +387 @@
         }
     }
+
+# if (TESTRESID)
+
+    pmConfig *config = psphotGetConfig();
+    char *output = psMetadataLookupStr (&status, config->arguments, "OUTPUT");
+    
+    // XXX // save the star coordinate vectors
+    // XXX 
+    // XXX psString filename = NULL;
+    // XXX psStringAppend (&filename, "%s.instar.dat", output);
+
+    char filename[1024];
+    sprintf (filename, "%s.instar.dat", output);
+    
+    FILE *fout = fopen (filename, "w");
+    for (int i = 0; i < xC->n; i++) {
+      fprintf (fout, "%d %f %f\n", i, xC->data.F32[i], yC->data.F32[i]);
+    }
+    fclose (fout);
+
+    _psphotSaveCube (output, "raw.raw.fits", sourceRawCube);
+    _psphotSaveCube (output, "sig.raw.fits", sourceSigCube);
+    _psphotSaveCube (output, "var.raw.fits", sourceVarCube);
+    _psphotSaveCube (output, "msk.raw.fits", sourceMskCube);
+
+    _psphotSaveCube (output, "sig.fits", inputSigCube);
+    _psphotSaveCube (output, "var.fits", inputVarCube);
+    _psphotSaveCube (output, "msk.fits", inputMskCube);
+
+    psFree (sourceRawCube);
+    psFree (sourceSigCube);
+    psFree (sourceVarCube);
+    psFree (sourceMskCube);
+
+    psFree (inputSigCube);
+    psFree (inputVarCube);
+    psFree (inputMskCube);
+
+# endif
 
     psFree (A);
@@ -343 +456 @@
     return (resid != NULL) ? true : false;
 }
+
+bool _psphotSaveCube (char *basename, char *extname, psArray *cube) {
+
+    psString filename = NULL;
+    psStringAppend (&filename, "%s.%s", basename, extname);
+
+    psFits *fits = psFitsOpen (filename, "w");
+    if (!psFitsWriteImageCube (fits, NULL, cube, NULL)) {
+      fprintf (stderr, "failed to write the cube %s\n", filename);
+    }
+    psFitsClose (fits);
+
+    psFree (filename);
+
+    return true;
+}

branches/eam_branches/ipp-20230313/psphot/src/psphotOutput.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+pmConfig *staticConfig = NULL;
+
+void psphotSaveConfig (pmConfig *config) {
+  staticConfig = config;
+}
+
+pmConfig *psphotGetConfig () {
+  return staticConfig;
+}
 
 // convert filerule to filerule.NUM and look up in the config->arguments metadata