IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 9527


Ignore:
Timestamp:
Oct 12, 2006, 4:24:34 PM (20 years ago)
Author:
rhl
Message:

Use CONSTANT_PHOTOMETRIC_WEIGHTS to control weighting of pmSourceCrossProduct; provide symbolic names for model parameters

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmSourcePhotometry.c

    r8815 r9527  
    33 *  @author EAM, IfA; GLG, MHPCC
    44 *
    5  *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2006-09-15 09:49:01 $
     5 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2006-10-13 02:24:34 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    9595    }
    9696
    97     if (model->dparams->data.F32[1] > 0) {
    98         SN = model->params->data.F32[1] / model->dparams->data.F32[1];
     97    if (model->dparams->data.F32[PM_PAR_FLUX] > 0) {
     98        SN = model->params->data.F32[PM_PAR_FLUX] / model->dparams->data.F32[PM_PAR_FLUX];
    9999        source->errMag = 1.0 / SN;
    100100    } else {
    … …  
    102102        source->errMag = 0.0;
    103103    }
    104     x = model->params->data.F32[2];
    105     y = model->params->data.F32[3];
     104    x = model->params->data.F32[PM_PAR_XPOS];
     105    y = model->params->data.F32[PM_PAR_YPOS];
    106106
    107107    // measure object model photometry
    … …  
    194194
    195195    if (DO_SKY) {
    196         sky = model->params->data.F32[0];
     196        sky = model->params->data.F32[PM_PAR_SKY];
    197197    } else {
    198198        sky = 0;
    … …  
    233233    // we only care about the value of the object model, not the local sky
    234234    if (DO_SKY) {
    235         sky = model->params->data.F32[0];
     235        sky = model->params->data.F32[PM_PAR_SKY];
    236236    } else {
    237237        sky = 0;
    … …  
    245245    psVector *params = model->params;
    246246
    247     Xo = params->data.F32[2];
    248     Yo = params->data.F32[3];
     247    Xo = params->data.F32[PM_PAR_XPOS];
     248    Yo = params->data.F32[PM_PAR_YPOS];
    249249
    250250    dX = Xo - image->col0;
    … …  
    297297}
    298298
    299 float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj)
     299double pmSourceCrossProduct (const pmSource *Mi,
     300                             const pmSource *Mj,
     301                             const bool unweighted_sum) // should the cross product be weighted?
    300302{
    301303
    … …  
    304306    int xIs, xJs, yIs, yJs;
    305307    int xIe, yIe;
    306     float flux, wt;
    307 
    308     psImage *Pi = Mi->pixels;
    309     psImage *Pj = Mj->pixels;
    310 
    311     psImage *Wi = Mi->weight;
    312 
    313     psImage *Ti = Mi->mask;
    314     psImage *Tj = Mj->mask;
     308    double flux, wt;
     309
     310    const psImage *Pi = Mi->pixels;
     311    const psImage *Pj = Mj->pixels;
     312
     313    const psImage *Wi = Mi->weight;
     314
     315    const psImage *Ti = Mi->mask;
     316    const psImage *Tj = Mj->mask;
    315317
    316318    Xs = PS_MAX (Pi->col0, Pj->col0);
    … …  
    337339            if (Tj->data.U8[yj][xj])
    338340                continue;
    339             wt = Wi->data.F32[yi][xi];
    340             if (wt > 0) {
    341                 flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
     341
     342            if (unweighted_sum) {
     343                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
     344            } else {
     345                wt = Wi->data.F32[yi][xi];
     346                if (wt > 0) {
     347                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
     348                }
    342349            }
    343350        }
    344351    }
    345     return (flux);
    346 }
    347 
    348 float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj)
     352    return flux;
     353}
     354
     355double pmSourceCrossWeight(const pmSource *Mi,
     356                           const pmSource *Mj,
     357                           const bool unweighted_sum) // should the cross product be weighted?
    349358{
    350359
    … …  
    353362    int xIs, xJs, yIs, yJs;
    354363    int xIe, yIe;
    355     float flux, wt;
    356 
    357     psImage *Pi = Mi->pixels;
    358     psImage *Pj = Mj->pixels;
    359 
    360     psImage *Wi = Mi->weight;
    361 
    362     psImage *Ti = Mi->mask;
    363     psImage *Tj = Mj->mask;
     364    double flux, wt;
     365
     366    const psImage *Pi = Mi->pixels;
     367    const psImage *Pj = Mj->pixels;
     368
     369    const psImage *Wi = Mi->weight;
     370
     371    const psImage *Ti = Mi->mask;
     372    const psImage *Tj = Mj->mask;
    364373
    365374    Xs = PS_MAX (Pi->col0, Pj->col0);
    … …  
    386395            if (Tj->data.U8[yj][xj])
    387396                continue;
    388             wt = Wi->data.F32[yi][xi];
    389             if (wt > 0) {
    390                 flux += 1.0 / wt;
     397
     398            if (unweighted_sum) {
     399                flux++;
     400            } else {
     401                wt = Wi->data.F32[yi][xi];
     402                if (wt > 0) {
     403                    flux += 1.0 / wt;
     404                }
    391405            }
    392406        }
    393407    }
    394     return (flux);
    395 }
    396 
     408    return flux;
     409}
     410

  • trunk/psModules/src/objects/pmSourcePhotometry.h

    r6872 r9527  
    33 *  @author EAM, IfA; GLG, MHPCC
    44 *
    5  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    6  *  @date $Date: 2006-04-17 18:01:05 $
     5 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     6 *  @date $Date: 2006-10-13 02:24:34 $
    77 *
    88 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4747bool pmSourceMagnitudesInit (psMetadata *config);
    4848bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode);
    49 float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj);
    50 float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj);
     49double pmSourceCrossProduct(const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
     50double pmSourceCrossWeight(const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
    5151bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *image, psImage *mask);
    5252

  • trunk/psphot/src/psphotEnsemblePSF.c

    r9270 r9527  
    3737    // psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    3838    if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined");
     39
     40    const bool CONSTANT_PHOTOMETRIC_WEIGHTS =
     41        psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
     42    if (!status) {
     43        psAbort(PS_FILE_LINE, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
     44    }
    3945
    4046    for (int i = 0; i < sources->n; i++) {
    … …  
    6874        pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type);
    6975        if (inSource->mode &  PM_SOURCE_MODE_SATSTAR) {
    70             modelEXT->params->data.F32[2] = inSource->moments->x;
    71             modelEXT->params->data.F32[3] = inSource->moments->y;
     76            modelEXT->params->data.F32[PM_PAR_XPOS] = inSource->moments->x;
     77            modelEXT->params->data.F32[PM_PAR_YPOS] = inSource->moments->y;
    7278        } else {
    7379            // peak-up on peak (for non-sat objects)
    … …  
    8187
    8288            psTrace ("psphotEnsemblePSF", 5, "peak coord: %f %f -> %f %f\n",
    83                      modelEXT->params->data.F32[2], modelEXT->params->data.F32[3], min.x + ix, min.y + iy);
     89                     modelEXT->params->data.F32[PM_PAR_XPOS], modelEXT->params->data.F32[PM_PAR_YPOS], min.x + ix, min.y + iy);
    8490           
    8591            // if min point is too deviant, keep the old value
    8692            if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
    87                 modelEXT->params->data.F32[2] = min.x + ix;
    88                 modelEXT->params->data.F32[3] = min.y + iy;
     93                modelEXT->params->data.F32[PM_PAR_XPOS] = min.x + ix;
     94                modelEXT->params->data.F32[PM_PAR_YPOS] = min.y + iy;
    8995            }
    9096            psFree (bicube);
    … …  
    96102
    97103        // save the original coords
    98         x = model->params->data.F32[2];
    99         y = model->params->data.F32[3];
     104        x = model->params->data.F32[PM_PAR_XPOS];
     105        y = model->params->data.F32[PM_PAR_YPOS];
    100106
    101107        // set the fit radius based on the object flux limit and the model
    … …  
    112118
    113119        // set model to unit peak, zero sky (we assume sky is subtracted)
    114         model->params->data.F32[0] = 0.0;
    115         model->params->data.F32[1] = 1.0;
     120        model->params->data.F32[PM_PAR_SKY] = 0.0;
     121        model->params->data.F32[PM_PAR_FLUX] = 1.0;
    116122
    117123        // fill in the model pixel values
    … …  
    143149
    144150        // scale by diagonal element (auto-cross-product)
    145         r = pmSourceCrossProduct (Mi, Mi);
     151        r = pmSourceCrossProduct (Mi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS);
    146152        weight->data.F32[i] = r;
    147153
    … …  
    149155
    150156        // find the image x model value
    151         f = pmSourceCrossProduct (Fi, Mi);
     157        f = pmSourceCrossProduct (Fi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS);
    152158        psSparseVectorElement (sparse, i, f / r);
    153159
    … …  
    163169           
    164170            // got an overlap; calculate cross-product and add to output array
    165             f = pmSourceCrossProduct (Mi, Mj);
     171            f = pmSourceCrossProduct (Mi, Mj, CONSTANT_PHOTOMETRIC_WEIGHTS);
    166172            psSparseMatrixElement (sparse, j, i, f / r);
    167173        }
    … …  
    194200            psAbort ("psphot", "ensemble source is nan");
    195201        }
    196         Fi->modelPSF->params->data.F32[1] = norm->data.F32[i];
    197         Fi->modelPSF->dparams->data.F32[1] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i]));
     202        Fi->modelPSF->params->data.F32[PM_PAR_FLUX] = norm->data.F32[i];
     203        Fi->modelPSF->dparams->data.F32[PM_PAR_FLUX] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i]));
    198204        // XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it
    199205
    … …  
    210216        pmModel *model = Fi->modelPSF;
    211217
    212         x = model->params->data.F32[2];
    213         y = model->params->data.F32[3];
     218        x = model->params->data.F32[PM_PAR_XPOS];
     219        y = model->params->data.F32[PM_PAR_YPOS];
    214220
    215221        psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK);
Note: See TracChangeset for help on using the changeset viewer.