IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 1, 2006, 12:24:13 PM (20 years ago)
Author:
eugene
Message:

added KAPA option and two summary plots; dropped old test code

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_01_branch/psphot/src/psphotSummaryPlots.c

    r10281 r10387  
    11# include "psphot.h"
    22
     3// the top portion of this file defines plotting functions which use kapa for plotting.
     4// if kapa is not available, these functions are defined in the bottom portion as stubs
     5// which perform NOP and return false (XXX should this be true??)
     6
     7# if (HAVE_KAPA)
     8
     9# include <kapa.h>
     10
     11// XXX not thread safe (perhaps not needed)
     12// to make this thread safe, we could check the thread ID and tie to it
     13static int kapa_fd = -1;
     14int psphotKapaOpen () {
     15
     16    if (kapa_fd == -1) {
     17        // XXX make -noX an option
     18        // XXX -noX is crashing kapa on point plotting
     19        kapa_fd = KapaOpen ("kapa", "psphot");
     20    }
     21    return kapa_fd;
     22}
     23
     24bool psphotKapaClose () {
     25   
     26    if (kapa_fd == -1) return true;
     27    KapaClose (kapa_fd);
     28    kapa_fd = -1;
     29    return true;
     30}
     31
    332// plot the sx, sy moments plane (faint and bright sources)
    4 bool psphotPlotMoments (psArray *sources) {
    5 
    6 
    7 
    8   Graphdata graphdata;
    9 
    10   // XXX make -noX an option
    11   int kapa = KiiOpen ("kapa", "psphot");
    12 
    13   KapaResize (kapa, 500, 500);
    14   KapaInitGraph (&graphdata);
    15 
    16   // examine sources to set data range
    17   graphdata.xmin = -0.05;
    18   graphdata.ymin = -0.05;
    19   graphdata.xmax = +2.05;
    20   graphdata.xmin = +2.05;
    21  
    22   KapaSetFont (kapa, "helvetica", 14);
    23   KapaBox (kapa, &graphdata);
    24   KapaSendLabel (kapa, "sigma_x|", KAPA_LABEL_XM);
    25   KapaSendLabel (kapa, "sigma_y|", KAPA_LABEL_YM);
     33bool psphotPlotMoments (pmConfig *config, pmFPAview *view, psArray *sources) {
     34
     35    // select model pixels (from output background model file, or create internal file)
     36    pmFPAfile *file = psMetadataLookupPtr (NULL, config->files, "PSPHOT.MOMENT.PLT");
     37    if (file == NULL) {
     38        psLogMsg ("psphot", 3, "skipping moments plot");
     39        return false;
     40    }
     41
     42    // pmFPAfileOpen defers disk I/O for KAPA files: just get the correct name
     43    pmFPAfileOpen (file, view, config);
     44
     45    Graphdata graphdata;
     46
     47    psLogMsg ("psphot", 3, "creating moments plot");
     48
     49    int kapa = psphotKapaOpen ();
     50    if (kapa == -1) {
     51        psError(PSPHOT_ERR_UNKNOWN, true, "failure to open kapa");
     52        return false;
     53    }
     54
     55    KapaResize (kapa, 500, 500);
     56    KapaInitGraph (&graphdata);
     57
     58    // examine sources to set data range
     59    graphdata.xmin = -0.05;
     60    graphdata.ymin = -0.05;
     61    graphdata.xmax = +2.05;
     62    graphdata.ymax = +2.05;
     63    KapaSetLimits (kapa, &graphdata);
     64 
     65    KapaSetFont (kapa, "helvetica", 14);
     66    KapaBox (kapa, &graphdata);
     67    KapaSendLabel (kapa, "&ss&h_x| (pixels)", KAPA_LABEL_XM);
     68    KapaSendLabel (kapa, "&ss&h_y| (pixels)", KAPA_LABEL_YM);
    2669               
    27   psVector *xBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    28   psVector *yBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    29   psVector *xFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    30   psVector *yFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    31 
    32   // construct the vectors
    33   int nB = 0;
    34   int nF = 0;
    35   for (int i = 0; i < sources->n; i++) {
    36     pmSource *source = sources->data[i];
    37     if (source->moments == NULL) continue;
    38    
    39     xFaint->data.F32[nF] = source->moments->Sx;
    40     yFaint->data.F32[nF] = source->moments->Sy;
    41     nF++;
    42    
    43     // XXX make this a user-defined cutoff
    44     if (source->moments->SN < 25) continue;
    45 
    46     xBright->data.F32[nB] = source->moments->Sx;
    47     yBright->data.F32[nB] = source->moments->Sy;
    48     nB++;
    49   }
    50   xFaint->n = nF;
    51   yFaint->n = nF;
    52 
    53   xBright->n = nB;
    54   yBright->n = nB;
    55  
    56   graphdata.color = KapaColorByName ("black");
    57   graphdata.ptype = 0;
    58   graphdata.size = 0.3;
    59   graphdata.style = 2;
    60   KapaPrepPlot (kapa, nF, &graphdata);
    61   KapaPlotVector (kapa, nF, xFaint->data.F32);
    62   KapaPlotVector (kapa, nF, yFaint->data.F32);
    63  
    64   graphdata.color = KapaColorByName ("red");
    65   graphdata.ptype = 2;
    66   graphdata.size = 0.5;
    67   graphdata.style = 2;
    68   KapaPrepPlot (kapa, nB, &graphdata);
    69   KapaPlotVector (kapa, nB, xBright->data.F32);
    70   KapaPlotVector (kapa, nB, yBright->data.F32);
    71 
    72   KapaPNG (kapa, "test.png");
    73   return true;
    74 }
     70    psVector *xBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     71    psVector *yBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     72    psVector *xFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     73    psVector *yFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     74
     75    // construct the vectors
     76    int nB = 0;
     77    int nF = 0;
     78    for (int i = 0; i < sources->n; i++) {
     79        pmSource *source = sources->data[i];
     80        if (source->moments == NULL) continue;
     81   
     82        xFaint->data.F32[nF] = source->moments->Sx;
     83        yFaint->data.F32[nF] = source->moments->Sy;
     84        nF++;
     85   
     86        // XXX make this a user-defined cutoff
     87        if (source->moments->SN < 25) continue;
     88
     89        xBright->data.F32[nB] = source->moments->Sx;
     90        yBright->data.F32[nB] = source->moments->Sy;
     91        nB++;
     92    }
     93    xFaint->n = nF;
     94    yFaint->n = nF;
     95
     96    xBright->n = nB;
     97    yBright->n = nB;
     98 
     99    graphdata.color = KapaColorByName ("black");
     100    graphdata.ptype = 0;
     101    graphdata.size = 0.3;
     102    graphdata.style = 2;
     103    KapaPrepPlot (kapa, nF, &graphdata);
     104    KapaPlotVector (kapa, nF, xFaint->data.F32);
     105    KapaPlotVector (kapa, nF, yFaint->data.F32);
     106 
     107    graphdata.color = KapaColorByName ("red");
     108    graphdata.ptype = 2;
     109    graphdata.size = 0.5;
     110    graphdata.style = 2;
     111    KapaPrepPlot (kapa, nB, &graphdata);
     112    KapaPlotVector (kapa, nB, xBright->data.F32);
     113    KapaPlotVector (kapa, nB, yBright->data.F32);
     114
     115    psLogMsg ("psphot", 3, "saving plot to %s", file->filename);
     116    KapaPNG (kapa, file->filename);
     117
     118    psFree (xBright);
     119    psFree (yBright);
     120    psFree (xFaint);
     121    psFree (yFaint);
     122
     123    return true;
     124}
     125
     126// plot the sx, sy, sxy as vector field,
     127// plot the PSF measured sx, sy, sxy as vector field
     128// pull the sources from the config / file?
     129bool psphotPlotPSFModel (pmConfig *config, pmFPAview *view, psArray *sources) {
     130
     131    // select model pixels (from output background model file, or create internal file)
     132    pmFPAfile *file = psMetadataLookupPtr (NULL, config->files, "PSPHOT.PSFMODEL.PLT");
     133    if (file == NULL) {
     134        psLogMsg ("psphot", 3, "skipping psf model plot");
     135        return false;
     136    }
     137
     138    // pmFPAfileOpen defers disk I/O for KAPA files: just get the correct name
     139    pmFPAfileOpen (file, view, config);
     140
     141    Graphdata graphdata;
     142
     143    psLogMsg ("psphot", 3, "creating psf model plot");
     144
     145    int kapa = psphotKapaOpen ();
     146    if (kapa == -1) {
     147        psError(PSPHOT_ERR_UNKNOWN, true, "failure to open kapa");
     148        return false;
     149    }
     150
     151    // XXX make the aspect-ratio match the image
     152    KapaResize (kapa, 800, 800);
     153    KapaInitGraph (&graphdata);
     154
     155    psVector *xMNT = psVectorAllocEmpty (2*sources->n, PS_TYPE_F32);
     156    psVector *yMNT = psVectorAllocEmpty (2*sources->n, PS_TYPE_F32);
     157    psVector *xPSF = psVectorAllocEmpty (2*sources->n, PS_TYPE_F32);
     158    psVector *yPSF = psVectorAllocEmpty (2*sources->n, PS_TYPE_F32);
     159    psVector *xMIN = psVectorAllocEmpty (2*sources->n, PS_TYPE_F32);
     160    psVector *yMIN = psVectorAllocEmpty (2*sources->n, PS_TYPE_F32);
     161
     162    // construct the plot vectors
     163    int nMNT = 0;
     164    int nPSF = 0;
     165    int nMIN = 0;
     166    float dx = 0;
     167    float dy = 0;
     168    float scale = 10;
     169    for (int i = 0; i < sources->n; i++) {
     170        pmSource *source = sources->data[i];
     171        if (source->moments == NULL) continue;
     172        if (source->moments->SN < 25) continue;
     173        if (source->type != PM_SOURCE_TYPE_STAR) continue;
     174   
     175        pmModel *model = source->modelPSF;
     176        if (model == NULL) continue;
     177
     178        psF32 *PAR = model->params->data.F32;
     179
     180        psEllipseMoments moments;
     181        moments.x2 = source->moments->Sx;
     182        moments.y2 = source->moments->Sy;
     183        moments.xy = source->moments->Sxy;
     184
     185        psEllipseShape shape;
     186        shape.sx  = PAR[PM_PAR_SXX] / sqrt(2.0);
     187        shape.sy  = PAR[PM_PAR_SYY] / sqrt(2.0);
     188        shape.sxy = PAR[PM_PAR_SXY];
     189
     190        // force the axis ratio to be < 20.0
     191        psEllipseAxes axes_mnt = psEllipseMomentsToAxes (moments, 20.0);
     192        psEllipseAxes axes_psf = psEllipseShapeToAxes (shape, 20.0);
     193
     194        // moments major axis
     195        dx = scale*axes_mnt.major*cos(axes_mnt.theta);
     196        dy = scale*axes_mnt.major*sin(axes_mnt.theta);
     197        xMNT->data.F32[nMNT] = PAR[PM_PAR_XPOS] - dx;
     198        yMNT->data.F32[nMNT] = PAR[PM_PAR_YPOS] - dy;
     199        nMNT++;
     200        xMNT->data.F32[nMNT] = PAR[PM_PAR_XPOS] + dx;
     201        yMNT->data.F32[nMNT] = PAR[PM_PAR_YPOS] + dy;
     202        nMNT++;
     203   
     204        // psf major axis
     205        dx = scale*axes_psf.major*cos(axes_psf.theta);
     206        dy = scale*axes_psf.major*sin(axes_psf.theta);
     207        xPSF->data.F32[nPSF] = PAR[PM_PAR_XPOS] - dx;
     208        yPSF->data.F32[nPSF] = PAR[PM_PAR_YPOS] - dy;
     209        nPSF++;
     210        xPSF->data.F32[nPSF] = PAR[PM_PAR_XPOS] + dx;
     211        yPSF->data.F32[nPSF] = PAR[PM_PAR_YPOS] + dy;
     212        nPSF++;
     213
     214        // minor axis (to show size)
     215        dy = +scale*axes_psf.minor*cos(axes_psf.theta);
     216        dx = -scale*axes_psf.minor*sin(axes_psf.theta);
     217        xMIN->data.F32[nMIN] = PAR[PM_PAR_XPOS] - dx;
     218        yMIN->data.F32[nMIN] = PAR[PM_PAR_YPOS] - dy;
     219        nMIN++;
     220        xMIN->data.F32[nMIN] = PAR[PM_PAR_XPOS] + dx;
     221        yMIN->data.F32[nMIN] = PAR[PM_PAR_YPOS] + dy;
     222        nMIN++;
     223
     224        graphdata.xmin = PS_MIN(graphdata.xmin, PAR[PM_PAR_XPOS]);
     225        graphdata.xmax = PS_MAX(graphdata.xmax, PAR[PM_PAR_XPOS]);
     226        graphdata.ymin = PS_MIN(graphdata.ymin, PAR[PM_PAR_YPOS]);
     227        graphdata.ymax = PS_MAX(graphdata.ymax, PAR[PM_PAR_YPOS]);
     228    }
     229    xMNT->n = yMNT->n = nMNT;
     230    xPSF->n = yPSF->n = nPSF;
     231    xMIN->n = yMIN->n = nMIN;
     232
     233    float range;
     234    range = graphdata.xmax - graphdata.xmin;
     235    graphdata.xmax += 0.05*range;
     236    graphdata.xmin -= 0.05*range;
     237    range = graphdata.ymax - graphdata.ymin;
     238    graphdata.ymax += 0.05*range;
     239    graphdata.ymin -= 0.05*range;
     240
     241    // XXX set the plot range to match the image
     242    KapaSetLimits (kapa, &graphdata);
     243 
     244    KapaSetFont (kapa, "helvetica", 14);
     245    KapaBox (kapa, &graphdata);
     246    KapaSendLabel (kapa, "x (pixels)", KAPA_LABEL_XM);
     247    KapaSendLabel (kapa, "y (pixels)", KAPA_LABEL_YM);
     248    KapaSendLabel (kapa, "vector is major axis (scale by 20) : black are moments, blue are psf model, red is psf minor axis", KAPA_LABEL_XP);
     249               
     250    graphdata.color = KapaColorByName ("black");
     251    graphdata.ptype = 100;
     252    graphdata.size = 0.3;
     253    graphdata.style = 2;
     254    KapaPrepPlot (kapa, nMNT, &graphdata);
     255    KapaPlotVector (kapa, nMNT, xMNT->data.F32);
     256    KapaPlotVector (kapa, nMNT, yMNT->data.F32);
     257 
     258    graphdata.color = KapaColorByName ("blue");
     259    graphdata.ptype = 100;
     260    graphdata.size = 0.5;
     261    graphdata.style = 2;
     262    KapaPrepPlot (kapa, nPSF, &graphdata);
     263    KapaPlotVector (kapa, nPSF, xPSF->data.F32);
     264    KapaPlotVector (kapa, nPSF, yPSF->data.F32);
     265
     266    graphdata.color = KapaColorByName ("red");
     267    graphdata.ptype = 100;
     268    graphdata.size = 0.5;
     269    graphdata.style = 2;
     270    KapaPrepPlot (kapa, nMIN, &graphdata);
     271    KapaPlotVector (kapa, nMIN, xMIN->data.F32);
     272    KapaPlotVector (kapa, nMIN, yMIN->data.F32);
     273
     274    psLogMsg ("psphot", 3, "saving plot to %s", file->filename);
     275    KapaPNG (kapa, file->filename);
     276
     277    psFree (xMNT);
     278    psFree (yMNT);
     279    psFree (xPSF);
     280    psFree (yPSF);
     281    psFree (xMIN);
     282    psFree (yMIN);
     283
     284    return true;
     285}
     286
     287# else
     288
     289bool psphotPlotMoments (pmConfig *config, pmFPAview *view, psArray *sources) {
     290    psLogMsg ("psphot", 3, "skipping moments plot");
     291    return true;
     292}
     293
     294bool psphotPlotPSFModel (pmConfig *config, pmFPAview *view, psArray *sources) {
     295    psLogMsg ("psphot", 3, "skipping psf model plot");
     296    return true;
     297}
     298
     299# endif
Note: See TracChangeset for help on using the changeset viewer.