IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34247


Ignore:
Timestamp:
Jul 31, 2012, 11:50:18 AM (14 years ago)
Author:
eugene
Message:

merge changes from trunk

Location:
branches/eam_branches/ipp-20120627/psphot
Files:
23 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20120627/psphot

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

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

    r34172 r34247  
    102102bool            psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    103103
    104 bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final);
    105 bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode);
     104bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final, bool skipNegativeFluxSources);
     105bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode, bool skipNegativeFluxSources);
    106106
    107107bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize);
    … …  
    265265bool            psphotEllipticalProfile (pmSource *source, bool RAW_RADIUS);
    266266bool            psphotEllipticalContour (pmSource *source);
     267bool            psphotLimitRadialApertures(psMetadata *recipe, long nSources);
    267268
    268269// psphotVisual functions
    … …  
    423424// bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index);
    424425
     426bool psphotStackSetInputsToSkip(pmConfig *config, const pmFPAview *view, const char *filerule, bool set) ;
     427
    425428bool psphotStackRenormaliseVariance(const pmConfig *config, pmReadout *readout);
    426429

  • branches/eam_branches/ipp-20120627/psphot/src/psphotChoosePSF.c

    r32348 r34247  
    4747    if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) {
    4848        psLogMsg ("psphot", PS_LOG_DETAIL, "psf model supplied for input file %d", index);
     49        return true;
     50    }
     51
     52    if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) {
     53        psLogMsg ("psphot", PS_LOG_DETAIL, "skipping choose PSF for input file %d", index);
    4954        return true;
    5055    }

  • branches/eam_branches/ipp-20120627/psphot/src/psphotCullPeaks.c

    r31154 r34247  
    6363# if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE)
    6464    psLogMsg ("psphot", PS_LOG_INFO, "Culling peaks from footprints using the smoothed image");
     65    bool setNaNsToZero = psMetadataLookupF32 (&status, recipe, "FOOTPRINT_SET_NANS_TO_ZERO");
     66    if (setNaNsToZero) {
     67        // Set NaN pixels to zero so that they do not considered in the footprint analysis
     68        // XXX: Currently the caller does nothing further with the significance readout so
     69        // modifiying the image does no harm
     70        // XXX: Note: signifR->weight is not used by pmFootprintCullPeaks so we don't
     71        // need to touch it
     72        psImageClipNaN(signifR->image, 0);
     73    }
     74       
     75
    6576# else
    6677    psLogMsg ("psphot", PS_LOG_INFO, "Culling peaks from footprints using the raw (unsmoothed) image");
    … …  
    7384
    7485        if (fp->npix > 30000) {
    75             fprintf (stderr, "big footprint: %f %f to %f %f (%d pix)\n", fp->bbox.x0, fp->bbox.y0, fp->bbox.x1, fp->bbox.y1, fp->npix);
     86            psLogMsg("psphot.cull.footprints", PS_LOG_WARN, "big footprint: %f %f to %f %f (%d pix)\n", fp->bbox.x0, fp->bbox.y0, fp->bbox.x1, fp->bbox.y1, fp->npix);
    7687        }
    7788        psTimerStart ("psphot.cull.footprints");
    … …  
    97108# endif
    98109
     110
    99111        float dtime = psTimerMark ("psphot.cull.footprints");
    100112        if (dtime > 1.0) {
    101             fprintf (stderr, "slow cull for %d (%f sec)\n", i, dtime);
     113            psLogMsg("psphot.cull.footprints", PS_LOG_WARN, "slow cull for %d (%f sec)\n", i, dtime);
    102114        }
    103115    }

  • branches/eam_branches/ipp-20120627/psphot/src/psphotEfficiency.c

    r34086 r34247  
    420420
    421421    // psphotFitSourcesLinearReadout subtracts the model fits
    422     if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) {
     422    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST, false)) {
    423423        psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
    424424        psFree(fakeSources);

  • branches/eam_branches/ipp-20120627/psphot/src/psphotExtendedSourceFits.c

    r34190 r34247  
    175175            psArrayAdd(job->args, 1, cells->data[j]); // sources
    176176            psArrayAdd(job->args, 1, models);
     177            // Allocate a metadata iterator here because psMetadataIteratorAlloc/Free are not thread safe
     178            psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
     179            psArrayAdd(job->args, 1, iter);
    177180            psArrayAdd(job->args, 1, AnalysisRegion); // XXX make a pointer
    178181
    … …  
    203206            }
    204207            psScalar *scalar = NULL;
    205             scalar = job->args->data[7];
     208            scalar = job->args->data[8];
    206209            Next += scalar->data.S32;
    207             scalar = job->args->data[8];
     210            scalar = job->args->data[9];
    208211            Nconvolve += scalar->data.S32;
    209             scalar = job->args->data[9];
     212            scalar = job->args->data[10];
    210213            NconvolvePass += scalar->data.S32;
    211             scalar = job->args->data[10];
     214            scalar = job->args->data[11];
    212215            Nplain += scalar->data.S32;
    213             scalar = job->args->data[11];
     216            scalar = job->args->data[12];
    214217            NplainPass += scalar->data.S32;
    215             scalar = job->args->data[12];
     218            scalar = job->args->data[13];
    216219            Nfaint += scalar->data.S32;
    217             scalar = job->args->data[13];
     220            scalar = job->args->data[14];
    218221            Nfail += scalar->data.S32;
     222            psFree(job->args->data[3]); // iterator allocated above
    219223            psFree(job);
    220224# endif
    … …  
    235239            } else {
    236240                psScalar *scalar = NULL;
    237                 scalar = job->args->data[7];
     241                scalar = job->args->data[8];
    238242                Next += scalar->data.S32;
    239                 scalar = job->args->data[8];
     243                scalar = job->args->data[9];
    240244                Nconvolve += scalar->data.S32;
    241                 scalar = job->args->data[9];
     245                scalar = job->args->data[10];
    242246                NconvolvePass += scalar->data.S32;
    243                 scalar = job->args->data[10];
     247                scalar = job->args->data[11];
    244248                Nplain += scalar->data.S32;
    245                 scalar = job->args->data[11];
     249                scalar = job->args->data[12];
    246250                NplainPass += scalar->data.S32;
    247                 scalar = job->args->data[12];
     251                scalar = job->args->data[13];
    248252                Nfaint += scalar->data.S32;
    249                 scalar = job->args->data[13];
     253                scalar = job->args->data[14];
    250254                Nfail += scalar->data.S32;
     255                psFree(job->args->data[3]); // metadata iterator allocated above
    251256            }
    252257            psFree(job);
    … …  
    285290    psArray *sources        = job->args->data[1];
    286291    psMetadata *models      = job->args->data[2];
    287     psRegion *region        = job->args->data[3];
    288     int psfSize             = PS_SCALAR_VALUE(job->args->data[4],S32);
    289     psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
    290     psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
     292    psMetadataIterator *iter = job->args->data[3];
     293    psRegion *region        = job->args->data[4];
     294    int psfSize             = PS_SCALAR_VALUE(job->args->data[5],S32);
     295    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
     296    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
    291297
    292298    // Define source fitting parameters for extended source fits
    … …  
    338344
    339345        // loop here over the models chosen for each source (exclude by S/N)
    340         psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
     346        // Reset the iterator
     347        psMetadataIteratorSet(iter, PS_LIST_HEAD);
    341348        psMetadataItem *item = NULL;
    342349        while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
    … …  
    428435          psFree (modelFit);
    429436        }
    430         psFree (iter);
    431437
    432438        // evaluate the relative quality of the models, choose one
    … …  
    489495
    490496    // change the value of a scalar on the array (wrap this and put it in psArray.h)
    491     scalar = job->args->data[7];
     497    scalar = job->args->data[8];
    492498    scalar->data.S32 = Next;
    493499
    494     scalar = job->args->data[8];
     500    scalar = job->args->data[9];
    495501    scalar->data.S32 = Nconvolve;
    496502
    497     scalar = job->args->data[9];
     503    scalar = job->args->data[10];
    498504    scalar->data.S32 = NconvolvePass;
    499505
    500     scalar = job->args->data[10];
     506    scalar = job->args->data[11];
    501507    scalar->data.S32 = Nplain;
    502508
    503     scalar = job->args->data[11];
     509    scalar = job->args->data[12];
    504510    scalar->data.S32 = NplainPass;
    505511
    506     scalar = job->args->data[12];
     512    scalar = job->args->data[13];
    507513    scalar->data.S32 = Nfaint;
    508514
    509     scalar = job->args->data[13];
     515    scalar = job->args->data[14];
    510516    scalar->data.S32 = Nfail;
    511517

  • branches/eam_branches/ipp-20120627/psphot/src/psphotFitSourcesLinear.c

    r34086 r34247  
    1313
    1414// for now, let's store the detections on the readout->analysis for each readout
    15 bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final)
     15bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final, bool skipNegativeFluxSources)
    1616{
    1717    bool status = true;
    … …  
    5151        psAssert (readout, "missing readout?");
    5252
     53        if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) {
     54            psLogMsg ("psphot", PS_LOG_DETAIL, "skipping fit sources for input file %d", i);
     55            continue;
     56        }
     57
     58
    5359        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    5460        psAssert (detections, "missing detections?");
    … …  
    6066        psAssert (psf, "missing psf?");
    6167
    62         if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) {
     68        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1, skipNegativeFluxSources)) {
    6369            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
    6470            return false;
    … …  
    7985
    8086            // rerun fit with correct fitVarMode
    81             if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarMode)) {
     87            // XXX: does skipNegativeFlux work here?
     88            if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarMode, skipNegativeFluxSources)) {
    8289                psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
    8390                return false;
    … …  
    130137}
    131138
    132 bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) {
     139bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode, bool skipNegativeFluxSources) {
    133140    bool status;
    134141    float x;
    … …  
    177184    float MIN_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MIN_VALID_FLUX");
    178185    if (!status) {
     186        MIN_VALID_FLUX = 0.0;
     187    }
     188    if (skipNegativeFluxSources && MIN_VALID_FLUX < 0) {
    179189        MIN_VALID_FLUX = 0.0;
    180190    }
    … …  
    252262        if (modelSum < 0.5) continue; // skip sources with no model constraint (somewhat arbitrary limit)
    253263        if (modelSum < 0.8) {
    254             fprintf (stderr, "low-sig model @ %f, %f (%f masked sum, %f sum, %f peak)\n",
     264            psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "low-sig model @ %f, %f (%f masked sum, %f sum, %f peak)\n",
    255265                     source->peak->xf, source->peak->yf, maskedSum, modelSum, source->peak->rawFlux);
    256266        }
    257267        if (maskedSum < 0.5) {
    258             fprintf (stderr, "worrying model @ %f, %f (%f masked sum, %f sum, %f peak)\n",
     268            psLogMsg ("psphot.ensemble",  PS_LOG_MINUTIA, "worrying model @ %f, %f (%f masked sum, %f sum, %f peak)\n",
    259269                     source->peak->xf, source->peak->yf, maskedSum, modelSum, source->peak->rawFlux);
    260270        }

  • branches/eam_branches/ipp-20120627/psphot/src/psphotForcedReadout.c

    r33140 r34247  
    6464
    6565    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    66     psphotFitSourcesLinear (config, view, filerule, false);
     66    psphotFitSourcesLinear (config, view, filerule, false, false);
    6767
    6868    // identify CRs and extended sources

  • branches/eam_branches/ipp-20120627/psphot/src/psphotRadialApertures.c

    r33089 r34247  
    7171    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    7272    psAssert (readout, "missing readout?");
     73   
     74    if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) {
     75        psLogMsg ("psphot", PS_LOG_DETAIL, "skipping radial aptertures for input file %d", index);
     76        return true;
     77    }
     78
    7379
    7480    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");

  • branches/eam_branches/ipp-20120627/psphot/src/psphotRadialBins.c

    r28013 r34247  
    202202}
    203203
     204// If the number of sources is large, edit the recipe to remove radial bins beyond
     205// a value specified in the recipe
     206bool psphotLimitRadialApertures(psMetadata *recipe, long nSources) {
     207
     208    bool status = false;
     209    long sourceLimit = psMetadataLookupS32 (&status, recipe, "RADIAL.NUM.SOURCES.LIMIT");
     210    if (!status) {
     211        sourceLimit = 50000;
     212    }
     213    if (nSources < sourceLimit) {
     214        return true;
     215    }
     216    psF32 maxRadius = psMetadataLookupF32 (&status, recipe, "RADIAL.SOURCES.OVER.LIMIT.RADIUS");
     217    if (!status) {
     218        maxRadius = 20.0;
     219    }
     220    psLogMsg ("psphot", PS_LOG_INFO, "Number of objects: %ld is greater than limit %ld. Limiting radial annular bins to %.3f pixels\n",
     221            nSources, sourceLimit, maxRadius);
     222
     223    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
     224    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
     225    if (!radMin || !radMin->n) {
     226        psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMin missing or empty)");
     227        return false;
     228    }
     229    if (!radMax || !radMax->n) {
     230        psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMax missing or empty)");
     231        return false;
     232    }
     233    if (radMax->n != radMin->n) {
     234        psError (PSPHOT_ERR_CONFIG, true, "length of radMin %ld and radMax %ld not equal)", radMin->n, radMax->n);
     235        return false;
     236    }
     237    int i = 0;
     238    psF32 lastRadius = 0;
     239    for (; i < radMax->n; i++) {
     240        if (radMax->data.F32[i] > maxRadius) {
     241            break;
     242        }
     243        lastRadius = radMax->data.F32[i];
     244    }
     245    if (i == radMax->n) {
     246        psLogMsg ("psphot", PS_LOG_INFO, "Radius of all bins is within the limit. lastRadius: %.3f\n", lastRadius);
     247        return true;
     248    }
     249    psLogMsg ("psphot", PS_LOG_INFO, "  radius limit exceeded at bin %d of %ld. New max radius is %.3f\n",
     250            i, radMax->n, lastRadius);
     251
     252    psVector *radMinNew = psVectorRealloc(radMin, i);
     253    if (radMinNew != radMin) {
     254        psMetadataAddVector(recipe, PS_LIST_TAIL, "RADIAL.ANNULAR.BINS.LOWER", PS_META_REPLACE, "", radMinNew);
     255        // XXX: I don't need this so I?
     256        // psFree(radMinNew);
     257    }
     258
     259    psVector *radMaxNew = psVectorRealloc(radMax, i);
     260    if (radMaxNew != radMax) {
     261        psMetadataAddVector(recipe, PS_LIST_TAIL, "RADIAL.ANNULAR.BINS.UPPER", PS_META_REPLACE, "", radMaxNew);
     262        // XXX: I don't need to free this do I?
     263        // psFree(radMaxNew);
     264    }
     265
     266    return true;
     267}
     268
    204269// the area-weighted mean radius is given by:
    205270

  • branches/eam_branches/ipp-20120627/psphot/src/psphotReadout.c

    r34086 r34247  
    185185
    186186    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    187     psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources)
     187    psphotFitSourcesLinear (config, view, filerule, false, false); // pass 1 (detections->allSources)
    188188
    189189    // measure the radial profiles to the sky
    … …  
    212212    // linear fit to include all sources (subtract again)
    213213    // NOTE : apply to ALL sources (extended + psf)
    214     psphotFitSourcesLinear (config, view, filerule, true); // pass 2 (detections->allSources)
     214    psphotFitSourcesLinear (config, view, filerule, true, false); // pass 2 (detections->allSources)
    215215
    216216    // if we only do one pass, skip to extended source analysis
    … …  
    261261
    262262        // NOTE: apply to ALL sources
    263         psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources)
     263        psphotFitSourcesLinear (config, view, filerule, true, false); // pass 3 (detections->allSources)
    264264    }
    265265
    … …  
    302302
    303303        // NOTE: apply to ALL sources
    304         psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources)
     304        psphotFitSourcesLinear (config, view, filerule, true, false); // pass 3 (detections->allSources)
    305305    }
    306306

  • branches/eam_branches/ipp-20120627/psphot/src/psphotReadoutForcedKnownSources.c

    r32996 r34247  
    4545
    4646    // linear PSF fit to source peaks
    47     psphotFitSourcesLinear (config, view, filerule, false);
     47    psphotFitSourcesLinear (config, view, filerule, false, false);
    4848
    4949    // calculate source magnitudes

  • branches/eam_branches/ipp-20120627/psphot/src/psphotReadoutKnownSources.c

    r33692 r34247  
    5757
    5858    // linear PSF fit to source peaks
    59     psphotFitSourcesLinear (config, view, filerule, false);
     59    psphotFitSourcesLinear (config, view, filerule, false, false);
    6060
    6161    // measure aperture photometry corrections

  • branches/eam_branches/ipp-20120627/psphot/src/psphotReadoutMinimal.c

    r34190 r34247  
    8181
    8282    // linear PSF fit to source peaks
    83     psphotFitSourcesLinear (config, view, filerule, false);
     83    psphotFitSourcesLinear (config, view, filerule, false, false);
    8484
    8585    // measure the radial profiles to the sky

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

    r33980 r34247  
    5555    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    5656    psAssert (readout, "missing readout?");
     57
     58    if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) {
     59        psLogMsg ("psphot", PS_LOG_DETAIL, "skipping replace all sources for input file %d", index);
     60        return true;
     61    }
    5762
    5863    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    … …  
    304309    psAssert (readout, "missing readout?");
    305310
     311    if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) {
     312        psLogMsg ("psphot", PS_LOG_DETAIL, "skipping reset models for input file %d", index);
     313        return true;
     314    }
     315
     316
    306317    // XXX the sources have already been copied (merge into here?)
    307318    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");

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

    r33964 r34247  
    4040    psFree(task);
    4141
    42     task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 14);
     42    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 15);
    4343    task->function = &psphotExtendedSourceFits_Threaded;
    4444    psThreadTaskAdd(task);

  • branches/eam_branches/ipp-20120627/psphot/src/psphotSourceSize.c

    r34179 r34247  
    1818    float sizeLimitCR;
    1919    float magLimitCR;
     20    int maxWindowCR;
    2021} psphotSourceSizeOptions;
    2122
    … …  
    2627bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
    2728bool psphotSourceSelectCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options);
    28 bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal);
     29bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal, int maxWindowCR);
    2930bool psphotMaskCosmicRayFootprintCheck (psArray *sources);
    3031int  psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh);
    … …  
    139140    if (!status) {
    140141        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CRMASK.APPLY is not defined.");
     142        return false;
     143    }
     144    options.maxWindowCR =  psMetadataLookupS32 (&status, recipe, "PSPHOT.CR.MAX.WINDOW");
     145    if (!status) {
     146        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CR.MAX.WINDOW is not defined.");
    141147        return false;
    142148    }
    … …  
    647653        psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf);
    648654        if (options->applyCRmask) {
    649             psphotMaskCosmicRay(readout, source, options->crMask);
     655            psphotMaskCosmicRay(readout, source, options->crMask, options->maxWindowCR);
    650656        } else {
    651657            source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    … …  
    687693// does no repair or recovery of the CR pixels, it only masks them out.  My test code can be
    688694// found at /data/ipp031.0/watersc1/psphot.20091209/algo_check.c
    689 bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal) {
     695bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal, int maxWindowCR) {
    690696
    691697    // Get the actual images and information about the peak.
    … …  
    699705    int ys = footprint->bbox.y0;
    700706    int ye = footprint->bbox.y1 + 1;
     707
     708    // We occasionally get very large footprints. When the footprint bounding box is large this function will
     709    // do an incredible amount of work for no benefit.
     710    // Limit the size of the area we examine.
     711    if (xe - xs > maxWindowCR) {
     712        xs = peak->x - maxWindowCR / 2;
     713        xe = xs + maxWindowCR;
     714    }
     715    if (ye - ys > maxWindowCR) {
     716        ys = peak->y - maxWindowCR / 2;
     717        ye = ys + maxWindowCR;
     718    }
    701719
    702720    LIMIT_XRANGE(xs, mask);
    … …  
    741759    for (int i = 0; i < footprint->spans->n; i++) {
    742760        pmSpan *sp = footprint->spans->data[i];
    743         for (int j = sp->x0; j <= sp->x1; j++) {
    744             int y = sp->y - ys;
    745             int x = j - xs;
     761        int y = sp->y - ys;
     762        if (y < 0 || y >= mymask->numRows) {
     763            // can we break if y >= numRows?
     764            continue;
     765        }
     766        for (int x = PS_MAX(sp->x0 - xs, 0); x <= PS_MIN(sp->x1 - xs, mymask->numCols-1); x++) {
    746767            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01;
    747768        }
    … …  
    883904
    884905bool psphotMaskCosmicRayFootprintCheck (psArray *sources) {
    885 
     906#ifdef CHECK_FOOTPRINTS
     907    // This gets really expensive for complex images
    886908    for (int i = 0; i < sources->n; i++) {
    887909        pmSource *source = sources->data[i];
    … …  
    894916        }
    895917    }
     918#endif
    896919    return true;
    897920}

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

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/eam_branches/ipp-20120627/psphot/src/psphotStackMatchPSFs.c

    r32868 r34247  
    8787    }
    8888
    89     // copy the header data from Src to Out
    90     // pmHDU *hduSrc = pmHDUFromReadout(readoutSrc);
    91     // psAssert (hduSrc, "input missing hdu?");
    92 
    93 
    9489    // set NAN pixels to 'SAT'
    9590    psImageMaskType maskSat = pmConfigMaskGet("SAT", config);
    … …  
    116111    // save the output fwhm values in the readout->analysis.  we may have / will have multiple output PSF sizes,
    117112    // so we save this in a vector.  if the vector is not yet defined, create it
     113    // Skip this if the readout deconvolution fraction was over the limit.
    118114    // NOTE: fwhmValues as defined here has 1 + nMatched PSF : 0 == unmatched
    119115    psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
    120116    psVectorAppend(fwhmValues, NAN); // XXX this corresponds to the unmatched image set
    121     for (int i = 0; i < options->targetSeeing->n; i++) {
    122         psVectorAppend(fwhmValues, options->targetSeeing->data.F32[i]);
     117
     118    bool overLimit = psMetadataLookupBool(NULL, readoutOut->analysis, "DECONV.OVERLIMIT");
     119    if (!overLimit) {
     120        for (int i = 0; i < options->targetSeeing->n; i++) {
     121            psVectorAppend(fwhmValues, options->targetSeeing->data.F32[i]);
     122        }
    123123    }
    124124    psMetadataAddVector(readoutSrc->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);

  • branches/eam_branches/ipp-20120627/psphot/src/psphotStackMatchPSFsNext.c

    r32996 r34247  
    77    int nRadialEntries = 0;
    88
    9     // find the currently selected readout
    10     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest
    11     psAssert (file, "missing file?");
     9    // find the numer of fwhmValues in the first non-skipped input
     10    int num = psphotFileruleCount(config, filerule);
     11    for (int i=0; i<num; i++) {
     12        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
     13        psAssert (file, "missing file?");
    1214   
    13     pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    14     psAssert (readout, "missing readout?");
    15    
    16     bool status = false;
    17     psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
    18     if (!fwhmValues) {
    19         return 1;
     15        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     16        psAssert (readout, "missing readout?");
     17        bool status = false;
     18        if (!psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) {
     19            psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
     20            if (fwhmValues) {
     21                nRadialEntries = fwhmValues->n;
     22            } else {
     23                nRadialEntries = 1;
     24            }
     25            break;
     26        }
    2027    }
    21    
    22     nRadialEntries = fwhmValues->n;
    2328    return nRadialEntries;
    2429}
    … …  
    6065    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    6166    psAssert (readout, "missing readout?");
     67
     68    if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) {
     69        psLogMsg("psphot", PS_LOG_INFO, "skipping smooth %d to next psf", index);
     70        return true;
     71    }
    6272
    6373    psLogMsg("psphot", PS_LOG_INFO, "smooth %d to next psf", index);

  • branches/eam_branches/ipp-20120627/psphot/src/psphotStackMatchPSFsUtils.c

    r33841 r34247  
    310310    float deconv = psMetadataLookupF32(NULL, readoutOut->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
    311311    if (deconv > deconvLimit) {
     312#if (1)
     313        // XXX: don't reject the image set
     314        psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image for matched psf analysis%d\n", deconv, deconvLimit, index);
     315        psMetadataAddBool(readoutOut->analysis, PS_LIST_TAIL, "DECONV.OVERLIMIT", PS_META_REPLACE, "", true);
     316#else
    312317        psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index);
    313318        goto escape;
     319#endif
     320    } else {
     321        psMetadataAddBool(readoutOut->analysis, PS_LIST_TAIL, "DECONV.OVERLIMIT", PS_META_REPLACE, "", false);
    314322    }
    315323
    … …  
    461469    return optWidths;
    462470}
     471
     472// Set input to be skipped if the decovolution fraction was overlimit. Use for radial apertures
     473// This interface can be potentiall be extended for other uses
     474bool psphotStackSetInputsToSkip(pmConfig *config, const pmFPAview *view, const char *filerule, bool set) {
     475    int num = psphotFileruleCount(config, filerule);
     476    bool status;
     477    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
     478    if (!status) chisqNum = -1;
     479    for (int i = 0; i < num; i++) {
     480        if (i == chisqNum) {
     481            continue;
     482        }
     483        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
     484        psAssert (file, "missing file?");
     485
     486        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     487        psAssert (readout, "missing readout?");
     488        if (set) {
     489            bool overLimit = psMetadataLookupBool(&status, readout->analysis, "DECONV.OVERLIMIT");
     490            psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSPHOT.SKIP.INPUT", PS_META_REPLACE, "Skip analysis", overLimit);
     491        } else {
     492            psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSPHOT.SKIP.INPUT", PS_META_REPLACE, "Skip analysis", false);
     493        }
     494    }
     495
     496    return true;
     497}

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

    r33958 r34247  
    11# include "psphotInternal.h"
     2
     3static bool psphotStackLoadWCS(pmConfig *config, const pmFPAview *view, const char *filerule);
    24
    35// we have 3 possible real filesets:
    … …  
    6769    char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
    6870    char *STACK_DET = STACK_RAW;
     71
     72    // load WCS
     73    if (!psphotStackLoadWCS(config, view, STACK_SRC)) {
     74        psError (PSPHOT_ERR_CONFIG, false, "trouble loading WCS for %s", STACK_SRC);
     75        return false;
     76    }
    6977
    7078    // set the photcode for each image
    … …  
    179187
    180188    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    181     psphotFitSourcesLinear (config, view, STACK_SRC, false);
     189    psphotFitSourcesLinear (config, view, STACK_SRC, false, false);
    182190    psphotStackVisualFilerule(config, view, STACK_SRC);
    183191
    … …  
    207215    // NOTE : apply to ALL sources (extended + psf)
    208216    // NOTE 2 : this function subtracts the models from the given filerule (SRC), not DET
    209     psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources)
     217    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 2 (detections->allSources)
    210218
    211219    // NOTE: possibly re-measure background model here with objects subtracted / or masked
    … …  
    278286    }
    279287
     288    // gcc doesn't like the label to refer to a declaration so we declare this here
     289    bool splitLinearFit;
     290
    280291pass1finish:
     292
     293    splitLinearFit = psMetadataLookupBool(NULL, recipe, "PSPHOT.STACK.SPLIT.LINEAR.FIT");
     294    if (splitLinearFit) {
     295        psLogMsg ("psphot", 3, "splitting fit of detected and matched soures\n");
     296        // Fit the detected sources separately from matched that wea are about to create.
     297        // NOTE: apply to ALL sources but only include sources with postitive flux in the fit
     298        psphotFitSourcesLinear (config, view, STACK_SRC, true, true); // pass 3 (detections->allSources)
     299    }
    281300
    282301    // generate the objects (objects unify the sources from the different images) NOTE: could
    … …  
    286305    psMemDump("matchsources");
    287306
     307    // check the source density. If it too high change the number of radial bins
     308    // in the recipe.
     309    psphotLimitRadialApertures(recipe, objects->n);
     310
    288311    // Construct an initial model for each object, set the radius to fitRadius, set circular
    289312    // fit mask.  NOTE: only applied to sources without guess models
    … …  
    294317    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
    295318
    296     // NOTE: apply to ALL sources
    297     psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)
     319    if (splitLinearFit) {
     320        // NOTE: apply to Matched sources. Since the sources that we fit above are subtracted, they will
     321        // not be included in this fit.
     322        psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
     323    } else {
     324        // Fit all sources together
     325        psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 3 (detections->allSources)
     326    }
    298327
    299328    // measure the radial profiles to the sky (only measures new objects)
    … …  
    344373            }
    345374        }
    346         psphotRadialApertures (config, view, STACK_CNV, 0); // XXX entry 0 == unmatched?
     375        psphotRadialApertures (config, view, STACK_CNV, 0); // entry 0 == unmatched
    347376        psMemDump("extmeas");
    348377
     378        // mark any inputs that we want to skip the matched apertures for
     379        psphotStackSetInputsToSkip(config, view, STACK_OUT, true);
     380
    349381        int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
    350 
    351382        for (int entry = 1; entry < nRadialEntries; entry++) {
    352383            // NOTE: entry 0 is the unmatched image set
    … …  
    359390
    360391            // this is necessary to get the right normalization for the new models
    361             psphotFitSourcesLinear (config, view, STACK_OUT, false);
     392            psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
    362393
    363394            // measure circular, radial apertures (objects sorted by S/N)
    364             // entry 0 == unmatched? pass entry + 1?
    365395            psphotRadialApertures (config, view, STACK_OUT, entry);
    366396
    … …  
    373403        }
    374404    }
     405    psphotStackSetInputsToSkip(config, view, STACK_OUT, false);
    375406
    376407    // measure aperture photometry corrections
    … …  
    413444    // create the exported-metadata and free local data
    414445    return psphotReadoutCleanup (config, view, STACK_SRC);
     446}
     447
     448static bool psphotStackLoadWCS(pmConfig *config, const pmFPAview *view, const char *filerule) {
     449
     450    int num = psphotFileruleCount(config, filerule);
     451
     452    for (int i=0; i<num; i++) {
     453        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
     454        psAssert (file, "missing file?");
     455   
     456        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     457        psAssert (readout, "missing readout?");
     458
     459        pmHDU *hdu = pmHDUFromReadout(readout);
     460        psAssert (hdu, "input missing hdu?");
     461
     462        if (!pmAstromReadWCS(file->fpa, readout->parent->parent, hdu->header, 1.0)) {
     463            psError (PSPHOT_ERR_UNKNOWN, false, "failed to read WCS from header for input %d", i);
     464            return false;
     465        }
     466    }
     467    return true;
    415468}
    416469
Note: See TracChangeset for help on using the changeset viewer.