IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25390


Ignore:
Timestamp:
Sep 15, 2009, 3:46:31 PM (17 years ago)
Author:
eugene
Message:

lots of updates to use moments and psf-moments(sum) to classify objects

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psphot/src/psphotSourceSize.c

    r25359 r25390  
    22# include <gsl/gsl_sf_gamma.h>
    33
    4 bool psphotSourceSizePSF (float *ApResid, float *ApSysErr, psArray *sources, pmPSF *psf, psImageMaskType maskVal);
    5 
    6 bool psphotSourceSizeExtended (FILE *output, pmSource *source, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr);
    7 
    8 bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSyserr);
    9 bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr);
    10 
    11 float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
    12                          psImageMaskType maskVal, const pmModel *model, float Ro);
    13 
    14 bool psphotMaskCosmicRay_Old (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
    15 bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
     4typedef struct {
     5    psImageMaskType maskVal;
     6    psImageMaskType crMask;
     7    float ApResid;
     8    float ApSysErr;
     9    float nSigmaApResid;
     10    float nSigmaMoments;
     11    float nSigmaCR;
     12    float soft;
     13    int grow;
     14} psphotSourceSizeOptions;
     15
     16bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf);
     17bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
     18bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
     19bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options);
     20bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
     21bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
    1622
    1723// we need to call this function after sources have been fitted to the PSF model and
    … …  
    2127// deviation from the psf model at the r = FWHM/2 position
    2228
     29// XXX use an internal flag to mark sources which have already been measured
    2330bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first)
    2431{
    2532    bool status;
     33    psphotSourceSizeOptions options;
    2634
    2735    psTimerStart ("psphot.size");
    2836
    2937    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    30     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    31     assert (maskVal);
     38    options.maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     39    assert (options.maskVal);
    3240
    3341    // bit to mask the cosmic-ray pixels
    34     psImageMaskType crMask  = pmConfigMaskGet("CR", config); // Mask value for cosmic rays
    35 
    36     float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT");
     42    options.crMask  = pmConfigMaskGet("CR", config); // Mask value for cosmic rays
     43
     44    options.nSigmaCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT");
    3745    assert (status);
    3846
    39     // float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT");
    40     // assert (status);
    41 
    42     int grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs
    43     if (!status || grow < 0) {
     47    // XXX recipe name is not great
     48    options.nSigmaApResid = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT");
     49    assert (status);
     50
     51    // XXX recipe name is not great
     52    options.nSigmaMoments = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.MOMENTS");
     53    assert (status);
     54
     55    options.grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs
     56    if (!status || options.grow < 0) {
    4457        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CR.GROW is not positive.");
    4558        return false;
    4659    }
    4760
    48     float soft = psMetadataLookupF32(&status, recipe, "PSPHOT.CR.NSIGMA.SOFTEN"); // Softening parameter
    49     if (!status || !isfinite(soft) || soft < 0.0) {
     61    options.soft = psMetadataLookupF32(&status, recipe, "PSPHOT.CR.NSIGMA.SOFTEN"); // Softening parameter
     62    if (!status || !isfinite(options.soft) || options.soft < 0.0) {
    5063        psWarning("PSPHOT.CR.NSIGMA.SOFTEN not set; defaulting to zero.");
    51         soft = 0.0;
     64        options.soft = 0.0;
    5265    }
    5366
    … …  
    5770    // XXX move this to the code that generates the PSF?
    5871    // XXX store the results on pmPSF?
    59     float ApResid, ApSysErr;
    60     psphotSourceSizePSF (&ApResid, &ApSysErr, sources, psf, maskVal);
     72    psphotSourceSizePSF (&options, sources, psf);
    6173
    6274    // classify the sources based on ApResid and Moments (extended sources)
    63     psphotSourceClass(readout, sources, recipe, psf, maskVal, ApResid, ApSysErr);
    64 
    65     // classify the sources based on the CR test (place this in a function?)
    66     for (int i = first; i < sources->n; i++) {
     75    psphotSourceClass(readout, sources, recipe, psf, &options);
     76
     77    psphotSourceSizeCR (readout, sources, &options);
     78
     79    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot.size"));
     80
     81    psphotVisualPlotSourceSize (recipe, sources);
     82    psphotVisualShowSourceSize (readout, sources);
     83
     84    return true;
     85}
     86
     87bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
     88
     89    // replace the source flux
     90    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     91    source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     92
     93    // flag this as a CR
     94    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     95    pmPeak *peak = source->peak;
     96    psAssert (peak, "NULL peak");
     97
     98    // grab the matching footprint
     99    pmFootprint *footprint = peak->footprint;
     100    if (!footprint) {
     101        // if we have not footprint, use the old code to mask by isophot
     102        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     103        return true;
     104    }
     105
     106    if (!footprint->spans) {
     107        // if we have no footprint, use the old code to mask by isophot
     108        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     109        return true;
     110    }
     111
     112    // mask all of the pixels covered by the spans of the footprint
     113    for (int j = 1; j < footprint->spans->n; j++) {
     114        pmSpan *span1 = footprint->spans->data[j];
     115
     116        int iy = span1->y;
     117        int xs = span1->x0;
     118        int xe = span1->x1;
     119
     120        for (int ix = xs; ix < xe; ix++) {
     121            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     122        }
     123    }
     124    return true;
     125}
     126
     127bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
     128
     129    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     130    pmPeak *peak = source->peak;
     131    psAssert (peak, "NULL peak");
     132
     133    psImage *mask   = source->maskView;
     134    psImage *pixels = source->pixels;
     135    psImage *variance = source->variance;
     136
     137    // XXX This should be a recipe variable
     138# define SN_LIMIT 5.0
     139
     140    int xo = peak->x - pixels->col0;
     141    int yo = peak->y - pixels->row0;
     142
     143    // mark the pixels in this row to the left, then the right
     144    for (int ix = xo; ix >= 0; ix--) {
     145        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
     146        if (SN > SN_LIMIT) {
     147            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
     148        }
     149    }
     150    for (int ix = xo + 1; ix < pixels->numCols; ix++) {
     151        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
     152        if (SN > SN_LIMIT) {
     153            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
     154        }
     155    }
     156
     157    // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
     158    // first go up:
     159    for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
     160        // mark the pixels in this row to the left, then the right
     161        for (int ix = 0; ix < pixels->numCols; ix++) {
     162            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     163            if (SN < SN_LIMIT) continue;
     164
     165            bool valid = false;
     166            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
     167            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
     168            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
     169
     170            if (!valid) continue;
     171            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     172        }
     173    }
     174    // next go down:
     175    for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
     176        // mark the pixels in this row to the left, then the right
     177        for (int ix = 0; ix < pixels->numCols; ix++) {
     178            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     179            if (SN < SN_LIMIT) continue;
     180
     181            bool valid = false;
     182            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
     183            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
     184            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
     185
     186            if (!valid) continue;
     187            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     188        }
     189    }
     190    return true;
     191}
     192
     193// model the apmifit distribution for the psf stars:
     194bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) {
     195
     196    // select the psf stars
     197    psArray *psfstars = psArrayAllocEmpty (100);
     198
     199    psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32);
     200    psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
     201   
     202    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
     203
     204    for (int i = 0; i < sources->n; i++) {
    67205        pmSource *source = sources->data[i];
     206        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
     207        psArrayAdd (psfstars, 100, source);
     208
     209        // XXX can we test if psfMag is set and calculate only if needed?
     210        pmSourceMagnitudes (source, psf, photMode, options->maskVal);
     211       
     212        float apMag = -2.5*log10(source->moments->Sum);
     213        float dMag = source->psfMag - apMag;
     214       
     215        psVectorAppend (Ap, 100, dMag);
     216        psVectorAppend (ApErr, 100, source->errMag);
     217    }
     218
     219    // model the distribution as a mean or median value and a systematic error from that value:
     220    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     221    psVectorStats (stats, Ap, NULL, NULL, 0);
     222
     223    psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32);
     224    for (int i = 0; i < Ap->n; i++) {
     225        dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian;
     226    }
     227
     228    options->ApResid = stats->robustMedian;
     229    options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05);
     230    fprintf (stderr, "psf-sum: %f +/- %f\n", options->ApResid, options->ApSysErr);
     231
     232    return true;
     233}
     234
     235// classify sources based on the combination of psf-mag, Mxx, Myy
     236bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
     237
     238    bool status;
     239    pmPSFClump psfClump;
     240    char regionName[64];
     241
     242    int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     243    for (int i = 0; i < nRegions; i ++) {
     244        snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
     245        psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     246        psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion");
     247
     248        psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION");
     249        psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion");
     250
     251        // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
     252        psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   psAssert (status, "missing PSF.CLUMP.X");
     253        psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   psAssert (status, "missing PSF.CLUMP.Y");
     254        psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  psAssert (status, "missing PSF.CLUMP.DX");
     255        psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  psAssert (status, "missing PSF.CLUMP.DY");
     256
     257        if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) {
     258            psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
     259            continue;
     260        }
     261       
     262        if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, options)) {
     263            psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
     264            continue;
     265        }
     266    }   
     267
     268    return true;
     269}
     270
     271bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
     272
     273    PS_ASSERT_PTR_NON_NULL(sources, false);
     274    PS_ASSERT_PTR_NON_NULL(recipe, false);
     275
     276    int Nsat  = 0;
     277    int Next  = 0;
     278    int Npsf  = 0;
     279    int Ncr   = 0;
     280    int Nmiss = 0;
     281    int Nskip = 0;
     282
     283    pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
     284    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
     285
     286    for (psS32 i = 0 ; i < sources->n ; i++) {
     287
     288        pmSource *source = (pmSource *) sources->data[i];
     289
     290        // psfClumps are found for image subregions:
     291        // skip sources not in this region
     292        if (source->peak->x <  region->x0) continue;
     293        if (source->peak->x >= region->x1) continue;
     294        if (source->peak->y <  region->y0) continue;
     295        if (source->peak->y >= region->y1) continue;
    68296
    69297        // skip source if it was already measured
    70         if (isfinite(source->crNsigma)) {
    71             psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since already measured\n");
     298        if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
     299            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
    72300            continue;
    73301        }
    … …  
    76304        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
    77305            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    78             psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since source is not subtracted\n");
    79             continue;
    80         }
    81 
    82         psF32 **resid  = source->pixels->data.F32;
    83         psF32 **variance = source->variance->data.F32;
    84         psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    85 
    86         // Integer position of peak
    87         int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
    88         int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
    89 
    90         // XXX for now, skip sources which are too close to a boundary
    91         // XXX raise a flag?
    92         if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
    93             yPeak < 1 || yPeak > source->pixels->numRows - 2) {
    94             source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    95             psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
    96             continue;
    97         }
    98 
    99         // XXX for now, just skip any sources with masked pixels
    100         // XXX raise a flag?
    101         bool keep = true;
    102         for (int iy = -1; (iy <= +1) && keep; iy++) {
    103             for (int ix = -1; (ix <= +1) && keep; ix++) {
    104                 if (mask[yPeak+iy][xPeak+ix] & maskVal) {
    105                     keep = false;
    106                 }
    107             }
    108         }
    109         if (!keep) {
    110             psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
    111             source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    112             continue;
    113         }
    114 
    115         // Compare the central pixel with those on either side, for the four possible lines through it.
    116 
    117         // Soften variances (add systematic error)
    118         float softening = soft * PS_SQR(source->peak->flux); // Softening for variances
    119 
    120         // Across the middle: y = 0
    121         float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
    122         float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
    123         float nX = cX / sqrtf(dcX + softening);
    124 
    125         // Up the centre: x = 0
    126         float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
    127         float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
    128         float nY = cY / sqrtf(dcY + softening);
    129 
    130         // Diagonal: x = y
    131         float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
    132         float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
    133         float nL = cL / sqrtf(dcL + softening);
    134 
    135         // Diagonal: x = - y
    136         float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
    137         float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
    138         float nR = cR / sqrtf(dcR + softening);
    139 
    140         // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
    141         // Ndof = 4 ? (four measurements, no free parameters)
    142         // XXX this value is going to be biased low because of systematic errors.
    143         // we need to calibrate it somehow
    144         // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
    145 
    146         // not strictly accurate: overcounts the chisq contribution from the center pixel (by
    147         // factor of 4); also biases a bit low if any pixels are masked
    148         // XXX I am not sure I want to keep this value...
    149         source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR);
    150 
    151         float fCR = 0.0;
    152         int nCR = 0;
    153         if (nX > 0.0) {
    154             fCR += nX;
    155             nCR ++;
    156         }
    157         if (nY > 0.0) {
    158             fCR += nY;
    159             nCR ++;
    160         }
    161         if (nL > 0.0) {
    162             fCR += nL;
    163             nCR ++;
    164         }
    165         if (nR > 0.0) {
    166             fCR += nR;
    167             nCR ++;
    168         }
    169         source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
    170         if (!isfinite(source->crNsigma)) {
    171             continue;
    172         }
    173 
    174         // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    175         if (source->crNsigma > CR_NSIGMA_LIMIT) {
    176             // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask);
    177             // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask);
    178         }
    179     }
    180 
    181     // pause and wait for user input:
    182     // continue, save (provide name), ??
    183     char key[10];
    184     fprintf (stdout, "[c]ontinue? ");
    185     if (!fgets(key, 8, stdin)) {
    186         psWarning("Unable to read option");
    187     }
    188  
    189     // now that we have masked pixels associated with CRs, we can grow the mask
    190     if (grow > 0) {
    191         bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask
    192         psImage *newMask = psImageConvolveMask(NULL, readout->mask, crMask, crMask, -grow, grow, -grow, grow);
    193         psImageConvolveSetThreads(oldThreads);
    194         if (!newMask) {
    195             psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
    196             return false;
    197         }
    198         psFree(readout->mask);
    199         readout->mask = newMask;
    200     }
    201 
    202     psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n",
    203               sources->n - first, psTimerMark ("psphot.size"));
    204 
    205     psphotVisualPlotSourceSize (sources);
    206     psphotVisualShowSourceSize (readout, sources);
     306            psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
     307            continue;
     308        }
     309
     310        // we are basically classifying by moments; use the default if not found
     311        psAssert (source->moments, "why is this source missing moments?");
     312        if (source->mode & noMoments) {
     313            Nskip ++;
     314            continue;
     315        }
     316
     317        psF32 Mxx = source->moments->Mxx;
     318        psF32 Myy = source->moments->Myy;
     319
     320        // XXX can we test if psfMag is set and calculate only if needed?
     321        pmSourceMagnitudes (source, psf, photMode, options->maskVal);
     322
     323        float apMag = -2.5*log10(source->moments->Sum);
     324        float dMag = source->psfMag - apMag;
     325        float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
     326
     327        source->extNsigma = nSigma;
     328        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     329
     330        // Anything within this region is a probably PSF-like object. Saturated stars may land
     331        // in this region, but are detected elsewhere on the basis of their peak value.
     332        bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY);
     333        if (isPSF) {
     334            Npsf ++;
     335            continue;
     336        }
     337
     338        // Defects may not always match CRs from peak curvature analysis
     339        // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
     340        if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
     341            source->mode |= PM_SOURCE_MODE_DEFECT;
     342            Ncr ++;
     343            continue;
     344        }
     345
     346        // saturated star (determined in PSF fit).  These may also be saturated galaxies, or
     347        // just large saturated regions.
     348        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     349            Nsat ++;
     350            continue;
     351        }
     352
     353        // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
     354        bool isEXT = (nSigma > options->nSigmaApResid) && (Mxx > psfClump->X) && (Myy > psfClump->Y);
     355        if (isEXT) {
     356            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
     357            Next ++;
     358            continue;
     359        }
     360
     361        Nmiss ++;
     362    }
     363
     364    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %d %d %d %d %d %d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip);
    207365
    208366    return true;
    … …  
    211369// given the PSF ellipse parameters, navigate around the 1sigma contour, return the total
    212370// deviation in sigmas.  This is measured on the residual image - should we ignore negative
    213 // deviations?
     371// deviations?  NOTE: This function was an early attempt to classify extended objects, and is
     372// no longer used by psphot.
    214373float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
    215374                         psImageMaskType maskVal, const pmModel *model, float Ro)
    … …  
    282441}
    283442
    284 bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
    285 
    286     // replace the source flux
    287     pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    288     source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    289 
    290     // flag this as a CR
    291     source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    292     pmPeak *peak = source->peak;
    293     psAssert (peak, "NULL peak");
    294 
    295     // grab the matching footprint
    296     pmFootprint *footprint = peak->footprint;
    297     if (!footprint) {
    298         // if we have not footprint, use the old code to mask by isophot
    299         psphotMaskCosmicRay_Old (source, maskVal, crMask);
    300         return true;
    301     }
    302 
    303     if (!footprint->spans) {
    304         // if we have not footprint, use the old code to mask by isophot
    305         psphotMaskCosmicRay_Old (source, maskVal, crMask);
    306         return true;
    307     }
    308 
    309     // mask all of the pixels covered by the spans of the footprint
    310     for (int j = 1; j < footprint->spans->n; j++) {
    311         pmSpan *span1 = footprint->spans->data[j];
    312 
    313         int iy = span1->y;
    314         int xs = span1->x0;
    315         int xe = span1->x1;
    316 
    317         for (int ix = xs; ix < xe; ix++) {
    318             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    319         }
    320     }
    321     return true;
    322 }
    323 
    324 bool psphotMaskCosmicRay_Old (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
    325 
    326     source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    327     pmPeak *peak = source->peak;
    328     psAssert (peak, "NULL peak");
    329 
    330     psImage *mask   = source->maskView;
    331     psImage *pixels = source->pixels;
    332     psImage *variance = source->variance;
    333 
    334     // XXX This should be a recipe variable
    335 # define SN_LIMIT 5.0
    336 
    337     int xo = peak->x - pixels->col0;
    338     int yo = peak->y - pixels->row0;
    339 
    340     // mark the pixels in this row to the left, then the right
    341     for (int ix = xo; ix >= 0; ix--) {
    342         float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
    343         if (SN > SN_LIMIT) {
    344             mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
    345         }
    346     }
    347     for (int ix = xo + 1; ix < pixels->numCols; ix++) {
    348         float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
    349         if (SN > SN_LIMIT) {
    350             mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
    351         }
    352     }
    353 
    354     // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
    355     // first go up:
    356     for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
    357         // mark the pixels in this row to the left, then the right
    358         for (int ix = 0; ix < pixels->numCols; ix++) {
    359             float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    360             if (SN < SN_LIMIT) continue;
    361 
    362             bool valid = false;
    363             valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
    364             valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
    365             valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
    366 
    367             if (!valid) continue;
    368             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    369         }
    370     }
    371     // next go down:
    372     for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
    373         // mark the pixels in this row to the left, then the right
    374         for (int ix = 0; ix < pixels->numCols; ix++) {
    375             float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    376             if (SN < SN_LIMIT) continue;
    377 
    378             bool valid = false;
    379             valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
    380             valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
    381             valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
    382 
    383             if (!valid) continue;
    384             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    385         }
    386     }
    387     return true;
    388 }
    389 
    390 // is this source extended or not?
    391 bool psphotSourceSizeExtended (FILE *output, pmSource *source, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr) {
    392 
    393     // XXX can we test if psfMag is set and calculate only if needed?
    394     pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
    395     pmSourceMagnitudes (source, psf, photMode, maskVal);
    396 
    397     float apMag = -2.5*log10(source->moments->Sum);
    398     float dMag = source->psfMag - apMag;
    399    
    400     float nSigma = (dMag - ApResid) / hypot(source->errMag, ApSysErr);
    401 
    402     fprintf (output, "%f %f : %f %f : %f %f : %f %f\n", source->peak->xf, source->peak->yf, nSigma, dMag, source->psfMag, source->errMag, source->moments->Mxx, source->moments->Myy);
    403 
    404     return true;
    405 }
    406 
    407 // model the apmifit distribution for the psf stars:
    408 bool psphotSourceSizePSF (float *ApResid, float *ApSysErr, psArray *sources, pmPSF *psf, psImageMaskType maskVal) {
    409 
    410     // select the psf stars
    411     psArray *psfstars = psArrayAllocEmpty (100);
    412 
    413     psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32);
    414     psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
    415    
    416     pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
    417 
     443bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
     444
     445    // classify the sources based on the CR test (place this in a function?)
     446    // XXX use an internal flag to mark sources which have already been measured
    418447    for (int i = 0; i < sources->n; i++) {
    419448        pmSource *source = sources->data[i];
    420         if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
    421         psArrayAdd (psfstars, 100, source);
    422 
    423         // XXX can we test if psfMag is set and calculate only if needed?
    424         pmSourceMagnitudes (source, psf, photMode, maskVal);
    425        
    426         float apMag = -2.5*log10(source->moments->Sum);
    427         float dMag = source->psfMag - apMag;
    428        
    429         psVectorAppend (Ap, 100, dMag);
    430         psVectorAppend (ApErr, 100, source->errMag);
    431     }
    432 
    433     // model the distribution as a mean or median value and a systematic error from that value:
    434     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    435     psVectorStats (stats, Ap, NULL, NULL, 0);
    436 
    437     psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32);
    438     for (int i = 0; i < Ap->n; i++) {
    439         dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian;
    440     }
    441 
    442     *ApResid = stats->robustMedian;
    443     *ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05);
    444     fprintf (stderr, "psf-sum: %f +/- %f\n", *ApResid, *ApSysErr);
    445 
    446     return true;
    447 }
    448 
    449 // classify sources based on the combination of psf-mag, Mxx, Myy
    450 bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr) {
    451 
    452     bool status;
    453     pmPSFClump psfClump;
    454     char regionName[64];
    455 
    456     int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
    457     for (int i = 0; i < nRegions; i ++) {
    458         snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
    459         psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
    460         psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion");
    461 
    462         psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION");
    463         psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion");
    464 
    465         // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
    466         psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   psAssert (status, "missing PSF.CLUMP.X");
    467         psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   psAssert (status, "missing PSF.CLUMP.Y");
    468         psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  psAssert (status, "missing PSF.CLUMP.DX");
    469         psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  psAssert (status, "missing PSF.CLUMP.DY");
    470 
    471         if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) {
    472             psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
    473             continue;
    474         }
    475        
    476         if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, maskVal, ApResid, ApSysErr)) {
    477             psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
    478             continue;
    479         }
    480     }   
    481 
    482     return true;
    483 }
    484 
    485 bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr) {
    486 
    487     PS_ASSERT_PTR_NON_NULL(sources, false);
    488     PS_ASSERT_PTR_NON_NULL(recipe, false);
    489 
    490     int Nsat  = 0;
    491     int Next  = 0;
    492     int Npsf  = 0;
    493     int Ncr   = 0;
    494     int Nmiss = 0;
    495     int Nskip = 0;
    496 
    497     pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
    498     pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
    499 
    500     for (psS32 i = 0 ; i < sources->n ; i++) {
    501 
    502         pmSource *source = (pmSource *) sources->data[i];
    503 
    504         // psfClumps are found for image subregions:
    505         // skip sources not in this region
    506         if (source->peak->x <  region->x0) continue;
    507         if (source->peak->x >= region->x1) continue;
    508         if (source->peak->y <  region->y0) continue;
    509         if (source->peak->y >= region->y1) continue;
    510 
    511         // we are basically classifying by moments; use the default if not found
    512         psAssert (source->moments, "why is this source missing moments?");
    513         if (source->mode & noMoments) {
    514             Nskip ++;
    515             continue;
    516         }
    517 
    518         psF32 Mxx = source->moments->Mxx;
    519         psF32 Myy = source->moments->Myy;
    520 
    521         // saturated star (determined in PSF fit)
    522         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    523             Nsat ++;
    524             continue;
    525         }
    526 
    527         // XXX can we test if psfMag is set and calculate only if needed?
    528         pmSourceMagnitudes (source, psf, photMode, maskVal);
    529 
    530         float apMag = -2.5*log10(source->moments->Sum);
    531         float dMag = source->psfMag - apMag;
    532         float nSigma = (dMag - ApResid) / hypot(source->errMag, ApSysErr);
    533 
    534         source->extNsigma = nSigma;
    535 
    536         // XXX these sigma cuts should be user-configured parameters
    537         bool isPSF = (fabs(nSigma) < 3.0) && (fabs(Mxx - psfClump->X) < 2.0*psfClump->dX) && (fabs(Myy - psfClump->Y) < 2.0*psfClump->dY);
    538         if (isPSF) {
    539             Npsf ++;
    540             continue;
    541         }
    542 
    543         // XXX get the sign on nSigma right here and above
    544         bool isEXT = (nSigma < 3.0) && (Mxx > psfClump->X) && (Myy > psfClump->Y);
    545         if (isEXT) {
    546             source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    547             Next ++;
    548             continue;
    549         }
    550 
    551        
    552         // XXX possible or likely?  need to compare my standard CR test with this
    553         // Mark in both cases?
    554         if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
    555             Ncr ++;
    556             continue;
    557         }
    558 
    559         Nmiss ++;
    560     }
    561 
    562     psLogMsg("psModules.objects", PS_LOG_INFO, "Rough classifications: %d %d %d %d %d %d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip);
    563 
    564     return true;
    565 }
     449
     450        // skip source if it was already measured
     451        if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
     452            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
     453            continue;
     454        }
     455
     456        // source must have been subtracted
     457        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
     458            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
     459            psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
     460            continue;
     461        }
     462
     463        psF32 **resid  = source->pixels->data.F32;
     464        psF32 **variance = source->variance->data.F32;
     465        psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     466
     467        // Integer position of peak
     468        int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
     469        int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
     470
     471        // Skip sources which are too close to a boundary.  These are mostly caught as DEFECT
     472        if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
     473            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
     474            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
     475            continue;
     476        }
     477
     478        // Skip sources with masked pixels.  These are mostly caught as DEFECT
     479        bool keep = true;
     480        for (int iy = -1; (iy <= +1) && keep; iy++) {
     481            for (int ix = -1; (ix <= +1) && keep; ix++) {
     482                if (mask[yPeak+iy][xPeak+ix] & options->maskVal) {
     483                    keep = false;
     484                }
     485            }
     486        }
     487        if (!keep) {
     488            psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
     489            continue;
     490        }
     491
     492        // Compare the central pixel with those on either side, for the four possible lines through it.
     493
     494        // Soften variances (add systematic error)
     495        float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances
     496
     497        // Across the middle: y = 0
     498        float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
     499        float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
     500        float nX = cX / sqrtf(dcX + softening);
     501
     502        // Up the centre: x = 0
     503        float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
     504        float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
     505        float nY = cY / sqrtf(dcY + softening);
     506
     507        // Diagonal: x = y
     508        float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
     509        float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
     510        float nL = cL / sqrtf(dcL + softening);
     511
     512        // Diagonal: x = - y
     513        float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
     514        float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
     515        float nR = cR / sqrtf(dcR + softening);
     516
     517        // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
     518        // Ndof = 4 ? (four measurements, no free parameters)
     519        // XXX this value is going to be biased low because of systematic errors.
     520        // we need to calibrate it somehow
     521        // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
     522
     523        // not strictly accurate: overcounts the chisq contribution from the center pixel (by
     524        // factor of 4); also biases a bit low if any pixels are masked
     525        // XXX I am not sure I want to keep this value...
     526        source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR);
     527
     528        float fCR = 0.0;
     529        int nCR = 0;
     530        if (nX > 0.0) {
     531            fCR += nX;
     532            nCR ++;
     533        }
     534        if (nY > 0.0) {
     535            fCR += nY;
     536            nCR ++;
     537        }
     538        if (nL > 0.0) {
     539            fCR += nL;
     540            nCR ++;
     541        }
     542        if (nR > 0.0) {
     543            fCR += nR;
     544            nCR ++;
     545        }
     546        source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
     547        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     548
     549        if (!isfinite(source->crNsigma)) {
     550            continue;
     551        }
     552
     553        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
     554        if (source->crNsigma > options->nSigmaCR) {
     555            source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     556            // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
     557            // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask);
     558        }
     559    }
     560
     561    // now that we have masked pixels associated with CRs, we can grow the mask
     562    if (options->grow > 0) {
     563        bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask
     564        psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow);
     565        psImageConvolveSetThreads(oldThreads);
     566        if (!newMask) {
     567            psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
     568            return false;
     569        }
     570        psFree(readout->mask);
     571        readout->mask = newMask;
     572    }
     573    return true;
     574}
Note: See TracChangeset for help on using the changeset viewer.