IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 9, 2025, 11:21:40 AM (15 months ago)
Author:
eugene
Message:

add psphotFindFeatures; enable forced-photometry between cameras; do NOT skip negative-flux sources in full-force photometry; make PSF residual component generation more robust: IRLS vs median; exclude pixels with poor values (scaled value > 1.5, relative to normalized PSF); exclude pixels with too few inputs (<5); allow PSF clump to use clipped stats or robust stats

Location:
trunk/psphot
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot

  • trunk/psphot/src/psphotMakeResiduals.c

    r35559 r42842  
    11# include "psphotInternal.h"
     2pmConfig *psphotGetConfig ();
     3bool _psphotSaveCube (char *basename, char *extname, psArray *cube);
    24
    35# define RESIDUAL_SOFTENING 0.005
    … …  
    9395    psVector *yC = psVectorAllocEmpty (100, PS_TYPE_F32);
    9496
     97# define TESTRESID 0
     98# if (TESTRESID)
     99    psArray *sourceRawCube = psArrayAllocEmpty (1);
     100    psArray *sourceSigCube = psArrayAllocEmpty (1);
     101    psArray *sourceVarCube = psArrayAllocEmpty (1);
     102    psArray *sourceMskCube = psArrayAllocEmpty (1);
     103
     104    psImage *refImage = NULL;
     105# endif
     106
    95107    // build (DATA - MODEL) [an image] for each psf star
    96108    psArray *input = psArrayAllocEmpty (100);
    … …  
    108120        psImage *mask     = psImageCopy (NULL, source->maskView ? source->maskView : source->maskObj,  PS_TYPE_IMAGE_MASK);
    109121        psImage *variance = psImageCopy (NULL, source->variance ? source->variance : source->pixels,   PS_TYPE_F32);
     122
     123# if (TESTRESID)
     124        psImage *tmpImage = psImageCopy(NULL, image, PS_TYPE_F32);
     125        if (refImage == NULL) {
     126          refImage = tmpImage;
     127          psArrayAdd (sourceRawCube, 1, tmpImage);
     128        } else {
     129          if (tmpImage->numCols != refImage->numCols) {
     130            fprintf (stderr, "mismatched image sizes (%d : %f, %f)\n", i, model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     131            continue;
     132          }
     133          if (tmpImage->numRows != refImage->numRows) {
     134            fprintf (stderr, "mismatched image sizes (%d : %f, %f)\n", i, model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     135            continue;
     136          }
     137          psArrayAdd (sourceRawCube, 1, tmpImage);
     138          psFree (tmpImage);
     139        }
     140# endif
     141
    110142        pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal);
    111143
    … …  
    114146        psBinaryOp (image, image, "/", psScalarAlloc(Io, PS_TYPE_F32));
    115147        psBinaryOp (variance, variance, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
     148
     149# if (TESTRESID)
     150        psArrayAdd (sourceSigCube, 1, image);
     151        psArrayAdd (sourceVarCube, 1, variance);
     152        psArrayAdd (sourceMskCube, 1, mask);
     153# endif
    116154
    117155        // we interpolate the image and variance - include the mask or not?
    … …  
    142180    psImageInit (resid->mask, 0);
    143181
     182
     183# if (TESTRESID)
     184    psArray *inputSigCube = psArrayAllocEmpty (1);
     185    psArray *inputVarCube = psArrayAllocEmpty (1);
     186    psArray *inputMskCube = psArrayAllocEmpty (1);
     187
     188    // below, we are generating the residual images with interpolation
     189    // to save these in an output cube, I need to generate place-holder images here
     190    for (int i = 0; i < input->n; i++) {
     191      psImage *inputSig = psImageAlloc (resid->Ro->numCols, resid->Ro->numRows, PS_TYPE_F32);
     192      psImage *inputVar = psImageAlloc (resid->Ro->numCols, resid->Ro->numRows, PS_TYPE_F32);
     193      psImage *inputMsk = psImageAlloc (resid->Ro->numCols, resid->Ro->numRows, PS_TYPE_IMAGE_MASK);
     194
     195      psImageInit (inputSig, NAN);
     196      psImageInit (inputVar, NAN);
     197      psImageInit (inputMsk, 0);
     198
     199      psArrayAdd (inputSigCube, 1, inputSig);
     200      psArrayAdd (inputVarCube, 1, inputVar);
     201      psArrayAdd (inputMskCube, 1, inputMsk);
     202
     203      // XXX free these here, right?
     204      psFree (inputSig);
     205      psFree (inputVar);
     206      psFree (inputMsk);
     207    }
     208# endif
     209
    144210    // x(resid) = (x(image) - Xo)*xBin + xCenter
    145211
    … …  
    148214    psVector *fmasks  = psVectorAlloc (input->n, PS_TYPE_VECTOR_MASK);
    149215
    150     // statistic to use to determine baseline for clipping
    151     psStats *fluxClip     = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    152     psStats *fluxClipDef  = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    153216    // statistic to use to determine output flux
    154217    // XXX make API to convert statOption for MEAN/MEDIAN in to corresponding STDEV?
    … …  
    156219    psStats *fluxStatsDef = psStatsAlloc (statOption | PS_STAT_SAMPLE_STDEV);
    157220
     221    // use this for the IRLS fitting below (defines niter)
     222    psStats *statsIRLS = psStatsAlloc(PS_STAT_CLIPPED_MEAN); statsIRLS->clipIter = 10; // max number of iterations
     223    psPolynomial2D *polyIRLS = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1);
     224
    158225    // Use psF64 to minimize overflow problems?
    159226    psImage *A = psImageAlloc(3, 3, PS_TYPE_F64); // Least-squares matrix
    … …  
    166233
    167234            int nGoodPixel = 0;              // pixel is off the image
     235
     236            // skip pixels outside of the requested radius
     237            float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
     238            if (radius > radiusMax) {
     239              resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
     240              continue;
     241            }
    168242
    169243            // build the vector of data values for this output pixel
    … …  
    194268                    nGoodPixel ++;
    195269                }
     270
     271# if (TESTRESID)
     272                psImage *inputSig = inputSigCube->data[i]; inputSig->data.F32[oy][ox] = flux;
     273                psImage *inputVar = inputVarCube->data[i]; inputVar->data.F32[oy][ox] = flux;
     274                psImage *inputMsk = inputMskCube->data[i]; inputMsk->data.PS_TYPE_IMAGE_MASK_DATA[oy][ox] = fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i];
     275# endif
    196276            }
     277
     278            if (nGoodPixel < 0.05*input->n) {
     279              fprintf (stderr, "warning: %d, %d : residual pixel with few input pixels: %d vs %d\n", ox, oy, nGoodPixel, (int) input->n);
     280            }
    197281
    198282            // skip pixels with insufficient data
    … …  
    207291
    208292            // measure the robust median to determine a baseline reference value
    209             *fluxClip = *fluxClipDef;
     293            psStats *fluxClip = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    210294            if (!psVectorStats (fluxClip, fluxes, NULL, fmasks, fmaskVal)) {
    211295                psError(PSPHOT_ERR_CONFIG, false, "Error calculating residual stats");
     296                psFree (fluxClip);
    212297                return false;
    213298            }
    214             if (isnan(fluxClip->robustMedian)) {
     299            if (isnan(fluxClip->sampleMedian)) {
    215300                resid->Ro->data.F32[oy][ox] = 0.0;
    216301                resid->Rx->data.F32[oy][ox] = 0.0;
    217302                resid->Ry->data.F32[oy][ox] = 0.0;
    218303                resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask;
     304                psFree (fluxClip);
    219305                continue;
     306            }
     307
     308            if (fabs(fluxClip->sampleMedian) > 1.5) {
     309              fprintf (stderr, "warning: %d, %d : residual pixel has funny initial median value: %f\n", ox, oy, fluxClip->sampleMedian);
    220310            }
    221311
    … …  
    223313            int nKeep = 0;
    224314            for (int i = 0; i < fluxes->n; i++) {
    225                 float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
     315                float delta = fluxes->data.F32[i] - fluxClip->sampleMedian;
    226316                float sigma = sqrt (dfluxes->data.F32[i]);
    227317                float swing = fabs(delta) / sigma;
    … …  
    233323                if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++;
    234324            }
     325            psFree (fluxClip);
     326
     327            if (nKeep < 5) {
     328              fprintf (stderr, "warning: %d, %d : residual pixel with few good pixels: %d vs %d\n", ox, oy, nKeep, (int) input->n);
     329            }
    235330
    236331            if (SPATIAL_ORDER == 0) {
    … …  
    242337                }
    243338
    244                 float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
    245                 if (radius > radiusMax) {
    246                   resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
    247                   continue;
    248                 }
     339                resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
    249340
    250341                resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
    251                 resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
    252 
    253342                if (isnan(resid->Ro->data.F32[oy][ox])) {
    254343                    resid->Ro->data.F32[oy][ox] = 0.0;
    255                     resid->Rx->data.F32[oy][ox] = 0.0;
    256                     resid->Ry->data.F32[oy][ox] = 0.0;
    257344                    resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask;
    258345                    continue;
    … …  
    265352                assert (SPATIAL_ORDER == 1);
    266353
    267                 float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
    268                 if (radius > radiusMax) {
    269                   resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
    270                   continue;
    271                 }
    272 
    273                 psImageInit(A, 0.0);
    274                 psVectorInit(B, 0.0);
    275                 for (int i = 0; i < fluxes->n; i++) {
    276                     if (fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
    277                     B->data.F64[0] += fluxes->data.F32[i]/dfluxes->data.F32[i];
    278                     B->data.F64[1] += fluxes->data.F32[i]*xC->data.F32[i]/dfluxes->data.F32[i];
    279                     B->data.F64[2] += fluxes->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
    280 
    281                     A->data.F64[0][0] += 1.0/dfluxes->data.F32[i];
    282                     A->data.F64[1][0] += xC->data.F32[i]/dfluxes->data.F32[i];
    283                     A->data.F64[2][0] += yC->data.F32[i]/dfluxes->data.F32[i];
    284 
    285                     A->data.F64[1][1] += PS_SQR(xC->data.F32[i])/dfluxes->data.F32[i];
    286                     A->data.F64[2][2] += PS_SQR(yC->data.F32[i])/dfluxes->data.F32[i];
    287                     A->data.F64[1][2] += xC->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
    288                 }
    289 
    290                 A->data.F64[0][1] = A->data.F64[1][0];
    291                 A->data.F64[0][2] = A->data.F64[2][0];
    292                 A->data.F64[2][1] = A->data.F64[1][2];
    293 
    294                 if (!psMatrixGJSolve(A, B)) {
    295                     resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
    296                     psWarning("Singular matrix solving for (y,x) = (%d,%d)'s residuals, masking", oy, ox);
    297                     continue;
    298                 }
    299 
    300                 resid->Ro->data.F32[oy][ox] = B->data.F64[0];
    301                 resid->Rx->data.F32[oy][ox] = B->data.F64[1];
    302                 resid->Ry->data.F32[oy][ox] = B->data.F64[2];
    303 
    304                 float dRo = sqrt(A->data.F32[0][0]);
     354                // we have the following vectors to describe this pixel:
     355                // fluxes, dfluxes, fmasks, xC, yC
     356
     357                psVectorIRLSFitPolynomial2D (polyIRLS, statsIRLS, fmasks, fmaskVal, fluxes, dfluxes, xC, yC);
     358
     359                resid->Ro->data.F32[oy][ox] = polyIRLS->coeff[0][0];
     360                resid->Rx->data.F32[oy][ox] = polyIRLS->coeff[1][0];
     361                resid->Ry->data.F32[oy][ox] = polyIRLS->coeff[0][1];
     362
     363                // is the error well-defined?
     364                float dRo = polyIRLS->coeffErr[0][0];
     365
     366                if (fabs(resid->Ro->data.F32[oy][ox]) > 1.5) {
     367                  fprintf (stderr, "warning: %d, %d : residual pixel has funny fit value: %f\n", ox, oy, resid->Ro->data.F32[oy][ox]);
     368                }
    305369
    306370                if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) {
     371                  // fprintf (stderr, "mask : %d, %d : %f +/- %f (%f for %d) : %d\n", ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]);
    307372                  resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
    308373                  resid->Ro->data.F32[oy][ox] = 0.0;
    309374                  resid->Rx->data.F32[oy][ox] = 0.0;
    310375                  resid->Ry->data.F32[oy][ox] = 0.0;
    311                 }
     376                } else {
     377                  // fprintf (stderr, "keep : %d, %d : %f +/- %f (%f for %d) : %d\n", ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]);
     378                }
    312379            }
    313380        }
    314381    }
     382
     383# if (TESTRESID)
     384
     385    pmConfig *config = psphotGetConfig();
     386    char *output = psMetadataLookupStr (&status, config->arguments, "OUTPUT");
     387   
     388    // XXX // save the star coordinate vectors
     389    // XXX
     390    // XXX psString filename = NULL;
     391    // XXX psStringAppend (&filename, "%s.instar.dat", output);
     392
     393    char filename[1024];
     394    sprintf (filename, "%s.instar.dat", output);
     395   
     396    FILE *fout = fopen (filename, "w");
     397    for (int i = 0; i < xC->n; i++) {
     398      fprintf (fout, "%d %f %f\n", i, xC->data.F32[i], yC->data.F32[i]);
     399    }
     400    fclose (fout);
     401
     402    _psphotSaveCube (output, "raw.raw.fits", sourceRawCube);
     403    _psphotSaveCube (output, "sig.raw.fits", sourceSigCube);
     404    _psphotSaveCube (output, "var.raw.fits", sourceVarCube);
     405    _psphotSaveCube (output, "msk.raw.fits", sourceMskCube);
     406
     407    _psphotSaveCube (output, "sig.fits", inputSigCube);
     408    _psphotSaveCube (output, "var.fits", inputVarCube);
     409    _psphotSaveCube (output, "msk.fits", inputMskCube);
     410
     411    psFree (sourceRawCube);
     412    psFree (sourceSigCube);
     413    psFree (sourceVarCube);
     414    psFree (sourceMskCube);
     415
     416    psFree (inputSigCube);
     417    psFree (inputVarCube);
     418    psFree (inputMskCube);
     419
     420# endif
    315421
    316422    psFree (A);
    … …  
    329435    psFree (fluxStats);
    330436    psFree (fluxStatsDef);
    331     psFree (fluxClip);
    332     psFree (fluxClipDef);
    333437
    334438    if (resid != NULL && psTraceGetLevel("psphot") > 5) {
    … …  
    340444    }
    341445
     446    psFree (statsIRLS);
     447    psFree (polyIRLS);
     448
    342449    psf->residuals = resid;
    343450    return (resid != NULL) ? true : false;
    344451}
     452
     453bool _psphotSaveCube (char *basename, char *extname, psArray *cube) {
     454
     455    psString filename = NULL;
     456    psStringAppend (&filename, "%s.%s", basename, extname);
     457
     458    psFits *fits = psFitsOpen (filename, "w");
     459    if (!psFitsWriteImageCube (fits, NULL, cube, NULL)) {
     460      fprintf (stderr, "failed to write the cube %s\n", filename);
     461    }
     462    psFitsClose (fits);
     463
     464    psFree (filename);
     465
     466    return true;
     467}
Note: See TracChangeset for help on using the changeset viewer.