IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5762


Ignore:
Timestamp:
Dec 12, 2005, 10:32:44 AM (21 years ago)
Author:
desonia
Message:

renamed the variable 'try' to 'psfTry' to avoid the conflict of keywords in
C++.

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

Legend:

Unmodified
Added
Removed
  • trunk/psModules/src/objects/pmPSFtry.c

    r5255 r5762  
    7979    int Npsf = 0;
    8080
    81     pmPSFtry *try
    82     = pmPSFtryAlloc (sources, modelName);
     81    pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName);
    8382
    8483    // stage 1:  fit an independent model (freeModel) to all sources
    8584    psTimerStart ("fit");
    86     for (int i = 0; i < try
    87                 ->sources->n; i++) {
    88 
    89             pmSource *source = try
    90                                    ->sources->data[i];
    91             pmModel  *model  = pmSourceModelGuess (source, try
    92                                                        ->psf->type);
    93             x = source->peak->x;
    94             y = source->peak->y;
    95 
    96             // set temporary object mask and fit object
    97             // fit model as FLT, not PSF
    98             psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
    99             status = pmSourceFitModel (source, model, false);
    100             psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
    101 
    102             // exclude the poor fits
    103             if (!status) {
    104                 try
    105                     ->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
    106                 psFree (model);
    107                 continue;
    108             }
    109             try
    110                 ->modelFLT->data[i] = model;
    111             Nflt ++;
    112         }
     85    for (int i = 0; i < psfTry->sources->n; i++) {
     86
     87        pmSource *source = psfTry->sources->data[i];
     88        pmModel  *model  = pmSourceModelGuess (source, psfTry->psf->type);
     89        x = source->peak->x;
     90        y = source->peak->y;
     91
     92        // set temporary object mask and fit object
     93        // fit model as FLT, not PSF
     94        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
     95        status = pmSourceFitModel (source, model, false);
     96        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
     97
     98        // exclude the poor fits
     99        if (!status) {
     100            psfTry->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
     101            psFree (model);
     102            continue;
     103        }
     104        psfTry->modelFLT->data[i] = model;
     105        Nflt ++;
     106    }
    113107    psLogMsg ("psphot.psftry", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    114108    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
    115109
    116110    // make this optional?
    117     // DumpModelFits (try->modelFLT, "modelsFLT.dat");
     111    // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat");
    118112
    119113    // stage 2: construct a psf (pmPSF) from this collection of model fits
    120     pmPSFFromModels (try
    121                      ->psf, try
    122                          ->modelFLT, try
    123                              ->mask);
     114    pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask);
    124115
    125116    // stage 3: refit with fixed shape parameters
    126117    psTimerStart ("fit");
    127     for (int i = 0; i < try
    128                 ->sources->n; i++) {
    129             // masked for: bad model fit, outlier in parameters
    130             if (try
    131                     ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    132 
    133             pmSource *source = try
    134                                    ->sources->data[i];
    135             pmModel  *modelFLT = try
    136                                      ->modelFLT->data[i];
    137 
    138             // set shape for this model based on PSF
    139             pmModel *modelPSF = pmModelFromPSF (modelFLT, try
    140                                                     ->psf);
    141             x = source->peak->x;
    142             y = source->peak->y;
    143 
    144             psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
    145             status = pmSourceFitModel (source, modelPSF, true);
    146 
    147             // skip poor fits
    148             if (!status) {
    149                 try
    150                     ->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL;
    151                 psFree (modelPSF);
    152                 goto next_source;
    153             }
    154 
    155             // otherwise, save the resulting model
    156             try
    157                 ->modelPSF->data[i] = modelPSF;
    158 
    159             // XXX : use a different aperture radius from the fit radius?
    160             // XXX : use a different estimator for the local sky?
    161             // XXX : pass 'source' as input?
    162             if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
    163                 try
    164                     ->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
    165                 goto next_source;
    166             }
    167 
    168             try
    169                 ->metric->data.F64[i] = obsMag - fitMag;
    170             try
    171                 ->fitMag->data.F64[i] = fitMag;
    172             Npsf ++;
     118    for (int i = 0; i < psfTry->sources->n; i++) {
     119        // masked for: bad model fit, outlier in parameters
     120        if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL)
     121            continue;
     122
     123        pmSource *source = psfTry->sources->data[i];
     124        pmModel  *modelFLT = psfTry->modelFLT->data[i];
     125
     126        // set shape for this model based on PSF
     127        pmModel *modelPSF = pmModelFromPSF (modelFLT, psfTry->psf);
     128        x = source->peak->x;
     129        y = source->peak->y;
     130
     131        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
     132        status = pmSourceFitModel (source, modelPSF, true);
     133
     134        // skip poor fits
     135        if (!status) {
     136            psfTry->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL;
     137            psFree (modelPSF);
     138            goto next_source;
     139        }
     140
     141        // otherwise, save the resulting model
     142        psfTry->modelPSF->data[i] = modelPSF;
     143
     144        // XXX : use a different aperture radius from the fit radius?
     145        // XXX : use a different estimator for the local sky?
     146        // XXX : pass 'source' as input?
     147        if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
     148            psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
     149            goto next_source;
     150        }
     151
     152        psfTry->metric->data.F64[i] = obsMag - fitMag;
     153        psfTry->fitMag->data.F64[i] = fitMag;
     154        Npsf ++;
    173155
    174156next_source:
    175             psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
    176 
    177         }
     157        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
     158
     159    }
    178160    psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    179161    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
    180162
    181163    // make this optional
    182     // DumpModelFits (try->modelPSF, "modelsPSF.dat");
     164    // DumpModelFits (psfTry->modelPSF, "modelsPSF.dat");
    183165
    184166    // XXX this function wants aperture radius for pmSourcePhotometry
    185     pmPSFtryMetric (try
    186                     , RADIUS);
     167    pmPSFtryMetric (psfTry, RADIUS);
    187168    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
    188               modelName, try
    189                   ->psf->ApResid, try
    190                       ->psf->dApResid, try
    191                           ->psf->skyBias);
    192 
    193     return (try
    194            );
     169              modelName,
     170              psfTry->psf->ApResid,
     171              psfTry->psf->dApResid,
     172              psfTry->psf->skyBias);
     173
     174    return (psfTry);
    195175}
    196176
    197177
    198 bool pmPSFtryMetric (pmPSFtry *try
    199                      , float RADIUS)
     178bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS)
    200179{
    201180
     
    203182    int   nKeep, nSkip;
    204183
    205     // the measured (aperture - fit) magnitudes (dA == try->metric)
     184    // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
    206185    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    207186    //     dA = dAo + dsky/flux
     
    213192
    214193    // rflux = ten(0.4*fitMag);
    215     psVector *rflux = psVectorAlloc (try
    216                                      ->sources->n, PS_TYPE_F64);
    217     for (int i = 0; i < try
    218                 ->sources->n; i++) {
    219             if (try
    220                     ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    221             rflux->data.F64[i] = pow(10.0, 0.4*try
    222                                      ->fitMag->data.F64[i]);
    223         }
     194    psVector *rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     195    for (int i = 0; i < psfTry->sources->n; i++) {
     196        if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL)
     197            continue;
     198        rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]);
     199    }
    224200
    225201    // find min and max of (1/flux):
    226202    psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    227     psVectorStats (stats, rflux, NULL, try
    228                        ->mask, PSFTRY_MASK_ALL);
     203    psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL);
    229204
    230205    // build binned versions of rflux, metric
     
    240215    for (int i = 0; i < daBin->n; i++) {
    241216
    242         psVector *tmp = psVectorAlloc (try
    243                                        ->sources->n, PS_TYPE_F64);
     217        psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    244218        tmp->n = 0;
    245219
    246220        // accumulate data within bin range
    247         for (int j = 0; j < try
    248                     ->sources->n; j++) {
    249                 // masked for: bad model fit, outlier in parameters
    250                 if (try
    251                         ->mask->data.U8[j] & PSFTRY_MASK_ALL) continue;
    252 
    253                 // skip points with extreme dA values
    254                 if (fabs(try
    255                          ->metric->data.F64[j]) > 0.5) continue;
    256 
    257                 // skip points outside of this bin
    258                 if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin)
    259                     continue;
    260                 if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin)
    261                     continue;
    262 
    263                 tmp->data.F64[tmp->n] = try
    264                                             ->metric->data.F64[j];
    265                 tmp->n ++;
    266             }
     221        for (int j = 0; j < psfTry->sources->n; j++) {
     222            // masked for: bad model fit, outlier in parameters
     223            if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL)
     224                continue;
     225
     226            // skip points with extreme dA values
     227            if (fabs(psfTry->metric->data.F64[j]) > 0.5)
     228                continue;
     229
     230            // skip points outside of this bin
     231            if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin)
     232                continue;
     233            if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin)
     234                continue;
     235
     236            tmp->data.F64[tmp->n] = psfTry->metric->data.F64[j];
     237            tmp->n ++;
     238        }
    267239
    268240        // is this a valid point?
     
    319291    stats = psVectorStats (stats, daResid, NULL, maskB, 1);
    320292
    321     try
    322         ->psf->ApResid = poly->coeff[0];
    323     try
    324         ->psf->dApResid = stats->clippedStdev;
    325     try
    326         ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
     293    psfTry->psf->ApResid = poly->coeff[0];
     294    psfTry->psf->dApResid = stats->clippedStdev;
     295    psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    327296
    328297    psFree (rflux);
  • trunk/psModules/src/objects/pmPSFtry.h

    r5255 r5762  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-10-10 19:53:40 $
     8 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-12-12 20:32:44 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    105105 */
    106106bool pmPSFtryMetric(
    107     pmPSFtry *try
    108     ,                      ///< Add comment.
    109     float RADIUS                        ///< Add comment.
     107    pmPSFtry *psfTry,                  ///< Add comment.
     108    float RADIUS                       ///< Add comment.
    110109);
    111110
Note: See TracChangeset for help on using the changeset viewer.