IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 17, 2006, 1:04:04 PM (20 years ago)
Author:
eugene
Message:

important code re-organization: split psphotEnsemblePSF into two functions, added general model capability, added sky offset fit

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotEnsemblePSF.c

    r10032 r10052  
    11# include "psphot.h"
    22
    3 // fit all reasonable sources with the linear PSF model
    4 bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
    5 
    6     bool  status;
     3// fit flux (and optionally sky model) to all reasonable sources
     4// with the linear fitting process.  sources must have an associated
     5// model with selected pixels, and the fit radius must be defined
     6
     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. 
     10// the analysis is performed wrt the simulated pixel values
     11
     12static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, bool fitSlope);
     13
     14bool psphotEnsemblePSF (pmReadout *readout, psArray *refSources, psMetadata *recipe, pmPSF *psf, bool final) {
     15
     16    bool isPSF;
     17    bool status;
    718    float x;
    819    float y;
    … …  
    1021    float r;
    1122
    12     // psRegion allArray = psRegionSet (0, 0, 0, 0);
    13 
    1423    psTimerStart ("psphot");
    1524
    1625    // source analysis is done in spatial order
    17     sources = psArraySort (sources, psphotSortByY);
    18 
    19     // radius for source moments if star is saturated
    20     // XXX we already do this in psphotSourceStats, is it still needed?
    21     float INNER_RADIUS     = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    22 
    23     // setup the PSF fit details
    24     psphotInitRadiusPSF (recipe, psf->type);
    25 
    26     // pre-calculate all model pixels
    27     psArray  *models = psArrayAllocEmpty (sources->n);
    28     psVector *index  = psVectorAllocEmpty (sources->n, PS_TYPE_U32);
     26    refSources = psArraySort (refSources, psphotSortByY);
     27
     28    // storage arrays for fitSources and sequence index
     29    psArray *fitSources = psArrayAllocEmpty (refSources->n);
     30    psVector *index = psVectorAllocEmpty (refSources->n, PS_TYPE_U32);
    2931
    3032    // option to limit analysis to a specific region
    … …  
    3234    psRegion AnalysisRegion = psRegionFromString (region);
    3335    AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion);
    34     // psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    3536    if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined");
    3637
    … …  
    4142    }
    4243
    43     for (int i = 0; i < sources->n; i++) {
    44         pmSource *inSource = sources->data[i];
     44    // construct source table with stand-alone fitSources
     45    for (int i = 0; i < refSources->n; i++) {
     46        pmSource *refSource = refSources->data[i];
    4547
    4648        // skip non-astronomical objects (very likely defects)
    47         // XXX EAM : should we try these anyway?
    48         if (inSource->type == PM_SOURCE_TYPE_DEFECT) continue;
    49         if (inSource->type == PM_SOURCE_TYPE_SATURATED) continue;
     49        if (refSource->type == PM_SOURCE_TYPE_DEFECT) continue;
     50        if (refSource->type == PM_SOURCE_TYPE_SATURATED) continue;
    5051        if (final) {
    51             if (inSource->mode &  PM_SOURCE_MODE_SUBTRACTED) continue;
     52            if (refSource->mode &  PM_SOURCE_MODE_SUBTRACTED) continue;
    5253        } else {
    53             if (inSource->mode &  PM_SOURCE_MODE_BLEND) continue;
     54            if (refSource->mode &  PM_SOURCE_MODE_BLEND) continue;
    5455        }
    5556
    56         if (inSource->moments->x < AnalysisRegion.x0) continue;
    57         if (inSource->moments->y < AnalysisRegion.y0) continue;
    58         if (inSource->moments->x > AnalysisRegion.x1) continue;
    59         if (inSource->moments->y > AnalysisRegion.y1) continue;
    60 
    61         pmSource *otSource = pmSourceAlloc ();
    62 
    63         // really saturated stars should be re-measured for a better centroid
    64         // XXX EAM : move this to a 'clear satstar function'
    65         // XXX : place this in pmSourceMoments or wherever it is called?
    66         if (inSource->mode &  PM_SOURCE_MODE_SATSTAR) {
    67             status = pmSourceMoments (inSource, INNER_RADIUS);
    68         }
    69 
    70         // XXX EAM : add option to use EXT or PSF form
    71         // use the source moments, etc to guess basic model parameters
    72         pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type);
    73         if (inSource->mode &  PM_SOURCE_MODE_SATSTAR) {
    74             modelEXT->params->data.F32[PM_PAR_XPOS] = inSource->moments->x;
    75             modelEXT->params->data.F32[PM_PAR_YPOS] = inSource->moments->y;
    76         } else {
    77             // peak-up on peak (for non-sat objects)
    78 
    79             // ix,iy must land on inSource->pixels
    80             int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x));
    81             int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y));
    82 
    83             psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
    84             psPlane min = psImageBicubeMin (bicube);
    85 
    86             psTrace ("psphot", 5, "peak coord: %f %f -> %f %f\n",
    87                      modelEXT->params->data.F32[PM_PAR_XPOS], modelEXT->params->data.F32[PM_PAR_YPOS], min.x + ix, min.y + iy);
    88 
    89             // if min point is too deviant, keep the old value
    90             if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
    91                 modelEXT->params->data.F32[PM_PAR_XPOS] = min.x + ix;
    92                 modelEXT->params->data.F32[PM_PAR_YPOS] = min.y + iy;
    93             }
    94             psFree (bicube);
    95         }
    96 
    97         // set PSF parameters for this model
    98         pmModel *model = pmModelFromPSF (modelEXT, psf);
    99         psFree (modelEXT);
     57        // which model to use?
     58        pmModel *model = pmSourceGetModel (&isPSF, refSource);
     59        if (model == NULL) continue;  // model must be defined
    10060
    10161        // save the original coords
    … …  
    10363        y = model->params->data.F32[PM_PAR_YPOS];
    10464
    105         // set the fit radius based on the object flux limit and the model
    106         psphotCheckRadiusPSF (readout, inSource, model);
     65        // is the source in the region of interest?
     66        if (x < AnalysisRegion.x0) continue;
     67        if (y < AnalysisRegion.y0) continue;
     68        if (x > AnalysisRegion.x1) continue;
     69        if (y > AnalysisRegion.y1) continue;
     70
     71        pmSource *fitSource = pmSourceAlloc ();
    10772
    10873        // make temporary copies of the image pixels and mask
    10974        // we need to have a copy which will not be modified by changes to its neighbor
    110         otSource->mask   = psImageCopy (NULL, inSource->mask,   PS_TYPE_U8);
    111         otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32);
    112         otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32);
    113 
    114         // build the model image
    115         psImage *flux = otSource->pixels;
    116         psImage *mask = otSource->mask;
     75        fitSource->mask   = psImageCopy (NULL, refSource->mask,   PS_TYPE_U8);
     76        fitSource->pixels = psImageCopy (NULL, refSource->pixels, PS_TYPE_F32);
     77        fitSource->weight = psImageCopy (NULL, refSource->weight, PS_TYPE_F32);
    11778
    11879        // set model to unit peak, zero sky (we assume sky is subtracted)
     80        // XXX do this with the actual model or use keep the original guess?
    11981        model->params->data.F32[PM_PAR_SKY] = 0.0;
    12082        model->params->data.F32[PM_PAR_I0] = 1.0;
    12183
    12284        // fill in the model pixel values
    123         psImageInit (flux, 0.0);
    124         psImageKeepCircle (mask, x, y, model->radiusFit, "OR", PM_MASK_MARK);
    125         pmModelAdd (flux, mask, model, false, false);
    126 
    127         // calculate nDOF (nPix - 1)
    128         // int Nmaskpix = psImageCountPixelMask (mask, allArray, PM_MASK_SAT);
    129         // model->nDOF  = mask->numCols*mask->numRows - Nmaskpix - 1;
     85        psImageInit (fitSource->pixels, 0.0);
     86        psImageKeepCircle (fitSource->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK);
     87        pmModelAdd (fitSource->pixels, fitSource->mask, model, false, false);
    13088
    13189        // save source in list
    132         otSource->modelPSF = model;
    133         index->data.U32[models->n] = i;
    134         psArrayAdd (models, 100, otSource);
    135         psFree (otSource);
    136     }
    137     psLogMsg ("psphot.emsemble", 4, "built models: %f (%ld objects)\n", psTimerMark ("psphot"), sources->n);
    138 
    139     // fill out the sparse matrix
    140     psSparse *sparse = psSparseAlloc (models->n, 100);
    141     psVector *weight = psVectorAlloc (models->n, PS_TYPE_F32);
    142     psVector *errors = psVectorAlloc (models->n, PS_TYPE_F32);
     90        if (isPSF) {
     91          fitSource->modelPSF = psMemIncrRefCounter (model);
     92        } else {
     93          fitSource->modelEXT = psMemIncrRefCounter (model);
     94        }
     95        index->data.U32[fitSources->n] = i;
     96        psArrayAdd (fitSources, 100, fitSource);
     97        psFree (fitSource);
     98    }
     99    psLogMsg (__func__, 4, "built fitSources: %f (%ld objects)\n", psTimerMark ("psphot"), refSources->n);
     100
     101    // vectors to store stats for each object
     102    psVector *weight = psVectorAlloc (fitSources->n, PS_TYPE_F32);
     103    psVector *errors = psVectorAlloc (fitSources->n, PS_TYPE_F32);
    143104
    144105    // create the border matrix (includes the sparse matrix)
    145106    // for just sky: 1 row; for x,y terms: 3 rows
     107    psSparse *sparse = psSparseAlloc (fitSources->n, 100);
    146108    psSparseBorder *border = psSparseBorderAlloc (sparse, 1);
    147109
    148     for (int i = 0; i < models->n; i++) {
     110    // fill out the sparse matrix elements and border elements (B)
     111    // Ri is the current refSource of interest (sci pixels)
     112    // Fi is the current fitSource of interest (fit pixels)
     113    for (int i = 0; i < fitSources->n; i++) {
    149114        int N = index->data.U32[i];
    150         pmSource *Fi = sources->data[N];
    151         pmSource *Mi = models->data[i];
     115        pmSource *Ri = refSources->data[N];
     116        pmSource *Fi = fitSources->data[i];
    152117
    153118        // scale by diagonal element (auto-cross-product)
    154         r = pmSourceCrossProduct (Mi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS);
     119        r = pmSourceCrossProduct (Fi, Fi, CONSTANT_PHOTOMETRIC_WEIGHTS);
    155120        weight->data.F32[i] = r;
    156121
    157122        // the formal error depends on the weighting scheme
    158123        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
    159             float var = pmSourceCrossProduct (Mi, Mi, false);
     124            float var = pmSourceCrossProduct (Fi, Fi, false);
    160125            errors->data.F32[i] = 1.0 / sqrt(var);
    161126        } else {
    … …  
    166131
    167132        // find the image x model value
    168         f = pmSourceCrossProduct (Fi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS);
     133        f = pmSourceCrossProduct (Ri, Fi, CONSTANT_PHOTOMETRIC_WEIGHTS);
    169134        psSparseVectorElement (sparse, i, f / r);
    170135
    171136        // add the per-source weights (border region)
    172         p = pmSourceWeight (Mi, CONST, CONSTANT_PHOTOMETRIC_WEIGHTS);
    173         // px = pmSourceWeight (Mi, XTERM, CONSTANT_PHOTOMETRIC_WEIGHTS);
    174         // py = pmSourceWeight (Mi, YTERM, CONSTANT_PHOTOMETRIC_WEIGHTS);
     137        float p = pmSourceWeight (Fi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
     138        // float px = pmSourceWeight (Fi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS);
     139        // float py = pmSourceWeight (Fi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS);
    175140        psSparseBorderElementB (border, i, 0, p);
    176141        // psSparseBorderElementB (border, i, 1, px);
    … …  
    178143
    179144        // loop over all other stars following this one
    180         for (int j = i + 1; j < models->n; j++) {
    181             pmSource *Mj = models->data[j];
     145        for (int j = i + 1; j < fitSources->n; j++) {
     146            pmSource *Fj = fitSources->data[j];
    182147
    183148            // skip over disjoint source images, break after last possible overlap
    184             if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break;
    185             if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue;
    186             if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue;
    187             if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue;
     149            if (Fi->pixels->row0 + Fi->pixels->numRows < Fj->pixels->row0) break;
     150            if (Fj->pixels->row0 + Fj->pixels->numRows < Fi->pixels->row0) continue;
     151            if (Fi->pixels->col0 + Fi->pixels->numCols < Fj->pixels->col0) continue;
     152            if (Fj->pixels->col0 + Fj->pixels->numCols < Fi->pixels->col0) continue;
    188153
    189154            // got an overlap; calculate cross-product and add to output array
    190             f = pmSourceCrossProduct (Mi, Mj, CONSTANT_PHOTOMETRIC_WEIGHTS);
     155            f = pmSourceCrossProduct (Fi, Fj, CONSTANT_PHOTOMETRIC_WEIGHTS);
    191156            psSparseMatrixElement (sparse, j, i, f / r);
    192157        }
    193158    }
    194     psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
     159    psSparseResort (sparse);
     160    psLogMsg (__func__, 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
     161
     162    // set the sky, sky_x, sky_y components of border matrix
     163    // XXX ignore sky slope terms for now
     164    SetBorderT (border, readout, refSources, false);
     165
     166    psSparseConstraint constraint;
     167    constraint.paramMin   = 0.0;
     168    constraint.paramMax   = 1e8;
     169    constraint.paramDelta = 1e8;
     170
     171    // solve for normalization terms (need include local sky?)
     172    psVector *norm = NULL;
     173    psVector *skyfit = NULL;
     174    psSparseBorderSolve (&norm, &skyfit, constraint, border, 3);
     175
     176    // adjust fitSources, set refSources and subtract
     177    for (int i = 0; i < fitSources->n; i++) {
     178        int N = index->data.U32[i];
     179        pmSource *Ri = refSources->data[N];
     180        pmModel *model = pmSourceGetModel (NULL, Ri);
     181
     182        // assign linearly-fitted normalization
     183        if (isnan(norm->data.F32[i])) {
     184            psAbort ("psphot", "ensemble source is nan");
     185        }
     186        model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
     187        // XXX note: is the value of 'errors' modified by the sky fit?
     188        model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
     189
     190        // subtract object
     191        pmModelSub (Ri->pixels, Ri->mask, model, false, false);
     192        Ri->mode |= PM_SOURCE_MODE_SUBTRACTED;
     193        if (!final) Ri->mode |= PM_SOURCE_MODE_TEMPSUB;
     194        // XXX not sure about the use of TEMPSUB
     195    }
     196
     197    // measure chisq for each source
     198    for (int i = 0; final && (i < fitSources->n); i++) {
     199        int N = index->data.U32[i];
     200        pmSource *Ri = refSources->data[N];
     201        pmModel *model = pmSourceGetModel (NULL, Ri);
     202
     203        x = model->params->data.F32[PM_PAR_XPOS];
     204        y = model->params->data.F32[PM_PAR_YPOS];
     205
     206        psImageKeepCircle (Ri->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK);
     207        pmSourceChisq (model, Ri->pixels, Ri->mask, Ri->weight);
     208        psImageKeepCircle (Ri->mask, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
     209    }
     210
     211    psFree (index);
     212    psFree (sparse);
     213    psFree (fitSources);
     214    psFree (norm);
     215    psFree (weight);
     216    psFree (errors);
     217
     218    psLogMsg (__func__, 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot"));
     219    return true;
     220}
     221
     222// calculate the weight terms for the sky fit component of the matrix
     223static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, bool fitSlope) {
    195224
    196225    // generate the image-wide weight terms
    197     // first, mask all object pixels
    198     for (int i = 0; i < models->n; i++) {
    199         int N = index->data.U32[i];
    200         pmSource *Fi = sources->data[N];
    201         pmSource *Mi = models->data[i];
    202 
    203         // XXX need to select the correct mode (check this against the selections below
    204         // - which model is appropriate for the given source?
    205         // - which radius is appropriate for the given source?
    206         // - where do x,y come from?
    207         psImageKeepCircle (Mi->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK);
     226    // turn on MARK for all image pixels
     227    psRegion fullArray = psRegionSet (0, 0, 0, 0);
     228    fullArray = psRegionForImage (readout->mask, fullArray);
     229    psImageMaskRegion (readout->mask, fullArray, "OR", PM_MASK_MARK);
     230
     231    // turn off MARK for all object pixels
     232    // we must use the refSources here since their pixels point at the science image
     233    for (int i = 0; i < refSources->n; i++) {
     234        pmSource *refSource = refSources->data[i];
     235        pmModel *model = pmSourceGetModel (NULL, refSource);
     236        float x = model->params->data.F32[PM_PAR_XPOS];
     237        float y = model->params->data.F32[PM_PAR_YPOS];
     238        psImageMaskCircle (refSource->mask, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
    208239    }   
    209240
    210241    // accumulate the image statistics from the masked regions
    211     for (int j = 0; j < numRows; j++) {
    212         for (int i = 0; i < numCols; i++) {
    213             w  +=     image->data.F32[j][i]/weight->data.F32[j][i];
    214             x  +=   i*image->data.F32[j][i]/weight->data.F32[j][i];
    215             y  +=   j*image->data.F32[j][i]/weight->data.F32[j][i];
    216             x2 += i*i*image->data.F32[j][i]/weight->data.F32[j][i];
    217             xy += i*j*image->data.F32[j][i]/weight->data.F32[j][i];
    218             y2 += j*j*image->data.F32[j][i]/weight->data.F32[j][i];
     242    psF32 **image = readout->image->data.F32;
     243    psF32 **weight = readout->weight->data.F32;
     244    psU8 **mask = readout->mask->data.U8;
     245
     246    double w, x, y, x2, xy, y2, xc, yc;
     247    w = x = y = x2 = xy = y2 = 0;
     248   
     249    int col0 = readout->image->col0;
     250    int row0 = readout->image->row0;
     251
     252    for (int j = 0; j < readout->image->numRows; j++) {
     253        for (int i = 0; i < readout->image->numCols; i++) {
     254            if (mask[j][i]) continue;
     255            w  +=       image[j][i]/weight[j][i];
     256            if (!fitSlope) continue;
     257            xc = i + col0;
     258            yc = j + row0;
     259            x  +=    xc*image[j][i]/weight[j][i];
     260            y  +=    yc*image[j][i]/weight[j][i];
     261            x2 += xc*xc*image[j][i]/weight[j][i];
     262            xy += xc*yc*image[j][i]/weight[j][i];
     263            y2 += yc*yc*image[j][i]/weight[j][i];
    219264        }
    220265    }
     266
     267    // turn off MARK for all image pixels
     268    psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_U8(PM_MASK_MARK));
     269
     270    // set the Border T elements
    221271    psSparseBorderElementT (border, 0, 0, w);
    222     if (0) {
     272    if (fitSlope) {
    223273        psSparseBorderElementT (border, 1, 0, x);
    224274        psSparseBorderElementT (border, 2, 0, y);
    … …  
    230280        psSparseBorderElementT (border, 2, 2, y2);
    231281    }   
    232 
    233     psSparseConstraint constraint;
    234     constraint.paramMin   = 0.0;
    235     constraint.paramMax   = 1e8;
    236     constraint.paramDelta = 1e8;
    237 
    238     // solve for normalization terms (need include local sky?)
    239     psSparseResort (sparse);
    240     psSparseBorderSolve (&xFit, &yFit, constraint, border, 3);
    241 
    242     // psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);
    243 
    244     // adjust models, set sources and subtract
    245     for (int i = 0; i < models->n; i++) {
    246         int N = index->data.U32[i];
    247         pmSource *Fi = sources->data[N];
    248         pmSource *Mi = models->data[i];
    249 
    250         // if we already have a PSF model, free it.
    251         psFree (Fi->modelPSF);
    252 
    253         // need to increment counter so we can free models here and sources above
    254         Fi->modelPSF = psMemIncrRefCounter (Mi->modelPSF);
    255 
    256         // assign linearly-fitted normalization
    257         if (isnan(norm->data.F32[i])) {
    258             psAbort ("psphot", "ensemble source is nan");
    259         }
    260         Fi->modelPSF->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
    261         // XXX I was using this before, but it is not strictly accurate
    262         // XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it
    263         // Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i]));
    264         Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
    265 
    266         // subtract object
    267         pmModelSub (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
    268         Fi->mode |= PM_SOURCE_MODE_SUBTRACTED;
    269         if (!final) Fi->mode |= PM_SOURCE_MODE_TEMPSUB;
    270     }
    271 
    272     // measure chisq for each source
    273     for (int i = 0; final && (i < models->n); i++) {
    274         int N = index->data.U32[i];
    275         pmSource *Fi = sources->data[N];
    276         pmModel *model = Fi->modelPSF;
    277 
    278         x = model->params->data.F32[PM_PAR_XPOS];
    279         y = model->params->data.F32[PM_PAR_YPOS];
    280 
    281         psImageKeepCircle (Fi->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK);
    282         pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight);
    283         psImageKeepCircle (Fi->mask, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
    284     }
    285 
    286     psFree (index);
    287     psFree (sparse);
    288     psFree (models);
    289     psFree (norm);
    290     psFree (weight);
    291     psFree (errors);
    292 
    293     psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot"));
    294282    return true;
    295283}
Note: See TracChangeset for help on using the changeset viewer.