IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 13407


Ignore:
Timestamp:
May 17, 2007, 3:35:34 AM (19 years ago)
Author:
rhl
Message:

Allow linear or constant order; check for singular matrix

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotMakeResiduals.c

    r13196 r13407  
    11# include "psphotInternal.h"
    2 # define ZERO_ORDER 0
    32
    43bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf) {
    … …  
    1211    if (!psMetadataLookupBool(&status, recipe, "PSF.RESIDUALS")) return true;
    1312
     13    int SPATIAL_ORDER = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.SPATIAL_ORDER");
     14#if 0                                   // True when PSF.RESIDUALS.SPATIAL_ORDER's in default recipe
     15    PS_ASSERT (status, false);
     16#else
     17    if (!status) {
     18        SPATIAL_ORDER = 1;              // i.e. linear
     19    }
     20#endif
     21    if (SPATIAL_ORDER != 0 && SPATIAL_ORDER != 1) {
     22        psError(PSPHOT_ERR_CONFIG, true, "PSF.RESIDUALS.SPATIAL_ORDER must be 0 or 1 (not %d)",
     23                SPATIAL_ORDER);
     24        return false;
     25    }
     26
    1427    int xBin = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.XBIN");
    1528    PS_ASSERT (status, false);
     29    psErrorClear(); PS_ASSERT (xBin != 0, false);
    1630
    1731    int yBin = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.YBIN");
    1832    PS_ASSERT (status, false);
     33    psErrorClear(); PS_ASSERT (yBin != 0, false);
    1934
    2035    float nSigma = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.NSIGMA");
    … …  
    122137    psStats *fluxStatsDef = psStatsAlloc (statOption | PS_STAT_SAMPLE_STDEV);
    123138
    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++) {
     139    psImage *A = psImageAlloc(3, 3, PS_TYPE_F64); // Least-squares matrix
     140    psVector *B = psVectorAlloc(3, PS_TYPE_F64); // Least-squares vector
     141
     142    // Solve MODEL = R + x R_x + y R_y in pixel-by-pixel a least-squares sense
     143    // (If SPATIAL_ORDER == 0, just solve MODEL = R)
     144    for (int oy = 0; resid != NULL && oy < resid->Ro->numRows; oy++) {
    128145        for (int ox = 0; ox < resid->Ro->numCols; ox++) {
    129146           
    … …  
    148165            *fluxClip = *fluxClipDef;
    149166            psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff);
     167            psErrorClear();             // clear (ignore) any outstanding errors
    150168
    151169            // mark input pixels which are more than N sigma from the median
    … …  
    161179            }               
    162180
    163             // measure the desired statistic on the unclipped pixels
    164             *fluxStats = *fluxStatsDef;
    165             psVectorStats (fluxStats, fluxes, NULL, fmasks, 0xff);
    166 
    167             resid->Ro->data.F32[oy][ox] = psStatsGetValue (fluxStats, statOption);
    168             resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
    169 
    170             // clear (ignore) any outstanding errors
    171             psErrorClear();
    172         }
    173     }
    174     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generated 0-th order residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
    175 
    176 # else
    177 
    178     psImage *A = psImageAlloc(3, 3, PS_TYPE_F64); // Least-squares matrix
    179     psVector *B = psVectorAlloc(3, PS_TYPE_F64); // Least-squares vector
    180 
    181     // build (x,y)*(DATA - MODEL - Ro) pixel-by-pixel
    182     for (int oy = 0; oy < resid->Ro->numRows; oy++) {
    183         for (int ox = 0; ox < resid->Ro->numCols; ox++) {
    184            
    185             // build the vector of data values for this output pixel
    186             // XXX this is identical to the pass above: we could cache the results for speed
    187             for (int i = 0; i < input->n; i++) {
    188 
    189                 psImageInterpolateOptions *interp = input->data[i];
     181            if (SPATIAL_ORDER == 0) {
     182                // measure the desired statistic on the unclipped pixels
     183                *fluxStats = *fluxStatsDef;
     184                psVectorStats (fluxStats, fluxes, NULL, fmasks, 0xff);
    190185               
    191                 // fractional image position
    192                 float ix = (ox + 0.5 - resid->xCenter) / (float) xBin + xC->data.F32[i] - interp->image->col0;
    193                 float iy = (oy + 0.5 - resid->yCenter) / (float) yBin + yC->data.F32[i] - interp->image->row0;
    194 
    195                 mflux = 0;
    196                 psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp);
    197                 fluxes->data.F32[i] = flux;
    198                 dfluxes->data.F32[i] = dflux;
    199                 fmasks->data.U8[i] = mflux;
    200                 // fprintf (stderr, "%f %f : %f %f (%d)\n", ix, iy, flux, dflux, fmasks->data.U8[i]);
     186                resid->Ro->data.F32[oy][ox] = psStatsGetValue (fluxStats, statOption);
     187                //resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
     188            } else {
     189                assert (SPATIAL_ORDER == 1);
     190                psImageInit(A, 0.0);
     191                psVectorInit(B, 0.0);
     192                for (int i = 0; i < fluxes->n; i++) {
     193                    if (fmasks->data.U8[i]) continue;
     194                    B->data.F64[0] += fluxes->data.F32[i]/dfluxes->data.F32[i];
     195                    B->data.F64[1] += fluxes->data.F32[i]*xC->data.F32[i]/dfluxes->data.F32[i];
     196                    B->data.F64[2] += fluxes->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
     197
     198                    A->data.F64[0][0] += 1.0/dfluxes->data.F32[i];
     199                    A->data.F64[1][0] += xC->data.F32[i]/dfluxes->data.F32[i];
     200                    A->data.F64[2][0] += yC->data.F32[i]/dfluxes->data.F32[i];
     201
     202                    A->data.F64[1][1] += PS_SQR(xC->data.F32[i])/dfluxes->data.F32[i];
     203                    A->data.F64[2][2] += PS_SQR(yC->data.F32[i])/dfluxes->data.F32[i];
     204                    A->data.F64[1][2] += xC->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
     205                }
     206
     207                A->data.F64[0][1] = A->data.F64[1][0];
     208                A->data.F64[0][2] = A->data.F64[2][0];
     209                A->data.F64[2][1] = A->data.F64[1][2];
     210
     211                if (!psMatrixGJSolve(A, B)) {
     212                    psError(PSPHOT_ERR_PSF, false, "Singular matrix solving for (r,x) = (%d,%d)'s residuals",
     213                            oy, ox);
     214                    psFree(resid); resid = NULL;
     215                    break;                 
     216                }
     217
     218                resid->Ro->data.F32[oy][ox] = B->data.F64[0];
     219                resid->Rx->data.F32[oy][ox] = B->data.F64[1];
     220                resid->Ry->data.F32[oy][ox] = B->data.F64[2];
    201221            }
    202 
    203             // measure the robust median to determine a baseline reference value
    204             *fluxClip = *fluxClipDef;
    205             psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff);
    206             psErrorClear();             // clear (ignore) any outstanding errors
    207 
    208             // mark input pixels which are more than N sigma from the median
    209             for (int i = 0; i < fluxes->n; i++) {
    210                 float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
    211                 float sigma = sqrt (dfluxes->data.F32[i]);
    212                 float swing = fabs(delta) / sigma;
    213 
    214                 // make this a user option
    215                 if (swing > nSigma) {
    216                     fmasks->data.U8[i] = 1;
    217                 }
    218             }               
    219 
    220             psImageInit(A, 0.0);
    221             psVectorInit(B, 0.0);
    222             for (int i = 0; i < fluxes->n; i++) {
    223                 if (fmasks->data.U8[i]) continue;
    224                 B->data.F64[0] += fluxes->data.F32[i]/dfluxes->data.F32[i];
    225                 B->data.F64[1] += fluxes->data.F32[i]*xC->data.F32[i]/dfluxes->data.F32[i];
    226                 B->data.F64[2] += fluxes->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
    227 
    228                 A->data.F64[0][0] += 1.0/dfluxes->data.F32[i];
    229                 A->data.F64[1][0] += xC->data.F32[i]/dfluxes->data.F32[i];
    230                 A->data.F64[2][0] += yC->data.F32[i]/dfluxes->data.F32[i];
    231 
    232                 A->data.F64[1][1] += PS_SQR(xC->data.F32[i])/dfluxes->data.F32[i];
    233                 A->data.F64[2][2] += PS_SQR(yC->data.F32[i])/dfluxes->data.F32[i];
    234                 A->data.F64[1][2] += xC->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
    235             }
    236 
    237             A->data.F64[0][1] = A->data.F64[1][0];
    238             A->data.F64[0][2] = A->data.F64[2][0];
    239             A->data.F64[2][1] = A->data.F64[1][2];
    240             psMatrixGJSolve(A, B);
    241 
    242             resid->Ro->data.F32[oy][ox] = B->data.F64[0];
    243             resid->Rx->data.F32[oy][ox] = B->data.F64[1];
    244             resid->Ry->data.F32[oy][ox] = B->data.F64[2];
    245222        }
    246223    }
    … …  
    248225    psFree (A);
    249226    psFree (B);
    250 
    251 # endif
    252227
    253228    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
    … …  
    266241    psFree (fluxClipDef);
    267242
    268     if (psTraceGetLevel("psphot") > 5) {
     243    if (resid != NULL && psTraceGetLevel("psphot") > 5) {
    269244      psphotSaveImage (NULL, resid->Ro,     "resid.ro.fits");
    270245      psphotSaveImage (NULL, resid->Rx,     "resid.rx.fits");
    … …  
    275250
    276251    psf->residuals = resid;
    277     return true;
     252    return (resid != NULL) ? true : false;
    278253}
Note: See TracChangeset for help on using the changeset viewer.