IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 17396 for trunk/psModules


Ignore:
Timestamp:
Apr 8, 2008, 8:36:06 AM (18 years ago)
Author:
eugene
Message:

merging from eam_branch_20080324 : psphot work on extended source fitting and related I/O functions

Location:
trunk/psModules/src/objects
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/models/pmModel_QGAUSS.c

    r15834 r17396  
    237237    if (!isfinite(shape.sxy)) return false;
    238238
    239     PAR[PM_PAR_SKY]  = moments->Sky;
    240     PAR[PM_PAR_I0]   = moments->Peak - moments->Sky;
     239    // XXX turn this off here for now PAR[PM_PAR_SKY]  = moments->Sky;
     240    PAR[PM_PAR_SKY]  = 0.0;
     241    PAR[PM_PAR_I0]   = moments->Peak;
    241242    PAR[PM_PAR_XPOS] = peak->x;
    242243    PAR[PM_PAR_YPOS] = peak->y;
    … …  
    457458
    458459    status = true;
    459     status &= (dP < 0.5);
     460//    status &= (dP < 0.5);
    460461    status &= (PAR[PM_PAR_I0] > 0);
    461462    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
     463
     464    fprintf (stderr, "QGAUSS status pars: dP: %f, I0: %f, S/N: %f\n",
     465             dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]));
    462466
    463467    return status;

  • trunk/psModules/src/objects/models/pmModel_SERSIC.c

    r15834 r17396  
    1919
    2020   note that a standard sersic model uses exp(-K*(z^(1/n) - 1).  the additional elements (K,
    21    the -1 offset) are absorbed in this model by the normalization, the exponenet, and the
     21   the -1 offset) are absorbed in this model by the normalization, the exponent, and the
    2222   radial scale.  We fit the elements in this form, then re-normalize them on output.
    2323   *****************************************************************************/
    … …  
    157157            break;
    158158        case PM_PAR_SXX:
    159             params_min =   0.25;
     159            params_min =   0.05;
    160160            break;
    161161        case PM_PAR_SYY:
    162             params_min =   0.25;
     162            params_min =   0.05;
    163163            break;
    164164        case PM_PAR_SXY:
    … …  
    166166            break;
    167167        case PM_PAR_7:
    168             params_min =   0.1;
     168            params_min =   0.05;
    169169            break;
    170170        default:
    … …  
    247247    if (!isfinite(shape.sxy)) return false;
    248248
    249     PAR[PM_PAR_SKY]  = moments->Sky;
    250     PAR[PM_PAR_I0]   = moments->Peak - moments->Sky;
     249    // XXX PAR[PM_PAR_SKY]  = moments->Sky;
     250    PAR[PM_PAR_SKY]  = 0.0;
     251    PAR[PM_PAR_I0]   = moments->Peak;
    251252    PAR[PM_PAR_XPOS] = peak->x;
    252253    PAR[PM_PAR_YPOS] = peak->y;
    … …  
    442443
    443444    status = true;
    444     status &= (dP < 0.5);
     445//    status &= (dP < 0.5);
    445446    status &= (PAR[PM_PAR_I0] > 0);
    446447    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
     448
     449    fprintf (stderr, "SERSIC status pars: dP: %f, I0: %f, S/N: %f\n",
     450             dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]));
    447451
    448452    return status;

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

    r15979 r17396  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2008-01-02 20:38:28 $
     8 *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2008-04-08 18:35:38 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    141141    tmp->xf = x;
    142142    tmp->yf = y;
    143 
     143    tmp->assigned = false;
    144144    tmp->type = type;
    145145

  • trunk/psModules/src/objects/pmPeaks.h

    r15984 r17396  
    1010 * @author GLG, MHPCC
    1111 *
    12  * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
    13  * @date $Date: 2008-01-02 20:42:46 $
     12 * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
     13 * @date $Date: 2008-04-08 18:35:38 $
    1414 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1515 */
    … …  
    5858    float flux;                         ///< level in unsmoothed sci image
    5959    float SN;                           ///< S/N implied by detection level
    60     pmPeakType type;   ///< Description of peak.
     60    bool assigned;                      ///< is peak assigned to a source?
     61    pmPeakType type;                    ///< Description of peak.
    6162}
    6263pmPeak;

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

    r17049 r17396  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.51 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2008-03-19 00:51:09 $
     8 *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2008-04-08 18:35:38 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    5050    psFree(tmp->modelPSF);
    5151    psFree(tmp->modelEXT);
    52     psFree(tmp->modelConv);
     52    psFree(tmp->modelFits);
     53    psFree(tmp->extpars);
    5354    psFree(tmp->blends);
    5455    psTrace("psModules.objects", 5, "---- end ----\n");
    … …  
    107108    source->modelPSF = NULL;
    108109    source->modelEXT = NULL;
    109     source->modelConv = NULL;
     110    source->modelFits = NULL;
    110111    source->type = PM_SOURCE_TYPE_UNKNOWN;
    111112    source->mode = PM_SOURCE_MODE_DEFAULT;
    … …  
    947948        return model;
    948949
    949 // XXX when should I return the modelConv ??
     950        // the 'best' extended model is saved in source->modelEXT (may be a pointer to one of
     951        // the elements of source->modelFits)
    950952      case PM_SOURCE_TYPE_EXTENDED:
    951         model = source->modelConv;
    952         if (!model) {
    953             model = source->modelEXT;
    954         }
     953        model = source->modelEXT;
    955954        if (!model && source->modelPSF) {
     955            // XXX raise an error or warning here?
    956956            if (isPSF) {
    957957                *isPSF = true;
    … …  
    10151015  if (!strcasecmp (name, "DEFECT"    )) return PM_SOURCE_MODE_DEFECT;
    10161016  if (!strcasecmp (name, "SATURATED" )) return PM_SOURCE_MODE_SATURATED;
    1017   if (!strcasecmp (name, "CRLIMIT"   )) return PM_SOURCE_MODE_CRLIMIT;
     1017  if (!strcasecmp (name, "CRLIMIT"   )) return PM_SOURCE_MODE_CR_LIMIT;
     1018  if (!strcasecmp (name, "EXTLIMIT"  )) return PM_SOURCE_MODE_EXT_LIMIT;
    10181019  if (!strcasecmp (name, "SUBTRACTED")) return PM_SOURCE_MODE_SUBTRACTED;
    10191020  return PM_SOURCE_MODE_DEFAULT;
    … …  
    10361037    case PM_SOURCE_MODE_DEFECT     : return psStringCopy ("DEFECT"    );
    10371038    case PM_SOURCE_MODE_SATURATED  : return psStringCopy ("SATURATED" );
    1038     case PM_SOURCE_MODE_CRLIMIT    : return psStringCopy ("CRLIMIT");
     1039    case PM_SOURCE_MODE_CR_LIMIT   : return psStringCopy ("CRLIMIT"   );
     1040    case PM_SOURCE_MODE_EXT_LIMIT  : return psStringCopy ("EXTLIMIT"  );
    10391041    case PM_SOURCE_MODE_SUBTRACTED : return psStringCopy ("SUBTRACTED");
    10401042    default:

  • trunk/psModules/src/objects/pmSource.h

    r16819 r17396  
    33 * @author EAM, IfA; GLG, MHPCC
    44 *
    5  * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
    6  * @date $Date: 2008-03-05 01:08:08 $
     5 * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
     6 * @date $Date: 2008-04-08 18:35:38 $
    77 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    88 */
    … …  
    2222 * source.
    2323 *
    24  * XXX: The values given below are currently illustrative and will require
    25  * some modification as the source classification code is developed. (TBD)
    26  *
    2724 */
    2825typedef enum {
    29     PM_SOURCE_TYPE_UNKNOWN,  ///< a cosmic-ray
    30     PM_SOURCE_TYPE_DEFECT,  ///< a cosmic-ray
    31     PM_SOURCE_TYPE_SATURATED,  ///< random saturated pixels
    32     PM_SOURCE_TYPE_STAR,  ///< a good-quality star
    33     PM_SOURCE_TYPE_EXTENDED,  ///< an extended object (eg, galaxy)
     26    PM_SOURCE_TYPE_UNKNOWN,             ///< not yet classified
     27    PM_SOURCE_TYPE_DEFECT,              ///< a cosmic-ray
     28    PM_SOURCE_TYPE_SATURATED,           ///< random saturated pixels (eg, bleed trails)
     29    PM_SOURCE_TYPE_STAR,                ///< a good-quality star (subtracted model is PSF)
     30    PM_SOURCE_TYPE_EXTENDED,            ///< an extended object (eg, galaxy) (subtracted model is EXT)
    3431} pmSourceType;
    3532
    … …  
    4845    PM_SOURCE_MODE_BADPSF     = 0x0400, ///< Failed to get good estimate of object's PSF
    4946    PM_SOURCE_MODE_DEFECT     = 0x0800, ///< Source is thought to be a defect
    50     PM_SOURCE_MODE_SATURATED  = 0x1000, ///< Source is thought to be saturation
    51     PM_SOURCE_MODE_CRLIMIT    = 0x2000, ///< Source has crNsigma above limit
    52     PM_SOURCE_MODE_SUBTRACTED = 0x4000, ///< XXX this flag is actually only used internally (move)
     47    PM_SOURCE_MODE_SATURATED  = 0x1000, ///< Source is thought to be saturated pixels (bleed trail)
     48    PM_SOURCE_MODE_CR_LIMIT   = 0x2000, ///< Source has crNsigma above limit
     49    PM_SOURCE_MODE_EXT_LIMIT  = 0x4000, ///< Source has extNsigma above limit
     50    PM_SOURCE_MODE_SUBTRACTED = 0x8000, ///< XXX this flag is actually only used internally (move)
    5351} pmSourceMode;
    5452
    … …  
    6058 *
    6159 *  XXX do I have to re-organize this (again!) to allow an arbitrary set of extended model fits??
     60 *  XXX put the Mag and Err inside the pmModel?
     61 *  XXX keep the modelEXT or add to the psArray
     62 * 
    6263 *
    6364 */
    … …  
    6566    const int id;                       ///< Unique ID for object
    6667    int seq;                            ///< ID for output (generated on write)
    67     pmPeak *peak;                       ///< Description of peak pixel.
     68    pmPeak  *peak;                      ///< Description of peak pixel.
    6869    psImage *pixels;                    ///< Rectangular region including object pixels.
    6970    psImage *weight;                    ///< Image variance.
    … …  
    7273    psImage *modelFlux;                 ///< cached copy of the best model for this source
    7374    psImage *psfFlux;                   ///< cached copy of the psf model for this source
    74     pmMoments *moments;                 ///< Basic moments measure for the object.
     75    pmMoments *moments;                 ///< Basic moments measured for the object.
    7576    pmModel *modelPSF;                  ///< PSF Model fit (parameters and type)
    76     pmModel *modelEXT;                  ///< EXT (floating) Model fit (parameters and type).
    77     pmModel *modelConv;                 ///< PSF-Convolved Model fit (parameters and type).
     77    pmModel *modelEXT;                  ///< EXT Model fit used for subtraction (parameters and type)
     78    psArray *modelFits;                 ///< collection of extended source models (best == modelEXT)
    7879    pmSourceType type;                  ///< Best identification of object.
    79     pmSourceMode mode;                  ///< Best identification of object.
    80     psArray *blends;
    81     float psfMag;                       ///< calculated from flux in modelPsf
     80    pmSourceMode mode;                  ///< analysis flags set for object.
     81    psArray *blends;                    ///< collection of sources thought to be confused with object
     82    float psfMag;                       ///< calculated from flux in modelPSF
    8283    float extMag;                       ///< calculated from flux in modelEXT
    8384    float errMag;                       ///< error in psfMag OR extMag (depending on type)
    8485    float apMag;                        ///< apMag corresponding to psfMag or extMag (depending on type)
    8586    float pixWeight;                    ///< model-weighted coverage of valid pixels
    86     float psfChisq;                      ///< probability of PSF
     87    float psfChisq;                     ///< probability of PSF
    8788    float crNsigma;                     ///< Nsigma deviation from PSF to CR
    8889    float extNsigma;                    ///< Nsigma deviation from PSF to EXT
     90    float sky, skyErr;                  ///< The sky and its error at the center of the object
    8991    psRegion region;                    ///< area on image covered by selected pixels
    90     float sky, skyErr;                  ///< The sky and its error at the center of the object
    9192    pmSourceExtendedPars *extpars;      ///< extended source parameters
    9293};

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

    r17010 r17396  
    33 *  @author EAM, IfA
    44 *
    5  *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2008-03-17 22:04:27 $
     5 *  @version $Revision: 1.56 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2008-04-08 18:35:38 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    327327
    328328        // if this is not TRUE, the output files only contain the psf measurements.
    329         bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "SAVE.XSRC");
    330         bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "SAVE.XFIT");
     329        bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS");
     330        bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS");
    331331
    332332        // define the EXTNAME values for the different data segments:
    … …  
    445445            if (xsrcname) {
    446446              if (!strcmp (exttype, "PS1_DEV_1")) {
    447                 status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname);
     447                status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe);
    448448              }
    449449            }
    … …  
    457457                psFree (headname);
    458458                psFree (dataname);
     459                psFree (xsrcname);
     460                psFree (xfitname);
    459461                psFree (outhead);
    460462                return false;
    … …  
    466468        psFree (headname);
    467469        psFree (dataname);
     470        psFree (xsrcname);
     471        psFree (xfitname);
    468472        psFree (outhead);
    469473        break;

  • trunk/psModules/src/objects/pmSourceIO.h

    r16819 r17396  
    44 * @author EAM, IfA; GLG, MHPCC
    55 *
    6  * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2008-03-05 01:08:08 $
     6 * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
     7 * @date $Date: 2008-04-08 18:35:38 $
    88 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    99 *
    … …  
    2727bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
    2828bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, char *xsrcname);
    29 bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname);
     29bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe);
    3030bool pmSourcesWrite_PS1_DEV_1_XFIT (psFits *fits, psArray *sources, char *extname);
    3131

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

    r16819 r17396  
    33 *  @author EAM, IfA
    44 *
    5  *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2008-03-05 01:08:08 $
     5 *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2008-04-08 18:35:38 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    251251}
    252252
    253 bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname)
     253bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
    254254{
    255255
     256    bool status;
    256257    psArray *table;
    257258    psMetadata *row;
    258     int i;
    259259    psF32 *PAR, *dPAR;
    260     psEllipseAxes axes;
    261260    psF32 xPos, yPos;
    262261    psF32 xErr, yErr;
    … …  
    273272    table = psArrayAllocEmpty (sources->n);
    274273
     274    // which extended source analyses should we perform?
     275    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     276    bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
     277    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
     278    bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
     279
     280    psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
     281    psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
     282    assert (radialBinsLower->n == radialBinsUpper->n);
     283
    275284    // we write out all sources, regardless of quality.  the source flags tell us the state
    276     for (i = 0; i < sources->n; i++) {
     285    for (int i = 0; i < sources->n; i++) {
    277286        // skip source if it is not a ext sourc
    278287        // XXX we have two places that extended source parameters are measured:
    … …  
    283292        pmSource *source = sources->data[i];
    284293
    285         // choose the convolved EXT model, if available, otherwise the simple one
    286         pmModel *model = source->modelConv;
    287         if (model == NULL) {
    288             model = source->modelEXT;
    289         }
    290         if (model == NULL) continue;
    291 
    292         // XXX do we need to do this?
     294        // skip sources without measurements
    293295        if (source->extpars == NULL) continue;
     296
     297        // we require a PSF model fit (ignore the real crud)
     298        pmModel *model = source->modelPSF;
     299        if (model == NULL) continue;
    294300
    295301        // XXX I need to split the extended models from the extended aperture measurements
    … …  
    301307        yErr = dPAR[PM_PAR_YPOS];
    302308
    303         // XXX for the aperture values, I probably can / should just refer to the psf position and not write any of this
    304         axes = pmPSF_ModelToAxes (PAR, 20.0);
    305 
    306309        row = psMetadataAlloc ();
     310
    307311        // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
    308312        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    … …  
    311315        psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG",        PS_DATA_F32, "Sigma in EXT x coordinate",                  xErr);
    312316        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    313         psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG",     PS_DATA_F32, "EXT fit instrumental magnitude",             source->extMag);
    314         psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    315 
    316         // XXX these should be major and minor, not 'x' and 'y'
    317         psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    PS_DATA_F32, "EXT width in x coordinate",                  axes.major);
    318         psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    PS_DATA_F32, "EXT width in y coordinate",                  axes.minor);
    319         psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA",        PS_DATA_F32, "EXT orientation angle",                      axes.theta);
    320 
    321         // XXX at the moment, this will be a fixed sized table, but perhaps have multiple output formats depending on what is measured?
    322 
    323         // other values that I need to report:
    324         // R, Mag Petrosian, .....
    325         if (source->extpars) {
    326 
    327             // Petrosian measurements
     317
     318        // Petrosian measurements
     319        // XXX insert header data: petrosian ref radius, flux ratio
     320        if (doPetrosian) {
    328321            pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
    329322            if (petrosian) {
    … …  
    332325                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
    333326                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
     327            } else {
     328                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
     329                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", NAN);
     330                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
     331                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
    334332            }
    335 
    336             // Kron measurements
     333        }
     334
     335        // Kron measurements
     336        if (doKron) {
    337337            pmSourceKronValues *kron = source->extpars->kron;
    338338            if (kron) {
    … …  
    341341                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
    342342                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
     343            } else {
     344                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       NAN);
     345                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", NAN);
     346                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          NAN);
     347                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    NAN);
    343348            }
    344 
    345             // Isophot measurements
    346             // XXX insert header data: isophotal level
     349        }
     350
     351        // Isophot measurements
     352        // XXX insert header data: isophotal level
     353        if (doIsophotal) {
    347354            pmSourceIsophotalValues *isophot = source->extpars->isophot;
    348355            if (isophot) {
    … …  
    351358                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
    352359                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
     360            } else {
     361                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       NAN);
     362                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", NAN);
     363                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          NAN);
     364                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    NAN);
    353365            }
    354 
     366        }
     367
     368        // Flux Annuli
     369        if (doAnnuli) {
    355370            pmSourceAnnuli *annuli = source->extpars->annuli;
    356371            if (annuli) {
    … …  
    367382                    sprintf (name, "FLUX_VAR_R_%02d", j);
    368383                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
    369                 }
     384                }
     385            } else {
     386                for (int j = 0; j < radialBinsLower->n; j++) {
     387                    char name[32];
     388                    sprintf (name, "FLUX_VAL_R_%02d", j);
     389                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
     390                    sprintf (name, "FLUX_ERR_R_%02d", j);
     391                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
     392                    sprintf (name, "FLUX_VAR_R_%02d", j);
     393                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
     394                }
    370395            }
    371396        }
    372397
    373         psArrayAdd (table, 100, row);
    374         psFree (row);
     398        psArrayAdd (table, 100, row);
     399        psFree (row);
    375400    }
    376401
    377402    if (table->n == 0) {
    378         psFitsWriteBlank (fits, outhead, extname);
    379         psFree (table);
    380         return true;
     403        psFitsWriteBlank (fits, outhead, extname);
     404        psFree (outhead);
     405        psFree (table);
     406        return true;
    381407    }
    382408
    383409    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    384410    if (!psFitsWriteTable (fits, outhead, table, extname)) {
    385         psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
    386         psFree(table);
    387         return false;
    388     }
     411        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
     412        psFree (outhead);
     413        psFree(table);
     414        return false;
     415    }
     416    psFree (outhead);
    389417    psFree (table);
    390418
    … …  
    397425    psArray *table;
    398426    psMetadata *row;
    399     int i;
    400427    psF32 *PAR, *dPAR;
    401428    psEllipseAxes axes;
    402429    psF32 xPos, yPos;
    403430    psF32 xErr, yErr;
     431    char name[64];
    404432
    405433    // create a header to hold the output data
    … …  
    412440    sources = psArraySort (sources, pmSourceSortBySN);
    413441
     442    // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams
     443    int nParamMax = 0;
     444    for (int i = 0; i < sources->n; i++) {
     445        pmSource *source = sources->data[i];
     446        if (source->modelFits == NULL) continue;
     447        for (int j = 0; j < source->modelFits->n; j++) {
     448            pmModel *model = source->modelFits->data[j];
     449            assert (model);
     450            nParamMax = PS_MAX (nParamMax, model->params->n);
     451        }
     452    }
     453
    414454    table = psArrayAllocEmpty (sources->n);
    415455
    416456    // we write out all sources, regardless of quality.  the source flags tell us the state
    417     for (i = 0; i < sources->n; i++) {
    418         // skip source if it is not a ext sourc
    419         // XXX we have two places that extended source parameters are measured:
    420         // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models,
    421         // psphotFitEXT, which does the simple extended source model fit (not psf-convolved)
    422         // should we require both?
     457    for (int i = 0; i < sources->n; i++) {
    423458
    424459        pmSource *source = sources->data[i];
    425460
    426         // XXX need to have an array of model fits; loop over the array and write out all
    427         // which we calculated
    428 
    429         // choose the convolved EXT model, if available, otherwise the simple one
    430         // XXX should not need to choose: write both out
    431         pmModel *model = source->modelConv;
    432         if (model == NULL) {
    433             model = source->modelEXT;
     461        // XXX if no model fits are saved, write out modelEXT?
     462        if (source->modelFits == NULL) continue;
     463
     464        // We have multiple sources : need to flag the one used to subtract the light (the 'best' model)
     465        for (int j = 0; j < source->modelFits->n; j++) {
     466
     467            // choose the convolved EXT model, if available, otherwise the simple one
     468            pmModel *model = source->modelFits->data[j];
     469            assert (model);
     470
     471            PAR = model->params->data.F32;
     472            dPAR = model->dparams->data.F32;
     473            xPos = PAR[PM_PAR_XPOS];
     474            yPos = PAR[PM_PAR_YPOS];
     475            xErr = dPAR[PM_PAR_XPOS];
     476            yErr = dPAR[PM_PAR_YPOS];
     477
     478            axes = pmPSF_ModelToAxes (PAR, 20.0);
     479
     480            row = psMetadataAlloc ();
     481
     482            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
     483            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
     484            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
     485            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT",            0, "EXT model y coordinate",                     yPos);
     486            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG",        0, "Sigma in EXT x coordinate",                  xErr);
     487            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG",        0, "Sigma in EXT y coordinate",                  yErr);
     488            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
     489            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
     490
     491            psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
     492            psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "name of model",                              pmModelClassGetName (model->type));
     493
     494            // XXX these should be major and minor, not 'x' and 'y'
     495            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width in x coordinate",                  axes.major);
     496            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width in y coordinate",                  axes.minor);
     497            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                      axes.theta);
     498
     499            // write out the other generic parameters
     500            for (int k = 0; k < nParamMax; k++) {
     501                if (k == PM_PAR_I0) continue;
     502                if (k == PM_PAR_SKY) continue;
     503                if (k == PM_PAR_XPOS) continue;
     504                if (k == PM_PAR_YPOS) continue;
     505                if (k == PM_PAR_SXX) continue;
     506                if (k == PM_PAR_SXY) continue;
     507                if (k == PM_PAR_SYY) continue;
     508
     509                snprintf (name, 64, "EXT_PAR_%02d", k);
     510
     511                if (k < model->params->n) {
     512                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
     513                } else {
     514                    psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
     515                }
     516            }
     517
     518            // XXX other parameters which may be set.
     519            // XXX flag / value to define the model
     520            // XXX write out the model type, fit status flags
     521
     522            psArrayAdd (table, 100, row);
     523            psFree (row);
    434524        }
    435         if (model == NULL) continue;
    436 
    437         PAR = model->params->data.F32;
    438         dPAR = model->dparams->data.F32;
    439         xPos = PAR[PM_PAR_XPOS];
    440         yPos = PAR[PM_PAR_YPOS];
    441         xErr = dPAR[PM_PAR_XPOS];
    442         yErr = dPAR[PM_PAR_YPOS];
    443 
    444         axes = pmPSF_ModelToAxes (PAR, 20.0);
    445 
    446         row = psMetadataAlloc ();
    447         // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
    448         psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    449         psMetadataAdd (row, PS_LIST_TAIL, "X_EXT",            PS_DATA_F32, "EXT model x coordinate",                     xPos);
    450         psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT",            PS_DATA_F32, "EXT model y coordinate",                     yPos);
    451         psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG",        PS_DATA_F32, "Sigma in EXT x coordinate",                  xErr);
    452         psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    453         psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG",     PS_DATA_F32, "EXT fit instrumental magnitude",             source->extMag);
    454         psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    455 
    456         // XXX these should be major and minor, not 'x' and 'y'
    457         psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    PS_DATA_F32, "EXT width in x coordinate",                  axes.major);
    458         psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    PS_DATA_F32, "EXT width in y coordinate",                  axes.minor);
    459         psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA",        PS_DATA_F32, "EXT orientation angle",                      axes.theta);
    460 
    461         // XXX other parameters which may be set.
    462         // XXX flag / value to define the model
    463 
    464         psArrayAdd (table, 100, row);
    465         psFree (row);
    466525    }
    467526
    468527    if (table->n == 0) {
    469         psFitsWriteBlank (fits, outhead, extname);
    470         psFree (table);
    471         return true;
     528        psFitsWriteBlank (fits, outhead, extname);
     529        psFree (outhead);
     530        psFree (table);
     531        return true;
    472532    }
    473533
    474534    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    475535    if (!psFitsWriteTable (fits, outhead, table, extname)) {
    476         psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
    477         psFree(table);
    478         return false;
    479     }
     536        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
     537        psFree (outhead);
     538        psFree(table);
     539        return false;
     540    }
     541    psFree (outhead);
    480542    psFree (table);
    481 
    482543    return true;
    483544}

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

    r17287 r17396  
    33 *  @author EAM, IfA; GLG, MHPCC
    44 *
    5  *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2008-04-02 22:40:36 $
     5 *  @version $Revision: 1.41 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2008-04-08 18:35:38 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    119119    }
    120120
    121     // measure EXT model photometry (do both modelEXT and modelConv or just the one?)
    122     status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
     121    // if we have a collection of model fits, one of them is a pointer to modelEXT?
     122    // XXX not guaranteed
     123    if (source->modelFits) {
     124      for (int i = 0; i < source->modelFits->n; i++) {
     125        pmModel *model = source->modelFits->data[i];
     126        status = pmSourcePhotometryModel (&model->mag, model);
     127      }
     128      if (source->modelEXT) {
     129        source->extMag = source->modelEXT->mag;
     130      }
     131    } else {
     132      if (source->modelEXT) {
     133        status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
     134      }
     135    }
    123136
    124137    // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
Note: See TracChangeset for help on using the changeset viewer.