IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 13806


Ignore:
Timestamp:
Jun 13, 2007, 2:32:20 PM (19 years ago)
Author:
Paul Price
Message:

Clear an error in the stats.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotMakeResiduals.c

    r13570 r13806  
    1414    PS_ASSERT (status, false);
    1515    if (SPATIAL_ORDER != 0 && SPATIAL_ORDER != 1) {
    16         psError(PSPHOT_ERR_CONFIG, true, "PSF.RESIDUALS.SPATIAL_ORDER must be 0 or 1 (not %d)",
    17                 SPATIAL_ORDER);
    18         return false;
     16        psError(PSPHOT_ERR_CONFIG, true, "PSF.RESIDUALS.SPATIAL_ORDER must be 0 or 1 (not %d)",
     17                SPATIAL_ORDER);
     18        return false;
    1919    }
    2020
    … …  
    3535    psImageInterpolateMode mode = psImageInterpolateModeFromString (modeString);
    3636    if (mode == PS_INTERPOLATE_NONE) {
    37         psError(PSPHOT_ERR_CONFIG, false, "invalid interpolation in psphot.config");
    38         return false;
     37        psError(PSPHOT_ERR_CONFIG, false, "invalid interpolation in psphot.config");
     38        return false;
    3939    }
    4040
    … …  
    4444    psStatsOptions statOption = psStatsOptionFromString (statString);
    4545    if (!statOption) {
    46         psError(PSPHOT_ERR_CONFIG, false, "invalid residual statistic in psphot.config");
    47         return false;
     46        psError(PSPHOT_ERR_CONFIG, false, "invalid residual statistic in psphot.config");
     47        return false;
    4848    }
    4949
    … …  
    5656    // - construct a residual image, renormalized
    5757    // - construct a renormalized weight image
    58     // - construct a new mask image 
     58    // - construct a new mask image
    5959
    6060    // construct the output residual table (Nx*DX,Ny*DY)
    … …  
    6666    // - set output pixel, weight, and mask
    6767
    68     const int badMask = 1;              // mask bits
    69     const int poorMask = 2;             //       from psImageInterpolate
    70     const int clippedMask = 4;          // mask bit set for clipped values
    71 
    72     bool offImage = false;              // pixel is off the image
     68    const int badMask = 1;              // mask bits
     69    const int poorMask = 2;             //       from psImageInterpolate
     70    const int clippedMask = 4;          // mask bit set for clipped values
     71
     72    bool offImage = false;              // pixel is off the image
    7373
    7474    // determine the maximum image size from the input sources
    … …  
    8787        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
    8888
    89         // which model to use?
    90         pmModel *model = pmSourceGetModel (&isPSF, source);
    91         if (model == NULL) continue;  // model must be defined
     89        // which model to use?
     90        pmModel *model = pmSourceGetModel (&isPSF, source);
     91        if (model == NULL) continue;  // model must be defined
    9292
    9393        psImage *image  = psImageCopy (NULL, source->pixels,   PS_TYPE_F32);
    … …  
    9595        psImage *mask   = psImageCopy (NULL, source->maskView, PS_TYPE_U8);
    9696        pmModelSub (image, mask, model, PM_MODEL_OP_FUNC);
    97        
    98         // re-normalize image and weight
    99         float Io = model->params->data.F32[PM_PAR_I0];
    100         psBinaryOp (image, image, "/", psScalarAlloc(Io, PS_TYPE_F32));
    101         psBinaryOp (weight, weight, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
    102 
    103         // we will interpolate the image and weight - include the mask or not?
    104         // XXX consider better values for the mask bits
    105         psImageInterpolateOptions *interp =
    106             psImageInterpolateOptionsAlloc(mode, image, weight, NULL, 0xff, 0.0, 0.0, badMask, poorMask, 0.0);
    107         psArrayAdd (input,  100, interp);
    108 
    109         // save the X,Y position for future reference
    110         xC->data.F32[xC->n] = model->params->data.F32[PM_PAR_XPOS];
    111         yC->data.F32[yC->n] = model->params->data.F32[PM_PAR_YPOS];
    112         psVectorExtend (xC, 100, 1);
    113         psVectorExtend (yC, 100, 1);
    114 
    115         xSize = PS_MAX (xSize, image->numCols);
    116         ySize = PS_MAX (ySize, image->numRows);
    117 
    118         // free up the excess references
    119         psFree (mask);
    120         psFree (image);
    121         psFree (weight);
    122         psFree (interp);
     97
     98        // re-normalize image and weight
     99        float Io = model->params->data.F32[PM_PAR_I0];
     100        psBinaryOp (image, image, "/", psScalarAlloc(Io, PS_TYPE_F32));
     101        psBinaryOp (weight, weight, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
     102
     103        // we will interpolate the image and weight - include the mask or not?
     104        // XXX consider better values for the mask bits
     105        psImageInterpolateOptions *interp =
     106            psImageInterpolateOptionsAlloc(mode, image, weight, NULL, 0xff, 0.0, 0.0, badMask, poorMask, 0.0);
     107        psArrayAdd (input,  100, interp);
     108
     109        // save the X,Y position for future reference
     110        xC->data.F32[xC->n] = model->params->data.F32[PM_PAR_XPOS];
     111        yC->data.F32[yC->n] = model->params->data.F32[PM_PAR_YPOS];
     112        psVectorExtend (xC, 100, 1);
     113        psVectorExtend (yC, 100, 1);
     114
     115        xSize = PS_MAX (xSize, image->numCols);
     116        ySize = PS_MAX (ySize, image->numRows);
     117
     118        // free up the excess references
     119        psFree (mask);
     120        psFree (image);
     121        psFree (weight);
     122        psFree (interp);
    123123    }
    124124    pmResiduals *resid = pmResidualsAlloc (xSize, ySize, xBin, yBin);
    125125
    126126    // x(resid) = (x(image) - Xo)*xBin + xCenter
    127    
     127
    128128    psVector *fluxes  = psVectorAlloc (input->n, PS_TYPE_F32);
    129129    psVector *dfluxes = psVectorAlloc (input->n, PS_TYPE_F32);
    … …  
    144144    // (If SPATIAL_ORDER == 0, just solve MODEL = R)
    145145    for (int oy = 0; resid != NULL && oy < resid->Ro->numRows; oy++) {
    146         for (int ox = 0; ox < resid->Ro->numCols; ox++) {
    147            
    148             // build the vector of data values for this output pixel
    149             for (int i = 0; i < input->n; i++) {
    150 
    151                 psImageInterpolateOptions *interp = input->data[i];
    152                
    153                 // fractional image position
    154                 float ix = (ox + 0.5 - resid->xCenter) / (float) xBin + xC->data.F32[i] - interp->image->col0;
    155                 float iy = (oy + 0.5 - resid->yCenter) / (float) yBin + yC->data.F32[i] - interp->image->row0;
    156 
    157                 mflux = 0;
    158                 offImage = false;
    159                 if (psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp) == PS_INTERPOLATE_STATUS_OFF) {
    160                     // This pixel is off the image
    161                     offImage = true;
    162                 }
    163                 fluxes->data.F32[i] = flux;
    164                 dfluxes->data.F32[i] = dflux;
    165                 fmasks->data.U8[i] = mflux;
    166                 // fprintf (stderr, "%f %f : %f %f (%d)\n", ix, iy, flux, dflux, fmasks->data.U8[i]);
    167             }
    168 
    169             // measure the robust median to determine a baseline reference value
    170             *fluxClip = *fluxClipDef;
    171             psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff);
    172             psErrorClear();             // clear (ignore) any outstanding errors
    173 
    174             // mark input pixels which are more than N sigma from the median
    175             for (int i = 0; i < fluxes->n; i++) {
    176                 float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
    177                 float sigma = sqrt (dfluxes->data.F32[i]);
    178                 float swing = fabs(delta) / sigma;
    179 
    180                 // make this a user option
    181                 if (swing > nSigma) {
    182                     fmasks->data.U8[i] = clippedMask;
    183                 }
    184             }               
    185 
    186             if (offImage || SPATIAL_ORDER == 0) {
    187                 // measure the desired statistic on the unclipped pixels
    188                 *fluxStats = *fluxStatsDef;
    189                 psVectorStats (fluxStats, fluxes, NULL, fmasks, 0xff);
    190                
    191                 resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
    192                 resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
    193                 //resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
    194             } else {
    195                 assert (SPATIAL_ORDER == 1);
    196                 psImageInit(A, 0.0);
    197                 psVectorInit(B, 0.0);
    198                 for (int i = 0; i < fluxes->n; i++) {
    199                     if (fmasks->data.U8[i]) continue;
    200                     B->data.F64[0] += fluxes->data.F32[i]/dfluxes->data.F32[i];
    201                     B->data.F64[1] += fluxes->data.F32[i]*xC->data.F32[i]/dfluxes->data.F32[i];
    202                     B->data.F64[2] += fluxes->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
    203 
    204                     A->data.F64[0][0] += 1.0/dfluxes->data.F32[i];
    205                     A->data.F64[1][0] += xC->data.F32[i]/dfluxes->data.F32[i];
    206                     A->data.F64[2][0] += yC->data.F32[i]/dfluxes->data.F32[i];
    207 
    208                     A->data.F64[1][1] += PS_SQR(xC->data.F32[i])/dfluxes->data.F32[i];
    209                     A->data.F64[2][2] += PS_SQR(yC->data.F32[i])/dfluxes->data.F32[i];
    210                     A->data.F64[1][2] += xC->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
    211                 }
    212 
    213                 A->data.F64[0][1] = A->data.F64[1][0];
    214                 A->data.F64[0][2] = A->data.F64[2][0];
    215                 A->data.F64[2][1] = A->data.F64[1][2];
    216 
    217                 if (!psMatrixGJSolve(A, B)) {
    218                     psError(PSPHOT_ERR_PSF, false, "Singular matrix solving for (y,x) = (%d,%d)'s residuals",
    219                             oy, ox);
    220                     psFree(resid); resid = NULL;
    221                     break;                 
    222                 }
    223 
    224                 resid->Ro->data.F32[oy][ox] = B->data.F64[0];
    225                 resid->Rx->data.F32[oy][ox] = B->data.F64[1];
    226                 resid->Ry->data.F32[oy][ox] = B->data.F64[2];
    227                 //resid->weight->data.F32[oy][ox] = XXX;
    228             }
    229         }
     146        for (int ox = 0; ox < resid->Ro->numCols; ox++) {
     147
     148            // build the vector of data values for this output pixel
     149            for (int i = 0; i < input->n; i++) {
     150
     151                psImageInterpolateOptions *interp = input->data[i];
     152
     153                // fractional image position
     154                float ix = (ox + 0.5 - resid->xCenter) / (float) xBin + xC->data.F32[i] - interp->image->col0;
     155                float iy = (oy + 0.5 - resid->yCenter) / (float) yBin + yC->data.F32[i] - interp->image->row0;
     156
     157                mflux = 0;
     158                offImage = false;
     159                if (psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp) == PS_INTERPOLATE_STATUS_OFF) {
     160                    // This pixel is off the image
     161                    offImage = true;
     162                }
     163                fluxes->data.F32[i] = flux;
     164                dfluxes->data.F32[i] = dflux;
     165                fmasks->data.U8[i] = mflux;
     166                // fprintf (stderr, "%f %f : %f %f (%d)\n", ix, iy, flux, dflux, fmasks->data.U8[i]);
     167            }
     168
     169            // measure the robust median to determine a baseline reference value
     170            *fluxClip = *fluxClipDef;
     171            psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff);
     172            psErrorClear();             // clear (ignore) any outstanding errors
     173
     174            // mark input pixels which are more than N sigma from the median
     175            for (int i = 0; i < fluxes->n; i++) {
     176                float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
     177                float sigma = sqrt (dfluxes->data.F32[i]);
     178                float swing = fabs(delta) / sigma;
     179
     180                // make this a user option
     181                if (swing > nSigma) {
     182                    fmasks->data.U8[i] = clippedMask;
     183                }
     184            }
     185
     186            if (offImage || SPATIAL_ORDER == 0) {
     187                // measure the desired statistic on the unclipped pixels
     188                *fluxStats = *fluxStatsDef;
     189                psVectorStats (fluxStats, fluxes, NULL, fmasks, 0xff);
     190                psErrorClear();         // clear (ignore) any outstanding errors
     191
     192                resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
     193                resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
     194                //resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
     195            } else {
     196                assert (SPATIAL_ORDER == 1);
     197                psImageInit(A, 0.0);
     198                psVectorInit(B, 0.0);
     199                for (int i = 0; i < fluxes->n; i++) {
     200                    if (fmasks->data.U8[i]) continue;
     201                    B->data.F64[0] += fluxes->data.F32[i]/dfluxes->data.F32[i];
     202                    B->data.F64[1] += fluxes->data.F32[i]*xC->data.F32[i]/dfluxes->data.F32[i];
     203                    B->data.F64[2] += fluxes->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
     204
     205                    A->data.F64[0][0] += 1.0/dfluxes->data.F32[i];
     206                    A->data.F64[1][0] += xC->data.F32[i]/dfluxes->data.F32[i];
     207                    A->data.F64[2][0] += yC->data.F32[i]/dfluxes->data.F32[i];
     208
     209                    A->data.F64[1][1] += PS_SQR(xC->data.F32[i])/dfluxes->data.F32[i];
     210                    A->data.F64[2][2] += PS_SQR(yC->data.F32[i])/dfluxes->data.F32[i];
     211                    A->data.F64[1][2] += xC->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
     212                }
     213
     214                A->data.F64[0][1] = A->data.F64[1][0];
     215                A->data.F64[0][2] = A->data.F64[2][0];
     216                A->data.F64[2][1] = A->data.F64[1][2];
     217
     218                if (!psMatrixGJSolve(A, B)) {
     219                    psError(PSPHOT_ERR_PSF, false, "Singular matrix solving for (y,x) = (%d,%d)'s residuals",
     220                            oy, ox);
     221                    psFree(resid); resid = NULL;
     222                    break;
     223                }
     224
     225                resid->Ro->data.F32[oy][ox] = B->data.F64[0];
     226                resid->Rx->data.F32[oy][ox] = B->data.F64[1];
     227                resid->Ry->data.F32[oy][ox] = B->data.F64[2];
     228                //resid->weight->data.F32[oy][ox] = XXX;
     229            }
     230        }
    230231    }
    231232
Note: See TracChangeset for help on using the changeset viewer.