IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27532


Ignore:
Timestamp:
Mar 30, 2010, 1:33:47 PM (16 years ago)
Author:
eugene
Message:

some API re-work to support psphotForced and detection efficiency measurements; better visualization; CR masking (optional); better source size analysis; some motion towards stack photometry

Location:
trunk/psphot
Files:
11 edited
2 copied

Legend:

Unmodified
Added
Removed

  • trunk/psphot/doc/notes.20100131.txt

    r26894 r27532  
     1
     22010.03.24
     3
     4 I have finished another iteration on source size analysis.  I feel
     5 fairly confident of the extended source and CR identifications.  They
     6 are not perfect, but they are in pretty good shape.  The CR IDs use
     7 the moments to find things smaller in one dimension than a PSF, while
     8 the EXT is set for things larger than a PSF in one dimention.
     9
     10 There still seems to be some trouble with SATSTARS : these are being
     11 set for sources which are fainter than the SAT limit -- I'm not sure
     12 if this is due to the quality of the PSF fit or if they really are
     13 SAT, but have integrated fluxes lower than expected (underfitted in
     14 PSF).
     15
     16 Additional Parameters needed from psphot for distinguishing dipoles:
     17
     18 -- note regarding CfA parameters --
     19
     20  Here are our settings.
     21 
     22  Note that Fpos is the flux which is 3-sigma above background, and
     23  Fneg is the flux which is 3-sigma below background.
     24 
     25  Npos is the # of pixels with flux which is 3-sigma above background.
     26 
     27  etc.
     28 
     29  # for positive objects:
     30  # Ngood = Npos, Nbad = Nneg, Fgood = Fpos, Fbad = Fneg
     31
     32  # for negative objects:
     33  # Ngood = Nneg, Nbad = Npos, Fgood = Fneg, Fbad = Fpos
     34  # FRATIO=Fgood/(Fgood+Fbad)
     35  # NRATIO_BAD=Ngood/(Ngood+Nbad)
     36  # NRATIO_MASK=Ngood/(Ngood+Nmask)
     37  # NRATIO_ALL=Ngood/(Ngood+Nbad+Nmask)
     38  DC_FRATIO_MIN         0.75   
     39  DC_NRATIO_BAD_MIN     0.70   
     40  DC_NRATIO_MASK_MIN    0.60
     41  DC_NRATIO_ALL_MIN     0.50
     42  DC_NGOOD_MIN          3
     43
     44  *** double-check on the use of fluxScale in pmSourcePhotometry.c:106
    145
    24620100131

  • trunk/psphot/src/psphot.h

    r26894 r27532  
    9292
    9393bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final);
    94 bool            psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final);
     94bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final);
    9595
    9696bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, bool getPSFsize);
    … …  
    119119
    120120bool            psphotMagnitudes (pmConfig *config, const pmFPAview *view);
    121 bool            psphotMagnitudesReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     121bool            psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf);
    122122bool            psphotMagnitudes_Threaded (psThreadJob *job);
    123123

  • trunk/psphot/src/psphotEfficiency.c

    r26894 r27532  
    44                    PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
    55
    6 //#define TESTING
    7 
    8 
    9 # if 0
     6# define TESTING 0
    107
    118// Calculate the limiting magnitude for an image
    … …  
    1916                     float *covarFactor,// Covariance factor
    2017                     const pmReadout *ro,     // Readout of interest
    21                      const pmPSF *psf,        // Point-spread function
     18                     pmPSF *psf,              // Point-spread function
    2219                     float thresh,            // Threshold for source identification
    2320                     float smoothSigma,       // Gaussian smoothing sigma
    … …  
    152149}
    153150
    154 # endif
    155 
    156151bool psphotEfficiency (pmConfig *config, const pmFPAview *view)
    157152{
    … …  
    178173bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
    179174{
    180 # if 0
    181175    bool status = true;
    182176
    … …  
    206200
    207201    // Collect recipe information
    208     float fwhmMajor = psMetadataLookupF32(NULL, recipe, "FWHM_MAJ"); // PSF size in x
    209     float fwhmMinor = psMetadataLookupF32(NULL, recipe, "FWHM_MIN"); // PSF size in y
     202    float smoothNsigma = psMetadataLookupF32(&status, recipe, "PEAKS_SMOOTH_NSIGMA"); // Smoothing limit
     203    psAssert (status && isfinite(smoothNsigma), "Unable to find PEAKS_SMOOTH_NSIGMA in recipe (or invalid value)");
     204
     205    float thresh = psMetadataLookupF32(&status, recipe, "PEAKS_NSIGMA_LIMIT_2");
     206    psAssert (status && isfinite(thresh), "Unable to find PEAKS_NSIGMA_LIMIT_2 in recipe (or invalid value)");
     207
     208    psVector *magOffsets = psMetadataLookupVector(&status, recipe, "EFF.MAG"); // Magnitude offsets
     209    psAssert (status, "Unable to find EFF.MAG F32 vector in recipe");
     210    psAssert (magOffsets->type.type == PS_TYPE_F32, "Unable to find EFF.MAG F32 vector in recipe");
     211
     212    int numSources = psMetadataLookupS32(&status, recipe, "EFF.NUM"); // Number of sources for each bin
     213    psAssert (status && (numSources > 0), "Unable to find EFF.NUM in recipe (or invalid value)");
     214
     215    float minGauss = psMetadataLookupF32(&status, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel
     216    psAssert (status && isfinite(minGauss), "PEAKS_MIN_GAUSS is not set in recipe (or invalid)");
     217
     218    // find the PSF size information (why is this not part of the psf structure?)
     219    float fwhmMajor = psMetadataLookupF32(NULL, readout->analysis, "FWHM_MAJ"); // PSF size in x
     220    float fwhmMinor = psMetadataLookupF32(NULL, readout->analysis, "FWHM_MIN"); // PSF size in y
    210221    if (!isfinite(fwhmMajor) || !isfinite(fwhmMinor) || fwhmMajor == 0.0 || fwhmMinor == 0.0) {
    211         psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_MAJ and FWHM_MIN in recipe");
    212         return false;
    213     }
    214     float smoothNsigma = psMetadataLookupF32(NULL, recipe, "PEAKS_SMOOTH_NSIGMA"); // Smoothing limit
    215     if (!isfinite(smoothNsigma)) {
    216         psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_SMOOTH_NSIGMA in recipe");
    217         return false;
    218     }
    219     float thresh = psMetadataLookupF32(NULL, recipe, "PEAKS_NSIGMA_LIMIT_2");
    220     if (!isfinite(thresh)) {
    221         psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_NSIGMA_LIMIT_2 in recipe");
    222         return false;
    223     }
    224     psVector *magOffsets = psMetadataLookupVector(NULL, recipe, "EFF.MAG"); // Magnitude offsets
    225     if (!magOffsets || magOffsets->type.type != PS_TYPE_F32) {
    226         psError(PSPHOT_ERR_CONFIG, false, "Unable to find EFF.MAG F32 vector in recipe");
    227         return NULL;
    228     }
    229     int numSources = psMetadataLookupS32(NULL, recipe, "EFF.NUM"); // Number of sources for each bin
    230     if (numSources == 0) {
    231         psError(PSPHOT_ERR_CONFIG, false, "Unable to find EFF.NUM in recipe");
    232         return NULL;
    233     }
    234     float minGauss = psMetadataLookupF32(NULL, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel
    235     if (!isfinite(minGauss)) {
    236         psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value");
    237         minGauss = 0.5;
     222        psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_MAJ and FWHM_MIN in readout->analysis");
     223        return false;
    238224    }
    239225
    … …  
    246232    // remove all sources, adding noise for subtracted sources
    247233    psphotRemoveAllSources(realSources, recipe);
    248     // psphotAddNoise(readout, realSources, recipe);
    249 
    250 
    251 #if defined(TESTING) && 0
     234
     235#if TESTING
    252236    {
    253237        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    … …  
    278262    }
    279263
    280 #ifdef TESTING
     264#if TESTING
    281265    psphotSaveImage(NULL, readout->image, "orig_image.fits");
    282266    psphotSaveImage(NULL, readout->variance, "orig_variance.fits");
    … …  
    293277    }
    294278
    295 #ifdef TESTING
     279#if TESTING
    296280    psphotSaveImage(NULL, readout->image, "fake_image.fits");
    297281    psphotSaveImage(NULL, readout->variance, "fake_variance.fits");
    … …  
    408392    psFree(significance);
    409393
    410     if (!psphotFitSourcesLinear(readout, fakeSourcesAll, recipe, psf, true)) {
     394    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
    411395        psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
    412396        psFree(fakeSources);
    … …  
    415399    }
    416400
    417     // Disable aperture corrections; casting away const!
     401    // Disable aperture corrections (save current value)
    418402    pmTrend2D *apTrend = psf->ApTrend;  // Aperture trend
    419     ((pmPSF*)psf)->ApTrend = NULL;
    420 
    421     if (!psphotMagnitudes(config, readout, view, fakeSourcesAll, psf)) {
     403    psf->ApTrend = NULL;
     404
     405    if (!psphotMagnitudesReadout(config, recipe, view, readout, fakeSourcesAll, psf)) {
    422406        psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources.");
    423407        psFree(fakeSources);
    424408        psFree(count);
    425         ((pmPSF*)psf)->ApTrend = apTrend; // Casting away const!
     409        psf->ApTrend = apTrend; // Casting away const!
    426410        return false;
    427411    }
    428412    psFree(fakeSourcesAll);
    429413
    430     // Re-enable aperture corrections; casting away const!
    431     ((pmPSF*)psf)->ApTrend = apTrend;
     414    // Replace aperture corrections
     415    psf->ApTrend = apTrend;
    432416
    433417    psVector *magDiffMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean difference in magnitude for each bin
    … …  
    443427        psVectorInit(magMask, 0);
    444428
    445 #ifdef TESTING
     429#if TESTING
    446430        psString name = NULL;
    447431        psStringAppend(&name, "fake_%d.dat", i);
    … …  
    459443            }
    460444
    461 #ifdef TESTING
     445#if TESTING
    462446            fprintf(file, "%f %f %f %f %f %f %f\n", source->peak->xf, source->peak->yf,
    463447                    source->modelPSF->params->data.F32[PM_PAR_XPOS],
    … …  
    495479        }
    496480
    497 #ifdef TESTING
     481#if TESTING
    498482        fclose(file);
    499483#endif
    … …  
    520504    psLogMsg("psphot", PS_LOG_INFO, "Detection efficiency: %lf sec\n", psTimerClear("psphot.fake"));
    521505
    522 # endif
    523 
    524506    return true;
    525507}

  • trunk/psphot/src/psphotFitSourcesLinear.c

    r26894 r27532  
    1717    bool status = true;
    1818
     19    // select the appropriate recipe information
     20    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     21    assert (recipe);
     22
    1923    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    2024    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    … …  
    2226    // loop over the available readouts
    2327    for (int i = 0; i < num; i++) {
    24         if (!psphotFitSourcesLinearReadout (config, view, "PSPHOT.INPUT", i, final)) {
     28
     29        // find the currently selected readout
     30        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest
     31        psAssert (file, "missing file?");
     32
     33        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     34        psAssert (readout, "missing readout?");
     35
     36        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     37        psAssert (detections, "missing detections?");
     38
     39        psArray *sources = detections->allSources;
     40        psAssert (sources, "missing sources?");
     41
     42        pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     43        psAssert (psf, "missing psf?");
     44
     45        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) {
    2546            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i);
    2647            return false;
    … …  
    3051}
    3152
    32 bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final) {
     53bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final) {
    3354
    3455    bool status;
    … …  
    3657    float y;
    3758    float f;
    38     // float r;
    39 
    40     // select the appropriate recipe information
    41     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    42     assert (recipe);
    43 
    44     // find the currently selected readout
    45     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
    46     psAssert (file, "missing file?");
    47 
    48     pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    49     psAssert (readout, "missing readout?");
    50 
    51     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    52     psAssert (detections, "missing detections?");
    53 
    54     psArray *sources = detections->allSources;
    55     psAssert (sources, "missing sources?");
    5659
    5760    if (!sources->n) {
    … …  
    5962        return true;
    6063    }
    61 
    62     pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    63     psAssert (sources, "missing psf?");
    6464
    6565    psTimerStart ("psphot.linear");

  • trunk/psphot/src/psphotFitSourcesLinearStack.c

    r25755 r27532  
    99// these are used to determine the simultaneous linear fit of fluxes.
    1010// the analysis is performed wrt the simulated pixel values
    11 
    12 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, psImageMaskType markVal);
    1311
    1412bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
    … …  
    194192    return true;
    195193}
    196 
    197 // Calculate the weight terms for the sky fit component of the matrix.  This function operates
    198 // on the pixels which correspond to all of the sources of interest.  These elements fill in
    199 // the border matrix components in the sparse matrix equation.
    200 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, psImageMaskType markVal) {
    201 
    202     // generate the image-wide weight terms
    203     // turn on MARK for all image pixels
    204     psRegion fullArray = psRegionSet (0, 0, 0, 0);
    205     fullArray = psRegionForImage (readout->mask, fullArray);
    206     psImageMaskRegion (readout->mask, fullArray, "OR", markVal);
    207 
    208     // turn off MARK for all object pixels
    209     for (int i = 0; i < sources->n; i++) {
    210         pmSource *source = sources->data[i];
    211         pmModel *model = pmSourceGetModel (NULL, source);
    212         if (model == NULL) continue;
    213         float x = model->params->data.F32[PM_PAR_XPOS];
    214         float y = model->params->data.F32[PM_PAR_YPOS];
    215         psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_IMAGE_MASK(markVal));
    216     }
    217 
    218     // accumulate the image statistics from the masked regions
    219     psF32 **image  = readout->image->data.F32;
    220     psF32 **variance = readout->variance->data.F32;
    221     psImageMaskType  **mask   = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA;
    222 
    223     double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy;
    224     w = x = y = x2 = xy = y2 = fo = fx = fy = 0;
    225 
    226     int col0 = readout->image->col0;
    227     int row0 = readout->image->row0;
    228 
    229     for (int j = 0; j < readout->image->numRows; j++) {
    230         for (int i = 0; i < readout->image->numCols; i++) {
    231             if (mask[j][i]) continue;
    232             if (constant_weights) {
    233                 wt = 1.0;
    234             } else {
    235                 wt = variance[j][i];
    236             }
    237             f = image[j][i];
    238             w   += 1/wt;
    239             fo  += f/wt;
    240         }
    241     }
    242 
    243     // turn off MARK for all image pixels
    244     psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_IMAGE_MASK(markVal));
    245 
    246     // set the Border T elements
    247     psSparseBorderElementG (border, 0, fo);
    248     psSparseBorderElementT (border, 0, 0, w);
    249 
    250     return true;
    251 }

  • trunk/psphot/src/psphotKernelFromPSF.c

    r17396 r27532  
    44
    55    assert (source);
    6     assert (source->psfFlux); // XXX build if needed?
     6    assert (source->psfImage); // XXX build if needed?
    77
    8     int x0 = source->peak->xf - source->psfFlux->col0;
    9     int y0 = source->peak->yf - source->psfFlux->row0;
     8    int x0 = source->peak->xf - source->psfImage->col0;
     9    int y0 = source->peak->yf - source->psfImage->row0;
    1010
    1111    // need to decide on the size: dynamically? statically?
    … …  
    1717    // if the realized PSF for this object does not cover the full kernel, give up for now
    1818    if (x0 + psf->xMin < 0) goto escape;
    19     if (x0 + psf->xMax >= source->psfFlux->numCols) goto escape;
     19    if (x0 + psf->xMax >= source->psfImage->numCols) goto escape;
    2020    if (y0 + psf->yMin < 0) goto escape;
    21     if (y0 + psf->yMax >= source->psfFlux->numRows) goto escape;
     21    if (y0 + psf->yMax >= source->psfImage->numRows) goto escape;
    2222
    2323    double sum = 0.0;
    2424    for (int j = psf->yMin; j <= psf->yMax; j++) {
    2525        for (int i = psf->xMin; i <= psf->xMax; i++) {
    26             double value = source->psfFlux->data.F32[y0 + j][x0 + i];
     26            double value = source->psfImage->data.F32[y0 + j][x0 + i];
    2727            psf->kernel[j][i] = value;
    2828            sum += value;

  • trunk/psphot/src/psphotMagnitudes.c

    r26894 r27532  
    1414    // loop over the available readouts
    1515    for (int i = 0; i < num; i++) {
    16         if (!psphotMagnitudesReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     16
     17        // find the currently selected readout
     18        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest
     19        psAssert (file, "missing file?");
     20
     21        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     22        psAssert (readout, "missing readout?");
     23
     24        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     25        psAssert (detections, "missing detections?");
     26
     27        psArray *sources = detections->allSources;
     28        psAssert (sources, "missing sources?");
     29
     30        pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     31        psAssert (psf, "missing psf?");
     32
     33        if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) {
    1734            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for PSPHOT.INPUT entry %d", i);
    1835            return false;
    … …  
    2239}
    2340
    24 bool psphotMagnitudesReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
     41bool psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
    2542
    2643    bool status = false;
    2744    int Nap = 0;
    2845
     46    if (!sources->n) {
     47        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source magnitudes");
     48        return true;
     49    }
     50       
    2951    psTimerStart ("psphot.mags");
    30 
    31     // find the currently selected readout
    32     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
    33     psAssert (file, "missing file?");
    34 
    35     pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    36     psAssert (readout, "missing readout?");
    37 
    38     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    39     psAssert (detections, "missing detections?");
    40 
    41     psArray *sources = detections->allSources;
    42     psAssert (sources, "missing sources?");
    43 
    44     if (!sources->n) {
    45         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
    46         return true;
    47     }
    48 
    49     pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    50     psAssert (psf, "missing psf?");
    5152
    5253    // determine the number of allowed threads
    … …  
    7980    }
    8081
    81     bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
     82    bool IGNORE_GROWTH  = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
    8283    bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
     84    bool DIFF_STATS     = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
    8385
    8486    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_APCORR | PM_SOURCE_PHOT_WEIGHT;
    8587    if (!IGNORE_GROWTH) photMode |= PM_SOURCE_PHOT_GROWTH;
    8688    if (INTERPOLATE_AP) photMode |= PM_SOURCE_PHOT_INTERP;
     89    if (DIFF_STATS)     photMode |= PM_SOURCE_PHOT_DIFFSTATS;
    8790
    8891    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)

  • trunk/psphot/src/psphotMakeResiduals.c

    r25755 r27532  
    305305                if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) {
    306306                  resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
     307                  resid->Ro->data.F32[oy][ox] = 0.0;
     308                  resid->Rx->data.F32[oy][ox] = 0.0;
     309                  resid->Ry->data.F32[oy][ox] = 0.0;
    307310                }
    308311            }

  • trunk/psphot/src/psphotReadout.c

    r26894 r27532  
    9191        return psphotReadoutCleanup (config, view);
    9292    }
    93     if (!strcasecmp (breakPt, "MOMENTS")) {
    94         return psphotReadoutCleanup(config, view);
    95     }
    96 
    9793    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
    9894    if (!psphotImageQuality (config, view)) { // pass 1
    9995        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
     96        return psphotReadoutCleanup(config, view);
     97    }
     98    if (!strcasecmp (breakPt, "MOMENTS")) {
    10099        return psphotReadoutCleanup(config, view);
    101100    }

  • trunk/psphot/src/psphotSourceSize.c

    r26894 r27532  
    1313    float soft;
    1414    int grow;
     15    int xtest, ytest;
     16    bool apply; // apply CR mask?
    1517} psphotSourceSizeOptions;
    1618
    … …  
    2224bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal);
    2325bool psphotMaskCosmicRayFootprintCheck (psArray *sources);
     26int  psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh);
    2427
    2528// we need to call this function after sources have been fitted to the PSF model and
    … …  
    101104    assert (status);
    102105
     106    // XXX recipe name is not great
     107    options.xtest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.XTEST");
     108    options.ytest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.YTEST");
     109    assert (status);
     110
    103111    options.grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs
    104112    if (!status || options.grow < 0) {
    … …  
    111119        psWarning("PSPHOT.CR.NSIGMA.SOFTEN not set; defaulting to zero.");
    112120        options.soft = 0.0;
     121    }
     122
     123    options.apply = psMetadataLookupBool(&status, recipe, "PSPHOT.CRMASK.APPLY"); // Growth size for CRs
     124    if (!status) {
     125        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CRMASK.APPLY is not defined.");
     126        return false;
    113127    }
    114128
    … …  
    239253            continue;
    240254        }
     255        // psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
    241256    }
    242257
    … …  
    244259}
    245260
     261# define SIZE_SN_LIM 10
    246262bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
    247263
    … …  
    285301        }
    286302
    287         // we are basically classifying by moments; use the default if not found
     303        // we are basically classifying by moments
    288304        psAssert (source->moments, "why is this source missing moments?");
    289305        if (source->mode & noMoments) {
    … …  
    292308        }
    293309
     310        // convert to Mmaj, Mmin:
    294311        psF32 Mxx = source->moments->Mxx;
    295312        psF32 Myy = source->moments->Myy;
    … …  
    315332        float apMag = -2.5*log10(source->moments->Sum);
    316333        float dMag = source->psfMag - apMag;
    317         float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
    318 
    319         source->extNsigma = nSigma;
    320         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     334
     335        // set nSigma to include both systematic and poisson error terms
     336        // XXX the 'poisson error' contribution for size is probably wrong...
     337        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
     338        float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag);
     339        float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag);
     340
     341        // partially-masked sources are more likely to be mis-measured PSFs
     342        float sizeBias = 1.0;
     343        if (source->pixWeight < 0.9) {
     344            sizeBias = 3.0;
     345        }
     346
     347        float minMxx = psfClump->X - sizeBias*options->nSigmaMoments*psfClump->dX;
     348        float minMyy = psfClump->Y - sizeBias*options->nSigmaMoments*psfClump->dY;
     349
     350        // include MAG, MXX, and MYY?
     351        source->extNsigma = nSigmaMAG;
     352
     353        // notes to clarify the source size classification rules:
     354        // * a defect should be functionally equivalent to a cosmic ray
     355        // * CR & defect should have a faintess limit (min S/N)
     356        // * SAT stars should not be faint, but defects may?
    321357
    322358        // Anything within this region is a probably PSF-like object. Saturated stars may land
    323359        // in this region, but are detected elsewhere on the basis of their peak value.
    324         bool isPSF = ((fabs(nSigma) < options->nSigmaApResid) &&
    325                       (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) &&
    326                       (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY));
     360        bool isPSF = (fabs(nSigmaMAG) < options->nSigmaApResid) && (fabs(nSigmaMXX) < sizeBias*options->nSigmaMoments) && (fabs(nSigmaMYY) < sizeBias*options->nSigmaMoments);
    327361        if (isPSF) {
    328           psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
    329                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    330                   options->nSigmaApResid,options->nSigmaMoments);
     362          psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
     363                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     364                  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
     365          source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    331366          Npsf ++;
    332367          continue;
    … …  
    336371        // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
    337372        // XXX this rule is not great
    338         if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
    339 
    340           psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
    341                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    342                   options->nSigmaApResid,options->nSigmaMoments);
     373        // XXX only accept brightish detections as CRs
     374        // (nSigmaMAG < -options->nSigmaApResid) ||
     375        bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy));
     376        if (isCR) {
     377            psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g %f\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
     378                    source->peak->xf,source->peak->yf,source->pixWeight,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     379                    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    343380            source->mode |= PM_SOURCE_MODE_DEFECT;
     381            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
    344382            Ncr ++;
    345383            continue;
    … …  
    349387        // just large saturated regions.
    350388        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    351           psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
    352                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    353                   options->nSigmaApResid,options->nSigmaMoments);
    354          
     389            psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
     390                    source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     391                    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
     392            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    355393            Nsat ++;
    356394            continue;
    … …  
    358396
    359397        // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
    360         bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
     398        bool isEXT = (nSigmaMAG > options->nSigmaApResid) || (Mxx > minMxx) || (Myy > minMyy);
    361399        if (isEXT) {
    362           psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
    363                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    364                   options->nSigmaApResid,options->nSigmaMoments);
     400          psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
     401                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     402                  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    365403         
    366404            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
     405            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    367406            Next ++;
    368407            continue;
    369408        }
    370         psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
    371                 source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    372                 options->nSigmaApResid,options->nSigmaMoments);
     409        psTrace("psphotSourceClassRegion.MISS",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
     410                source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     411                options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    373412       
    374         psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
     413        // sources that reach here are probably too faint for a reasonable source size measurement
     414        // psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigmaMAG);
     415        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    375416        Nmiss ++;
    376417    }
    … …  
    385426bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
    386427
     428    psTimerStart ("psphot.cr");
     429
     430    int nMasked = 0;
    387431    for (int i = 0; i < sources->n; i++) {
    388432        pmSource *source = sources->data[i];
    … …  
    393437            continue;
    394438        }
     439
     440        // only check candidates marked above
     441        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) {
     442            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
     443            continue;
     444        }
     445
     446        // skip unless this source is thought to be a cosmic ray.  flag the detection and mask the pixels
     447        // XXX this may be degenerate with the above test
     448        if (!(source->mode & PM_SOURCE_MODE_DEFECT)) continue;
    395449
    396450        // Integer position of peak
    … …  
    402456            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
    403457            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
    404             //      psTrace("psphot.czw", 2, "Not calculating crNsigma due to edge\n");
    405458            continue;
    406459        }
    407460       
    408         // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    409         if (source->mode & PM_SOURCE_MODE_DEFECT) {
    410             // XXX this is running slowly and is too agressive, but it more-or-less works
    411             // psphotMaskCosmicRayCZW(readout, source, options->crMask);
    412            
    413             // XXX these are the old versions which we are not using
    414             // psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
    415             // psphotMaskCosmicRay_Old (source, maskVal, crMask);
    416         }
     461        // XXX for testing, only CRMASK a single source:
     462        if (options->xtest && (fabs(source->peak->xf - options->xtest) > 5)) continue;
     463        if (options->ytest && (fabs(source->peak->yf - options->ytest) > 5)) continue;
     464
     465        // replace object in image
     466        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     467            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
     468        }
     469
     470        // XXX this is running slowly and is too agressive, but it more-or-less works
     471        psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf);
     472        if (options->apply) {
     473            psphotMaskCosmicRay(readout, source, options->crMask);
     474        } else {
     475            source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     476        }
     477        nMasked ++;
     478
     479        // re-subtract the object, leave local sky
     480        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
    417481    }
    418482   
    … …  
    430494    }
    431495
     496    psLogMsg ("psphot.cr", PS_LOG_INFO, "mask CR: %d masked in %f sec\n", nMasked, psTimerMark ("psphot.cr"));
     497
    432498    // XXX test : save the mask image
    433499    if (0) {
    434500        psphotSaveImage (NULL, readout->mask,   "mask.fits");
    435501    }
     502
    436503    return true;
    437504}
     505
     506# define DUMPPICS 0
     507# define LIMIT_XRANGE(X, IMAGE) { X = PS_MIN(PS_MAX(0, X), IMAGE->numCols); }
     508# define LIMIT_YRANGE(Y, IMAGE) { Y = PS_MIN(PS_MAX(0, Y), IMAGE->numRows); }
    438509
    439510// Comments by CZW 20091209 : Mechanics of how to identify CR pixels taken from "Cosmic-Ray
    … …  
    448519    pmFootprint *footprint = peak->footprint;
    449520
    450     int xm = footprint->bbox.x0;
    451     int xM = footprint->bbox.x1;
    452     int ym = footprint->bbox.y0;
    453     int yM = footprint->bbox.y1;
    454 
    455     if (xm < 0) { xm = 0; }
    456     if (ym < 0) { ym = 0; }
    457     if (xM > mask->numCols) { xM = mask->numCols; }
    458     if (yM > mask->numRows) { yM = mask->numRows; }
    459     int dx = xM - xm;
    460     int dy = yM - ym;
    461 
    462521    // Bounding boxes are inclusive of final pixel
    463     // XXX ensure dx,dy do not become too large here
    464     dx++;
    465     dy++;
     522    int xs = footprint->bbox.x0;
     523    int xe = footprint->bbox.x1 + 1;
     524    int ys = footprint->bbox.y0;
     525    int ye = footprint->bbox.y1 + 1;
     526
     527    LIMIT_XRANGE(xs, mask);
     528    LIMIT_XRANGE(xe, mask);
     529    LIMIT_YRANGE(ys, mask);
     530    LIMIT_YRANGE(ye, mask);
     531
     532    int dx = xe - xs;
     533    int dy = ye - ys;
    466534
    467535    psImage *image= readout->image;
    468536    psImage *variance = readout->variance;
    469537     
    470     int binning = 1;
    471     float sigma_thresh = 2.0;
    472     int iteration = 0;
    473     int max_iter = 5;
    474     float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
     538    int binning = 2;
     539    float sigma_thresh = 3.0;
     540    int max_iter = 1; // XXX with isophot masking, we only want to do a single pass
    475541
    476542    // Temporary images.
    477543    psImage *mypix  = psImageAlloc(dx,dy,image->type.type);
     544    psImage *myfix  = psImageAlloc(dx,dy,image->type.type);
    478545    psImage *myvar  = psImageAlloc(dx,dy,image->type.type);
    479546    psImage *binned = psImageAlloc(dx * binning,dy * binning,image->type.type);
    … …  
    482549    psImage *mymask = psImageAlloc(dx,dy,PS_TYPE_IMAGE_MASK);
    483550     
    484     int x,y;
    485551    // Load my copy of things.
    486     for (x = 0; x < dx; x++) {
    487         for (y = 0; y < dy; y++) {
    488             psImageSet(mypix,x,y,psImageGet(image,x+xm,y+ym));
    489             psImageSet(myvar,x,y,psImageGet(variance,x+xm,y+ym));
     552    for (int y = 0; y < dy; y++) {
     553        for (int x = 0; x < dx; x++) {
     554            mypix->data.F32[y][x] = image->data.F32[y+ys][x+xs];
     555            myvar->data.F32[y][x] = variance->data.F32[y+ys][x+xs];
    490556            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00;
    491557        }
    … …  
    495561        pmSpan *sp = footprint->spans->data[i];
    496562        for (int j = sp->x0; j <= sp->x1; j++) {
    497             y = sp->y - ym;
    498             x = j - xm;
     563            int y = sp->y - ys;
     564            int x = j - xs;
    499565            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01;
    500566        }
    501567    }
    502568
    503     int CRpix_count = 0;     
    504     do {
    505         CRpix_count = 0;
    506         // Zero out my temp images.
    507         for (x = 0; x < binning * dx; x++) {
    508             for (y = 0; y < binning * dy; y++) {
    509                 psImageSet(binned,x,y,0.0);
    510                 psImageSet(conved,x,y,0.0);
    511                 psImageSet(edges,x/binning,y/ binning,0.0);
     569    int nCRpix = 1; // force at least one pass...
     570    for (int iteration = 0; (iteration < max_iter) && (nCRpix > 0); iteration++) {
     571        nCRpix = 0;
     572        psImageInit (binned, 0.0);
     573        psImageInit (conved, 0.0);
     574        psImageInit (edges, 0.0);
     575
     576        // Make subsampled image. Maybe this should be called "unbinned" or something
     577        for (int y = 0; y < binning * dy; y++) {
     578            int yraw = y / binning;
     579            for (int x = 0; x < binning * dx; x++) {
     580                int xraw = x / binning;
     581                binned->data.F32[y][x] = mypix->data.F32[yraw][xraw];
    512582            }
    513         }     
    514         // Make subsampled image. Maybe this should be called "unbinned" or something
    515         for (x = 0; x < binning * dx; x++) {
    516             for (y = 0; y < binning * dy; y++) {
    517                 psImageSet(binned,x,y,psImageGet(mypix,x / binning,y / binning));
     583        }
     584
     585        // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
     586        for (int y = 1; y < binning * dy - 1; y++) {
     587            for (int x = 1; x < binning * dx - 1; x++) {
     588                float value = binned->data.F32[y][x] - 0.25 *
     589                    (binned->data.F32[y+0][x-1] + binned->data.F32[y+0][x+1] +
     590                     binned->data.F32[y-1][x+0] + binned->data.F32[y+1][x+0]);
     591                value = PS_MAX(0.0, value);
     592
     593                conved->data.F32[y][x] = value;
    518594            }
    519595        }
    520         // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
    521         for (x = 1; x < dx - 1; x++) {
    522             for (y = 1; y < dy - 1; y++) {
    523                 psImageSet(conved,x,y,psImageGet(binned,x,y) - 0.25 *
    524                            (psImageGet(binned,x-1,y) + psImageGet(binned,x+1,y) +
    525                             psImageGet(binned,x,y-1) + psImageGet(binned,x,y+1)));
    526                 if (psImageGet(conved,x,y) < 0.0) {
    527                     psImageSet(conved,x,y,0.0);
    528                 }
     596
     597        // Create an edge map by rebinning
     598        for (int y = 0; y < binning * dy; y++) {
     599            int yraw = y / binning;
     600            for (int x = 0; x < binning * dx; x++) {
     601                int xraw = x / binning;
     602                edges->data.F32[yraw][xraw] += conved->data.F32[y][x];
    529603            }
    530604        }
    531         // Create an edge map by rebinning
    532         for (x = 0; x < binning * dx; x++) {
    533             for (y = 0; y < binning * dy; y++) {
    534                 psImageSet(edges,x / binning, y / binning,
    535                            psImageGet(edges, x / binning, y / binning) +
    536                            psImageGet(conved,x,y));
    537             }
    538         }
    539         // Modify my mask if we're above the significance threshold
    540         for (x = 0; x < dx; x++) {
    541             for (y = 0; y < dy; y++) {
    542                 if ( psImageGet(edges,x,y) / (binning * sqrt(noise_factor * psImageGet(myvar,x,y))) > sigma_thresh ) {
    543                     if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01) {
    544                         mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x40;
    545                         CRpix_count++;
    546                     }
    547                 }
    548             }
    549         }
    550 
     605
     606        // coordinate of peak in subimage pixels:
     607        int xPeak = peak->x - xs;
     608        int yPeak = peak->y - ys;
     609
     610        // Modify my mask if we're above the significance threshold, but only for connected pixels
     611        nCRpix = psphotMaskCosmicRayConnected (xPeak, yPeak, mymask, myvar, edges, binning, sigma_thresh);
     612       
     613# if DUMPPICS
     614        psphotSaveImage (NULL, mypix,   "crmask.pix.fits");
     615# endif
     616       
     617// XXX do not repair the pixels in isophot version
     618# if 0
    551619        // "Repair" Masked pixels for the next round.
    552         for (x = 1; x < dx - 1; x++) {
    553             for (y = 1; y < dy - 1; y++) {
    554                 if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
    555                     psImageSet(mypix,x,y,0.25 *
    556                                (psImageGet(mypix,x-1,y) + psImageGet(mypix,x+1,y) +
    557                                 psImageGet(mypix,x,y-1) + psImageGet(mypix,x,y+1)));
    558                 }
    559             }
    560         }
    561        
    562 # if 0
    563         if ((psTraceGetLevel("psphot.czw") >= 2)&&(abs(xm - 2770) < 5)&&(abs(ym - 2581) < 5)&&(iteration == 0)) {
    564           psTrace("psphot.czw",2,"TRACEMOTRON %d %d %d %d %d\n",xm,ym,dx,dy,iteration);
    565           psphotSaveImage (NULL, mypix,   "czw.pix.fits");
    566           psphotSaveImage (NULL, myvar,   "czw.var.fits");
    567           psphotSaveImage (NULL, binned,  "czw.binn.fits");
    568           psphotSaveImage (NULL, conved,  "czw.conv.fits");
    569           psphotSaveImage (NULL, edges,   "czw.edge.fits");
    570           psphotSaveImage (NULL, mymask,  "czw.mask.fits");
    571         }
    572 # endif
    573 
    574         psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration,CRpix_count);
    575         iteration++;
    576     } while ((iteration < max_iter)&&(CRpix_count > 0));
    577 
    578     // A solitary masked pixel is likely a lie. Remove those.
    579     for (x = 0; x < dx; x++) {
    580         for (y = 0; y < dy; y++) {
    581             if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
    582                 if ((x-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
     620        for (int y = 1; y < dy - 1; y++) {
     621            for (int x = 1; x < dx - 1; x++) {
     622                if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) {
     623                    myfix->data.F32[y][x] = mypix->data.F32[y][x];
    583624                    continue;
    584625                }
    585                 if ((y-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
    586                     continue;
    587                 }
    588                 if ((x+1 < dx)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
    589                     continue;
    590                 }
    591                 if ((y+1 < dy)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
    592                     continue;
    593                 }
    594                 mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
     626                myfix->data.F32[y][x] = 0.25 *
     627                    (mypix->data.F32[y+0][x-1] + mypix->data.F32[y+0][x+1] +
     628                     mypix->data.F32[y-1][x+0] + mypix->data.F32[y+1][x+0]);
    595629            }
    596630        }
    597     }
    598 
    599     // transfer temporary mask to real mask
    600     for (x = 0; x < dx; x++) {
    601         for (y = 0; y < dy; y++) {
     631       
     632        // "Repair" Masked pixels for the next round.
     633        for (int y = 1; y < dy - 1; y++) {
     634            for (int x = 1; x < dx - 1; x++) {
     635                mypix->data.F32[y][x] = myfix->data.F32[y][x];
     636            }
     637        }
     638# endif
     639
     640# if DUMPPICS
     641        fprintf (stderr, "CRMASK %d %d %d %d %d\n", xs, ys, dx, dy, iteration);
     642        psphotSaveImage (NULL, mypix,   "crmask.fix.fits");
     643        psphotSaveImage (NULL, myvar,   "crmask.var.fits");
     644        psphotSaveImage (NULL, binned,  "crmask.binn.fits");
     645        psphotSaveImage (NULL, conved,  "crmask.conv.fits");
     646        psphotSaveImage (NULL, edges,   "crmask.edge.fits");
     647        psphotSaveImage (NULL, mymask,  "crmask.mask.fits");
     648# endif
     649        psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration, nCRpix);
     650    }
     651
     652# if 0
     653    // A solitary masked pixel is likely a lie. Remove those
     654    // XXX can't we use nCRpix == 1 to test for these?
     655    for (int x = 0; x < dx; x++) {
     656        for (int y = 0; y < dy; y++) {
     657            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) continue;
     658            if ((x-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
     659                continue;
     660            }
     661            if ((y-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
     662                continue;
     663            }
     664            if ((x+1 < dx) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
     665                continue;
     666            }
     667            if ((y+1 < dy) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
     668                continue;
     669            }
     670            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
     671        }
     672    }
     673# endif
     674
     675    // transfer temporary mask to real mask & count masked pixels
     676    nCRpix = 0;
     677    for (int x = 0; x < dx; x++) {
     678        for (int y = 0; y < dy; y++) {
    602679            if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
    603                 mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ym+mask->row0][x+xm+mask->col0] |= maskVal;
     680                mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ys+mask->row0][x+xs+mask->col0] |= maskVal;
     681                nCRpix ++;
    604682            }
    605683        }
    … …  
    607685     
    608686    // XXX if we decide this REALLY is a cosmic ray, set the CR_LIMIT bit
    609     source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     687    if (nCRpix > 1) {
     688        source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     689        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     690    }
     691    // fprintf (stderr, "CRMASK %d %d %d %d %d\n", peak->x, peak->y, dx, dy, nCRpix);
    610692
    611693    psFree(mypix);
     694    psFree(myfix);
    612695    psFree(myvar);
    613696    psFree(binned);
    … …  
    680763    }
    681764    return true;
     765}
     766
     767# define VERBOSE 0
     768int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh) {
     769   
     770    int xLo, xRo;
     771    int nCRpix = 0;
     772
     773    float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
     774   
     775    // mark the pixels in this row to the left, then the right. stay within footprint
     776    int xL = xPeak; // find the range of valid pixels in this row
     777    int xR = xPeak;
     778    for (int ix = xPeak; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix--) {
     779        float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
     780        float value = edges->data.F32[yPeak][ix] / noise;
     781        if (value < sigma_thresh ) break;
     782        mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
     783        xL = ix;
     784        nCRpix ++;
     785        if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
     786    }
     787    for (int ix = xPeak; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix++) {
     788        float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
     789        float value = edges->data.F32[yPeak][ix] / noise;
     790        if (value < sigma_thresh ) break;
     791        mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
     792        xR = ix;
     793        nCRpix ++;
     794        if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
     795    }
     796    // xL and xR mark the first and last valid pixel in the row
     797
     798    // for each of the neighboring rows, mark the high pixels if they touch the range xL to xR
     799    xLo = PS_MAX(xL - 1, 0);
     800    xRo = PS_MIN(xR + 1, mymask->numCols);
     801
     802    // first go down:
     803    for (int iy = yPeak - 1; iy >= 0; iy--) {
     804
     805        int xLn = -1;
     806        int xRn = -1;
     807        int newPix = 0;
     808
     809        // mark the pixels in the good range
     810        for (int ix = xLo; ix < xRo; ix++) {
     811            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
     812            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
     813            float value = edges->data.F32[iy][ix] / noise;
     814            if (value < sigma_thresh ) continue;
     815            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
     816            if (xLn == -1) xLn = ix; // first valid pixel in this row
     817            xRn = ix;                // last valid pixel in this row
     818            nCRpix ++;
     819            newPix ++;
     820            if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
     821        }
     822
     823        // mark the pixels to the left of the good range
     824        for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
     825            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
     826            float value = edges->data.F32[iy][ix] / noise;
     827            if (value < sigma_thresh ) break;
     828            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
     829            if (xRn == -1) xRn = ix; // last valid pixel in this row
     830            xLn = ix;
     831            nCRpix ++;
     832            newPix ++;
     833            if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
     834        }
     835
     836        // mark the pixels to the right of the good range
     837        for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
     838            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
     839            float value = edges->data.F32[iy][ix] / noise;
     840            if (value < sigma_thresh ) break;
     841            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
     842            if (xLn == -1) xLn = ix; // first valid pixel in this row
     843            xRn = ix;
     844            nCRpix ++;
     845            newPix ++;
     846            if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
     847        }
     848        if (newPix == 0) break;
     849        xLo = PS_MAX(xLn - 1, 0);
     850        xRo = PS_MIN(xRn + 1, mymask->numCols);
     851    }
     852
     853    xLo = PS_MAX(xL - 1, 0);
     854    xRo = PS_MIN(xR + 1, mymask->numCols);
     855
     856    // next go up:
     857    for (int iy = yPeak + 1; iy < mymask->numRows; iy++) {
     858
     859        int xLn = -1;
     860        int xRn = -1;
     861        int newPix = 0;
     862
     863        // mark the pixels in the good range
     864        for (int ix = xLo; ix < xRo; ix++) {
     865            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
     866            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
     867            float value = edges->data.F32[iy][ix] / noise;
     868            if (value < sigma_thresh ) continue;
     869            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
     870            if (xLn == -1) xLn = ix; // first valid pixel in this row
     871            xRn = ix;                // last valid pixel in this row
     872            nCRpix ++;
     873            newPix ++;
     874            if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
     875        }
     876
     877        // mark the pixels to the left of the good range
     878        for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
     879            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
     880            float value = edges->data.F32[iy][ix] / noise;
     881            if (value < sigma_thresh ) break;
     882            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
     883            if (xRn == -1) xRn = ix; // last valid pixel in this row
     884            xLn = ix;
     885            nCRpix ++;
     886            newPix ++;
     887            if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
     888        }
     889
     890        // mark the pixels to the right of the good range
     891        for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
     892            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
     893            float value = edges->data.F32[iy][ix] / noise;
     894            if (value < sigma_thresh ) break;
     895            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
     896            if (xLn == -1) xLn = ix; // first valid pixel in this row
     897            xRn = ix;
     898            nCRpix ++;
     899            newPix ++;
     900            if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
     901        }
     902        if (newPix == 0) break;
     903        xLo = PS_MAX(xLn - 1, 0);
     904        xRo = PS_MIN(xRn + 1, mymask->numCols);
     905    }
     906
     907    return nCRpix;
    682908}
    683909

  • trunk/psphot/src/psphotVisual.c

    r26894 r27532  
    10991099    psFree (outsat);
    11001100    return true;
     1101}
     1102
     1103static void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1)
     1104{
     1105    float x[2], y[2];
     1106    x[0] = x0;
     1107    x[1] = x1;
     1108    y[0] = y0;
     1109    y[1] = y1;
     1110    KapaPrepPlot   (myKapa, 2, graphdata);
     1111    KapaPlotVector (myKapa, 2, x, "x");
     1112    KapaPlotVector (myKapa, 2, y, "y");
    11011113}
    11021114
    … …  
    11381150    }
    11391151
     1152    // generate model profiles (major and minor axis):
     1153    // create a model with theta = 0.0 so major and minor axes are equiv to x and y:
     1154    psEllipseShape rawShape, rotShape;
     1155
     1156    rawShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
     1157    rawShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
     1158    rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
     1159
     1160    psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0);
     1161
     1162    axes.theta = 0.0;
     1163
     1164    rotShape = psEllipseAxesToShape (axes);
     1165
     1166    psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32);
     1167    for (int i = 0; i < source->modelPSF->params->n; i++) {
     1168        params->data.F32[i] = source->modelPSF->params->data.F32[i];
     1169    }
     1170    params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2;
     1171    params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2;
     1172    params->data.F32[PM_PAR_SXY] = rotShape.sxy;
     1173    params->data.F32[PM_PAR_XPOS] = 0.0;
     1174    params->data.F32[PM_PAR_YPOS] = 0.0;
     1175
     1176    psVector *rmod = psVectorAlloc(300, PS_TYPE_F32);
     1177    psVector *fmaj = psVectorAlloc(300, PS_TYPE_F32);
     1178    psVector *fmin = psVectorAlloc(300, PS_TYPE_F32);
     1179
     1180    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     1181
     1182    float r = 0.0;
     1183    for (int i = 0; i < rmod->n; i++) {
     1184        r = i*0.1;
     1185        rmod->data.F32[i] = r;
     1186
     1187        coord->data.F32[1] = r;
     1188        coord->data.F32[0] = 0.0;
     1189        fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
     1190
     1191        coord->data.F32[0] = r;
     1192        coord->data.F32[1] = 0.0;
     1193        fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
     1194    }
     1195    psFree (coord);
     1196    psFree (params);
     1197
     1198    float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
     1199    float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
     1200    if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR);
     1201
     1202    psEllipseMoments emoments;
     1203    emoments.x2 = source->moments->Mxx;
     1204    emoments.xy = source->moments->Mxy;
     1205    emoments.y2 = source->moments->Myy;
     1206    axes = psEllipseMomentsToAxes (emoments, 20.0);
     1207    float MOMENTS_MAJOR = 2.355*axes.major;
     1208    float MOMENTS_MINOR = 2.355*axes.minor;
     1209
     1210    float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
     1211
    11401212    // reset source Add/Sub state to recorded
    11411213    if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    … …  
    11741246    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
    11751247
     1248    graphdata.color = KapaColorByName ("blue");
     1249    graphdata.ptype = 0;
     1250    graphdata.size = 0.0;
     1251    graphdata.style = 0;
     1252    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1253    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1254    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
     1255    plotline (myKapa, &graphdata, 0.0, logHM, 30.0, logHM);
     1256    plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, 0.0, 0.5*FWHM_MINOR, 5.0);
     1257    graphdata.ltype = 1;
     1258    plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, 0.0, 0.5*MOMENTS_MINOR, 5.0);
     1259    graphdata.ltype = 0;
     1260       
     1261    graphdata.color = KapaColorByName ("green");
     1262    graphdata.ptype = 0;
     1263    graphdata.size = 0.0;
     1264    graphdata.style = 0;
     1265    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1266    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1267    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
     1268    plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, 0.0, 0.5*FWHM_MAJOR, 5.0);
     1269    graphdata.ltype = 1;
     1270    plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, 0.0, 0.5*MOMENTS_MAJOR, 5.0);
     1271    graphdata.ltype = 0;
     1272       
     1273    for (int i = 0; i < rmod->n; i++) {
     1274        rmod->data.F32[i] = log10(rmod->data.F32[i]);
     1275    }
     1276
    11761277    // ** loglog **
    11771278    KapaSelectSection (myKapa, "loglog");
    … …  
    11821283    graphdata.ymin = -0.05;
    11831284    graphdata.ymax = +5.05;
     1285    graphdata.color = KapaColorByName ("black");
    11841286    KapaSetLimits (myKapa, &graphdata);
    11851287
    … …  
    12041306    KapaPlotVector (myKapa, nb, Rb->data.F32, "x");
    12051307    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
     1308
     1309    graphdata.color = KapaColorByName ("blue");
     1310    graphdata.ptype = 0;
     1311    graphdata.size = 0.0;
     1312    graphdata.style = 0;
     1313    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1314    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1315    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
     1316
     1317    graphdata.color = KapaColorByName ("green");
     1318    graphdata.ptype = 0;
     1319    graphdata.size = 0.0;
     1320    graphdata.style = 0;
     1321    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1322    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1323    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
     1324
     1325    psFree (rmod);
     1326    psFree (fmin);
     1327    psFree (fmaj);
    12061328
    12071329    psFree (rg);
    … …  
    13921514        if (source == NULL) continue;
    13931515
    1394         // if (source->type != type) continue;
    13951516        if (mode) {
    13961517            if (keep) {
    … …  
    15161637    psVector *sDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    15171638
     1639    psVector *xLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1640    psVector *yLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1641    psVector *mLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1642    psVector *sLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1643
    15181644    psVector *xCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    15191645    psVector *yCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    … …  
    15261652    int nPSF = 0;
    15271653    int nDEF = 0;
     1654    int nLOW = 0;
    15281655    int nCR  = 0;
    15291656    for (int i = 0; i < sources->n; i++) {
    15301657        pmSource *source = sources->data[i];
    15311658        if (source->moments == NULL) continue;
     1659
     1660        // only plot the measured sources...
     1661        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) continue;
    15321662
    15331663        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
    … …  
    15611691            continue;
    15621692        }
    1563         if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
     1693        if (source->errMag > 0.1) {
     1694            xLOW->data.F32[nLOW] = source->moments->Mxx;
     1695            yLOW->data.F32[nLOW] = source->moments->Myy;
     1696            mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum);
     1697            sLOW->data.F32[nLOW] = source->extNsigma;
     1698            nLOW++;
    15641699            continue;
    15651700        }
    … …  
    15701705        nPSF++;
    15711706    }
     1707
    15721708    xSAT->n = nSAT;
    15731709    ySAT->n = nSAT;
    … …  
    15941730    mDEF->n = nDEF;
    15951731    sDEF->n = nDEF;
     1732
     1733    xLOW->n = nLOW;
     1734    yLOW->n = nLOW;
     1735    mLOW->n = nLOW;
     1736    sLOW->n = nLOW;
    15961737
    15971738    // four sections: MxxMyy, MagMxx, MagMyy, MagSigma
    … …  
    16571798    KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y");
    16581799
     1800    graphdata.color = KapaColorByName ("black");
     1801    graphdata.ptype = 7;
     1802    graphdata.size = 1.0;
     1803    graphdata.style = 2;
     1804    KapaPrepPlot   (myKapa, nLOW, &graphdata);
     1805    KapaPlotVector (myKapa, nLOW, xLOW->data.F32, "x");
     1806    KapaPlotVector (myKapa, nLOW, yLOW->data.F32, "y");
     1807
    16591808    // second section: MagMyy
    16601809    section.dx = 0.75;
    … …  
    16681817    graphdata.color = KapaColorByName ("black");
    16691818    graphdata.xmin = -17.1;
    1670     graphdata.xmax =  -7.9;
     1819    graphdata.xmax =  -6.9;
    16711820    graphdata.ymin = Ymin;
    16721821    graphdata.ymax = Ymax;
    … …  
    17301879    graphdata.xmin = Xmin;
    17311880    graphdata.xmax = Xmax;
    1732     graphdata.ymin =  -7.9;
     1881    graphdata.ymin =  -6.9;
    17331882    graphdata.ymax = -17.1;
    17341883    KapaSetLimits (myKapa, &graphdata);
    … …  
    17891938
    17901939    graphdata.color = KapaColorByName ("black");
    1791     graphdata.xmax =  -7.9;
     1940    graphdata.xmax =  -6.9;
    17921941    graphdata.xmin = -17.1;
    17931942    graphdata.ymin = -20.1;
Note: See TracChangeset for help on using the changeset viewer.