IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 6571


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

Location:
trunk/psphot/src
Files:
2 added
5 deleted
18 edited

Legend:

Unmodified
Added
Removed
  • trunk/psphot/src/Makefile.am

    r6481 r6571  
    1212
    1313psphot_SOURCES =                \
     14        psphot.c                \
     15        psphotModelGroupInit.c  \
     16        psphotArguments.c       \
     17        psphotParseCamera.c     \
     18        testPSphotLoop.c        \
     19        psphotMaskCell.c        \
    1420        psphotReadout.c         \
    15         psphotImageStats.c      \
    16         psphotImageBackground.c \
     21        psphotImageMedian.c     \
    1722        psphotFindPeaks.c       \
    1823        psphotSourceStats.c     \
     
    2429        psphotReplaceUnfit.c    \
    2530        psphotApResid.c         \
    26         psphotOutput.c          \
    27         psphotLoadPixels.c      \
    28         psphotModelGroupInit.c  \
    29         psphotGrowthCurve.c     \
    30         psphotCleanup.c         \
     31        psphotMagnitudes.c      \
     32        psphotSkyReplace.c      \
    3133        psphotEvalPSF.c         \
    3234        psphotEvalFLT.c         \
    33         psphotSourceFits.c      \
     35        psphotSourceFits.c     
     36
     37junk =  psphotGrowthCurve.c     \
     38        psphotCleanup.c         \
    3439        psphotSortBySN.c        \
    3540        psphotDefinePixels.c    \
    36         psphotMagnitudes.c      \
    3741        psphotRadiusChecks.c    \
    38         psphotImageMedian.c     \
    39         psphotSkyReplace.c      \
    40         psLine.c                \
    41         psSparse.c              \
    42         psPolynomialUtils.c     \
    4342        psModulesUtils.c        \
    4443        pmSourceContour.c       \
    4544        pmModelFitSet.c         \
    46         pmCellSetMask.c         \
    4745        pmSourceFitSet.c        \
    48         psphot.c                \
    49         psphotArguments.c       \
    50         psphotParseCamera.c     \
    51         psphotModelTest.c       \
    52         psphotImageLoop.c       
     46        psphotModelTest.c       
    5347
    5448noinst_HEADERS =                \
     
    6862#       psphotApplyPSF.c       
    6963#       psphotFitGalaxies.c
     64#       psphotImageStats.c     
     65#       psphotImageBackground.c
     66#       psphotOutput.c         
     67#       psphotLoadPixels.c     
  • trunk/psphot/src/psModulesUtils.c

    r6522 r6571  
    11# include "psModulesUtils.h"
    2 
    3 // extract config informatin from config data or from image header, if specified
    4 // XXX EAM : this function should be replaced with Paul's image/header/metadata load scheme
    5 psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) {
    6 
    7     char *source;
    8     char *keyword;
    9     psF32 value;
    10     psMetadataItem *item;
    11 
    12     // look for the entry in the config collection
    13     item = psMetadataLookup (config, name);
    14     if (item == NULL) {
    15         psLogMsg ("pmConfigLookupF32", 2, "no key %s in config data, trying as header keyword", name);
    16         value = psMetadataLookupF32 (status, header, name);
    17         if (*status == false) {
    18             psAbort ("pmConfigLookupF32", "no key %s in header", name);
    19         }
    20         *status = true;
    21         return (value);
    22     }   
    23 
    24     // I'm either expecting a string, with the name "HD:keyword"...
    25     if (item->type == PS_DATA_STRING) {
    26         source = item->data.V;
    27         if (!strncasecmp (source, "HD:", 3)) {
    28             keyword = &source[3];
    29             value = psMetadataLookupF32 (status, header, keyword);
    30             if (*status == false) {
    31                 psAbort ("pmConfigLookupF32", "no key %s in config", name);
    32             }
    33             *status = true;
    34             // psFree (item);
    35             return (value);
    36         }       
    37     }
    38 
    39     //  ... or a value (F32?)
    40     if (item->type == PS_DATA_F32) {
    41         value = item->data.F32;
    42         // psFree (item);
    43         return (value);
    44     }
    45 
    46     *status = false;
    47     psError(PS_ERR_UNKNOWN, true, "invalid item");
    48     return (0);
    49 }
    50 
    51 bool pmModelFitStatus (pmModel *model) {
    52 
    53     bool status;
    54 
    55     pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type);
    56     status = statusFunc (model);
    57 
    58     return (status);
    59 }
    60 
    61 pmModel *pmModelSelect (pmSource *source) {
    62     switch (source->type) {
    63       case PM_SOURCE_STAR:
    64         return source->modelPSF;
    65        
    66       case PM_SOURCE_EXTENDED:
    67         return source->modelEXT;
    68        
    69       default:
    70         return NULL;
    71     }
    72     return NULL;
    73 }
    74 
    75 bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) {
    76 
    77     float obsSum = 0;
    78     float fitSum = 0;
    79     float sky = model->params->data.F32[0];
    80 
    81     *fitMag = 99.0;
    82     *obsMag = 99.0;
    83 
    84     pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
    85     fitSum = modelFluxFunc (model->params);
    86     if (fitSum <= 0) return false;
    87     if (!isfinite(fitSum)) return false;
    88     *fitMag = -2.5*log10(fitSum);
    89 
    90     for (int ix = 0; ix < image->numCols; ix++) {
    91         for (int iy = 0; iy < image->numRows; iy++) {
    92             if (mask->data.U8[iy][ix]) continue;
    93             obsSum += image->data.F32[iy][ix] - sky;
    94         }
    95     }
    96     if (obsSum <= 0) return false;
    97     *obsMag = -2.5*log10(obsSum);
    98 
    99     return (true);
    100 }
    101 
    102 // force the fitted sky to meet the source flux at outer radius
    103 bool pmModelSkyOffset (pmModel *model, float x, float y, float radius) {
    104 
    105     float flux;
    106 
    107     psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    108 
    109     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    110     psVector *params = model->params;
    111  
    112     coord->data.F32[0] = x + radius;
    113     coord->data.F32[1] = y;
    114     flux = modelFunc (NULL, params, coord);
    115     params->data.F32[0] = flux;
    116 
    117     psFree (coord);
    118 
    119     return true;
    120 }
    121 
    122 pmModel *pmModelCopy (pmModel *model) {
    123 
    124     pmModel *new = pmModelAlloc (model->type);
    125    
    126     new->chisq  = model->chisq;
    127     new->nDOF   = model->nDOF;
    128     new->nIter  = model->nIter;
    129     new->status = model->status;
    130     new->radius = model->radius;
    131 
    132     for (int i = 0; i < new->params->n; i++) {
    133         new->params->data.F32[i]  = model->params->data.F32[i];
    134         new->dparams->data.F32[i] = model->dparams->data.F32[i];
    135     }
    136    
    137     return (new);
    138 }
    139 
    140 float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj) {
    141 
    142     int Xs, Xe, Ys, Ye;
    143     int xi, xj, yi, yj;
    144     int xIs, xJs, yIs, yJs;
    145     int xIe, yIe;
    146     float flux, wt;
    147 
    148     psImage *Pi = Mi->pixels;
    149     psImage *Pj = Mj->pixels;
    150 
    151     psImage *Wi = Mi->weight;
    152 
    153     psImage *Ti = Mi->mask;
    154     psImage *Tj = Mj->mask;
    155 
    156     Xs = PS_MAX (Pi->col0, Pj->col0);
    157     Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
    158    
    159     Ys = PS_MAX (Pi->row0, Pj->row0);
    160     Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
    161    
    162     xIs = Xs - Pi->col0;
    163     xJs = Xs - Pj->col0;
    164     yIs = Ys - Pi->row0;
    165     yJs = Ys - Pj->row0;
    166 
    167     xIe = Xe - Pi->col0;
    168     yIe = Ye - Pi->row0;
    169 
    170     // note that this is addressing the same image pixels,
    171     // though only if both are source not model images
    172     flux = 0;
    173     for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    174         for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
    175             if (Ti->data.U8[yi][xi]) continue;
    176             if (Tj->data.U8[yj][xj]) continue;
    177             wt = Wi->data.F32[yi][xi];
    178             if (wt > 0) {
    179                 flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
    180             }
    181         }
    182     }
    183     return (flux);
    184 }
    185 
    186 float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj) {
    187 
    188     int Xs, Xe, Ys, Ye;
    189     int xi, xj, yi, yj;
    190     int xIs, xJs, yIs, yJs;
    191     int xIe, yIe;
    192     float flux, wt;
    193 
    194     psImage *Pi = Mi->pixels;
    195     psImage *Pj = Mj->pixels;
    196 
    197     psImage *Wi = Mi->weight;
    198 
    199     psImage *Ti = Mi->mask;
    200     psImage *Tj = Mj->mask;
    201 
    202     Xs = PS_MAX (Pi->col0, Pj->col0);
    203     Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
    204    
    205     Ys = PS_MAX (Pi->row0, Pj->row0);
    206     Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
    207    
    208     xIs = Xs - Pi->col0;
    209     xJs = Xs - Pj->col0;
    210     yIs = Ys - Pi->row0;
    211     yJs = Ys - Pj->row0;
    212 
    213     xIe = Xe - Pi->col0;
    214     yIe = Ye - Pi->row0;
    215 
    216     // note that this is addressing the same image pixels,
    217     // though only if both are source not model images
    218     flux = 0;
    219     for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    220         for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
    221             if (Ti->data.U8[yi][xi]) continue;
    222             if (Tj->data.U8[yj][xj]) continue;
    223             wt = Wi->data.F32[yi][xi];
    224             if (wt > 0) {
    225                 flux += 1.0 / wt;
    226             }
    227         }
    228     }
    229     return (flux);
    230 }
    2312
    2323static void ppConfigFree (ppConfig *config) {
  • trunk/psphot/src/psModulesUtils.h

    r6522 r6571  
    2323    psMetadata *recipe;                 // The recipe (i.e., specific setups)
    2424    psMetadata *arguments;              // Command-line arguments
    25     psMetadata *options;                // Command-line recipe options
    2625    psDB       *database;               // Database handle
    2726} ppConfig;
  • trunk/psphot/src/psphot.c

    r6481 r6571  
    1010
    1111    // load command-line arguments, options, and system config data
    12     ppConfig *config = psphotArguments (&argc, argv);
     12    pmConfig *config = psphotArguments (&argc, argv);
    1313
    1414    // load input data (config and images (signal, noise, mask)
    15     ppFile *input = psphotParseCamera (config);
    16 
    17     // check on output mode,format, setup basic filename, etc
    18     psphotOutputPrep (input, config);
     15    psphotParseCamera (config);
    1916
    2017    // call psphot for each readout
    21     psphotImageLoop (input, config);
     18    psphotImageLoop (config);
    2219
    2320    psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete"));
    2421
    25     psFree (input);
    2622    psFree (config);
    2723    psphotCleanup ();
  • trunk/psphot/src/psphot.h

    r6481 r6571  
    1010# include <pmPSFtry.h>
    1111# include <pmModelGroup.h>
    12 # include "psLibUtils.h"
    13 # include "psModulesUtils.h"
    14 # include "psSparse.h"
     12// # include "psLibUtils.h"
     13// # include "psModulesUtils.h"
     14// # include "psSparse.h"
    1515# include "pmFPAConstruct.h"
    1616# include "pmConcepts.h"
     
    2121
    2222// top-level psphot functions
    23 ppConfig       *psphotArguments (int *argc, char **argv);
    24 ppFile         *psphotParseCamera (ppConfig *config);
    25 bool            psphotImageLoop (ppFile *input, ppConfig *config);
     23pmConfig       *psphotArguments (int *argc, char **argv);
     24bool            psphotParseCamera (pmConfig *config);
     25bool            psphotImageLoop (pmConfig *config);
    2626
    2727bool            psphotModelTest (pmReadout *readout, psMetadata *arguments, psMetadata *recipe);
    2828bool            psphotReadout (pmReadout *readout, psMetadata *config);
    29 bool            ppImageLoadPixels (ppFile *input, psDB *db, int chipNum, int cellNum);
    3029void            psphotCleanup (void);
     30
     31// bool            ppImageLoadPixels (ppFile *input, psDB *db, int chipNum, int cellNum);
    3132
    3233// psphotReadout functions
     
    4748bool            psphotApResid (pmReadout *readout, psArray *sources, psMetadata *config, pmPSF *psf);
    4849
    49 // XXX deprecate these?
    50 // bool            psphotFullFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
    51 // bool            psphotApplyPSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
    52 // bool            psphotFitExtended (pmReadout *readout, psMetadata *config, psArray *sources, psStats *skyStats);
    53 
    5450// basic support functions
    55 pmModel        *pmModelCopy (pmModel *model);
    56 pmModel        *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius);
    57 float           pmSourceCrossProduct (pmSource *Mi, pmSource *Mj);
    58 float           pmSourceCrossWeight (pmSource *Mi, pmSource *Mj);
    59 psArray        *pmSourceContour_EAM (psImage *image, int x, int y, float threshold);
    6051void            psphotModelGroupInit (void);
    6152int             psphotSortBySN (const void **a, const void **b);
     
    6354bool            psphotGrowthCurve (pmReadout *readout, pmPSF *psf);
    6455void            psphotTestArguments (int *argc, char **argv);
    65 bool            pmCellSetMask (pmCell *cell, psMetadata *recipe);
     56bool            psphotMaskCell (pmCell *cell, psMetadata *recipe);
    6657bool            psphotBackgroundNames (psMetadata *arguments);
    6758bool            psphotSkyReplace (pmReadout *readout, psImage *background);
     
    7364bool            psphotInitRadiusEXT (psMetadata *config, pmModelType type);
    7465bool            psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model);
    75 bool            psphotDefinePixels (pmSource *mySource, const pmReadout *readout, psF32 x, psF32 y, psF32 Radius);
    76 bool            psphotRedefinePixels (pmSource *mySource, const pmReadout *readout, psF32 x, psF32 y, psF32 Radius);
    7766
    7867// output functions
     68# if (0)
    7969bool            pmSourcesWriteSX   (psArray *sources, char *filename);
    8070bool            pmSourcesWriteOBJ  (psArray *sources, char *filename);
     
    9383int             psphotSaveImage (psMetadata *header, psImage *image, char *filename);
    9484bool            psphotUpdateHeader (psMetadata *header, psMetadata *config);
    95 int             pmSourcesDophotType (pmSource *source);
    96 bool            psMetadataItemTransfer (psMetadata *out, psMetadata *in, char *key);
     85# endif
    9786
    98 char           *psphotSplitName (psMetadata *header);
    99 void            psphotOutputPrep (ppFile *file, ppConfig *config);
    100 void            psphotOutputCleanup (void);
    101 char           *psphotNameSubstitute (char *input, char *replace, char *key);
     87// void            psphotOutputPrep (ppFile *file, ppConfig *config);
     88// void            psphotOutputCleanup (void);
    10289
    10390bool            psphotMagnitudes (psMetadata *config, psArray *sources, pmPSF *psf);
     
    10794bool            psphotEvalDBL (pmSource *source, pmModel *model);
    10895bool            psphotEvalEXT (pmSource *source, pmModel *model);
    109 
    110 // functions to support simultaneous multi-source fitting
    111 bool            psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF);
    112 bool            pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF);
    113 psF32           pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x);
    114 bool            pmModelFitSetInit (pmModelType type);
    115 void            pmModelFitSetClear (void);
    11696
    11797//  functions to support the source fitting process
     
    124104psArray        *psphotFitDBL (pmReadout *readout, pmSource *source);
    125105
    126 // bicubic interpolation
    127 psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y);
    128 psPlane         psImageBicubeMin (psPolynomial2D *poly);
    129 
    130 bool psImageJpegColormap (char *name);
    131 bool psImageJpeg (psImage *image, char *filename, float zero, float scale);
     106// XXX these can probably be dropped:
    132107
    133108// optional mode for clip fit?
     
    143118    const psVector *z,
    144119    const psVector *t);
     120
     121// functions to support simultaneous multi-source fitting
     122bool            psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF);
     123bool            pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF);
     124psF32           pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x);
     125bool            pmModelFitSetInit (pmModelType type);
     126void            pmModelFitSetClear (void);
     127
     128// bicubic interpolation
     129psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y);
     130psPlane         psImageBicubeMin (psPolynomial2D *poly);
     131
     132bool psImageJpegColormap (char *name);
     133bool psImageJpeg (psImage *image, char *filename, float zero, float scale);
     134
     135
     136// XXX deprecate
     137
     138// char           *psphotNameSubstitute (char *input, char *replace, char *key);
     139// char           *psphotSplitName (psMetadata *header);
     140// pmModel        *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius);
     141// float           pmSourceCrossProduct (pmSource *Mi, pmSource *Mj);
     142// float           pmSourceCrossWeight (pmSource *Mi, pmSource *Mj);
     143// psArray        *pmSourceContour_EAM (psImage *image, int x, int y, float threshold);
  • trunk/psphot/src/psphotArguments.c

    r6522 r6571  
    11# include "psphot.h"
     2# include <glob.h>
    23
    34static void usage (void);
    45
    5 ppConfig *psphotArguments (int *argc, char **argv) {
     6pmConfig *psphotArguments (int *argc, char **argv) {
    67
    78    int N;
     
    1213    psLogSetFormat ("M");
    1314
    14     ppConfig *config = ppConfigAlloc ();
     15    // these other options override the PSPHOT recipe options
     16    psMetadata *options = psMetadataAlloc ();
    1517
    16     // load config data from default locations
    17     if (!pmConfigRead(&config->site, &config->camera, &config->recipe, argc, argv, PSPHOT_RECIPE)) {
    18         psErrorStackPrint(stderr, "Can't find site configuration!\n");
    19         exit (EXIT_FAILURE);
     18    // run the test model (requires X,Y coordinate)
     19    if ((N = psArgumentGet (*argc, argv, "-modeltest"))) {
     20        psMetadataAddBool (options, PS_LIST_TAIL, "TEST_FIT",   0, "", true);
     21        psMetadataAddF32  (options, PS_LIST_TAIL, "TEST_FIT_X", 0, "", atof(argv[N+1]));
     22        psMetadataAddF32  (options, PS_LIST_TAIL, "TEST_FIT_Y", 0, "", atof(argv[N+2]));
     23
     24        psArgumentRemove (N, argc, argv);
     25        psArgumentRemove (N, argc, argv);
     26        psArgumentRemove (N, argc, argv);
     27
     28        // specify the modeltest model
     29        if ((N = psArgumentGet (*argc, argv, "-model"))) {
     30            psArgumentRemove (N, argc, argv);
     31            psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODEL", 0, "", argv[N]);
     32            psArgumentRemove (N, argc, argv);
     33        }
     34
     35        // specify the test fit mode
     36        if ((N = psArgumentGet (*argc, argv, "-fitmode"))) {
     37            psArgumentRemove (N, argc, argv);
     38            psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODE", 0, "", argv[N]);
     39            psArgumentRemove (N, argc, argv);
     40        }
     41        if ((N = psArgumentGet (*argc, argv, "-fitset"))) {
     42            psArgumentRemove (N, argc, argv);
     43            psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_SET", 0, "", argv[N]);
     44            psArgumentRemove (N, argc, argv);
     45        }
    2046    }
    2147
    22     // Parse other command-line arguments
    23     config->arguments = psMetadataAlloc ();
    24 
    25     // arguments (must be supplied for each run, or not used)
    26     // mask image (input) is used by psphotImageLoop
    27     if ((N = psArgumentGet (*argc, argv, "-mask"))) {
     48    // photcode : used in output to supplement header data (argument or recipe?)
     49    if ((N = psArgumentGet (*argc, argv, "-photcode"))) {
    2850        psArgumentRemove (N, argc, argv);
    29         psMetadataAddStr (config->arguments, PS_LIST_TAIL, "MASK_IMAGE", PS_META_REPLACE, "", argv[N]);
     51        psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", argv[N]);
    3052        psArgumentRemove (N, argc, argv);
    3153    }
    3254
    33     // weight image (input) is used by psphotImageLoop
    34     if ((N = psArgumentGet (*argc, argv, "-weight"))) {
     55    // break : used from recipe throughout psphotReadout
     56    if ((N = psArgumentGet (*argc, argv, "-break"))) {
    3557        psArgumentRemove (N, argc, argv);
    36         psMetadataAddStr (config->arguments, PS_LIST_TAIL, "WEIGHT_IMAGE", PS_META_REPLACE, "", argv[N]);
     58        psMetadataAddStr (options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]);
    3759        psArgumentRemove (N, argc, argv);
    3860    }
    3961
    40     // residual image (output) is used by psphotOutput
    41     if ((N = psArgumentGet (*argc, argv, "-resid"))) {
     62    // fitmode : used from recipe throughout psphotReadout
     63    if ((N = psArgumentGet (*argc, argv, "-fitmode"))) {
    4264        psArgumentRemove (N, argc, argv);
    43         psMetadataAddStr (config->arguments, PS_LIST_TAIL, "RESID_IMAGE", PS_META_REPLACE, "", argv[N]);
     65        psMetadataAddStr (options, PS_LIST_TAIL, "FITMODE", PS_META_REPLACE, "", argv[N]);
    4466        psArgumentRemove (N, argc, argv);
    4567    }
    4668
    47     // background image (output) is used by psphotOutput
    48     if ((N = psArgumentGet (*argc, argv, "-background"))) {
     69    // analysis region : overrides recipe value, used in psphotReadout/psphotEnsemblePSF
     70    if ((N = psArgumentGet (*argc, argv, "-region"))) {
    4971        psArgumentRemove (N, argc, argv);
    50         psMetadataAddStr (config->arguments, PS_LIST_TAIL, "BACKGROUND_IMAGE", PS_META_REPLACE, "", argv[N]);
     72        psMetadataAddStr (options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", argv[N]);
    5173        psArgumentRemove (N, argc, argv);
    5274    }
    5375
    54     // background model (output) is used by psphotOutput
    55     if ((N = psArgumentGet (*argc, argv, "-backmodel"))) {
     76    // other arbitrary recipe values: -D key value (all added as string)
     77    while ((N = psArgumentGet (*argc, argv, "-D"))) {
    5678        psArgumentRemove (N, argc, argv);
    57         psMetadataAddStr (config->arguments, PS_LIST_TAIL, "BACKGROUND_MODEL", PS_META_REPLACE, "", argv[N]);
     79        psMetadataAddStr (options, PS_LIST_TAIL, argv[N], 0, "", argv[N+1]);
     80        psArgumentRemove (N, argc, argv);
    5881        psArgumentRemove (N, argc, argv);
    5982    }
     83
     84    // other arbitrary recipe values: -Df key value (all added as float)
     85    while ((N = psArgumentGet (*argc, argv, "-Df"))) {
     86        psArgumentRemove (N, argc, argv);
     87        psMetadataAddF32 (options, PS_LIST_TAIL, argv[N], 0, "", atof(argv[N+1]));
     88        psArgumentRemove (N, argc, argv);
     89        psArgumentRemove (N, argc, argv);
     90    }
     91
     92    // other arbitrary recipe values: -Di key value (all added as int)
     93    while ((N = psArgumentGet (*argc, argv, "-Di"))) {
     94        psArgumentRemove (N, argc, argv);
     95        psMetadataAddS32 (options, PS_LIST_TAIL, argv[N], 0, "", atoi(argv[N+1]));
     96        psArgumentRemove (N, argc, argv);
     97        psArgumentRemove (N, argc, argv);
     98    }
     99
     100    // load config data from default locations
     101    pmConfig *config = pmConfigRead(argc, argv, PSPHOT_RECIPE);
     102
     103    // Storage for other command-line arguments
     104    config->arguments = psMetadataAlloc ();
     105
     106    // save these recipe options until we have loaded the options
     107    psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "PSPHOT.OPTIONS",  PS_META_REPLACE, "", options);
    60108
    61109    // chip selection is used to limit chips to be processed
     
    66114    }
    67115
    68     // run the test model (requires X,Y coordinate)
    69     if ((N = psArgumentGet (*argc, argv, "-modeltest"))) {
    70         psMetadataAddBool (config->arguments, PS_LIST_TAIL, "TEST_FIT",   0, "", true);
    71         psMetadataAddF32  (config->arguments, PS_LIST_TAIL, "TEST_FIT_X", 0, "", atof(argv[N+1]));
    72         psMetadataAddF32  (config->arguments, PS_LIST_TAIL, "TEST_FIT_Y", 0, "", atof(argv[N+2]));
     116    // we load all input files onto a psArray, to be parsed later
     117    psArray *input = psArrayAlloc (16);
    73118
     119    // load the list of filenames the supplied file (may be a glob: "file*.fits")
     120    if ((N = psArgumentGet (*argc, argv, "-file"))) {
     121        glob_t globList;
    74122        psArgumentRemove (N, argc, argv);
    75         psArgumentRemove (N, argc, argv);
    76         psArgumentRemove (N, argc, argv);
    77 
    78         // specify the modeltest model
    79         if ((N = psArgumentGet (*argc, argv, "-model"))) {
    80             psArgumentRemove (N, argc, argv);
    81             psMetadataAddStr (config->arguments, PS_LIST_TAIL, "TEST_FIT_MODEL", 0, "", argv[N]);
    82             psArgumentRemove (N, argc, argv);
     123        globList.gl_offs = 0;
     124        glob (argv[N], 0, NULL, &globList);
     125        for (int i = 0; i < globList.gl_pathc; i++) {
     126            char *filename = psStringCopy (globList.gl_pathv[i]);
     127            psArrayAdd (input, 16, filename);
    83128        }
    84 
    85         // specify the test fit mode
    86         if ((N = psArgumentGet (*argc, argv, "-fitmode"))) {
    87             psArgumentRemove (N, argc, argv);
    88             psMetadataAddStr (config->arguments, PS_LIST_TAIL, "TEST_FIT_MODE", 0, "", argv[N]);
    89             psArgumentRemove (N, argc, argv);
    90         }
    91         if ((N = psArgumentGet (*argc, argv, "-fitset"))) {
    92             psArgumentRemove (N, argc, argv);
    93             psMetadataAddStr (config->arguments, PS_LIST_TAIL, "TEST_FIT_SET", 0, "", argv[N]);
    94             psArgumentRemove (N, argc, argv);
    95         }
    96     }
    97 
    98     // these other options override the recipe options
    99     config->options = psMetadataAlloc ();
    100 
    101     // PSF : optional psf file to be loaded
    102     if ((N = psArgumentGet (*argc, argv, "-psf"))) {
    103         psArgumentRemove (N, argc, argv);
    104         psMetadataAddStr (config->options, PS_LIST_TAIL, "PSF_INPUT_FILE", PS_META_REPLACE, "", argv[N]);
    105129        psArgumentRemove (N, argc, argv);
    106130    }
    107131
    108     // PSF : optional psf file to be loaded
    109     if ((N = psArgumentGet (*argc, argv, "-psfout"))) {
     132    // load the list from the supplied text file
     133    if ((N = psArgumentGet (*argc, argv, "-list"))) {
     134        int nItems;
     135        char line[1024]; // XXX limits the list lines to 1024 chars
     136        char word[1024];
     137        char *filename;
     138
    110139        psArgumentRemove (N, argc, argv);
    111         psMetadataAddStr (config->options, PS_LIST_TAIL, "PSF_OUTPUT_FILE", PS_META_REPLACE, "", argv[N]);
     140        FILE *f = fopen (argv[N], "r");
     141        if (f == NULL) {
     142            psAbort ("psphot", "unable to open specified list file");
     143        }
     144        while (fgets (line, 1024, f) != NULL) {
     145            nItems = sscanf (line, "%s", word);
     146            switch (nItems) {
     147              case 0:
     148                break;
     149              case 1:
     150                filename = psStringCopy (word);
     151                psArrayAdd (input, 16, filename);
     152                break;
     153              default:
     154                // rigid format, no comments allowed?
     155                psAbort ("psphot", "error parsing input list file");
     156                break;
     157            }
     158        }
    112159        psArgumentRemove (N, argc, argv);
    113     }
     160    }           
     161    if (input->n < 1) usage ();
    114162
    115     // photcode : used in output to supplement header data (argument or recipe?)
    116     if ((N = psArgumentGet (*argc, argv, "-photcode"))) {
    117         psArgumentRemove (N, argc, argv);
    118         psMetadataAddStr (config->options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", argv[N]);
    119         psArgumentRemove (N, argc, argv);
    120     }
     163    // input list gets places as an array on the config->arguements list
     164    psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "INPUT",  PS_META_REPLACE, "", input);
    121165
    122     // break : used from recipe throughout psphotReadout
    123     if ((N = psArgumentGet (*argc, argv, "-break"))) {
    124         psArgumentRemove (N, argc, argv);
    125         psMetadataAddStr (config->options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]);
    126         psArgumentRemove (N, argc, argv);
    127     }
    128 
    129     // fitmode : used from recipe throughout psphotReadout
    130     if ((N = psArgumentGet (*argc, argv, "-fitmode"))) {
    131         psArgumentRemove (N, argc, argv);
    132         psMetadataAddStr (config->options, PS_LIST_TAIL, "FITMODE", PS_META_REPLACE, "", argv[N]);
    133         psArgumentRemove (N, argc, argv);
    134     }
    135 
    136     // analysis region : overrides recipe value, used in psphotReadout/psphotEnsemblePSF
    137     if ((N = psArgumentGet (*argc, argv, "-region"))) {
    138         psArgumentRemove (N, argc, argv);
    139         psMetadataAddStr (config->options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", argv[N]);
    140         psArgumentRemove (N, argc, argv);
    141     }
    142 
    143     // other arbitrary recipe values: -D key value (all added as string)
    144     while ((N = psArgumentGet (*argc, argv, "-D"))) {
    145         psArgumentRemove (N, argc, argv);
    146         psMetadataAddStr (config->options, PS_LIST_TAIL, argv[N], 0, "", argv[N+1]);
    147         psArgumentRemove (N, argc, argv);
    148         psArgumentRemove (N, argc, argv);
    149     }
    150 
    151     // other arbitrary recipe values: -Df key value (all added as float)
    152     while ((N = psArgumentGet (*argc, argv, "-Df"))) {
    153         psArgumentRemove (N, argc, argv);
    154         psMetadataAddF32 (config->options, PS_LIST_TAIL, argv[N], 0, "", atof(argv[N+1]));
    155         psArgumentRemove (N, argc, argv);
    156         psArgumentRemove (N, argc, argv);
    157     }
    158 
    159     // other arbitrary recipe values: -Di key value (all added as int)
    160     while ((N = psArgumentGet (*argc, argv, "-Di"))) {
    161         psArgumentRemove (N, argc, argv);
    162         psMetadataAddS32 (config->options, PS_LIST_TAIL, argv[N], 0, "", atoi(argv[N+1]));
    163         psArgumentRemove (N, argc, argv);
    164         psArgumentRemove (N, argc, argv);
    165     }
    166 
    167     if (*argc != 3) usage ();
     166    if (*argc != 2) usage ();
    168167
    169168    // input and output positions are fixed
    170     psMetadataAddStr (config->arguments, PS_LIST_TAIL, "INPUT_FILE",  PS_META_REPLACE, "", argv[1]);
    171     psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT_ROOT", PS_META_REPLACE, "", argv[2]);
     169    psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", PS_META_REPLACE, "", argv[1]);
    172170
    173171    psTrace(__func__, 1, "Done with psphotArguments...\n");
     
    176174
    177175static void usage (void) {
    178     fprintf (stderr, "USAGE: psphot (image) (output)\n");
     176    fprintf (stderr, "USAGE: psphot [-file image(s)] [-list imagelist] (output)\n");
    179177    exit (2);
    180178}
  • trunk/psphot/src/psphotChoosePSF.c

    r6481 r6571  
    33// 2006.02.07 : no leaks!
    44// try PSF models and select best option
    5 pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources) {
     5pmPSF *psphotChoosePSF (psArray *sources, psMetadata *config) {
    66
    77    bool            status;
  • trunk/psphot/src/psphotCleanup.c

    r6399 r6571  
    99    psTimeFinalize ();
    1010    pmConceptsDone ();
    11     fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot");
     11    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot");
    1212    return;
    1313}
  • trunk/psphot/src/psphotEnsemblePSF.c

    r6481 r6571  
    44// 2006.02.07 : no leaks!
    55// fit all reasonable sources with the linear PSF model
    6 bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, bool final) {
     6bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *config, pmPSF *psf, bool final) {
    77
    88    bool  status;
  • trunk/psphot/src/psphotImageMedian.c

    r6481 r6571  
    99static int MAX_SAMPLE_PIXELS;
    1010
    11 static char *backImage = NULL;
    12 static char *backModel = NULL;
    13 
    14 bool psphotBackgroundNames (psMetadata *arguments) {
    15 
    16     bool status;
    17 
    18     backImage = psMetadataLookupStr (&status, arguments, "BACKGROUND_IMAGE");
    19     backModel = psMetadataLookupStr (&status, arguments, "BACKGROUND_MODEL");
    20 
    21     return true;
    22 }
    23 
    2411// generate the median in NxN boxes, clipping heavily
    2512// linear interpolation to generate full-scale model
    26 psImage *psphotImageMedian (pmReadout *readout, psMetadata *config)
     13psImage *psphotImageMedian (psMetadata *config, pmFPAview *view)
    2714{
    2815    bool status;
    2916    psRegion region;
     17    psImage *model = NULL;
     18
     19    psTimerStart ("psphot");
     20
     21    // find the currently selected readout
     22    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
     23    pmFPA *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
     24    pmReadout  *readout = pmFPAviewThisReadout (view, input);
     25
     26    psImage *image = readout->image;
     27    psImage *mask  = readout->mask;
    3028
    3129    rnd = psRandomAlloc (PS_RANDOM_TAUS, 0);
    3230    MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, config, "IMSTATS_NPIX");
    3331    if (!status) MAX_SAMPLE_PIXELS = 1000;
    34 
    35     psTimerStart ("psphot");
    36 
    37     psImage *image = readout->image;
    38     psImage *mask  = readout->mask;
    3932
    4033    // dimensions of input & output image
     
    5649    int ny = (Ny % DY) ? (int)(Ny / DY) + 1 : Ny / DY;
    5750
    58     // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    59     psImage *model = psImageAlloc (nx, ny, PS_TYPE_F32);
     51    // select model pixels, from output background model file, or create
     52    model = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32);
    6053
    6154    // measure clipped median for subimages
     
    7669
    7770            // XXX psImageStats is still very inefficient and poorly coded...
     71            // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    7872            // stats = psImageStats (stats, subset, maskset, 0xff);
    79             // model->data.F32[iy][ix] = stats->clippedMean;
     73            // backMdl->data.F32[iy][ix] = stats->clippedMean;
    8074        }
    8175    }
     
    8377
    8478    // linear interpolation to full-scale
    85    
    86     psImage *background = psImageAlloc (Nx, Ny, PS_TYPE_F32);
     79
     80    // select background pixels, from output background file, or create
     81    background = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
    8782
    8883    // XXX this code skips the initial pixels
     
    216211    psLogMsg ("psphot", 3, "build resampled image: %f sec\n", psTimerMark ("psphot"));
    217212
     213    // back-sub image pixels, from output background file (don't create if not requested)
     214    backSub = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKSUB", 0, 0, PS_TYPE_F32);
     215
    218216    // subtract the background model
    219217    for (int j = 0; j < image->numRows; j++) {
     
    221219            if (!mask->data.U8[j][i]) {
    222220                image->data.F32[j][i] -= background->data.F32[j][i];
    223             }
    224         }
    225     }
    226 
    227     // optionally save background
    228     if (backImage != NULL) psphotSaveImage (NULL, background, backImage);
    229     if (backModel != NULL) psphotSaveImage (NULL, model, backModel);
     221                if (backSub) {
     222                    backSub->data.F32[j][i] = image->data.F32[j][i];
     223                }
     224            }
     225        }
     226    }
    230227
    231228    psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
    232229    psFree (rnd);
     230
     231    // it is safe to free all of these: either they are saved on config->files or are local
    233232    psFree (model);
    234     return (background);
     233    psFree (background);
     234    psFree (backSub);
     235    return true;
    235236}
    236237
  • trunk/psphot/src/psphotMagnitudes.c

    r6495 r6571  
    11# include "psphot.h"
    2 
    3 // XXX EAM : the apMag should only be calculated for the brighter sources?
    4 // XXX EAM : SN limit set by user?
    5 // XXX EAM : masked region should be (optionally) elliptical
    6 pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius) {
    7 
    8     int status;
    9     bool isPSF;
    10     float x, y;
    11     float rflux;
    12     float radius;
    13     pmModel *model;
    14 
    15     switch (source->type) {
    16       case PM_SOURCE_STAR:
    17         model = source->modelPSF;
    18         if (model == NULL) return NULL;
    19         radius = (apRadius > 0) ? apRadius : model->radius;
    20         isPSF = true;
    21         break;
    22 
    23       case PM_SOURCE_EXTENDED:
    24         model = source->modelEXT;
    25         if (model == NULL) return NULL;
    26         radius = model->radius;
    27         isPSF = false;
    28         break;
    29            
    30       default:
    31         return NULL;
    32     }
    33 
    34     x = model->params->data.F32[2];
    35     y = model->params->data.F32[3];
    36 
    37     // replace source flux
    38     pmSourceAddModel (source->pixels, source->mask, model, false, false);
    39 
    40     // set aperture mask circle of PSF_FIT_RADIUS
    41     psImageKeepCircle (source->mask, x, y, radius, "OR", PSPHOT_MASK_MARKED);
    42 
    43     // measure object photometry
    44     status = pmSourcePhotometry (&source->fitMag, &source->apMag, model, source->pixels, source->mask);
    45 
    46     // for PSFs, correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS
    47     if (isPSF && (psf != NULL)) {
    48       if (psf->growth != NULL) {
    49           source->apMag += pmGrowthCurveCorrect (psf->growth, model->radius);
    50       }
    51 
    52       rflux   = pow (10.0, 0.4*source->fitMag);
    53       source->apMag  -= PS_SQR(model->radius)*rflux * psf->skyBias + psf->skySat / rflux;
    54       source->fitMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0);
    55     }
    56 
    57     // unmask aperture
    58     psImageKeepCircle (source->mask, x, y, radius, "AND", ~PSPHOT_MASK_MARKED);
    59 
    60     // subtract object, leave local sky
    61     pmSourceSubModel (source->pixels, source->mask, model, false, false);
    62 
    63     if (!status) return NULL;
    64     return model;
    65 }
    66 
    67 /*
    68 aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y)
    69 (aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y)
    70 (aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y)
    71 
    72 */
    732
    743bool psphotMagnitudes (psMetadata *config, psArray *sources, pmPSF *psf) {
     
    8817}
    8918
     19// XXX should this function use RADIUS or should it
  • trunk/psphot/src/psphotModelGroupInit.c

    r5718 r6571  
    11# include "psphot.h"
     2
     3// Add locally-defined models here.  As these mature, they can be moved to
     4// psModule/src/objects/models
    25
    36# include "models/pmModel_TAUSS.c"
     
    69    {"PS_MODEL_TAUSS", 9, pmModelFunc_TAUSS,  pmModelFlux_TAUSS,  pmModelRadius_TAUSS,  pmModelLimits_TAUSS,  pmModelGuess_TAUSS, pmModelFromPSF_TAUSS, pmModelFitStatus_TAUSS},
    710};
    8 
    911
    1012void psphotModelGroupInit (void)
  • trunk/psphot/src/psphotModelTest.c

    r6481 r6571  
    33static char DEFAULT_MODE[] = "EXT";
    44
    5 bool psphotModelTest (pmReadout *readout, psMetadata *arguments, psMetadata *recipe) {
     5bool psphotModelTest (pmReadout *readout, recipe) {
    66
    77    bool status;
     
    1515
    1616    // run model fitting tests on a single source
    17     if (!psMetadataLookupBool (&status, arguments, "TEST_FIT")) return false;
     17    if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false;
    1818
    1919    // what fitting mode to use?
    20     char *psfModeWord = psMetadataLookupStr (&status, arguments, "TEST_FIT_MODE");
     20    char *psfModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE");
    2121    if (!status) {
    2222        psfModeWord = DEFAULT_MODE;
     
    3333    } else {
    3434        // find the model: supplied by user or first in the PSF_MODEL list
    35         char *modelName  = psMetadataLookupStr (&status, arguments, "TEST_FIT_MODEL");
     35        char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
    3636        if (modelName == NULL) {
    3737            // get the list pointers for the PSF_MODEL entries
     
    5757
    5858    // find the fitting parameters (try test values first)
    59     float INNER = psMetadataLookupF32 (&status, arguments, "TEST_FIT_INNER_RADIUS");
     59    float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS");
    6060    if (!status) {
    6161        INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    6262    }
    6363
    64     float OUTER = psMetadataLookupF32 (&status, arguments, "TEST_FIT_OUTER_RADIUS");
     64    float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS");
    6565    if (!status) {
    6666        OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    6767    }
    6868
    69     float RADIUS = psMetadataLookupF32 (&status, arguments, "TEST_FIT_RADIUS");
     69    float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS");
    7070    if (!status) {
    7171        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
    7272    }
    7373
    74     float mRADIUS = psMetadataLookupF32 (&status, arguments, "TEST_MOMENTS_RADIUS");
     74    float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS");
    7575    if (!status) {
    7676        mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     
    7878
    7979    // define the source of interest
    80     float xObj     = psMetadataLookupF32 (&status, arguments, "TEST_FIT_X");
    81     float yObj     = psMetadataLookupF32 (&status, arguments, "TEST_FIT_Y");
     80    float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
     81    float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
    8282
    8383    // construct the source structures
     
    120120        }
    121121        sprintf (name, "TEST_FIT_PAR%d", i);
    122         value = psMetadataLookupF32 (&status, arguments, name);
     122        value = psMetadataLookupF32 (&status, recipe, name);
    123123        if (status) {
    124124            params[i] = value;
     
    154154    psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PSPHOT_MASK_MARKED);
    155155
    156     char *fitset = psMetadataLookupStr (&status, arguments, "TEST_FIT_SET");
     156    char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET");
    157157    if (status) {
    158158        status = psphotFitSet (source, model, fitset, psfMode);
  • 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
  • trunk/psphot/src/psphotParseCamera.c

    r6379 r6571  
    88      } }
    99
    10 // 2006.02.07 : no leaks!
    11 ppFile *psphotParseCamera (ppConfig *config) {
     10bool psphotParseCamera (pmConfig *config) {
    1211
    13     ppFile *input = ppFileAlloc ();
     12    // psphot is supplied with a list of input images (may be only one image)
     13    psArray *infiles = psMetadataLookupPtr(NULL, config->arguments, "INPUT");
     14    if (infiles->n < 1) psAbort (__func__, "empty input list");
    1415
    15     input->filename = psMetadataLookupStr(NULL, config->arguments, "INPUT_FILE");
    16     psMemCopy (input->filename); // keep for external use
    17 
    18     // Open the input image
    19     psLogMsg("psphot", PS_LOG_INFO, "Opening input image: %s\n", input->filename);
    20     input->fits = psFitsOpen (input->filename, "r"); // File handle for FITS file
    21     if (! input->fits) {
    22         psErrorStackPrint(stderr, "Can't open input image: %s\n", input->filename);
    23         exit(EXIT_FAILURE);
     16    // if no camera has been specified, use the first image as a template for the rest.
     17    if (config->camera == NULL) {
     18        psFits *fits = psFitsOpen (infiles->data[0], "r");
     19        psMetadata *phu = psFitsReadHeader (NULL, fits);
     20        config->camera = pmConfigCameraFromHeader (config->site, phu);
     21        psFitsClose (fits);
     22        psFree (phu);
    2423    }
    25     input->phu = psFitsReadHeader(NULL, input->fits); // FITS header
    26 
    27     // Get camera configuration from header if not already defined
     24    // There's no point in continuing if we can't recognize what we've got
    2825    if (! config->camera) {
    29         config->camera = pmConfigCameraFromHeader(config->site, input->phu);
    30         if (! config->camera) {
    31              // There's no point in continuing if we can't recognize what we've got
    32             psErrorStackPrint(stderr, "Can't find camera configuration!\n");
    33             exit(EXIT_FAILURE);
    34         }
    35    } else if (! pmConfigValidateCamera(config->camera, input->phu)) {
    36        // There's no point in continuing if what we've got doesn't match what we've been told
    37         psError(PS_ERR_IO, true, "%s does not seem to be from the specified camera.\n",
    38                 input->filename);
    39         exit(EXIT_FAILURE);
     26        psErrorStackPrint(stderr, "Can't find camera configuration!\n");
     27        exit(EXIT_FAILURE);
    4028    }
    4129
    42     // Determine the correct recipe to use (from camera or from user)
    43     // if config->recipe is not NULL, it is a user-supplied recipe
    44     if (!config->recipe) {
    45         config->recipe = pmConfigRecipeFromCamera(config->camera, PSPHOT_RECIPE);
    46         if (config->recipe == NULL) {
    47             psErrorStackPrint(stderr, "Can't find recipe configuration!\n");
    48             exit(EXIT_FAILURE);
    49         }
     30    // XXX place the "data" element on the pmConfig structure?
     31    config->files = psMetadataAlloc ();
     32
     33    // build the template fpa, set up the basic view
     34    pmFPA *input = pmFPAConstruct (camera);
     35
     36    // assign the I/O files (potentially) needed by psphot
     37    pmFPAfile *file = pmFPAfileDefine (config->files, camera, input, "PSPHOT.INPUT");
     38
     39    pmFPAview *view = pmFPAviewAlloc (input, config->camera, 0);
     40    for (int i = 0; i < infiles->n; i++) {
     41        // XXX save the phu from above?
     42        psFits *fits = psFitsOpen (infiles->data[i], "r");
     43        psMetadata *phu = psFitsReadHeader (NULL, fits);
     44        pmConfigValidateCamera (config->camera, phu);
     45
     46        // set the view to the corresponding entry for this phu
     47        // XXX should I save the phu data at the appropriate location in the fpa?
     48        pmFPAsetView (view, phu);
     49
     50        // XXX is this the correct psMD to save the filename?
     51        name = pmFPAviewNameFromRule (file->filextra, view);
     52        psMetadataAddStr (file->names, PS_LIST_TAIL, name, 0, "", infiles->data[i]);
    5053    }
     54    pmFPAfileDefine (config->files, camera, input, "PSPHOT.OUTPUT");
     55    pmFPAfileDefine (config->files, camera, input, "PSPHOT.RESID");
     56    pmFPAfileDefine (config->files, camera, input, "PSPHOT.PSF_INPUT");
     57    pmFPAfileDefine (config->files, camera, input, "PSPHOT.PSF_OUTPUT");
     58
     59    // build the template fpa, set up the basic view
     60    // supply the backgnd with a different camera?
     61    pmFPAfileDefine (config->files, camera, NULL, "PSPHOT.BACKSUB");
     62    pmFPAfileDefine (config->files, camera, NULL, "PSPHOT.BACKGND");
     63    pmFPAfileDefine (config->files, camera, NULL, "PSPHOT.BACKMDL");
     64    pmFPAfileDefine (config->files, camera, NULL, "PSPHOT.PSF_SAMPLE");
    5165
    5266    // recipe override values (command-line options):
    53     psMetadataItem *item = NULL;
    54     psMetadataIterator *iter = psMetadataIteratorAlloc (config->options, PS_LIST_HEAD, NULL);
    55     while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
    56         psMetadataAddItem (config->recipe, item, PS_LIST_TAIL, PS_META_REPLACE);
     67    psMetadata *options = psMetadataLookupPtr (status, config->arguments, "PSPHOT.OPTIONS");
     68    psMetadata *recipe = psMetadataLookupPtr (status, config->recipes, "PSPHOT");
     69    psMetadataIterator *iter = psMetadataIteratorAlloc (options, PS_LIST_HEAD, NULL);
     70    while ((psMetadataItem *item = psMetadataGetAndIncrement (iter)) != NULL) {
     71        psMetadataAddItem (recipe, item, PS_LIST_TAIL, PS_META_REPLACE);
    5772    }
    5873    psFree (iter);
    5974
    6075    // set default recipe values here
    61     METADATA_ADD_DEFAULT (config->recipe, Str, "FITMODE", "NONE", "");
    62     METADATA_ADD_DEFAULT (config->recipe, Str, "PHOTCODE", "NONE", "");
    63     METADATA_ADD_DEFAULT (config->recipe, Str, "BREAK_POINT", "NONE", "");
    64     METADATA_ADD_DEFAULT (config->recipe, Str, "OUTPUT_FORMAT", "CMP", "");
    65     METADATA_ADD_DEFAULT (config->recipe, Str, "OUTPUT_MODE", "SPLIT", "");
     76    // XXX this needs re-thinking...
     77    METADATA_ADD_DEFAULT (recipe, Str, "FITMODE",     "NONE", "");
     78    METADATA_ADD_DEFAULT (recipe, Str, "PHOTCODE",    "NONE", "");
     79    METADATA_ADD_DEFAULT (recipe, Str, "BREAK_POINT", "NONE", "");
    6680
    6781    // Chip selection: if we are using a single chip, select it for each FPA
     
    88102        }
    89103    }
    90 
    91     // Construct cameras in preparation for reading
    92     input->fpa = pmFPAConstruct(config->camera);
    93 
    94104    psTrace(__func__, 1, "Done with psphotParseCamera...\n");
    95     return input;
     105    return true;
    96106}
  • trunk/psphot/src/psphotReadout.c

    r6522 r6571  
    22
    33// XXX 2006.02.07 : no leaks!
    4 bool psphotReadout (pmReadout *readout, psMetadata *config) {
     4// XXX change 'config' to 'recipe' eventually
     5bool psphotReadout (pmConfig *config, pmFPAview *view) {
    56
    67    psArray     *sources  = NULL;
    78    psArray     *peaks    = NULL;
    8     psImage     *skymodel = NULL;
    99    pmPSF       *psf      = NULL;
    1010    bool         status;
    1111
     12    // find the currently selected readout
     13    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
     14    pmFPA *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
     15    pmReadout  *readout = pmFPAviewThisReadout (view, input);
     16
    1217    // generate a background model (median, smoothed image)
    13     skymodel = psphotImageMedian (readout, config);
    14 
    15     // optional output of background-subtracted image
    16     psphot
     18    psphotImageMedian (config, view);
    1719
    1820    // find the peaks in the image
    19     peaks = psphotFindPeaks (readout, config);
     21    peaks = psphotFindPeaks (readout, recipe);
    2022
    2123    // construct sources and measure basic stats
    2224    // limit moments analysis by S/N?
    23     sources = psphotSourceStats (readout, config, peaks);
     25    sources = psphotSourceStats (readout, recipe, peaks);
    2426
    2527    // classify sources based on moments, brightness
    2628    // faint sources not classified?
    27     psphotRoughClass (sources, config);
     29    psphotRoughClass (sources, recipe);
    2830
    2931    // mark blended peaks PS_SOURCE_BLEND
    30     psphotBasicDeblend (sources, config);
     32    psphotBasicDeblend (sources, recipe);
    3133
    3234    // use bright stellar objects to measure PSF
    33     psf = psphotChoosePSF (config, sources);
     35    psf = psphotChoosePSF (sources, recipe);
    3436
    3537    // linear PSF fit to peaks
    36     psphotEnsemblePSF (readout, config, sources, psf, FALSE);
     38    psphotEnsemblePSF (readout, sources, recipe, psf, FALSE);
    3739
    3840    // non-linear PSF and EXT fit to brighter sources
    39     psphotBlendFit (readout, config, sources, psf);
     41    psphotBlendFit (readout, sources, recipe, psf);
    4042
    4143    // replace fitted sources
     
    4648
    4749    // linear PSF fit to remaining peaks
    48     psphotEnsemblePSF (readout, config, sources, psf, TRUE);
     50    psphotEnsemblePSF (readout, sources, recipe, psf, TRUE);
    4951
    5052    // measure aperture photometry corrections
    51     psphotApResid (readout, sources, config, psf);
     53    psphotApResid (readout, sources, recipe, psf);
    5254
    5355    // calculate source magnitudes
    54     psphotMagnitudes (config, sources, psf);
     56    psphotMagnitudes (sources, recipe, psf);
    5557
    5658    // update the header with stats results
     59    // XXX need to do this conditionally?
    5760    psMetadata *header = pmReadoutGetHeader (readout);
    5861    psphotUpdateHeader (header, config);
     
    6063    // XXX make this an option?
    6164    // replace background in residual image
    62     psphotSkyReplace (readout, skymodel);
     65    psphotSkyReplace (config, view);
    6366
    6467    // need to do something with the sources, psf, and sky
     
    6770    // XXX : replace? status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SKY.MEAN",  PS_DATA_F32,     "psphot sky mean", sky->sampleMean);
    6871    // XXX : replace? status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SKY.SIGMA", PS_DATA_F32,     "psphot sky stdev", sky->sampleStdev);
    69     // XXX : what should I do with the skymodel image?
    7072
    7173    // free up the local copies of the data
     
    7375    psFree (sources);
    7476    psFree (peaks);
    75     psFree (skymodel);
    7677    return true;
    7778}
  • trunk/psphot/src/psphotSkyReplace.c

    r6481 r6571  
    11# include "psphot.h"
    22
    3 // generate the median in NxN boxes, clipping heavily
    4 // linear interpolation to generate full-scale model
    5 bool psphotSkyReplace (pmReadout *readout, psImage *background) {
     3// in order to  successfully replace the sky, we must define a corresponding file...
     4bool psphotSkyReplace (psMetadata *config, pmFPAview *view) {
    65
     6    // find the currently selected readout
     7    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
     8    pmFPA *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
     9    pmReadout  *readout = pmFPAviewThisReadout (view, input);
     10
     11    // select the corresponding images
    712    psImage *image = readout->image;
    813    psImage *mask  = readout->mask;
     14
     15    // select background pixels, from output background file, or create
     16    background = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKGND", 0, 0, PS_TYPE_F32);
     17    if (background == NULL) {
     18        return false;
     19    }
    920
    1021    // replace the background model
     
    1728    }
    1829   
     30    psFree (background);
    1931    return true;
    2032}
  • trunk/psphot/src/testPSphotLoop.c

    r6522 r6571  
    11# include "psphot.h"
    22
    3 bool psphotLoop (psphotData *data, ppConfig *config) {
     3bool psphotImageLoop (pmConfig *config) {
    44
    5     psRegion region = {0,0,0,0};
    6     pmFPA *fpa = data->input->fpa;
    7     pmFPAiterator *fpi = pmFPAiteratorAlloc (fpa, region);
     5    psMetadata *recipe = psMetadataLookupPts (NULL, config->recipes, PSPHOT_RECIPE);
     6    if (!recipe) {
     7        psErrorStackPrint(stderr, "Can't find recipe configuration!\n");
     8        exit(EXIT_FAILURE);
     9    }
     10
     11    pmFPA *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
     12    pmFPAview *view = pmFPAviewAlloc (0);
    813
    914    // files associated with the science image
    10     pmFileDefine (fpi, config->camera, "PPIMAGE.INPUT");
    11     pmFileDefine (fpi, config->camera, "PPIMAGE.OUTPUT.HEAD");
    12     pmFileDefine (fpi, config->camera, "PPIMAGE.OUTPUT.DATA");
    13     pmFileDefine (fpi, config->camera, "PPIMAGE.BACKGND");
    14     pmFileDefine (fpi, config->camera, "PPIMAGE.BACKSUB");
    15     pmFileDefine (fpi, config->camera, "PPIMAGE.RESID");
    16     pmFileDefine (fpi, config->camera, "PPIMAGE.PSF");
    17     pmFileIOChecks (fpi);
     15    pmFPAfileReadChecks (config->files, view);
    1816
    19     while ((chip = pmChipNext (fpi)) != NULL) {
     17    while ((chip = pmChipNext (view, input)) != NULL) {
    2018
    2119        psLogMsg ("psphot", 4, "Chip %d: %x %x\n", i, chip->exists, chip->process);
    2220        if (! chip->process) { continue; }
    23         pmFileIOChecks (fpi);
     21        pmFPAfileReadChecks (config->files, view);
    2422
    25         while ((cell = pmCellNext (fpi)) != NULL) {
     23        while ((cell = pmCellNext (view, input)) != NULL) {
    2624
    2725            psLogMsg ("psphot", 4, "Cell %d: %x %x\n", j, cell->exists, cell->process);
    2826            if (! cell->process) { continue; }
    29             pmFileIOChecks (fpi);
     27            pmFPAfileReadChecks (config->files, view);
     28
     29            // XXX optional mask and weight input image should be loaded here?
     30            // this sets the weight map and basic mask applying CELL.BAD and CELL.SATURATION
     31            pmCellSetWeights(cell);
     32
     33            // I have a valid mask, now mask in the analysis region of interest
     34            pmCellSetMask (cell, recipe);
    3035
    3136            // process each of the readouts
    32             while ((readout = pmReadoutNext (fpi)) != NULL) {
    33                 pmFileIOChecks (fpi);
     37            while ((readout = pmReadoutNext (view, input)) != NULL) {
     38                pmFPAfileReadChecks (config->files, view);
    3439               
     40                // run a single-model test if desired
     41                // XXX move this to psphotReadout??
     42                // psphotModelTest (readout, recipe);
     43
    3544                // run the actual photometry analysis
    36                 psphotReadout (readout, config);
     45                psphotReadout (config, view);
     46
     47                pmFPAfileWriteChecks (config->files, view);
    3748
    3849                // write out the desired output dataset(s)
    39                 psphotOutput (readout, config);
     50                // psphotOutput (view, recipe);
    4051            }
     52            pmFPAfileWriteChecks (config->files, view);
    4153        }
     54        pmFPAfileWriteChecks (config->files, view);
    4255    }
    43     pmFileClose (fpi);
     56    pmFPAfileWriteChecks (config->files, view);
    4457    return true;
    4558}
     59
     60// I/O files related to psphot:
     61// PSPHOT.INPUT   : input image file(s)
     62// PSPHOT.RESID   : residual image
     63// PSPHOT.OUTPUT  : output object tables (object)
     64
     65// PSPHOT.BACKSUB : background subtracted image
     66// PSPHOT.BACKGND : background model (full-scale image?)
     67// PSPHOT.BACKMDL : background model (binned image?)
     68// PSPHOT.PSF     : sample PSF images
     69
     70
     71/**
     72
     73filename | pmFPAfile  | pmFPA  | pmFPAview
     74chip00.f | PSPHOT.IN  | input  | view
     75chip01.f | PSPHOT.IN  | input  | view
     76chip02.f | PSPHOT.IN  | input  | view
     77chip03.f | PSPHOT.IN  | input  | view
     78
     79out00.f  | PSPHOT.OUT | input  | view
     80out01.f  | PSPHOT.OUT | input  | view
     81out02.f  | PSPHOT.OUT | input  | view
     82out03.f  | PSPHOT.OUT | input  | view
     83
     84obj00.f  | PSPHOT.OBJ | input  | view
     85obj01.f  | PSPHOT.OBJ | input  | view
     86obj02.f  | PSPHOT.OBJ | input  | view
     87obj03.f  | PSPHOT.OBJ | input  | view
     88
     89bin00.f  | PSPHOT.BIN | binned | view
     90bin01.f  | PSPHOT.BIN | binned | view
     91bin02.f  | PSPHOT.BIN | binned | view
     92bin03.f  | PSPHOT.BIN | binned | view
     93
     94mosaic.f | PSPHOT.MOS | mosaic | view
     95
     96**/
Note: See TracChangeset for help on using the changeset viewer.