IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 7, 2014, 5:38:26 AM (12 years ago)
Author:
eugene
Message:

merge changes from trunk

Location:
branches/eam_branches/ipp-20140423/psphot/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20140423/psphot/src

  • branches/eam_branches/ipp-20140423/psphot/src/psphot.h

    r36623 r36828  
    432432
    433433typedef struct {
    434     float fRmajorMin;
    435     float fRmajorMax;
    436     float fRmajorDel;
    437     float fRminorMin;
    438     float fRminorMax;
    439     float fRminorDel;
     434    float   Q;
     435    float   NSigma;
     436    float   clampSN;
     437    int     extModelType;
    440438} psphotGalaxyShapeOptions;
    441439
    … …  
    537535bool psphotGalaxyShape_Threaded (psThreadJob *job);
    538536bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize);
    539 bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize, bool saveResults);
     537bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, pmSourceGalaxyFits *galaxyFits, psImageMaskType maskVal, int psfSize, bool saveResults);
    540538psphotGalaxyShapeOptions *psphotGalaxyShapeOptionsAlloc();
     539bool psphotGalaxyShapeOptionsSet(pmSource *source, pmSourceGalaxyFits *galaxyFits, psphotGalaxyShapeOptions *defaultOptions);
    541540
    542541bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule);

  • branches/eam_branches/ipp-20140423/psphot/src/psphotFullForceSummary.c

    r36666 r36828  
    186186    psFitsClose(fits);
    187187
    188     // Convert to an array indexed by ID
    189     psArray *sortedTable = psArrayAlloc(4*inTable->n);
     188    // Convert to a set of arrays indexed by model type + 1
     189    // which contain pointers to arrays indexed by ID
     190#define MAX_MODEL_TYPE 10
     191    psArray *sortedTables = psArrayAlloc(MAX_MODEL_TYPE+1);
    190192    for (int i=0; i<inTable->n; i++) {
    191193        psMetadata *row = inTable->data[i];
    192194        psS32 ID = psMetadataLookupS32(&status, row, "ID");
     195        psS32 modelType = psMetadataLookupS32(&status, row, "MODEL_TYPE");
     196        // XXX: need to use the lookup table functions to be ready for changes in the model type numbers
     197        if (modelType+1 >= MAX_MODEL_TYPE) {
     198            psError(PSPHOT_ERR_CONFIG, false, "found modelType %d max allowed is %d", modelType, MAX_MODEL_TYPE);
     199            return false;
     200        }
     201        psArray *sortedTable = sortedTables->data[modelType+1];
     202        if (!sortedTable) {
     203            sortedTable = psArrayAlloc(4*inTable->n);
     204            sortedTables->data[modelType+1] = sortedTable;
     205            // dont' free sortedTable the array of tables gets our reference
     206        }
    193207        if (ID >= sortedTable->n) {
    194208            sortedTable = psArrayRealloc(sortedTable, 2*ID);
    195209            // Why doesn't psArrayRealloc do this?????
    196             sortedTable->n = sortedTable->nalloc;
     210                sortedTable->n = sortedTable->nalloc;
    197211        }
    198212        if (sortedTable->data[ID]) {
    … …  
    202216        sortedTable->data[ID] = psMemIncrRefCounter(row);
    203217    }
    204     psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLE", PS_META_REPLACE, "cff table", sortedTable);
     218    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLES", PS_META_REPLACE, "cff tables", sortedTables);
    205219    psFree(inTable);
    206     psFree(sortedTable);
     220    psFree(sortedTables);
    207221
    208222    return true;

  • branches/eam_branches/ipp-20140423/psphot/src/psphotFullForceSummaryReadout.c

    r36664 r36828  
    33
    44typedef struct {
     5    psF32   Q;
     6    psF32   NSigma;
     7    psF32   clampSN;
     8#ifdef notdef
     9    // these are per object now
    510    int     numTrials;
    6     psF32   fRmajorMin;
    7     psF32   fRmajorMax;
    8     psF32   fRmajorDel;
    9     psF32   fRminorMin;
    10     psF32   fRminorMax;
    11     psF32   fRminorDel;
    12     psVector    *rMajor;
    13     psVector    *rMinor;
    14     psArray *zeroPt;
    15     psArray *exptime;
    16     psArray *cffTable;
    17 } galaxyShapeOptions;
    18 
    19 static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options);
     11    psF64   fRmajorMin;
     12    psF64   fRmajorMax;
     13    psF64   fRmajorDel;
     14    psF64   fRminorMin;
     15    psF64   fRminorMax;
     16    psF64   fRminorDel;
     17    psVector    *fRmajor;
     18    psVector    *fRminor;
     19#endif
     20    psArray *zeroPt;    // zero points for each input used with exptime to scale flux
     21    psArray *exptime;   // exposure times for each input
     22    psArray *cffTables; // one for each model type index is (model_type + 1) entry 0 is no extended model (star)
     23} galaxyShapeSummaryOptions;
     24
     25
     26
     27static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeSummaryOptions *options);
    2028static pmPhotObj *findObjectForSource(psArray **pObjects, pmSource *source);
    2129
    22 static bool setOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe, bool saveVectors);
    23 static bool checkOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe);
     30static bool setOptions(galaxyShapeSummaryOptions *options, pmReadout *readout, psMetadata *recipe, bool saveVectors);
     31static bool checkOptions(galaxyShapeSummaryOptions *options, pmReadout *readout, psMetadata *recipe);
    2432
    2533
    … …  
    4755    psMetadataAddF32(outputCell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)",
    4856        outputExptime);
    49     // don't think this one matters
    50 //    psMetadataAddF32(output->fpa->concepts, PS_LIST_TAIL, "FPA.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)",
    51  //       outputExptime);
    5257
    5358    // Create objects from the various input's sources
    5459    // loop over the available readouts
    55 //    psArray *readouts = psArrayAlloc(num);
     60
    5661    psVector *fluxScaleFactor = psVectorAlloc(num, PS_TYPE_F32);
    57     galaxyShapeOptions options;
     62    galaxyShapeSummaryOptions options;
    5863    psArray *objects = NULL;
    5964    for (int index = 0; index < num; index++) {
    … …  
    7277                return false;
    7378            }
    74             options.cffTable = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLE");
    75             if (!options.cffTable) {
     79            options.cffTables = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLES");
     80            if (!options.cffTables) {
    7681                psError (PS_ERR_UNKNOWN, true, "Cannot find cff table in arguments.");
    7782                return false;
    … …  
    9196        psF32 scaleFactor = pow(10, 0.4 * (outputZeroPoint - zero_point)) * outputExptime / exptime;
    9297        fluxScaleFactor->data.F32[index] = scaleFactor;
    93 
    94 //        readouts->data[index] = readout;
    9598
    9699        // find detections
    … …  
    177180}
    178181
    179 static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options) {
     182static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeSummaryOptions *options) {
    180183
    181184    pmSource *outSrc = NULL;
    182185
    183     psVector *sumWeightedFlux = NULL;
    184     psVector *sumInvSig2 = NULL;
    185     psVector *numerator   = NULL;
    186 //    psVector *tmp = NULL;
    187     psF32   totalNPix = 0;
    188     long    vectorLength = 0;
    189 
     186
     187    pmSource *source0 = obj->sources->data[0];
     188    int ID = source0->seq;
     189
     190    // Find the parameters from the cff tables. nModels is the number of extended model types that
     191    // this source has an entry for
     192    int nModels = 0;
     193    // skip table 0 which is modelType -1 not extended
     194    for (int i=1; i<options->cffTables->n; i++) {
     195        psArray *table = options->cffTables->data[i];
     196        if (!table) continue;
     197        if (ID >= table->n) continue;
     198        if (table->data[ID]) {
     199            nModels++;
     200        }
     201    }
     202
     203    // Loop over sources for this object. Start output source for first input that has
     204    // good pixel frac
    190205    for (int i=0; i < obj->sources->n; i++) {
    191206        pmSource *source = obj->sources->data[i];
    … …  
    212227                outSrc->extpars =  psMemIncrRefCounter(source->extpars);
    213228            }
    214             if (source->modelEXT) {
    215                 outSrc->modelEXT =  psMemIncrRefCounter(source->modelEXT);
    216             }
    217             if (source->modelFits) {
    218                 outSrc->modelFits =  psMemIncrRefCounter(source->modelFits);
    219             }
    220         }
    221 
    222         if (source->galaxyFits && isfinite(source->galaxyFits->nPix) && source->galaxyFits->chisq->n) {
    223             if (numerator == NULL) {
    224                 vectorLength = source->galaxyFits->chisq->n;
    225                 // tmp = psVectorAlloc(vectorLength, PS_TYPE_F32);
    226                 sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32);
    227                 psVectorInit(sumWeightedFlux, 0.0);
    228                 sumInvSig2 = psVectorAlloc(vectorLength, PS_TYPE_F32);
    229                 psVectorInit(sumInvSig2, 0.0);
    230                 numerator   = psVectorAlloc(vectorLength, PS_TYPE_F32);
    231                 psVectorInit(numerator, 0.0);
    232             }
    233 
    234             psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d",
    235                 vectorLength, options->numTrials);
    236             psAssert(source->galaxyFits->chisq->n == vectorLength, "length of chisq vectors do not match %ld %ld",
    237                     source->galaxyFits->chisq->n, vectorLength);
    238 
    239 
    240             psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID];
     229        }
     230        if (nModels && source->modelFits && source->modelFits->n) {
     231            if (!outSrc->modelFits) {
     232                outSrc->modelFits = psArrayAllocEmpty(nModels);
     233            }
     234            for (int iModel = 0; iModel < source->modelFits->n; iModel++) {
     235                pmModel *newModel = source->modelFits->data[iModel];
     236                pmModel *outModel = NULL;
     237                for (int jModel = 0; jModel < outSrc->modelFits->n; jModel++) {
     238                    outModel = outSrc->modelFits->data[jModel];
     239                    if (newModel->type == outModel->type)  {
     240                        // already have an entry for this type
     241                        break;
     242                    }
     243                    outModel = NULL;
     244                }
     245                if (!outModel) {
     246                    // no model of this type in output source yet. Add this one.
     247                    // Note: the parameters that we need (position, angle, index, if applicable)
     248                    // are the same for all sources so copying any of them is fine
     249                    psArrayAdd(outSrc->modelFits, 1, newModel);
     250                }
     251            }
     252        }
     253    }
     254
     255    if (!outSrc) {
     256        // no good measurements for this source (how?)
     257        return NULL;
     258    }
     259
     260    // now loop over the model fits and summarize the galaxy shape measurements for each
     261    for (int iModel=0; outSrc->modelFits && iModel < outSrc->modelFits->n; iModel++) {
     262        pmModel *model = outSrc->modelFits->data[iModel];
     263        int modelType = model->type;
     264        psAssert(modelType >= 0 && modelType+1 < options->cffTables->n, "model type out of range");
     265
     266        psArray *cffTable = options->cffTables->data[modelType+1];
     267        psAssert(cffTable, "missing cff table?");
     268
     269        psMetadata *cffRow = cffTable->data[outSrc->seq];
     270        if (!cffRow) {
     271            psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d model %d", outSrc->seq, modelType);
     272            psFree(outSrc);
     273            return NULL;
     274        }
     275        bool mdok;
     276        psF32 R_MAJOR = psMetadataLookupF32(&mdok, cffRow, "R_MAJOR");
     277        if (!mdok) {
     278            psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d type %d", outSrc->seq, modelType);
     279            psFree(outSrc);
     280            return NULL;
     281        }
     282        psF32 R_MINOR = psMetadataLookupF32(&mdok, cffRow, "R_MINOR");
     283        if (!mdok) {
     284            psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d type %d", outSrc->seq, modelType);
     285            psFree(outSrc);
     286            return NULL;
     287        }
     288
     289        psVector *sumWeightedFlux = NULL;
     290        psVector *sumInvSig2 = NULL;
     291        psVector *numerator   = NULL;
     292        psVector *fRmajorVec = NULL;
     293        psVector *fRminorVec = NULL;
     294        psF32   totalNPix = 0;
     295        long    vectorLength = 0;
     296        psF32   fRmajorMin = 0;
     297        psF32   fRmajorMax = 0;
     298        psF32   fRmajorDel = 0;
     299        psF32   fRminorMin = 0;
     300        psF32   fRminorMax = 0;
     301        psF32   fRminorDel = 0;
     302        for (int i=0; i < obj->sources->n; i++) {
     303            pmSource *source = obj->sources->data[i];
     304
     305            if (source->pixWeightNotPoor < .9) continue;
     306
     307            // The only parameters that we are summarizing currently are galaxy shapes.
     308            if (!source->modelFits) continue;
     309            if (!source->galaxyFits) continue;
     310
     311            // XXX: put this into a function
     312            pmSourceGalaxyFits *galaxyFits = NULL;
     313            for (int j=0; j<source->galaxyFits->n; j++) {
     314                galaxyFits = source->galaxyFits->data[j];
     315                if (galaxyFits && galaxyFits->modelType == modelType) break;
     316                galaxyFits = NULL;
     317            }
     318
     319            if (galaxyFits && isfinite(galaxyFits->nPix) && galaxyFits->chisq->n) {
     320                if (numerator == NULL) {
     321                    // first source with galaxyFits allocate accumulators
     322                    vectorLength = galaxyFits->chisq->n;
     323                    sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32);
     324                    psVectorInit(sumWeightedFlux, 0.0);
     325                    sumInvSig2 = psVectorAlloc(vectorLength, PS_TYPE_F32);
     326                    psVectorInit(sumInvSig2, 0.0);
     327                    numerator   = psVectorAlloc(vectorLength, PS_TYPE_F32);
     328                    psVectorInit(numerator, 0.0);
     329
     330                    // Initialize the fractional radius vectors
     331                    // save these to initialize the output galaxy fits
     332                    fRmajorMin = galaxyFits->fRmajorMin;
     333                    fRmajorMax = galaxyFits->fRmajorMax;
     334                    fRmajorDel = galaxyFits->fRmajorDel;
     335                    fRminorMin = galaxyFits->fRminorMin;
     336                    fRminorMax = galaxyFits->fRminorMax;
     337                    fRminorDel = galaxyFits->fRminorDel;
     338                    fRmajorVec = psVectorAlloc(vectorLength, PS_TYPE_F32);
     339                    fRminorVec = psVectorAlloc(vectorLength, PS_TYPE_F32);
     340                    int k = 0;
     341                    for (float fRmajor = fRmajorMin; fRmajor < fRmajorMax + 0.5*fRmajorDel; fRmajor += fRmajorDel) {
     342                        for (float fRminor = fRminorMin; fRminor < fRminorMax + 0.5*fRminorDel; fRminor += fRminorDel) {
     343                            fRminorVec->data.F32[k] = fRminor;
     344                            fRmajorVec->data.F32[k] = fRmajor;
     345                            k++;
     346                        }
     347                    }
     348                    psAssert(k == vectorLength, "mismatched vectors");
     349                }
     350
     351                // Die if the lengths of the vectors is not the same in all sources
     352#ifdef notdef
     353                psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d",
     354                    vectorLength, options->numTrials);
     355#endif
     356                psAssert(galaxyFits->chisq->n == vectorLength, "length of chisq vectors do not match %ld %ld",
     357                             galaxyFits->chisq->n, vectorLength);
     358
     359                psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID];
     360
     361                totalNPix += galaxyFits->nPix;
     362
     363                for (int k = 0; k < vectorLength; k++) {
     364                    psF32 chisq = galaxyFits->chisq->data.F32[k];
     365                    psF32 flux  = galaxyFits->Flux->data.F32[k]  * scaleFactor;
     366                    psF32 dFlux = galaxyFits->dFlux->data.F32[k] * scaleFactor;
     367
     368                    numerator->data.F32[k] += chisq * galaxyFits->nPix;
     369
     370                    psF32 invSig2 = 1.0 / (dFlux * dFlux);
     371                    sumInvSig2->data.F32[k] += invSig2;
     372
     373                    sumWeightedFlux->data.F32[k] += flux * invSig2;
     374                }
     375            }
     376        }
     377
     378        if (vectorLength) {
     379            // allocate galaxyFits for the output source
     380            if (!outSrc->galaxyFits) {
     381                outSrc->galaxyFits = psArrayAllocEmpty(1);
     382            }
     383            pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc();
     384            psArrayAdd(outSrc->galaxyFits, 1, galaxyFits);
     385            psFree(galaxyFits);
     386
     387            galaxyFits->nPix  = totalNPix;
     388            galaxyFits->modelType = modelType;
     389            psVector *fluxVec =
     390                galaxyFits->Flux  = psVectorRecycle(galaxyFits->Flux,  vectorLength, PS_TYPE_F32);
     391            psVector *dFluxVec =
     392                galaxyFits->dFlux = psVectorRecycle(galaxyFits->dFlux, vectorLength, PS_TYPE_F32);
     393            psVector *chisqVec =
     394                galaxyFits->chisq = psVectorRecycle(galaxyFits->chisq, vectorLength, PS_TYPE_F32);
     395
     396            galaxyFits->fRmajorMin = fRmajorMin;
     397            galaxyFits->fRmajorMax = fRmajorMax;
     398            galaxyFits->fRmajorDel = fRmajorDel;
     399            galaxyFits->fRminorMin = fRminorMin;
     400            galaxyFits->fRminorMax = fRminorMax;
     401            galaxyFits->fRminorDel = fRminorDel;
     402
     403            // fill the summary galaxyFits vectors and find the trial with the minimum value for chisq
     404
     405            int min_k = -1;
     406            psF32 minChisq = NAN;
    241407
    242408            for (int k = 0; k < vectorLength; k++) {
    243 
    244                 psF32 chisq = source->galaxyFits->chisq->data.F32[k];
    245                 psF32 flux  = source->galaxyFits->Flux->data.F32[k]  * scaleFactor;
    246                 psF32 dFlux = source->galaxyFits->dFlux->data.F32[k] * scaleFactor;
    247                 // for chisq
    248                 // Numerator = sum( nPix[i] * chisq[i] )
    249                 // denominator will be sum( nPix )
    250 
    251 
    252 //                tmp = (psVector *) psBinaryOp(tmp,
    253  //                   source->galaxyFits->chisq, "*", psScalarAlloc(source->galaxyFits->nPix, PS_TYPE_F32));
    254 //                numerator = (psVector *) psBinaryOp(numerator, numerator, "+", tmp);
    255 
    256                 numerator->data.F32[k] += chisq * source->galaxyFits->nPix;
    257                 totalNPix += source->galaxyFits->nPix;
    258 
    259                 // sumInvSig2 = sum( 1 / dFlux**2 )
    260 //                tmp = (psVector *) psBinaryOp(tmp, source->galaxyFits->dFlux, "*", source->galaxyFits->dFlux);
    261 //                tmp = (psVector *) psBinaryOp(tmp, psScalarAlloc(1.0, PS_TYPE_F32), "/", tmp);
    262 //                sumInvSig2 = (psVector *) psBinaryOp(sumInvSig2, sumInvSig2, "+", tmp);
    263 //
    264                 psF32 invSig2 = 1.0 / (dFlux * dFlux);
    265                 sumInvSig2->data.F32[k] += invSig2;
    266 
    267                 // sumWeightedFlux = sum ( Flux / dFlux**2 )
    268 //                tmp = (psVector *) psBinaryOp(tmp, source->galaxyFits->Flux, "*", tmp);
    269 //                sumWeightedFlux = (psVector *) psBinaryOp(sumWeightedFlux, sumWeightedFlux, "+", tmp);
    270                 sumWeightedFlux->data.F32[k] += flux * invSig2;
    271             }
    272         }
    273     }
    274 
    275     if (outSrc) {
    276         if (vectorLength) {
    277             outSrc->galaxyFits = pmSourceGalaxyFitsAlloc();
    278             outSrc->galaxyFits->nPix = totalNPix;
    279             outSrc->galaxyFits->Flux  = psVectorRecycle(outSrc->galaxyFits->Flux, vectorLength, PS_TYPE_F32);
    280             outSrc->galaxyFits->dFlux = psVectorRecycle(outSrc->galaxyFits->dFlux, vectorLength, PS_TYPE_F32);
    281             outSrc->galaxyFits->chisq = psVectorRecycle(outSrc->galaxyFits->chisq, vectorLength, PS_TYPE_F32);
    282             for (int k = 0; k < vectorLength; k++) {
    283                 outSrc->galaxyFits->Flux->data.F32[k]  = sumWeightedFlux->data.F32[k] / sumInvSig2->data.F32[k];
    284                 outSrc->galaxyFits->dFlux->data.F32[k] = 1.0 / sumInvSig2->data.F32[k];
    285                 outSrc->galaxyFits->chisq->data.F32[k] = numerator->data.F32[k] / totalNPix;
    286 
    287 #ifdef nodef
    288                 outSrc->galaxyFits->Flux = (psVector *) psBinaryOp(outSrc->galaxyFits->Flux,
    289                         sumWeightedFlux, "/", sumInvSig2);
    290 
    291                 outSrc->galaxyFits->dFlux = (psVector *) psBinaryOp(outSrc->galaxyFits->dFlux,
    292                         psScalarAlloc(1.0, PS_TYPE_F32), "/", sumInvSig2);
    293 
    294                 outSrc->galaxyFits->chisq = (psVector *) psBinaryOp(outSrc->galaxyFits->chisq,
    295                         numerator, "/", psScalarAlloc(totalNPix, PS_TYPE_F32));
    296 #endif
     409                fluxVec->data.F32[k]  = sumWeightedFlux->data.F32[k] / sumInvSig2->data.F32[k];
     410                dFluxVec->data.F32[k] = 1.0 / sumInvSig2->data.F32[k];
     411
     412                psF32 thischisq = chisqVec->data.F32[k] = numerator->data.F32[k] / totalNPix;
     413
     414                if (isfinite(thischisq)  && (!isfinite(minChisq) || thischisq < minChisq)) {
     415                    min_k = k;
     416                    minChisq = thischisq;
     417                }
    297418            }
    298419
    … …  
    300421            psFree(sumInvSig2);
    301422            psFree(sumWeightedFlux);
    302 //            psFree(tmp);
    303 
    304             int min_j = -1;
    305             psF32 minChisq = NAN;
    306             psVector *chisq = outSrc->galaxyFits->chisq;
    307             for (int j=0; j < chisq->n; j++) {
    308                 psF32 thischisq = chisq->data.F32[j];
    309                 if (isfinite(thischisq)  && (!isfinite(minChisq) || thischisq < minChisq)) {
    310                     min_j = j;
    311                     minChisq = thischisq;
     423
     424            if (min_k >= 0 && isfinite(minChisq)) {
     425                // copy the best fit params to the model
     426                // fractional radii with the lowest chisq
     427                psEllipseAxes axes = pmPSF_ModelToAxes(model->params->data.F32, model->type);
     428
     429                // examine the params for the trial with minimum chisq.
     430                bool fitMajor = true;
     431                bool fitMinor = true;
     432                psF64 fRmajorBest = fRmajorVec->data.F32[min_k];
     433                if ((fabs(fRmajorBest - galaxyFits->fRmajorMin) < galaxyFits->fRmajorDel) ||
     434                    (fabs(fRmajorBest - galaxyFits->fRmajorMax) < galaxyFits->fRmajorDel)) {
     435                    fitMajor = false;
    312436                }
    313             }
    314             if (min_j >= 0 && isfinite(minChisq) && outSrc->modelEXT) {
    315                 // copy the best fit params to the model
    316                 psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type);
    317                 if (outSrc->seq >= options->cffTable->n) {
    318                     psError(PS_ERR_PROGRAMMING, true, "object sequence number %d is larger than cff table length %ld",
    319                         outSrc->seq, options->cffTable->n);
    320                     return NULL;
     437                psF64 fRminorBest = fRminorVec->data.F32[min_k];
     438                if ((fabs(fRminorBest - galaxyFits->fRminorMin) < galaxyFits->fRminorDel) ||
     439                    (fabs(fRminorBest - galaxyFits->fRminorMax) < galaxyFits->fRminorDel)) {
     440                    fitMinor = false;
    321441                }
    322                 psMetadata *row = options->cffTable->data[outSrc->seq];
    323                 if (!row) {
    324                     psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d", outSrc->seq);
    325                     psFree(outSrc);
    326                     return NULL;
     442                // If either major or minor is at one of the limits do not fit report the minimum value
     443                bool useFit = fitMajor && fitMinor;
     444
     445                // save the flux and dFlux values from entry with lowest chisq
     446                psF64 fluxBest = fluxVec->data.F32[min_k];
     447                psF64 dFluxBest = dFluxVec->data.F32[min_k];
     448                psF64 dFlux0   = NAN;
     449                psF64 flux0    = NAN;
     450
     451                if (useFit) {
     452                    #define NUM_TRIALS_INIT 9
     453                    psVector *major = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64);
     454                    psVector *minor = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64);
     455                    psVector *chisq = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64);
     456                    psVector *flux = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64);
     457                    psVector *dFlux = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64);
     458
     459                    // XXX: use a recipe parameter instead of 2.2
     460                    psF64 maxDeltaMaj = 2.2 * galaxyFits->fRmajorDel;
     461                    psF64 maxDeltaMin = 2.2 * galaxyFits->fRminorDel;
     462
     463                    psF64 majorMin = NAN;
     464                    psF64 majorMax = NAN;
     465                    psF64 minorMin = NAN;
     466                    psF64 minorMax = NAN;
     467                    psF64 chisqMin = NAN;
     468                    psF64 chisqMax = NAN;
     469                    for (int k = 0; k < vectorLength; k++) {
     470
     471                        if (fabs(fRmajorVec->data.F32[k] - fRmajorBest) < maxDeltaMaj &&
     472                            fabs(fRminorVec->data.F32[k] - fRminorBest) < maxDeltaMin) {
     473
     474                            if (isfinite(chisqVec->data.F32[k]) &&
     475                                isfinite(fluxVec->data.F32[k])  &&
     476                                isfinite(dFluxVec->data.F32[k])) {
     477
     478                                // compute major and minor radius vectors from nominal and trial fractions
     479                                // also find the ranges in the vectors
     480                                psF64 thisMajor = R_MAJOR * fRmajorVec->data.F32[k];
     481                                if (!isfinite(majorMin) || thisMajor < majorMin) {
     482                                    majorMin = thisMajor;
     483                                }
     484                                if (!isfinite(majorMax) || thisMajor > majorMax) {
     485                                    majorMax = thisMajor;
     486                                }
     487                                psVectorAppend(major, thisMajor);
     488
     489                                psF64 thisMinor = R_MINOR * fRminorVec->data.F32[k];
     490                                if (!isfinite(minorMin) || thisMinor < minorMin) {
     491                                    minorMin = thisMinor;
     492                                }
     493                                if (!isfinite(minorMax) || thisMinor > minorMax) {
     494                                    minorMax = thisMinor;
     495                                }
     496                                psVectorAppend(minor, thisMinor);
     497
     498                                psF64 thisChisq = chisqVec->data.F32[k];
     499                                if (!isfinite(chisqMin) || thisChisq < chisqMin) {
     500                                    chisqMin = thisChisq;
     501                                }
     502                                if (!isfinite(chisqMax) || thisChisq > chisqMax) {
     503                                    chisqMax = thisChisq;
     504                                }
     505                                psVectorAppend(chisq, thisChisq);
     506
     507                                psVectorAppend(flux,  fluxVec->data.F32[k]);
     508                                psVectorAppend(dFlux, dFluxVec->data.F32[k]);
     509                            }
     510                        }
     511                    }
     512
     513
     514                    // see if we ever get too few good values (haven't seen this happen)
     515                    if (major->n < NUM_TRIALS_INIT) {
     516                        fprintf(stderr, "only found %ld good points wanted %d for seq: %d\n",
     517                            major->n, NUM_TRIALS_INIT, outSrc->seq);
     518                    }
     519
     520                    // rescale data - this helps avoid precision errors.
     521                    psF64 majorRange = majorMax - majorMin;
     522                    if (majorRange == 0) {
     523                        majorRange = 1;
     524                    }
     525                    psF64 minorRange = minorMax - minorMin;
     526                    if (minorRange == 0) {
     527                        minorRange = 1;
     528                    }
     529                    psF64 chisqRange = chisqMax - chisqMin;
     530                    if (chisqRange == 0) {
     531                        chisqRange = 1;
     532                    }
     533                    for (int k = 0; k < major->n; k++) {
     534                        major->data.F64[k] = (major->data.F64[k] - majorMin) / majorRange;
     535                        minor->data.F64[k] = (minor->data.F64[k] - minorMin) / minorRange;
     536                        chisq->data.F64[k] = (chisq->data.F64[k] - chisqMin) / chisqRange;
     537                    }
     538
     539                    // Fit chisq versus rMajor and rMinor
     540                    psPolynomial2D *chisqFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2);
     541                    psPolynomial2D *fluxFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2);
     542                    psPolynomial2D *fluxErrorFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2);
     543                    chisqFit->coeffMask[2][2] = PS_POLY_MASK_SET;
     544                    chisqFit->coeffMask[2][1] = PS_POLY_MASK_SET;
     545                    chisqFit->coeffMask[1][2] = PS_POLY_MASK_SET;
     546                    fluxFit->coeffMask[2][2] = PS_POLY_MASK_SET;
     547                    fluxFit->coeffMask[2][1] = PS_POLY_MASK_SET;
     548                    fluxFit->coeffMask[1][2] = PS_POLY_MASK_SET;
     549                    fluxErrorFit->coeffMask[2][2] = PS_POLY_MASK_SET;
     550                    fluxErrorFit->coeffMask[2][1] = PS_POLY_MASK_SET;
     551                    fluxErrorFit->coeffMask[1][2] = PS_POLY_MASK_SET;
     552
     553                    bool goodFit = psVectorFitPolynomial2D (chisqFit, NULL, 0xff, chisq, NULL, major, minor);
     554
     555                    if (goodFit) {
     556                        // find the axes for the minimum chisq from the formula
     557                        psF64 **coeff = chisqFit->coeff;
     558                        psF64 D = 4*coeff[2][0]*coeff[0][2] - coeff[1][1]*coeff[1][1];
     559                        psF64 major0 = ((coeff[1][1]*coeff[0][1] - 2*coeff[0][2]*coeff[1][0]) / D) * majorRange + majorMin;
     560                        psF64 minor0 = ((coeff[1][1]*coeff[1][0] - 2*coeff[2][0]*coeff[0][1]) / D) * minorRange + minorMin;
     561
     562                        // estimated chisq at the minimum
     563                        psF64 chisq0 = psPolynomial2DEval(chisqFit, major0, minor0) * chisqRange + chisqMin;
     564
     565                        // now fit the flux ...
     566                        bool goodFluxFit = psVectorFitPolynomial2D(fluxFit, NULL, 0xFF, flux, NULL, major, minor);
     567
     568                        // .. and compute flux at the minimum chisq position
     569                        if (goodFluxFit) {
     570                            flux0  = psPolynomial2DEval(fluxFit, major0, minor0);
     571                        } else {
     572                            flux0 = fluxBest;
     573                        }
     574
     575                        // .. and compute dFlux at the minimum chisq position
     576                        bool goodFluxErrorFit = psVectorFitPolynomial2D(fluxErrorFit, NULL, 0xFF, dFlux, NULL, major, minor);
     577                        if (goodFluxErrorFit) {
     578                            dFlux0  = psPolynomial2DEval(fluxErrorFit, major0, minor0);
     579                        } else {
     580                            dFlux0 = dFluxBest;
     581                        }
     582
     583#ifdef PRINTVALS
     584                        fprintf (stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %4.2f |  %6.3f %6.3f %6.3f %4.2f | %7.4f %7.4f %7.1f %7.1f\n",
     585                            outSrc->seq, min_k, major->n, R_MAJOR, fRmajorBest*R_MAJOR, major0, major0/R_MAJOR, R_MINOR, fRminorBest*R_MINOR, minor0, minor0/R_MINOR, minChisq, chisq0, fluxBest, flux0);
     586#endif
     587
     588
     589                        axes.major = major0;
     590                        axes.minor = minor0;
     591                        model->chisq = chisq0;
     592#ifdef DUMPVECTORS
     593                        char fn[80];
     594
     595                        sprintf(fn, "vectors/%s.%05d.txt", goodFit ? "g" : "b", outSrc->seq);
     596
     597                        FILE *f = fopen(fn, "w");
     598                        fprintf(f, "#major minor chisq i\n");
     599                        for (int k = 0; k < chisq->n; k++) {
     600                            fprintf(f, "%10.6f %10.6f %10.6f %5d\n",
     601                                major->data.F64[k]/R_MAJOR, minor->data.F64[k]/R_MINOR, chisq->data.F64[k], k);
     602                        }
     603                        fclose(f);
     604#endif
     605                    } else {
     606                        model->flags |= PM_MODEL_STATUS_NONCONVERGE;
     607
     608#ifdef PRINTVALS
     609                        fprintf(stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %6.3f bad fit\n",
     610                            outSrc->seq, min_k, major->n, R_MAJOR, fRmajorBest, R_MINOR, fRminorBest);
     611#endif
     612                        psErrorClear();
     613                        // psFree(outSrc->galaxyFits);
     614
     615                    }
     616                    psFree(chisqFit);
     617                    psFree(fluxFit);
     618                    psFree(fluxErrorFit);
     619                    psFree(major);
     620                    psFree(minor);
     621                    psFree(chisq);
     622                    psFree(flux);
     623                    psFree(dFlux);
     624                } else {
     625                    // No fit
     626                    // Set the values based on the trial with the best chisq
     627                    axes.major = R_MAJOR * fRmajorBest;
     628                    axes.minor = R_MINOR * fRminorBest;
     629                    model->chisq = minChisq;
     630                    flux0 = fluxBest;
     631                    dFlux0 = dFluxBest;
     632#ifdef PRINTVALS
     633                    fprintf(stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %6.3f skip fit\n",
     634                            outSrc->seq, min_k, 0L, R_MAJOR, fRmajorBest, R_MINOR, fRminorBest);
     635#endif
    327636                }
    328                 bool mdok;
    329                 psF32 rMajor = psMetadataLookupF32(&mdok, row, "R_MAJOR");
    330                 if (!mdok) {
    331                     psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d", outSrc->seq);
    332                     psFree(outSrc);
    333                     return NULL;
    334                 }
    335                 psF32 rMinor = psMetadataLookupF32(&mdok, row, "R_MINOR");
    336                 if (!mdok) {
    337                     psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d", outSrc->seq);
    338                     psFree(outSrc);
    339                     return NULL;
    340                 }
    341                 axes.major = rMajor * options->rMajor->data.F32[min_j];
    342                 axes.minor = rMinor * options->rMinor->data.F32[min_j];
    343                 pmPSF_AxesToModel (outSrc->modelEXT->params->data.F32, axes, outSrc->modelEXT->type);
    344                 outSrc->modelEXT->chisq = minChisq;
    345                 outSrc->modelEXT->mag = -2.5 * log10(outSrc->galaxyFits->Flux->data.F32[min_j]);
    346                 outSrc->modelEXT->magErr = 1.0 /
    347                     (outSrc->galaxyFits->Flux->data.F32[min_j]/outSrc->galaxyFits->dFlux->data.F32[min_j]); // 1 / SN
    348             }
    349         }
     637                // now save the model parameters in the model structure
     638
     639                pmPSF_AxesToModel (model->params->data.F32, axes, modelType);
     640
     641                model->mag = -2.5 * log10(flux0);
     642                model->magErr = dFlux0 / flux0; // 1 / SN
     643                // XXX: should there be a different flag if we didn't do a fit of the chisq
     644                model->flags |= PM_MODEL_STATUS_FITTED;
     645            }
     646        }
     647        psFree(fRminorVec);
     648        psFree(fRmajorVec);
    350649    }
    351650
    … …  
    360659    }
    361660
    362 static bool setOptions(galaxyShapeOptions *opt, pmReadout *readout, psMetadata *recipe, bool makeVectors) {
     661static bool setOptions(galaxyShapeSummaryOptions *opt, pmReadout *readout, psMetadata *recipe, bool makeVectors) {
    363662    bool status;
    364     bool useAnalysis;
    365 
     663    bool useAnalysis;   // fall back to recipe if we dont' find values in analysis. Probably should no longer do this
     664
     665#ifdef notdef
    366666    psMetadataLookupF32(&useAnalysis, readout->analysis, "GALAXY_SHAPES_FR_MAJOR_MIN");
     667#endif
     668    psMetadataLookupF32(&useAnalysis, readout->analysis, "GALAXY_SHAPES_Q");
    367669    psMetadata *md = useAnalysis ? readout->analysis : recipe;
    368670
     671    GETVAL(Q, "GALAXY_SHAPES_Q");
     672    GETVAL(NSigma, "GALAXY_SHAPES_NSIGMA");
     673    GETVAL(clampSN, "GALAXY_SHAPES_CLAMP_SN");
     674
     675#ifdef notdef
     676    // these are per object now
    369677    GETVAL(fRmajorMin, "GALAXY_SHAPES_FR_MAJOR_MIN");
    370678    GETVAL(fRmajorMax, "GALAXY_SHAPES_FR_MAJOR_MAX");
    … …  
    378686
    379687    if (makeVectors) {
    380         opt->rMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
    381         opt->rMajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
     688        opt->fRminor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
     689        opt->fRmajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
    382690        int i = 0;
    383691        for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel;
    … …  
    385693            for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel;
    386694                    fRminor += opt->fRminorDel) {
    387                 opt->rMinor->data.F32[i] = fRminor;
    388                 opt->rMajor->data.F32[i] = fRmajor;
     695                opt->fRminor->data.F32[i] = fRminor;
     696                opt->fRmajor->data.F32[i] = fRmajor;
    389697                i++;
    390698            }
    … …  
    392700        psAssert(i == opt->numTrials, "Something's wrong with my loop got %d entries expected %d", i, opt->numTrials);
    393701    } else {
    394         opt->rMinor = NULL;
    395         opt->rMajor = NULL;
    396     }
     702        opt->fRminor = NULL;
     703        opt->fRmajor = NULL;
     704    }
     705#endif
    397706       
    398707    return true;
    … …  
    405714    }
    406715
    407 static bool checkOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe) {
    408     galaxyShapeOptions thisReadoutsOptions;
     716static bool checkOptions(galaxyShapeSummaryOptions *options, pmReadout *readout, psMetadata *recipe) {
     717    galaxyShapeSummaryOptions thisReadoutsOptions;
    409718   
    410719    if (!setOptions(&thisReadoutsOptions, readout, recipe, false)) {
    … …  
    412721        return false;
    413722    }
     723    CHECKVAL(options, thisReadoutsOptions, Q, "mismatched Q");
     724    CHECKVAL(options, thisReadoutsOptions, NSigma, "mismatched NSIGMA");
     725    CHECKVAL(options, thisReadoutsOptions, clampSN, "mismatched cleampSN");
     726
     727#ifdef notdef
     728    // these are per object now
    414729    CHECKVAL(options, thisReadoutsOptions, numTrials, "mismatched number of trials")
    415730    CHECKVAL(options, thisReadoutsOptions, fRmajorMin, "mismatched fRmajorMin")
    … …  
    419734    CHECKVAL(options, thisReadoutsOptions, fRminorMax, "mismatched fRminorMax")
    420735    CHECKVAL(options, thisReadoutsOptions, fRminorDel, "mismatched fRminorDel")
     736#endif
    421737
    422738    return true;

  • branches/eam_branches/ipp-20140423/psphot/src/psphotGalaxyShape.c

    r36633 r36828  
    8080    assert (markVal);
    8181
    82     // maskVal is used to test for rejected pixels, and must include markVal
    8382    maskVal |= markVal;
    8483
    85     // bit-mask to mark pixels not used in analysis
    8684    psphotGalaxyShapeOptions *opt = psphotGalaxyShapeOptionsAlloc();
     85    opt->Q = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_Q");
     86    psAssert (status, "missing GALAXY_SHAPES_Q");
     87    opt->NSigma = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_NSIGMA");
     88    psAssert (status, "missing GALAXY_SHAPES_NSIGMA");
     89    opt->clampSN = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_CLAMP_SN");
     90    psAssert (status, "missing GALAXY_SHAPES_NSIGMA");
     91    psString modelTypeToSave = psMetadataLookupStr(&status, recipe, "EXT_MODEL_TYPE_FORCE");
     92    if (modelTypeToSave && strcmp(modelTypeToSave, "ALL") && strcmp(modelTypeToSave, "BEST")) {
     93        opt->extModelType = pmModelClassGetType(modelTypeToSave);
     94    } else {
     95        opt->extModelType = pmModelClassGetType("PS_MODEL_SERSIC");
     96    }
     97
     98#ifdef notdef
    8799    opt->fRmajorMin = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MIN"); if (!status) opt->fRmajorMin = 0.5;
    88100    opt->fRmajorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MAX"); if (!status) opt->fRmajorMax = 2.0;
    … …  
    91103    opt->fRminorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_MAX"); if (!status) opt->fRminorMax = 2.0;
    92104    opt->fRminorDel = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_DEL"); if (!status) opt->fRminorDel = 0.1;
     105#endif
    93106
    94107    // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
    … …  
    238251bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize) {
    239252
    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;
     253
     254    for (int iModel = 0 ; iModel < source->modelFits->n; iModel++) {
     255        pmModel *model = source->modelFits->data[iModel];
     256        if (!model) return false;
     257
     258        pmModelType modelType = model->type;
     259        model->flags = PM_MODEL_STATUS_NONE;
     260
     261        // we are using fitOptions->mode : be sure this makes sense
     262        pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
     263        if (!pcm) return false;
     264
     265        // we are fitting only PM_PAR_I0; the shape elements are generated from a grid
     266        psF32 *PAR = pcm->modelConv->params->data.F32;
     267
     268        pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc();
     269
     270        // set the options for this source
     271        if (!psphotGalaxyShapeOptionsSet(source, galaxyFits, opt)) {
     272            // this source won't work
     273            psFree(pcm);
     274            psFree(galaxyFits);
     275            continue;
     276        }
     277
     278        if (!source->galaxyFits) {
     279            source->galaxyFits = psArrayAllocEmpty(1);
     280        }
     281        psArrayAdd(source->galaxyFits, 1, galaxyFits);
     282        psFree(galaxyFits);
     283
     284        // I have some source guess (e0, e1, e2)
     285        psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType);
     286
     287        float fRmajorBest = NAN;
     288        float fRminorBest = NAN;
     289        float chisqBest = NAN;
     290        for (float fRmajor = galaxyFits->fRmajorMin; fRmajor < galaxyFits->fRmajorMax + 0.5*galaxyFits->fRmajorDel; fRmajor += galaxyFits->fRmajorDel) {
     291            for (float fRminor = galaxyFits->fRminorMin; fRminor < galaxyFits->fRminorMax + 0.5*galaxyFits->fRminorDel; fRminor += galaxyFits->fRminorDel) {
     292     
     293                psEllipseAxes testAxes = guessAxes;
     294                testAxes.major = guessAxes.major * fRmajor;
     295                testAxes.minor = guessAxes.minor * fRminor;
     296               
     297                pmPSF_AxesToModel (PAR, testAxes, modelType);
     298               
     299                psphotGalaxyShapeSource (pcm, source, galaxyFits, maskVal, psfSize, true);
     300
     301                int i = galaxyFits->chisq->n - 1;
     302                float flux = galaxyFits->Flux->data.F32[i];
     303                if (isfinite(flux)) {
     304                    float thisChisq = galaxyFits->chisq->data.F32[i];
     305                    if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
     306                        chisqBest = thisChisq;
     307                        fRmajorBest = fRmajor;
     308                        fRminorBest = fRminor;
     309                    }
    284310                }
     311                // reset I0 to avoid potential problems on the next iteration
     312                PAR[PM_PAR_I0] = 1.0;
    285313            }
    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;
     314        }
     315
     316        if (isfinite(chisqBest)) {
     317            // now save the best fitting model as the source's extended model ...
     318            psEllipseAxes testAxes = guessAxes;
     319
     320            // ... unless this macro is defined
     321#ifndef SAVE_NOMINAL_MODEL
     322            testAxes.major = guessAxes.major * fRmajorBest;
     323            testAxes.minor = guessAxes.minor * fRminorBest;
    300324#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);
     325           
     326            pmPSF_AxesToModel (PAR, testAxes, modelType);
     327               
     328            psphotGalaxyShapeSource (pcm, source, galaxyFits, maskVal, psfSize, false);
     329
     330            // Replace modelEXT with the best model from the first of the model fits if one of them is good
     331            if (isfinite(PAR[PM_PAR_I0]) && modelType == opt->extModelType) {
     332                psFree (source->modelEXT);
     333
     334                source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
     335                source->type = PM_SOURCE_TYPE_EXTENDED;
     336                source->mode |= PM_SOURCE_MODE_EXTMODEL;
     337                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
     338
     339                // cache the model flux
     340                pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
     341            }
     342        }
     343
     344        psFree (pcm);
     345    }
    325346    return true;
    326347}
    … …  
    328349// fit the given model to the source and find chisq & normalization
    329350// 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) {
     351bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, pmSourceGalaxyFits *galaxyFits, psImageMaskType maskVal, int psfSize, bool saveResults) {
    331352
    332353    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    384405
    385406    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;
     407        psVectorAppend (galaxyFits->Flux, flux);
     408        psVectorAppend (galaxyFits->dFlux, dflux);
     409        psVectorAppend (galaxyFits->chisq, Chisq);
     410        galaxyFits->nPix = nPix;
     411        galaxyFits->modelType = pcm->modelConv->type;
    390412    }
    391413
    … …  
    404426    psMemSetDeallocator(opt, (psFreeFunc) psphotGalaxyShapeOptionsFree);
    405427   
    406     opt->fRmajorMin = 0.5;
    407     opt->fRmajorMax = 2.0;
    408     opt->fRmajorDel = 0.1;
    409     opt->fRminorMin = 0.5;
    410     opt->fRminorMax = 2.0;
    411     opt->fRminorDel = 0.1;
    412 
    413428    return opt;
    414429}
    415430
     431bool psphotGalaxyShapeOptionsSet(pmSource *source, pmSourceGalaxyFits *galaxyFits, psphotGalaxyShapeOptions *options)
     432{
     433    // XXX: put these in recipe
     434    // doesn't make sense to use too small of a spacing
     435    float clampedSN = source->extSN < options->clampSN ? source->extSN : options->clampSN;
     436
     437    float f_del = options->Q / clampedSN;
     438    float f_min = 1 - options->NSigma * f_del;
     439    float f_max = 1 + options->NSigma * f_del;
     440
     441    // if f_min goes negative skip this object
     442    if (f_min < 0) return false;
     443
     444    galaxyFits->fRmajorMin = f_min;
     445    galaxyFits->fRmajorMax = f_max;
     446    galaxyFits->fRmajorDel = f_del;
     447    galaxyFits->fRminorMin = f_min;
     448    galaxyFits->fRminorMax = f_max;
     449    galaxyFits->fRminorDel = f_del;
     450
     451    return true;
     452}
     453

  • branches/eam_branches/ipp-20140423/psphot/src/psphotKronIterate.c

    r36375 r36828  
    532532            source->moments->KronFlux  = NAN;
    533533            source->moments->KronFluxErr  = NAN;
     534            source->moments->KronFinner = NAN;
     535            source->moments->KronFouter = NAN;
    534536            return false;
    535537        }

  • branches/eam_branches/ipp-20140423/psphot/src/psphotStackImageLoop.c

    • Property svn:mergeinfo changed (with no actual effect on merging)
Note: See TracChangeset for help on using the changeset viewer.