IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 11, 2006, 5:27:13 PM (20 years ago)
Author:
eugene
Message:

major revisions to use the new pmFPAfile,view I/O handling

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotOutput.c

    r6522 r6571  
    11# include "psphot.h"
    2 
    3 static int   extNumber  = -1;
    4 
    5 static char *outputName = NULL;
    6 static char *outputRoot = NULL;
    7 static char *outputMode = NULL;
    8 static char *outputFormat = NULL;
    9 
    10 static char *extNumberFormat = NULL;
    11 static char *extNameKey = NULL;
    12 static char *extRoot    = NULL;
    13 static char *psfFile    = NULL;
    14 static char *psfSample  = NULL;
    15 
    16 static psFits *outputFits = NULL;
    17 static FILE   *outputFile = NULL;
    18 
    19 bool psphotDataIO (psphotData *data, ppConfig *config, int chip, int cell) {
    20 
    21     // load input data, if depth is appropriate
    22     ppImageLoadDepth loadDepth = ppImageCheckDepth (config->recipe, "LOAD.DEPTH");
    23     if ((chip == -1) && (cell == -1) && (loadDepth == PP_LOAD_FPA)) {
    24         psTrace(__func__, 1, "Loading pixels for FPA...\n");
    25         ppImageLoadPixels(data->input, config->database, -1, -1);
    26     }
    27     if ((chip != -1) && (cell == -1) && (loadDepth == PP_LOAD_CHIP)) {
    28         psTrace(__func__, 1, "Loading pixels for chip %d...\n", i);
    29         ppImageLoadPixels(data->input, config->database, chip, -1);
    30     }
    31     if ((chip != -1) && (cell != -1) && (loadDepth == PP_LOAD_CELL)) {
    32         psTrace(__func__, 1, "Loading pixels for chip %d, cell %d...\n", i, j);
    33         ppImageLoadPixels(data->input, config->database, chip, cell);
    34     }
    35 
    36 
    37 
    38     if (saveDepth == PP_LOAD_FPA) {
    39         psTrace(__func__, 1, "Saving objects for FPA...\n");
    40         psphotOutputOpen (file->fpa->phu);
    41     }
    42 
    43     ppImageLoadDepth saveDepth = ppImageCheckDepth (config->recipe, "SAVE.DEPTH");
    44    
    45 
    46 }
    47 
    48 void psphotOutputPrep (ppFile *file, ppConfig *config) {
    49 
    50     bool status;
    51 
    52     outputRoot      = psMetadataLookupStr (&status, config->arguments, "OUTPUT_ROOT");
    53     if (!status) psAbort ("psphot", "output file not specified");
    54 
    55     outputName      = psMetadataLookupStr (&status, config->recipe, "OUTPUT_NAME");
    56     outputMode      = psMetadataLookupStr (&status, config->recipe, "OUTPUT_MODE");
    57     outputFormat    = psMetadataLookupStr (&status, config->recipe, "OUTPUT_FORMAT");
    58 
    59     psfFile         = psMetadataLookupStr (&status, config->recipe, "PSF_OUTPUT_FILE");
    60     psfSample       = psMetadataLookupStr (&status, config->recipe, "PSF_SAMPLE_FILE");
    61 
    62     // rules to construct output names
    63     extNumberFormat = psMetadataLookupStr (&status, config->recipe, "EXTNUMBER_FORMAT");
    64     extNameKey      = psMetadataLookupStr (&status, config->recipe, "EXTNAME");
    65     extRoot         = psMetadataLookupStr (&status, config->recipe, "EXTROOT");
    66 
    67     headExtname     = psMetadataLookupStr (&status, config->recipe, "HEAD_EXTNAME");
    68     dataExtname     = psMetadataLookupStr (&status, config->recipe, "DATA_EXTNAME");
    69 
    70     if (extNumberFormat == NULL) {
    71         extNumberFormat = psStringCopy ("%02d");
    72     }
    73 
    74     psStringStrip (outputName);
    75     psStringStrip (extNameKey);
    76     psStringStrip (extRoot);
    77 
    78     // validate the outputMode and outputFormat
    79     status = false;
    80     status |= !strcasecmp (outputFormat, "TEXT");
    81     status |= !strcasecmp (outputFormat, "OBJ");
    82     status |= !strcasecmp (outputFormat, "SX");
    83     status |= !strcasecmp (outputFormat, "CMP");
    84     status |= !strcasecmp (outputFormat, "CMF");
    85     if (!status) {
    86         psLogMsg ("psphot", PS_LOG_ERROR, "invalid output format %s", outputFormat);
    87         exit (1);
    88     }
    89 
    90     // validate the outputMode and outputFormat
    91     status = false;
    92     status |= !strcasecmp (outputMode, "MEF");
    93     status |= !strcasecmp (outputMode, "SPLIT");
    94     if (!status) {
    95         psLogMsg ("psphot", PS_LOG_ERROR, "invalid output mode %s", outputMode);
    96         exit (1);
    97     }
    98 
    99     // for MEF output, we only allow certain output formats
    100     if (!strcasecmp (outputMode, "MEF")) {
    101         if (!strcasecmp (outputFormat, "CMF")) goto valid_format;
    102         psLogMsg ("psphot", PS_LOG_ERROR, "MEF output mode unavailable");
    103         exit (1);
    104     }
    105 valid_format:
    106 
    107     // register background image-file names
    108     psphotBackgroundNames (config->arguments);
    109 
    110     // save so freeing config will not drop these references
    111     psMemCopy (outputRoot);
    112     psMemCopy (outputName);
    113     psMemCopy (outputMode);
    114     psMemCopy (outputFormat);
    115     psMemCopy (extNumberFormat);
    116     psMemCopy (extNameKey);
    117     psMemCopy (extRoot);
    118     psMemCopy (psfFile);
    119     psMemCopy (psfSample);
    120 
    121     return;
    122 }
    123 
    124 bool psphotOutputOpen (psMetadata *phu){
    125 
    126     /* outputFile is used by text-format outputs */
    127 
    128     psFree (outputFile);
    129     outputFile = psphotOutputName (phu, outputName);
    130 
    131     /* for fits-format output, open outputFits */
    132     if (!strcasecmp (outputFormat, "CMF")) {
    133         if (outputFits != NULL) {
    134             psFitsClose (outputFits);
    135         }
    136         outputFits = psFitsOpen (outputFile, "w");
    137 
    138         if (!strcasecmp (outputMode, "MEF")) {
    139             psFitsWriteHeader (outputFits, phu);
    140         }
    141     }
    142     return true;
    143 }
    144 
    145 bool psphotOutputClose (void) {
    146 
    147     if (outputFile != NULL) {
    148         psFitsClose (outputFile);
    149     }
    150     return true;
    151 }
    1522
    1533// output functions: we have several fixed modes (sx, obj, cmp)
    … …  
    1599    psMetadata *header = pmReadoutGetHeader (readout);
    16010
    161     psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
    16211    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    163 
    164     // starting value is -1
    165     extNumber ++;
    166 
    167     // for SPLIT output, we need to open a file for each output call
    168     if (!strcasecmp (outputMode, "SPLIT")) {
    169         outputFile = psphotSplitName (header);
    170     }
    171 
    172     if (!strcasecmp (outputMode, "MEF")) {
    173         headExt = psphotOutputName (header, headExtname);
    174         dataExt = psphotOutputName (header, dataExtname);
    175     }
    176     if (!strcasecmp (outputMode, "SPLIT")) {
    177         dataExt = psphotOutputName (header, dataExtname);
    178     }
    179 
    180     char *residImage = psMetadataLookupStr (&status, arguments, "RESID_IMAGE");
     12    // sample PSF images??
    18113    if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample);
    182     if (residImage != NULL) psphotSaveImage (header, readout->image, residImage);
     14    if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample);
    18315
    18416    if (psfFile != NULL) {
    … …  
    18719        psFree (psfData);
    18820    }
    189     if (outputFile == NULL) {
    190         psLogMsg ("output", 3, "no data output file selected");
    191         return;
    192     }
    193     if (outputFormat == NULL) {
    194         psLogMsg ("output", 3, "no data output format selected");
    195         psFree (outputFile);
    196         return;
    197     }
    198     if (!strcasecmp (outputFormat, "SX")) {
    199         pmSourcesWriteSX (sources, outputFile);
    200         psFree (outputFile);
    201         return;
    202     }
    203     if (!strcasecmp (outputFormat, "OBJ")) {
    204         pmSourcesWriteOBJ (sources, outputFile);
    205         psFree (outputFile);
    206         return;
    207     }
    208     if (!strcasecmp (outputFormat, "CMP")) {
    209         pmSourcesWriteCMP (sources, outputFile, header);
    210         psFree (outputFile);
    211         return;
    212     }
    213     if (!strcasecmp (outputFormat, "CMF")) {
    214         pmSourcesWriteCMF (sources, header, headExt, dataExt);
    215         psFree (outputFile);
    216         return;
    217     }
    218     if (!strcasecmp (outputFormat, "TEXT")) {
    219         pmSourcesWriteText (sources, outputFile);
    220         psFree (outputFile);
    221         return;
    222     }
    223 
    224     psAbort ("psphot", "unknown output mode %s", outputFormat);
    225 }
    226 
    227 // dophot-style output list with fixed line width
    228 bool pmSourcesWriteOBJ (psArray *sources, char *filename) {
    229 
    230     int type;
    231     psF32 *PAR, *dPAR;
    232     float dmag, apResid;
    233 
    234     psTimerStart ("string");
    235 
    236     psLine *line = psLineAlloc (104);  // 104 is dophot-defined line length
    237 
    238     FILE *f = fopen (filename, "w");
    239     if (f == NULL) {
    240         psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename);
    241         return false;
    242     }
    243 
    244     // write sources with models
    245     for (int i = 0; i < sources->n; i++) {
    246         pmSource *source = (pmSource *) sources->data[i];
    247         pmModel *model = pmModelSelect (source);
    248         if (model == NULL) continue;
    249 
    250         PAR = model->params->data.F32;
    251         dPAR = model->dparams->data.F32;
    252 
    253         dmag = dPAR[1] / PAR[1];
    254         type = pmSourceDophotType (source);
    255         apResid = source->apMag - source->fitMag;
    256 
    257         psLineInit (line);
    258         psLineAdd (line, "%3d",   type);
    259         psLineAdd (line, "%8.2f", PAR[2]);
    260         psLineAdd (line, "%8.2f", PAR[3]);
    261         psLineAdd (line, "%8.3f", source->fitMag);
    262         psLineAdd (line, "%6.3f", dmag);
    263         psLineAdd (line, "%9.2f", PAR[0]);
    264         psLineAdd (line, "%9.3f", PAR[4]);
    265         psLineAdd (line, "%9.3f", PAR[5]);
    266         psLineAdd (line, "%7.2f", PAR[6]);
    267         psLineAdd (line, "%8.3f", 99.999);
    268         psLineAdd (line, "%8.3f", source->apMag);
    269         psLineAdd (line, "%8.2f\n", apResid);
    270         fwrite (line->line, 1, line->Nline, f);
    271     }
    272     fclose (f);
    273     psFree (line);
    274     fprintf (stderr, "%f seconds for %d objects with psLine\n", psTimerMark ("string"), (int)sources->n);
    275 
    276     psTimerStart ("string");
    277 
    278     f = fopen ("test.obj", "w");
    279     if (f == NULL) {
    280         psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", "test.obj");
    281         return false;
    282     }
    283 
    284     char *string;
    285     // write sources with models
    286     for (int i = 0; i < sources->n; i++) {
    287         pmSource *source = (pmSource *) sources->data[i];
    288         pmModel *model = pmModelSelect (source);
    289         if (model == NULL) continue;
    290 
    291         PAR = model->params->data.F32;
    292         dPAR = model->dparams->data.F32;
    293 
    294         dmag = dPAR[1] / PAR[1];
    295         type = pmSourceDophotType (source);
    296         apResid = source->apMag - source->fitMag;
    297 
    298         string = NULL;
    299         psStringAppend (&string, "%3d",   type);
    300         psStringAppend (&string, "%8.2f", PAR[2]);
    301         psStringAppend (&string, "%8.2f", PAR[3]);
    302         psStringAppend (&string, "%8.3f", source->fitMag);
    303         psStringAppend (&string, "%6.3f", dmag);
    304         psStringAppend (&string, "%9.2f", PAR[0]);
    305         psStringAppend (&string, "%9.3f", PAR[4]);
    306         psStringAppend (&string, "%9.3f", PAR[5]);
    307         psStringAppend (&string, "%7.2f", PAR[6]);
    308         psStringAppend (&string, "%8.3f", 99.999);
    309         psStringAppend (&string, "%8.3f", source->apMag);
    310         psStringAppend (&string, "%8.2f\n", apResid);
    311         fwrite (string, 1, strlen(string), f);
    312         psFree (string);
    313     }
    314     fclose (f);
    315     fprintf (stderr, "%f seconds for %d objects with psString\n", psTimerMark ("string"), (int)sources->n);
    316 
    317     return true;
    318 }
    319 
    320 // elixir-mode / sextractor-style output list with fixed line width
    321 bool pmSourcesWriteSX (psArray *sources, char *filename) {
    322 
    323     psF32 *PAR, *dPAR;
    324     float dmag;
    325 
    326     psLine *line = psLineAlloc (110);  // 110 is sextractor line length
    327 
    328     FILE *f = fopen (filename, "w");
    329     if (f == NULL) {
    330         psLogMsg ("WriteSourceSX", 3, "can't open output file for output %s\n", filename);
    331         return false;
    332     }
    333 
    334     // write sources with models
    335     for (int i = 0; i < sources->n; i++) {
    336         pmSource *source = (pmSource *) sources->data[i];
    337         pmModel *model = pmModelSelect (source);
    338         if (model == NULL) continue;
    339 
    340         PAR = model->params->data.F32;
    341         dPAR = model->dparams->data.F32;
    342 
    343         // pmSourceSextractType (source, &type, &flags); XXX EAM : implement me...
    344 
    345         dmag = dPAR[1] / PAR[1];
    346 
    347         psLineInit (line);
    348         psLineAdd (line, "%5.2f", 0.0); // should be type
    349         psLineAdd (line, "%11.3f", PAR[2]);
    350         psLineAdd (line, "%11.3f", PAR[3]);
    351         psLineAdd (line, "%9.4f", source->fitMag);
    352         psLineAdd (line, "%9.4f", dmag);
    353         psLineAdd (line, "%13.4f", PAR[0]);
    354         psLineAdd (line, "%9.2f", 0.0); // should be FWHMx
    355         psLineAdd (line, "%9.2f", 0.0); // should be FWHMy
    356         psLineAdd (line, "%6.1f", 0.0); // should be Theta
    357         psLineAdd (line, "%9.4f", 99.999); // should be MAG_ISO
    358         psLineAdd (line, "%9.4f", source->apMag);
    359         psLineAdd (line, "%4d\n", 0); // should be flags
    360         fwrite (line->line, 1, line->Nline, f);
    361     }
    362     fclose (f);
    363     return true;
    364 }
    365 
    366 // elixir-style pseudo FITS table (header + ascii list)
    367 bool pmSourcesWriteCMP (psArray *sources, char *filename, psMetadata *header) {
    368 
    369     int i, type;
    370     psMetadataItem *mdi;
    371     psF32 *PAR, *dPAR;
    372     float dmag, lsky;
    373     bool status;
    374 
    375     // find config information for output header
    376     float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT");
    377     if (!status) ZERO_POINT = 25.0;
    378 
    379     // write necessary information to output header
    380     psMetadataAdd (header, PS_LIST_TAIL, "NSTARS",   PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n);
    381 
    382     // create file, write-out header
    383     psFits *fits = psFitsOpen (filename, "w");
    384 
    385     // set NAXIS to 0 : CFITSIO requires isolated header to have NAXIS = 0
    386     mdi = psMetadataLookup (header, "NAXIS");
    387     mdi->data.S32 = 0;
    388     mdi->type = PS_DATA_S32;
    389 
    390     // psMetadataAdd (header, PS_LIST_HEAD, "NAXIS", PS_DATA_S32 | PS_META_REPLACE, "", 0);
    391     psFitsWriteHeader (header, fits);
    392     psFitsClose (fits);
    393 
    394     // re-open, add data to end of file
    395     FILE *f = fopen (filename, "a+");
    396     if (f == NULL) {
    397         psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename);
    398         return false;
    399     }
    400     fseek (f, 0, SEEK_END);
    401 
    402     psLine *line = psLineAlloc (67);  // 66 is imclean-defined line length
    403 
    404     // write sources with models first
    405     for (i = 0; i < sources->n; i++) {
    406         pmSource *source = (pmSource *) sources->data[i];
    407         pmModel *model = pmModelSelect (source);
    408         if (model == NULL) continue;
    409 
    410         PAR = model->params->data.F32;
    411         dPAR = model->dparams->data.F32;
    412 
    413         dmag = dPAR[1] / PAR[1];
    414         type = pmSourceDophotType (source);
    415         lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]);
    416 
    417         psLineInit (line);
    418         psLineAdd (line, "%6.1f ", PAR[2]);
    419         psLineAdd (line, "%6.1f ", PAR[3]);
    420         psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->fitMag + ZERO_POINT));
    421         psLineAdd (line, "%03d ",  PS_MIN (999, (int)(1000*dmag)));
    422         psLineAdd (line, "%2d ",   type);
    423         psLineAdd (line, "%3.1f ", lsky);
    424         psLineAdd (line, "%6.3f ", 99.999); // should be 'Mgal
    425         psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->apMag + ZERO_POINT));
    426         psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x'
    427         psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y'
    428         psLineAdd (line, "%5.1f\n", 0); // should be theta
    429         fwrite (line->line, 1, line->Nline, f);
    430     }
    431     fclose (f);
    432     return true;
    433 }
    434 
    435 // elixir-style FITS table output (header + table in 1st extension)
    436 // this format consists of a header derived from the image header
    437 // followed by a zero-size matrix, followed by the table data
    438 bool pmSourcesWriteCMF (psArray *sources, char *filename, psMetadata *header) {
    439 
    440     psArray *table;
    441     psMetadataItem *mdi;
    442     psMetadata *row;
    443     psMetadata *theader;
    444     int i, type;
    445     psF32 *PAR, *dPAR;
    446     float dmag, lsky;
    447     bool status;
    448 
    449     // find config information for output header
    450     float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT");
    451 
    452     table = psArrayAlloc (sources->n);
    453     table->n = 0;
    454 
    455     for (i = 0; i < sources->n; i++) {
    456         pmSource *source = (pmSource *) sources->data[i];
    457         pmModel *model = pmModelSelect (source);
    458         if (model == NULL) continue;
    459 
    460         PAR = model->params->data.F32;
    461         dPAR = model->dparams->data.F32;
    462 
    463         dmag = dPAR[1] / PAR[1];
    464         type = pmSourceDophotType (source);
    465         lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]);
    466 
    467         row = psMetadataAlloc ();
    468         psMetadataAdd (row, PS_LIST_TAIL, "X_PIX",   PS_DATA_F32, "", PAR[2]);
    469         psMetadataAdd (row, PS_LIST_TAIL, "Y_PIX",   PS_DATA_F32, "", PAR[3]);
    470         psMetadataAdd (row, PS_LIST_TAIL, "MAG_RAW", PS_DATA_F32, "", source->fitMag + ZERO_POINT);
    471         psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", (int)(1000*dmag));
    472         psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", 32.0);
    473         psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP",  PS_DATA_F32, "", source->apMag + ZERO_POINT);
    474         psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", lsky);
    475         psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X",  PS_DATA_F32, "", type);
    476         psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y",  PS_DATA_F32, "", PAR[4]);
    477         psMetadataAdd (row, PS_LIST_TAIL, "THETA",   PS_DATA_F32, "", PAR[5]);
    478         psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT",  PS_DATA_STRING, "", "0");
    479         psMetadataAdd (row, PS_LIST_TAIL, "DUMMY",   PS_DATA_STRING, "", "123");
    480    
    481         psArrayAdd (table, 100, row);
    482         psFree (row);
    483     }
    484 
    485     // write the correct number of sources in the header
    486     psMetadataAdd (header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", table->n);
    487     psMetadataAdd (header, PS_LIST_TAIL, "EXTEND", PS_DATA_BOOL | PS_META_REPLACE, "this file has extensions", true);
    488 
    489     // set NAXIS to 0 (we don't write out the data array)
    490     mdi = psMetadataLookup (header, "NAXIS");
    491     mdi->data.S32 = 0;
    492     mdi->type = PS_DATA_S32;
    493 
    494     // create the basic table header
    495     theader = psMetadataAlloc ();
    496     psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE");
    497 
    498     // write out the header and table
    499     psFits *fits = psFitsOpen (filename, "w");
    500     psFitsWriteHeader (header, fits);
    501     psFitsWriteTable (fits, theader, table);
    502     psFitsClose (fits);
    503 
    504     return true;
    505 }
    506 
    507 /***** Text Output Methods *****/
    508 
    509 bool pmSourcesWriteText (psArray *sources, char *filename) {
    510 
    511     char *name = (char *) psAlloc (strlen(filename) + 10);
    512 
    513     sprintf (name, "%s.psf.dat", filename);
    514     pmModelWritePSFs (sources, name);
    515 
    516     sprintf (name, "%s.ext.dat", filename);
    517     pmModelWriteEXTs (sources, name);
    518 
    519     sprintf (name, "%s.nul.dat", filename);
    520     pmModelWriteNULLs (sources, name);
    521 
    522     sprintf (name, "%s.mnt.dat", filename);
    523     pmMomentsWriteText (sources, name);
    524 
    525     psFree (name);
    526     return true;
    527 }
    528 
    529 // write the PSF sources to an output file
    530 bool pmModelWritePSFs (psArray *sources, char *filename) {
    531 
    532     double dPos, dMag;
    533     int i, j;
    534     FILE *f;
    535     psF32 *PAR, *dPAR;
    536     pmModel  *model;
    537 
    538     f = fopen (filename, "w");
    539     if (f == NULL) {
    540         psLogMsg ("pmModelWritePSFs", 3, "can't open output file for moments%s\n", filename);
    541         return false;
    542     }
    543 
    544     // write sources with models first
    545     for (i = 0; i < sources->n; i++) {
    546         pmSource *source = (pmSource *) sources->data[i];
    547         if (source->type != PM_SOURCE_STAR) continue;
    548         model = source->modelPSF;
    549         if (model == NULL) continue;
    550 
    551         PAR  = model->params->data.F32;
    552         dPAR = model->dparams->data.F32;
    553        
    554         // dPos is positional error, dMag is mag error
    555         dPos = hypot (dPAR[2], dPAR[3]);
    556         dMag = dPAR[1] / PAR[1];
    557 
    558         fprintf (f, "%7.1f %7.1f  %7.1f %8.4f  %7.4f %7.4f  ",
    559                  PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos);
    560 
    561         for (j = 4; j < model->params->n; j++) {
    562             fprintf (f, "%9.6f ", PAR[j]);
    563         }
    564         fprintf (f, " : ");
    565         for (j = 4; j < model->params->n; j++) {
    566             fprintf (f, "%9.6f ", dPAR[j]);
    567         }
    568         fprintf (f, ": %8.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",
    569                  source[0].apMag, source[0].type, source[0].mode,
    570                  log10(model[0].chisq/model[0].nDOF),
    571                  source[0].moments->SN,
    572                  model[0].radius,
    573                  model[0].nDOF,
    574                  model[0].nIter);
    575     }
    576     fclose (f);
    577     return true;
    578 }
    579 
    580 // dump the sources to an output file
    581 bool pmModelWriteEXTs (psArray *sources, char *filename) {
    582 
    583     double dPos, dMag;
    584     int i, j;
    585     FILE *f;
    586     psF32 *PAR, *dPAR;
    587     pmModel  *model;
    588 
    589     f = fopen (filename, "w");
    590     if (f == NULL) {
    591         psLogMsg ("pmModelWriteEXTs", 3, "can't open output file for moments%s\n", filename);
    592         return false;
    593     }
    594 
    595     // write sources with models first
    596     for (i = 0; i < sources->n; i++) {
    597         pmSource *source = (pmSource *) sources->data[i];
    598 
    599         if (source->type != PM_SOURCE_EXTENDED) continue;
    600         model = pmSourceMagnitudes (source, NULL, 0.0);
    601         if (model == NULL) continue;
    602 
    603         PAR  = model->params->data.F32;
    604         dPAR = model->dparams->data.F32;
    605        
    606         // dPos is shape error
    607         // XXX these are hardwired for SGAUSS
    608         dPos = hypot ((dPAR[4] / PAR[4]), (dPAR[5] / PAR[5]));
    609         dMag = dPAR[1] / PAR[1];
    610 
    611         fprintf (f, "%7.1f %7.1f  %7.1f %8.4f  %7.4f %7.4f  ",
    612                  PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos);
    613 
    614         for (j = 4; j < model->params->n; j++) {
    615             fprintf (f, "%9.6f ", PAR[j]);
    616         }
    617         fprintf (f, " : ");
    618         for (j = 4; j < model->params->n; j++) {
    619             fprintf (f, "%9.6f ", dPAR[j]);
    620         }
    621         fprintf (f, ": %7.4f  %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",
    622                  source->apMag,
    623                  source[0].type, source[0].mode,
    624                  log10(model[0].chisq/model[0].nDOF),
    625                  source[0].moments->SN,
    626                  model[0].radius,
    627                  model[0].nDOF,
    628                  model[0].nIter);
    629     }
    630     fclose (f);
    631     return true;
    632 }
    633 
    634 // dump the sources to an output file
    635 bool pmModelWriteNULLs (psArray *sources, char *filename)
    636 {
    637 
    638     int i;
    639     FILE *f;
    640     pmMoments *moment = NULL;
    641     pmSource *source = NULL;
    642 
    643     f = fopen (filename, "w");
    644     if (f == NULL) {
    645         psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename);
    646         return false;
    647     }
    648 
    649     pmMoments *empty = pmMomentsAlloc ();
    650 
    651     // write sources with models first
    652     for (i = 0; i < sources->n; i++) {
    653         source = sources->data[i];
    654 
    655         // skip these sources (in PSF or EXT)
    656         if (source->type == PM_SOURCE_STAR) continue;
    657         if (source->type == PM_SOURCE_EXTENDED) continue;
    658            
    659         if (source->moments == NULL) {
    660             moment = empty;
    661         } else {
    662             moment = source->moments;
    663         }
    664 
    665         fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n",
    666                  source->peak->x, source->peak->y, source->peak->counts,
    667                  source->moments->x, source->moments->y,
    668                  source->moments->Sx, source->moments->Sy,
    669                  source->moments->Sum, source->moments->Peak,
    670                  source->moments->Sky, source->moments->SN,
    671                  source->moments->nPixels, source->type);
    672     }
    673     fclose (f);
    674     psFree (empty);
    675     return true;
    676 }
    677 
    678 // write the moments to an output file
    679 bool pmMomentsWriteText (psArray *sources, char *filename)
    680 {
    681 
    682     int i;
    683     FILE *f;
    684     pmSource *source = NULL;
    685 
    686     f = fopen (filename, "w");
    687     if (f == NULL) {
    688         psLogMsg ("pmMomentsWriteText", 3, "can't open output file for moments%s\n", filename);
    689         return false;
    690     }
    691 
    692     for (i = 0; i < sources->n; i++) {
    693         source = sources->data[i];
    694         if (source->moments == NULL)
    695             continue;
    696         fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %10.1f %7.1f %7.1f %7.1f  %4d %2d %#5x\n",
    697                  source->peak->x, source->peak->y, source->peak->counts,
    698                  source->moments->x, source->moments->y,
    699                  source->moments->Sx, source->moments->Sy,
    700                  source->moments->Sum, source->moments->Peak,
    701                  source->moments->Sky, source->moments->SN,
    702                  source->moments->nPixels, source->type, source->mode);
    703     }
    704     fclose (f);
    705     return true;
    706 }
    707 
    708 // write the peaks to an output file
    709 bool pmPeaksWriteText (psArray *peaks, char *filename) {
    710 
    711     int i;
    712     FILE *f;
    713        
    714     f = fopen (filename, "w");
    715     if (f == NULL) {
    716         psLogMsg ("pmPeaksWriteText", 3, "can't open output file for peaks%s\n", filename);
    717         return false;
    718     }
    719 
    720     for (i = 0; i < peaks->n; i++) {
    721         pmPeak *peak = peaks->data[i];
    722         if (peak == NULL) continue;
    723         fprintf (f, "%5d %5d  %7.1f\n",
    724                  peak->x, peak->y, peak->counts);
    725     }
    726     fclose (f);
    727     return true;
    728 }
    729 
    730 // translations between psphot object types and dophot object types
    731 int pmSourceDophotType (pmSource *source) {
    732 
    733     switch (source->type) {
    734 
    735       case PM_SOURCE_DEFECT:
    736       case PM_SOURCE_SATURATED:
    737         return (8);
    738 
    739       case PM_SOURCE_STAR:
    740         if (source->mode & PM_SOURCE_SATSTAR) return (10);
    741         if (source->mode & PM_SOURCE_POOR) return (7);
    742         if (source->mode & PM_SOURCE_FAIL) return (4);
    743         return (1);
    744 
    745       case PM_SOURCE_EXTENDED:
    746         return (2);
    747 
    748       default:
    749         return (0);
    750     }
    751     return (0);
    75221}
    75322
Note: See TracChangeset for help on using the changeset viewer.