Changeset 5828 for trunk/psphot/src/psphotOutput.c

Timestamp:

Dec 22, 2005, 4:20:27 PM (21 years ago)

Author:

eugene

Message:

various cleanups, fixes, additions...

File:

• trunk/psphot/src/psphotOutput.c (modified) (7 diffs)

trunk/psphot/src/psphotOutput.c

@@ -5 +5 @@
 
     bool status;
+
+    psTimerStart ("psphot");
 
     char *outputMode = psMetadataLookupPtr (&status, config, "OUTPUT_MODE");
     char *outputFile = psMetadataLookupPtr (&status, config, "OUTPUT_FILE");
     char *residImage = psMetadataLookupPtr (&status, config, "RESID_IMAGE");
+    char *psfFile    = psMetadataLookupPtr (&status, config, "PSF_OUTPUT_FILE");
 
     if (residImage != NULL) psphotSaveImage (imdata->header, imdata->image, residImage);
+    if (psfFile != NULL) psphotWritePSF (psf, psfFile);
 
     if (outputFile == NULL) {
@@ -44 +48 @@
         return;
     }
+
     psAbort ("psphot", "unknown output mode %s", outputMode);
 }
@@ -287 +292 @@
 
     // set NAXIS to 0
-    fprintf (stderr, "setting naxis\n"); 
     mdi = psMetadataLookup (imdata->header, "NAXIS");
     mdi->data.S32 = 0;
@@ -540 +544 @@
     }
 
