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/psphotFitSourcesLinear.c

    r13375 r13900  
    55// model with selected pixels, and the fit radius must be defined
    66
    7 // given the set of sources, each of which points to the pixels in the 
    8 // science image, we construct a set of simulated sources with their own pixels. 
    9 // these are used to determine the simultaneous linear fit of fluxes. 
     7// given the set of sources, each of which points to the pixels in the
     8// science image, we construct a set of simulated sources with their own pixels.
     9// these are used to determine the simultaneous linear fit of fluxes.
    1010// the analysis is performed wrt the simulated pixel values
    1111
    1212static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER);
    1313
    14 bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
     14bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final, psMaskType maskVal) {
    1515
    1616    bool status;
    … …  
    4141    int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
    4242    if (!status) {
    43         SKY_FIT_ORDER = 0;
     43        SKY_FIT_ORDER = 0;
    4444    }
    4545    bool SKY_FIT_LINEAR = psMetadataLookupBool(&status, recipe, "SKY_FIT_LINEAR");
    4646    if (!status) {
    47         SKY_FIT_LINEAR = false;
     47        SKY_FIT_LINEAR = false;
    4848    }
    4949
    … …  
    5555        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    5656        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    57         if (source->type == PM_SOURCE_TYPE_STAR &&
     57        if (source->type == PM_SOURCE_TYPE_STAR &&
    5858            source->mode & PM_SOURCE_MODE_SATSTAR) continue;
    5959        if (final) {
    … …  
    6767        y = source->peak->yf;
    6868
    69         // is the source in the region of interest?
     69        // is the source in the region of interest?
    7070        if (x < AnalysisRegion.x0) continue;
    7171        if (y < AnalysisRegion.y0) continue;
    … …  
    9797        psSparseMatrixElement (sparse, i, i, f);
    9898
    99         // the formal error depends on the weighting scheme
    100         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
    101             float var = pmSourceModelDotModel (SRCi, SRCi, false);
    102             errors->data.F32[i] = 1.0 / sqrt(var);
    103         } else {
    104             errors->data.F32[i] = 1.0 / sqrt(f);
    105         }
     99        // the formal error depends on the weighting scheme
     100        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
     101            float var = pmSourceModelDotModel (SRCi, SRCi, false);
     102            errors->data.F32[i] = 1.0 / sqrt(var);
     103        } else {
     104            errors->data.F32[i] = 1.0 / sqrt(f);
     105        }
    106106
    107107
    … …  
    110110        psSparseVectorElement (sparse, i, f);
    111111
    112         // add the per-source weights (border region)
    113         switch (SKY_FIT_ORDER) {
    114           case 1:
    115             f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS);
    116             psSparseBorderElementB (border, i, 1, f);
    117             f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS);
    118             psSparseBorderElementB (border, i, 2, f);
    119 
    120           case 0:
    121             f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
    122             psSparseBorderElementB (border, i, 0, f);
    123             break;
    124 
    125           default:
    126             psAbort("Invalid SKY_FIT_ORDER %d\n", SKY_FIT_ORDER);
    127             break;
    128         }
     112        // add the per-source weights (border region)
     113        switch (SKY_FIT_ORDER) {
     114          case 1:
     115            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS);
     116            psSparseBorderElementB (border, i, 1, f);
     117            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS);
     118            psSparseBorderElementB (border, i, 2, f);
     119
     120          case 0:
     121            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
     122            psSparseBorderElementB (border, i, 0, f);
     123            break;
     124
     125          default:
     126            psAbort("Invalid SKY_FIT_ORDER %d\n", SKY_FIT_ORDER);
     127            break;
     128        }
    129129
    130130        // loop over all other stars following this one
    … …  
    159159    psVector *skyfit = NULL;
    160160    if (SKY_FIT_LINEAR) {
    161         psSparseBorderSolve (&norm, &skyfit, constraint, border, 5);
    162         fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]);
     161        psSparseBorderSolve (&norm, &skyfit, constraint, border, 5);
     162        fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]);
    163163    } else {
    164         norm = psSparseSolve (NULL, constraint, sparse, 5);
    165         skyfit = NULL;
     164        norm = psSparseSolve (NULL, constraint, sparse, 5);
     165        skyfit = NULL;
    166166    }
    167167
    … …  
    169169    for (int i = 0; i < fitSources->n; i++) {
    170170        pmSource *source = fitSources->data[i];
    171         pmModel *model = pmSourceGetModel (NULL, source);
     171        pmModel *model = pmSourceGetModel (NULL, source);
    172172
    173173        // assign linearly-fitted normalization
    … …  
    177177        model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
    178178        model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
    179         // XXX is the value of 'errors' modified by the sky fit?
     179        // XXX is the value of 'errors' modified by the sky fit?
    180180
    181181        // subtract object
    182         pmSourceSub (source, PM_MODEL_OP_FULL);
     182        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    183183        source->mode |= PM_SOURCE_MODE_SUBTRACTED;
    184184        if (!final) source->mode |= PM_SOURCE_MODE_TEMPSUB;
    185         // XXX not sure about the use of TEMPSUB
     185        // XXX not sure about the use of TEMPSUB
    186186    }
    187187
    … …  
    190190        pmSource *source = fitSources->data[i];
    191191        pmModel *model = pmSourceGetModel (NULL, source);
    192         pmSourceChisq (model, source->pixels, source->maskObj, source->weight);
     192        pmSourceChisq (model, source->pixels, source->maskObj, source->weight, maskVal);
    193193    }
    194194
    … …  
    219219    for (int i = 0; i < sources->n; i++) {
    220220        pmSource *source = sources->data[i];
    221         pmModel *model = pmSourceGetModel (NULL, source);
    222         if (model == NULL) continue;
     221        pmModel *model = pmSourceGetModel (NULL, source);
     222        if (model == NULL) continue;
    223223        float x = model->params->data.F32[PM_PAR_XPOS];
    224224        float y = model->params->data.F32[PM_PAR_YPOS];
    225         psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
    226     }   
     225        psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
     226    }
    227227
    228228    // accumulate the image statistics from the masked regions
    … …  
    233233    double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy;
    234234    w = x = y = x2 = xy = y2 = fo = fx = fy = 0;
    235    
     235
    236236    int col0 = readout->image->col0;
    237237    int row0 = readout->image->row0;
    238238
    239239    for (int j = 0; j < readout->image->numRows; j++) {
    240         for (int i = 0; i < readout->image->numCols; i++) {
    241             if (mask[j][i]) continue;
    242             if (constant_weights) {
    243                 wt = 1.0;
    244             } else {
    245                 wt = weight[j][i];
    246             }
    247             f = image[j][i];
    248             w   += 1/wt;
    249             fo  += f/wt;
    250             if (SKY_FIT_ORDER == 0) continue;
    251 
    252             xc  = i + col0;
    253             yc  = j + row0;
    254             x  +=    xc/wt;
    255             y  +=    yc/wt;
    256             x2 += xc*xc/wt;
    257             xy += xc*yc/wt;
    258             y2 += yc*yc/wt;
    259             fx +=  f*xc/wt;
    260             fy +=  f*yc/wt;
    261         }
     240        for (int i = 0; i < readout->image->numCols; i++) {
     241            if (mask[j][i]) continue;
     242            if (constant_weights) {
     243                wt = 1.0;
     244            } else {
     245                wt = weight[j][i];
     246            }
     247            f = image[j][i];
     248            w   += 1/wt;
     249            fo  += f/wt;
     250            if (SKY_FIT_ORDER == 0) continue;
     251
     252            xc  = i + col0;
     253            yc  = j + row0;
     254            x  +=    xc/wt;
     255            y  +=    yc/wt;
     256            x2 += xc*xc/wt;
     257            xy += xc*yc/wt;
     258            y2 += yc*yc/wt;
     259            fx +=  f*xc/wt;
     260            fy +=  f*yc/wt;
     261        }
    262262    }
    263263
    … …  
    269269    psSparseBorderElementT (border, 0, 0, w);
    270270    if (SKY_FIT_ORDER > 0) {
    271         psSparseBorderElementG (border, 0, fx);
    272         psSparseBorderElementG (border, 0, fy);
    273         psSparseBorderElementT (border, 1, 0, x);
    274         psSparseBorderElementT (border, 2, 0, y);
    275         psSparseBorderElementT (border, 0, 1, x);
    276         psSparseBorderElementT (border, 1, 1, x2);
    277         psSparseBorderElementT (border, 2, 1, xy);
    278         psSparseBorderElementT (border, 0, 2, y);
    279         psSparseBorderElementT (border, 1, 2, xy);
    280         psSparseBorderElementT (border, 2, 2, y2);
    281     }   
     271        psSparseBorderElementG (border, 0, fx);
     272        psSparseBorderElementG (border, 0, fy);
     273        psSparseBorderElementT (border, 1, 0, x);
     274        psSparseBorderElementT (border, 2, 0, y);
     275        psSparseBorderElementT (border, 0, 1, x);
     276        psSparseBorderElementT (border, 1, 1, x2);
     277        psSparseBorderElementT (border, 2, 1, xy);
     278        psSparseBorderElementT (border, 0, 2, y);
     279        psSparseBorderElementT (border, 1, 2, xy);
     280        psSparseBorderElementT (border, 2, 2, y2);
     281    }
    282282    return true;
    283283}
Note: See TracChangeset for help on using the changeset viewer.