IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 7, 2014, 5:40:37 AM (12 years ago)
Author:
eugene
Message:

merge changes from trunk

Location:
branches/eam_branches/ipp-20140423/psModules/src/objects
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20140423/psModules/src/objects/pmSource.c

    r36623 r36829  
    162162    source->sky              = NAN;
    163163    source->skyErr           = NAN;   
     164    source->extSN            = NAN;   
    164165
    165166    source->region = psRegionSet(NAN, NAN, NAN, NAN);

  • branches/eam_branches/ipp-20140423/psModules/src/objects/pmSource.h

    r36623 r36829  
    113113    float sky;                          ///< The sky at the center of the object
    114114    float skyErr;                       ///< The sky error at the center of the object
     115    float extSN;                        ///< for externally supplied source the kron signal to noise (used by full force)
    115116
    116117    psRegion region;                    ///< area on image covered by selected pixels
    … …  
    119120    pmSourceExtendedPars *extpars;      ///< extended source parameters
    120121    pmSourceDiffStats *diffStats;       ///< extra parameters for difference detections
    121     pmSourceGalaxyFits *galaxyFits;     ///< fits to galaxy models (psphotFullForce only)
     122    psArray *galaxyFits;                ///< fits to galaxy models (psphotFullForce only)
    122123    pmSourceLensing *lensingOBJ;        ///< lensing moments parameters (per object)
    123124    pmSourceLensing *lensingPSF;        ///< lensing moments parameters (psf, interpolated)

  • branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceExtendedPars.h

    r36623 r36829  
    8383
    8484typedef struct {
     85  int       modelType;
    8586  psVector *Flux;
    8687  psVector *dFlux;
    8788  psVector *chisq;
    88   int nPix;
     89  int       nPix;
     90  bool      reducedTrials;
     91  float     fRmajorMin;
     92  float     fRmajorMax;
     93  float     fRmajorDel;
     94  float     fRminorMin;
     95  float     fRminorMax;
     96  float     fRminorDel;
    8997} pmSourceGalaxyFits;
    9098

  • branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO.c

    r36623 r36829  
    12021202        break;
    12031203
    1204       case PM_FPA_FILE_CFF:
     1204      case PM_FPA_FILE_CFF: {
     1205        // determine the output table format
     1206        psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSPHOT");
     1207        if (!status) {
     1208            psError(PS_ERR_UNKNOWN, true, "missing recipe PSPHOT in config data");
     1209            return false;
     1210        }
    12051211        // read in header, if not yet loaded
    12061212        hdu = pmFPAviewThisHDU (view, file->fpa);
    … …  
    12441250        }
    12451251
    1246         sources = pmSourcesRead_CFF(file->fits, hdu->header);
     1252        sources = pmSourcesRead_CFF(file->fits, hdu->header, recipe);
    12471253
    12481254        psTrace("psModules.objects", 6, "read CMF table from %s : %s : %s", file->filename, headname, dataname);
    … …  
    12501256        psFree (dataname);
    12511257        psFree (tableHeader);
     1258        }
    12521259        break;
    12531260

  • branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO.h

    r36623 r36829  
    5555
    5656psArray *pmSourcesReadCMP (char *filename, psMetadata *header);
    57 psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header);
     57psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header, psMetadata *recipe);
    5858bool pmSourcesWrite_CFF (pmReadout *readout, psFits *fits, psArray *sources, psMetadata *header, psMetadata *recipe);
    5959

  • branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO_CFF.c

    r36633 r36829  
    5353#include "pmSourceOutputs.h"
    5454
    55 // read in a readout from the fits file
    56 psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header)
     55// read in sources readout from a cff fits file
     56psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header, psMetadata *recipe)
    5757{
    5858    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    7575    PS_ASSERT_INT_NONNEGATIVE(modelType, NULL);
    7676
     77    pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC");
     78    pmModelType devModelType    = pmModelClassGetType("PS_MODEL_DEV");
     79
     80    psString modelForce = psMetadataLookupStr(&status, recipe, "EXT_MODEL_TYPE_FORCE");
     81    psF32 forceDevSersicMin = NAN;
     82    bool forceAll = false;
     83    int modelTypeForce = 0;
     84    if (!strcmp(modelForce, "ALL")) {
     85        forceAll = true;
     86    } else if (!strcmp(modelForce, "PS_MODEL_SERSIC")) {
     87        modelTypeForce = sersicModelType;
     88        forceDevSersicMin = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_DEV_SERSIC_MIN");
     89    } else {
     90        modelTypeForce = pmModelClassGetType(modelForce);
     91        PS_ASSERT_INT_NONNEGATIVE(modelTypeForce, NULL);
     92    }
     93
     94    // skip extended model types for likely stars
     95    // max value of KronMag - psfMag to keep ...
     96    psF32 starCut = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_MAGDIFF_MAX");
     97    if (!status) {
     98        starCut = 0;
     99    }
     100    // ... unless SN is less than this value
     101    psF32 SNMinForCut = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_CUT_SN_MIN");
     102    if (!status) {
     103        SNMinForCut = 10;
     104    }
     105
    77106    // We get the size of the table, and allocate the array of sources first because the table
    78107    // is large and ephemeral --- when the table gets blown away, whatever is allocated after
    79108    // the table is read blocks the free.  In fact, it's better to read the table row by row.
    80     long numSources = psFitsTableSize(fits); // Number of sources in table
    81     psArray *sources = psArrayAlloc(numSources); // Array of sources, to return
    82 
    83     // convert the table to the pmSource entriesa
    84     for (int i = 0; i < numSources; i++) {
     109    long numRows = psFitsTableSize(fits); // Number of rows in table
     110    psArray *sources = psArrayAllocEmpty(numRows); // Array of sources, to return
     111
     112    // convert the table to the pmSource entries
     113    pmSource *source = NULL;
     114    for (int i = 0; i < numRows; i++) {
    85115        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
    86116        if (!row) {
    … …  
    108138        float kronRadius = psMetadataLookupF32 (&status, row, "KRON_RADIUS");
    109139        float petRadius  = psMetadataLookupF32 (&status, row, "PETRO_RADIUS");
     140        float SN         = psMetadataLookupF32 (&status, row, "SN");
    110141        bool fitGalaxy   = psMetadataLookupU8 (&status, row, "FIT_GALAXY");
    111142        bool psfStar     = psMetadataLookupU8 (&status, row, "PSF_STAR");
    … …  
    114145        float Rminor     = psMetadataLookupF32 (&status, row, "R_MINOR");
    115146        float theta      = psMetadataLookupF32 (&status, row, "THETA");
     147        float chisq      = psMetadataLookupF32 (&status, row, "CHISQ");
     148        float nDOF       = psMetadataLookupF32 (&status, row, "NDOF");
     149        float magDiff    = psMetadataLookupF32 (&status, row, "MAG_DIFF");
     150        psS16 modelFlags = psMetadataLookupS32 (&status, row, "MODEL_FLAGS");
    116151
    117152        int   galaxyModelType = psMetadataLookupS32(&status, row, "MODEL_TYPE");
    118         if (status) {
     153        if (status && galaxyModelType >= 0) {
    119154            galaxyModelType = pmModelClassGetLocalType(galaxyModelType);
    120155        } else {
    … …  
    123158        float Sindex     = psMetadataLookupF32 (&status, row, "INDEX"); // Should this be PAR_07 not sersic index
    124159
    125         pmSource *source = pmSourceAlloc ();
    126         pmModel *model = pmModelAlloc (modelType);
    127         source->modelPSF  = model;
    128 //        RoughClass wants source type to be unknown
    129 //        source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags
    130         source->type = PM_SOURCE_TYPE_UNKNOWN;
    131 
    132         // XXX we can set this in general, but for a specific image, we need to weed out SATSTARS and
    133         // stars that are masked
    134         if (psfStar) {
    135             source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR;
    136         }
    137 
    138         // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models.
    139         psF32 *PAR = model->params->data.F32;
    140         psF32 *dPAR = model->dparams->data.F32;
    141 
    142         source->seq       = ID;
    143         PAR[PM_PAR_XPOS]  = X;
    144         PAR[PM_PAR_YPOS]  = Y;
    145 
    146         dPAR[PM_PAR_XPOS] = 0.0;
    147         dPAR[PM_PAR_YPOS] = 0.0;
    148 
    149         PAR[PM_PAR_SKY]   = 0.0;
    150         dPAR[PM_PAR_SKY]  = 0.0;
    151 
    152         PAR[PM_PAR_I0]    = 1.0;
    153         dPAR[PM_PAR_I0]   = 0.0;
    154 
    155         source->sky       = PAR[PM_PAR_SKY];
    156         source->skyErr    = dPAR[PM_PAR_SKY];
    157 
    158         source->psfMag    = 0.0;
    159         source->psfMagErr = 0.0;
    160         source->apMag     = 0.0;
    161         source->apRadius  = apRadius;
    162 
    163         // we generate a somewhat fake PSF model here --
    164         // in most (all?) contexts, we will replace this with a measured psf model
    165         // elsewhere
    166         axes.major        = 1.0;
    167         axes.minor        = 1.0;
    168         axes.theta        = 0.0;
    169         pmPSF_AxesToModel (PAR, axes, modelType);
    170 
    171         // peak->detValue, rawFlux, smoothFlux are all set to the flux argument which is counts per second
    172         source->peak      = pmPeakAlloc(X, Y, flux, PM_PEAK_LONE);
    173         source->peak->xf  = X; // pmPeakAlloc converts X,Y to int, so reset here
    174         source->peak->yf  = Y; // pmPeakAlloc converts X,Y to int, so reset here
    175         source->peak->dx  = 0.0;
    176         source->peak->dy  = 0.0;
    177 
    178         source->moments = pmMomentsAlloc ();
    179         source->moments->Mx = X;
    180         source->moments->My = Y;
    181         source->moments->Mrf = kronRadius * 0.4; // kronRadius is 2.5 * first radial moment
    182 
    183         // Don't mark the moments as measured because that causes many fields to be left blank.
    184         // The moments code knows not to change the position or the Mrf for external sources
    185         // source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;
    186 
    187         if (isfinite(petRadius)) {
    188             source->extpars = pmSourceExtendedParsAlloc ();
    189             source->extpars->petrosianRadius = petRadius;
    190         }
    191 
     160        if (!source || ID != source->seq) {
     161            if (source) {
     162                psArrayAdd (sources, 1, source);
     163                psFree(source);
     164            }
     165            source = pmSourceAlloc ();
     166            pmModel *model = pmModelAlloc (modelType);
     167            source->modelPSF  = model;
     168            //        RoughClass wants source type to be unknown
     169            //        source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags
     170            source->type = PM_SOURCE_TYPE_UNKNOWN;
     171
     172            // XXX we can set this in general, but for a specific image, we need to weed out SATSTARS and
     173            // stars that are masked
     174            if (psfStar) {
     175                source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR;
     176            }
     177
     178            // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models.
     179            psF32 *PAR = model->params->data.F32;
     180            psF32 *dPAR = model->dparams->data.F32;
     181
     182            source->seq       = ID;
     183
     184            PAR[PM_PAR_XPOS]  = X;
     185            PAR[PM_PAR_YPOS]  = Y;
     186
     187            dPAR[PM_PAR_XPOS] = 0.0;
     188            dPAR[PM_PAR_YPOS] = 0.0;
     189
     190            PAR[PM_PAR_SKY]   = 0.0;
     191            dPAR[PM_PAR_SKY]  = 0.0;
     192
     193            PAR[PM_PAR_I0]    = 1.0;
     194            dPAR[PM_PAR_I0]   = 0.0;
     195
     196            source->sky       = PAR[PM_PAR_SKY];
     197            source->skyErr    = dPAR[PM_PAR_SKY];
     198
     199            source->psfMag    = 0.0;
     200            source->psfMagErr = 0.0;
     201            source->apMag     = 0.0;
     202            source->apRadius  = apRadius;
     203
     204            // we generate a somewhat fake PSF model here --
     205            // in most (all?) contexts, we will replace this with a measured psf model
     206            // elsewhere
     207            axes.major        = 1.0;
     208            axes.minor        = 1.0;
     209            axes.theta        = 0.0;
     210            pmPSF_AxesToModel (PAR, axes, modelType);
     211
     212            // peak->detValue, rawFlux, smoothFlux are all set to the flux argument which is counts per second
     213            source->peak      = pmPeakAlloc(X, Y, flux, PM_PEAK_LONE);
     214            source->peak->xf  = X; // pmPeakAlloc converts X,Y to int, so reset here
     215            source->peak->yf  = Y; // pmPeakAlloc converts X,Y to int, so reset here
     216            source->peak->dx  = 0.0;
     217            source->peak->dy  = 0.0;
     218
     219            source->extSN     = SN;
     220
     221            source->moments = pmMomentsAlloc ();
     222            source->moments->Mx = X;
     223            source->moments->My = Y;
     224            source->moments->Mrf = kronRadius * 0.4; // kronRadius is 2.5 * first radial moment
     225
     226            // Don't mark the moments as measured because that causes many fields to be left blank.
     227            // The moments code knows not to change the position or the Mrf for external sources
     228            // source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;
     229
     230            if (isfinite(petRadius)) {
     231                source->extpars = pmSourceExtendedParsAlloc ();
     232                source->extpars->petrosianRadius = petRadius;
     233            }
     234
     235        }
     236        bool saveExtModelParams = false;
    192237        if (fitGalaxy && galaxyModelType >= 0) {
    193             source->modelFits = psArrayAllocEmpty (1);
    194             pmModel *model = pmModelAlloc(galaxyModelType);
    195             psF32 *xPAR = model->params->data.F32;
    196 
    197             xPAR[PM_PAR_SKY]  = 0.0;
    198             xPAR[PM_PAR_I0]   = 1.0;
    199             xPAR[PM_PAR_XPOS] = X;
     238            // skip likely stars
     239            if (magDiff < starCut || SN < SNMinForCut) {
     240                if (forceAll) {
     241                    saveExtModelParams = true;
     242                } else if (galaxyModelType == modelTypeForce) {
     243                    // This is the model type that we are looking for
     244                    // proceed
     245                    saveExtModelParams = true;
     246                } else if (modelTypeForce == sersicModelType && galaxyModelType == devModelType) {
     247                    // We're doing sersic models, if sersic index is greater than the recipe's minimum index
     248                    // do dev model as well
     249                    if (isfinite(forceDevSersicMin) || Sindex >= forceDevSersicMin) {
     250                        saveExtModelParams = true;
     251                    }
     252                } else {
     253                    // not interested in this model
     254                }
     255            }
     256        }
     257
     258        if (saveExtModelParams) {
     259            if (!source->modelFits) {
     260                source->modelFits = psArrayAllocEmpty (1);
     261            }
     262            pmModel *model = pmModelAlloc(galaxyModelType);
     263            psF32 *xPAR = model->params->data.F32;
     264
     265            xPAR[PM_PAR_SKY]  = 0.0;
     266            xPAR[PM_PAR_I0]   = 1.0;
     267            xPAR[PM_PAR_XPOS] = X;
    200268            xPAR[PM_PAR_YPOS] = Y;
    201269           
    … …  
    210278            }
    211279
     280            model->chisq = chisq;
     281            model->nDOF = nDOF;
     282            model->flags = modelFlags;
     283
    212284            psArrayAdd (source->modelFits, 1, model);
    213285
    214 #ifdef notyet
    215             // XXX: set source->modelEXT to this model and flag as extended so that we have a better
    216             // shot at subtracting extended sources?
    217 
    218             // This doesn't work right. Need to do some more work to flesh out the model before it can be
    219             // used in subtaction. Better idea might be to do 2 passes in psphotFullForceReadout to
    220             // First find the best extended models and then do everything else.
    221             source->modelEXT = psMemIncrRefCounter(model);
    222             // is this safe? (no see above)
    223             source->type = PM_SOURCE_TYPE_EXTENDED;
    224             source->mode |= PM_SOURCE_MODE_EXTMODEL;
    225 #endif
    226 
    227286            psFree (model);
    228287        }
    229288
    230         sources->data[i] = source;
    231289        psFree(row);
     290    }
     291    if (source) {
     292        // close out last source
     293        psArrayAdd (sources, 1, source);
     294        psFree(source);
    232295    }
    233296
    … …  
    250313
    251314    pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC");
     315    pmModelType devModelType    = pmModelClassGetType("PS_MODEL_DEV");
     316    pmModelType selectedModelType = -1;
     317    bool chooseBest = false;
     318    bool chooseAll = false;
    252319
    253320    psString modelToChoose = psMetadataLookupStr(&mdok, recipe, "EXT_MODEL_TYPE_FOR_CFF");
    254     pmModelType selectedModelType = -1;
     321
    255322    if (mdok && modelToChoose != NULL) {
    256         if (strcmp(modelToChoose, "BEST")) {
     323        if (!strcmp(modelToChoose, "BEST")) {
     324            chooseBest = true;
     325        } else if (!strcmp(modelToChoose, "ALL")) {
     326            chooseAll = true;
     327        } else if (strcmp(modelToChoose, "PS_MODEL_SERSIC")) {
     328            // We have selected a model type other than Sersic.
     329            // Save it's type for use below.  Sersic is handled specially
    257330            selectedModelType = pmModelClassGetType(modelToChoose);
    258331        }
    259332    }
     333
     334    // minimum sersic index to force devModel
     335    psF32 sersicMinDev = psMetadataLookupF32(&mdok, recipe, "EXT_MODEL_FORCE_DEV_SERSIC_MIN");
     336    if (!mdok) {
     337        sersicMinDev = NAN;
     338    }
     339
     340    sources = psArraySort (sources, pmSourceSortBySeq);
    260341
    261342    for (int i = 0; i < sources->n; i++) {
    … …  
    263344        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    264345
    265         psF32 xPos, yPos, flux, rMajor, rMinor, theta;
    266         psS32 modelType = 0;
     346        #define MAX_ROWS_PER_SRC 10
     347        psF32 xPos[MAX_ROWS_PER_SRC], yPos[MAX_ROWS_PER_SRC], flux[MAX_ROWS_PER_SRC];
     348        psF32 rMajor[MAX_ROWS_PER_SRC], rMinor[MAX_ROWS_PER_SRC], theta[MAX_ROWS_PER_SRC];
     349        psF32 chisq[MAX_ROWS_PER_SRC], nDOF[MAX_ROWS_PER_SRC];
     350        psS32 modelFlags[MAX_ROWS_PER_SRC];
     351        psS32 modelType[MAX_ROWS_PER_SRC];
     352        psF32 sersicIndex = NAN;
    267353        bool fitGalaxy = false;
    268354        bool psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false;
    269         psF32 sersicIndex = 0;
    270         // For now only perform galaxy fits on extended objects
    271         if (source->modelEXT == NULL) {
    272             pmModel *model = source->modelPSF;
    273             if (model == NULL) continue;
    274             psF32 *PAR = model->params->data.F32;
    275             if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY])  || !isfinite(PAR[PM_PAR_SXY]) ||
    276                 !isfinite(source->psfFlux)) {
    277                 continue;
    278             }
    279 
    280             xPos = model->params->data.F32[PM_PAR_XPOS];
    281             yPos = model->params->data.F32[PM_PAR_YPOS];
    282             flux = source->psfFlux;
    283             rMajor = 0;
    284             rMinor = 0;
    285             theta = 0;
    286         } else {
    287             //   Find the model with the selected type. If selected type is -1 choose the one selected ad
    288             //   modelEXT which was the best
    289             int iModel = -1;
    290             if (source->modelEXT) {
    291                 pmModelType ext_model_type =  selectedModelType != -1  ? selectedModelType : source->modelEXT->type;
     355        int n_rows = 0;
     356
     357        psF32 kronFlux = source->moments->KronFlux;
     358        psF32 SN = NAN;
     359        psF32 magDiff = NAN;
     360        if (isfinite(kronFlux) && isfinite(source->moments->KronFluxErr) && isfinite(source->psfMag)) {
     361            SN = kronFlux/source->moments->KronFluxErr;
     362            // kronMag - psfMag for use as star/glaxy separator
     363            magDiff = -2.5 * log10(kronFlux) - source->psfMag ;
     364        }
     365
     366        // start with psf model
     367        pmModel *model = source->modelPSF;
     368        if (model == NULL) continue;
     369        psF32 *PAR = model->params->data.F32;
     370        if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY])  || !isfinite(PAR[PM_PAR_SXY]) ||
     371            !isfinite(source->psfFlux)) {
     372            continue;
     373        }
     374
     375        xPos[0] = model->params->data.F32[PM_PAR_XPOS];
     376        yPos[0] = model->params->data.F32[PM_PAR_YPOS];
     377        flux[0] = source->psfFlux;
     378        rMajor[0] = 0;
     379        rMinor[0] = 0;
     380        theta[0] = 0;
     381        modelType[0] = -1;
     382        sersicIndex = NAN;
     383        chisq[0] = NAN;
     384        nDOF[0] = NAN;
     385        modelFlags[0] = 0;
     386
     387        if (source->modelFits != NULL) {
     388            // figure out which models to use based on recipe paramters
     389            if (chooseAll) {
     390                // Save parameters for all valid extended models
     391
     392                // but make sure we aren't going to overflow our arrays
     393                assert (source->modelFits->n < MAX_ROWS_PER_SRC);
     394
    292395                for (int j=0; j<source->modelFits->n; j++) {
    293                     pmModel *aModel = source->modelFits->data[j];
    294                     if (aModel->type == ext_model_type) {
    295                         iModel = j;
    296                         break;
     396                    pmModel *model = source->modelFits->data[j];
     397                    psF32 *PAR = model->params->data.F32;
     398
     399                    if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY])  && isfinite(PAR[PM_PAR_SXY]) &&
     400                        isfinite(model->mag)) {
     401
     402                        xPos[n_rows] = PAR[PM_PAR_XPOS];
     403                        yPos[n_rows] = PAR[PM_PAR_YPOS];
     404
     405                        psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
     406                        rMajor[n_rows] = axes.major;
     407                        rMinor[n_rows] = axes.minor;
     408                        theta[n_rows]  = axes.theta*PS_DEG_RAD;
     409                        flux[n_rows] = pow(10.0, -0.4*model->mag);
     410                        modelType[n_rows] = model->type;
     411                        chisq[n_rows] = model->chisq;
     412                        nDOF[n_rows] = model->nDOF;
     413                        modelFlags[n_rows] = model->flags;
     414                        fitGalaxy = true;
     415                        if (model->params->n == 8) {
     416                            // this will save the sersic index in all models but that might
     417                            // be useful
     418                            sersicIndex = 0.5 / PAR[PM_PAR_7];
     419                        } else {
     420                            sersicIndex = NAN;
     421                        }
     422
     423                        n_rows++;
     424                    }
     425                }
     426            } else {
     427                int jModelSersic = -1;
     428                int jModelDev = -1;
     429                int jModelSelected = -1;
     430                psF32 minChisq = NAN;
     431                for (int j=0; j<source->modelFits->n; j++) {
     432                    pmModel *model = source->modelFits->data[j];
     433                    if (chooseBest) {
     434                        // choose the model with lowest chisq
     435                        if (isfinite(model->chisq) && (!isfinite(minChisq) || model->chisq < minChisq)) {
     436                            jModelSelected = j;
     437                            minChisq = model->chisq;
     438                        }
     439                    } else {
     440                        // find the index of models of interest
     441                        if (model->type == selectedModelType) {
     442                            jModelSelected = j;
     443                        } else if (model->type == sersicModelType) {
     444                            jModelSersic = j;
     445                        } else if (model->type == devModelType) {
     446                            jModelDev = j;
     447                        }
     448                    }
     449                }
     450                if (jModelSelected >= 0 || jModelSersic >= 0) {
     451                    // If a specific non-sersic model we take paramers from that one.
     452                    // Otherwise we do the sersic model.
     453                    pmModel *model = jModelSelected >= 0 ? source->modelFits->data[jModelSelected] :
     454                                                           source->modelFits->data[jModelSersic];
     455                    psF32 *PAR = model->params->data.F32;
     456
     457                    if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY])  && isfinite(PAR[PM_PAR_SXY]) &&
     458                        isfinite(model->mag)) {
     459
     460                        xPos[0] = PAR[PM_PAR_XPOS];
     461                        yPos[0] = PAR[PM_PAR_YPOS];
     462
     463                        psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
     464                        rMajor[0] = axes.major;
     465                        rMinor[0] = axes.minor;
     466                        theta[0]  = axes.theta*PS_DEG_RAD;
     467                        flux[0] = pow(10.0, -0.4*model->mag);
     468                        modelType[0] = model->type;
     469                        chisq[0] = model->chisq;
     470                        nDOF[0] = model->nDOF;
     471                        modelFlags[0] = model->flags;
     472                        fitGalaxy = true;
     473                        if (model->type == sersicModelType) {
     474                            PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false);
     475                            sersicIndex = 0.5 / PAR[PM_PAR_7];
     476                        }
     477
     478                        n_rows = 1;
     479
     480                        // Unless a specific non-sersic model type was selected do dev model for sources with
     481                        // sersic index above the recipe limit.
     482                        if (jModelSelected == -1 && jModelDev >= 0 && isfinite(sersicMinDev) &&
     483                            sersicIndex > sersicMinDev) {
     484
     485                            model = source->modelFits->data[jModelDev];
     486                            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SYY])  && isfinite(PAR[PM_PAR_SXY]) &&
     487                                isfinite(model->mag)) {
     488
     489                                xPos[1] = PAR[PM_PAR_XPOS];
     490                                yPos[1] = PAR[PM_PAR_YPOS];
     491
     492                                psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
     493                                rMajor[1] = axes.major;
     494                                rMinor[1] = axes.minor;
     495                                theta[1]  = axes.theta*PS_DEG_RAD;
     496                                flux[1] = pow(10.0, -0.4*model->mag);
     497                                modelType[1] = model->type;
     498                                sersicIndex = NAN;
     499                                chisq[1] = model->chisq;
     500                                nDOF[1] = model->nDOF;
     501                                modelFlags[1] = model->flags;
     502                                n_rows = 2;
     503                            }
     504                        }
    297505                    }
    298506                }
    299507            }
    300             if (iModel == -1) {
    301                 // Can this happen? perhaps if model type is qgauss (extended but below s/n for model fits?
    302                 // XXX: Should this be an assert?
    303                 continue;
    304             }
    305             pmModel *model = source->modelFits->data[iModel];
    306             psF32 *PAR = model->params->data.F32;
    307             xPos = PAR[PM_PAR_XPOS];
    308             yPos = PAR[PM_PAR_YPOS];
    309             if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) {
    310                 // XXX: do we need to handle this type
    311                 continue;
    312             } else {
    313                 if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY])  || !isfinite(PAR[PM_PAR_SXY]) ||
    314                     !isfinite(model->mag)) {
    315                     // bad model
    316                     continue;
    317                 }
    318                 psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
    319                 rMajor = axes.major;
    320                 rMinor = axes.minor;
    321                 theta  = axes.theta*PS_DEG_RAD;
    322                 flux = pow(10.0, -0.4*model->mag);
    323                 fitGalaxy = true;
    324             }
    325             modelType = model->type;
    326             if (modelType == sersicModelType) {
    327                 PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false);
    328                 sersicIndex = 0.5 / PAR[PM_PAR_7];
    329             }
    330         }
    331         psMetadata *row = psMetadataAlloc();
    332         psMetadataAddU32 (row, PS_LIST_TAIL, "ID",         0,   "IPP detection identifier",  source->seq);
    333         psMetadataAddF32 (row, PS_LIST_TAIL, "X",                0, "x coordinate",          xPos);
    334         psMetadataAddF32 (row, PS_LIST_TAIL, "Y",                0, "y coordinate",          yPos);
    335         psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX",             0, "flux per second",       flux/exptime);
    336         psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS",        0, "aperture radius",       source->apRadius);
    337         psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS",      0, "Kron radius",           source->moments->Mrf * 2.5);
    338         psMetadataAddF32 (row, PS_LIST_TAIL, "PETRO_RADIUS",     0, "Petrosian Radius",      source->extpars ? source->extpars->petrosianRadius : NAN);
    339         psMetadataAddBool (row, PS_LIST_TAIL, "FIT_GALAXY",      0, "source has xfit",       fitGalaxy);
    340         psMetadataAddBool (row, PS_LIST_TAIL, "PSF_STAR",        0, "source was psf star",   psfStar);
    341         psMetadataAddF32 (row, PS_LIST_TAIL, "R_MAJOR",          0, "radius of major axis",  rMajor);
    342         psMetadataAddF32 (row, PS_LIST_TAIL, "R_MINOR",          0, "radius of minor axis",  rMinor);
    343         psMetadataAddF32 (row, PS_LIST_TAIL, "THETA",            0, "theta",                 theta);
    344         psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "model type",            modelType);
    345         psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX",            0, "sersic index",          sersicIndex);
    346 
    347         psArrayAdd(table, 100, row);
    348         psFree(row);
     508        }
     509        if (n_rows == 0) {
     510            // didn't find a suitable extended model we have one row with the psf parameters
     511            n_rows = 1;
     512        }
     513
     514        for (int j = 0; j < n_rows; j++) {
     515            psMetadata *row = psMetadataAlloc();
     516            psMetadataAddU32 (row, PS_LIST_TAIL, "ID",         0,   "IPP detection identifier",  source->seq);
     517            psMetadataAddF32 (row, PS_LIST_TAIL, "X",                0, "x coordinate",          xPos[j]);
     518            psMetadataAddF32 (row, PS_LIST_TAIL, "Y",                0, "y coordinate",          yPos[j]);
     519            psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX",             0, "flux per second",       flux[j]/exptime);
     520            psMetadataAddF32 (row, PS_LIST_TAIL, "SN",               0, "kron flux signal to noise", SN);
     521            psMetadataAddF32 (row, PS_LIST_TAIL, "MAG_DIFF",         0, "psf mag - kron mag",    magDiff);
     522            psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS",        0, "aperture radius",       source->apRadius);
     523            psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS",      0, "Kron radius",           source->moments->Mrf * 2.5);
     524            psMetadataAddF32 (row, PS_LIST_TAIL, "PETRO_RADIUS",     0, "Petrosian Radius",      source->extpars ? source->extpars->petrosianRadius : NAN);
     525            psMetadataAddBool (row, PS_LIST_TAIL, "FIT_GALAXY",      0, "source has xfit",       fitGalaxy);
     526            psMetadataAddBool (row, PS_LIST_TAIL, "PSF_STAR",        0, "source was psf star",   psfStar);
     527            psMetadataAddF32 (row, PS_LIST_TAIL, "R_MAJOR",          0, "radius of major axis",  rMajor[j]);
     528            psMetadataAddF32 (row, PS_LIST_TAIL, "R_MINOR",          0, "radius of minor axis",  rMinor[j]);
     529            psMetadataAddF32 (row, PS_LIST_TAIL, "THETA",            0, "theta",                 theta[j]);
     530            psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "model type",            modelType[j]);
     531            psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX",            0, "sersic index",          sersicIndex);
     532            psMetadataAddF32 (row, PS_LIST_TAIL, "CHISQ",            0, "chisq",                 chisq[j]);
     533            psMetadataAddF32 (row, PS_LIST_TAIL, "NDOF",             0, "n degrees of freedom",  nDOF[j]);
     534            psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_FLAGS",      0, "model flags",           modelFlags[j]);
     535
     536            psArrayAdd(table, 100, row);
     537            psFree(row);
     538        }
    349539    }
    350540

  • branches/eam_branches/ipp-20140423/psModules/src/objects/pmSourceIO_CMF.c.in

    r36623 r36829  
    482482    }
    483483
     484#ifdef SORT_EXTENSIONS_BY_FLUX
    484485    // let's write these out in S/N order
    485486    sources = psArraySort (sources, pmSourceSortByFlux);
     487#endif
    486488
    487489    table = psArrayAllocEmpty (sources->n);
    … …  
    821823    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname);
    822824
     825#ifdef SORT_EXTENSIONS_BY_FLUX
    823826    // let's write these out in S/N order
    824827    sources = psArraySort (sources, pmSourceSortByFlux);
     828#endif
    825829
    826830    float magOffset;
    … …  
    942946
    943947            // EAM : adding for PV2 outputs:
    944             @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", source->modelEXT ? source->modelEXT->flags : 0);
     948            @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags);
    945949
    946950            @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE",   0, "position angle at source (degrees)",         posAngle);
    … …  
    11091113        model->chisq = psMetadataLookupF32(&status, row, "EXT_CHISQ");
    11101114        model->nDOF = psMetadataLookupF32(&status, row, "EXT_NDOF");
     1115        model->flags = psMetadataLookupS16(&status, row, "EXT_FLAGS");
    11111116
    11121117        // EXT_MODEL_TYPE gives the model chosen by psphot as the best.
    … …  
    12251230    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
    12261231
     1232#ifdef SORT_EXTENSIONS_BY_FLUX
    12271233    // let's write these out in S/N order
    12281234    sources = psArraySort (sources, pmSourceSortByFlux);
     1235#endif
    12291236
    12301237    table = psArrayAllocEmpty (sources->n);
    … …  
    14391446    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "galaxy table extension", extname);
    14401447
    1441     psMetadataAddStr (outhead, PS_LIST_TAIL, "HI", PS_META_REPLACE, "does this get through?", "THERE");
    1442 
    1443     // let's write these out in S/N order
    1444     sources = psArraySort (sources, pmSourceSortByFlux);
     1448    psF32 Q = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_Q");
     1449    psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_Q", PS_META_REPLACE, "", Q);
     1450
     1451    psF32 NSigma = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_NSIGMA");
     1452    psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_NSIGMA", PS_META_REPLACE, "", NSigma);
     1453
     1454    psF32 clampSN = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_CLAMP_SN");
     1455    psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_CLAMP_SN", PS_META_REPLACE, "", clampSN);
     1456
     1457    // They are probably already in this order but ...
     1458    sources = psArraySort (sources, pmSourceSortBySeq);
    14451459
    14461460    psArray *table = psArrayAllocEmpty (sources->n);
    … …  
    14591473        if (source->galaxyFits == NULL) continue;
    14601474
    1461         pmModel *model = source->modelFits->data[0];
    1462         if (!model) return false;
    1463 
    1464         // X,Y coordinates are stored with the model parameters
    1465         psF32 *PAR = model->params->data.F32;
    1466 
    1467         psMetadata *row = psMetadataAlloc ();
    1468 
    1469         // we write out the x,y positions so people can link to the psf either way (position or ID)
    1470         psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index", source->seq);
    1471         psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT",            0, "model x coordinate",             PAR[PM_PAR_XPOS]);
    1472         psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT",            0, "model y coordinate",             PAR[PM_PAR_YPOS]);
    1473         psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX",             0, "number of pixels for fits",      source->galaxyFits->nPix);
    1474 
    1475         psVector *Flux = source->galaxyFits->Flux;
    1476         psVector *dFlux = source->galaxyFits->dFlux;
    1477         psVector *chisq = source->galaxyFits->chisq;
    1478 
    1479         psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX",     PS_META_REPLACE, "normalization for galaxy flux", Flux);
    1480         psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX_ERR", PS_META_REPLACE, "error on normalization", dFlux);
    1481         psMetadataAddVector (row, PS_LIST_TAIL, "GAL_CHISQ",    PS_META_REPLACE, "galaxy fit chisq", chisq);
    1482 
    1483         psArrayAdd (table, 100, row);
    1484         psFree (row);
     1475        for (int iModel = 0; iModel < source->modelFits->n; iModel++) {
     1476            pmModel *model = source->modelFits->data[iModel];
     1477            if (!model) continue;
     1478
     1479            pmSourceGalaxyFits *galaxyFits = NULL;
     1480            for (int iFit = 0; iFit < source->galaxyFits->n; iFit++) {
     1481                galaxyFits = source->galaxyFits->data[iFit];
     1482                if (model->type == galaxyFits->modelType) break;
     1483                galaxyFits = NULL;
     1484            }
     1485
     1486            if (!galaxyFits) continue;
     1487
     1488            // X,Y coordinates are stored with the model parameters
     1489            psF32 *PAR = model->params->data.F32;
     1490
     1491            psMetadata *row = psMetadataAlloc ();
     1492
     1493            // we write out the x,y positions so people can link to the psf either way (position or ID)
     1494            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index", source->seq);
     1495            psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "model type",                     galaxyFits->modelType);
     1496            psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT",            0, "model x coordinate",             PAR[PM_PAR_XPOS]);
     1497            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT",            0, "model y coordinate",             PAR[PM_PAR_YPOS]);
     1498            psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX",             0, "number of pixels for fits",      galaxyFits->nPix);
     1499            // psMetadataAddS32 (row, PS_LIST_TAIL, "REDUCED_TRIALS",   0, "source has reduced number of trials",      galaxyFits->reducedTrials);
     1500
     1501            psVector *Flux = galaxyFits->Flux;
     1502            psVector *dFlux = galaxyFits->dFlux;
     1503            psVector *chisq = galaxyFits->chisq;
     1504
     1505            psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX",     PS_META_REPLACE, "normalization for galaxy flux", Flux);
     1506            psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX_ERR", PS_META_REPLACE, "error on normalization", dFlux);
     1507            psMetadataAddVector (row, PS_LIST_TAIL, "GAL_CHISQ",    PS_META_REPLACE, "galaxy fit chisq", chisq);
     1508
     1509            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_MIN",    0, "fractional major axis min",      galaxyFits->fRmajorMin);
     1510            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_MAX",    0, "fractional major axis max",      galaxyFits->fRmajorMax);
     1511            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_DEL",    0, "fractional major axis max",      galaxyFits->fRmajorDel);
     1512            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_MIN",    0, "fractional minor axis min",      galaxyFits->fRminorMin);
     1513            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_MAX",    0, "fractional minor axis max",      galaxyFits->fRminorMax);
     1514            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_DEL",    0, "fractional minor axis max",      galaxyFits->fRminorDel);
     1515
     1516            psArrayAdd (table, 100, row);
     1517            psFree (row);
     1518        }
    14851519    }
    14861520
    … …  
    15201554        return false;
    15211555    }
     1556    if (!readout->analysis) {
     1557        readout->analysis = psMetadataAlloc();
     1558    }
     1559    psF32 Q = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_Q");
     1560    // XXX: turn this into an assert once we're done
     1561    if (status) {
     1562        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_Q", PS_META_REPLACE, "", Q);
     1563
     1564        psF32 NSigma = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_NSIGMA");
     1565        psAssert(status, "missing GALAXY_SHAPES_NSIGMA");
     1566        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_NSIGMA", PS_META_REPLACE, "", NSigma);
     1567
     1568        psF32 clampSN = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_CLAMP_SN");
     1569        psAssert(status, "missing GALAXY_SHAPES_CLAMP_SN");
     1570        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_CLAMP_SN", PS_META_REPLACE, "", clampSN);
     1571    }
    15221572
    15231573    for (long i = 0; i < numSources; i++) {
    … …  
    15411591        }
    15421592
     1593        int modelType = psMetadataLookupS32 (&status,    row, "MODEL_TYPE");
    15431594        psVector *Flux  = psMetadataLookupVector(&status, row, "GAL_FLUX");
    15441595        psVector *dFlux = psMetadataLookupVector(&status, row, "GAL_FLUX_ERR");
    … …  
    15461597
    15471598        if (Flux && Flux->n > 0) {
    1548             psFree(source->galaxyFits);
    1549             source->galaxyFits = pmSourceGalaxyFitsAlloc();
    1550             source->galaxyFits->nPix = psMetadataLookupF32(&status, row, "NPIX");
    1551 
    1552             psFree(source->galaxyFits->Flux);
    1553             source->galaxyFits->Flux  = psMemIncrRefCounter(Flux);
    1554             psFree(source->galaxyFits->dFlux);
    1555             source->galaxyFits->dFlux = psMemIncrRefCounter(dFlux);
    1556             psFree(source->galaxyFits->chisq);
    1557             source->galaxyFits->chisq = psMemIncrRefCounter(chisq);
     1599            if (!source->galaxyFits) {
     1600                source->galaxyFits = psArrayAllocEmpty(1);
     1601            }
     1602
     1603            pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc();
     1604
     1605            psArrayAdd(source->galaxyFits, 1, galaxyFits);
     1606
     1607            psFree(galaxyFits);
     1608            galaxyFits->modelType = modelType;
     1609            galaxyFits->nPix = psMetadataLookupF32(&status, row, "NPIX");
     1610
     1611            galaxyFits->fRmajorMin = psMetadataLookupF32(&status, row, "FR_MAJOR_MIN");
     1612            galaxyFits->fRmajorMax = psMetadataLookupF32(&status, row, "FR_MAJOR_MAX");
     1613            galaxyFits->fRmajorDel = psMetadataLookupF32(&status, row, "FR_MAJOR_DEL");
     1614            galaxyFits->fRminorMin = psMetadataLookupF32(&status, row, "FR_MINOR_MIN");
     1615            galaxyFits->fRminorMax = psMetadataLookupF32(&status, row, "FR_MINOR_MAX");
     1616            galaxyFits->fRminorDel = psMetadataLookupF32(&status, row, "FR_MINOR_DEL");
     1617
     1618            psFree(galaxyFits->Flux);
     1619            galaxyFits->Flux  = psMemIncrRefCounter(Flux);
     1620            psFree(galaxyFits->dFlux);
     1621            galaxyFits->dFlux = psMemIncrRefCounter(dFlux);
     1622            psFree(galaxyFits->chisq);
     1623            galaxyFits->chisq = psMemIncrRefCounter(chisq);
    15581624        }
    15591625
Note: See TracChangeset for help on using the changeset viewer.