IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 5, 2011, 11:02:53 AM (15 years ago)
Author:
eugene
Message:

merge updates from eam branch 20110404

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmPSF_IO.c

    r30031 r31451  
    6363
    6464bool pmPSFmodelReadPSFClump (psMetadata *analysis, psMetadata *header);
     65bool pmPSFmodelRead_ApTrend (pmPSF *psf, pmFPAfile *file);
     66bool pmPSFmodelWrite_ApTrend (pmFPAfile *file, pmPSF *psf);
     67
     68bool pmPSFmodelRead_GrowthCurve (pmPSF *psf, pmFPAfile *file);
     69bool pmPSFmodelWrite_GrowthCurve (pmFPAfile *file, pmPSF *psf);
    6570
    6671bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file)
    … …  
    197202        psError(psErrorCodeLast(), false, "Failed to write PSF for chip");
    198203        return false;
     204    }
     205    return true;
     206}
     207
     208// XXX we save the model term identifiers (item) as S32, but they probably should be more flexible
     209bool pmTrend2DtoTable (psArray *table, pmTrend2D *trend, char *label, int item) {
     210
     211    if (trend == NULL) return true;
     212
     213    if (trend->mode == PM_TREND_MAP) {
     214        // write the image components into a table: this is needed because they may each be a different size
     215        psImageMap *map = trend->map;
     216        for (int ix = 0; ix < map->map->numCols; ix++) {
     217            for (int iy = 0; iy < map->map->numRows; iy++) {
     218                psMetadata *row = psMetadataAlloc ();
     219                psMetadataAddS32 (row, PS_LIST_TAIL, label,        0, "", item);
     220                psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
     221                psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
     222                psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
     223                psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[iy][ix]);
     224                psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
     225
     226                psArrayAdd (table, 100, row);
     227                psFree (row);
     228            }
     229        }
     230    } else {
     231        // write the polynomial components into a table
     232        psPolynomial2D *poly = trend->poly;
     233        for (int ix = 0; ix <= poly->nX; ix++) {
     234            for (int iy = 0; iy <= poly->nY; iy++) {
     235                psMetadata *row = psMetadataAlloc ();
     236                psMetadataAddS32 (row, PS_LIST_TAIL, label,        0, "", item);
     237                psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
     238                psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
     239                psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
     240                psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
     241                psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->coeffMask[ix][iy]);
     242
     243                psArrayAdd (table, 100, row);
     244                psFree (row);
     245            }
     246        }
     247    }
     248    return true;
     249}
     250
     251// extra trend2D elements from a row
     252bool pmTrend2DfromTableRow (pmTrend2D *trend, psMetadata *row) {
     253
     254    bool status = false;
     255
     256    int xPow = psMetadataLookupS32 (&status, row, "X_POWER");
     257    int yPow = psMetadataLookupS32 (&status, row, "Y_POWER");
     258
     259    if (trend->mode == PM_TREND_MAP) {
     260        psImageMap *map = trend->map;
     261        assert (map);
     262        assert (map->map);
     263        assert (map->error);
     264        assert (xPow >= 0);
     265        assert (yPow >= 0);
     266        assert (xPow < map->map->numCols);
     267        assert (yPow < map->map->numRows);
     268        map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
     269        map->error->data.F32[yPow][xPow]  = psMetadataLookupF32 (&status, row, "ERROR");
     270    } else {
     271        psPolynomial2D *poly = trend->poly;
     272        assert (poly);
     273        assert (xPow >= 0);
     274        assert (yPow >= 0);
     275        assert (xPow <= poly->nX);
     276        assert (yPow <= poly->nY);
     277        poly->coeff[xPow][yPow]     = psMetadataLookupF32 (&status, row, "VALUE");
     278        poly->coeffErr[xPow][yPow]  = psMetadataLookupF32 (&status, row, "ERROR");
     279        poly->coeffMask[xPow][yPow] = psMetadataLookupU8  (&status, row, "MASK");
    199280    }
    200281    return true;
    … …  
    403484
    404485            if (trend->mode == PM_TREND_MAP) {
    405               nX = trend->map->map->numCols;
    406               nY = trend->map->map->numRows;
     486                nX = trend->map->map->numCols;
     487                nY = trend->map->map->numRows;
    407488            } else {
    408               nX = trend->poly->nX;
    409               nY = trend->poly->nY;
     489                nX = trend->poly->nX;
     490                nY = trend->poly->nY;
    410491            }
    411492            snprintf (name, 9, "PAR%02d_NX", i);
    … …  
    428509        psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
    429510
     511        float PSF_APERTURE =  psMetadataLookupF32(&status, roAnalysis, "PSF_APERTURE");
     512        if (status) {
     513            psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_APERTURE", PS_DATA_F32, "aperture for psf objects", PSF_APERTURE);
     514        }
     515        float PSF_FIT_RADIUS =  psMetadataLookupF32(&status, roAnalysis, "PSF_FIT_RADIUS");
     516        if (status) {
     517            psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_DATA_F32, "aperture for psf objects", PSF_FIT_RADIUS);
     518        }
     519
    430520        // build a FITS table of the PSF parameters
    431521        psArray *psfTable = psArrayAllocEmpty (100);
    432522        for (int i = 0; i < nPar; i++) {
    433523            pmTrend2D *trend = psf->params->data[i];
    434             if (trend == NULL) continue; // skip unset parameters (eg, XPOS)
    435 
    436             if (trend->mode == PM_TREND_MAP) {
    437                 // write the image components into a table: this is needed because they may each be a different size
    438                 psImageMap *map = trend->map;
    439                 for (int ix = 0; ix < map->map->numCols; ix++) {
    440                     for (int iy = 0; iy < map->map->numRows; iy++) {
    441                         psMetadata *row = psMetadataAlloc ();
    442                         psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
    443                         psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
    444                         psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
    445                         psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
    446                         psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[iy][ix]);
    447                         psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
    448 
    449                         psArrayAdd (psfTable, 100, row);
    450                         psFree (row);
    451                     }
    452                 }
    453             } else {
    454                 // write the polynomial components into a table
    455                 psPolynomial2D *poly = trend->poly;
    456                 for (int ix = 0; ix <= poly->nX; ix++) {
    457                     for (int iy = 0; iy <= poly->nY; iy++) {
    458                         psMetadata *row = psMetadataAlloc ();
    459                         psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
    460                         psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
    461                         psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
    462                         psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
    463                         psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
    464                         psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->coeffMask[ix][iy]);
    465 
    466                         psArrayAdd (psfTable, 100, row);
    467                         psFree (row);
    468                     }
    469                 }
    470             }
     524            pmTrend2DtoTable (psfTable, trend, "MODEL_TERM", i);
    471525        }
    472526
    … …  
    547601    }
    548602
     603    if (!pmPSFmodelWrite_ApTrend(file, psf)) {
     604        psError(psErrorCodeLast(), false, "Unable to write PSF ApTrend");
     605        return false;
     606    }
     607
     608    if (!pmPSFmodelWrite_GrowthCurve(file, psf)) {
     609        psError(psErrorCodeLast(), false, "Unable to write PSF Growth Curve");
     610        return false;
     611    }
     612
    549613    // write a representation of the psf model
    550614    {
    551615        psMetadata *header = psMetadataAlloc ();
    552 
    553         if (0) {
    554             // set some header keywords to make it clear there are no residuals?
    555             if (!psFitsWriteBlank (file->fits, header, residName)) {
    556                 psError(psErrorCodeLast(), false, "Unable to write blank PSF residuals.");
    557                 psFree(residName);
    558                 psFree(header);
    559                 return false;
    560             }
    561             psFree (residName);
    562             psFree (header);
    563             return true;
    564         }
    565616
    566617        int DX = 65;
    … …  
    651702}
    652703
     704
     705
    653706// if this file needs to have a PHU written out, write one
    654707bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config)
    … …  
    934987    psf->skyBias   = psMetadataLookupF32 (&status, header, "SKY_BIAS");
    935988
     989    if (roAnalysis) {
     990        float PSF_APERTURE =  psMetadataLookupF32(&status, header, "PSF_APERTURE");
     991        if (status) {
     992            psMetadataAddF32 (roAnalysis, PS_LIST_TAIL, "PSF_APERTURE", PS_DATA_F32, "aperture for psf objects", PSF_APERTURE);
     993        }
     994        float PSF_FIT_RADIUS =  psMetadataLookupF32(&status, header, "PSF_FIT_RADIUS");
     995        if (status) {
     996            psMetadataAddF32 (roAnalysis, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_DATA_F32, "aperture for psf objects", PSF_FIT_RADIUS);
     997        }
     998    } else {
     999        psWarning ("unable to read PSF_APERTURE or PSF_FIT_RADIUS");
     1000    }
     1001
    9361002    // read the raw table data
    9371003    psArray *table = psFitsReadTable (file->fits);
    … …  
    9451011    for (int i = 0; i < table->n; i++) {
    9461012        psMetadata *row = table->data[i];
     1013
    9471014        int iPar = psMetadataLookupS32 (&status, row, "MODEL_TERM");
    948         int xPow = psMetadataLookupS32 (&status, row, "X_POWER");
    949         int yPow = psMetadataLookupS32 (&status, row, "Y_POWER");
    9501015
    9511016        pmTrend2D *trend = psf->params->data[iPar];
    … …  
    9551020        }
    9561021
    957         if (trend->mode == PM_TREND_MAP) {
    958             psImageMap *map = trend->map;
    959             assert (map);
    960             assert (map->map);
    961             assert (map->error);
    962             assert (xPow >= 0);
    963             assert (yPow >= 0);
    964             assert (xPow < map->map->numCols);
    965             assert (yPow < map->map->numRows);
    966             map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
    967             map->error->data.F32[yPow][xPow]  = psMetadataLookupF32 (&status, row, "ERROR");
    968         } else {
    969             psPolynomial2D *poly = trend->poly;
    970             assert (poly);
    971             assert (xPow >= 0);
    972             assert (yPow >= 0);
    973             assert (xPow <= poly->nX);
    974             assert (yPow <= poly->nY);
    975             poly->coeff[xPow][yPow]     = psMetadataLookupF32 (&status, row, "VALUE");
    976             poly->coeffErr[xPow][yPow]  = psMetadataLookupF32 (&status, row, "ERROR");
    977             poly->coeffMask[xPow][yPow] = psMetadataLookupU8  (&status, row, "MASK");
    978         }
     1022        pmTrend2DfromTableRow(trend, row);
    9791023    }
    9801024    psFree (header);
    … …  
    10631107    }
    10641108
     1109    if (!pmPSFmodelRead_ApTrend (psf, file)) {
     1110        psError(psErrorCodeLast(), false, "Unable to read PSF ApTrend data.");
     1111        return false;
     1112    }
     1113
     1114    if (!pmPSFmodelRead_GrowthCurve(psf, file)) {
     1115        psError(psErrorCodeLast(), false, "Unable to read PSF Growth Curve");
     1116        return false;
     1117    }
     1118
    10651119    psMetadataAdd (chipAnalysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
    10661120    psFree (psf);
    … …  
    10701124    psFree (header);
    10711125
     1126    return true;
     1127}
     1128
     1129// write aperture trend to a FITS table
     1130bool pmPSFmodelWrite_ApTrend (pmFPAfile *file, pmPSF *psf) {
     1131
     1132    pmTrend2D *trend = psf->ApTrend;
     1133    if (trend == NULL) {
     1134        psWarning ("no PSF ApTrend to write out, skipping");
     1135        return true;
     1136    }
     1137
     1138    // we need to write a header for the table,
     1139    psMetadata *header = psMetadataAlloc();
     1140
     1141    int nX = 0, nY = 0;
     1142    if (trend->mode == PM_TREND_MAP) {
     1143        nX = trend->map->map->numCols;
     1144        nY = trend->map->map->numRows;
     1145    } else {
     1146        nX = trend->poly->nX;
     1147        nY = trend->poly->nY;
     1148    }
     1149    psMetadataAddS32 (header, PS_LIST_TAIL, "TREND_NX", 0, "", nX);
     1150    psMetadataAddS32 (header, PS_LIST_TAIL, "TREND_NY", 0, "", nY);
     1151    char *modeName = pmTrend2DModeToString (trend->mode);
     1152    psMetadataAddStr (header, PS_LIST_TAIL, "TREND_MD", 0, "", modeName);
     1153    psFree (modeName);
     1154
     1155    // build a FITS table of the ApTrend (only 1)
     1156    psArray *table = psArrayAllocEmpty (100);
     1157    pmTrend2DtoTable (table, trend, "APTREND", 0);
     1158
     1159    // write an empty FITS segment if we have no PSF information
     1160    if (table->n == 0) {
     1161        psError(PM_ERR_PROG, true, "No PSF data to write.");
     1162        psFree(table);
     1163        psFree(header);
     1164        return false;
     1165    }
     1166
     1167    psTrace ("pmFPAfile", 5, "writing psf ApTrend data %s\n", "AP_TREND");
     1168    if (!psFitsWriteTable(file->fits, header, table, "AP_TREND")) {
     1169        psError(psErrorCodeLast(), false, "Error writing psf table data %s\n", "AP_TREND");
     1170        psFree(table);
     1171        psFree(header);
     1172        return false;
     1173    }
     1174
     1175    psFree (table);
     1176    psFree (header);
     1177    return true;
     1178}
     1179
     1180// read aperture trend to a FITS table
     1181bool pmPSFmodelRead_ApTrend (pmPSF *psf, pmFPAfile *file) {
     1182
     1183    bool status;
     1184
     1185    // move fits pointer to AP_TREND section
     1186    // advance to the table data extension
     1187    if (!psFitsMoveExtNameClean (file->fits, "AP_TREND")) {
     1188        psWarning ("no Aperture Trend data in PSF file, skipping");
     1189        return true;
     1190    }
     1191
     1192    psMetadata *header = psFitsReadHeader (NULL, file->fits);
     1193    if (!header) {
     1194        psError(psErrorCodeLast(), false, "Unable to read AP_TREND header.");
     1195        return false;
     1196    }
     1197       
     1198    // read the raw table data
     1199    psArray *table = psFitsReadTable (file->fits);
     1200    if (!table) {
     1201        psError(psErrorCodeLast(), false, "Unable to read AP_TREND table.");
     1202        psFree(header);
     1203        return false;
     1204    }
     1205
     1206    // XXX allow user to set this optionally?
     1207    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     1208
     1209    psImageBinning *binning = psImageBinningAlloc();
     1210    binning->nXfine = psf->fieldNx;
     1211    binning->nYfine = psf->fieldNy;
     1212    binning->nXruff = psMetadataLookupS32 (&status, header, "TREND_NX");
     1213    binning->nYruff = psMetadataLookupS32 (&status, header, "TREND_NY");
     1214    psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
     1215    char *modeName  = psMetadataLookupStr (&status, header, "TREND_MD");
     1216    if (!status) {
     1217        psError(PM_ERR_PROG, true, "inconsistent PSF header: NX & NY defined for AP TREND, but not MD");
     1218        psFree (header);
     1219        psFree (stats);
     1220        psFree (table);
     1221        return false;
     1222    }
     1223    pmTrend2DMode psfTrendMode = pmTrend2DModeFromString (modeName);
     1224    if (psfTrendMode == PM_TREND_NONE) {
     1225        psfTrendMode = PM_TREND_POLY_ORD;
     1226    }
     1227
     1228    // measure Trend2D for the current spatial scale
     1229    pmTrend2D *apTrend = pmTrend2DNoImageAlloc (PM_TREND_MAP, binning, stats);
     1230
     1231    // fill in the matching psf->params entries
     1232    for (int i = 0; i < table->n; i++) {
     1233        psMetadata *row = table->data[i];
     1234        pmTrend2DfromTableRow(apTrend, row);
     1235    }
     1236    psf->ApTrend = apTrend;
     1237
     1238    psFree (binning);
     1239    psFree (header);
     1240    psFree (stats);
     1241    psFree (table);
     1242    return true;
     1243}
     1244
     1245// write aperture trend to a FITS table
     1246bool pmPSFmodelWrite_GrowthCurve (pmFPAfile *file, pmPSF *psf) {
     1247
     1248    pmGrowthCurve *growth = psf->growth;
     1249    if (growth == NULL) {
     1250        psWarning ("no PSF Growth Curve to write out, skipping");
     1251        return true;
     1252    }
     1253
     1254    // we need to write a header for the table,
     1255    psMetadata *header = psMetadataAlloc();
     1256
     1257    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_MIN_RAD", 0, "", growth->radius->data.F32[0]);
     1258    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_MAX_RAD", 0, "", growth->maxRadius);
     1259    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_REF_RAD", 0, "", growth->refRadius);
     1260    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_AP_LOSS", 0, "", growth->apLoss);
     1261    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_AP_REF",  0, "", growth->apRef);
     1262    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_FIT_MAG", 0, "", growth->fitMag);
     1263
     1264    // build a FITS table of the ApTrend (only 1)
     1265    psArray *table = psArrayAllocEmpty (100);
     1266    for (int i = 0; i < growth->apMag->n; i++) {
     1267        psMetadata *row = psMetadataAlloc ();
     1268        psMetadataAddF32 (row, PS_LIST_TAIL, "RADIUS", 0, "", growth->radius->data.F32[i]);
     1269        psMetadataAddF32 (row, PS_LIST_TAIL, "AP_MAG", 0, "", growth->apMag->data.F32[i]);
     1270        psArrayAdd (table, 100, row);
     1271        psFree (row);
     1272    }
     1273
     1274    // write an empty FITS segment if we have no PSF information
     1275    if (table->n == 0) {
     1276        psError(PM_ERR_PROG, true, "No PSF data to write.");
     1277        psFree(table);
     1278        psFree(header);
     1279        return false;
     1280    }
     1281
     1282    psTrace ("pmFPAfile", 5, "writing psf Growth Curve data %s\n", "GROWTH_CURVE");
     1283    if (!psFitsWriteTable(file->fits, header, table, "GROWTH_CURVE")) {
     1284        psError(psErrorCodeLast(), false, "Error writing psf table data %s\n", "GROWTH_CURVE");
     1285        psFree(table);
     1286        psFree(header);
     1287        return false;
     1288    }
     1289
     1290    psFree (table);
     1291    psFree (header);
     1292    return true;
     1293}
     1294
     1295// read aperture trend to a FITS table
     1296bool pmPSFmodelRead_GrowthCurve (pmPSF *psf, pmFPAfile *file) {
     1297
     1298    bool status;
     1299
     1300    // move fits pointer to AP_TREND section
     1301    // advance to the table data extension
     1302    if (!psFitsMoveExtNameClean (file->fits, "GROWTH_CURVE")) {
     1303        psWarning ("no Growth Curve data in PSF file, skipping");
     1304        return true;
     1305    }
     1306
     1307    psMetadata *header = psFitsReadHeader (NULL, file->fits);
     1308    if (!header) {
     1309        psError(psErrorCodeLast(), false, "Unable to read GROWTH_CURVE header.");
     1310        return false;
     1311    }
     1312       
     1313    // read the raw table data
     1314    psArray *table = psFitsReadTable (file->fits);
     1315    if (!table) {
     1316        psError(psErrorCodeLast(), false, "Unable to read GROWTH_CURVE table.");
     1317        psFree(header);
     1318        return false;
     1319    }
     1320
     1321    float minRadius = psMetadataLookupF32 (&status, header, "GROWTH_MIN_RAD"); if (!status) return false;
     1322    float maxRadius = psMetadataLookupF32 (&status, header, "GROWTH_MAX_RAD"); if (!status) return false;
     1323    float refRadius = psMetadataLookupF32 (&status, header, "GROWTH_REF_RAD"); if (!status) return false;
     1324
     1325    psf->growth = pmGrowthCurveAlloc(minRadius, maxRadius, refRadius);
     1326
     1327    psf->growth->apLoss = psMetadataLookupF32 (&status, header, "GROWTH_AP_LOSS"); if (!status) return false;
     1328    psf->growth->apRef  = psMetadataLookupF32 (&status, header, "GROWTH_AP_REF"); if (!status) return false;
     1329    psf->growth->fitMag = psMetadataLookupF32 (&status, header, "GROWTH_FIT_MAG"); if (!status) return false;
     1330
     1331    // fill in the matching psf->params entries
     1332    for (int i = 0; i < table->n; i++) {
     1333        psMetadata *row = table->data[i];
     1334        psf->growth->apMag->data.F32[i] = psMetadataLookupF32 (&status, row, "AP_MAG"); if (!status) return false;
     1335    }
     1336
     1337    psFree (header);
     1338    psFree (table);
    10721339    return true;
    10731340}
Note: See TracChangeset for help on using the changeset viewer.