IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36718 for trunk/psphot/src


Ignore:
Timestamp:
May 1, 2014, 4:36:26 PM (12 years ago)
Author:
bills
Message:

Modify psphotFullForce and psphotFullForceSummary to support forcing
multiple extended model fits per source.

Location:
trunk/psphot/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphot.h

    r36623 r36718  
    537537bool psphotGalaxyShape_Threaded (psThreadJob *job);
    538538bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize);
    539 bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize, bool saveResults);
     539bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, pmSourceGalaxyFits *galaxyFits, psImageMaskType maskVal, int psfSize, bool saveResults);
    540540psphotGalaxyShapeOptions *psphotGalaxyShapeOptionsAlloc();
    541541

  • trunk/psphot/src/psphotFullForceSummary.c

    r36666 r36718  
    187187
    188188    // Convert to an array indexed by ID
    189     psArray *sortedTable = psArrayAlloc(4*inTable->n);
     189#define MAX_MODEL_TYPE 10
     190    psArray *sortedTables = psArrayAlloc(MAX_MODEL_TYPE+1);
    190191    for (int i=0; i<inTable->n; i++) {
    191192        psMetadata *row = inTable->data[i];
    192193        psS32 ID = psMetadataLookupS32(&status, row, "ID");
     194        psS32 modelType = psMetadataLookupS32(&status, row, "MODEL_TYPE");
     195        // XXX: need to use the lookup table functions to be ready for changes in the model type numbers
     196        if (modelType+1 >= MAX_MODEL_TYPE) {
     197            psError(PSPHOT_ERR_CONFIG, false, "found modelType %d max allowed is %d", modelType, MAX_MODEL_TYPE);
     198            return false;
     199        }
     200        psArray *sortedTable = sortedTables->data[modelType+1];
     201        if (!sortedTable) {
     202            sortedTable = psArrayAlloc(4*inTable->n);
     203            sortedTables->data[modelType+1] = sortedTable;
     204            // dont' free sortedTable the array of tables gets our reference
     205        }
    193206        if (ID >= sortedTable->n) {
    194207            sortedTable = psArrayRealloc(sortedTable, 2*ID);
    195208            // Why doesn't psArrayRealloc do this?????
    196             sortedTable->n = sortedTable->nalloc;
     209                sortedTable->n = sortedTable->nalloc;
    197210        }
    198211        if (sortedTable->data[ID]) {
    … …  
    202215        sortedTable->data[ID] = psMemIncrRefCounter(row);
    203216    }
    204     psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLE", PS_META_REPLACE, "cff table", sortedTable);
     217    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLES", PS_META_REPLACE, "cff tables", sortedTables);
    205218    psFree(inTable);
    206     psFree(sortedTable);
     219    psFree(sortedTables);
    207220
    208221    return true;

  • trunk/psphot/src/psphotFullForceSummaryReadout.c

    r36701 r36718  
    1414    psArray *zeroPt;
    1515    psArray *exptime;
    16     psArray *cffTable;
     16    psArray *cffTables;
    1717} galaxyShapeOptions;
    1818
    … …  
    7171                return false;
    7272            }
    73             options.cffTable = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLE");
    74             if (!options.cffTable) {
     73            options.cffTables = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLES");
     74            if (!options.cffTables) {
    7575                psError (PS_ERR_UNKNOWN, true, "Cannot find cff table in arguments.");
    7676                return false;
    … …  
    178178    pmSource *outSrc = NULL;
    179179
    180     psVector *sumWeightedFlux = NULL;
    181     psVector *sumInvSig2 = NULL;
    182     psVector *numerator   = NULL;
    183     psF32   totalNPix = 0;
    184     long    vectorLength = 0;
     180
     181    pmSource *source0 = obj->sources->data[0];
     182    int ID = source0->seq;
     183
     184    // Find the parameters from the cff tables
     185    int nModels = 0;
     186    // skip table 0 which is modelType -1 not extended
     187    for (int i=1; i<options->cffTables->n; i++) {
     188        psArray *table = options->cffTables->data[i];
     189        if (!table) continue;
     190        if (ID >= table->n) continue;
     191        if (table->data[ID]) {
     192            nModels++;
     193        }
     194    }
    185195
    186196    // Loop over sources for this object. Start output source for first input that has
    … …  
    210220                outSrc->extpars =  psMemIncrRefCounter(source->extpars);
    211221            }
    212             if (source->modelEXT) {
    213                 outSrc->modelEXT =  psMemIncrRefCounter(source->modelEXT);
    214             }
    215             if (source->modelFits) {
    216                 outSrc->modelFits =  psMemIncrRefCounter(source->modelFits);
    217             }
    218         }
    219 
    220         // The only parameters that we are summarizing currently are galaxy shapes.
    221         // Continue for inputs that do not have modelEXT which contains the model paramters
    222         if (!source->modelEXT) continue;
    223 
    224         if (source->galaxyFits && isfinite(source->galaxyFits->nPix) && source->galaxyFits->chisq->n) {
    225             if (numerator == NULL) {
    226                 // first source with galaxyFits allocate accumulators
    227                 vectorLength = source->galaxyFits->chisq->n;
    228                 sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32);
    229                 psVectorInit(sumWeightedFlux, 0.0);
    230                 sumInvSig2 = psVectorAlloc(vectorLength, PS_TYPE_F32);
    231                 psVectorInit(sumInvSig2, 0.0);
    232                 numerator   = psVectorAlloc(vectorLength, PS_TYPE_F32);
    233                 psVectorInit(numerator, 0.0);
    234             }
    235 
    236             // Die if the lengths of the vectors is not the same in all sources
    237             psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d",
    238                 vectorLength, options->numTrials);
    239             psAssert(source->galaxyFits->chisq->n == vectorLength, "length of chisq vectors do not match %ld %ld",
    240                     source->galaxyFits->chisq->n, vectorLength);
    241 
    242             psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID];
    243 
    244             totalNPix += source->galaxyFits->nPix;
    245 
    246             for (int k = 0; k < vectorLength; k++) {
    247 
    248                 psF32 chisq = source->galaxyFits->chisq->data.F32[k];
    249                 psF32 flux  = source->galaxyFits->Flux->data.F32[k]  * scaleFactor;
    250                 psF32 dFlux = source->galaxyFits->dFlux->data.F32[k] * scaleFactor;
    251 
    252                 numerator->data.F32[k] += chisq * source->galaxyFits->nPix;
    253 
    254                 psF32 invSig2 = 1.0 / (dFlux * dFlux);
    255                 sumInvSig2->data.F32[k] += invSig2;
    256 
    257                 sumWeightedFlux->data.F32[k] += flux * invSig2;
    258             }
    259         }
    260     }
    261 
    262     if (outSrc) {
     222            outSrc->modelEXT = psFree(outSrc->modelEXT);
     223        }
     224        if (nModels && !outSrc->modelFits && source->modelFits && source->modelFits->n == nModels) {
     225            outSrc->modelFits =  psMemIncrRefCounter(source->modelFits);
     226        }
     227    }
     228
     229    if (!outSrc) {
     230        // no good measurements for this source (how?)
     231        return NULL;
     232    }
     233
     234    // now loop over the model fits and summarize the galaxy shape measurements for each
     235    for (int iModel=0; outSrc->modelFits && iModel < outSrc->modelFits->n; iModel++) {
     236        pmModel *model = outSrc->modelFits->data[iModel];
     237        int modelType = model->type;
     238        psAssert(modelType >= 0 && modelType+1 < options->cffTables->n, "model type out of range");
     239
     240        psArray *cffTable = options->cffTables->data[modelType+1];
     241        psAssert(cffTable, "missing cff table?");
     242
     243        psMetadata *cffRow = cffTable->data[outSrc->seq];
     244        if (!cffRow) {
     245            psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d model %d", outSrc->seq, modelType);
     246            psFree(outSrc);
     247            return NULL;
     248        }
     249        bool mdok;
     250        psF32 R_MAJOR = psMetadataLookupF32(&mdok, cffRow, "R_MAJOR");
     251        if (!mdok) {
     252            psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d type %d", outSrc->seq, modelType);
     253            psFree(outSrc);
     254            return NULL;
     255        }
     256        psF32 R_MINOR = psMetadataLookupF32(&mdok, cffRow, "R_MINOR");
     257        if (!mdok) {
     258            psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d type %d", outSrc->seq, modelType);
     259            psFree(outSrc);
     260            return NULL;
     261        }
     262
     263        psVector *sumWeightedFlux = NULL;
     264        psVector *sumInvSig2 = NULL;
     265        psVector *numerator   = NULL;
     266        psF32   totalNPix = 0;
     267        long    vectorLength = 0;
     268        for (int i=0; i < obj->sources->n; i++) {
     269            pmSource *source = obj->sources->data[i];
     270
     271            if (source->pixWeightNotPoor < .9) continue;
     272
     273            // The only parameters that we are summarizing currently are galaxy shapes.
     274            if (!source->modelFits) continue;
     275            if (!source->galaxyFits) continue;
     276
     277            // XXX: put this into a function
     278            pmSourceGalaxyFits *galaxyFits = NULL;
     279            for (int j=0; j<source->galaxyFits->n; j++) {
     280                galaxyFits = source->galaxyFits->data[j];
     281                if (galaxyFits && galaxyFits->modelType == modelType) break;
     282                galaxyFits = NULL;
     283            }
     284            // end of proposed function
     285
     286            if (galaxyFits && isfinite(galaxyFits->nPix) && galaxyFits->chisq->n) {
     287                if (numerator == NULL) {
     288                    // first source with galaxyFits allocate accumulators
     289                    vectorLength = galaxyFits->chisq->n;
     290                    sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32);
     291                    psVectorInit(sumWeightedFlux, 0.0);
     292                    sumInvSig2 = psVectorAlloc(vectorLength, PS_TYPE_F32);
     293                    psVectorInit(sumInvSig2, 0.0);
     294                    numerator   = psVectorAlloc(vectorLength, PS_TYPE_F32);
     295                    psVectorInit(numerator, 0.0);
     296                }
     297
     298                // Die if the lengths of the vectors is not the same in all sources
     299                psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d",
     300                    vectorLength, options->numTrials);
     301                psAssert(galaxyFits->chisq->n == vectorLength, "length of chisq vectors do not match %ld %ld",
     302                             galaxyFits->chisq->n, vectorLength);
     303
     304                psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID];
     305
     306                totalNPix += galaxyFits->nPix;
     307
     308                for (int k = 0; k < vectorLength; k++) {
     309                    psF32 chisq = galaxyFits->chisq->data.F32[k];
     310                    psF32 flux  = galaxyFits->Flux->data.F32[k]  * scaleFactor;
     311                    psF32 dFlux = galaxyFits->dFlux->data.F32[k] * scaleFactor;
     312
     313                    numerator->data.F32[k] += chisq * galaxyFits->nPix;
     314
     315                    psF32 invSig2 = 1.0 / (dFlux * dFlux);
     316                    sumInvSig2->data.F32[k] += invSig2;
     317
     318                    sumWeightedFlux->data.F32[k] += flux * invSig2;
     319                }
     320            }
     321        }
     322
    263323        if (vectorLength) {
    264324            // allocate galaxyFits for the output source
    265             outSrc->galaxyFits = pmSourceGalaxyFitsAlloc();
    266             outSrc->galaxyFits->nPix  = totalNPix;
     325            if (!outSrc->galaxyFits) {
     326                outSrc->galaxyFits = psArrayAllocEmpty(1);
     327            }
     328            pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc();
     329            psArrayAdd(outSrc->galaxyFits, 1, galaxyFits);
     330            psFree(galaxyFits);
     331
     332            galaxyFits->nPix  = totalNPix;
     333            galaxyFits->modelType = modelType;
    267334            psVector *fluxVec =
    268                 outSrc->galaxyFits->Flux  = psVectorRecycle(outSrc->galaxyFits->Flux,  vectorLength, PS_TYPE_F32);
     335                galaxyFits->Flux  = psVectorRecycle(galaxyFits->Flux,  vectorLength, PS_TYPE_F32);
    269336            psVector *dFluxVec =
    270                 outSrc->galaxyFits->dFlux = psVectorRecycle(outSrc->galaxyFits->dFlux, vectorLength, PS_TYPE_F32);
     337                galaxyFits->dFlux = psVectorRecycle(galaxyFits->dFlux, vectorLength, PS_TYPE_F32);
    271338            psVector *chisqVec =
    272                 outSrc->galaxyFits->chisq = psVectorRecycle(outSrc->galaxyFits->chisq, vectorLength, PS_TYPE_F32);
    273 
    274             // Find the nominal radii paramters in the cff table
    275             if (outSrc->seq >= options->cffTable->n) {
    276                 psError(PS_ERR_PROGRAMMING, true, "object sequence number %d is larger than cff table length %ld",
    277                     outSrc->seq, options->cffTable->n);
    278                 return NULL;
    279             }
    280             psMetadata *cffRow = options->cffTable->data[outSrc->seq];
    281             if (!cffRow) {
    282                 psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d", outSrc->seq);
    283                 psFree(outSrc);
    284                 return NULL;
    285             }
    286             bool mdok;
    287             psF32 R_MAJOR = psMetadataLookupF32(&mdok, cffRow, "R_MAJOR");
    288             if (!mdok) {
    289                 psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d", outSrc->seq);
    290                 psFree(outSrc);
    291                 return NULL;
    292             }
    293             psF32 R_MINOR = psMetadataLookupF32(&mdok, cffRow, "R_MINOR");
    294             if (!mdok) {
    295                 psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d", outSrc->seq);
    296                 psFree(outSrc);
    297                 return NULL;
    298             }
     339                galaxyFits->chisq = psVectorRecycle(galaxyFits->chisq, vectorLength, PS_TYPE_F32);
     340
    299341
    300342            // fill the summary galaxyFits vectors and find the trial with the minimum value for chisq
    … …  
    319361            psFree(sumWeightedFlux);
    320362
    321             if (min_k >= 0 && isfinite(minChisq) && outSrc->modelEXT) {
     363            if (min_k >= 0 && isfinite(minChisq)) {
    322364                // copy the best fit params to the model
    323365                // fractional radii with the lowest chisq
    324                 psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type);
     366                psEllipseAxes axes = pmPSF_ModelToAxes(model->params->data.F32, model->type);
    325367
    326368                // examine the params for the trial with minimum chisq.
    … …  
    472514                        axes.major = major0;
    473515                        axes.minor = minor0;
    474                         outSrc->modelEXT->chisq = chisq0;
     516                        model->chisq = chisq0;
    475517#ifdef DUMPVECTORS
    476518                        char fn[80];
    … …  
    487529#endif
    488530                    } else {
    489                         outSrc->modelEXT->flags |= PM_MODEL_STATUS_NONCONVERGE;
     531                        model->flags |= PM_MODEL_STATUS_NONCONVERGE;
    490532
    491533#ifdef PRINTVALS
    … …  
    496538                        // psFree(outSrc->galaxyFits);
    497539
    498                         // XXX drop the model since it is bogus? Perhaps we do this or use the nominal values
    499                         // and add a flag
    500 #ifdef notnow
    501                         psFree(outSrc->modelEXT);
    502                         outSrc->modelEXT = 0;
    503                         psFree(outSrc->modelFits);
    504                         outSrc->modelFits = 0;
    505 #endif
    506540                    }
    507541                    psFree(chisqFit);
    … …  
    516550                    axes.major = R_MAJOR * fRMajorBest;
    517551                    axes.minor = R_MINOR * fRMinorBest;
    518                     outSrc->modelEXT->chisq = minChisq;
     552                    model->chisq = minChisq;
    519553                    flux0 = fluxBest;
    520554#ifdef PRINTVALS
    … …  
    525559                // now save the model parameters in the model structure
    526560
    527                 pmPSF_AxesToModel (outSrc->modelEXT->params->data.F32, axes, outSrc->modelEXT->type);
    528 
    529                 outSrc->modelEXT->mag = -2.5 * log10(flux0);
    530                 outSrc->modelEXT->magErr = dFlux0 / flux0; // 1 / SN
     561                pmPSF_AxesToModel (model->params->data.F32, axes, modelType);
     562
     563                model->mag = -2.5 * log10(flux0);
     564                model->magErr = dFlux0 / flux0; // 1 / SN
    531565                // XXX: should there be a different flag if we didn't do a fit of the chisq
    532                 outSrc->modelEXT->flags |= PM_MODEL_STATUS_FITTED;
     566                model->flags |= PM_MODEL_STATUS_FITTED;
    533567            }
    534568        }

  • trunk/psphot/src/psphotGalaxyShape.c

    r36633 r36718  
    238238bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize) {
    239239
    240     // we use the 0th model to define the initial guess shape
    241     pmModel *model = source->modelFits->data[0];
    242     if (!model) return false;
    243 
    244     pmModelType modelType = model->type;
    245 
    246     // we are using fitOptions->mode : be sure this makes sense
    247     pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
    248     if (!pcm) return false;
    249 
    250     // we are fitting only PM_PAR_I0; the shape elements are generated from a grid
    251     psF32 *PAR = pcm->modelConv->params->data.F32;
    252 
    253     // double check that the guess is carried along...
    254 
    255     // I have some source guess (e0, e1, e2)
    256     psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType);
    257 
    258     if (!source->galaxyFits) {
    259         source->galaxyFits = pmSourceGalaxyFitsAlloc();
    260     }
    261 
    262     float fRmajorBest = NAN;
    263     float fRminorBest = NAN;
    264     float chisqBest = NAN;
    265     for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; fRmajor += opt->fRmajorDel) {
    266         for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel; fRminor += opt->fRminorDel) {
    267  
    268             psEllipseAxes testAxes = guessAxes;
    269             testAxes.major = guessAxes.major * fRmajor;
    270             testAxes.minor = guessAxes.minor * fRminor;
    271            
    272             pmPSF_AxesToModel (PAR, testAxes, modelType);
    273            
    274             psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, true);
    275 
    276             int i = source->galaxyFits->chisq->n - 1;
    277             float flux = source->galaxyFits->Flux->data.F32[i];
    278             if (isfinite(flux)) {
    279                 float thisChisq = source->galaxyFits->chisq->data.F32[i];
    280                 if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
    281                     chisqBest = thisChisq;
    282                     fRmajorBest = fRmajor;
    283                     fRminorBest = fRminor;
     240    for (int iModel = 0 ; iModel < source->modelFits->n; iModel++) {
     241        pmModel *model = source->modelFits->data[iModel];
     242        if (!model) return false;
     243
     244        pmModelType modelType = model->type;
     245
     246        // we are using fitOptions->mode : be sure this makes sense
     247        pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
     248        if (!pcm) return false;
     249
     250        // we are fitting only PM_PAR_I0; the shape elements are generated from a grid
     251        psF32 *PAR = pcm->modelConv->params->data.F32;
     252
     253        // double check that the guess is carried along...
     254
     255        // I have some source guess (e0, e1, e2)
     256        psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType);
     257
     258        if (!source->galaxyFits) {
     259            source->galaxyFits = psArrayAllocEmpty(1);
     260            // source->galaxyFits = pmSourceGalaxyFitsAlloc();
     261        }
     262        pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc();
     263        psArrayAdd(source->galaxyFits, 1, galaxyFits);
     264        psFree(galaxyFits);
     265
     266        float fRmajorBest = NAN;
     267        float fRminorBest = NAN;
     268        float chisqBest = NAN;
     269        for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; fRmajor += opt->fRmajorDel) {
     270            for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel; fRminor += opt->fRminorDel) {
     271     
     272                psEllipseAxes testAxes = guessAxes;
     273                testAxes.major = guessAxes.major * fRmajor;
     274                testAxes.minor = guessAxes.minor * fRminor;
     275               
     276                pmPSF_AxesToModel (PAR, testAxes, modelType);
     277               
     278                psphotGalaxyShapeSource (pcm, source, galaxyFits, maskVal, psfSize, true);
     279
     280                int i = galaxyFits->chisq->n - 1;
     281                float flux = galaxyFits->Flux->data.F32[i];
     282                if (isfinite(flux)) {
     283                    float thisChisq = galaxyFits->chisq->data.F32[i];
     284                    if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
     285                        chisqBest = thisChisq;
     286                        fRmajorBest = fRmajor;
     287                        fRminorBest = fRminor;
     288                    }
    284289                }
     290                // reset I0 to avoid potential problems on the next iteration
     291                PAR[PM_PAR_I0] = 1.0;
    285292            }
    286             // reset I0 to avoid potential problems on the next iteration
    287             PAR[PM_PAR_I0] = 1.0;
    288         }
    289     }
    290 
    291 #define SAVE_BEST_MODEL
    292 #ifdef SAVE_BEST_MODEL
    293     // Save model with smallest chisq
    294     if (isfinite(chisqBest)) {
    295 #else
    296     // Save model with nominal parameters
    297     {
    298         fRmajorBest = 1;
    299         fRminorBest = 1;
     293        }
     294
     295        if (isfinite(chisqBest)) {
     296            // now save the best fitting model as the source's extended model ...
     297            psEllipseAxes testAxes = guessAxes;
     298
     299            // ... unless this macro is defined
     300#ifndef SAVE_NOMINAL_MODEL
     301            testAxes.major = guessAxes.major * fRmajorBest;
     302            testAxes.minor = guessAxes.minor * fRminorBest;
    300303#endif
    301         // now save the best fitting model as the source's extended model
    302         psEllipseAxes testAxes = guessAxes;
    303         testAxes.major = guessAxes.major * fRmajorBest;
    304         testAxes.minor = guessAxes.minor * fRminorBest;
    305        
    306         pmPSF_AxesToModel (PAR, testAxes, modelType);
    307        
    308         psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, false);
    309 
    310         // Replace modelEXT with this model, if the model is good
    311         if (isfinite(PAR[PM_PAR_I0])) {
    312             psFree (source->modelEXT);
    313 
    314             source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
    315             source->type = PM_SOURCE_TYPE_EXTENDED;
    316             source->mode |= PM_SOURCE_MODE_EXTMODEL;
    317             source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    318 
    319             // cache the model flux
    320             pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
    321         }
    322     }
    323 
    324     psFree (pcm);
     304           
     305            pmPSF_AxesToModel (PAR, testAxes, modelType);
     306               
     307            psphotGalaxyShapeSource (pcm, source, galaxyFits, maskVal, psfSize, false);
     308
     309            // Replace modelEXT with the best model from the first of the model fits if one of them is good
     310            if (isfinite(PAR[PM_PAR_I0]) && iModel == 0) {
     311                psFree (source->modelEXT);
     312
     313                source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
     314                source->type = PM_SOURCE_TYPE_EXTENDED;
     315                source->mode |= PM_SOURCE_MODE_EXTMODEL;
     316                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
     317
     318                // cache the model flux
     319                pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
     320            }
     321        }
     322
     323        psFree (pcm);
     324    }
    325325    return true;
    326326}
    … …  
    328328// fit the given model to the source and find chisq & normalization
    329329// XXX is this a single-component model? sersic with a supplied index, Reff, axis ratio, and theta?
    330 bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize, bool saveResults) {
     330bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, pmSourceGalaxyFits *galaxyFits, psImageMaskType maskVal, int psfSize, bool saveResults) {
    331331
    332332    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    384384
    385385    if (saveResults) {
    386         psVectorAppend (source->galaxyFits->Flux, flux);
    387         psVectorAppend (source->galaxyFits->dFlux, dflux);
    388         psVectorAppend (source->galaxyFits->chisq, Chisq);
    389         source->galaxyFits->nPix = nPix;
     386        psVectorAppend (galaxyFits->Flux, flux);
     387        psVectorAppend (galaxyFits->dFlux, dflux);
     388        psVectorAppend (galaxyFits->chisq, Chisq);
     389        galaxyFits->nPix = nPix;
     390        galaxyFits->modelType = pcm->modelConv->type;
    390391    }
    391392
Note: See TracChangeset for help on using the changeset viewer.