IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34570


Ignore:
Timestamp:
Oct 24, 2012, 3:56:18 PM (14 years ago)
Author:
eugene
Message:

merge changes from trunk

Location:
branches/eam_branches/ipp-20120905/psphot
Files:
20 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20120905/psphot

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

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

    r34416 r34570  
    372372bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    373373bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS);
    374 bool psphotDropBadMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
     374bool psphotFilterMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
    375375
    376376bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final);
    … …  
    525525bool psphotSourceMemoryReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    526526
     527const char * psphotGetFilerule(const char *baseRule);
     528extern bool psphotINpsphotStack;
     529
    527530#endif

  • branches/eam_branches/ipp-20120905/psphot/src/psphotAddNoise.c

    r34415 r34570  
    3434static int Nmasked = 0;
    3535
     36// the return state indicates if any sources were actually replaced
    3637bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add) {
    3738
    … …  
    5152
    5253    psArray *sources = detections->allSources;
    53     //psAssert (sources, "missing sources?");
    54     // if no work, should just return true
    55     if (!sources) return true;
     54    // if no work to do, should just return true
     55    if (!sources) return false;
    5656
    5757    psTimerStart ("psphot.noise");
    … …  
    7474
    7575    if (SIZE <= 0) {
    76        return true;
     76       return false;
    7777    }
    7878

  • branches/eam_branches/ipp-20120905/psphot/src/psphotDeblendSatstars.c

    r34416 r34570  
    544544    float logRdel = 0.1;
    545545
    546     float Xc = source->satstar->Xo - source->pixels->col0 - 0.5;
    547     float Yc = source->satstar->Yo - source->pixels->row0 - 0.5;
    548     psVector *logRmodel = source->satstar->logRmodel;
    549     psVector *logFmodel = source->satstar->logFmodel;
     546    pmSourceSatstar *satstar = source->satstar ? source->satstar : (source->parent ? source->parent->satstar : NULL);
     547    psAssert (satstar, "null satstar");
     548    float Xc = satstar->Xo - source->pixels->col0 - 0.5;
     549    float Yc = satstar->Yo - source->pixels->row0 - 0.5;
     550    psVector *logRmodel = satstar->logRmodel;
     551    psVector *logFmodel = satstar->logFmodel;
    550552
    551553    for (int iy = 0; iy < source->pixels->numRows; iy++) {
    … …  
    10481050}
    10491051
     1052// the return state indicates if any sources were actually subtracted
    10501053bool psphotAddOrSubSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex, psMetadata *recipe, bool add) {
    10511054
    10521055    bool status;
     1056    bool modified = false;
    10531057
    10541058    psTimerStart ("psphot.deblend.sat");
    … …  
    10651069
    10661070    psArray *sources = detections->allSources;
    1067     psAssert (sources, "missing sources?");
     1071    // if no work to do, should just return false
     1072    if (!sources) return false;
    10681073
    10691074    if (!sources->n) {
    10701075        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend");
    1071         return true;
     1076        return false;
    10721077    }
    10731078
    … …  
    10821087        if (!(source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE)) continue;
    10831088
     1089        // tell the calling function that we modified the image
     1090        modified = true;
     1091
    10841092        if (!psphotSatstarProfileOp (source, maskVal, 1.0, 0, add)) continue;
    10851093    }
    10861094
    10871095    psLogMsg ("psphot", PS_LOG_DETAIL, "satstar op: %f sec\n", psTimerMark ("psphot.deblend.sat"));
    1088     return true;
     1096    return modified;
    10891097}
    10901098

  • branches/eam_branches/ipp-20120905/psphot/src/psphotDefineFiles.c

    r33963 r34570  
    11# include "psphotInternal.h"
     2
     3bool psphotINpsphotStack = false;
    24
    35// List of output files
    … …  
    181183    return;
    182184}
     185
     186// psphotGetFilerule
     187// Since psphotStack processes multipe FPAs at a time it has a different file rule structure than regular psphot.
     188// For the background output files we define a function psphotGetFilerule which given a base psphot file rule
     189// returns the corresponding psphotStack rule *if* the program is psphotStack. That is indicated by a global
     190// boolean which defaults to false, and psphotStack only sets to true
     191
     192const char *psphotGetFilerule(const char *psphotRule) {
     193    const char *rule = psphotRule;
     194    if (psphotINpsphotStack) {
     195        if (!strcmp(psphotRule, "PSPHOT.BACKMDL")) {
     196            rule =  "PSPHOT.STACK.BACKMDL";
     197        } else if (!strcmp(psphotRule, "PSPHOT.BACKMDL.STDEV")) {
     198            rule = "PSPHOT.STACK.BACKMDL.STDEV";
     199        } else if (!strcmp(psphotRule, "PSPHOT.BACKSUB")) {
     200            rule = "PSPHOT.STACK.BACKSUB";
     201        } else if (!strcmp(psphotRule, "PSPHOT.BACKGND")) {
     202            rule = "PSPHOT.STACK.BACKGND";
     203        } else {
     204            psAssert(0, "unsupported file rule %s", psphotRule);
     205        }
     206    }
     207    return rule;
     208}

  • branches/eam_branches/ipp-20120905/psphot/src/psphotFindDetections.c

    r34415 r34570  
    7171    }
    7272
     73    // first pass vs other: if this is the first pass, the code will use PEAKS_NSIGMA_LIMIT and
     74    // only attempt to detect PEAKS_NMAX entries.  If 'firstPass' is false, the code will
     75    // attempt to replace the subtracted sources in order to measure the footprints.  After
     76    // replacement, it is necessary to regenerate the significance image.  If no sources are
     77    // available, the code will skip the significance image regeneration step.
     78
    7379    bool replaceSourcesForFootprints = false;
    7480    if (firstPass) {
    … …  
    126132    if (useFootprints) {
    127133        if (replaceSourcesForFootprints) {
     134            bool modified = false;
    128135            // subtract the noise for all sources including satstars
    129             psphotAddOrSubNoiseReadout(config, view, filerule, index, recipe, false);
    130             psphotReplaceAllSourcesReadout (config, view, filerule, index, recipe, false);
     136            modified |= psphotAddOrSubNoiseReadout(config, view, filerule, index, recipe, false);
     137            modified |= psphotReplaceAllSourcesReadout (config, view, filerule, index, recipe, false);
    131138
    132139            // add in the satstars
    133             psphotAddOrSubSatstarsReadout (config, view, filerule, index, recipe, true);
     140            modified |= psphotAddOrSubSatstarsReadout (config, view, filerule, index, recipe, true);
    134141
    135             psFree (significance);
    136             significance = psphotSignificanceImage (readout, recipe, maskVal);
     142            if (modified) {
     143                psFree (significance);
     144                significance = psphotSignificanceImage (readout, recipe, maskVal);
     145            }
    137146
    138147            // display the significance image

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

    r34415 r34570  
    11# include "psphotInternal.h"
    22
    3 bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool oldWindow);
    4 bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, psImage *smoothedPixels);
    5 
    6 
    7 bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass)
    8 {
     3bool psphotKronRadiusMeasure (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal, psImage *smoothedPixels);
     4bool psphotKronFluxMeasure (pmSource *source, psImageMaskType maskVal);
     5
     6
     7bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass) {
    98    bool status = true;
    109
    … …  
    9897    }
    9998
    100     bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
    101     if (!status) {
    102         KRON_APPLY_WEIGHT = true;
    103     }
    104 
    105     bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
    106     if (!status) {
    107         KRON_APPLY_WINDOW = false;
    108     }
    10999    bool KRON_SMOOTH = psMetadataLookupBool (&status, recipe, "KRON_SMOOTH");
    110100    if (!status) {
    … …  
    170160    }
    171161
    172     // generate the window image: multiply the flux by this to downweight neighbors
    173     // XXX: we don't need this image if we aren't going to apply the window
    174     psImage *kronWindow = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
    175     psImageInit (kronWindow, 1.0);
    176 
    177     // start with the currently known moments (Mxx, Myy, Mxy) and generate a window image
    178     for (int i = 0; i < sources->n; i++) {
    179 
    180         pmSource *source = sources->data[i];
    181 
    182         // set a window function for each source based on the moments
    183         // (this skips really bad sources (no peak, no moments, DEFECT)
    184         psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW);
     162    // pass 3 we only measure fluxes for matched sources so we don't need smooth pixels or multiple iterations
     163    if (pass == 3) {
     164        KRON_SMOOTH = false;
     165        KRON_ITERATIONS = 1;
    185166    }
    186167
    … …  
    235216
    236217            psArrayAdd(job->args, 1, readout);
    237             psArrayAdd(job->args, 1, kronWindow);
    238218            psArrayAdd(job->args, 1, cells->data[j]); // sources
    239219            psArrayAdd(job->args, 1, smoothedImage);
    … …  
    243223            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS,    PS_TYPE_F32);
    244224            PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS,    PS_TYPE_S32);
    245             PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32);
    246             PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32);
    247225            PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA,  PS_TYPE_F32);
    248226            PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32);
    … …  
    281259    }
    282260    psFree (cellGroups);
    283     psFree (kronWindow);
    284261    if (KRON_SMOOTH) {
    285262        for (int i = 0; i < sources->n; i++) {
    … …  
    298275
    299276    pmReadout *readout              = job->args->data[0];
    300     psImage *kronWindow             = job->args->data[1];
    301     psArray *sources                = job->args->data[2];
    302     psImage *smoothedImage          = job->args->data[3];
    303     psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    304     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
    305     float RADIUS                    = PS_SCALAR_VALUE(job->args->data[6],F32);
    306     float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[7],F32);
    307     int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[8],S32);
    308     bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[9],S32);
    309     bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[10],S32);
    310     float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[11],F32);
    311     float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[12],F32);
    312     int pass                        = PS_SCALAR_VALUE(job->args->data[13],S32);
    313 #ifndef REVERT_ON_BAD_MEASUREMENT
    314     (void) pass;
    315 #endif
     277    psArray *sources                = job->args->data[1];
     278    psImage *smoothedImage          = job->args->data[2];
     279    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
     280    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
     281    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[5],F32);
     282    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
     283    int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[7],S32);
     284    float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[8],F32);
     285    float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[9],F32);
     286    int pass                        = PS_SCALAR_VALUE(job->args->data[10],S32);
     287
     288    bool measureRadius = true;
     289    if (pass == 3) {
     290        measureRadius = false;
     291    }
    316292
    317293    for (int j = 0; j < KRON_ITERATIONS; j++) {
    … …  
    320296            pmSource *source = sources->data[i];
    321297            if (!source->peak) continue; // XXX how can we have a peak-less source?
     298            if (!source->moments) continue;
     299            if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     300            if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     301            // skip saturated stars modeled with a radial profile
     302            if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
     303
     304            // in pass 3 we only measure the flux for matched sources
     305            if ((pass == 3) && !(source->mode2 & PM_SOURCE_MODE2_MATCHED)) continue;
    322306
    323307# if (0)
    324308# define TEST_X 653
    325309# define TEST_Y 466
    326         if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    327             fprintf (stderr, "test object\n");
    328         }
     310            if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     311                fprintf (stderr, "test object\n");
     312            }
    329313# undef TEST_X
    330314# undef TEST_Y
    331315# endif
    332316
    333             // check status of this source's moments
    334             if (!source->moments) continue;
    335             if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
    336             if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
     317            if (measureRadius) {
     318                // check status of this source's moments
     319                // XXX: I don't think that we have to apply these restrictions since we dropped the window function
     320                if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
     321                if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
     322            }
     323
    337324            if (!isfinite(source->moments->Mrf)) {
    338325                // Once we save a bad Mrf measurement we give up on this source
    339                 // checking here allows us to avoid adding and subtracting the model
     326                // XXX: is this the right thing to do?
    340327                continue;
    341328            }
    342 
    343             // skip saturated stars modeled with a radial profile
    344             if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
    345329
    346330            // replace object in image
    … …  
    358342            // On first iteration set window radius to sky radius (if valid). We also use this on subsequent
    359343            // iterations if we cannot find a better limit
    360             float maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
     344            float maxWindow;
     345            if (measureRadius) {
     346                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
     347            } else {
     348                maxWindow = source->moments->Mrf;
     349            }
    361350            if (j > 0) {
    362351                // on subsequent iterations we use a factor times the previous radial moment value
    … …  
    380369            float windowRadius = PS_MAX(RADIUS, maxWindow);
    381370
    382 #ifdef REVERT_ON_BAD_MEASURMENT
    383             // save previous measurements. We might revert back to them if this round fails
    384             float MrfPrior = source->moments->Mrf;
    385             float KronFluxPrior = source->moments->KronFlux;
    386             float KronFluxErrPrior = source->moments->KronFluxErr;
    387 #endif
    388 
    389371            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    390372            bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    391             psAssert (source->pixels, "WTF?");
     373            psAssert (source->pixels, "redefine pixels failed?");
    392374            if (extend && smoothedPixels) {
    393375                psFree(source->tmpPtr);
    394376                smoothedPixels = psImageSubset(smoothedImage, source->region);
    395                 psAssert (smoothedPixels, "WTF?");
     377                psAssert (smoothedPixels, "redefine smoothed pixels failed?");
    396378                source->tmpPtr = (psPtr) smoothedPixels ;
    397379            }
    398380
    399 
    400             // clear the window function for this source based on the moments
    401             // Note: this function also applies cuts on the source and returns false if it
    402             // does not meet the requirements for measuring the Kron Radius or Magnitude.
    403             // Note: that it performs the cuts even if KRON_APPLY_WINDOW is false
    404             if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) {
    405 
    406                 // this function populates moments->Mrf,KronFlux,KronFluxErr
    407                 psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, smoothedPixels);
    408                 psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    409 
    410                 // set a window function for each source based on the moments
    411                 psphotKronWindowSetSource (source, kronWindow, true, true, KRON_APPLY_WINDOW);
    412             }
    413 
    414 #ifdef REVERT_ON_BAD_MEASUREMENT
    415             // on pass 2 if we get an invalid measurement on a pass 1 source where we had a good one previously
    416             // in pass 1 keep that measurement
    417             bool reverted = false;
    418             if (pass > 1 && !isfinite(source->moments->Mrf) && (source->mode2 & PM_SOURCE_MODE2_PASS1_SRC)) {
    419                 source->moments->Mrf = MrfPrior;    // This is finite otherwise we wouldn't have gotten here
    420                 source->moments->KronFlux = KronFluxPrior;
    421                 source->moments->KronFluxErr = KronFluxErrPrior;
    422                 reverted = true;
    423             }
    424 #endif
     381            bool measureFlux = true;
     382            if (measureRadius) {
     383                measureFlux = psphotKronRadiusMeasure (source, windowRadius, MIN_KRON_RADIUS, maskVal, smoothedPixels);
     384            }
     385
     386            if (measureFlux) {
     387                // Make sure the sources images are large enough for the measured Kron Radius
     388                bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y,
     389                    2.5 * source->moments->Mrf + 2);
     390
     391                psAssert (source->pixels, "redefine pixels failed?");
     392                if (extend && smoothedPixels) {
     393                    psFree(source->tmpPtr);
     394                    smoothedPixels = psImageSubset(smoothedImage, source->region);
     395                    psAssert (smoothedPixels, "redefine smoothed pixels failed?");
     396                    source->tmpPtr = (psPtr) smoothedPixels ;
     397                }
     398
     399                // this function populates moments->Mrf,KronFlux,KronFluxErr, KronFinner, and KronFouter
     400                psphotKronFluxMeasure (source, maskVal);
     401            }
     402
     403            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     404
    425405#ifdef DUMP_KRS
    426 #ifndef REVERT_ON_BAD_MEASUREMENT
    427             bool reverted = false;
    428 #endif
    429             fprintf(dumpFile, "%7d %1d %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %2d\n", source->id, reverted, source->moments->Mrf, MrfPrior, maxWindow, windowRadius, source->peak->xf, source->peak->yf, source->imageID);
     406            fprintf(dumpFile, "%7d %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %2d\n", source->id, source->moments->Mrf, MrfPrior, maxWindow, windowRadius, source->peak->xf, source->peak->yf, source->imageID);
    430407#endif
    431408
    … …  
    443420}
    444421
    445 bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal,
    446     bool applyWeight, psImage *smoothedPixels) {
     422bool psphotKronRadiusMeasure (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal,
     423    psImage *smoothedPixels) {
    447424
    448425    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    452429
    453430    psF32 R2 = PS_SQR(radius);
    454     float rsigma2 =  applyWeight ? 0.5 / R2 : 0;
    455431
    456432    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    465441    // Xn  = SUM (x - xc)^n * (z - sky)
    466442
    467 
    468443    psF32 RF = 0.0;
    469444    psF32 RS = 0.0;
    … …  
    481456    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
    482457
    483     int Xwo = source->pixels->col0;
    484     int Ywo = source->pixels->row0;
    485 
    486458    psF32 **vPix;
    487459    if (smoothedPixels) {
    … …  
    490462        vPix = source->pixels->data.F32;
    491463    }
    492     psF32 **vWin = kronWindow->data.F32;
     464   
     465    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     466
     467    {
     468        for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     469
     470            psF32 yDiff = row - yCM;
     471            if (fabs(yDiff) > radius) continue;
     472
     473            // coordinate of mirror pixel
     474            int yFlip = yCM - yDiff;
     475            if (yFlip < 0) continue;
     476            if (yFlip >= source->pixels->numRows) continue;
     477
     478            for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     479                // check mask and value for this pixel
     480                if (vMsk && (vMsk[row][col] & maskVal)) continue;
     481                if (isnan(vPix[row][col])) continue;
     482
     483                psF32 xDiff = col - xCM;
     484                if (fabs(xDiff) > radius) continue;
     485
     486                // coordinate of mirror pixel
     487                int xFlip = xCM - xDiff;
     488                if (xFlip < 0) continue;
     489                if (xFlip >= source->pixels->numCols) continue;
     490
     491                // check mask and value for mirror pixel
     492                if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
     493                if (isnan(vPix[yFlip][xFlip])) continue;
     494
     495                // radius is just a function of (xDiff, yDiff)
     496                psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     497                if (r2 > R2) continue;
     498
     499                float fDiff1 = vPix[row][col];
     500                float fDiff2 = vPix[yFlip][xFlip];
     501
     502                float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
     503
     504                // Kron Flux uses the 1st radial moment
     505                psF32 rf = pDiff * sqrt(r2);
     506                psF32 rs = 0.5 * (fDiff1 + fDiff2);
     507
     508                RF  += rf;
     509                RS  += rs;
     510            }
     511        }
     512
     513        float MrfTry = RF/RS;
     514        if (RF <= 0. || RS <= 0 || !isfinite(MrfTry)) {
     515            // We did not get a good measurement
     516            source->moments->Mrf = NAN;
     517            source->moments->KronFlux  = NAN;
     518            source->moments->KronFluxErr  = NAN;
     519            return false;
     520        }
     521
     522        float Mrf = MAX(minKronRadius, MrfTry);
     523        // Saturate the 1st radial moment
     524        if (sqrt(source->peak->detValue) < 10.0) {
     525            Mrf = MIN (radius, Mrf);
     526        }
     527        source->moments->Mrf = Mrf;
     528    }
     529    return true;
     530}
     531
     532bool psphotKronFluxMeasure(pmSource *source, psImageMaskType maskVal) {
     533
     534    PS_ASSERT_PTR_NON_NULL(source, false);
     535    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     536    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
     537
     538    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
     539    // image coordinates.  the source->pixels image has an offset relative to its parent of
     540    // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in
     541    // this subimage.  we subtract off the peak coordinates, adjusted to this subimage, to have
     542    // minimal round-off error in the sums.  since these values are subtracted just to minimize
     543    // the dynamic range and are added back below, the exact value does not matter. these are
     544    // (int) so they can be used in the image index below.
     545
     546    // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
     547    // Xn  = SUM (x - xc)^n * (z - sky)
     548
     549    // the peak position is less accurate but less subject to extreme deviations
     550    float dX = source->moments->Mx - source->peak->xf;
     551    float dY = source->moments->My - source->peak->yf;
     552    float dR = hypot(dX, dY);
     553    float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
     554    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
     555
     556    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     557    // xCM, yCM from pixel coords to pixel index here.
     558    psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
     559    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
     560
     561    psF32 **vPix = source->pixels->data.F32;
    493562    psF32 **vWgt = source->variance->data.F32;
    494563   
    495564    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    496565
     566    source->moments->KronFlux  = NAN;
     567    source->moments->KronFluxErr  = NAN;
     568    source->moments->KronFinner = NAN;
     569    source->moments->KronFouter = NAN;
     570
     571    // Calculate the Kron fluxes
     572    float radKinner = 1.0 * source->moments->Mrf;
     573    float radKron   = 2.5 * source->moments->Mrf;
     574    float radKouter = 4.0 * source->moments->Mrf;
     575
     576    float limitRadius = MIN (radKouter, source->windowRadius);
     577    if (radKouter > source->windowRadius) {
     578        // This happens but the measurement isn't important enough to allocate the extra pixels
     579        // psWarning ("outer kron radius: %f is larger than windowRadius: %f for %d\n",
     580        //     radKouter, source->windowRadius, source->id);
     581        limitRadius = MIN (radKron, source->windowRadius);
     582    }
     583    if (radKron > source->windowRadius) {
     584        // caller should have prevented this from happening
     585        psWarning ("kron radius: %f is larger than windowRadius: %f for %d\n",
     586            radKron, source->windowRadius, source->id);
     587        return false;
     588    }
     589
     590    float Sum = 0.0;
     591    float Var = 0.0;
     592    float SumInner = 0.0;
     593    float SumOuter = 0.0;
     594
     595    // set vPix to the source pixels (it may have been set to the smoothed image above)
     596    vPix = source->pixels->data.F32;
     597
    497598    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    498599
    499600        psF32 yDiff = row - yCM;
    500         if (fabs(yDiff) > radius) continue;
    501 
    502         // coordinate of mirror pixel
    503         int yFlip = yCM - yDiff;
    504         if (yFlip < 0) continue;
    505         if (yFlip >= source->pixels->numRows) continue;
     601        if (fabs(yDiff) > limitRadius) continue;
    506602
    507603        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    … …  
    511607
    512608            psF32 xDiff = col - xCM;
    513             if (fabs(xDiff) > radius) continue;
    514 
    515             // coordinate of mirror pixel
    516             int xFlip = xCM - xDiff;
    517             if (xFlip < 0) continue;
    518             if (xFlip >= source->pixels->numCols) continue;
    519 
    520             // check mask and value for mirror pixel
    521             if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
    522             if (isnan(vPix[yFlip][xFlip])) continue;
    523 
    524             // radius is just a function of (xDiff, yDiff)
     609            if (fabs(xDiff) > limitRadius) continue;
     610
    525611            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    526             if (r2 > R2) continue;
    527 
    528             // flux * window
    529             float z = r2 * rsigma2;
    530             assert (z >= 0.0);
    531 
    532             // weight by window image and wide Gaussian
    533             float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
    534             float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
    535 
    536             float fDiff1 = vPix[row][col]*weight1;
    537             float fDiff2 = vPix[yFlip][xFlip]*weight2;
    538 
    539             float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
    540 
    541             // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
    542             psF32 rf = pDiff * sqrt(r2);
    543             psF32 rs = 0.5 * (fDiff1 + fDiff2);
    544 
    545             RF  += rf;
    546             RS  += rs;
     612            psF32 r = sqrt(r2);
     613
     614            float pDiff = vPix[row][col];
     615            psF32 wDiff = vWgt[row][col];
     616
     617            if (r > radKinner && r < radKron) {
     618                SumInner += pDiff;
     619            }
     620            if (r < radKron) {
     621                Sum += pDiff;
     622                Var += wDiff;
     623            }
     624            if (r > radKron && r < radKouter) {
     625                SumOuter += pDiff;
     626            }
    547627        }
    548628    }
    549629
    550     float MrfTry = RF/RS;
    551     if (RF <= 0. || RS <= 0 || !isfinite(MrfTry)) {
    552         // We did not get a good measurement
    553         source->moments->Mrf = NAN;
    554         source->moments->KronFlux  = NAN;
    555         source->moments->KronFluxErr  = NAN;
    556         return false;
    557     }
    558 
    559     float Mrf = MAX(minKronRadius, MrfTry);
    560     // Saturate the 1st radial moment
    561     if (sqrt(source->peak->detValue) < 10.0) {
    562         Mrf = MIN (radius, Mrf);
    563     }
    564 
    565     // Calculate the Kron magnitude (make this block optional?)
    566     float radKron  = 2.5*Mrf;
    567     float radKron2 = radKron*radKron;
    568 
    569     int nKronPix = 0;
    570     float Sum = 0.0;
    571     float Var = 0.0;
    572     float Win = 0.0;
    573 
    574     // set vPix to the source pixels (it may have been set to the
    575     // smoothed image above)
    576     vPix = source->pixels->data.F32;
    577 
    578 
    579     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    580 
    581         psF32 yDiff = row - yCM;
    582         if (fabs(yDiff) > radKron) continue;
    583 
    584         for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    585             // check mask and value for this pixel
    586             if (vMsk && (vMsk[row][col] & maskVal)) continue;
    587             if (isnan(vPix[row][col])) continue;
    588 
    589             psF32 xDiff = col - xCM;
    590             if (fabs(xDiff) > radKron) continue;
    591 
    592             // radKron is just a function of (xDiff, yDiff)
    593             psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    594             if (r2 > radKron2) continue;
    595 
    596             float weight1  = vWin[row+Ywo][col+Xwo];
    597             float fDiff1 = vPix[row][col]*weight1;
    598 
    599             float pDiff = fDiff1;
    600             psF32 wDiff = vWgt[row][col] * weight1;
    601 
    602             Sum += pDiff;
    603             Var += wDiff;
    604             Win += weight1;
    605             nKronPix ++;
    606         }
    607     }
    608 
    609     source->moments->Mrf         = Mrf;
    610630    source->moments->KronFlux    = Sum;
    611631    source->moments->KronFluxErr = sqrt(Var);
     632    source->moments->KronFinner  = SumInner;
     633
     634    // only save radKouter if the radius is inside the integration radius limit
     635    if (radKouter <= limitRadius) {
     636        source->moments->KronFouter = SumOuter;
     637    }
    612638
    613639    return true;
    614640}
    615 
    616 bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) {
    617 
    618     if (!source) return false;
    619     if (!source->peak) return false; // XXX how can we have a peak-less source?
    620     if (!source->moments) return false;
    621     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
    622     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
    623     if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) return false;
    624     if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false;
    625     psAssert(kronWindow, "need a window");
    626 
    627     // XXX: If we are not applying the window then we don't need to check for valid Mrf here.
    628     // We should give the this module a chance to measure a good value.
    629     // However experiments show that it hardly ever succeeds in getting a better value
    630     if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
    631 
    632     if (!applyWindow) {
    633         return true;
    634     }
    635 
    636     // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has
    637     // these values (and peak of 1.0)
    638 
    639     int Nx = kronWindow->numCols;
    640     int Ny = kronWindow->numRows;
    641 
    642     float Xo = source->moments->Mx;
    643     float Yo = source->moments->My;
    644 
    645     psEllipseMoments moments;
    646     moments.x2 = source->moments->Mxx;
    647     moments.y2 = source->moments->Myy;
    648     moments.xy = source->moments->Mxy;
    649 
    650     psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
    651     if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) {
    652         // Shall we log a proper warning? This happens often with matched sources (forced photometry)
    653         // fprintf(dump, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta);
    654         return false;
    655     }
    656 
    657     // Why this factor of 0.5 ?
    658     float scale = 0.5 * source->moments->Mrf / axes.major;
    659     axes.major *= scale;
    660     axes.minor *= scale;
    661 
    662     psEllipseShape shape = psEllipseAxesToShape(axes);
    663     if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) {
    664         // Shall we log a proper warning? This happens often with matched sources (forced photometry)
    665         // fprintf(dump, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy);
    666         return false;
    667     }
    668 
    669     float Smajor = axes.major;
    670 
    671     int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1);
    672     int maxX = PS_MIN(PS_MAX(Xo + 5*Smajor, 0), Nx - 1);
    673     int minY = PS_MIN(PS_MAX(Yo - 5*Smajor, 0), Ny - 1);
    674     int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1);
    675 
    676     float rMxx = 0.5 / PS_SQR(shape.sx);
    677     float rMyy = 0.5 / PS_SQR(shape.sy);
    678     float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
    679                                     // implementation. XXX: Is this correct?
    680 
    681     for (int iy = minY; iy < maxY; iy++) {
    682         for (int ix = minX; ix < maxX; ix++) {
    683 
    684             float dX = (ix + 0.5 - Xo);
    685             float dY = (iy + 0.5 - Yo);
    686 
    687             float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
    688 
    689             float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z));
    690 
    691             kronWindow->data.F32[iy][ix] *= f;
    692         }
    693     }
    694 
    695     return true;
    696 }

  • branches/eam_branches/ipp-20120905/psphot/src/psphotMaskBackground.c

    r31154 r34570  
    1717
    1818    // find the currently selected readout (XXX note that the model is saved on PSPHOT.BACKMDL regardless of 'filename'
    19     pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest
     19    pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, psphotGetFilerule("PSPHOT.BACKMDL"), index); // File of interest
    2020    assert (modelFile);
    2121
    22     pmFPAfile *stdevFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL.STDEV", index);
     22    pmFPAfile *stdevFile = pmFPAfileSelectSingle(config->files, psphotGetFilerule("PSPHOT.BACKMDL.STDEV"), index);
    2323    assert (stdevFile);
    2424

  • branches/eam_branches/ipp-20120905/psphot/src/psphotModelBackground.c

    r33963 r34570  
    405405
    406406    psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters
    407     pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning, index);
    408     pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning, index);
     407    pmReadout *model = pmFPAGenerateReadout(config, view, psphotGetFilerule("PSPHOT.BACKMDL"), inFPA, binning, index);
     408    pmReadout *modelStdev = pmFPAGenerateReadout(config, view, psphotGetFilerule("PSPHOT.BACKMDL.STDEV"), inFPA, binning, index);
    409409
    410410    if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config, false)) {

  • branches/eam_branches/ipp-20120905/psphot/src/psphotOutput.c

    r32695 r34570  
    3434   
    3535
    36     pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL");
     36    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, psphotGetFilerule("PSPHOT.BACKMDL"));
    3737    if (!file) return NULL;
    3838
    … …  
    4545    bool status;
    4646
    47     pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL.STDEV");
     47    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, psphotGetFilerule("PSPHOT.BACKMDL.STDEV"));
    4848    if (!file) return NULL;
    4949

  • branches/eam_branches/ipp-20120905/psphot/src/psphotReplaceUnfit.c

    r34136 r34570  
    4343}
    4444
     45// the return state indicates if any sources were actually replaced
    4546bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) {
    4647
    … …  
    5859    if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) {
    5960        psLogMsg ("psphot", PS_LOG_DETAIL, "skipping replace all sources for input file %d", index);
    60         return true;
     61        return false;
    6162    }
    6263
    … …  
    6566
    6667    psArray *sources = detections->allSources;
    67     //psAssert (sources, "missing sources?");
    68     if (!sources) return true;
     68    // if no work to do, should just return false
     69    if (!sources) return false;
    6970
    7071    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)

  • branches/eam_branches/ipp-20120905/psphot/src/psphotSetThreads.c

    r34317 r34570  
    3030    psFree(task);
    3131
    32     task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 14);
     32    task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 11);
    3333    task->function = &psphotKronIterate_Threaded;
    3434    psThreadTaskAdd(task);

  • branches/eam_branches/ipp-20120905/psphot/src/psphotSkyReplace.c

    r31673 r34570  
    3636
    3737    // select background pixels, from output background file, or create
    38     pmReadout *background = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKGND");
     38    pmReadout *background = pmFPAfileThisReadout (config->files, view, psphotGetFilerule("PSPHOT.BACKGND"));
    3939    if (background == NULL) psAbort("background not defined");
    4040

  • branches/eam_branches/ipp-20120905/psphot/src/psphotSourceMatch.c

    r34354 r34570  
    559559}
    560560
    561 bool psphotDropBadMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
     561bool psphotFilterMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
    562562
    563563    bool status = false;
    564564
    565     psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Drop Bad Matched Sources ---");
     565    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Filter Matched Sources ---");
    566566
    567567    psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    … …  
    615615    psFree(dropped);
    616616
     617    // find the "best" Mrf from the detected sources.
     618    // Currently we use the smallest positive value
     619    for (int i=0; i< objects->n; i++) {
     620        pmPhotObj *obj = objects->data[i];
     621
     622        float minMrf = 1000.;
     623        bool hasMatched = false;
     624        for (int j = 0; j < obj->sources->n; j++) {
     625            pmSource *source = obj->sources->data[j];
     626            if (source->mode2 & PM_SOURCE_MODE2_MATCHED) {
     627                hasMatched = true;
     628                continue;
     629            }
     630            float Mrf = source->moments->Mrf;
     631            if (isfinite(Mrf) && Mrf < minMrf && Mrf > 0) {
     632                minMrf = Mrf;
     633            }
     634        }
     635
     636        if (!hasMatched || minMrf > 120.) {
     637            continue;
     638        }
     639
     640        // set Mrf for matched sources to the value found above
     641        for (int j = 0; j < obj->sources->n; j++) {
     642            pmSource *source = obj->sources->data[j];
     643            if (source->mode2 & PM_SOURCE_MODE2_MATCHED) {
     644                source->moments->Mrf = minMrf;
     645            }
     646        }
     647    }
     648
    617649    return true;
    618650}

  • branches/eam_branches/ipp-20120905/psphot/src/psphotStack.c

    r34283 r34570  
    22
    33int main (int argc, char **argv) {
     4
     5    // Set flag to tell certain library functions that we are in the psphotStack program.
     6    // (This is defined in psphotDefineFiles.c)
     7    psphotINpsphotStack = true;
    48
    59    // uncomment to turn on memory dumps (move this to an option)

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

    r34466 r34570  
    4343    if (!needConvolved) {
    4444        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.INPUT.CNV");
     45        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.MASK.CNV");
     46        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.VARIANCE.CNV");
     47        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.PSF.CNV");
    4548    }
    4649
    … …  
    97100            // drop all versions of the internal files
    98101            status = true;
    99             status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
    100             status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV");
    101             status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND");
     102            status &= pmFPAfileDropInternal (config->files, psphotGetFilerule("PSPHOT.BACKMDL"));
     103            status &= pmFPAfileDropInternal (config->files, psphotGetFilerule("PSPHOT.BACKMDL.STDEV"));
     104            status &= pmFPAfileDropInternal (config->files, psphotGetFilerule("PSPHOT.BACKGND"));
    102105            if (!status) {
    103106                psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files");

  • branches/eam_branches/ipp-20120905/psphot/src/psphotStackParseCamera.c

    r33913 r34570  
    1414        return false;
    1515    }
     16
     17    // select the appropriate recipe information
     18    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     19    bool savePSF = psMetadataLookupBool(&status, recipe, "SAVE.PSF");
     20    bool saveBackgroundModel = psMetadataLookupBool(&status, recipe, "SAVE.BACKMDL");
     21    bool saveBackground = psMetadataLookupBool(&status, recipe, "SAVE.BACKGND");
     22    bool saveBackSub = psMetadataLookupBool(&status, recipe, "SAVE.BACKSUB");
     23    bool saveResid = psMetadataLookupBool(&status, recipe, "SAVE.RESID");
     24
     25    bool saveCnv = psMetadataLookupBool(&status, recipe, "SAVE.CNV");
     26    bool saveChisq = psMetadataLookupBool(&status, recipe, "SAVE.CHISQ");
     27    bool useRaw = psMetadataLookupBool(&status, recipe, "PSPHOT.STACK.USE.RAW");
     28
    1629
    1730    int nRaw = 0;
    … …  
    120133            }
    121134        }
     135        pmFPAfile *inputTemplate;
     136        if (useRaw) {
     137            inputTemplate = rawInputFile;
     138        } else {
     139            inputTemplate = cnvInputFile;
     140        }
     141        if (!inputTemplate) {
     142            psError(PS_ERR_UNKNOWN, true, "cannot determinte inputTemplate: USE.RAW: %d\n", useRaw);
     143            return false;
     144        }
    122145
    123146        psString sources = psMetadataLookupStr(&status, input, "SOURCES"); // Name of mask
    … …  
    135158        }
    136159        // generate an pmFPAimage for the output convolved image
    137         // XXX output of these files should be optional
    138160        {
    139161            // pmFPAfile *srcInputFile = rawInputFile ? rawInputFile : cnvInputFile;
    … …  
    143165                return false;
    144166            }
    145             outputImage->save = true;
     167            outputImage->save = saveCnv;
    146168            outputImage->fileID = stack_id;             // this is used to generate output names
    147169
    … …  
    155177                return NULL;
    156178            }
    157             outputMask->save = true;
     179            outputMask->save = saveCnv;
    158180            outputMask->fileID = stack_id;              // this is used to generate output names
    159181
    … …  
    167189                return NULL;
    168190            }
    169             outputVariance->save = true;
     191            outputVariance->save = saveCnv;
    170192            outputVariance->fileID = stack_id;          // this is used to generate output names
    171193
    … …  
    178200            outsources->save = true;
    179201            outsources->fileID = stack_id;              // this is used to generate output names
     202
     203            if (savePSF) {
     204                pmFPAfile *output = pmFPAfileDefineOutputFromFile (config, inputTemplate, "PSPHOT.STACK.PSF.SAVE");
     205                if (!output) {
     206                    psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.PSF.SAVE");
     207                    return false;
     208                }
     209                output->save = true;
     210                output->fileID = stack_id;
     211            }
     212            if (saveBackgroundModel) {
     213                int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN");
     214                int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN");
     215                pmFPAfile *output = pmFPAfileDefineFromFile (config, inputTemplate, DX, DY, "PSPHOT.STACK.BACKMDL");
     216                if (!output) {
     217                    psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.BACKMDL");
     218                    return false;
     219                }
     220
     221                output->save = true;
     222                output->fileID = stack_id;
     223            }
     224            if (saveBackground) {
     225                pmFPAfile *output = pmFPAfileDefineFromFile (config, inputTemplate, 1, 1, "PSPHOT.STACK.BACKGND");
     226                if (!output) {
     227                    psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.STACK.BACKGND");
     228                    return false;
     229                }
     230                output->save = true;
     231                output->fileID = stack_id;
     232            }
     233            if (saveBackSub) {
     234                pmFPAfile *output = pmFPAfileDefineFromFile (config, inputTemplate, 1, 1, "PSPHOT.STACK.BACKSUB");
     235                if (!output) {
     236                    psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.STACK.BACKSUB");
     237                    return false;
     238                }
     239                output->save = true;
     240                output->fileID = stack_id;
     241            }
     242            if (saveResid) {
     243                pmFPAfile *output = pmFPAfileDefineOutputFromFile (config, inputTemplate, "PSPHOT.STACK.RESID");
     244                if (!output) {
     245                    psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.RESID");
     246                    return false;
     247                }
     248                output->save = true;
     249                output->fileID = stack_id;
     250            }
    180251        }
    181252    }
    … …  
    191262    }
    192263
    193     // select the appropriate recipe information
    194     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    195     bool saveCnv = psMetadataLookupBool(&status, recipe, "SAVE.CNV");
    196     bool saveChisq = psMetadataLookupBool(&status, recipe, "SAVE.CHISQ");
    197 
    198     // loop over the available readouts
    199     for (int i = 0; i < nInputs; i++) {
    200         pmFPAfile *file = NULL;
    201 
    202         file = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i);
    203         file->save = saveCnv;
    204 
    205         file = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.MASK", i);
    206         file->save = saveCnv;
    207 
    208         file = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.VARIANCE", i);
    209         file->save = saveCnv;
    210     }
    211 
    212264    // generate an pmFPAimage for the chisqImage
    213     // XXX output of these files should be optional
    214265    {
    215266        pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE");

  • branches/eam_branches/ipp-20120905/psphot/src/psphotStackReadout.c

    r34416 r34570  
    369369    psMemDump("psfstats");
    370370
    371     // drop matched sources without any useful measurements
    372     psphotDropBadMatchedSources (config, view, STACK_SRC, objects);
     371    // drop matched sources without any useful measurements and set kron radii for the ones
     372    // we decide to keep
     373    psphotFilterMatchedSources (config, view, STACK_SRC, objects);
     374
     375    // measure kron fluxes for the matched sources only
     376    psphotKronIterate(config, view, STACK_SRC, 3);
    373377
    374378    // measure elliptical apertures, petrosians (objects sorted by S/N)

  • branches/eam_branches/ipp-20120905/psphot/src/psphotSubtractBackground.c

    r32348 r34570  
    2121
    2222    // find the currently selected readout (XXX note that the model is saved on PSPHOT.BACKMDL regardless of 'filename'
    23     pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest
     23    pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, psphotGetFilerule("PSPHOT.BACKMDL"), index); // File of interest
    2424    assert (modelFile);
    2525
    … …  
    3636    // select background pixels, from output background file, or create
    3737    // XXX for now, we will only allow a single background image to be generated
    38     file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND");
     38    file = psMetadataLookupPtr (&status, config->files, psphotGetFilerule("PSPHOT.BACKGND"));
    3939    if (file) {
    4040        // we are using PSPHOT.BACKGND as an I/O file: select readout or create
    … …  
    5050        }
    5151    } else {
    52         background = pmFPAfileDefineInternal (config->files, "PSPHOT.BACKGND", image->numCols, image->numRows, PS_TYPE_F32);
     52        background = pmFPAfileDefineInternal (config->files, psphotGetFilerule("PSPHOT.BACKGND"), image->numCols, image->numRows, PS_TYPE_F32);
    5353    }
    5454    psF32 **backData = background->image->data.F32;
Note: See TracChangeset for help on using the changeset viewer.