IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 19862


Ignore:
Timestamp:
Oct 2, 2008, 11:00:57 AM (18 years ago)
Author:
eugene
Message:

adding visualizations

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branch_20080926/psphot/src/psphotVisual.c

    r19762 r19862  
    11# include "psphotInternal.h"
     2bool psphotVisualPlotMoments (pmConfig *config, const pmFPAview *view, psArray *sources);
     3
     4// this function displays representative images as the psphot analysis progresses:
     5// 0 : image, 1 : weight
     6// 0 : backsub, 1 : weight, 2 : backgnd
     7// 0 : backsub, 1 : weight, 2 : signif
     8// (overlay peaks on images)
     9// (overlay footprints on images)
     10// (overlay moments on images)
     11// (overlay rough class on images)
     12// 0 : backsub, 1 : psfpos, 2: psfsub
     13// 0 : backsub, 1 : lin_resid, 2: psfsub
    214
    315# if (HAVE_KAPA)
    … …  
    921static bool isVisual = false;
    1022static int kapa = -1;
     23static int kapa2 = -1;
     24static int kapa3 = -1;
    1125
    1226bool psphotSetVisual (bool mode) {
    1327
    14   isVisual = mode;
    15   return true;
    16 }
    17 
    18 bool psphotVisualShowImage (pmReadout *readout) {
    19 
    20   KiiImage image;
    21   KapaImageData data;
    22   Coords coords;
    23 
    24   if (!isVisual) return true;
    25 
    26   if (kapa == -1) {
    27     kapa = pmKapaOpen (true, "psphot");
     28    isVisual = mode;
     29    return true;
     30}
     31
     32bool psphotVisualScaleImage (int kapaFD, psImage *inImage, const char *name, int channel) {
     33
     34    KiiImage image;
     35    KapaImageData data;
     36    Coords coords;
     37
     38    strcpy (coords.ctype, "RA---TAN");
     39
     40    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     41    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
     42    if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) {
     43        fprintf (stderr, "failed to get background values\n");
     44        return false;
     45    }
     46
     47    image.data2d = inImage->data.F32;
     48    image.Nx = inImage->numCols;
     49    image.Ny = inImage->numRows;
     50
     51    strcpy (data.name, name);
     52    strcpy (data.file, name);
     53    data.zero = stats->robustMedian - stats->robustStdev;
     54    data.range = 5*stats->robustStdev;
     55    data.logflux = 0;
     56 
     57    KiiSetChannel (kapaFD, channel);
     58    KiiNewPicture2D (kapaFD, &image, &data, &coords);
     59
     60    psFree (stats);
     61    psFree (rng);
     62
     63    return true;
     64}
     65
     66bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max) {
     67
     68    KiiImage image;
     69    KapaImageData data;
     70    Coords coords;
     71
     72    strcpy (coords.ctype, "RA---TAN");
     73
     74    image.data2d = inImage->data.F32;
     75    image.Nx = inImage->numCols;
     76    image.Ny = inImage->numRows;
     77
     78    strcpy (data.name, name);
     79    strcpy (data.file, name);
     80    data.zero = min;
     81    data.range = max - min;
     82    data.logflux = 0;
     83 
     84    KiiSetChannel (kapaFD, channel);
     85    KiiNewPicture2D (kapaFD, &image, &data, &coords);
     86
     87    return true;
     88}
     89
     90bool psphotVisualShowImage (pmConfig *config, pmReadout *readout) {
     91
     92    if (!isVisual) return true;
     93
    2894    if (kapa == -1) {
    29       fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    30       isVisual = false;
    31       return false;
    32     }
    33   } 
    34 
    35   // XXX get image stats
    36   image.data2d = readout->image->data.F32;
    37   image.Nx = readout->image->numCols;
    38   image.Ny = readout->image->numRows;
    39 
    40   strcpy (data.name, "image");
    41   strcpy (data.file, "image");
    42   data.zero = 0.0;
    43   data.range = 1000.0;
    44   data.logflux = 0;
     95        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
     96        if (kapa == -1) {
     97            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     98            isVisual = false;
     99            return false;
     100        }
     101    } 
     102
     103    // psphotVisualShowMask (kapa, readout->mask, "mask", 2);
     104    psphotVisualScaleImage (kapa, readout->weight, "weight", 1);
     105    psphotVisualScaleImage (kapa, readout->image, "image", 0);
     106
     107    // pause and wait for user input:
     108    // continue, save (provide name), ??
     109    char key[10];
     110    fprintf (stdout, "[c]ontinue? ");
     111    if (!fgets(key, 8, stdin)) {
     112        psWarning("Unable to read option");
     113    }
     114    return true;
     115}
     116
     117bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout) {
     118
     119    pmReadout *backgnd;
     120
     121    if (!isVisual) return true;
     122
     123    if (kapa == -1) {
     124        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
     125        if (kapa == -1) {
     126            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     127            isVisual = false;
     128            return false;
     129        }
     130    } 
     131
     132    bool status = false;
     133    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND");
     134
     135    if (file->mode == PM_FPA_MODE_INTERNAL) {
     136        backgnd = file->readout;
     137    } else {
     138        backgnd = pmFPAviewThisReadout (view, file->fpa);
     139    }
     140
     141    psphotVisualScaleImage (kapa, backgnd->image, "backgnd", 2);
     142    psphotVisualScaleImage (kapa, readout->image, "backsub", 0);
     143
     144    // pause and wait for user input:
     145    // continue, save (provide name), ??
     146    char key[10];
     147    fprintf (stdout, "[c]ontinue? ");
     148    if (!fgets(key, 8, stdin)) {
     149        psWarning("Unable to read option");
     150    }
     151    return true;
     152}
     153
     154bool psphotVisualShowSignificance (psImage *image) {
     155
     156    if (!isVisual) return true;
     157
     158    if (kapa == -1) {
     159        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
     160        if (kapa == -1) {
     161            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     162            isVisual = false;
     163            return false;
     164        }
     165    } 
     166
     167    // XXX test: image->data.F32[10][10] = 10000;
     168    psphotVisualRangeImage (kapa, image, "signif", 2, -1.0, 25.0*25.0);
     169
     170    // pause and wait for user input:
     171    // continue, save (provide name), ??
     172    char key[10];
     173    fprintf (stdout, "[c]ontinue? ");
     174    if (!fgets(key, 8, stdin)) {
     175        psWarning("Unable to read option");
     176    }
     177    return true;
     178}
     179
     180bool psphotVisualShowPeaks (pmConfig *config, const pmFPAview *view, pmDetections *detections) {
     181
     182    int Noverlay;
     183    KiiOverlay *overlay;
    45184 
    46   KiiSetChannel (kapa, 0);
    47   KiiNewPicture2D (kapa, &image, &data, &coords);
     185    if (!isVisual) return true;
     186
     187    if (kapa == -1) {
     188        fprintf (stderr, "kapa not opened, skipping\n");
     189        return false;
     190    }
     191
     192    psArray *peaks = detections->peaks;
     193
     194    // note: this uses the Ohana allocation tools:
     195    // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1);
     196    ALLOCATE (overlay, KiiOverlay, peaks->n);
     197
     198    Noverlay = 0;
     199    for (int i = 0; i < peaks->n; i++) {
     200
     201        pmPeak *peak = peaks->data[i];
     202        if (peak == NULL) continue;
     203
     204        overlay[Noverlay].type = KII_OVERLAY_BOX;
     205        overlay[Noverlay].x = peak->xf;
     206        overlay[Noverlay].y = peak->yf;
     207        overlay[Noverlay].dx = 2.0;
     208        overlay[Noverlay].dy = 2.0;
     209        overlay[Noverlay].angle = 0.0;
     210        overlay[Noverlay].text = NULL;
     211        Noverlay ++;
     212
     213# if (0)
     214        overlay[Noverlay].type = KII_OVERLAY_BOX;
     215        overlay[Noverlay].x = peak->x;
     216        overlay[Noverlay].y = peak->y;
     217        overlay[Noverlay].dx = 1.0;
     218        overlay[Noverlay].dy = 1.0;
     219        overlay[Noverlay].angle = 0.0;
     220        overlay[Noverlay].text = NULL;
     221        Noverlay ++;
     222
     223        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     224        overlay[Noverlay].x = peak->xf;
     225        overlay[Noverlay].y = peak->yf;
     226        overlay[Noverlay].dx = 2.0;
     227        overlay[Noverlay].dy = 2.0;
     228        overlay[Noverlay].angle = 0.0;
     229        overlay[Noverlay].text = NULL;
     230        Noverlay ++;
     231# endif
     232    }
     233
     234# if (0)
     235    overlay[Noverlay].type = KII_OVERLAY_BOX;
     236    overlay[Noverlay].x = 10.0;
     237    overlay[Noverlay].y = 10.0;
     238    overlay[Noverlay].dx = 0.5;
     239    overlay[Noverlay].dy = 0.5;
     240    overlay[Noverlay].angle = 0.0;
     241    overlay[Noverlay].text = NULL;
     242    Noverlay ++;
     243# endif
     244
     245    KiiLoadOverlay (kapa, overlay, Noverlay, "red");
     246    FREE (overlay);
     247
     248    // pause and wait for user input:
     249    // continue, save (provide name), ??
     250    char key[10];
     251    fprintf (stdout, "[c]ontinue? ");
     252    if (!fgets(key, 8, stdin)) {
     253        psWarning("Unable to read option");
     254    }
     255    return true;
     256}
     257
     258// XXX this is just wrong: outlines entire image
     259bool psphotVisualShowFootprints (pmConfig *config, const pmFPAview *view, pmDetections *detections) {
     260
     261    int Noverlay;
     262    KiiOverlay *overlay;
    48263 
    49   // display the weight
    50   image.data2d = readout->weight->data.F32;
    51   image.Nx = readout->weight->numCols;
    52   image.Ny = readout->weight->numRows;
    53 
    54   strcpy (data.name, "weight");
    55   strcpy (data.file, "weight");
    56   data.zero = 0.0;
    57   data.range = 1000.0;
    58   data.logflux = 0;
     264    if (!isVisual) return true;
     265
     266    if (kapa == -1) {
     267        fprintf (stderr, "kapa not opened, skipping\n");
     268        return false;
     269    }
     270
     271    psArray *footprints = detections->footprints;
     272    if (!footprints) return true;
     273
     274    // note: this uses the Ohana allocation tools:
     275    int NOVERLAY = footprints->n;
     276    ALLOCATE (overlay, KiiOverlay, NOVERLAY);
     277
     278    Noverlay = 0;
     279    for (int i = 0; i < footprints->n; i++) {
     280
     281        pmSpan *span = NULL;
     282
     283        pmFootprint *footprint = footprints->data[i];
     284        if (footprint == NULL) continue;
     285        if (footprint->spans == NULL) continue;
     286        if (footprint->spans->n < 1) continue;
     287
     288        // draw the top
     289        span = footprint->spans->data[0];
     290        overlay[Noverlay].type = KII_OVERLAY_LINE;
     291        overlay[Noverlay].x = span->x0;
     292        overlay[Noverlay].y = span->y;
     293        overlay[Noverlay].dx = span->x1 - span->x0;
     294        overlay[Noverlay].dy = 0;
     295        overlay[Noverlay].angle = 0.0;
     296        overlay[Noverlay].text = NULL;
     297        Noverlay ++;
     298        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     299
     300        int ys = span->y;
     301        int x0s = span->x0;
     302        int x1s = span->x1;
     303
     304        // draw the outer span edges
     305        for (int j = 1; j < footprint->spans->n; j++) {
     306            pmSpan *span1 = footprint->spans->data[j];
     307
     308            int ye = span1->y;
     309            int x0e = span1->x0;
     310            int x1e = span1->x1;
     311
     312            // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right)
     313            // find all of the spans in this row and generate x0e, x01:
     314            for (int k = j + 1; k < footprint->spans->n; k++) {
     315                pmSpan *span2 = footprint->spans->data[k];
     316                if (span2->y > span1->y) break;
     317                x0e = PS_MIN (x0e, span2->x0);
     318                x1e = PS_MAX (x1e, span2->x1);
     319                j++;
     320            }
     321
     322            overlay[Noverlay].type = KII_OVERLAY_LINE;
     323            overlay[Noverlay].x = x0s;
     324            overlay[Noverlay].y = ys;
     325            overlay[Noverlay].dx = x0e - x0s;
     326            overlay[Noverlay].dy = ye - ys;
     327            overlay[Noverlay].angle = 0.0;
     328            overlay[Noverlay].text = NULL;
     329            Noverlay ++;
     330            CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     331
     332            overlay[Noverlay].type = KII_OVERLAY_LINE;
     333            overlay[Noverlay].x = x1s;
     334            overlay[Noverlay].y = ys;
     335            overlay[Noverlay].dx = x1e - x1s;
     336            overlay[Noverlay].dy = ye - ys;
     337            overlay[Noverlay].angle = 0.0;
     338            overlay[Noverlay].text = NULL;
     339            Noverlay ++;
     340            CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     341
     342            ys = ye;
     343            x0s = x0e;
     344            x1s = x1e;
     345        }
     346
     347        // draw the bottom
     348        span = footprint->spans->data[footprint->spans->n - 1];
     349        overlay[Noverlay].type = KII_OVERLAY_LINE;
     350        overlay[Noverlay].x = span->x0;
     351        overlay[Noverlay].y = span->y;
     352        overlay[Noverlay].dx = span->x1 - span->x0;
     353        overlay[Noverlay].dy = 0;
     354        overlay[Noverlay].angle = 0.0;
     355        overlay[Noverlay].text = NULL;
     356        Noverlay ++;
     357        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     358    }
     359
     360    KiiLoadOverlay (kapa, overlay, Noverlay, "blue");
     361    FREE (overlay);
     362
     363    // pause and wait for user input:
     364    // continue, save (provide name), ??
     365    char key[10];
     366    fprintf (stdout, "[c]ontinue? ");
     367    if (!fgets(key, 8, stdin)) {
     368        psWarning("Unable to read option");
     369    }
     370    return true;
     371}
     372
     373bool psphotVisualShowMoments (pmConfig *config, const pmFPAview *view, psArray *sources) {
     374
     375    int Noverlay;
     376    KiiOverlay *overlay;
    59377 
    60   KiiSetChannel (kapa, 0);
    61   KiiNewPicture2D (kapa, &image, &data, &coords);
     378    psEllipseMoments emoments;
     379    psEllipseAxes axes;
     380
     381    if (!isVisual) return true;
     382
     383    if (kapa == -1) {
     384        fprintf (stderr, "kapa not opened, skipping\n");
     385        return false;
     386    }
     387
     388    // note: this uses the Ohana allocation tools:
     389    ALLOCATE (overlay, KiiOverlay, sources->n);
     390
     391    Noverlay = 0;
     392    for (int i = 0; i < sources->n; i++) {
     393
     394        pmSource *source = sources->data[i];
     395        if (source == NULL) continue;
     396
     397        pmMoments *moments = source->moments;
     398        if (moments == NULL) continue;
     399
     400        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     401        overlay[Noverlay].x = moments->x;
     402        overlay[Noverlay].y = moments->y;
     403
     404        emoments.x2 = moments->Sx;
     405        emoments.y2 = moments->Sy;
     406        emoments.xy = moments->Sxy;
     407
     408        axes = psEllipseMomentsToAxes (emoments, 20.0);
     409
     410        overlay[Noverlay].dx = 2.0*axes.major;
     411        overlay[Noverlay].dy = 2.0*axes.minor;
     412        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
     413        overlay[Noverlay].text = NULL;
     414        Noverlay ++;
     415    }
     416
     417    KiiLoadOverlay (kapa, overlay, Noverlay, "yellow");
     418    FREE (overlay);
     419
     420    // pause and wait for user input:
     421    // continue, save (provide name), ??
     422    char key[10];
     423    fprintf (stdout, "[c]ontinue? ");
     424    if (!fgets(key, 8, stdin)) {
     425        psWarning("Unable to read option");
     426    }
     427
     428    return true;
     429}
     430
     431bool psphotVisualPlotMoments (pmConfig *config, const pmFPAview *view, psArray *sources) {
     432
     433    bool status;
     434    Graphdata graphdata;
     435
     436    if (!isVisual) return true;
     437
     438    if (kapa3 == -1) {
     439        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
     440        if (kapa3 == -1) {
     441            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     442            isVisual = false;
     443            return false;
     444        }
     445    } 
     446
     447    // select the current recipe
     448    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
     449    if (!recipe) {
     450        fprintf (stderr, "missing recipe, skipping moments plot\n");
     451        return false;
     452    }
     453
     454    KapaClearPlots (kapa3);
     455    KapaInitGraph (&graphdata);
     456
     457    float psfX  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X");
     458    float psfY  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.Y");
     459    float psfdX = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DX");
     460    float psfdY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY");
     461    float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
     462    float Rx = psfdX * PSF_CLUMP_NSIGMA;
     463    float Ry = psfdY * PSF_CLUMP_NSIGMA;
     464
     465    // examine sources to set data range
     466    graphdata.xmin = -0.05;
     467    graphdata.ymin = -0.05;
     468    graphdata.xmax = 2.0*psfX;
     469    graphdata.ymax = 2.0*psfY;
     470    KapaSetLimits (kapa3, &graphdata);
     471
     472    KapaSetFont (kapa3, "helvetica", 14);
     473    KapaBox (kapa3, &graphdata);
     474    KapaSendLabel (kapa3, "&ss&h_x| (pixels)", KAPA_LABEL_XM);
     475    KapaSendLabel (kapa3, "&ss&h_y| (pixels)", KAPA_LABEL_YM);
     476
     477    psVector *xBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     478    psVector *yBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     479    psVector *xFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     480    psVector *yFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     481
     482    // construct the vectors
     483    int nB = 0;
     484    int nF = 0;
     485    for (int i = 0; i < sources->n; i++) {
     486        pmSource *source = sources->data[i];
     487        if (source->moments == NULL)
     488            continue;
     489
     490        xFaint->data.F32[nF] = source->moments->Sx;
     491        yFaint->data.F32[nF] = source->moments->Sy;
     492        nF++;
     493
     494        // XXX make this a user-defined cutoff
     495        if (source->moments->SN < 50)
     496            continue;
     497
     498        xBright->data.F32[nB] = source->moments->Sx;
     499        yBright->data.F32[nB] = source->moments->Sy;
     500        nB++;
     501    }
     502    xFaint->n = nF;
     503    yFaint->n = nF;
     504
     505    xBright->n = nB;
     506    yBright->n = nB;
     507
     508    graphdata.color = KapaColorByName ("black");
     509    graphdata.ptype = 0;
     510    graphdata.size = 0.3;
     511    graphdata.style = 2;
     512    KapaPrepPlot (kapa3, nF, &graphdata);
     513    KapaPlotVector (kapa3, nF, xFaint->data.F32, "x");
     514    KapaPlotVector (kapa3, nF, yFaint->data.F32, "y");
     515
     516    graphdata.color = KapaColorByName ("red");
     517    graphdata.ptype = 2;
     518    graphdata.size = 0.5;
     519    graphdata.style = 2;
     520    KapaPrepPlot (kapa3, nB, &graphdata);
     521    KapaPlotVector (kapa3, nB, xBright->data.F32, "x");
     522    KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
     523
     524    // draw a circle centered on psfX,Y with size of the psf limit
     525    psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
     526    psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
     527    for (int i = 0; i < xLimit->n; i++) {
     528        xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
     529        yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
     530    }
     531    graphdata.color = KapaColorByName ("blue");
     532    graphdata.style = 0;
     533    KapaPrepPlot (kapa3, xLimit->n, &graphdata);
     534    KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x");
     535    KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y");
     536    psFree (xLimit);
     537    psFree (yLimit);
     538
     539# if (0)
     540    // *** make a histogram of the source counts in the x and y directions
     541    psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0);
     542    psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0);
     543    psVectorHistogram (nX, xFaint, NULL, NULL, 0);
     544    psVectorHistogram (nY, yFaint, NULL, NULL, 0);
     545    psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
     546    psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32);
     547    psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
     548    psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
     549    for (int i = 0; i < nX->nums->n; i++) {
     550        dX->data.F32[i] = nX->nums->data.S32[i];
     551        vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
     552    }
     553    for (int i = 0; i < nY->nums->n; i++) {
     554        dY->data.F32[i] = nY->nums->data.S32[i];
     555        vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
     556    }
     557
     558    graphdata.color = KapaColorByName ("black");
     559    graphdata.ptype = 0;
     560    graphdata.size = 0.0;
     561    graphdata.style = 0;
     562    KapaPrepPlot (kapa3, dX->n, &graphdata);
     563    KapaPlotVector (kapa3, dX->n, dX->data.F32, "x");
     564    KapaPlotVector (kapa3, vX->n, vX->data.F32, "y");
     565
     566    psFree (nX);
     567    psFree (dX);
     568    psFree (vX);
     569
     570    psFree (nY);
     571    psFree (dY);
     572    psFree (vY);
     573# endif
     574
     575    psFree (xBright);
     576    psFree (yBright);
     577    psFree (xFaint);
     578    psFree (yFaint);
     579
     580    // pause and wait for user input:
     581    // continue, save (provide name), ??
     582    char key[10];
     583    fprintf (stdout, "[c]ontinue? ");
     584    if (!fgets(key, 8, stdin)) {
     585        psWarning("Unable to read option");
     586    }
     587    return true;
     588}
     589
     590// assumes 'kapa' value is checked and set
     591bool psphotVisualShowRoughClass_Single (pmConfig *config, const pmFPAview *view, psArray *sources, pmSourceType type, pmSourceMode mode, char *color) {
     592
     593    int Noverlay;
     594    KiiOverlay *overlay;
    62595 
    63   // display the weight
    64   image.data2d = readout->weight->data.F32;
    65   image.Nx = readout->weight->numCols;
    66   image.Ny = readout->weight->numRows;
    67 
    68   strcpy (data.name, "mask");
    69   strcpy (data.file, "maak");
    70   data.zero = 0.0;
    71   data.range = 1000.0;
    72   data.logflux = 0;
    73  
    74   KiiSetChannel (kapa, 0);
    75   KiiNewPicture2D (kapa, &image, &data, &coords);
    76 }
    77 
     596    psEllipseMoments emoments;
     597    psEllipseAxes axes;
     598
     599    // note: this uses the Ohana allocation tools:
     600    ALLOCATE (overlay, KiiOverlay, sources->n);
     601
     602    Noverlay = 0;
     603    for (int i = 0; i < sources->n; i++) {
     604
     605        pmSource *source = sources->data[i];
     606        if (source == NULL) continue;
     607
     608        if (source->type != type) continue;
     609        if (mode && !(source->mode & mode)) continue;
     610
     611        pmMoments *moments = source->moments;
     612        if (moments == NULL) continue;
     613
     614        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     615        overlay[Noverlay].x = moments->x;
     616        overlay[Noverlay].y = moments->y;
     617
     618        emoments.x2 = moments->Sx;
     619        emoments.y2 = moments->Sy;
     620        emoments.xy = moments->Sxy;
     621
     622        axes = psEllipseMomentsToAxes (emoments, 20.0);
     623
     624        overlay[Noverlay].dx = 2.0*axes.major;
     625        overlay[Noverlay].dy = 2.0*axes.minor;
     626        overlay[Noverlay].angle = axes.theta * PS_DEG_RAD;
     627        overlay[Noverlay].text = NULL;
     628        Noverlay ++;
     629    }
     630
     631    KiiLoadOverlay (kapa, overlay, Noverlay, color);
     632    FREE (overlay);
     633
     634    return true;
     635}
     636
     637bool psphotVisualShowRoughClass (pmConfig *config, const pmFPAview *view, psArray *sources) {
     638
     639    psphotVisualPlotMoments (config, view, sources);
     640
     641    if (!isVisual) return true;
     642
     643    if (kapa == -1) {
     644        fprintf (stderr, "kapa not opened, skipping\n");
     645        return false;
     646    }
     647
     648    KiiEraseOverlay (kapa, "yellow"); // moments
     649
     650    psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_STAR, 0, "red");
     651    psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_EXTENDED, 0, "blue");
     652    psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_DEFECT, 0, "blue");
     653    psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_SATURATED, 0, "red");
     654    psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_PSFSTAR, "yellow");
     655    psphotVisualShowRoughClass_Single (config, view, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_SATSTAR, "green");
     656
     657    // pause and wait for user input:
     658    // continue, save (provide name), ??
     659    char key[10];
     660    fprintf (stdout, "red: STAR or SAT AREA; blue: EXTENDED or DEFECT; green: SATSTAR; yellow: PSFSTAR\n");
     661    fprintf (stdout, "[c]ontinue? ");
     662    if (!fgets(key, 8, stdin)) {
     663        psWarning("Unable to read option");
     664    }
     665    return true;
     666}
     667
     668bool psphotVisualShowPSFStars (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources) {
     669
     670    bool status;
     671
     672    if (!isVisual) return true;
     673
     674    if (kapa2 == -1) {
     675        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars");
     676        if (kapa2 == -1) {
     677            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     678            isVisual = false;
     679            return false;
     680        }
     681    } 
     682
     683    // select the current recipe
     684    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
     685    if (!recipe) {
     686        psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE);
     687        return false;
     688    }
     689
     690    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
     691    psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     692    assert (maskVal);
     693
     694    // the source images are written to an image 10x the size of a PSF object
     695    // float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
     696    // PS_ASSERT (status, false);
     697
     698    int DX = 21;
     699    int DY = 21;
     700
     701    // examine PSF sources in S/N order (brightest first)
     702    sources = psArraySort (sources, pmSourceSortBySN);
     703
     704    // counters to track the size of the image and area used in a row
     705    int dX = 0;                         // starting corner of next box
     706    int dY = 0;                         // height of row so far
     707    int NX = 20*DX;                     // full width of output image
     708    int NY = 0;                         // total height of output image
     709
     710    // first, examine the PSF stars:
     711    // - determine bounding boxes for summary image
     712    for (int i = 0; i < sources->n; i++) {
     713
     714        pmSource *source = sources->data[i];
     715
     716        bool keep = false;
     717        keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
     718        if (!keep) continue;
     719
     720        // how does this subimage get placed into the output image?
     721        // DX = source->pixels->numCols
     722        // DY = source->pixels->numRows
     723
     724        if (dX + DX > NX) {
     725            // too wide for the rest of this row
     726            if (dX == 0) {
     727                // alone on this row
     728                NY += DY;
     729                dX = 0;
     730                dY = 0;
     731            } else {
     732                // start the next row
     733                NY += dY;
     734                dX = DX;
     735                dY = DY;
     736            }
     737        } else {
     738            // extend this row
     739            dX += DX;
     740            dY = PS_MAX (dY, DY);
     741        }
     742    }
     743    NY += DY;
     744
     745    // allocate output image
     746    psImage *outpos = psImageAlloc (NX, NY, PS_TYPE_F32);
     747    psImage *outsub = psImageAlloc (NX, NY, PS_TYPE_F32);
     748
     749    int Xo = 0;                         // starting corner of next box
     750    int Yo = 0;                         // starting corner of next box
     751    dY = 0;                             // height of row so far
     752
     753    int nPSF = 0;
     754
     755    // next, examine the PSF stars:
     756    // - create output image array
     757    for (int i = 0; i < sources->n; i++) {
     758
     759        pmSource *source = sources->data[i];
     760
     761        bool keep = false;
     762        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
     763            nPSF ++;
     764            keep = true;
     765        }
     766        if (!keep) continue;
     767
     768        if (Xo + DX > NX) {
     769            // too wide for the rest of this row
     770            if (Xo == 0) {
     771                // place source alone on this row
     772                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     773                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
     774
     775                psphotSubWithTest (source, false, maskVal); // remove source (force)
     776                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
     777                psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded
     778
     779                Yo += DY;
     780                Xo = 0;
     781                dY = 0;
     782            } else {
     783                // start the next row
     784                Yo += dY;
     785                Xo = 0;
     786
     787                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     788                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
     789
     790                psphotSubWithTest (source, false, maskVal); // remove source (force)
     791                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
     792                psphotSetState (source, false, maskVal); // replace source (has been subtracted)
     793
     794                Xo = DX;
     795                dY = DY;
     796            }
     797        } else {
     798            // extend this row
     799            psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     800            psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
     801
     802            psphotSubWithTest (source, false, maskVal); // remove source (force)
     803            psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
     804            psphotSetState (source, false, maskVal); // replace source (has been subtracted)
     805
     806            Xo += DX;
     807            dY = PS_MAX (dY, DY);
     808        }
     809    }
     810
     811    psphotVisualScaleImage (kapa2, outpos, "psfpos", 0);
     812    psphotVisualScaleImage (kapa2, outsub, "psfsub", 1);
     813
     814    // pause and wait for user input:
     815    // continue, save (provide name), ??
     816    char key[10];
     817    fprintf (stdout, "[c]ontinue? ");
     818    if (!fgets(key, 8, stdin)) {
     819        psWarning("Unable to read option");
     820    }
     821
     822    psFree (outpos);
     823    psFree (outsub);
     824    return true;
     825}
     826
     827bool psphotVisualShowSatStars (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources) {
     828
     829    bool status;
     830
     831    if (!isVisual) return true;
     832
     833    if (kapa2 == -1) {
     834        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:images");
     835        if (kapa2 == -1) {
     836            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     837            isVisual = false;
     838            return false;
     839        }
     840    } 
     841
     842    // select the current recipe
     843    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
     844    if (!recipe) {
     845        psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE);
     846        return false;
     847    }
     848
     849    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
     850    psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     851    assert (maskVal);
     852
     853    // the source images are written to an image 10x the size of a PSF object
     854    // float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
     855    // PS_ASSERT (status, false);
     856
     857    int DX = 41;
     858    int DY = 41;
     859
     860    // examine PSF sources in S/N order (brightest first)
     861    sources = psArraySort (sources, pmSourceSortBySN);
     862
     863    // counters to track the size of the image and area used in a row
     864    int dX = 0;                         // starting corner of next box
     865    int dY = 0;                         // height of row so far
     866    int NX = 10*DX;                     // full width of output image
     867    int NY = 0;                         // total height of output image
     868
     869    // first, examine the PSF and SAT stars:
     870    // - determine bounding boxes for summary image
     871    for (int i = 0; i < sources->n; i++) {
     872
     873        pmSource *source = sources->data[i];
     874
     875        bool keep = false;
     876        keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
     877        if (!keep) continue;
     878
     879        // how does this subimage get placed into the output image?
     880        // DX = source->pixels->numCols
     881        // DY = source->pixels->numRows
     882
     883        if (dX + DX > NX) {
     884            // too wide for the rest of this row
     885            if (dX == 0) {
     886                // alone on this row
     887                NY += DY;
     888                dX = 0;
     889                dY = 0;
     890            } else {
     891                // start the next row
     892                NY += dY;
     893                dX = DX;
     894                dY = DY;
     895            }
     896        } else {
     897            // extend this row
     898            dX += DX;
     899            dY = PS_MAX (dY, DY);
     900        }
     901    }
     902    NY += DY;
     903
     904    // allocate output image
     905    psImage *outsat = psImageAlloc (NX, NY, PS_TYPE_F32);
     906
     907    int Xo = 0;                         // starting corner of next box
     908    int Yo = 0;                         // starting corner of next box
     909    dY = 0;                             // height of row so far
     910
     911    int nSAT = 0;
     912
     913    // next, examine the SAT stars:
     914    // - create output image array
     915    for (int i = 0; i < sources->n; i++) {
     916
     917        pmSource *source = sources->data[i];
     918
     919        bool keep = false;
     920        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     921            nSAT ++;
     922            keep = true;
     923        }           
     924        if (!keep) continue;
     925
     926        if (Xo + DX > NX) {
     927            // too wide for the rest of this row
     928            if (Xo == 0) {
     929                // place source alone on this row
     930                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     931                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
     932                psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded
     933
     934                Yo += DY;
     935                Xo = 0;
     936                dY = 0;
     937            } else {
     938                // start the next row
     939                Yo += dY;
     940                Xo = 0;
     941                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     942                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
     943                psphotSetState (source, true, maskVal); // replace source (has been subtracted)
     944
     945                Xo = DX;
     946                dY = DY;
     947            }
     948        } else {
     949            // extend this row
     950            psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     951            psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
     952            psphotSetState (source, true, maskVal); // replace source (has been subtracted)
     953
     954            Xo += DX;
     955            dY = PS_MAX (dY, DY);
     956        }
     957    }
     958
     959    psphotVisualScaleImage (kapa2, outsat, "satstar", 2);
     960
     961    // pause and wait for user input:
     962    // continue, save (provide name), ??
     963    char key[10];
     964    fprintf (stdout, "[c]ontinue? ");
     965    if (!fgets(key, 8, stdin)) {
     966        psWarning("Unable to read option");
     967    }
     968
     969    psFree (outsat);
     970    return true;
     971}
     972
     973bool psphotVisualShowPSFModel (pmConfig *config, pmReadout *readout, pmPSF *psf) {
     974
     975    if (!isVisual) return true;
     976
     977    if (kapa2 == -1) {
     978        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars");
     979        if (kapa2 == -1) {
     980            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     981            isVisual = false;
     982            return false;
     983        }
     984    } 
     985
     986    int DX = 64;
     987    int DY = 64;
     988
     989    psImage *psfMosaic = psImageAlloc (3*DX, 3*DY, PS_TYPE_F32);
     990    psImageInit (psfMosaic, 0.0);
     991   
     992    pmModel *modelRef = pmModelAlloc(psf->type);
     993
     994# if (1)
     995    // generate a fake model at each of the 3x3 image grid positions
     996    for (int x = -1; x <= +1; x ++) {
     997        for (int y = -1; y <= +1; y ++) {
     998            // use the center of the center pixel of the image
     999            float xc = (0.5 + 0.45*x)*readout->image->numCols + readout->image->col0;
     1000            float yc = (0.5 + 0.45*y)*readout->image->numRows + readout->image->row0;
     1001
     1002            // assign the x and y coords to the image center
     1003            // create an object with center intensity of 1000
     1004            modelRef->params->data.F32[PM_PAR_SKY] = 0;
     1005            modelRef->params->data.F32[PM_PAR_I0] = 1000;
     1006            modelRef->params->data.F32[PM_PAR_XPOS] = xc;
     1007            modelRef->params->data.F32[PM_PAR_YPOS] = yc;
     1008   
     1009            // create modelPSF from this model
     1010            pmModel *model = pmModelFromPSF (modelRef, psf);
     1011
     1012            // place the reference object in the image center
     1013            // no need to mask the source here
     1014            // XXX should we measure this for the analytical model only or the full model?
     1015            pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, -x*DX, -y*DY);
     1016        }
     1017    }
     1018# else
     1019    {
     1020        // use the center of the center pixel of the image
     1021        float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
     1022        float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
     1023
     1024        // assign the x and y coords to the image center
     1025        // create an object with center intensity of 1000
     1026        modelRef->params->data.F32[PM_PAR_SKY] = 0;
     1027        modelRef->params->data.F32[PM_PAR_I0] = 1000;
     1028        modelRef->params->data.F32[PM_PAR_XPOS] = xc;
     1029        modelRef->params->data.F32[PM_PAR_YPOS] = yc;
     1030   
     1031        // create modelPSF from this model
     1032        pmModel *model = pmModelFromPSF (modelRef, psf);
     1033
     1034        // place the reference object in the image center
     1035        // no need to mask the source here
     1036        // XXX should we measure this for the analytical model only or the full model?
     1037        pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, 32.0, 16.0);
     1038    }
     1039# endif
     1040
     1041    psImage *psfLogFlux = (psImage *) psUnaryOp (NULL, psfMosaic, "log");
     1042    psphotVisualRangeImage (kapa2, psfLogFlux, "psf_mosaic", 1, -2.0, 3.0);
     1043
     1044    // pause and wait for user input:
     1045    // continue, save (provide name), ??
     1046    char key[10];
     1047    fprintf (stdout, "[c]ontinue? ");
     1048    if (!fgets(key, 8, stdin)) {
     1049        psWarning("Unable to read option");
     1050    }
     1051    return true;
     1052}
     1053
     1054bool psphotVisualShowLinearFit (pmConfig *config, pmReadout *readout) {
     1055
     1056    if (!isVisual) return true;
     1057
     1058    if (kapa == -1) {
     1059        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
     1060        if (kapa == -1) {
     1061            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     1062            isVisual = false;
     1063            return false;
     1064        }
     1065    } 
     1066
     1067    psphotVisualScaleImage (kapa, readout->image, "lin_resid", 1);
     1068
     1069    // pause and wait for user input:
     1070    // continue, save (provide name), ??
     1071    char key[10];
     1072    fprintf (stdout, "[c]ontinue? ");
     1073    if (!fgets(key, 8, stdin)) {
     1074        psWarning("Unable to read option");
     1075    }
     1076    return true;
     1077}
     1078
     1079# else
     1080
     1081bool psphotSetVisual (bool mode){}
     1082bool psphotVisualShowImage (pmConfig *config, pmReadout *readout){}
     1083bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout){}
     1084bool psphotVisualShowSignificance (psImage *image){}
     1085bool psphotVisualShowPeaks (pmConfig *config, const pmFPAview *view, pmDetections *detections){}
     1086bool psphotVisualShowFootprints (pmConfig *config, const pmFPAview *view, pmDetections *detections){}
     1087bool psphotVisualShowMoments (pmConfig *config, const pmFPAview *view, psArray *sources){}
     1088bool psphotVisualShowRoughClass (pmConfig *config, const pmFPAview *view, psArray *sources){}
     1089bool psphotVisualShowPSF (pmConfig *config, const pmFPAview *view, pmPSF *psf, psArray *sources){}
     1090bool psphotVisualShowLinearFit (pmConfig *config, pmReadout *readout){}
     1091
     1092# endif
Note: See TracChangeset for help on using the changeset viewer.