IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 12665


Ignore:
Timestamp:
Mar 28, 2007, 5:20:02 PM (19 years ago)
Author:
eugene
Message:

adding source plots

Location:
trunk/psphot/src
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/Makefile.am

    r12633 r12665  
    4949        psphotMergeSources.c     \
    5050        psphotReadoutCleanup.c   \
     51        psphotSourcePlots.c      \
     52        psphotRadialPlot.c       \
     53        psphotMosaicSubimage.c   \
    5154        psphotAddNoise.c         \
    5255        psphotCleanup.c         
    5356
    54 #       psphotSourcePlots.c     
    55 #       psphotRadialPlot.c       
    5657
    5758include_HEADERS = \

  • trunk/psphot/src/psphot.h

    r12627 r12665  
    104104bool psphotSetHeaderNstars (psMetadata *recipe, psArray *sources);
    105105bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add);
    106 bool psphotRadialPlot (const char *filename, pmSource *source);
     106bool psphotRadialPlot (int *kapa, char *filename, pmSource *source);
    107107bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe);
     108bool psphotMosaicSubimage (psImage *outImage, pmSource *source, int Xo, int Yo, int DX, int DY);
     109
     110bool psphotAddWithTest (pmSource *source, bool useState);
     111bool psphotSubWithTest (pmSource *source, bool useState);
     112bool psphotSetState (pmSource *source, bool curState);

  • trunk/psphot/src/psphotChoosePSF.c

    r11702 r12665  
    7474
    7575    // select the candidate PSF stars (pointers to original sources)
    76     for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) {
     76    for (int i = 0; i < sources->n; i++) {
    7777        pmSource *source = sources->data[i];
    78         if (source->mode & PM_SOURCE_MODE_PSFSTAR) psArrayAdd (stars, 200, source);
     78        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
     79            // keep NSTARS PSF stars, unmark the rest
     80            if (stars->n < NSTARS) {
     81                psArrayAdd (stars, 200, source);
     82            } else {
     83                source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
     84            }
     85        }
    7986    }
    8087    psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n);

  • trunk/psphot/src/psphotRadialPlot.c

    r12627 r12665  
    66# include <kapa.h>
    77
    8 bool psphotRadialPlot (const char *filename, pmSource *source) {
     8static int nCount = 0;
     9
     10bool psphotRadialPlot (int *kapa, char *filename, pmSource *source) {
     11
     12    Graphdata graphdata;
     13
     14    // only plot 50 stars for now...
     15    if (nCount > 00) {
     16        if (*kapa != 0) {
     17            KiiClose (*kapa);
     18            *kapa = 0;
     19        }
     20        return true;
     21    }
    922
    1023    // XXX get the 'showWindow' option from the recipes somewhere
    11     int kapa = pmKapaOpen (false);
    12     if (kapa == -1) {
     24    // XXX 'showWindow = false' is broken
     25    if (*kapa == 0) {
     26        *kapa = pmKapaOpen (true);
     27        KapaResize (*kapa, 500, 500);
     28        unlink (filename);
     29    }
     30    if (*kapa == -1) {
    1331        psError(PSPHOT_ERR_UNKNOWN, true, "failure to open kapa");
    1432        return false;
    1533    }
    1634
    17     KapaResize (kapa, 500, 500);
    1835    KapaInitGraph (&graphdata);
     36    KapaClear(*kapa, true);
    1937
    2038    // examine sources to set data range
    21     graphdata.xmin = -0.05;
     39    graphdata.xmin =  -0.05;
     40    graphdata.xmax = +30.05;
    2241    graphdata.ymin = -0.05;
    23     graphdata.xmax = +2.05;
    24     graphdata.ymax = +2.05;
    25     KapaSetLimits (kapa, &graphdata);
     42    graphdata.ymax = +5.05;
     43    KapaSetLimits (*kapa, &graphdata);
    2644 
    27     KapaSetFont (kapa, "helvetica", 14);
    28     KapaBox (kapa, &graphdata);
    29     KapaSendLabel (kapa, "&ss&h_x| (pixels)", KAPA_LABEL_XM);
    30     KapaSendLabel (kapa, "&ss&h_y| (pixels)", KAPA_LABEL_YM);
     45    KapaSetFont (*kapa, "helvetica", 14);
     46    KapaBox (*kapa, &graphdata);
     47    KapaSendLabel (*kapa, "radius (pixels)", KAPA_LABEL_XM);
     48    KapaSendLabel (*kapa, "log flux (counts)", KAPA_LABEL_YM);
    3149               
    3250    int nPts = source->pixels->numRows * source->pixels->numCols;
    33     psVector *radius = psVectorAllocEmpty (nPts, PS_TYPE_F32);
    34     psVector *lflux = psVectorAllocEmpty (nPts, PS_TYPE_F32);
    35     psVector *rmask = psVectorAllocEmpty (nPts, PS_TYPE_F32);
    36     psVector *lfmask = psVectorAllocEmpty (nPts, PS_TYPE_F32);
     51    psVector *rg = psVectorAllocEmpty (nPts, PS_TYPE_F32);
     52    psVector *fg = psVectorAllocEmpty (nPts, PS_TYPE_F32);
     53    psVector *rb = psVectorAllocEmpty (nPts, PS_TYPE_F32);
     54    psVector *fb = psVectorAllocEmpty (nPts, PS_TYPE_F32);
    3755
    38     int n = 0;
    39     int nm = 0;
     56    int ng = 0;
     57    int nb = 0;
     58    float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
     59    float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
    4060    for (int iy = 0; iy < source->pixels->numRows; iy++) {
    4161        for (int ix = 0; ix < source->pixels->numCols; ix++) {
    42            
    4362            if (source->mask->data.U8[iy][ix]) {
    44                 rmask->data.F32[nm] = RADIUS;
    45                 lfmask->data.F32[nm] = log10(source->pixels->data.F32[iy][ix]);
    46                 nm++;
     63                rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
     64                fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]);
     65                nb++;
    4766            } else {
    48                 radius->data.F32[n] = RADIUS;
    49                 lflux->data.F32[n] = log10(source->pixels->data.F32[iy][ix]);
    50                 n++;
     67                rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
     68                fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]);
     69                ng++;
    5170            }
    5271        }
    … …  
    5675
    5776    graphdata.color = KapaColorByName ("black");
     77    graphdata.ptype = 2;
     78    graphdata.size = 0.5;
     79    graphdata.style = 2;
     80    KapaPrepPlot (*kapa, ng, &graphdata);
     81    KapaPlotVector (*kapa, ng, rg->data.F32);
     82    KapaPlotVector (*kapa, ng, fg->data.F32);
     83
     84    graphdata.color = KapaColorByName ("red");
    5885    graphdata.ptype = 0;
    5986    graphdata.size = 0.3;
    6087    graphdata.style = 2;
    61     KapaPrepPlot (kapa, nm, &graphdata);
    62     KapaPlotVector (kapa, nm, rmask->data.F32);
    63     KapaPlotVector (kapa, nm, lfmask->data.F32);
     88    KapaPrepPlot (*kapa, nb, &graphdata);
     89    KapaPlotVector (*kapa, nb, rb->data.F32);
     90    KapaPlotVector (*kapa, nb, fb->data.F32);
    6491 
    65     graphdata.color = KapaColorByName ("red");
    66     graphdata.ptype = 2;
    67     graphdata.size = 0.5;
    68     graphdata.style = 2;
    69     KapaPrepPlot (kapa, n, &graphdata);
    70     KapaPlotVector (kapa, n, radius->data.F32);
    71     KapaPlotVector (kapa, n, lflux->data.F32);
     92    psLogMsg ("psphot", 3, "saving plot to %s", filename);
    7293
    73     psLogMsg ("psphot", 3, "saving plot to %s", file->filename);
    74     KapaPS (kapa, false, KAPA_PS_NEWPAGE, filename, "none");
     94    char pagename[16];
     95    sprintf (pagename, "%02d", nCount);
     96    if (nCount == 0) {
     97        KiiPS (*kapa, false, KAPA_PS_NEWPLOT, filename, pagename);
     98    } else {
     99        KiiPS (*kapa, false, KAPA_PS_NEWPAGE, filename, pagename);
     100    }
    75101
    76     psFree (radius);
    77     psFree (lflux);
    78     psFree (rmask);
    79     psFree (lfmask);
     102    psFree (rg);
     103    psFree (fg);
     104    psFree (rb);
     105    psFree (fb);
    80106
     107    nCount ++;
     108    return true;
    81109}
     110
     111# else
     112
     113bool psphotRadialPlot (int *kapa, const char *filename, pmSource *source) {
     114    psLogMsg ("psphot", 3, "skipping source radial plots");
     115    return true;
     116}
     117
     118# endif

  • trunk/psphot/src/psphotReadout.c

    r12632 r12665  
    108108    }
    109109
     110    // plot positive + negative sources (replace, then remove)
     111    psphotSourcePlots (readout, sources, recipe);
     112
    110113    // non-linear PSF and EXT fit to brighter sources
    111114    psphotBlendFit (readout, sources, recipe, psf);
    … …  
    113116    // replace all sources
    114117    psphotReplaceAll (sources);
     118
     119    // plot positive sources
     120    // psphotSourcePlots (readout, sources, recipe);
    115121
    116122    // linear PSF fit to remaining peaks

  • trunk/psphot/src/psphotReplaceUnfit.c

    r12625 r12665  
    5353    return true;
    5454}
     55
     56// add or sub source replace or if the source has
     57bool psphotAddWithTest (pmSource *source, bool useState) {
     58
     59    // what is current state? (true : add; false : sub)
     60    bool state = !(source->mode & PM_SOURCE_MODE_SUBTRACTED);
     61    if (state && useState) return true;
     62
     63    // select appropriate model
     64    pmModel *model = pmSourceGetModel (NULL, source);
     65    if (model == NULL) return false;  // model must be defined
     66   
     67    psTrace ("psphot", 3, "replacing object at %f,%f\n",
     68             model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     69   
     70    // replace the model if 1) state says it is missing or 2) useState is false (just do it)
     71    if (!state || !useState) {
     72        pmModelAdd (source->pixels, source->mask, model, false, false);
     73    }
     74    return true;
     75}
     76
     77// add or sub source replace or if the source has
     78bool psphotSubWithTest (pmSource *source, bool useState) {
     79
     80    // what is current state? (true : sub; false : add)
     81    bool state = (source->mode & PM_SOURCE_MODE_SUBTRACTED);
     82    if (state && useState) return true;
     83
     84    // select appropriate model
     85    pmModel *model = pmSourceGetModel (NULL, source);
     86    if (model == NULL) return false;  // model must be defined
     87   
     88    psTrace ("psphot", 3, "replacing object at %f,%f\n",
     89             model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     90   
     91    // replace the model if 1) state says it is missing or 2) useState is false (just do it)
     92    if (!state || !useState) {
     93        pmModelSub (source->pixels, source->mask, model, false, false);
     94    }
     95    return true;
     96}
     97
     98// add or sub source replace or if the source has
     99bool psphotSetState (pmSource *source, bool curState) {
     100
     101    // what is desired state? (true : add; false : sub)
     102    bool newState = !(source->mode & PM_SOURCE_MODE_SUBTRACTED);
     103    if (curState == newState) return true;
     104
     105    // select appropriate model
     106    pmModel *model = pmSourceGetModel (NULL, source);
     107    if (model == NULL) return false;  // model must be defined
     108   
     109    psTrace ("psphot", 3, "replacing object at %f,%f\n",
     110             model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
     111   
     112    if (curState && !newState) {
     113        pmModelSub (source->pixels, source->mask, model, false, false);
     114    }
     115    if (newState && !curState) {
     116        pmModelAdd (source->pixels, source->mask, model, false, false);
     117    }
     118    return true;
     119}

  • trunk/psphot/src/psphotSourcePlots.c

    r12627 r12665  
    33bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe) {
    44
    5     bool status = false;
     5    // bool status = false;
    66    psTimerStart ("psphot");
    77
    88    // the source images are written to an image 10x the size of a PSF object
    9     float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    10     PS_ASSERT (status);
     9    // float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
     10    // PS_ASSERT (status, false);
     11
     12    int DX = 21;
     13    int DY = 21;
    1114
    1215    // examine PSF sources in S/N order (brightest first)
    … …  
    1619    int dX = 0;                         // starting corner of next box
    1720    int dY = 0;                         // height of row so far
    18     int NX = 10*(2*OUTER+1);            // full width of output image
     21    int NX = 20*DX;                     // full width of output image
    1922    int NY = 0;                         // total height of output image
    2023
    2124    // first, examine the PSF and SAT stars:
    22     // - generate radial plots (PS plots)
    2325    // - determine bounding boxes for summary image
    2426    for (int i = 0; i < sources->n; i++) {
    … …  
    3234
    3335        // how does this subimage get placed into the output image?
    34         if (dX + source->pixels->numCols >= NX) {
     36        // DX = source->pixels->numCols
     37        // DY = source->pixels->numRows
     38
     39        if (dX + DX > NX) {
    3540            // too wide for the rest of this row
    3641            if (dX == 0) {
    3742                // alone on this row
    38                 NY += source->pixels->numRows;
     43                NY += DY;
    3944                dX = 0;
    4045                dY = 0;
    … …  
    4247                // start the next row
    4348                NY += dY;
    44                 dX = source->pixels->numCols;
    45                 dY = source->pixels->numRows;
     49                dX = DX;
     50                dY = DY;
    4651            }
    4752        } else {
    4853            // extend this row
    49             dX += source->pixels->numCols;
    50             dY = PS_MAX (dY, source->pixels->numRows);
     54            dX += DX;
     55            dY = PS_MAX (dY, DY);
    5156        }
    52         psphotRadialPlot ("radial.plots.ps", source);
    5357    }
    5458
    … …  
    5963    int Xo = 0;                         // starting corner of next box
    6064    int Yo = 0;                         // starting corner of next box
    61     int dY = 0;                         // height of row so far
     65    dY = 0;                             // height of row so far
     66
     67    int nPSF = 0;
     68    int nSAT = 0;
     69    int kapa = 0;                       // file descriptor for plotting routine
    6270
    6371    // first, examine the PSF and SAT stars:
    6472    // - generate radial plots (PS plots)
    65     // - determine bounding boxes for summary image
     73    // - create output image array
    6674    for (int i = 0; i < sources->n; i++) {
    6775
    … …  
    6977
    7078        bool keep = false;
    71         keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
    72         keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
     79        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
     80            nPSF ++;
     81            keep = true;
     82        }
     83        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     84            nSAT ++;
     85            keep = true;
     86        }           
    7387        if (!keep) continue;
    7488
    7589        // how does this subimage get placed into the output image?
    76         if (Xo + source->pixels->numCols >= NX) {
     90        // DX = source->pixels->numCols
     91        // DY = source->pixels->numRows
     92
     93        if (Xo + DX > NX) {
    7794            // too wide for the rest of this row
    7895            if (Xo == 0) {
    79                 // alone on this row
    80                 PlaceImage (source, Xo, Yo);
    81                 Yo += source->pixels->numRows;
     96                // place source alone on this row
     97                psphotAddWithTest (source, true); // replace source if subtracted
     98                psphotRadialPlot (&kapa, "radial.plots.ps", source);
     99                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY);
     100
     101                psphotSubWithTest (source, false); // remove source (force)
     102                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY);
     103
     104                psphotSetState (source, false); // replace source (has been subtracted)
     105                Yo += DY;
    82106                Xo = 0;
    83107                dY = 0;
    … …  
    86110                Yo += dY;
    87111                Xo = 0;
    88                 PlaceImage (source, Xo, Yo);
    89                 Xo = source->pixels->numCols;
    90                 dY = source->pixels->numRows;
     112                psphotAddWithTest (source, true); // replace source if subtracted
     113                psphotRadialPlot (&kapa, "radial.plots.ps", source);
     114                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY);
     115
     116                psphotSubWithTest (source, false); // remove source (force)
     117                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY);
     118                psphotSetState (source, false); // replace source (has been subtracted)
     119
     120                Xo = DX;
     121                dY = DY;
    91122            }
    92123        } else {
    93124            // extend this row
    94             PlaceImage (source, Xo, Yo);
    95             Xo += source->pixels->numCols;
    96             dY = PS_MAX (dY, source->pixels->numRows);
     125            psphotAddWithTest (source, true); // replace source if subtracted
     126            psphotRadialPlot (&kapa, "radial.plots.ps", source);
     127            psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY);
     128
     129            psphotSubWithTest (source, false); // remove source (force)
     130            psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY);
     131            psphotSetState (source, false); // replace source (has been subtracted)
     132
     133            Xo += DX;
     134            dY = PS_MAX (dY, DY);
    97135        }
    98136    }
    99137
    100     psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d failed: %f sec\n", sources->n, Nmoments, Nfail, psTimerMark ("psphot"));
     138    psphotSaveImage (NULL, outpos, "outpos.fits");
     139    psphotSaveImage (NULL, outsub, "outsub.fits");
     140    psLogMsg ("psphot", PS_LOG_INFO, "plotted %d sources (%d psf, %d sat): %f sec\n", nPSF + nSAT, nPSF, nSAT, psTimerMark ("psphot"));
     141
     142    psFree (outpos);
     143    psFree (outsub);
    101144    return (sources);
    102145}
Note: See TracChangeset for help on using the changeset viewer.