-    fprintf (stderr, "writing out moments\n");
     for (i = 0; i < sources->n; i++) {
         source = sources->data[i];
@@ -604 +607 @@
 
 
-bool psphotSamplePSFs (pmPSF *psf, psImage *image) {
-
-  // make sample PSFs for 4 corners and the center
+bool psphotSamplePSFs (psMetadata *config, pmPSF *psf, psImage *image) {
+
+    bool status;
+
+    // make sample PSFs for 4 corners and the center
     psImage *sample;
+
+    // optional dump of all rough source data
+    char *output = psMetadataLookupPtr (&status, config, "PSF_SAMPLE_FILE");
+    if (!status) return false;
+    if (output == NULL) return false;
+    if (output[0] == 0) return false;
 
     pmModel *modelFLT = pmModelAlloc (psf->type);
@@ -613 +624 @@
     modelFLT->params->data.F32[1] = 1;
 
-    psFits *fits = psFitsOpen ("sample.psf.fits", "w");
+    psFits *fits = psFitsOpen (output, "w");
 
     sample = pmModelPSFatXY (image, modelFLT, psf, 25, 25, 25, 25);
@@ -630 +641 @@
 }
 
+psPolynomial2D *psPolynomial2DfromMD (psMetadata *folder) {
+
+    bool status;
+    char keyword[80];
+
+    // get polynomial orders
+    // XXX add status failures tests
+    int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
+    int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
+
+    psPolynomial2D *poly = psPolynomial2DAlloc (nXorder, nYorder, PS_POLYNOMIAL_ORD);
+
+    for (int nx = 0; nx < poly->nX + 1; nx++) {
+        for (int ny = 0; ny < poly->nY + 1; ny++) {
+            sprintf (keyword, "VAL_X%02d_Y%02d", nx, ny);
+            poly->coeff[nx][ny] = psMetadataLookupF32 (&status, folder, keyword);
+            if (!status) poly->mask[nx][ny] = 1;
+
+            sprintf (keyword, "ERR_X%02d_Y%02d", nx, ny);
+            poly->coeffErr[nx][ny] = psMetadataLookupF32 (&status, folder, keyword);
+        }
+    }
+    return (poly);
+}
+
+psPolynomial3D *psPolynomial3DfromMD (psMetadata *folder) {
+
+    bool status;
+    char keyword[80];
+
+    // get polynomial orders
+    // XXX add status failures tests
+    int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
+    int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
+    int nZorder = psMetadataLookupS32 (&status, folder, "NORDER_Z");
+
+    psPolynomial3D *poly = psPolynomial3DAlloc (nXorder, nYorder, nZorder, PS_POLYNOMIAL_ORD);
+
+    for (int nx = 0; nx < poly->nX + 1; nx++) {
+        for (int ny = 0; ny < poly->nY + 1; ny++) {
+            for (int nz = 0; nz < poly->nZ + 1; nz++) {
+                sprintf (keyword, "VAL_X%02d_Y%02d_Z%02d", nx, ny, nz);
+                poly->coeff[nx][ny][nz] = psMetadataLookupF32 (&status, folder, keyword);
+                if (!status) poly->mask[nx][ny][nz] = 1;
+
+                sprintf (keyword, "ERR_X%02d_Y%02d_Z%02d", nx, ny, nz);
+                poly->coeffErr[nx][ny][nz] = psMetadataLookupF32 (&status, folder, keyword);
+            }
+        }
+    }
+    return (poly);
+}
+
+// XXX : these may need F64, or %g format for output
+bool psPolynomial2DtoMD (psMetadata *md, psPolynomial2D *poly, char *format, ...) {
+
+    int Nbyte;
+    char tmp;
+    char *root;
+    va_list argp;  
+
+    va_start (argp, format);
+    Nbyte = vsnprintf (&tmp, 0, format, argp);
+    va_end (argp);
+
+    if (!Nbyte) return false;
+
+    va_start (argp, format);
+    root = (char *) psAlloc (Nbyte + 1);
+    memset (root, 0, Nbyte + 1);
+    vsnprintf (root, Nbyte + 1, format, argp);
+    va_end (argp);
+
+    psMetadata *folder = psMetadataAlloc ();
+    psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 2D polynomial", folder);
+    psFree (root);
+
+    // specify the polynomial orders
+    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
+    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
+
+    // place polynomial entries on folder
+    for (int nx = 0; nx < poly->nX + 1; nx++) {
+        for (int ny = 0; ny < poly->nY + 1; ny++) {
+            if (poly->mask[nx][ny]) continue;
+            psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny], nx, ny);
+            psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d", PS_DATA_F32, "polynomial coefficient error", poly->coeffErr[nx][ny], nx, ny);
+        }
+    }
+    return true;
+}
+
+bool psPolynomial3DtoMD (psMetadata *md, psPolynomial3D *poly, char *format, ...) {
+
+    int Nbyte;
+    char tmp;
+    char *root;
+    va_list argp;  
+
+    va_start (argp, format);
+    Nbyte = vsnprintf (&tmp, 0, format, argp);
+    va_end (argp);
+
+    if (!Nbyte) return false;
+
+    va_start (argp, format);
+    root = (char *) psAlloc (Nbyte + 1);
+    memset (root, 0, Nbyte + 1);
+    vsnprintf (root, Nbyte + 1, format, argp);
+    va_end (argp);
+
+    psMetadata *folder = psMetadataAlloc ();
+    psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 3D polynomial", folder);
+    psFree (root);
+
+    // specify the polynomial orders
+    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
+    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
+    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Z", PS_DATA_S32, "number of z orders", poly->nZ);
+
+    // place polynomial entries on folder
+    for (int nx = 0; nx < poly->nX + 1; nx++) {
+        for (int ny = 0; ny < poly->nY + 1; ny++) {
+            for (int nz = 0; nz < poly->nZ + 1; nz++) {
+                if (poly->mask[nx][ny][nz]) continue;
+                psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d_Z%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny][nz], nx, ny, nz);
+                psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d_Z%02d", PS_DATA_F32, "polynomial coeffficient error", poly->coeffErr[nx][ny][nz], nx, ny, nz);
+            }
+        }
+    }
+    return true;
+}
+
+bool psphotWritePSF (pmPSF *psf, char *filename) {
+
+    psMetadata *psfdata = psMetadataAlloc ();
+
+    char *modelName = pmModelGetType (psf->type);
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName);
+
+    int nPar = pmModelParameterCount (psf->type)    ;
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
+
+    for (int i = 0; i < nPar - 4; i++) {
+        psPolynomial2D *poly = psf->params->data[i];
+        psPolynomial2DtoMD (psfdata, poly, "PSF_PAR%02d", i);
+    }
+    psPolynomial3DtoMD (psfdata, psf->ApTrend, "APTREND");
+
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid);
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid);
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
+
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq);
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars);
+    
+    psMetadataConfigWrite (psfdata, filename);
+    psFree (psfdata);
+
+    return true;
+}
+
+pmPSF *psphotReadPSF (char *filename) {
+
+    bool status;
+    unsigned int Nfail;
+    char keyword[80];
+
+    psMetadata *psfdata = psMetadataConfigParse (NULL, &Nfail, filename, FALSE);
+
+    char *modelName = psMetadataLookupPtr (&status, psfdata, "PSF_MODEL_NAME");
+    pmModelType type = pmModelSetType (modelName);
+
+    pmPSF *psf = pmPSFAlloc (type);
+
+    int nPar = psMetadataLookupS32 (&status, psfdata, "PSF_MODEL_NPAR");
+    if (nPar != pmModelParameterCount (psf->type)) psAbort ("read PSF" , "mismatch model par count");
+
+    for (int i = 0; i < nPar - 4; i++) {
+        sprintf (keyword, "PSF_PAR%02d", i);
+        psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword);
+        psPolynomial2D *poly = psPolynomial2DfromMD (folder);
+        psFree (psf->params->data[i]);
+        psf->params->data[i] = poly;
+    }
+    sprintf (keyword, "APTREND");
+    psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword);
+    psPolynomial3D *poly = psPolynomial3DfromMD (folder);
+    psFree (psf->ApTrend);
+    psf->ApTrend = poly;
+
+    psf->ApResid = psMetadataLookupF32 (&status, psfdata, "PSF_AP_RESID");
+    psf->dApResid = psMetadataLookupF32 (&status, psfdata, "PSF_dAP_RESID");
+    psf->skyBias = psMetadataLookupF32 (&status, psfdata, "PSF_SKY_BIAS");
+
+    psf->chisq = psMetadataLookupF32 (&status, psfdata, "PSF_CHISQ");
+    psf->nPSFstars = psMetadataLookupS32 (&status, psfdata, "PSF_NSTARS");
+
+    psFree (psfdata);
+    return (psf);
+}
+
 bool psphotDumpMoments (psMetadata *config, psArray *sources) {