IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 20, 2010, 7:21:53 PM (16 years ago)
Author:
eugene
Message:

update PCM fitting (not quite ready)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitPCM.c

    r28687 r28692  
    1 # include "psphotInternal.h"
    2 # define USE_DELTA_PSF 0
     1/* @file  pmSourceFitPCM.c
     2 * structures and functions to support PSF-convolved model fitting
     3 *
     4 * @author EAM, IfA
     5 *
     6 * @version $Revision: 1.29 $
     7 * @date $Date: 2009-02-16 22:30:50 $
     8 * Copyright 2010 Institute for Astronomy, University of Hawaii
     9 */
    310
    4 // save as static values so they may be set externally
    5 static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
    6 static psF32 PM_SOURCE_FIT_MODEL_MIN_TOL = 0.1;
    7 static psF32 PM_SOURCE_FIT_MODEL_MAX_TOL = 2.0;
     11#ifdef HAVE_CONFIG_H
     12#include <config.h>
     13#endif
     14
     15#include <stdio.h>
     16#include <math.h>
     17#include <string.h>
     18#include <strings.h>
     19#include <pslib.h>
     20#include "pmHDU.h"
     21#include "pmFPA.h"
     22#include "pmFPAMaskWeight.h"
     23
     24#include "pmTrend2D.h"
     25#include "pmResiduals.h"
     26#include "pmGrowthCurve.h"
     27#include "pmSpan.h"
     28#include "pmFootprintSpans.h"
     29#include "pmFootprint.h"
     30#include "pmPeaks.h"
     31#include "pmMoments.h"
     32#include "pmModelFuncs.h"
     33#include "pmModel.h"
     34#include "pmModelUtils.h"
     35#include "pmModelClass.h"
     36#include "pmSourceMasks.h"
     37#include "pmSourceExtendedPars.h"
     38#include "pmSourceDiffStats.h"
     39#include "pmSource.h"
     40#include "pmSourceFitModel.h"
     41#include "pmPCMdata.h"
     42
     43# define FACILITY "psModules.objects"
    844
    945// input source has both modelPSF and modelEXT.  on successful exit, we set the
    … …  
    1147// convolved model image.
    1248
    13 // XXX need to generalize this -- number of fitted parameters must be flexible based on the fitOptions
    14 
    15 pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
     49bool pmSourceFitPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
    1650   
    17     // maskVal is used to test for rejected pixels, and must include markVal
    18     maskVal |= markVal;
    19 
    20     // make sure we save a cached copy of the psf flux
    21     pmSourceCachePSF (source, maskVal);
    22 
    23     // convert the cached cached psf model for this source to a psKernel
    24     psKernel *psf = psphotKernelFromPSF (source, psfSize);
    25     if (!psf) return NULL;
    26 
    27 # if (USE_DELTA_PSF)
    28     psImageInit (psf->image, 0.0);
    29     psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0;
    30 # endif
    31 
    32     // generate copy of the model
    33     // XXX we could modify the parameter values or even the model
    34     // here based on the observed seeing (some lookup table...)
    35 
    36     // use the source moments, etc to guess basic model parameters
    37     pmModel *modelConv = pmSourceModelGuess (source, modelType);
    38     if (!modelConv) {
    39         psFree (psf);
    40         return NULL;
    41     }
    42 
    43     // adjust the pixels based on the footprint
    44     float radius = psphotSetRadiusEXT (readout, source, markVal);
    45     if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false;
    46 
    47     // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
    48     psEllipseShape psfShape;
    49     psfShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
    50     psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
    51     psfShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
    52     psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0);
    53 
    54     // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
    55     psEllipseShape extShape;
    56     extShape.sx  = modelConv->params->data.F32[PM_PAR_SXX] / M_SQRT2;
    57     extShape.sxy = modelConv->params->data.F32[PM_PAR_SXY];
    58     extShape.sy  = modelConv->params->data.F32[PM_PAR_SYY] / M_SQRT2;
    59     psEllipseAxes extAxes = psEllipseShapeToAxes (extShape, 20.0);
    60 
    61     // decrease the initial guess ellipse by psf_minor axis:
    62     psEllipseAxes extAxesMod;
    63     extAxesMod.major = sqrt (PS_MAX (1.0, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor)));
    64     extAxesMod.minor = sqrt (PS_MAX (1.0, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor)));
    65     extAxesMod.theta = extAxes.theta;
    66 
    67     psEllipseShape extShapeMod = psEllipseAxesToShape (extAxesMod);
    68     modelConv->params->data.F32[PM_PAR_SXX] = extShapeMod.sx * M_SQRT2;
    69     modelConv->params->data.F32[PM_PAR_SXY] = extShapeMod.sxy;
    70     modelConv->params->data.F32[PM_PAR_SYY] = extShapeMod.sy * M_SQRT2;
    71 
    72     // increase the initial guess central intensity by 2pi r^2:
    73     modelConv->params->data.F32[PM_PAR_I0] *= (1.0 + PS_SQR(psfAxes.minor) / PS_SQR(extAxesMod.minor));
    74 
    75     psVector *params  = modelConv->params;
    76     psVector *dparams = modelConv->dparams;
    77 
    78     // create the minimization constraints
    79     psMinConstraint *constraint = psMinConstraintAlloc();
    80     constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK);
    81     constraint->checkLimits = modelConv->modelLimits;
    82 
    83     // set parameter mask based on fitting mode
    84     // we fit a model without a floating sky term
    85     int nParams = params->n - 1;
    86     psVectorInit (constraint->paramMask, 0);
    87     constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
     51    psVector *params  = pcm->modelConv->params;
     52    psVector *dparams  = pcm->modelConv->dparams;
    8853
    8954    // force the floating parameters to fall within the contraint ranges
    9055    for (int i = 0; i < params->n; i++) {
    91         modelConv->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
    92         modelConv->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
     56        pcm->modelConv->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
     57        pcm->modelConv->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
    9358    }
    9459
    9560    // set up the minimization process
    96     psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_MIN_TOL, PM_SOURCE_FIT_MODEL_MAX_TOL);
     61    psMinimization *myMin = psMinimizationAlloc (fitOptions->nIter, fitOptions->minTol, fitOptions->maxTol);
    9762
    9863    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    9964
    100     bool fitStatus = psphotModelWithPSF_LMM (myMin, covar, params, constraint, source, psf, modelConv->modelFunc);
     65    bool fitStatus = pmPCM_MinimizeChisq (myMin, covar, params, source, pcm);
    10166    for (int i = 0; i < dparams->n; i++) {
    102         if (psTraceGetLevel("psphot") >= 4) {
    103             fprintf (stderr, "%f ", params->data.F32[i]);
    104         }
    105         if ((constraint->paramMask != NULL) && constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i])
     67        if ((pcm->constraint->paramMask != NULL) && pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i])
    10668            continue;
    10769        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
     70        psTrace ("psModules.objects", 4, "%f +/- %f", params->data.F32[i], dparams->data.F32[i]);
    10871    }
    10972    psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
    … …  
    11881
    11982    // save the resulting chisq, nDOF, nIter
    120     modelConv->chisq = myMin->value;
    121     modelConv->nIter = myMin->iter;
    122 
    123     // XXX I actually need to count the number of unmasked pixels here
    124     modelConv->nDOF  = source->pixels->numCols*source->pixels->numRows  -  nParams;
    125 
    126     modelConv->flags |= PM_MODEL_STATUS_FITTED;
    127     if (!fitStatus) modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE;
     83    pcm->modelConv->chisq = myMin->value;
     84    pcm->modelConv->nIter = myMin->iter;
     85    pcm->modelConv->nPix = pcm->nPix;
     86    pcm->modelConv->nDOF = pcm->nDOF;
     87    pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF;
     88    pcm->modelConv->flags |= PM_MODEL_STATUS_FITTED;
     89    if (!fitStatus) pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE;
    12890
    12991    // models can go insane: reject these
    … …  
    13395    onPic &= (params->data.F32[PM_PAR_YPOS] >= source->pixels->row0);
    13496    onPic &= (params->data.F32[PM_PAR_YPOS] <  source->pixels->row0 + source->pixels->numRows);
    135     if (!onPic) modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE;
     97    if (!onPic) pcm->modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE;
    13698
    13799    source->mode |= PM_SOURCE_MODE_FITTED; // XXX is this needed?
    138100
    139     psFree(psf);
    140101    psFree(myMin);
    141102    psFree(covar);
    142     psFree(constraint);
    143103
    144     return modelConv;
     104    return true;
    145105}
     106
     107bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) {
     108
     109    pcm->modelConv->modelGuess(pcm->modelConv, source);
     110
     111    // generate copy of the model
     112    // XXX we could modify the parameter values or even the model
     113    // here based on the observed seeing (some lookup table...)
     114
     115    // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
     116    psEllipseShape psfShape;
     117    psfShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
     118    psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
     119    psfShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
     120    psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0);
     121
     122    // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
     123    psEllipseShape extShape;
     124    extShape.sx  = pcm->modelConv->params->data.F32[PM_PAR_SXX] / M_SQRT2;
     125    extShape.sxy = pcm->modelConv->params->data.F32[PM_PAR_SXY];
     126    extShape.sy  = pcm->modelConv->params->data.F32[PM_PAR_SYY] / M_SQRT2;
     127    psEllipseAxes extAxes = psEllipseShapeToAxes (extShape, 20.0);
     128
     129    // decrease the initial guess ellipse by psf_minor axis:
     130    psEllipseAxes extAxesMod;
     131    extAxesMod.major = sqrt (PS_MAX (1.0, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor)));
     132    extAxesMod.minor = sqrt (PS_MAX (1.0, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor)));
     133    extAxesMod.theta = extAxes.theta;
     134
     135    psEllipseShape extShapeMod = psEllipseAxesToShape (extAxesMod);
     136    pcm->modelConv->params->data.F32[PM_PAR_SXX] = extShapeMod.sx * M_SQRT2;
     137    pcm->modelConv->params->data.F32[PM_PAR_SXY] = extShapeMod.sxy;
     138    pcm->modelConv->params->data.F32[PM_PAR_SYY] = extShapeMod.sy * M_SQRT2;
     139
     140    // increase the initial guess central intensity by 2pi r^2:
     141    pcm->modelConv->params->data.F32[PM_PAR_I0] *= (1.0 + PS_SQR(psfAxes.minor) / PS_SQR(extAxesMod.minor));
     142
     143    return true;
     144}
Note: See TracChangeset for help on using the changeset viewer.