IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 12948


Ignore:
Timestamp:
Apr 20, 2007, 2:18:15 PM (19 years ago)
Author:
eugene
Message:

added spatial variations

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_02_branch/psphot/src/psphotMakeResiduals.c

    r12933 r12948  
    11# include "psphotInternal.h"
     2# define ZERO_ORDER 0
    23
    34bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf) {
    … …  
    910    psTimerStart ("residuals");
    1011
     12    if (!psMetadataLookupBool(&status, recipe, "PSF.RESIDUALS")) return true;
     13
    1114    int xBin = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.XBIN");
    1215    PS_ASSERT (status, false);
    … …  
    1417    int yBin = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.YBIN");
    1518    PS_ASSERT (status, false);
     19
     20    float nSigma = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.NSIGMA");
     21    PS_ASSERT (status, false);
     22
     23    char *modeString = psMetadataLookupStr(&status, recipe, "PSF.RESIDUALS.INTERPOLATION");
     24    PS_ASSERT (status, false);
     25
     26    psImageInterpolateMode mode = psImageInterpolateModeFromString (modeString);
     27    if (mode == PS_INTERPOLATE_NONE) {
     28        psError(PSPHOT_ERR_CONFIG, false, "invalid interpolation in psphot.config");
     29        return false;
     30    }
     31
     32    char *statString = psMetadataLookupStr(&status, recipe, "PSF.RESIDUALS.STATISTIC");
     33    PS_ASSERT (status, false);
     34
     35    psStatsOptions statOption = psStatsOptionFromString (statString);
     36    if (!statOption) {
     37        psError(PSPHOT_ERR_CONFIG, false, "invalid residual statistic in psphot.config");
     38        return false;
     39    }
    1640
    1741    // user parameters:
    … …  
    4064    psVector *yC = psVectorAllocEmpty (100, PS_TYPE_F32);
    4165
     66    // build (DATA - MODEL) [an image] for each psf star
    4267    psArray *input = psArrayAllocEmpty (100);
    4368    for (int i = 0; i < sources->n; i++) {
    … …  
    5580        psImage *mask   = psImageCopy (NULL, source->mask,   PS_TYPE_U8);
    5681        pmModelSub (image, mask, model, false, false);
    57         psFree (mask);
    5882       
    5983        // re-normalize image and weight
    … …  
    6286        psBinaryOp (weight, weight, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
    6387
    64         psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(
    65             PS_INTERPOLATE_BILINEAR,
    66             image, weight, NULL, 0, 0.0, 0.0, 1, 0, 0.0);
     88        // we will interpolate the image and weight - include the mask or not?
     89        // XXX consider better values for the mask bits
     90        psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(mode, image, weight, NULL, 0xff, 0.0, 0.0, 1, 2, 0.0);
    6791        psArrayAdd (input,  100, interp);
    6892
    … …  
    77101
    78102        // free up the excess references
     103        psFree (mask);
    79104        psFree (image);
    80105        psFree (weight);
    … …  
    89114    psVector *fmasks  = psVectorAlloc (input->n, PS_TYPE_U8);
    90115
     116    // statistic to use to determine baseline for clipping
    91117    psStats *fluxClip     = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    92118    psStats *fluxClipDef  = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    93     psStats *fluxStats    = psStatsAlloc (PS_STAT_SAMPLE_MEAN   | PS_STAT_SAMPLE_STDEV);
    94     psStats *fluxStatsDef = psStatsAlloc (PS_STAT_SAMPLE_MEAN   | PS_STAT_SAMPLE_STDEV);
    95 
    96     // build the residual image pixel-by-pixel
    97     for (int oy = 0; oy < resid->image->numRows; oy++) {
     119    // statistic to use to determine output flux
     120    // XXX make API to convert statOption for MEAN/MEDIAN in to corresponding STDEV?
     121    psStats *fluxStats    = psStatsAlloc (statOption | PS_STAT_SAMPLE_STDEV);
     122    psStats *fluxStatsDef = psStatsAlloc (statOption | PS_STAT_SAMPLE_STDEV);
     123
     124// this section builds just the 0th order term Ro
     125# if (ZERO_ORDER)
     126    // build Ro = DATA - MODEL (rebinned image) pixel-by-pixel
     127    for (int oy = 0; oy < resid->Ro->numRows; oy++) {
    98128        fprintf (stderr, ".");
    99         for (int ox = 0; ox < resid->image->numCols; ox++) {
     129        for (int ox = 0; ox < resid->Ro->numCols; ox++) {
    100130           
     131            // build the vector of data values for this output pixel
    101132            for (int i = 0; i < input->n; i++) {
    102133
    … …  
    104135               
    105136                // fractional image position
    106                 float ix = (ox - resid->xCenter) / (float) xBin + xC->data.F32[i] - interp->image->col0;
    107                 float iy = (oy - resid->yCenter) / (float) yBin + yC->data.F32[i] - interp->image->row0;               
     137                float ix = (ox + 0.5 - resid->xCenter) / (float) xBin + xC->data.F32[i] - interp->image->col0;
     138                float iy = (oy + 0.5 - resid->yCenter) / (float) yBin + yC->data.F32[i] - interp->image->row0;
    108139
    109140                mflux = 0;
    … …  
    115146            }
    116147
     148            // measure the robust median to determine a baseline reference value
    117149            *fluxClip = *fluxClipDef;
    118150            psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff);
    119             psErrorClear();
    120151
    121152            // mark input pixels which are more than N sigma from the median
    … …  
    123154                float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
    124155                float sigma = sqrt (dfluxes->data.F32[i]);
    125                
    126156                float swing = fabs(delta) / sigma;
    127157
    128158                // make this a user option
    129                 if (swing > 3.0) {
     159                if (swing > nSigma) {
    130160                    fmasks->data.U8[i] = 1;
    131161                }
    132162            }               
    133163
     164            // measure the desired statistic on the unclipped pixels
    134165            *fluxStats = *fluxStatsDef;
    135166            psVectorStats (fluxStats, fluxes, NULL, fmasks, 0xff);
     167
     168            resid->Ro->data.F32[oy][ox] = psStatsGetValue (fluxStats, statOption);
     169            resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
     170
     171            // clear (ignore) any outstanding errors
    136172            psErrorClear();
    137 
    138             resid->image->data.F32[oy][ox] = fluxStats->sampleMean;
    139             resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
    140173        }
    141174    }
    142 
    143     psphotSaveImage (NULL, resid->image, "resid.im.fits");
    144     psphotSaveImage (NULL, resid->weight, "resid.wt.fits");
    145     psphotSaveImage (NULL, resid->mask, "resid.mk.fits");
    146 
    147     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
     175    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generated 0-th order residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
     176
     177# else
     178
     179    psImage *A = psImageAlloc(3, 3, PS_TYPE_F64); // Least-squares matrix
     180    psVector *B = psVectorAlloc(3, PS_TYPE_F64); // Least-squares vector
     181
     182    // build (x,y)*(DATA - MODEL - Ro) pixel-by-pixel
     183    for (int oy = 0; oy < resid->Ro->numRows; oy++) {
     184        fprintf (stderr, ".");
     185        for (int ox = 0; ox < resid->Ro->numCols; ox++) {
     186           
     187            // build the vector of data values for this output pixel
     188            // XXX this is identical to the pass above: we could cache the results for speed
     189            for (int i = 0; i < input->n; i++) {
     190
     191                psImageInterpolateOptions *interp = input->data[i];
     192               
     193                // fractional image position
     194                float ix = (ox + 0.5 - resid->xCenter) / (float) xBin + xC->data.F32[i] - interp->image->col0;
     195                float iy = (oy + 0.5 - resid->yCenter) / (float) yBin + yC->data.F32[i] - interp->image->row0;
     196
     197                mflux = 0;
     198                psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp);
     199                fluxes->data.F32[i] = flux;
     200                dfluxes->data.F32[i] = dflux;
     201                fmasks->data.U8[i] = mflux;
     202                // fprintf (stderr, "%f %f : %f %f (%d)\n", ix, iy, flux, dflux, fmasks->data.U8[i]);
     203            }
     204
     205            // measure the robust median to determine a baseline reference value
     206            *fluxClip = *fluxClipDef;
     207            psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff);
     208            psErrorClear();             // clear (ignore) any outstanding errors
     209
     210            // mark input pixels which are more than N sigma from the median
     211            for (int i = 0; i < fluxes->n; i++) {
     212                float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
     213                float sigma = sqrt (dfluxes->data.F32[i]);
     214                float swing = fabs(delta) / sigma;
     215
     216                // make this a user option
     217                if (swing > nSigma) {
     218                    fmasks->data.U8[i] = 1;
     219                }
     220            }               
     221
     222            psImageInit(A, 0.0);
     223            psVectorInit(B, 0.0);
     224            for (int i = 0; i < fluxes->n; i++) {
     225                if (fmasks->data.U8[i]) continue;
     226                B->data.F64[0] += fluxes->data.F32[i]/dfluxes->data.F32[i];
     227                B->data.F64[1] += fluxes->data.F32[i]*xC->data.F32[i]/dfluxes->data.F32[i];
     228                B->data.F64[2] += fluxes->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
     229
     230                A->data.F64[0][0] += 1.0/dfluxes->data.F32[i];
     231                A->data.F64[1][0] += xC->data.F32[i]/dfluxes->data.F32[i];
     232                A->data.F64[2][0] += yC->data.F32[i]/dfluxes->data.F32[i];
     233
     234                A->data.F64[1][1] += PS_SQR(xC->data.F32[i])/dfluxes->data.F32[i];
     235                A->data.F64[2][2] += PS_SQR(yC->data.F32[i])/dfluxes->data.F32[i];
     236                A->data.F64[1][2] += xC->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
     237            }
     238
     239            A->data.F64[0][1] = A->data.F64[1][0];
     240            A->data.F64[0][2] = A->data.F64[2][0];
     241            A->data.F64[2][1] = A->data.F64[1][2];
     242            psMatrixGJSolve(A, B);
     243
     244            resid->Ro->data.F32[oy][ox] = B->data.F64[0];
     245            resid->Rx->data.F32[oy][ox] = B->data.F64[1];
     246            resid->Ry->data.F32[oy][ox] = B->data.F64[2];
     247        }
     248    }
     249
     250    psFree (A);
     251    psFree (B);
     252
     253# endif
    148254
    149255    psFree (xC);
    … …  
    160266    psFree (fluxClipDef);
    161267
     268    psphotSaveImage (NULL, resid->Ro,     "resid.ro.fits");
     269    psphotSaveImage (NULL, resid->Rx,     "resid.rx.fits");
     270    psphotSaveImage (NULL, resid->Ry,     "resid.ry.fits");
     271    psphotSaveImage (NULL, resid->weight, "resid.wt.fits");
     272    psphotSaveImage (NULL, resid->mask,   "resid.mk.fits");
     273
     274    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
     275
    162276    psf->residuals = resid;
    163277    return true;
Note: See TracChangeset for help on using the changeset viewer.