IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 4115


Ignore:
Timestamp:
Jun 4, 2005, 5:22:13 PM (21 years ago)
Author:
eugene
Message:

running sources in brightness order, using a mask

Location:
trunk/psphot
Files:
4 added
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/doc/psphot.txt

    r4114 r4115  
    11
    2 2005.04.12 : Notes on psphot
     2psphot is a stand-alone program which implements the Pan-STARRS IPP
     3object detection and analysis functions for using with individual
     4images. 
    35
    4 - psPeak, psSource, etc: should be pmPeak, pmSource, etc
    5 - pmCullPeaks: should be pmCullImagePeaks
    6 
    7 - does pmSourceMoments need an image argument?
    8 
    9 - psLibInit / p_psTimeInit
    10 
    11   we should not have to know about astronomy time concepts (ie, have a
    12   timeConfig file) in order to use the time functions for basic work.
    13   This means having a PS_TIME_UNIX which just works with the UNIX
    14   clock and avoids all initialization issues.
    15 
    16 - psImageStats is much too slow for basic median
    17   - this function does a complete sort on all pixel values.  we need
    18   to use the histogram method (spec a different type of median)
    19 
    20 - psImageStats, ROBUST_MEDIAN fails
    21 - there are errors in psVectorStats for ROBUST MEAN, STDEV (& median?)
    22   * this is ill-defined.  revise the ADD: too sensitive to the
    23   binning, smoothing scales
    24 
    25 - pmFindImagePeaks is not finding any peaks:
    26   typo in interior row section: see bug 359
    27   after this fix: PASS
    28 
    29 * p_psGetRowVectorFromImage : PASS
    30 * psGetRowVectorFromImage : PASS
    31 - both have about the same time (~15 us +/- 5us) for 2048 pix on alala
    32 
    33 * pmFindVectorPeaks : PASS (returned row data in data.U32)
    34 
    35 - # include <string.h> : needed in psMetadata.c
    36 
    37 * MyListAddPeak / pmPeakAlloc are inconsistent wrt col,row vs x,y
    38 
    39 - why does psMetadataItemAllocV allocate a string of length
    40   MAX_STRING_LENGTH for the comment, then use strncpy?
    41 
    42   should be:
    43 
    44     // set metadata item comment
    45     if (comment == NULL) {
    46         // Per SDRS, null isn't allowed, must use "" instead
    47         metadataItem->comment = psStringCopy ("");
    48     } else {
    49         metadataItem->comment = psStringCopy (comment);
    50     }
    51 
    52 psMetadataAdd -> psMetadataAddItem is all confused:
    53  - we end up with two psMetadataItemAlloc calls, when only one should
    54    be used
    55  - psMetadataItemAlloc was doing the wrong thing if the incoming data
    56    was NULL
    57  - I think the logic in this block is wrong as well. 
    58 
    59  * I fixed this in psMetadata.c
    60 
    61 2005.04.14 :
    62 
    63  - psImageStats is still much too slow
    64 
    65  - pmObjects : changed to using psArray to carry peaks / sources,
    66    rather than psList
    67 
    68  - added function pmPeaksSubset to replace pmCullPeaks
    69 
    70  - pmSourceLocalSky: added min/max fncs to force subimage to stay in
    71    image
    72 
    73  - Sx, Sy can go negative: have forced limit to 0,0
    74 
    75  - pmSourceRoughClass: implemented sigmaX, sigmaY search for stellar
    76    clump (uses pmComparePeakDescend)
    77 
    78    - the search ignores the 0,0 pixel
    79 
    80    - pmSourceRoughClass uses all detected peaks.  allow an exclusion
    81      for peak > max, peak < min?
    82 
    83    - need to push clump stats on the metadata
    84 
    85    - the source classes are assigned so they are mutually exclusive
    86 
    87    - change PS_SOURCE_OTHER to STAR?
    88 
    89  - pmSourceSetPixelCircle:
    90 
    91    - name should be pmSourceSetPixelsCircle (geometry could be an
    92      option...)
    93 
    94    - I cleaned the defined region to match the convensions of
    95      pmSourceLocalSky
    96 
    97    - truncate and force subimage to lie on image
    98 
    99    - CheckRadius2 seems like an inefficient coding
    100 
    101    - valid pixels have mask value of 1 (inverse of other funcs)
    102 
    103  - pmSourceFitModel:
    104 
    105    - num iterations is much too large
    106 
    107    -
    108 
    109 mrq:
    110 
    111  - allocate alpha, beta arrays
    112  - mrq2dcof (pars, alpha, beta)
    113  - make guess pars', alpha', beta'
    114  - mrq2dcof (pars', alpha', beta')
    115  - if (chisq < old chisq), keep new guess
    116 
    117 
    118 2005.04.19:
    119 
    120 - psMinimizeLMChi2 has some serious problems.  I re-wrote it starting
    121   from the mrqmin example in ohana.   after a variety little bugs, it
    122   seems to work quite well.  the program 'modeltest.c' runs two
    123   different gaussians through psMinimizeLMChi2, and gets the right
    124   answer quite quickly.  The first pass at the implementation had a
    125   variety of problems.  to avoid any worry about errors in the
    126   psMatrix code, I implemented psGaussJordan, basing the code on the
    127   elixir gaussj.c code.  this worked fine, and let me find the errors
    128   elsewhere in the code. the LUD version now works just the same as
    129   the GaussJ version.
    130 
    131 
    132 2005.04.21
    133 
    134 - psMinimizeLMChi2 does an excellent job now in most cases.  I am able
    135   to run pmSourceFitModel on the PSF stars quite well.  I added in the
    136   (nearly) correct errors (actually, just using sqrt(N), limiting to
    137   N>=1).  I see one problem object, which is quite faint, and does not
    138   converge well: I get very large values for x,y, Sx, Sy, etc.  I need
    139   to think about constraining these fits from running way off scale.
    140 
    141   fit:   0.123879 sec for 277 stars (gauss)
    142   fit:   0.155940 sec for 277 stars (pgauss) (!)
    143   fit:   0.153292 sec for 277 stars
    144   (0.5 msec per star)

  • trunk/psphot/src/apply_psf_model.c

    r4114 r4115  
    2222    keep           = psRegionForImage (keep, image, keep);
    2323
    24     // set the object surfact-brightness limit for fitted pixels
     24    float shapeNsigma = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
     25    float snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
     26
     27    // set the object surface-brightness limit for fitted pixels
    2528    float FLUX_LIMIT  = FIT_NSIGMA * sky->sampleStdev;
    2629    psLogMsg ("psphot.apply_psf_model", 3, "fitting pixels with at least %f object counts\n", FLUX_LIMIT);
    … …  
    3134    for (int i = 0; i < sources->n; i++) {
    3235        psSource *source = sources->data[i];
     36
     37        // skip the existing PSF stars
     38        // XXX drop this -- refit?
    3339        if (source->modelPSF != NULL) continue;
    3440
    … …  
    5561        pmSourceMaskSaturated (source, SATURATE);
    5662
    57         // fit as PSF, not FLT (drop poor fits)
     63        // fit as PSF, not FLT (skip poor fits)
    5864        if (!pmSourceFitModel (source, model, true)) continue;
    5965        source->modelPSF = model;
    6066        Niter += model[0].nIter;
    6167        Nfit ++;
     68
     69        mark_psf_source (source);
     70        subtract_psf_source (source);
    6271    }
    6372    psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);

  • trunk/psphot/src/choose_psf_model.c

    r4114 r4115  
    6868    for (int i = 0; i < test->sources->n; i++) {
    6969        psSource *source = test->sources->data[i];
    70         source->modelPSF = test->modelPSF->data[i];
     70        // drop masked sources from PSFSTAR list
     71        if (test->mask->data.U8[i]) {
     72          source->type = PS_SOURCE_OTHER;
     73          // XXX is this the right type to go to?
     74        } else {
     75          source->modelPSF = test->modelPSF->data[i];
     76        }
    7177    }
    7278    return (test->psf);

  • trunk/psphot/src/find_peaks.c

    r4114 r4115  
    99
    1010    // smooth the image
    11     // should we also subtract a super-binned image? (as an option)
     11    // should we also subtract a super-binned image? (as an option?)
    1212
    1313    psTimerStart ("psphot");
    … …  
    2424
    2525    psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot"));
    26 
    2726
    2827    // find the peaks in the smoothed image

  • trunk/psphot/src/image_stats.c

    r4114 r4115  
    1313    // pass on the stats
    1414    {
    15         bool   status    = false;
     15        bool     status  = false;
     16        int      Npixels = image->numRows*image->numCols;
     17        int      Nsubset = PS_MIN (Npixels, psMetadataLookupF32 (&status, config, "NSUBSET"));
     18        psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32);
    1619        psRandom *rnd    = psRandomAlloc (PS_RANDOM_TAUS, 0);
    17         int    Nsubset   = psMetadataLookupF32 (&status, config, "NSUBSET");
    18         double fSubset   = (double) Nsubset / (image->numRows*image->numCols);
    19         psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32);
    2020
    21         subset->n = 0;
    22         for (int i = 0; i < image->numRows; i++) {
    23             for (int j = 0; j < image->numCols; j++) {
    24                 double frnd = psRandomUniform (rnd);
    25                 if ((fSubset < 1.0) && (fSubset < frnd)) continue;
    26                 subset->data.F32[subset->n] = image->data.F32[j][i];
    27                 subset->n ++;
    28                 if (subset->n == Nsubset) goto got_subset;
    29             }
     21        // choose Nsubset points between 0 and Nx*Ny, convert to coords
     22        Npixels = image->numRows*image->numCols;
     23        for (int i = 0; i < Nsubset; i++) {
     24          double frnd = psRandomUniform (rnd);
     25          int pixel = Npixels * frnd;
     26          int ix = pixel / image->numCols;
     27          int iy = pixel % image->numCols;
     28          subset->data.F32[i] = image->data.F32[iy][ix];
    3029        }
    3130
    32     got_subset:
    3331        // this should use ROBUST not SAMPLE median
    3432        // robust median is broken in pslib (0.5)

  • trunk/psphot/src/psPolynomials.c

    r4114 r4115  
    191191// XXX EAM : test version of 1d fitting
    192192psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly,
    193         const psVector* x,
    194         const psVector* y,
    195         const psVector* yErr)
     193                                             psVector *mask,
     194                                             const psVector *x,
     195                                             const psVector *y,
     196                                             const psVector *yErr)
    196197{
    197198    // I think this is 1 dimension down
    … …  
    232233    // Build the B and A data structs.
    233234    for (int k = 0; k < x->n; k++) {
     235        if ((mask != NULL) && mask->U8[k]) continue;
    234236        xSums = BuildSums1D(xSums, x->data.F64[k], nTerm);
    235237     
    … …  
    291293// XXX EAM : test version of 2d fitting
    292294psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly,
    293         const psVector* x,
    294         const psVector* y,
    295         const psVector* z,
    296         const psVector* zErr)
     295                                             psVector* mask,
     296                                             const psVector* x,
     297                                             const psVector* y,
     298                                             const psVector* z,
     299                                             const psVector* zErr)
    297300{
    298301    // I think this is 1 dimension down
    … …  
    322325    // Build the B and A data structs.
    323326    for (int k  = 0; k < x->n; k++) {
     327      if ((mask != NULL) && mask->data.U8[i]) continue;
    324328        Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
    325329     
    … …  
    391395}   
    392396
    393 psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
     397psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly,
    394398                            const psVector* x,
    395399                            const psVector* y,
    … …  
    413417    for (int N = 0; N < 3; N++) {
    414418        // XXX EAM : this would be better defined with an element mask
    415         poly   = VectorFitPolynomial2DOrd_EAM (poly, X, Y, Z, dZ);
     419        poly   = VectorFitPolynomial2DOrd_EAM (poly, NULL, X, Y, Z, dZ);
    416420        zFit   = Polynomial2DEvalVectorD (poly, x, y);
    417421        zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit);
    … …  
    424428        int n = 0;
    425429        for (int i = 0; i < zResid->n; i++) {
    426             if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) continue;
     430            if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
     431              continue;
     432            }
    427433            X->data.F64[n]  =  x->data.F64[i];
    428434            Y->data.F64[n]  =  y->data.F64[i];
    … …  
    435441        Z->n = n;
    436442        dZ->n = n;
     443    }
     444    return (poly);
     445}
     446
     447// XXX EAM : be careful here with F32 vs F64 vectors
     448psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
     449                            psVector* mask,
     450                            const psVector* x,
     451                            const psVector* y,
     452                            const psVector* z,
     453                            const psVector* dz)
     454{
     455    psVector *zFit   = NULL;
     456    psVector *zResid = NULL;
     457    psStats  *stats  = NULL;
     458
     459    for (int N = 0; N < 3; N++) {
     460        poly   = VectorFitPolynomial2DOrd_EAM (poly, mask, x, y, z, dz);
     461        zFit   = Polynomial2DEvalVectorD (poly, x, y);
     462        zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit);
     463
     464        stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     465        stats  = psVectorStats (stats, zResid, NULL, mask, 1);
     466        psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g (%d pts)\n",
     467                 stats->clippedMean, stats->clippedStdev, stats->clippedNvalues);
     468
     469        // set mask if pts are not valid
     470        int n = 0;
     471        for (int i = 0; i < zResid->n; i++) {
     472            if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
     473              mask->data.U8[i] = 1;
     474              continue;
     475            } else {
     476              mask->data.U8[i] = 0;
     477            }         
     478            n++;
     479        }
    437480    }
    438481    return (poly);

  • trunk/psphot/src/psphot.c

    r4114 r4115  
    2424    peaks = find_peaks (image, config, sky);
    2525
     26    // construct sources and measure basic stats
    2627    sources = source_moments (image, config, peaks);
     28
     29    // source analysis is done in S/N order (brightest first)
     30    sources = psArraySort (sources, by_SN);
    2731
    2832    // use stellar objects SN > PSF_SN_LIM
    … …  
    3539    // XXX it has some problems
    3640    // test_psf_scatter (sources);
    37 
    38     // identify PSF-like objects (< PSF_SHAPE_NSIGMA from dSx,dSy model)
    39     // SN <= FAINT_SN_LIM : FAINTSTAR
    40     // SN >  FAINT_SN_LIM : BRIGHTSTAR
    41     mark_psf_sources (sources, config);
    42 
    43     // subtract BRIGHTSTAR, PSFSTAR, SATSTAR
    44     subtract_psf_sources (sources);
    4541
    4642    // fit_galaxies (image, config, sources);

  • trunk/psphot/src/psphot.h

    r4114 r4115  
    1313  psArray    *modelFLT;     // model fits, floating parameters
    1414  psArray    *modelPSF;     // model fits, PSF parameters
     15  psVector   *mask;
    1516  psVector   *metric;
    1617  psStats    *metricStats;

  • trunk/psphot/src/pspsf.c

    r4114 r4115  
    3636    test->metricStats = NULL;
    3737    test->metric      = psVectorAlloc (sources->n, PS_TYPE_F64);
     38    test->mask        = psVectorAlloc (sources->n, PS_TYPE_U8);
    3839
    3940    for (int i = 0; i < test->modelFLT->n; i++) {
    4041        test->modelFLT->data[i] = NULL;
    4142        test->modelPSF->data[i] = NULL;
     43        test->mask->data.U8[i]  = 0;
    4244    }   
    4345    return (test);
    … …  
    5759        psModel  *model  = pmSourceModelGuess (source, test->modelType);
    5860
    59         // fit model as FLT, not PSF (skip poor fits)
    60         if (!pmSourceFitModel (source, model, false)) continue;
     61        // fit model as FLT, not PSF (mask & skip poor fits)
     62        if (!pmSourceFitModel (source, model, false)) {
     63          test->mask->data.U8[i] = 1;
     64          continue;
     65        }
    6166        test->modelFLT->data[i] = model;
    6267        Nflt ++;
    … …  
    6671
    6772    // stage 2: construct a psf (pmPSF) from this collection of model fits
    68     pmPSFFromModels (test->psf, test->modelFLT);
     73    pmPSFFromModels (test->psf, test->modelFLT, test->mask);
    6974
    7075    // stage 3: refit with fixed shape parameters
    … …  
    7378        psSource *source = test->sources->data[i];
    7479        psModel  *modelFLT = test->modelFLT->data[i];
    75         if (modelFLT == NULL) continue;
     80
     81        // masked for: bad model fit, outlier in parameters
     82        if (test->mask->data.U8[i]) continue;
     83        if (modelFLT == NULL) {
     84          psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
     85          continue;
     86        }
    7687        psModel  *modelPSF = psModelFromPSF (modelFLT, test->psf); // set shape for this model
    7788
    7889        // fit model as PSF, not FLT (skip poor fits)
    79         if (!pmSourceFitModel (source, modelPSF, true)) continue;
     90        if (!pmSourceFitModel (source, modelPSF, true)) {
     91          test->mask->data.U8[i] = 1;
     92          continue;
     93        }
    8094        test->modelPSF->data[i] = modelPSF;
    8195        Npsf ++;
    … …  
    92106        float fitSum = 0;
    93107     
     108        // is this metricMask redundant with test->mask ?
    94109        metricMask->data.U8[i] = 1;
    95110        test->metric->data.F64[i] = 0;
    96111
     112        if (test->mask->data.U8[i]) continue;
     113
    97114        psModel *model = test->modelPSF->data[i];
    98         if (model == NULL) continue;
     115        if (model == NULL) {
     116          psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
     117          continue;
     118        }
    99119
    100120        psImage *image = ((psSource *)test->sources->data[i])->pixels;
    101121        psImage *mask  = ((psSource *)test->sources->data[i])->mask;
    102122
     123        // this metric is Ap-Fit
    103124        float sky = model->params->data.F32[0];
    104 
    105125        for (int ix = 0; ix < image->numCols; ix++) {
    106126            for (int iy = 0; iy < image->numRows; iy++) {
    … …  
    110130            }
    111131        }
    112         // this metric is Ap-Fit
    113132        // fprintf (stderr, "%d %f %f  %f\n", i, obsSum, fitSum, test->metric->data.F64[i]);
     133        // keep the good metrics
    114134        if ((fitSum > 0) && (obsSum > 0)) {
    115135            metricMask->data.U8[i] = 0;
    116136            test->metric->data.F64[i] = -2.5*log10(obsSum/fitSum);
    117         }       
     137        } else {
     138          test->mask->data.U8[i] = 1;
     139        }
    118140    }
    119141
    … …  
    128150// input: an array of psModels, pre-allocated psf
    129151// some of the array entries may be NULL, ignore them
    130 bool pmPSFFromModels (pmPSF *psf, psArray *models) {
     152bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) {
    131153
    132154    int n;
    133155
    134     // construct the x and y vectors from the collection of objects
    135     // count the number of valid model measurements
     156    // construct the fit vectors from the collection of objects
     157    // use the mask to ignore missing fits
    136158    psVector *x = psVectorAlloc (models->n, PS_TYPE_F64);
    137159    psVector *y = psVectorAlloc (models->n, PS_TYPE_F64);
    138     n = 0;
     160    psVector *z = psVectorAlloc (models->n, PS_TYPE_F64);
     161    psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
     162
    139163    for (int i = 0; i < models->n; i++) {
    140164        psModel *model = models->data[i];
    … …  
    142166
    143167        // XXX EAM : this is fragile: x and y must be stored consistently at 2,3
    144         x->data.F64[n] = model->params->data.F32[2];
    145         y->data.F64[n] = model->params->data.F32[3];
    146         n++;
    147     }
    148     x->n = y->n = n;
     168        x->data.F64[i] = model->params->data.F32[2];
     169        y->data.F64[i] = model->params->data.F32[3];
     170    }
    149171
    150172    // we are doing a robust fit.  after each pass, we drop points which are
    151     // more deviant than three sigma.
    152     psVector *z = psVectorAlloc (x->n, PS_TYPE_F64);
    153     psVector *dz = psVectorAlloc (x->n, PS_TYPE_F64);
     173    // more deviant than three sigma. 
     174    // mask is currently updated for each pass. this is inconsistent
     175
    154176    for (int i = 0; i < psf->params->n; i++) {
    155         n = 0;
    156177        for (int j = 0; j < models->n; j++) {
    157178            psModel *model = models->data[j];
    158179            if (model == NULL) continue;
    159             z->data.F64[n] = model->params->data.F32[i + 4];
    160             dz->data.F64[n] = 1;
    161             n++;
     180            z->data.F64[j] = model->params->data.F32[i + 4];
     181            dz->data.F64[j] = 1;
    162182            // XXX EAM : need to use actual errors?
    163183            // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
    164184        }
    165         if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n");
    166         z->n = n;
    167         dz->n = n;
    168 
    169         psf->params->data[i] = RobustFit2D (psf->params->data[i], x, y, z, dz);
     185        // if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n");
     186
     187        psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
    170188        // psPolynomial2DDump (psf->params->data[i]);
    171189    }

  • trunk/psphot/src/subtract_psf_sources.c

    r4114 r4115  
    66bool subtract_psf_sources (psArray *sources)
    77{
     8
     9    sources = psArraySort (sources, by_SN);
     10
    811    for (int i = 0; i < sources->n; i++) {
    912        float sky;

  • trunk/psphot/src/test_psf_scatter.c

    r4114 r4115  
    7171
    7272    psPolynomial1D *dsxLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD);
    73     psVectorFitPolynomial1D (dsxLine, snBin, dsxBin, NULL);
     73    psVectorFitPolynomial1D (dsxLine, NULL, snBin, dsxBin, NULL);
    7474    psPolynomial1D *dsyLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD);
    75     psVectorFitPolynomial1D (dsyLine, snBin, dsyBin, NULL);
     75    psVectorFitPolynomial1D (dsyLine, NULL, snBin, dsyBin, NULL);
    7676    // these should have a slope of 1.0
    7777    psPolynomial1DDump (dsxLine);
Note: See TracChangeset for help on using the changeset viewer.