IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 6, 2010, 1:36:51 PM (16 years ago)
Author:
eugene
Message:

working on stackphot

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/stackphot.20100406/psphot/src/psphotFitSourcesLinearStack.c

    r27547 r27625  
    11# include "psphotInternal.h"
    2 
    3 // for now, let's store the detections on the readout->analysis for each readout
    4 bool psphotFitSourcesLinearStack (pmConfig *config, const pmFPAview *view, bool final)
    5 {
    6     bool status = true;
    7 
    8     int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    9     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    10 
    11     // loop over the available readouts
    12     for (int i = 0; i < num; i++) {
    13         if (!psphotFitSourcesLinearReadoutStack (config, view, "PSPHOT.INPUT", i, final)) {
    14             psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i);
    15             return false;
    16         }
    17     }
    18     return true;
    19 }
    20 
    21 bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final) {
     2// XXX need array of covar factors for each image
     3// XXX define the 'good' / 'bad' flags?
     4
     5bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final) {
    226
    237    bool status;
    … …  
    259    float y;
    2610    float f;
    27     // float r;
    2811
    2912    // select the appropriate recipe information
    … …  
    3114    assert (recipe);
    3215
    33     // find the currently selected readout
    34     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
    35     psAssert (file, "missing file?");
    36 
    37     pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    38     psAssert (readout, "missing readout?");
    39 
    40     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    41     psAssert (detections, "missing detections?");
    42 
    43     psArray *sources = detections->allSources;
    44     psAssert (sources, "missing sources?");
    45 
    46     if (!sources->n) {
    47         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit");
    48         return true;
    49     }
    50 
    51     pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    52     psAssert (sources, "missing psf?");
    53 
    5416    psTimerStart ("psphot.linear");
    5517
    … …  
    6527    maskVal |= markVal;
    6628
    67     // source analysis is done in spatial order
    68     sources = psArraySort (sources, pmSourceSortByY);
     29    // analysis is done in spatial order (to speed up overlap search)
     30    // sort by first element in each source list
     31    objects = psArraySort (objects, pmPhotObjSortByY);
    6932
    7033    // storage array for fitSources
    71     psArray *fitSources = psArrayAllocEmpty (sources->n);
    72 
    73     // option to limit analysis to a specific region
    74     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    75     psRegion AnalysisRegion = psRegionFromString (region);
    76     AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion);
    77     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
    78 
    79     bool CONSTANT_PHOTOMETRIC_WEIGHTS =
    80         psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
    81     if (!status) {
    82         psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
    83     }
    84     int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
    85     if (!status) {
    86         SKY_FIT_ORDER = 0;
    87     }
    88     bool SKY_FIT_LINEAR = psMetadataLookupBool(&status, recipe, "SKY_FIT_LINEAR");
    89     if (!status) {
    90         SKY_FIT_LINEAR = false;
    91     }
    92 
    93     // XXX test: choose a larger-than expected radius:
    94     float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
    95     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor);
    96 
    97     // XXX do not apply covarFactor for the moment...
    98     // covarFactor = 1.0;
     34    psArray *fitSources = psArrayAllocEmpty (objects->n);
     35
     36    bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
     37    psAssert (status, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
     38
     39    // XXX store a local static array of covar factors for each of the images (by image ID)
     40    // float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
     41    // psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor);
    9942
    10043    // select the sources which will be used for the fitting analysis
    … …  
    10447        if (!object->sources) continue;
    10548
    106         // check an element of the group to see if we should use it
    107         if (!object->flags & PM_PHOT_OBJ_BAD) continue;
     49        // XXX check an element of the group to see if we should use it
     50        // if (!object->flags & PM_PHOT_OBJ_BAD) continue;
    10851
    10952        for (int j = 0; j < object->sources->n; j++) {
    … …  
    162105            pmSource *SRCj = fitSources->data[j];
    163106
    164             // XXX I need to know if this source is on the same image as SRCi --
    165             if (!sameImge) { continue; }
     107            // we only need to generate dot terms for source on the same image
     108            if (SRCj->imageID != SRCi->imageID) { continue; }
    166109
    167110            // skip over disjoint source images, break after last possible overlap
    … …  
    190133    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
    191134
    192     // XXXX **** philosophical question:
    193     // we measure bright objects in three passes: 1) linear fit; 2) non-linear fit; 3) linear fit:
    194     // should retain the chisq and errors from the intermediate non-linear fit?
    195     // the non-linear fit provides better values for the position errors, and for
    196     // extended sources, the shape errors
    197 
    198135    // adjust I0 for fitSources and subtract
    199136    for (int i = 0; i < fitSources->n; i++) {
    … …  
    208145        model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
    209146        model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
    210         // XXX is the value of 'errors' modified by the sky fit?
    211147
    212148        // subtract object
    … …  
    229165    psFree (fitSources);
    230166    psFree (norm);
    231     psFree (skyfit);
    232167    psFree (errors);
    233     psFree (border);
    234168
    235169    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
    236 
    237     psphotVisualShowResidualImage (readout);
    238     psphotVisualPlotChisq (sources);
    239     // psphotVisualShowFlags (sources);
    240 
    241     // We have to place this visualization here because the models are not realized until
    242     // psphotGuessModels or fitted until psphotFitSourcesLinear.
    243     psphotVisualShowPSFStars (recipe, psf, sources);
    244170
    245171    return true;
    246172}
    247173
    248 // XXX do we need this?
    249 // XXX disallow the simultaneous sky fit and remove this code...
    250 
    251 // Calculate the weight terms for the sky fit component of the matrix.  This function operates
    252 // on the pixels which correspond to all of the sources of interest.  These elements fill in
    253 // the border matrix components in the sparse matrix equation.
    254 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal) {
    255 
    256     // generate the image-wide weight terms
    257     // turn on MARK for all image pixels
    258     psRegion fullArray = psRegionSet (0, 0, 0, 0);
    259     fullArray = psRegionForImage (readout->mask, fullArray);
    260     psImageMaskRegion (readout->mask, fullArray, "OR", markVal);
    261 
    262     // turn off MARK for all object pixels
    263     for (int i = 0; i < sources->n; i++) {
    264         pmSource *source = sources->data[i];
    265         pmModel *model = pmSourceGetModel (NULL, source);
    266         if (model == NULL) continue;
    267         float x = model->params->data.F32[PM_PAR_XPOS];
    268         float y = model->params->data.F32[PM_PAR_YPOS];
    269         psImageMaskCircle (source->maskView, x, y, model->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal));
    270     }
    271 
    272     // accumulate the image statistics from the masked regions
    273     psF32 **image  = readout->image->data.F32;
    274     psF32 **variance = readout->variance->data.F32;
    275     psImageMaskType  **mask   = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA;
    276 
    277     double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy;
    278     w = x = y = x2 = xy = y2 = fo = fx = fy = 0;
    279 
    280     int col0 = readout->image->col0;
    281     int row0 = readout->image->row0;
    282 
    283     for (int j = 0; j < readout->image->numRows; j++) {
    284         for (int i = 0; i < readout->image->numCols; i++) {
    285             if (mask[j][i]) continue;
    286             if (constant_weights) {
    287                 wt = 1.0;
    288             } else {
    289                 wt = variance[j][i];
    290             }
    291             f = image[j][i];
    292             w   += 1/wt;
    293             fo  += f/wt;
    294             if (SKY_FIT_ORDER == 0) continue;
    295 
    296             xc  = i + col0;
    297             yc  = j + row0;
    298             x  +=    xc/wt;
    299             y  +=    yc/wt;
    300             x2 += xc*xc/wt;
    301             xy += xc*yc/wt;
    302             y2 += yc*yc/wt;
    303             fx +=  f*xc/wt;
    304             fy +=  f*yc/wt;
    305         }
    306     }
    307 
    308     // turn off MARK for all image pixels
    309     psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_IMAGE_MASK(markVal));
    310 
    311     // set the Border T elements
    312     psSparseBorderElementG (border, 0, fo);
    313     psSparseBorderElementT (border, 0, 0, w);
    314     if (SKY_FIT_ORDER > 0) {
    315         psSparseBorderElementG (border, 0, fx);
    316         psSparseBorderElementG (border, 0, fy);
    317         psSparseBorderElementT (border, 1, 0, x);
    318         psSparseBorderElementT (border, 2, 0, y);
    319         psSparseBorderElementT (border, 0, 1, x);
    320         psSparseBorderElementT (border, 1, 1, x2);
    321         psSparseBorderElementT (border, 2, 1, xy);
    322         psSparseBorderElementT (border, 0, 2, y);
    323         psSparseBorderElementT (border, 1, 2, xy);
    324         psSparseBorderElementT (border, 2, 2, y2);
    325     }
    326 
    327     return true;
     174// sort by Y (ascending)
     175int pmPhotObjSortBySN (const void **a, const void **b)
     176{
     177    pmPhotObj *objA = *(pmPhotObj **)a;
     178    pmPhotObj *objB = *(pmPhotObj **)b;
     179
     180    psAssert (objA->sources, "missing sources?");
     181    psAssert (objB->sources, "missing sources?");
     182
     183    psAssert (objA->sources->n, "missing sources");
     184    psAssert (objB->sources->n, "missing sources");
     185
     186    psAssert (objA->sources->data[0], "missing sources");
     187    psAssert (objB->sources->data[0], "missing sources");
     188
     189    pmSource *A = objA->sources->data[0];
     190    pmSource *B = objB->sources->data[0];
     191
     192    psF32 fA = (A->peak == NULL) ? 0 : A->peak->y;
     193    psF32 fB = (B->peak == NULL) ? 0 : B->peak->y;
     194
     195    psF32 diff = fA - fB;
     196    if (diff > FLT_EPSILON) return (+1);
     197    if (diff < FLT_EPSILON) return (-1);
     198    return (0);
    328199}
     200
Note: See TracChangeset for help on using the changeset viewer.