IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25528


Ignore:
Timestamp:
Sep 23, 2009, 5:31:22 PM (17 years ago)
Author:
eugene
Message:

use moments sigma if previously determinedpsphot/src/psphotReadoutFindPSF.c

Location:
branches/eam_branches/20090715/psphot/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r23442 r25528  
    4141
    4242    // construct sources and measure basic stats
    43     psArray *sources = psphotSourceStats (config, readout, detections);
     43    psArray *sources = psphotSourceStats (config, readout, detections, true);
    4444    if (!sources) return false;
    4545

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

    r25409 r25528  
    5454
    5555    // construct sources and measure basic stats
    56     psArray *sources = psphotSourceStats (config, readout, detections);
     56    psArray *sources = psphotSourceStats (config, readout, detections, true);
    5757    if (!sources) return false;
    5858

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

    r25389 r25528  
    11# include "psphotInternal.h"
    22
    3 psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections) {
     3bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources);
     4
     5psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow) {
    46
    57    bool status = false;
    … …  
    6567    }
    6668
    67     // XXX *** this section is an attempt to iteratively determine the best value for sigma of the moments weighting Gaussian
    68     {
    69         float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
    70         if (!status) return false;
    71 
    72         float sigma[4] = {1.0, 2.0, 3.0, 4.5};
    73         float Sout[4];
    74 
    75         // XXX check that this sorts by peak->SN
    76         sources = psArraySort (sources, pmSourceSortBySN);
    77 
    78         // loop over radii:
    79         for (int i = 0; i < 4; i++) {
    80 
    81             for (int j = 0; (j < sources->n) && (j < 400); j++) {
    82  
    83                 pmSource *source = sources->data[j];
    84                 psAssert (source->moments, "force moments to exist");
    85                 source->moments->nPixels = 0;
    86 
    87                 // skip faint sources for moments measurement
    88                 if (source->peak->SN < MIN_SN) {
    89                     continue;
    90                 }
    91 
    92                 // measure basic source moments
    93                 status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0);
    94             }
    95 
    96             // choose a grid scale that is a fixed fraction of the psf sigma^2
    97             psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i]));
    98             psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
    99             psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
    100 
    101             // determine the PSF parameters from the source moment values
    102             pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe);
    103             psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
    104 
    105             psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
    106             psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
    107             if (!regionMD) {
    108                 regionMD = psMetadataAlloc();
    109                 psMetadataAddMetadata (recipe, PS_LIST_TAIL, "PSF.CLUMP.REGION.000", PS_META_REPLACE, "psf clump region", regionMD);
    110                 psFree (regionMD);
    111             }
    112             psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
    113             psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
    114             psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
    115             psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
    116            
    117             // psphotVisualPlotMoments (recipe, sources);
    118 
    119             Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
     69    if (setWindow) {
     70        if (!psphotSetMomentsWindow(recipe, sources)) {
     71            psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
     72            return NULL;
    12073        }
    121 
    122         // we are looking for sigma for which Sout = 0.65 (or some other value)
    123 
    124         float Sigma = NAN;
    125         float minS = Sout[0];
    126         float maxS = Sout[0];
    127         for (int i = 0; i < 4; i++) {
    128             minS = PS_MIN(Sout[i], minS);
    129             maxS = PS_MAX(Sout[i], maxS);
    130         }
    131         if (minS > 0.65) Sigma = sigma[3];
    132         if (maxS < 0.65) Sigma = sigma[0];
    133 
    134         for (int i = 0; i < 3; i++) {
    135             if ((Sout[i] > 0.65) && (Sout[i+1] > 0.65)) continue;
    136             if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue;
    137             Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
    138         }
    139         psAssert (isfinite(Sigma), "did we miss a case?");
    140        
    141         // choose a grid scale that is a fixed fraction of the psf sigma^2
    142         psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));
    143         psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
    144         psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
    145         psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
    146         psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
    147 
    148         psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);
    14974    }
    15075
    … …  
    318243}
    319244
     245// this function attempts to iteratively determine the best value for sigma of the moments weighting Gaussian
     246bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources) {
     247
     248    bool status;
     249
     250    float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
     251    if (!status) return false;
     252
     253    // XXX move this to a config file?
     254    float sigma[4] = {1.0, 2.0, 3.0, 4.5};
     255    float Sout[4];
     256
     257    // this sorts by peak->SN
     258    sources = psArraySort (sources, pmSourceSortBySN);
     259
     260    // loop over radii:
     261    for (int i = 0; i < 4; i++) {
     262
     263        // XXX move max source number to config
     264        for (int j = 0; (j < sources->n) && (j < 400); j++) {
     265 
     266            pmSource *source = sources->data[j];
     267            psAssert (source->moments, "force moments to exist");
     268            source->moments->nPixels = 0;
     269
     270            // skip faint sources for moments measurement
     271            if (source->peak->SN < MIN_SN) {
     272                continue;
     273            }
     274
     275            // measure basic source moments
     276            status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0);
     277        }
     278
     279        // choose a grid scale that is a fixed fraction of the psf sigma^2
     280        psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i]));
     281        psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
     282        psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
     283
     284        // determine the PSF parameters from the source moment values
     285        pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe);
     286        psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
     287
     288        psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
     289        psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
     290        if (!regionMD) {
     291            regionMD = psMetadataAlloc();
     292            psMetadataAddMetadata (recipe, PS_LIST_TAIL, "PSF.CLUMP.REGION.000", PS_META_REPLACE, "psf clump region", regionMD);
     293            psFree (regionMD);
     294        }
     295        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
     296        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
     297        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
     298        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
     299           
     300        // psphotVisualPlotMoments (recipe, sources);
     301
     302        Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
     303    }
     304
     305    // we are looking for sigma for which Sout = 0.65 (or some other value)
     306
     307    float Sigma = NAN;
     308    float minS = Sout[0];
     309    float maxS = Sout[0];
     310    for (int i = 0; i < 4; i++) {
     311        minS = PS_MIN(Sout[i], minS);
     312        maxS = PS_MAX(Sout[i], maxS);
     313    }
     314    if (minS > 0.65) Sigma = sigma[3];
     315    if (maxS < 0.65) Sigma = sigma[0];
     316
     317    for (int i = 0; i < 3; i++) {
     318        if ((Sout[i] > 0.65) && (Sout[i+1] > 0.65)) continue;
     319        if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue;
     320        Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
     321    }
     322    psAssert (isfinite(Sigma), "did we miss a case?");
     323       
     324    // choose a grid scale that is a fixed fraction of the psf sigma^2
     325    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));
     326    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
     327    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
     328    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
     329    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
     330
     331    psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);
     332    return true;
     333}
     334
    320335# if (0)
    321336bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe) {
Note: See TracChangeset for help on using the changeset viewer.