IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27532 for trunk/psphot/src


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/src
Files:
10 edited
2 copied

Legend:

Unmodified
Added
Removed

  • 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.