IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14338


Ignore:
Timestamp:
Jul 20, 2007, 10:34:30 AM (19 years ago)
Author:
eugene
Message:

finished PSFConvModel code; added structure to carry convolved model data

Location:
trunk/psphot/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphot.h

    r13983 r14338  
    1616psString        psphotVersionLong(void);
    1717
    18 bool            psphotModelTest (pmReadout *readout, psMetadata *recipe, psMaskType maskVal, psMaskType mark);
     18bool            psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe, psMaskType maskVal, psMaskType mark);
    1919bool            psphotReadout (pmConfig *config, const pmFPAview *view);
    2020bool            psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources);
    … …  
    116116bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal);
    117117bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal);
    118 psKernel *psphotKernelFromPSF (pmSource *source);
     118psKernel *psphotKernelFromPSF (pmSource *source, int nPix);
    119119
    120120bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal);
    … …  
    122122bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal);
    123123bool psphotKron (pmSource *source, psMetadata *recipe, psMaskType maskVal);
     124
     125// structures & functions to support psf-convolved model fitting
     126
     127// pmPCMData : PSF Convolved Model data storage structure
     128typedef struct {
     129    psImage *model;
     130    psArray *dmodels;
     131    psImage *modelConv;
     132    psArray *dmodelsConv;
     133} pmPCMData;
     134
    124135
    125136// psf-convolved model fitting
    … …  
    129140    psVector *params,
    130141    psMinConstraint *constraint,
    131     const pmSource *source,
     142    pmSource *source,
    132143    const psKernel *psf,
    133144    psMinimizeLMChi2Func func);
    … …  
    138149    const psVector *params,
    139150    const psVector *paramMask,
     151    pmPCMData *pcm,
    140152    const pmSource *source,
    141153    const psKernel *psf,
    142154    psMinimizeLMChi2Func func);
    143155
     156pmPCMData *pmPCMDataAlloc (
     157    const psVector *params,
     158    const psVector *paramMask,
     159    pmSource *source);
     160
     161psImage *pmPCMDataSaveImage (pmPCMData *pcm);
    144162
    145163#endif

  • trunk/psphot/src/psphotKernelFromPSF.c

    r13983 r14338  
    11# include "psphot.h"
    22
    3 psKernel *psphotKernelFromPSF (pmSource *source) {
     3psKernel *psphotKernelFromPSF (pmSource *source, int nPix) {
    44
    55  assert (source);
    … …  
    1010
    1111  // need to decide on the size: dynamically? statically?
    12   psKernel *psf = psKernelAlloc (-1, +1, -1, +1);
     12  psKernel *psf = psKernelAlloc (-nPix, +nPix, -nPix, +nPix);
    1313
    1414  for (int j = psf->yMin; j <= psf->yMax; j++) {

  • trunk/psphot/src/psphotModelTest.c

    r14326 r14338  
    11# include "psphotInternal.h"
    2 static char DEFAULT_MODE[] = "EXT";
    3 
    4 // XXX consider this function : add more test information?
    5 bool psphotModelTest (pmReadout *readout, psMetadata *recipe, psMaskType maskVal, psMaskType mark) {
     2# define PM_SOURCE_FIT_PSF_X_EXT PM_SOURCE_FIT_PSF_AND_SKY
     3
     4// XXX add more test information?
     5bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe, psMaskType maskVal, psMaskType mark) {
    66
    77    bool status;
    88    int modelType;
    9     unsigned int Nfail;
    109    float obsMag, fitMag, value;
    1110    char name[64];
    … …  
    1312    pmSourceFitMode fitMode;
    1413
    15     psMetadataItem *item  = NULL;
     14    // run model fitting tests on a single source?
     15    if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false;
     16
     17    // find the currently selected readout
     18    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
     19    PS_ASSERT_PTR_NON_NULL (readout, false);
    1620
    1721    // use poissonian errors or local-sky errors
    … …  
    2024    pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS);
    2125
    22     // run model fitting tests on a single source
    23     if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false;
     26    // find the various fitting parameters (try test values first)
     27    float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS");
     28    if (!status || !isfinite(INNER)) {
     29        INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
     30    }
     31    float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS");
     32    if (!status || !isfinite(OUTER)) {
     33        OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
     34    }
     35    float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS");
     36    if (!status || !isfinite(RADIUS)) {
     37        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
     38    }
     39    float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS");
     40    if (!status || !isfinite(mRADIUS)) {
     41        mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     42    }
     43
     44    // define the source of interest
     45    float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
     46    float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
     47    if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined");
    2448
    2549    // what fitting mode to use?
     50    fitMode = PM_SOURCE_FIT_EXT;
    2651    char *fitModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE");
    27     if (!status || !strcasecmp (fitModeWord, "DEFAULT")) {
    28         fitModeWord = DEFAULT_MODE;
    29     }
    30     fitMode = PM_SOURCE_FIT_EXT;
    31     if (!strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF;
     52    if (fitModeWord && !strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF;
     53    if (fitModeWord && !strcasecmp (fitModeWord, "CONV")) fitMode = PM_SOURCE_FIT_PSF_X_EXT;
     54    if (fitModeWord && !strcasecmp (fitModeWord, "DEFAULT")) fitMode = PM_SOURCE_FIT_EXT;
    3255
    3356    // in fitMode, psf sets the model type
    3457    if (fitMode == PM_SOURCE_FIT_PSF) {
    35         // XXX load psf using psphotLoadPSF
    36         char *psfFile = psMetadataLookupStr (&status, recipe, "PSF_INPUT_FILE");
    37         if (!status) psAbort("PSF_INPUT_FILE not supplied");
    38         psMetadata *psfData = psMetadataConfigRead(NULL, &Nfail, psfFile, FALSE);
    39         psf = pmPSFfromMetadata (psfData);
     58        psf = psphotLoadPSF (config, view, recipe);
     59        if (!psf) psAbort("PSF_INPUT_FILE not supplied");
    4060        modelType = psf->type;
    41     } else {
     61    }
     62    if (fitMode == PM_SOURCE_FIT_EXT) {
    4263        // find the model: supplied by user or first in the PSF_MODEL list
    4364        char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
    … …  
    5778
    5879            // take the first list element
    59             item = psListGet (list, PS_LIST_HEAD);
     80            psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
    6081            modelName = item->data.V;
    6182        }
    … …  
    6384        if (modelType < 0) psAbort("unknown model %s", modelName);
    6485    }
    65 
    66     // find the fitting parameters (try test values first)
    67     float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS");
    68     if (!status || !isfinite(INNER)) {
    69         INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    70     }
    71 
    72     float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS");
    73     if (!status || !isfinite(OUTER)) {
    74         OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    75     }
    76 
    77     float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS");
    78     if (!status || !isfinite(RADIUS)) {
    79         RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
    80     }
    81 
    82     float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS");
    83     if (!status || !isfinite(mRADIUS)) {
    84         mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    85     }
    86 
    87     // define the source of interest
    88     float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
    89     float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
     86    if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
     87        // we need to load BOTH a psf and an ext model
     88        psf = psphotLoadPSF (config, view, recipe);
     89        if (!psf) psAbort("PSF_INPUT_FILE not supplied");
     90
     91        // find the model: supplied by user or first in the PSF_MODEL list
     92        char *modelName  = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL");
     93        if (!status || !strcasecmp (modelName, "DEFAULT")) {
     94            // get the list pointers for the PSF_MODEL entries
     95
     96            psList *list = NULL;
     97            psMetadataItem *mdi = psMetadataLookup (recipe, "PSF_MODEL");
     98            if (mdi == NULL) psAbort("missing PSF_MODEL selection");
     99            if (mdi->type == PS_DATA_STRING) {
     100                list = psListAlloc(NULL);
     101                psListAdd (list, PS_LIST_HEAD, mdi);
     102            } else {
     103                if (mdi->type != PS_DATA_METADATA_MULTI) psAbort("missing PSF_MODEL selection");
     104                list = psMemIncrRefCounter(mdi->data.list);
     105            }
     106
     107            // take the first list element
     108            psMetadataItem *item = psListGet (list, PS_LIST_HEAD);
     109            modelName = item->data.V;
     110        }
     111        modelType = pmModelSetType (modelName);
     112        if (modelType < 0) psAbort("unknown model %s", modelName);
     113    }
    90114
    91115    // construct the source structures
    … …  
    100124    // get the source moments
    101125    status = pmSourceMoments (source, mRADIUS);
    102     if (!status) psAbort("pmSourceLocalSky error");
     126    if (!status) psAbort("psSourceMoments error");
    103127    source->peak->value = source->moments->Peak;
    104128
    … …  
    116140    // get the initial model parameter guess
    117141    pmModel *model = pmSourceModelGuess (source, modelType);
    118     // if any parameters are defined, use those values
     142
     143    // if any parameters are defined by the user, take those values
    119144    int nParams = pmModelParameterCount (modelType);
    120145    psF32 *params = model->params->data.F32;
     146    params[PM_PAR_XPOS] = xObj; // XXX use the user-supplied value,
     147    params[PM_PAR_YPOS] = yObj; // XXX or use the centroid
    121148    for (int i = 0; i < nParams; i++) {
    122         if (i == 2) {
    123             params[i] = xObj;
    124             continue;
    125         }
    126         if (i == 3) {
    127             params[i] = yObj;
    128             continue;
    129         }
     149        if (i == PM_PAR_XPOS) continue;
     150        if (i == PM_PAR_YPOS) continue;
     151
    130152        sprintf (name, "TEST_FIT_PAR%d", i);
    131153        value = psMetadataLookupF32 (&status, recipe, name);
    … …  
    138160    fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
    139161
     162    // for PSF fitting, set the shape parameters based on the PSF & source position
    140163    if (fitMode == PM_SOURCE_FIT_PSF) {
    141164        pmModel *modelPSF = pmModelFromPSF (model, psf);
    … …  
    154177    fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
    155178
    156 
    157179    fprintf (stderr, "input parameters: \n");
    158180    for (int i = 0; i < nParams; i++) {
    … …  
    162184    // define the pixels used for the fit
    163185    psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", mark);
     186    psphotSaveImage (NULL, source->maskObj, "mask1.fits");
    164187
    165188    char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET");
    … …  
    169192    }
    170193
    171     status = pmSourceFitModel (source, model, fitMode, maskVal);
     194    if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
     195        // build the psf for the object
     196        source->modelPSF = pmModelFromPSF (model, psf);
     197        source->modelEXT = model;
     198        status = psphotPSFConvModel (source, recipe, maskVal);
     199        model = source->modelConv;
     200        params = model->params->data.F32;
     201    } else {
     202        status = pmSourceFitModel (source, model, fitMode, maskVal);
     203    }
    172204
    173205    // measure the source mags

  • trunk/psphot/src/psphotModelWithPSF.c

    r14327 r14338  
    66    psVector *params,
    77    psMinConstraint *constraint,
    8     const pmSource *source,
     8    pmSource *source,
    99    const psKernel *psf,
    1010    psMinimizeLMChi2Func func)
    1111{
    12     psTrace("psLib.math", 3, "---- begin ----\n");
     12    psTrace("psphot", 3, "---- begin ----\n");
    1313    PS_ASSERT_PTR_NON_NULL(min, false);
    1414    PS_ASSERT_VECTOR_NON_NULL(params, false);
    … …  
    4646    psF32 dLinear = 0.0;
    4747
     48    // generate PCM data storage structure
     49    pmPCMData *pcm = pmPCMDataAlloc (params, paramMask, source);
     50
    4851    // calculate initial alpha and beta, set chisq (min->value)
    49     min->value = psphotModelWithPSF_SetABX(alpha, beta, params, paramMask, source, psf, func);
     52    min->value = psphotModelWithPSF_SetABX(alpha, beta, params, paramMask, pcm, source, psf, func);
    5053    if (isnan(min->value)) {
    5154        min->iter = min->maxIter;
    … …  
    5356    }
    5457    // dump some useful info if trace is defined
    55     if (psTraceGetLevel("psLib.math") >= 6) {
     58    if (psTraceGetLevel("psphot") >= 6) {
    5659        p_psImagePrint(psTraceGetDestination(), alpha, "alpha guess (0)");
    5760        p_psVectorPrint(psTraceGetDestination(), beta, "beta guess (0)");
    5861    }
    59     if (psTraceGetLevel("psLib.math") >= 5) {
     62    if (psTraceGetLevel("psphot") >= 5) {
    6063        p_psVectorPrint(psTraceGetDestination(), params, "params guess (0)");
    6164    }
    … …  
    6366    // iterate until the tolerance is reached, or give up
    6467    while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) {
    65         psTrace("psLib.math", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
    66         psTrace("psLib.math", 5, "Last delta is %f.  Min->tol is %f.\n", min->lastDelta, min->tol);
     68        psTrace("psphot", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
     69        psTrace("psphot", 5, "Last delta is %f.  Min->tol is %f.\n", min->lastDelta, min->tol);
    6770
    6871
    … …  
    7578
    7679        // dump some useful info if trace is defined
    77         if (psTraceGetLevel("psLib.math") >= 6) {
     80        if (psTraceGetLevel("psphot") >= 6) {
    7881            p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (1)");
    7982            p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (1)");
    8083        }
    81         if (psTraceGetLevel("psLib.math") >= 5) {
     84        if (psTraceGetLevel("psphot") >= 5) {
    8285            p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)");
    8386        }
    8487
    8588        // calculate Chisq for new guess, update Alpha & Beta
    86         Chisq = psphotModelWithPSF_SetABX(Alpha, Beta, Params, paramMask, source, psf, func);
     89        Chisq = psphotModelWithPSF_SetABX(Alpha, Beta, Params, paramMask, pcm, source, psf, func);
    8790        if (isnan(Chisq)) {
    8891            min->iter ++;
    … …  
    97100        psF32 rho = (min->value - Chisq) / dLinear;
    98101
    99         psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value,
     102        psTrace("psphot", 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value,
    100103                Chisq, min->lastDelta, rho);
    101104
    102105        // dump some useful info if trace is defined
    103         if (psTraceGetLevel("psLib.math") >= 6) {
     106        if (psTraceGetLevel("psphot") >= 6) {
    104107            p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (2)");
    105108            p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (2)");
    … …  
    114117            params = psVectorCopy(params, Params, PS_TYPE_F32);
    115118            lambda *= 0.1;
     119
     120            // save the new convolved model image
     121            psFree (source->modelFlux);
     122            source->modelFlux = pmPCMDataSaveImage(pcm);
    116123        } else {
    117124            lambda *= 10.0;
    … …  
    119126        min->iter++;
    120127    }
    121     psTrace("psLib.math", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
     128    psTrace("psphot", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
    122129
    123130    // construct & return the covariance matrix (if requested)
    124131    if (covar != NULL) {
    125132        if (!psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
    126             psTrace ("psLib.math", 5, "failure to calculate covariance matrix\n");
     133            psTrace ("psphot", 5, "failure to calculate covariance matrix\n");
    127134        }
    128135    }
    … …  
    134141    psFree(Beta);
    135142    psFree(Params);
     143    psFree(pcm);
    136144
    137145    if (min->iter == min->maxIter) {
    138         psTrace("psLib.math", 3, "---- end (false) ----\n");
     146        psTrace("psphot", 3, "---- end (false) ----\n");
    139147        return(false);
    140148    }
    141149
    142     psTrace("psLib.math", 3, "---- end (true) ----\n");
     150    psTrace("psphot", 3, "---- end (true) ----\n");
    143151    return(true);
    144152}
    … …  
    149157    const psVector *params,
    150158    const psVector *paramMask,
     159    pmPCMData *pcm,
    151160    const pmSource *source,
    152161    const psKernel *psf,
    … …  
    168177    }
    169178
     179    // generate the model and derivative images for this parameter set
     180
     181    // storage for model derivatives
    170182    psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
    171 
    172     // generate the model and derivative images for this parameter set
    173     psImage *model = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
    174     psArray *dmodels = psArrayAlloc (params->n);
    175     for (psS32 n = 0; n < params->n; n++) {
    176       if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
    177       psImage *dmodel = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
    178       dmodels->data[n] = dmodel;
    179     }
    180183
    181184    // working vector to store local coordinate
    … …  
    206209            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    207210
    208             model->data.F32[i][j] = func (deriv, params, coord);
     211            pcm->model->data.F32[i][j] = func (deriv, params, coord);
    209212
    210213            for (int n = 0; n < params->n; n++) {
    211214              if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
    212               psImage *dmodel = dmodels->data[n];
     215              psImage *dmodel = pcm->dmodels->data[n];
    213216              dmodel->data.F32[i][j] = deriv->data.F32[n];
    214217            }
    215218        }
    216219    }
    217     psFree (coord);
     220    psFree(coord);
     221    psFree(deriv);
     222
     223    psphotSaveImage (NULL, pcm->model, "model1.fits");
    218224
    219225    // convolve model and dmodel arrays with PSF
    220     psImage *modelConv = psImageConvolveDirect (model, psf);
    221     psArray *dmodelsConv = psArrayAlloc (params->n);
    222     for (int n = 0; n < params->n; n++) {
    223       if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
    224       psImage *dmodel = dmodels->data[n];
    225       psImage *dmodelConv = psImageConvolveDirect (dmodel, psf);
    226       dmodelsConv->data[n] = dmodelConv;
    227     }
     226    psImageConvolveDirect (pcm->modelConv, pcm->model, psf);
     227    for (int n = 0; n < pcm->dmodels->n; n++) {
     228      if (pcm->dmodels->data[n] == NULL) continue;
     229      psImage *dmodel = pcm->dmodels->data[n];
     230      psImage *dmodelConv = pcm->dmodelsConv->data[n];
     231      psImageConvolveDirect (dmodelConv, dmodel, psf);
     232    }
     233
     234    // XXX TEST : SAVE IMAGES
     235    psphotSaveImage (NULL, pcm->model, "model.fits");
     236    psphotSaveImage (NULL, pcm->modelConv, "modelConv.fits");
     237    psphotSaveImage (NULL, source->pixels, "obj.fits");
     238    psphotSaveImage (NULL, source->maskObj, "mask.fits");
     239    psphotSaveImage (NULL, source->weight, "weight.fits");
     240    exit (0);
    228241
    229242    // zero alpha and beta for summing below
    … …  
    248261            }
    249262
    250             float ymodel  = modelConv->data.F32[i][j];
     263            float ymodel  = pcm->modelConv->data.F32[i][j];
    251264            float yweight = 1.0 / source->weight->data.F32[i][j];
    252265            float delta = ymodel - source->pixels->data.F32[i][j];
    … …  
    263276                continue;
    264277              }
    265               psImage *dmodel = dmodelsConv->data[n1];
     278              psImage *dmodel = pcm->dmodelsConv->data[n1];
    266279              float weight = dmodel->data.F32[i][j] * yweight;
    267280              for (psS32 n2 = 0; n2 <= n1; n2++) {
    … …  
    269282                  continue;
    270283                }
    271                 dmodel = dmodelsConv->data[n2];
     284                dmodel = pcm->dmodelsConv->data[n2];
    272285                alpha->data.F32[n1][n2] += weight * dmodel->data.F32[i][j];
    273286              }
    … …  
    294307    }
    295308
    296     psFree (model);
    297     psFree (dmodels);
    298     psFree (modelConv);
    299     psFree (dmodelsConv);
    300     psFree(deriv);
    301 
    302309    return(chisq);
    303310}
    304311
    305 
     312static void pmPCMDataFree (pmPCMData *pcm) {
     313
     314    if (pcm == NULL) return;
     315
     316    psFree (pcm->model);
     317    psFree (pcm->modelConv);
     318    psFree (pcm->dmodels);
     319    psFree (pcm->dmodelsConv);
     320    return;
     321}
     322
     323pmPCMData *pmPCMDataAlloc (
     324    const psVector *params,
     325    const psVector *paramMask,
     326    pmSource *source) {
     327
     328    pmPCMData *pcm = (pmPCMData *) psAlloc(sizeof(pmPCMData));
     329    psMemSetDeallocator(pcm, (psFreeFunc) pmPCMDataFree);
     330
     331    // Allocate storage images for raw model and derivative images
     332    pcm->model = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
     333    pcm->dmodels = psArrayAlloc (params->n);
     334    for (psS32 n = 0; n < params->n; n++) {
     335      if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
     336      pcm->dmodels->data[n] = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
     337    }
     338
     339    // Allocate storage images for convolved model and derivative images
     340    pcm->modelConv = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
     341    pcm->dmodelsConv = psArrayAlloc (params->n);
     342    for (psS32 n = 0; n < params->n; n++) {
     343      if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
     344      pcm->dmodelsConv->data[n] = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
     345    }
     346
     347    return pcm;
     348}
     349
     350psImage *pmPCMDataSaveImage (pmPCMData *pcm) {
     351
     352    psImage *model = psImageCopy (NULL, pcm->model, PS_TYPE_F32);
     353    return model;
     354}
    306355
    307356/*

  • trunk/psphot/src/psphotPSFConvModel.c

    r13983 r14338  
    77// static bool  PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;
    88
     9// input source has both modelPSF and modelEXT.  on successful exit, we set the
     10// modelConv to contain the fitted parameters, and the modelFlux to contain the
     11// convolved model image.
    912bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
    1013   
    11     // XXX make sure we save a cached copy of the source psf
     14    // make sure we save a cached copy of the psf flux
     15    pmSourceCachePSF (source, maskVal);
    1216
    1317    // convert the cached cached psf model for this source to a psKernel
    14     psKernel *psf = psphotKernelFromPSF (source);
     18    // XXX for the moment, hard-wire the kernel to be 5x5 (2 pix radius)
     19    psKernel *psf = psphotKernelFromPSF (source, 2);
    1520
    1621    // generate copy of the model
    … …  
    4853    }
    4954
     55    // set up the minimization process
    5056    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
    5157
    … …  
    5460    bool fitStatus = psphotModelWithPSF_LMM (myMin, covar, params, constraint, source, psf, modelFunc);
    5561    for (int i = 0; i < dparams->n; i++) {
    56         if (psTraceGetLevel("psModules.objects") >= 4) {
     62        if (psTraceGetLevel("psphot") >= 4) {
    5763            fprintf (stderr, "%f ", params->data.F32[i]);
    5864        }
    … …  
    6167        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
    6268    }
    63     psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
     69    psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
    6470
    6571    // save the resulting chisq, nDOF, nIter
    … …  
    9197
    9298    bool retval = (onPic && fitStatus);
    93     psTrace("psModules.objects", 5, "---- %s(%d) end ----\n", __func__, retval);
     99    psTrace("psphot", 5, "---- %s(%d) end ----\n", __func__, retval);
    94100    return(retval);
    95101}
Note: See TracChangeset for help on using the changeset viewer.