IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 13086


Ignore:
Timestamp:
Apr 30, 2007, 3:51:22 PM (19 years ago)
Author:
eugene
Message:

added psphotSourceSize, not fully integrated

Location:
trunk/psphot/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/Makefile.am

    r13035 r13086  
    5656        psphotMosaicSubimage.c   \
    5757        psphotMakeResiduals.c    \
     58        psphotSourceSize.c       \
    5859        psphotAddNoise.c
    5960

  • trunk/psphot/src/psphot.h

    r13052 r13086  
    108108bool            psphotSetState (pmSource *source, bool curState);
    109109bool            psphotDeblendSatstars (psArray *sources, psMetadata *recipe);
     110bool            psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe);
    110111
    111112bool            psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf);

  • trunk/psphot/src/psphotReadout.c

    r13035 r13086  
    128128    psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
    129129    if (dump) psphotSaveImage (NULL, readout->image,  "image.v1.fits");
     130
     131    // XXX test the CR/EXT measurement
     132    psphotSourceSize (readout, sources, recipe);
     133
    130134    if (!strcasecmp (breakPt, "ENSEMBLE")) {
    131135        goto finish;

  • trunk/psphot/src/psphotSourceSize.c

    r13081 r13086  
    11# include "psphotInternal.h"
     2# include <gsl/gsl_sf_gamma.h>
    23
    3 // we need to call this function after sources have been fitted to the PSF
    4 // model and subtracted.  this function examines the 9 pixels in the 3x3 square containing the peak
    5 // and compares the observed flux to the model.  For the 4 pixel lines which contain the peak, we
    6 // calculate the value
     4// we need to call this function after sources have been fitted to the PSF model and
     5// subtracted.  this function examines the 9 pixels in the 3x3 square containing the peak and
     6// compares the observed flux to the model.
    77
    8 /*
    9  *
    10  *
    11  *
    12  *
    13  *
    14  */
     8bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe) {
    159
    16 bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add) {
     10    FILE *f = fopen ("srcsize.dat", "w");
    1711
    1812    // loop over all source
    … …  
    2317        if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
    2418
    25         psF32 *resid = source->pixels->data.F32;
    26         psF32 *weight = source->weight->data.F32;
    27         psU8 *mask = source->maskObj->data.U8;
     19        psF32 **resid = source->pixels->data.F32;
     20        psF32 **weight = source->weight->data.F32;
     21        psU8 **mask = source->maskObj->data.U8;
    2822
    29         int xPeak = source->peak->xf - source->pixels->col0;
    30         int yPeak = source->peak->yf - source->pixels->row0;
     23        int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
     24        int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
    3125
    32         // XXX need to deal with edge peaks...
     26        // XXX for now, just skip any sources with masked pixels
     27        bool keep = true;
     28        for (int iy = -1; (iy <= +1) && keep; iy++) {
     29            for (int ix = -1; (ix <= +1) && keep; ix++) {
     30                if (mask[yPeak+iy][xPeak+ix]) { keep &= false; }
     31            }
     32        }
     33        if (!keep) continue;
     34
     35        // XXX need to deal with edge peaks... and mask
    3336        float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1];
    3437        float cY = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak+0] - resid[yPeak+1][xPeak+0];
    … …  
    3639        float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1];
    3740
    38         float dcX = 2*weight[yPeak][xPeak] - weight[yPeak+0][xPeak-1] - weight[yPeak+0][xPeak+1];
    39         float dcY = 2*weight[yPeak][xPeak] - weight[yPeak+1][xPeak+0] - weight[yPeak+1][xPeak+0];
    40         float dcL = 2*weight[yPeak][xPeak] - weight[yPeak-1][xPeak-1] - weight[yPeak+1][xPeak+1];
    41         float dcR = 2*weight[yPeak][xPeak] - weight[yPeak+1][xPeak-1] - weight[yPeak-1][xPeak+1];
     41        float dcX = 4*weight[yPeak][xPeak] + weight[yPeak+0][xPeak-1] + weight[yPeak+0][xPeak+1];
     42        float dcY = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak+0] + weight[yPeak+1][xPeak+0];
     43        float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1];
     44        float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1];
     45
     46        float nX = cX / sqrt(dcX);
     47        float nY = cY / sqrt(dcY);
     48        float nL = cL / sqrt(dcL);
     49        float nR = cR / sqrt(dcR);
     50
     51        float chisq = PS_SQR (nX) + PS_SQR (nY) + PS_SQR (nX) + PS_SQR (nR);
     52
     53        // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
     54        // Ndof = 4 ? (four measurements, no free parameters)
     55        float Ppsf = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
     56
     57        float fCR = 0.0;
     58        float fEXT = 0.0;
     59        int nCR = 0;
     60        int nEXT = 0;
     61        if (nX > 0.0) {
     62            fCR += nX;
     63            nCR ++;
     64        } else {
     65            fEXT += nX;
     66            nEXT ++;
     67        }
     68        if (nY > 0.0) {
     69            fCR += nY;
     70            nCR ++;
     71        } else {
     72            fEXT += nY;
     73            nEXT ++;
     74        }
     75        if (nL > 0.0) {
     76            fCR += nL;
     77            nCR ++;
     78        } else {
     79            fEXT += nL;
     80            nEXT ++;
     81        }
     82        if (nR > 0.0) {
     83            fCR += nR;
     84            nCR ++;
     85        } else {
     86            fEXT += nR;
     87            nEXT ++;
     88        }
     89        fCR = (nCR > 0) ? fCR / nCR : 0.0;
     90        fEXT = (nEXT > 0) ? fEXT / nEXT : 0.0;
    4291
    4392        // XXX plot this as a function of flux / magnitude?
    4493        // XXX overlay outliers on an image?
    4594       
    46         // some possible classes:
    47         // all 4 curvatures are highly negative : CR
    48         // all 4 curvatures are highly positive : EXT
     95        fprintf (f, "%f %f %f   %f %f   %f %f\n", source->peak->xf, source->peak->yf, source->peak->flux,
     96                 chisq, Ppsf, fCR, fEXT);
     97    }
     98    fclose (f);
     99    return true;
     100}
    49101
    50         // at least 2 are significantly negative, none are significantly positive : CR
    51         // at least 2 are significantly positive, none are significantly negative : EXT
    52102
    53         // any are significantly negative, some may be significantly positive : CR
    54         // any are significantly positive, none may be significantly positive : EXT
     103// some possible classes:
     104// all 4 curvatures are highly negative : CR
     105// all 4 curvatures are highly positive : EXT
     106
     107// at least 2 are significantly negative, none are significantly positive : CR
     108// at least 2 are significantly positive, none are significantly negative : EXT
     109
     110// any are significantly negative, some may be significantly positive : CR
     111// any are significantly positive, none may be significantly positive : EXT
    55112       
    56         /* Nn  Np  No
    57            4   0   0   CR_1
    58            3   1   0   CR_1
    59            3   0   1   CR_1
    60            2   2   0   CR_2
    61            2   1   1   CR_2
    62            2   0   2   CR_2
    63            1   3   0   CR_3
    64            1   2   1   CR_3
    65            1   1   2   CR_3
    66            1   0   3   CR_3
    67            0   4   0   EXT
    68            0   3   1   EXT
    69            0   2   2   EXT
    70            0   1   3   PSF
    71            0   0   4   PSF
    72         */
     113/* Nn  Np  No
     114   4   0   0   CR_1
     115   3   1   0   CR_1
     116   3   0   1   CR_1
     117   2   2   0   CR_2
     118   2   1   1   CR_2
     119   2   0   2   CR_2
     120   1   3   0   CR_3
     121   1   2   1   CR_3
     122   1   1   2   CR_3
     123   1   0   3   CR_3
     124   0   4   0   EXT
     125   0   3   1   EXT
     126   0   2   2   EXT
     127   0   1   3   PSF
     128   0   0   4   PSF
     129*/
    73130
    74         /* I can write the formal probability that the 4 measurements are consistent with a PSF
    75          * based on how large the error values are (Nsigma -> erf -> P(Nsigma)).
    76          * I should examine this value as a function of flux and also examine the distribution of
    77          * the probabilities for the PSF sources.
    78          */
     131/* Alternatively, write a f(CR) = Sum(nX,etc if >0) */
     132         
     133
     134/* I can write the formal probability that the 4 measurements are consistent with a PSF
     135 * based on how large the error values are (Nsigma -> erf -> P(Nsigma)).
     136 * I should examine this value as a function of flux and also examine the distribution of
     137 * the probabilities for the PSF sources.
     138 */
     139
Note: See TracChangeset for help on using the changeset viewer.