IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33953


Ignore:
Timestamp:
May 29, 2012, 5:16:09 PM (14 years ago)
Author:
eugene
Message:

changes from trunk

Location:
branches/eam_branches/ipp-20120405/psphot
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20120405/psphot

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

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

    r33946 r33953  
    189189// used by psphotFindDetections
    190190pmReadout      *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal);
    191 psArray        *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax);
     191psArray        *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks);
    192192bool            psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal);
    193193psErrorCode     psphotCullPeaks(const pmReadout *readout, const pmReadout *signifRO, const psMetadata *recipe, psArray *footprints);
    … …  
    472472
    473473bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule);
    474 bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf);
     474bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index);
    475475bool psphotKronIterate_Threaded (psThreadJob *job);
    476476

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

    r32348 r33953  
    4848    // Use the new pmFootprints approach?
    4949    const bool useFootprints = psMetadataLookupBool(NULL, recipe, "USE_FOOTPRINTS");
     50    const bool footprintUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED");
    5051
    5152    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    … …  
    6465    }
    6566
     67    bool replaceSourcesForFootprints = false;
    6668    if (firstPass) {
    6769        pass = 1;
    … …  
    7072    } else {
    7173        pass = 2;
     74        replaceSourcesForFootprints = footprintUseUnsubtracted;   
    7275        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); PS_ASSERT (status, NULL);
    7376        NMAX = 0; // unlimited number of peaks in final pass: allow a limit (PEAKS_NMAX_2) ?
    … …  
    9699
    97100    // detect the peaks in the significance image
    98     detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX);
     101    int totalPeaks = 0;
     102    detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX, &totalPeaks);
    99103    psMetadataAddF32  (readout->analysis, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold);
    100104    if (!detections->peaks) {
    … …  
    105109        return false;
    106110    }
     111    // hard limit on number of peaks we will accept. (To avoid memory overload in psphotStack)
     112    int maxPeaks = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX_TOTAL"); PS_ASSERT (status, NULL);
     113    if (maxPeaks && (totalPeaks > maxPeaks)) {
     114        psFree (detections);
     115        psError (PSPHOT_ERR_DATA, true, "Too many peaks %d found PEAKS_NMAX_TOTAL: %d", totalPeaks, maxPeaks);
     116        return false;
     117    }
    107118
    108119    // optionally merge peaks into footprints
    109120    if (useFootprints) {
     121        if (replaceSourcesForFootprints) {
     122            psphotAddOrSubNoiseReadout(config, view, filerule, index, recipe, false);
     123            psphotReplaceAllSourcesReadout (config, view, filerule, index, recipe, false);
     124            psFree (significance);
     125            significance = psphotSignificanceImage (readout, recipe, maskVal);
     126        }
     127
    110128        psphotFindFootprints (detections, significance, readout, recipe, threshold, pass, maskVal);
     129
     130        if (replaceSourcesForFootprints) {
     131            psphotRemoveAllSourcesReadout (config, view, filerule, index, recipe, false);
     132        }
    111133    }
    112134

  • branches/eam_branches/ipp-20120405/psphot/src/psphotFindPeaks.c

    r32348 r33953  
    44// image must be constructed to represent (S/N)^2.  If nMax is non-zero, only return a maximum
    55// of nMax peaks
    6 psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax) {
     6psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks) {
    77
    88    bool status = false;
    … …  
    1818        psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
    1919        return NULL;
     20    }
     21    // return the total number of peaks found before the nMax limit is applied
     22    if (totalPeaks) {
     23        *totalPeaks = peaks->n;
    2024    }
    2125

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

    r33946 r33953  
    11# include "psphotInternal.h"
    2 # ifndef ROUND
    3 # define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
    4 # endif
    5 
    6 bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert);
    7 bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal);
     2
     3bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool oldWindow);
     4bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, psImage *smoothedPixels);
     5
    86
    97bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule)
    … …  
    119    bool status = true;
    1210
    13     // return true;
    14 
    1511    fprintf (stdout, "\n");
    1612    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Kron Iterate ---");
    … …  
    4642        // psAssert (psf, "missing psf?");
    4743
    48         if (!psphotKronIterateReadout (config, recipe, view, readout, sources, psf)) {
     44        if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i)) {
    4945            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
    5046            return false;
    … …  
    5854bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
    5955
    60 bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
     56bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index) {
    6157
    6258    bool status = false;
    … …  
    6965    psTimerStart ("psphot.kron");
    7066
     67
    7168    // determine the number of allowed threads
    7269    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    … …  
    8784    int KRON_ITERATIONS = psMetadataLookupS32 (&status, recipe, "KRON_ITERATIONS");
    8885    if (!status) {
    89         KRON_ITERATIONS = 2;
     86        KRON_ITERATIONS = 1;
     87    }
     88
     89    bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
     90    if (!status) {
     91        KRON_APPLY_WEIGHT = true;
     92    }
     93
     94    bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
     95    if (!status) {
     96        KRON_APPLY_WINDOW = false;
     97    }
     98    bool KRON_SMOOTH = psMetadataLookupBool (&status, recipe, "KRON_SMOOTH");
     99    if (!status) {
     100        KRON_SMOOTH = false;
     101    }
     102    float KRON_SMOOTH_SIGMA = psMetadataLookupF32 (&status, recipe, "KRON_SMOOTH_SIGMA");
     103    if (!status) {
     104        KRON_SMOOTH_SIGMA = 1.7;
     105    }
     106    float KRON_SMOOTH_NSIGMA = psMetadataLookupS32 (&status, recipe, "KRON_SMOOTH_NSIGMA");
     107    if (!status) {
     108        KRON_SMOOTH_NSIGMA = 2;
    90109    }
    91110
    … …  
    119138        // set a window function for each source based on the moments
    120139        // (this skips really bad sources (no peak, no moments, DEFECT)
    121         psphotKronWindowSetSource (source, kronWindow, false, true);
     140        psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW);
     141    }
     142
     143    // We measure the Kron Radius on a smoothed copy of the readout image
     144    psImage *smoothedImage = NULL;
     145    if (KRON_SMOOTH) {
     146        // Build the smoothed source image
     147        // Replace the subtracted sources
     148        psphotReplaceAllSourcesReadout(config, view, filerule, index, recipe, false);
     149        // Copy the image and smooth
     150        psTimerStart ("psphot.kron.smooth");
     151        smoothedImage = psImageCopy(NULL, readout->image, PS_TYPE_F32);
     152        psImageSmooth(smoothedImage, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA);
     153        psLogMsg ("psphot.kron", PS_LOG_INFO, "smoothed image %f sec\n", psTimerMark ("psphot.kron.smooth"));
     154
     155        // remove the sources
     156        psphotRemoveAllSourcesReadout( config, view, filerule, index, recipe, false );
     157        // Now subtract smooth versions of the sources from the smoothed image
     158        psTimerStart ("psphot.kron.smooth.sources");
     159        for (int i=0; i< sources->n; i++) {
     160            pmSource *source = sources->data[i];
     161            // If source has been subtracted from the readout image subtract a "smoothed" version from the smoothedImage
     162            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     163                // cache copy of smoothedPixels in the source
     164                // tmpPtr is for use by a single "module" and must be null otherwise
     165                psAssert(source->tmpPtr == NULL, "source->tmpPtr is not null!");
     166
     167                psImage *smoothedPixels = psImageSubset(smoothedImage, source->region);
     168                source->tmpPtr = (psPtr) smoothedPixels;
     169                pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0);
     170            }
     171        }
     172        psLogMsg ("psphot.kron", PS_LOG_INFO, "removed %ld smoothed sources %f sec\n", sources->n, psTimerMark ("psphot.kron.smooth.sources"));
     173
    122174    }
    123175
    … …  
    146198            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS,    PS_TYPE_F32);
    147199            PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS,    PS_TYPE_S32);
     200            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32);
     201            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32);
     202            PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA,  PS_TYPE_F32);
     203            psArrayAdd(job->args, 1, smoothedImage);
    148204
    149205// set this to 0 to run without threading
    150 # if (0)
     206# if (1)
    151207            if (!psThreadJobAddPending(job)) {
    152208                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    … …  
    156212            if (!psphotKronIterate_Threaded(job)) {
    157213                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     214                // psFree(AnalysisRegion);
    158215                return false;
    159216            }
    … …  
    179236    psFree (cellGroups);
    180237    psFree (kronWindow);
     238    psFree (smoothedImage);
    181239
    182240    psLogMsg ("psphot.kron", PS_LOG_WARN, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
    … …  
    194252    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
    195253    int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[7],S32);
    196 
    197     // XXX TEST : set iteration to 1
     254    bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[8],S32);
     255    bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[9],S32);
     256    float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[10],F32);
     257    psImage *smoothedImage          = job->args->data[11];
     258
     259    // psImage *smoothedPixels = NULL;
    198260    for (int j = 0; j < KRON_ITERATIONS; j++) {
    199261        for (int i = 0; i < sources->n; i++) {
    … …  
    202264            if (!source->peak) continue; // XXX how can we have a peak-less source?
    203265
    204             // allocate space for moments
     266            // check status of this source's moments
    205267            if (!source->moments) continue;
     268            if (!source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue;
     269            if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
    206270
    207271            // replace object in image
    208272            bool reSubtract = false;
     273            psImage *smoothedPixels = NULL;
    209274            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    210                 pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     275                pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     276                smoothedPixels = (psImage *) source->tmpPtr;
     277                if (smoothedPixels) {
     278                    // psFree(source->tmpPtr);
     279                    // smoothedPixels = psImageSubset(smoothedImage, source->region);
     280                    pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, true, maskVal, 0, 0);
     281                    // source->tmpPtr = smoothedPixels;
     282                }
    211283                reSubtract = true;
    212284            }
    … …  
    223295            // XXX float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow);
    224296
    225             // Sextractor apparently takes something like the skyRadius window, measures the moments inside that window,
    226             // then uses an elliptical outline based on the 2nd moments, but 6x the size
    227 
    228297            // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS
    229298            // if we lack the skyRadius (eg MATCHED sources), go to the default value
    230             float maxWindow, windowRadius;
    231            
    232             if (j == 0) {
    233                 maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
    234                 windowRadius = PS_MAX(RADIUS, maxWindow);
    235             } else {
    236                 maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS;
    237                 windowRadius = PS_MAX(RADIUS, maxWindow);
    238             }
    239 
    240 #ifdef notdef
    241             fprintf(stderr, "Redefining pixels for source: %d %4d %4d new radius: %f\n",
    242                                         i, source->peak->x, source->peak->y, windowRadius+2);
    243 #endif
     299            float maxWindow;
     300            if (j == 0) {
     301                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
     302            } else {
     303                maxWindow = isfinite(source->moments->Mrf) ?  6.0*source->moments->Mrf : RADIUS;
     304            }
     305            float windowRadius = PS_MAX(RADIUS, maxWindow);
     306
    244307            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    245             pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
     308            bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    246309            psAssert (source->pixels, "WTF?");
     310            if (extend && smoothedPixels) {
     311                psFree(source->tmpPtr);
     312                smoothedPixels = psImageSubset(smoothedImage, source->region);
     313                psAssert (smoothedPixels, "WTF?");
     314                source->tmpPtr = (psPtr) smoothedPixels ;
     315            }
     316
    247317
    248318            // clear the window function for this source based on the moments
    249             psphotKronWindowSetSource (source, kronWindow, (j > 0), false);
    250 
    251             // this function populates moments->Mrf,KronFlux,KronFluxErr
    252             psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
    253             psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    254 
    255             // set a window function for each source based on the moments
    256             psphotKronWindowSetSource (source, kronWindow, true, true);
     319            // Note this function also applies cuts on the source and returns false if it
     320            // does not meet the requirements for measuring the Kron Radius or Magnitude.
     321            // Not that it performs that function even if KRON_APPLY_WINDOW is false
     322            if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) {
     323
     324                // this function populates moments->Mrf,KronFlux,KronFluxErr
     325                psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, smoothedPixels);
     326                psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     327
     328                // set a window function for each source based on the moments
     329                psphotKronWindowSetSource (source, kronWindow, true, true, KRON_APPLY_WINDOW);
     330            }
    257331
    258332            // if we subtracted it above, re-subtract the object, leave local sky
    259333            if (reSubtract) {
    260                 pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     334                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     335                if (smoothedPixels) {
     336                    pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0);
     337                    if (j + 1 == KRON_ITERATIONS) {
     338                        // We're done with the smoothedPixels
     339                        psFree(source->tmpPtr);
     340                        source->tmpPtr = NULL;
     341                        smoothedPixels = NULL;
     342                    }
     343                }
    261344            }
    262345        }
    263346    }
     347    // psFree(smoothedPixels);
    264348    return true;
    265349}
    266350
    267 bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal) {
     351bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal,
     352    bool applyWeight, psImage *smoothedPixels) {
    268353
    269354    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    273358
    274359    psF32 R2 = PS_SQR(radius);
    275     float rsigma2 = 0.5 / PS_SQR(radius/2.0);
     360    float rsigma2 =  applyWeight ? 0.5 / R2 : 0;
    276361
    277362    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    286371    // Xn  = SUM (x - xc)^n * (z - sky)
    287372
     373
    288374    psF32 RF = 0.0;
    289375    psF32 RS = 0.0;
    … …  
    304390    int Ywo = source->pixels->row0;
    305391
    306     psF32 **vPix = source->pixels->data.F32;
     392    psF32 **vPix;
     393    if (smoothedPixels) {
     394        vPix = smoothedPixels->data.F32;
     395    } else {
     396        vPix = source->pixels->data.F32;
     397    }
    307398    psF32 **vWin = kronWindow->data.F32;
    308399    psF32 **vWgt = source->variance->data.F32;
    … …  
    377468    float Var = 0.0;
    378469    float Win = 0.0;
     470
     471    // set vPix to the source pixels (it may have been set to the
     472    // smoothed image above)
     473    vPix = source->pixels->data.F32;
     474
    379475
    380476    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    415511}
    416512
    417 bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert) {
     513bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) {
    418514
    419515    if (!source) return false;
    … …  
    422518    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
    423519    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
     520    if (!source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) return false;
     521    if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false;
    424522    psAssert(kronWindow, "need a window");
     523
     524    if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
     525
     526    if (!applyWindow) {
     527        return true;
     528    }
    425529
    426530    // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has
    … …  
    433537    float Yo = source->moments->My;
    434538
    435     float Mxx = source->moments->Mxx;
    436     float Mxy = source->moments->Mxy;
    437     float Myy = source->moments->Myy;
    438 
    439     float Mmajor = 0.5*(Mxx + Myy) + 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
    440     float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
    441 
    442     // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width
    443     float scale = PS_SQR(0.5 * source->moments->Mrf) / Mmajor;
    444 
    445     float Sxx = scale * Mmajor * Mminor / Myy; // sigma_x^2
    446     float Sxy = Mxy / (scale * Mmajor * Mminor);
    447     float Syy = scale * Mmajor * Mminor / Mxx; // sigma_y^2
    448 
    449     float Smajor = sqrt(Mmajor);
     539    psEllipseMoments moments;
     540    moments.x2 = source->moments->Mxx;
     541    moments.y2 = source->moments->Myy;
     542    moments.xy = source->moments->Mxy;
     543
     544    psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
     545    if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) {
     546        // Shall we log a proper warning? This happens often with matched sources (forced photometry)
     547        // fprintf(dump, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta);
     548        return false;
     549    }
     550
     551    // Why this factor of 0.5 ?
     552    float scale = 0.5 * source->moments->Mrf / axes.major;
     553    axes.major *= scale;
     554    axes.minor *= scale;
     555
     556    psEllipseShape shape = psEllipseAxesToShape(axes);
     557    if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) {
     558        // Shall we log a proper warning? This happens often with matched sources (forced photometry)
     559        // fprintf(dump, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy);
     560        return false;
     561    }
     562
     563    float Smajor = axes.major;
    450564
    451565    int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1);
    … …  
    454568    int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1);
    455569
    456     float rMxx = 0.5 / Sxx;
    457     float rMyy = 0.5 / Syy;
    458    
     570    float rMxx = 0.5 / PS_SQR(shape.sx);
     571    float rMyy = 0.5 / PS_SQR(shape.sy);
     572    float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
     573                                    // implementation. XXX: Is this correct?
     574
    459575    for (int iy = minY; iy < maxY; iy++) {
    460576        for (int ix = minX; ix < maxX; ix++) {
    … …  
    467583            float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z));
    468584
    469             kronWindow->data.F32[iy][ix] *= f;
     585            kronWindow->data.F32[iy][ix] *= f;
    470586        }
    471587    }
     588
    472589    return true;
    473590}

  • branches/eam_branches/ipp-20120405/psphot/src/psphotMergeSources.c

    r33444 r33953  
    744744    psAssert (readoutOut, "missing readout?");
    745745
     746    pmDetections *detectionsOutOld = psMetadataLookupPtr (&status, readoutOut->analysis, "PSPHOT.DETECTIONS");
     747    psArray *oldFootprints = detectionsOutOld ? detectionsOutOld->footprints : NULL;
     748
    746749    // replace any existing DETECTION container on readoutOut->analysis with the new one
    747750    pmDetections *detectionsOut = pmDetectionsAlloc();
     751    if (oldFootprints) {
     752        // ... but hang on to any existing footprints so that they can be merged with new footprints in pass 2
     753        detectionsOut->footprints = psMemIncrRefCounter(oldFootprints);
     754    }
    748755    detectionsOut->allSources = psArrayAllocEmpty (100);
    749756    if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detectionsOut)) {

  • branches/eam_branches/ipp-20120405/psphot/src/psphotRadialProfileWings.c

    r33411 r33953  
    5555static float MIN_RADIUS = NAN;
    5656static float SKY_STDEV  = NAN;
     57static float SKY_SLOPE_MIN = NAN;
    5758// static FILE *file = NULL;
    5859
    … …  
    8182
    8283    MIN_RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS");
    83     if (!status) {
     84    if (!status || (MIN_RADIUS > MAX_RADIUS)) {
    8485        MIN_RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    8586    }
    … …  
    8990    if (!status) {
    9091        SKY_STDEV = 1.0; // a crude default value (why would this not exist?)
     92    }
     93
     94    // SKY_SLOPE_MIN is the sigma of the sky model (ie, smoothed on large scales)
     95    SKY_SLOPE_MIN = psMetadataLookupF32 (&status, recipe, "SKY_SLOPE_MIN");
     96    if (!status) {
     97        SKY_SLOPE_MIN = 3.0;
    9198    }
    9299
    … …  
    302309            limitSlope = slope;
    303310        }
    304         if (!limit && isfinite(slope) && (fabs(slope) < 3.0)) {
     311        if (!limit && isfinite(slope) && (fabs(slope) < SKY_SLOPE_MIN)) {
    305312            // SB no longer changing.       
    306313            limit = true;
    307314            // linearly interpolate to the radius at which we hit the sky, using the last flux and the limiting slope
    308315            if (isfinite(lastFlux)) {
    309                 float interpolatedRadius = lastRadius + lastFlux / 3.0;
     316                float interpolatedRadius = lastRadius + lastFlux / SKY_SLOPE_MIN;
    310317                if (interpolatedRadius < MAX_RADIUS) {
    311318                    limitRadius = interpolatedRadius;

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

    r33089 r33953  
    185185        // remove noise for subtracted objects (ie, return to normal noise level)
    186186        // NOTE: this needs to operate only on the OLD sources
    187         psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources)
     187        bool footprintUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED");
     188        // Note: if footprintUseUnsubtracted is true the noise was already subtracted in psphotFindDetections()
     189        if (!footprintUseUnsubtracted) {
     190            psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources)
     191        }
    188192
    189193        // define new sources based on only the new peaks & measure moments

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

    r32996 r33953  
    6464    // Construct an initial model for each object, set the radius to fitRadius, set circular
    6565    // fit mask.  NOTE: only applied to sources without guess models
    66     psphotGuessModels (config, view, filerule);
     66    if (!psphotGuessModels (config, view, filerule)) {
     67        psLogMsg ("psphot", 3, "failure to Guess Model");
     68        return psphotReadoutCleanup (config, view, filerule);
     69    }
    6770
    6871    // linear PSF fit to source peaks

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

    r33946 r33953  
    3030    psFree(task);
    3131
    32     task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 8);
     32    task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 12);
    3333    task->function = &psphotKronIterate_Threaded;
    3434    psThreadTaskAdd(task);

  • branches/eam_branches/ipp-20120405/psphot/src/psphotSourceStats.c

    r33089 r33953  
    605605    }
    606606    psAssert (isfinite(Sigma), "did we miss a case?");
     607    if (!isfinite(minKronRadius)) {
     608        minKronRadius = Sigma;
     609    }
    607610
    608611    // choose a grid scale that is a fixed fraction of the psf sigma^2

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

    r33946 r33953  
    4242    // XXX for now, just load the full set of images up front except for EXPNUM which we defer
    4343    pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.RAW");
     44    pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.CNV");
    4445    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
    4546

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

    r32348 r33953  
    105105    psAssert(psphotRecipe, "Need PSPHOT recipe");
    106106
     107#if (0)
     108    // DON'T OVERWRITE PSPHOT recipe value
    107109    psString maskBadStr = psMetadataLookupStr(NULL, ppStackRecipe, "MASK.BAD");// Name of bits to mask for bad
    108110    psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    109111
    110112    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    111     psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad);
     113    psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad);
     114#endif
    112115
    113116    psImage *binned = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model

  • branches/eam_branches/ipp-20120405/psphot/src/psphotStackObjects.c

    r32633 r33953  
    119119            }
    120120            if (!skipSourcePetro) keepObjectPetro = true;
    121 
    122             keepObjectPetro = true;
    123121        }
    124122

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

    r33690 r33953  
    129129        }
    130130
     131        psS64 stack_id = psMetadataLookupS64(&status, input, "STACK_ID");
     132        if (!status) {
     133            psError(PS_ERR_UNKNOWN, false, "Unable to find STACK_ID from sources %d", i);
     134            return false;
     135        }
    131136        // generate an pmFPAimage for the output convolved image
    132137        // XXX output of these files should be optional
    … …  
    139144            }
    140145            outputImage->save = true;
    141             outputImage->fileID = i;            // this is used to generate output names
     146            outputImage->fileID = stack_id;             // this is used to generate output names
    142147
    143148            pmFPAfile *outputMask = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.MASK");
    … …  
    151156            }
    152157            outputMask->save = true;
    153             outputMask->fileID = i;             // this is used to generate output names
     158            outputMask->fileID = stack_id;              // this is used to generate output names
    154159
    155160            pmFPAfile *outputVariance = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.VARIANCE");
    … …  
    163168            }
    164169            outputVariance->save = true;
    165             outputVariance->fileID = i;         // this is used to generate output names
     170            outputVariance->fileID = stack_id;          // this is used to generate output names
    166171
    167172            // the output sources are carried on the outputImage->fpa structures
    … …  
    172177            }
    173178            outsources->save = true;
    174             outsources->fileID = i;             // this is used to generate output names
     179            outsources->fileID = stack_id;              // this is used to generate output names
    175180        }
    176181    }

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

    r32996 r33953  
    8989        return psphotReadoutCleanup (config, view, STACK_SRC);
    9090    }
    91     if (!psphotSubtractBackground (config, view, STACK_SRC)) {
    92         return psphotReadoutCleanup (config, view, STACK_SRC);
     91    if (strcmp(STACK_SRC, STACK_DET)) {
     92#define MODEL_BACKGROUND_SRC 1
     93#ifdef MODEL_BACKGROUND_SRC
     94        // work around the fact that the background levels on the convolved
     95        // and unconvolved stacks can be different
     96        if (!psphotModelBackground (config, view, STACK_SRC)) {
     97            return psphotReadoutCleanup (config, view, STACK_SRC);
     98        }
     99#endif
     100        if (!psphotSubtractBackground (config, view, STACK_SRC)) {
     101            return psphotReadoutCleanup (config, view, STACK_SRC);
     102        }
    93103    }
    94104    if (!strcasecmp (breakPt, "BACKMDL")) {
    … …  
    190200    psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
    191201
     202
    192203    // if we only do one pass, skip to extended source analysis
    193204    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
    … …  
    211222            //  subtract all sources from DET (this will subtract using the psf model for SRC, which
    212223            //  will somewhat oversubtract the sources -- this is OK
    213             psphotRemoveAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
     224            psphotRemoveAllSources (config, view, STACK_DET, false); // do not ignore subtraction state for sources
    214225        }
    215226
    … …  
    223234        // remove noise for subtracted objects (ie, return to normal noise level)
    224235        // NOTE: this needs to operate only on the OLD sources
    225         psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
     236        // NOTE: if fooprintsUseUnsubtracted, the noise has already been removed by psphotFindDetections
     237        bool footprintsUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED");
     238        if (!footprintsUseUnsubtracted) {
     239            psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
     240        }
    226241
    227242        // if DET and SRC are different images, copy the detections from DET to SRC
    … …  
    303318
    304319    // create source children for the OUT filerule (for radial aperture photometry)
    305     psArray *objectsRadial = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
    306     if (!objectsRadial) {
     320    // (These are not just for radial aperture photomoetry. If they aren't defined we
     321    // get no sources output
     322    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
     323    if (!objectsOut) {
    307324        psFree(objects);
    308325        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    … …  
    310327    }
    311328
    312     // measure circular, radial apertures (objects sorted by S/N)
    313     // this forces photometry on the undetected sources from other images
    314     psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched?
    315     psMemDump("extmeas");
    316 
    317     int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
    318 
    319     for (int entry = 1; entry < nRadialEntries; entry++) {
    320         // NOTE: entry 0 is the unmatched image set
    321 
    322         // re-measure the PSF for the smoothed image (using entries in 'allSources')
    323         psphotChoosePSF (config, view, STACK_OUT, false);
    324 
    325         // this is necessary to update the models based on the new PSF
    326         psphotResetModels (config, view, STACK_OUT);
    327 
    328         // this is necessary to get the right normalization for the new models
    329         psphotFitSourcesLinear (config, view, STACK_OUT, false);
    330 
    331         // measure circular, radial apertures (objects sorted by S/N)
    332         // entry 0 == unmatched? pass entry + 1?
    333         psphotRadialApertures (config, view, STACK_OUT, entry);
    334 
    335         // replace the flux in the image so it is returned to its original state
    336         psphotReplaceAllSources (config, view, STACK_OUT, false);
    337 
    338         // smooth to the next FWHM, or set 'smoothAgain' to false if no more
    339         psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
    340         psMemDump("matched");
     329
     330    bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES");
     331    if (radial_apertures) {
     332        // measure circular, radial apertures (objects sorted by S/N)
     333        // this forces photometry on the undetected sources from other images
     334        psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched?
     335        psMemDump("extmeas");
     336
     337        int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
     338
     339        for (int entry = 1; entry < nRadialEntries; entry++) {
     340            // NOTE: entry 0 is the unmatched image set
     341
     342            // re-measure the PSF for the smoothed image (using entries in 'allSources')
     343            psphotChoosePSF (config, view, STACK_OUT, false);
     344
     345            // this is necessary to update the models based on the new PSF
     346            psphotResetModels (config, view, STACK_OUT);
     347
     348            // this is necessary to get the right normalization for the new models
     349            psphotFitSourcesLinear (config, view, STACK_OUT, false);
     350
     351            // measure circular, radial apertures (objects sorted by S/N)
     352            // entry 0 == unmatched? pass entry + 1?
     353            psphotRadialApertures (config, view, STACK_OUT, entry);
     354
     355            // replace the flux in the image so it is returned to its original state
     356            psphotReplaceAllSources (config, view, STACK_OUT, false);
     357
     358            // smooth to the next FWHM, or set 'smoothAgain' to false if no more
     359            psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
     360            psMemDump("matched");
     361        }
    341362    }
    342363
    … …  
    344365    if (!psphotApResid (config, view, STACK_SRC)) {
    345366        psFree (objects);
     367        psFree (objectsOut);
    346368        psLogMsg ("psphot", 3, "failed on psphotApResid");
    347369        return psphotReadoutCleanup (config, view, STACK_SRC);
    … …  
    375397
    376398    psFree (objects);
    377     psFree (objectsRadial);
     399    psFree (objectsOut);
    378400
    379401    // create the exported-metadata and free local data
Note: See TracChangeset for help on using the changeset viewer.