IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 19, 2007, 4:40:44 PM (19 years ago)
Author:
Paul Price
Message:

Extensive changes to APIs to allow use of a nominated value to mask
against (the maskVal). Previously, the mask values were either
hard-coded (e.g., PM_MASK_SAT) or taken as anything non-zero. The
code is tested with psModules (which has similar changes) and does not
crash, but neither is it successful in marking all bad pixels. For
this reason, I have left the "gutter" pixels (cell gaps) set to 0
instead of NAN in pmFPAMosaic.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotChoosePSF.c

    r13834 r13900  
    22
    33// try PSF models and select best option
    4 pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) {
     4pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal, psMaskType mark) {
    55
    66    bool status;
    … …  
    4343    psPolynomial2D *psfTrendMask = psPolynomial2DfromMetadata (md);
    4444    if (!psfTrendMask) {
    45         psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe");
    46         return NULL;
     45        psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe");
     46        return NULL;
    4747    }
    4848
    … …  
    5353        pmSource *source = sources->data[i];
    5454        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
    55             // keep NSTARS PSF stars, unmark the rest
    56             if (stars->n < NSTARS) {
    57                 psArrayAdd (stars, 200, source);
    58             } else {
    59                 source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
    60             }
    61         }
     55            // keep NSTARS PSF stars, unmark the rest
     56            if (stars->n < NSTARS) {
     57                psArrayAdd (stars, 200, source);
     58            } else {
     59                source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
     60            }
     61        }
    6262    }
    6363    psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n);
    … …  
    6565    if (stars->n == 0) {
    6666        psLogMsg ("psphot.choosePSF", PS_LOG_WARN, "Failed to find any PSF candidates");
    67         psFree (stars);
    68         psFree (psfTrendMask);
     67        psFree (stars);
     68        psFree (psfTrendMask);
    6969        return NULL;
    7070    }
    … …  
    9191        psMetadataItem *item = psListGetAndIncrement (iter);
    9292        char *modelName = item->data.V;
    93         models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS);
     93        models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS, maskVal, mark);
    9494    }
    9595
    … …  
    127127    // print/dump psf parameters
    128128    if (psTraceGetLevel("psphot") >= 5) {
    129         for (int i = PM_PAR_SXX; i < try->psf->params_NEW->n; i++) {
    130             psPolynomial2D *poly = try->psf->params_NEW->data[i];
    131             for (int nx = 0; nx <= poly->nX; nx++) {
    132                 for (int ny = 0; ny <= poly->nY; ny++) {
    133                     if (poly->mask[nx][ny]) continue;
    134                     fprintf (stderr, "%g x^%d y^%d\n", poly->coeff[nx][ny], nx, ny);
    135                 }
    136             }
    137         }
     129        for (int i = PM_PAR_SXX; i < try->psf->params_NEW->n; i++) {
     130            psPolynomial2D *poly = try->psf->params_NEW->data[i];
     131            for (int nx = 0; nx <= poly->nX; nx++) {
     132                for (int ny = 0; ny <= poly->nY; ny++) {
     133                    if (poly->mask[nx][ny]) continue;
     134                    fprintf (stderr, "%g x^%d y^%d\n", poly->coeff[nx][ny], nx, ny);
     135                }
     136            }
     137        }
    138138    }
    139139
    … …  
    147147    psVector *dSN = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
    148148    for (int i = 0; i < try->sources->n; i++) {
    149         // masked for: bad model fit, outlier in parameters
    150         if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
    151             continue;
    152 
    153         pmSource *source = try->sources->data[i];
    154         Sx->data.F32[Sx->n] = source->modelPSF->params->data.F32[PM_PAR_SXX];
    155         Sy->data.F32[Sy->n] = source->modelPSF->params->data.F32[PM_PAR_SYY];
    156         dSN->data.F32[dSN->n] = source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0];
    157         dSx->data.F32[dSx->n] = source->modelPSF->dparams->data.F32[PM_PAR_SXX];
    158         dSy->data.F32[dSy->n] = source->modelPSF->dparams->data.F32[PM_PAR_SYY];
    159         Sx->n ++;
    160         Sy->n ++;
    161         dSN->n ++;
    162         dSx->n ++;
    163         dSy->n ++;
     149        // masked for: bad model fit, outlier in parameters
     150        if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
     151            continue;
     152
     153        pmSource *source = try->sources->data[i];
     154        Sx->data.F32[Sx->n] = source->modelPSF->params->data.F32[PM_PAR_SXX];
     155        Sy->data.F32[Sy->n] = source->modelPSF->params->data.F32[PM_PAR_SYY];
     156        dSN->data.F32[dSN->n] = source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0];
     157        dSx->data.F32[dSx->n] = source->modelPSF->dparams->data.F32[PM_PAR_SXX];
     158        dSy->data.F32[dSy->n] = source->modelPSF->dparams->data.F32[PM_PAR_SYY];
     159        Sx->n ++;
     160        Sy->n ++;
     161        dSN->n ++;
     162        dSx->n ++;
     163        dSy->n ++;
    164164    }
    165165    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    … …  
    174174
    175175    for (int i = 0; i < dSN->n; i++) {
    176         dSx->data.F32[i] = (dSx->data.F32[i] - dSxo) / PS_MAX (PSF_MIN_DS, Sx->data.F32[i]*dSN->data.F32[i]);
    177         dSy->data.F32[i] = (dSy->data.F32[i] - dSyo) / PS_MAX (PSF_MIN_DS, Sy->data.F32[i]*dSN->data.F32[i]);
     176        dSx->data.F32[i] = (dSx->data.F32[i] - dSxo) / PS_MAX (PSF_MIN_DS, Sx->data.F32[i]*dSN->data.F32[i]);
     177        dSy->data.F32[i] = (dSy->data.F32[i] - dSyo) / PS_MAX (PSF_MIN_DS, Sy->data.F32[i]*dSN->data.F32[i]);
    178178    }
    179179
    … …  
    206206
    207207    // build a PSF residual image
    208     if (!psphotMakeResiduals (try->sources, recipe, try->psf)) {
    209         psError(PSPHOT_ERR_PSF, false, "Unable to construct residual table for PSF");
    210         psFree (models);
    211         return NULL;
     208    if (!psphotMakeResiduals (try->sources, recipe, try->psf, maskVal)) {
     209        psError(PSPHOT_ERR_PSF, false, "Unable to construct residual table for PSF");
     210        psFree (models);
     211        return NULL;
    212212    }
    213213
    214214    // XXX test dump of psf star data and psf-subtracted image
    215     if (psTraceGetLevel("psphot.psfstars") > 5) { 
    216         psphotSaveImage (NULL, readout->image,  "rawstars.fits");
    217 
    218         for (int i = 0; i < try->sources->n; i++) {
    219             // masked for: bad model fit, outlier in parameters
    220             if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
    221                 continue;
    222 
    223             pmSource *source = try->sources->data[i];
    224             float x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    225             float y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
    226 
    227             // set the mask and subtract the PSF model
    228             // XXX should we be using maskObj? should we be unsetting the mask?
    229             // use pmModelSub because modelFlux has not been generated
    230             assert (source->maskObj);
    231             psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK);
    232             pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL);
    233             psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
    234         }
    235 
    236         FILE *f = fopen ("shapes.dat", "w");
    237         for (int i = 0; i < try->sources->n; i++) {
    238             psF32 inPar[10];
    239 
    240             // masked for: bad model fit, outlier in parameters
    241             if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    242 
    243             pmSource *source = try->sources->data[i];
    244             psF32 *outPar = source->modelEXT->params->data.F32;
    245 
    246             psEllipseShape shape;
    247 
    248             shape.sx  = outPar[PM_PAR_SXX] / M_SQRT2;
    249             shape.sy  = outPar[PM_PAR_SYY] / M_SQRT2;
    250             shape.sxy = outPar[PM_PAR_SXY];
    251 
    252             psEllipsePol pol = pmPSF_ModelToFit (outPar);
    253             inPar[PM_PAR_E0] = pol.e0;
    254             inPar[PM_PAR_E1] = pol.e1;
    255             inPar[PM_PAR_E2] = pol.e2;
    256             pmPSF_FitToModel (inPar, 0.1);
    257 
    258             psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0);
    259             psEllipseAxes axes2;
    260             (void)psEllipsePolToAxes(pol, 0.1, &axes2);
    261             psEllipsePol pol2 = psEllipseAxesToPol (axes1);
    262 
    263             fprintf (f, "%3d  %7.2f %7.2f  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f  :  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n",
    264                      i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS],
    265                      outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY],
    266                      pol.e0, pol.e1, pol.e2,
    267                      pol2.e0, pol2.e1, pol2.e2,
    268                      inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY],
    269                      axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD,
    270                      axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD
    271                      );
    272         }
    273         fclose (f);
    274 
    275         psphotSaveImage (NULL, readout->image,  "psfstars.fits");
    276         pmSourcesWritePSFs (try->sources, "psfstars.dat");
    277         pmSourcesWriteEXTs (try->sources, "extstars.dat", false);
    278         psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
    279         psMetadataConfigWrite (psfData, "psfmodel.dat");
    280         psFree (psfData);
    281         psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n");
    282         exit (0);
     215    if (psTraceGetLevel("psphot.psfstars") > 5) {
     216        psphotSaveImage (NULL, readout->image,  "rawstars.fits");
     217
     218        for (int i = 0; i < try->sources->n; i++) {
     219            // masked for: bad model fit, outlier in parameters
     220            if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
     221                continue;
     222
     223            pmSource *source = try->sources->data[i];
     224            float x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
     225            float y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
     226
     227            // set the mask and subtract the PSF model
     228            // XXX should we be using maskObj? should we be unsetting the mask?
     229            // use pmModelSub because modelFlux has not been generated
     230            assert (source->maskObj);
     231            psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK);
     232            pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal);
     233            psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
     234        }
     235
     236        FILE *f = fopen ("shapes.dat", "w");
     237        for (int i = 0; i < try->sources->n; i++) {
     238            psF32 inPar[10];
     239
     240            // masked for: bad model fit, outlier in parameters
     241            if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
     242
     243            pmSource *source = try->sources->data[i];
     244            psF32 *outPar = source->modelEXT->params->data.F32;
     245
     246            psEllipseShape shape;
     247
     248            shape.sx  = outPar[PM_PAR_SXX] / M_SQRT2;
     249            shape.sy  = outPar[PM_PAR_SYY] / M_SQRT2;
     250            shape.sxy = outPar[PM_PAR_SXY];
     251
     252            psEllipsePol pol = pmPSF_ModelToFit (outPar);
     253            inPar[PM_PAR_E0] = pol.e0;
     254            inPar[PM_PAR_E1] = pol.e1;
     255            inPar[PM_PAR_E2] = pol.e2;
     256            pmPSF_FitToModel (inPar, 0.1);
     257
     258            psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0);
     259            psEllipseAxes axes2;
     260            (void)psEllipsePolToAxes(pol, 0.1, &axes2);
     261            psEllipsePol pol2 = psEllipseAxesToPol (axes1);
     262
     263            fprintf (f, "%3d  %7.2f %7.2f  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f  :  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n",
     264                     i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS],
     265                     outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY],
     266                     pol.e0, pol.e1, pol.e2,
     267                     pol2.e0, pol2.e1, pol2.e2,
     268                     inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY],
     269                     axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD,
     270                     axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD
     271                     );
     272        }
     273        fclose (f);
     274
     275        psphotSaveImage (NULL, readout->image,  "psfstars.fits");
     276        pmSourcesWritePSFs (try->sources, "psfstars.dat");
     277        pmSourcesWriteEXTs (try->sources, "extstars.dat", false);
     278        psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
     279        psMetadataConfigWrite (psfData, "psfmodel.dat");
     280        psFree (psfData);
     281        psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n");
     282        exit (0);
    283283    }
    284284
    … …  
    290290
    291291    if (!psphotPSFstats (readout, recipe, psf)) {
    292         psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms");
    293         psFree(psf);
    294         return NULL;
     292        psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms");
     293        psFree(psf);
     294        return NULL;
    295295    }
    296296
    … …  
    327327    pmModel *modelPSF = pmModelFromPSF (modelEXT, psf);
    328328    if (modelPSF == NULL) {
    329         psError(PSPHOT_ERR_PSF, false, "Failed to estimate PSF model at image centre");
    330         psFree(modelEXT);
    331         return false;
     329        psError(PSPHOT_ERR_PSF, false, "Failed to estimate PSF model at image centre");
     330        psFree(modelEXT);
     331        return false;
    332332    }
    333333
    … …  
    346346    psF64 FWHM_Y = FWHM_X * (axes.minor / axes.major);
    347347
    348     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_X",   PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X);
    349     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_Y",   PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y);
    350     psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
     348    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_X",   PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X);
     349    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_Y",   PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y);
     350    psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
    351351    psMetadataAddS32 (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars);
    352352    psMetadataAddBool(recipe, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true);
    … …  
    386386        moments.xy = source->moments->Sxy;
    387387
    388         // limit axis ratio < 20.0
     388        // limit axis ratio < 20.0
    389389        axes = psEllipseMomentsToAxes (moments, 20.0);
    390390
Note: See TracChangeset for help on using the changeset viewer.