IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 5, 2009, 5:03:33 PM (17 years ago)
Author:
Paul Price
Message:

Merging pap_branch_20090128. Resolved a small number of conflicts. Compiles, but not tested in detail.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotVisual.c

    r21183 r21366  
    22
    33// this function displays representative images as the psphot analysis progresses:
    4 // 0 : image, 1 : weight
    5 // 0 : backsub, 1 : weight, 2 : backgnd
    6 // 0 : backsub, 1 : weight, 2 : signif
     4// 0 : image, 1 : variance
     5// 0 : backsub, 1 : variance, 2 : backgnd
     6// 0 : backsub, 1 : variance, 2 : signif
    77// (overlay peaks on images)
    88// (overlay footprints on images)
    99// (overlay moments on images)
    10 // (overlay rough class on images) 
     10// (overlay rough class on images)
    1111// 0 : backsub, 1 : psfpos, 2: psfsub
    1212// 0 : backsub, 1 : lin_resid, 2: psfsub
    … …  
    4040    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
    4141    if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) {
    42         fprintf (stderr, "failed to get background values\n");
    43         return false;
     42        fprintf (stderr, "failed to get background values\n");
     43        return false;
    4444    }
    4545
    … …  
    4949    ALLOCATE (image.data2d, float *, image.Ny);
    5050    for (int iy = 0; iy < image.Ny; iy++) {
    51         ALLOCATE (image.data2d[iy], float, image.Nx);
    52         for (int ix = 0; ix < image.Nx; ix++) {
    53             image.data2d[iy][ix] = inImage->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix];
    54         }
     51        ALLOCATE (image.data2d[iy], float, image.Nx);
     52        for (int ix = 0; ix < image.Nx; ix++) {
     53            image.data2d[iy][ix] = inImage->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix];
     54        }
    5555    }
    5656
    … …  
    6060    data.range = 32;
    6161    data.logflux = 0;
    62  
     62
    6363    KiiSetChannel (kapaFD, channel);
    6464    KiiNewPicture2D (kapaFD, &image, &data, &coords);
    6565
    6666    for (int iy = 0; iy < image.Ny; iy++) {
    67         free (image.data2d[iy]);
     67        free (image.data2d[iy]);
    6868    }
    6969    free (image.data2d);
    … …  
    8686    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
    8787    if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) {
    88         fprintf (stderr, "failed to get background values\n");
    89         return false;
     88        fprintf (stderr, "failed to get background values\n");
     89        return false;
    9090    }
    9191
    … …  
    9999    data.range = 5*stats->robustStdev;
    100100    data.logflux = 0;
    101  
     101
    102102    KiiSetChannel (kapaFD, channel);
    103103    KiiNewPicture2D (kapaFD, &image, &data, &coords);
    … …  
    126126    data.range = max - min;
    127127    data.logflux = 0;
    128  
     128
    129129    KiiSetChannel (kapaFD, channel);
    130130    KiiNewPicture2D (kapaFD, &image, &data, &coords);
    … …  
    139139    if (kapa == -1) {
    140140        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
    141         if (kapa == -1) {
    142             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    143             isVisual = false;
    144             return false;
    145         }
    146     } 
     141        if (kapa == -1) {
     142            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     143            isVisual = false;
     144            return false;
     145        }
     146    }
    147147
    148148    // psphotVisualShowMask (kapa, readout->mask, "mask", 2);
    149     psphotVisualScaleImage (kapa, readout->weight, "weight", 1);
     149    psphotVisualScaleImage (kapa, readout->variance, "variance", 1);
    150150    psphotVisualScaleImage (kapa, readout->image, "image", 0);
    151151
    … …  
    155155    fprintf (stdout, "[c]ontinue? ");
    156156    if (!fgets(key, 8, stdin)) {
    157         psWarning("Unable to read option");
     157        psWarning("Unable to read option");
    158158    }
    159159    return true;
    … …  
    168168    if (kapa == -1) {
    169169        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
    170         if (kapa == -1) {
    171             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    172             isVisual = false;
    173             return false;
    174         }
    175     } 
     170        if (kapa == -1) {
     171            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     172            isVisual = false;
     173            return false;
     174        }
     175    }
    176176
    177177    bool status = false;
    … …  
    179179
    180180    if (file->mode == PM_FPA_MODE_INTERNAL) {
    181         backgnd = file->readout;
     181        backgnd = file->readout;
    182182    } else {
    183         backgnd = pmFPAviewThisReadout (view, file->fpa);
     183        backgnd = pmFPAviewThisReadout (view, file->fpa);
    184184    }
    185185
    … …  
    192192    fprintf (stdout, "[c]ontinue? ");
    193193    if (!fgets(key, 8, stdin)) {
    194         psWarning("Unable to read option");
     194        psWarning("Unable to read option");
    195195    }
    196196    return true;
    … …  
    203203    if (kapa == -1) {
    204204        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
    205         if (kapa == -1) {
    206             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    207             isVisual = false;
    208             return false;
    209         }
    210     } 
     205        if (kapa == -1) {
     206            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     207            isVisual = false;
     208            return false;
     209        }
     210    }
    211211
    212212    // XXX test: image->data.F32[10][10] = 10000;
    … …  
    218218    fprintf (stdout, "[c]ontinue? ");
    219219    if (!fgets(key, 8, stdin)) {
    220         psWarning("Unable to read option");
     220        psWarning("Unable to read option");
    221221    }
    222222    return true;
    … …  
    227227    int Noverlay;
    228228    KiiOverlay *overlay;
    229  
     229
    230230    if (!isVisual) return true;
    231231
    232232    if (kapa == -1) {
    233         fprintf (stderr, "kapa not opened, skipping\n");
    234         return false;
     233        fprintf (stderr, "kapa not opened, skipping\n");
     234        return false;
    235235    }
    236236
    … …  
    244244    for (int i = 0; i < peaks->n; i++) {
    245245
    246         pmPeak *peak = peaks->data[i];
    247         if (peak == NULL) continue;
    248 
    249         overlay[Noverlay].type = KII_OVERLAY_BOX;
    250         overlay[Noverlay].x = peak->xf;
    251         overlay[Noverlay].y = peak->yf;
    252         overlay[Noverlay].dx = 2.0;
    253         overlay[Noverlay].dy = 2.0;
    254         overlay[Noverlay].angle = 0.0;
    255         overlay[Noverlay].text = NULL;
    256         Noverlay ++;
     246        pmPeak *peak = peaks->data[i];
     247        if (peak == NULL) continue;
     248
     249        overlay[Noverlay].type = KII_OVERLAY_BOX;
     250        overlay[Noverlay].x = peak->xf;
     251        overlay[Noverlay].y = peak->yf;
     252        overlay[Noverlay].dx = 2.0;
     253        overlay[Noverlay].dy = 2.0;
     254        overlay[Noverlay].angle = 0.0;
     255        overlay[Noverlay].text = NULL;
     256        Noverlay ++;
    257257
    258258# if (0)
    259         overlay[Noverlay].type = KII_OVERLAY_BOX;
    260         overlay[Noverlay].x = peak->x;
    261         overlay[Noverlay].y = peak->y;
    262         overlay[Noverlay].dx = 1.0;
    263         overlay[Noverlay].dy = 1.0;
    264         overlay[Noverlay].angle = 0.0;
    265         overlay[Noverlay].text = NULL;
    266         Noverlay ++;
    267 
    268         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    269         overlay[Noverlay].x = peak->xf;
    270         overlay[Noverlay].y = peak->yf;
    271         overlay[Noverlay].dx = 2.0;
    272         overlay[Noverlay].dy = 2.0;
    273         overlay[Noverlay].angle = 0.0;
    274         overlay[Noverlay].text = NULL;
    275         Noverlay ++;
     259        overlay[Noverlay].type = KII_OVERLAY_BOX;
     260        overlay[Noverlay].x = peak->x;
     261        overlay[Noverlay].y = peak->y;
     262        overlay[Noverlay].dx = 1.0;
     263        overlay[Noverlay].dy = 1.0;
     264        overlay[Noverlay].angle = 0.0;
     265        overlay[Noverlay].text = NULL;
     266        Noverlay ++;
     267
     268        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     269        overlay[Noverlay].x = peak->xf;
     270        overlay[Noverlay].y = peak->yf;
     271        overlay[Noverlay].dx = 2.0;
     272        overlay[Noverlay].dy = 2.0;
     273        overlay[Noverlay].angle = 0.0;
     274        overlay[Noverlay].text = NULL;
     275        Noverlay ++;
    276276# endif
    277277    }
    … …  
    296296    fprintf (stdout, "[c]ontinue? ");
    297297    if (!fgets(key, 8, stdin)) {
    298         psWarning("Unable to read option");
     298        psWarning("Unable to read option");
    299299    }
    300300    return true;
    … …  
    305305    int Noverlay;
    306306    KiiOverlay *overlay;
    307  
     307
    308308    if (!isVisual) return true;
    309309
    310310    if (kapa == -1) {
    311         fprintf (stderr, "kapa not opened, skipping\n");
    312         return false;
     311        fprintf (stderr, "kapa not opened, skipping\n");
     312        return false;
    313313    }
    314314
    … …  
    323323    for (int i = 0; i < footprints->n; i++) {
    324324
    325         pmSpan *span = NULL;
    326 
    327         pmFootprint *footprint = footprints->data[i];
    328         if (footprint == NULL) continue;
    329         if (footprint->spans == NULL) continue;
    330         if (footprint->spans->n < 1) continue;
    331 
    332         // draw the top
    333         span = footprint->spans->data[0];
    334         overlay[Noverlay].type = KII_OVERLAY_LINE;
    335         overlay[Noverlay].x = span->x0;
    336         overlay[Noverlay].y = span->y;
    337         overlay[Noverlay].dx = span->x1 - span->x0;
    338         overlay[Noverlay].dy = 0;
    339         overlay[Noverlay].angle = 0.0;
    340         overlay[Noverlay].text = NULL;
    341         Noverlay ++;
    342         CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    343 
    344         int ys = span->y;
    345         int x0s = span->x0;
    346         int x1s = span->x1;
    347 
    348         // draw the outer span edges
    349         for (int j = 1; j < footprint->spans->n; j++) {
    350             pmSpan *span1 = footprint->spans->data[j];
    351 
    352             int ye = span1->y;
    353             int x0e = span1->x0;
    354             int x1e = span1->x1;
    355 
    356             // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right)
    357             // find all of the spans in this row and generate x0e, x01:
    358             for (int k = j + 1; k < footprint->spans->n; k++) {
    359                 pmSpan *span2 = footprint->spans->data[k];
    360                 if (span2->y > span1->y) break;
    361                 x0e = PS_MIN (x0e, span2->x0);
    362                 x1e = PS_MAX (x1e, span2->x1);
    363                 j++;
    364             }
    365 
    366             overlay[Noverlay].type = KII_OVERLAY_LINE;
    367             overlay[Noverlay].x = x0s;
    368             overlay[Noverlay].y = ys;
    369             overlay[Noverlay].dx = x0e - x0s;
    370             overlay[Noverlay].dy = ye - ys;
    371             overlay[Noverlay].angle = 0.0;
    372             overlay[Noverlay].text = NULL;
    373             Noverlay ++;
    374             CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    375 
    376             overlay[Noverlay].type = KII_OVERLAY_LINE;
    377             overlay[Noverlay].x = x1s;
    378             overlay[Noverlay].y = ys;
    379             overlay[Noverlay].dx = x1e - x1s;
    380             overlay[Noverlay].dy = ye - ys;
    381             overlay[Noverlay].angle = 0.0;
    382             overlay[Noverlay].text = NULL;
    383             Noverlay ++;
    384             CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    385 
    386             ys = ye;
    387             x0s = x0e;
    388             x1s = x1e;
    389         }
    390 
    391         // draw the bottom
    392         span = footprint->spans->data[footprint->spans->n - 1];
    393         overlay[Noverlay].type = KII_OVERLAY_LINE;
    394         overlay[Noverlay].x = span->x0;
    395         overlay[Noverlay].y = span->y;
    396         overlay[Noverlay].dx = span->x1 - span->x0;
    397         overlay[Noverlay].dy = 0;
    398         overlay[Noverlay].angle = 0.0;
    399         overlay[Noverlay].text = NULL;
    400         Noverlay ++;
    401         CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     325        pmSpan *span = NULL;
     326
     327        pmFootprint *footprint = footprints->data[i];
     328        if (footprint == NULL) continue;
     329        if (footprint->spans == NULL) continue;
     330        if (footprint->spans->n < 1) continue;
     331
     332        // draw the top
     333        span = footprint->spans->data[0];
     334        overlay[Noverlay].type = KII_OVERLAY_LINE;
     335        overlay[Noverlay].x = span->x0;
     336        overlay[Noverlay].y = span->y;
     337        overlay[Noverlay].dx = span->x1 - span->x0;
     338        overlay[Noverlay].dy = 0;
     339        overlay[Noverlay].angle = 0.0;
     340        overlay[Noverlay].text = NULL;
     341        Noverlay ++;
     342        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     343
     344        int ys = span->y;
     345        int x0s = span->x0;
     346        int x1s = span->x1;
     347
     348        // draw the outer span edges
     349        for (int j = 1; j < footprint->spans->n; j++) {
     350            pmSpan *span1 = footprint->spans->data[j];
     351
     352            int ye = span1->y;
     353            int x0e = span1->x0;
     354            int x1e = span1->x1;
     355
     356            // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right)
     357            // find all of the spans in this row and generate x0e, x01:
     358            for (int k = j + 1; k < footprint->spans->n; k++) {
     359                pmSpan *span2 = footprint->spans->data[k];
     360                if (span2->y > span1->y) break;
     361                x0e = PS_MIN (x0e, span2->x0);
     362                x1e = PS_MAX (x1e, span2->x1);
     363                j++;
     364            }
     365
     366            overlay[Noverlay].type = KII_OVERLAY_LINE;
     367            overlay[Noverlay].x = x0s;
     368            overlay[Noverlay].y = ys;
     369            overlay[Noverlay].dx = x0e - x0s;
     370            overlay[Noverlay].dy = ye - ys;
     371            overlay[Noverlay].angle = 0.0;
     372            overlay[Noverlay].text = NULL;
     373            Noverlay ++;
     374            CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     375
     376            overlay[Noverlay].type = KII_OVERLAY_LINE;
     377            overlay[Noverlay].x = x1s;
     378            overlay[Noverlay].y = ys;
     379            overlay[Noverlay].dx = x1e - x1s;
     380            overlay[Noverlay].dy = ye - ys;
     381            overlay[Noverlay].angle = 0.0;
     382            overlay[Noverlay].text = NULL;
     383            Noverlay ++;
     384            CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     385
     386            ys = ye;
     387            x0s = x0e;
     388            x1s = x1e;
     389        }
     390
     391        // draw the bottom
     392        span = footprint->spans->data[footprint->spans->n - 1];
     393        overlay[Noverlay].type = KII_OVERLAY_LINE;
     394        overlay[Noverlay].x = span->x0;
     395        overlay[Noverlay].y = span->y;
     396        overlay[Noverlay].dx = span->x1 - span->x0;
     397        overlay[Noverlay].dy = 0;
     398        overlay[Noverlay].angle = 0.0;
     399        overlay[Noverlay].text = NULL;
     400        Noverlay ++;
     401        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    402402    }
    403403
    … …  
    410410    fprintf (stdout, "[c]ontinue? ");
    411411    if (!fgets(key, 8, stdin)) {
    412         psWarning("Unable to read option");
     412        psWarning("Unable to read option");
    413413    }
    414414    return true;
    … …  
    419419    int Noverlay;
    420420    KiiOverlay *overlay;
    421  
     421
    422422    psEllipseMoments emoments;
    423423    psEllipseAxes axes;
    … …  
    426426
    427427    if (kapa == -1) {
    428         fprintf (stderr, "kapa not opened, skipping\n");
    429         return false;
     428        fprintf (stderr, "kapa not opened, skipping\n");
     429        return false;
    430430    }
    431431
    … …  
    436436    for (int i = 0; i < sources->n; i++) {
    437437
    438         pmSource *source = sources->data[i];
    439         if (source == NULL) continue;
    440 
    441         pmMoments *moments = source->moments;
    442         if (moments == NULL) continue;
    443 
    444         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    445         overlay[Noverlay].x = moments->Mx;
    446         overlay[Noverlay].y = moments->My;
    447 
    448         emoments.x2 = moments->Mxx;
    449         emoments.xy = moments->Mxy;
    450         emoments.y2 = moments->Myy;
    451 
    452         axes = psEllipseMomentsToAxes (emoments, 20.0);
    453 
    454         overlay[Noverlay].dx = 2.0*axes.major;
    455         overlay[Noverlay].dy = 2.0*axes.minor;
    456         overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;  // XXXXXXXX the axes angle is negative to display of object on kapa
    457         overlay[Noverlay].text = NULL;
    458         Noverlay ++;
     438        pmSource *source = sources->data[i];
     439        if (source == NULL) continue;
     440
     441        pmMoments *moments = source->moments;
     442        if (moments == NULL) continue;
     443
     444        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     445        overlay[Noverlay].x = moments->Mx;
     446        overlay[Noverlay].y = moments->My;
     447
     448        emoments.x2 = moments->Mxx;
     449        emoments.xy = moments->Mxy;
     450        emoments.y2 = moments->Myy;
     451
     452        axes = psEllipseMomentsToAxes (emoments, 20.0);
     453
     454        overlay[Noverlay].dx = 2.0*axes.major;
     455        overlay[Noverlay].dy = 2.0*axes.minor;
     456        overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;  // XXXXXXXX the axes angle is negative to display of object on kapa
     457        overlay[Noverlay].text = NULL;
     458        Noverlay ++;
    459459    }
    460460
    … …  
    467467    fprintf (stdout, "[c]ontinue? ");
    468468    if (!fgets(key, 8, stdin)) {
    469         psWarning("Unable to read option");
     469        psWarning("Unable to read option");
    470470    }
    471471
    … …  
    482482    if (kapa3 == -1) {
    483483        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    484         if (kapa3 == -1) {
    485             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    486             isVisual = false;
    487             return false;
    488         }
    489     } 
     484        if (kapa3 == -1) {
     485            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     486            isVisual = false;
     487            return false;
     488        }
     489    }
    490490
    491491    KapaClearPlots (kapa3);
    … …  
    494494    // there are N regions: use the first (guaranteed to exist) to get the overall limits
    495495    psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
    496    
     496
    497497    float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    498498    float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
    … …  
    558558
    559559    // XXX draw N circles to outline the clumps
    560     { 
    561         // draw a circle centered on psfX,Y with size of the psf limit
    562         psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
    563         psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
    564 
    565         int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
    566         float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
    567 
    568         graphdata.color = KapaColorByName ("blue");
    569         graphdata.style = 0;
    570 
    571         for (int n = 0; n < nRegions; n++) {
    572 
    573             char regionName[64];
    574             snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
    575             psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
    576    
    577             float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    578             float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
    579             float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
    580             float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
    581             float Rx = psfdX * PSF_CLUMP_NSIGMA;
    582             float Ry = psfdY * PSF_CLUMP_NSIGMA;
    583 
    584             for (int i = 0; i < xLimit->n; i++) {
    585                 xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
    586                 yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
    587             }
    588             KapaPrepPlot (kapa3, xLimit->n, &graphdata);
    589             KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x");
    590             KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y");
    591         }
    592         psFree (xLimit);
    593         psFree (yLimit);
     560    {
     561        // draw a circle centered on psfX,Y with size of the psf limit
     562        psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
     563        psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
     564
     565        int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     566        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
     567
     568        graphdata.color = KapaColorByName ("blue");
     569        graphdata.style = 0;
     570
     571        for (int n = 0; n < nRegions; n++) {
     572
     573            char regionName[64];
     574            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
     575            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     576
     577            float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     578            float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
     579            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
     580            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
     581            float Rx = psfdX * PSF_CLUMP_NSIGMA;
     582            float Ry = psfdY * PSF_CLUMP_NSIGMA;
     583
     584            for (int i = 0; i < xLimit->n; i++) {
     585                xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
     586                yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
     587            }
     588            KapaPrepPlot (kapa3, xLimit->n, &graphdata);
     589            KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x");
     590            KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y");
     591        }
     592        psFree (xLimit);
     593        psFree (yLimit);
    594594    }
    595595
    596596# if (0)
    597     // *** make a histogram of the source counts in the x and y directions 
     597    // *** make a histogram of the source counts in the x and y directions
    598598    psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0);
    599599    psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0);
    … …  
    605605    psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
    606606    for (int i = 0; i < nX->nums->n; i++) {
    607         dX->data.F32[i] = nX->nums->data.S32[i];
    608         vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
     607        dX->data.F32[i] = nX->nums->data.S32[i];
     608        vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
    609609    }
    610610    for (int i = 0; i < nY->nums->n; i++) {
    611         dY->data.F32[i] = nY->nums->data.S32[i];
    612         vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
     611        dY->data.F32[i] = nY->nums->data.S32[i];
     612        vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
    613613    }
    614614
    … …  
    640640    fprintf (stdout, "[c]ontinue? ");
    641641    if (!fgets(key, 8, stdin)) {
    642         psWarning("Unable to read option");
     642        psWarning("Unable to read option");
    643643    }
    644644    return true;
    … …  
    650650    int Noverlay;
    651651    KiiOverlay *overlay;
    652  
     652
    653653    psEllipseMoments emoments;
    654654    psEllipseAxes axes;
    … …  
    660660    for (int i = 0; i < sources->n; i++) {
    661661
    662         pmSource *source = sources->data[i];
    663         if (source == NULL) continue;
    664 
    665         if (source->type != type) continue;
    666         if (mode && !(source->mode & mode)) continue;
    667 
    668         pmMoments *moments = source->moments;
    669         if (moments == NULL) continue;
    670 
    671         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    672         overlay[Noverlay].x = moments->Mx;
    673         overlay[Noverlay].y = moments->My;
    674 
    675         emoments.x2 = moments->Mxx;
    676         emoments.y2 = moments->Myy;
    677         emoments.xy = moments->Mxy;
    678 
    679         axes = psEllipseMomentsToAxes (emoments, 20.0);
    680 
    681         overlay[Noverlay].dx = 2.0*axes.major;
    682         overlay[Noverlay].dy = 2.0*axes.minor;
    683         overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;
    684         overlay[Noverlay].text = NULL;
    685         Noverlay ++;
     662        pmSource *source = sources->data[i];
     663        if (source == NULL) continue;
     664
     665        if (source->type != type) continue;
     666        if (mode && !(source->mode & mode)) continue;
     667
     668        pmMoments *moments = source->moments;
     669        if (moments == NULL) continue;
     670
     671        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     672        overlay[Noverlay].x = moments->Mx;
     673        overlay[Noverlay].y = moments->My;
     674
     675        emoments.x2 = moments->Mxx;
     676        emoments.y2 = moments->Myy;
     677        emoments.xy = moments->Mxy;
     678
     679        axes = psEllipseMomentsToAxes (emoments, 20.0);
     680
     681        overlay[Noverlay].dx = 2.0*axes.major;
     682        overlay[Noverlay].dy = 2.0*axes.minor;
     683        overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;
     684        overlay[Noverlay].text = NULL;
     685        Noverlay ++;
    686686    }
    687687
    … …  
    697697
    698698    if (kapa == -1) {
    699         fprintf (stderr, "kapa not opened, skipping\n");
    700         return false;
     699        fprintf (stderr, "kapa not opened, skipping\n");
     700        return false;
    701701    }
    702702
    … …  
    716716    fprintf (stdout, "[c]ontinue? ");
    717717    if (!fgets(key, 8, stdin)) {
    718         psWarning("Unable to read option");
     718        psWarning("Unable to read option");
    719719    }
    720720    return true;
    … …  
    727727    if (kapa2 == -1) {
    728728        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars");
    729         if (kapa2 == -1) {
    730             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    731             isVisual = false;
    732             return false;
    733         }
    734     } 
     729        if (kapa2 == -1) {
     730            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     731            isVisual = false;
     732            return false;
     733        }
     734    }
    735735
    736736    int DX = 64;
    … …  
    739739    psImage *psfMosaic = psImageAlloc (5*DX, 5*DY, PS_TYPE_F32);
    740740    psImageInit (psfMosaic, 0.0);
    741    
     741
    742742    psImage *funMosaic = psImageAlloc (5*DX, 5*DY, PS_TYPE_F32);
    743743    psImageInit (funMosaic, 0.0);
    … …  
    750750    // generate a fake model at each of the 3x3 image grid positions
    751751    for (int x = -2; x <= +2; x ++) {
    752         for (int y = -2; y <= +2; y ++) {
    753             // use the center of the center pixel of the image
    754             float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5;
    755             float yc = (int)((0.5 + 0.225*y)*readout->image->numRows) + readout->image->row0 + 0.5;
    756 
    757             // assign the x and y coords to the image center
    758             // create an object with center intensity of 1000
    759             modelRef->params->data.F32[PM_PAR_SKY] = 0;
    760             modelRef->params->data.F32[PM_PAR_I0] = 1000;
    761             modelRef->params->data.F32[PM_PAR_XPOS] = xc;
    762             modelRef->params->data.F32[PM_PAR_YPOS] = yc;
    763    
    764             // create modelPSF from this model
    765             pmModel *model = pmModelFromPSF (modelRef, psf);
    766             if (!model) continue;
    767 
    768             // place the reference object in the image center
    769             // no need to mask the source here
    770             // XXX should we measure this for the analytical model only or the full model?
    771             pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
    772             pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
    773             pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
    774             psFree (model);
    775         }
     752        for (int y = -2; y <= +2; y ++) {
     753            // use the center of the center pixel of the image
     754            float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5;
     755            float yc = (int)((0.5 + 0.225*y)*readout->image->numRows) + readout->image->row0 + 0.5;
     756
     757            // assign the x and y coords to the image center
     758            // create an object with center intensity of 1000
     759            modelRef->params->data.F32[PM_PAR_SKY] = 0;
     760            modelRef->params->data.F32[PM_PAR_I0] = 1000;
     761            modelRef->params->data.F32[PM_PAR_XPOS] = xc;
     762            modelRef->params->data.F32[PM_PAR_YPOS] = yc;
     763
     764            // create modelPSF from this model
     765            pmModel *model = pmModelFromPSF (modelRef, psf);
     766            if (!model) continue;
     767
     768            // place the reference object in the image center
     769            // no need to mask the source here
     770            // XXX should we measure this for the analytical model only or the full model?
     771            pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
     772            pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
     773            pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
     774            psFree (model);
     775        }
    776776    }
    777777
    … …  
    792792    fprintf (stdout, "[c]ontinue? ");
    793793    if (!fgets(key, 8, stdin)) {
    794         psWarning("Unable to read option");
     794        psWarning("Unable to read option");
    795795    }
    796796    return true;
    … …  
    805805    if (kapa2 == -1) {
    806806        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars");
    807         if (kapa2 == -1) {
    808             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    809             isVisual = false;
    810             return false;
    811         }
    812     } 
     807        if (kapa2 == -1) {
     808            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     809            isVisual = false;
     810            return false;
     811        }
     812    }
    813813
    814814    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    827827
    828828    // counters to track the size of the image and area used in a row
    829     int dX = 0;                         // starting corner of next box
    830     int dY = 0;                         // height of row so far
    831     int NX = 20*DX;                     // full width of output image
    832     int NY = 0;                         // total height of output image
     829    int dX = 0;                         // starting corner of next box
     830    int dY = 0;                         // height of row so far
     831    int NX = 20*DX;                     // full width of output image
     832    int NY = 0;                         // total height of output image
    833833
    834834    // first, examine the PSF stars:
    … …  
    838838        pmSource *source = sources->data[i];
    839839
    840         bool keep = false;
     840        bool keep = false;
    841841        keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
    842         if (!keep) continue;
    843 
    844         // how does this subimage get placed into the output image?
    845         // DX = source->pixels->numCols
    846         // DY = source->pixels->numRows
    847 
    848         if (dX + DX > NX) {
    849             // too wide for the rest of this row
    850             if (dX == 0) {
    851                 // alone on this row
    852                 NY += DY;
    853                 dX = 0;
    854                 dY = 0;
    855             } else {
    856                 // start the next row
    857                 NY += dY;
    858                 dX = DX;
    859                 dY = DY;
    860             }
    861         } else {
    862             // extend this row
    863             dX += DX;
    864             dY = PS_MAX (dY, DY);
    865         }
     842        if (!keep) continue;
     843
     844        // how does this subimage get placed into the output image?
     845        // DX = source->pixels->numCols
     846        // DY = source->pixels->numRows
     847
     848        if (dX + DX > NX) {
     849            // too wide for the rest of this row
     850            if (dX == 0) {
     851                // alone on this row
     852                NY += DY;
     853                dX = 0;
     854                dY = 0;
     855            } else {
     856                // start the next row
     857                NY += dY;
     858                dX = DX;
     859                dY = DY;
     860            }
     861        } else {
     862            // extend this row
     863            dX += DX;
     864            dY = PS_MAX (dY, DY);
     865        }
    866866    }
    867867    NY += DY;
    … …  
    873873    psImageInit (outsub, 0.0);
    874874
    875     int Xo = 0;                         // starting corner of next box
    876     int Yo = 0;                         // starting corner of next box
    877     dY = 0;                             // height of row so far
     875    int Xo = 0;                         // starting corner of next box
     876    int Yo = 0;                         // starting corner of next box
     877    dY = 0;                             // height of row so far
    878878
    879879    int nPSF = 0;
    … …  
    885885        pmSource *source = sources->data[i];
    886886
    887         bool keep = false;
     887        bool keep = false;
    888888        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
    889             nPSF ++;
    890             keep = true;
    891         }
    892         if (!keep) continue;
    893 
    894         if (Xo + DX > NX) {
    895             // too wide for the rest of this row
    896             if (Xo == 0) {
    897                 // place source alone on this row
    898                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    899                 psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    900 
    901                 psphotSubWithTest (source, false, maskVal); // remove source (force)
    902                 psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    903                 psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded
    904 
    905                 Yo += DY;
    906                 Xo = 0;
    907                 dY = 0;
    908             } else {
    909                 // start the next row
    910                 Yo += dY;
    911                 Xo = 0;
    912 
    913                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    914                 psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    915 
    916                 psphotSubWithTest (source, false, maskVal); // remove source (force)
    917                 psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    918                 psphotSetState (source, false, maskVal); // replace source (has been subtracted)
    919 
    920                 Xo = DX;
    921                 dY = DY;
    922             }
    923         } else {
    924             // extend this row
    925             psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    926             psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    927 
    928             psphotSubWithTest (source, false, maskVal); // remove source (force)
    929             psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    930             psphotSetState (source, false, maskVal); // replace source (has been subtracted)
    931 
    932             Xo += DX;
    933             dY = PS_MAX (dY, DY);
    934         }
     889            nPSF ++;
     890            keep = true;
     891        }
     892        if (!keep) continue;
     893
     894        if (Xo + DX > NX) {
     895            // too wide for the rest of this row
     896            if (Xo == 0) {
     897                // place source alone on this row
     898                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     899                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
     900
     901                psphotSubWithTest (source, false, maskVal); // remove source (force)
     902                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
     903                psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded
     904
     905                Yo += DY;
     906                Xo = 0;
     907                dY = 0;
     908            } else {
     909                // start the next row
     910                Yo += dY;
     911                Xo = 0;
     912
     913                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     914                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
     915
     916                psphotSubWithTest (source, false, maskVal); // remove source (force)
     917                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
     918                psphotSetState (source, false, maskVal); // replace source (has been subtracted)
     919
     920                Xo = DX;
     921                dY = DY;
     922            }
     923        } else {
     924            // extend this row
     925            psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     926            psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
     927
     928            psphotSubWithTest (source, false, maskVal); // remove source (force)
     929            psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
     930            psphotSetState (source, false, maskVal); // replace source (has been subtracted)
     931
     932            Xo += DX;
     933            dY = PS_MAX (dY, DY);
     934        }
    935935    }
    936936
    … …  
    943943    fprintf (stdout, "[c]ontinue? ");
    944944    if (!fgets(key, 8, stdin)) {
    945         psWarning("Unable to read option");
     945        psWarning("Unable to read option");
    946946    }
    947947
    … …  
    965965    if (kapa2 == -1) {
    966966        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:images");
    967         if (kapa2 == -1) {
    968             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    969             isVisual = false;
    970             return false;
    971         }
    972     } 
     967        if (kapa2 == -1) {
     968            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     969            isVisual = false;
     970            return false;
     971        }
     972    }
    973973
    974974    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    987987
    988988    // counters to track the size of the image and area used in a row
    989     int dX = 0;                         // starting corner of next box
    990     int dY = 0;                         // height of row so far
    991     int NX = 10*DX;                     // full width of output image
    992     int NY = 0;                         // total height of output image
     989    int dX = 0;                         // starting corner of next box
     990    int dY = 0;                         // height of row so far
     991    int NX = 10*DX;                     // full width of output image
     992    int NY = 0;                         // total height of output image
    993993
    994994    // first, examine the PSF and SAT stars:
    … …  
    998998        pmSource *source = sources->data[i];
    999999
    1000         bool keep = false;
     1000        bool keep = false;
    10011001        keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
    1002         if (!keep) continue;
    1003 
    1004         // how does this subimage get placed into the output image?
    1005         // DX = source->pixels->numCols
    1006         // DY = source->pixels->numRows
    1007 
    1008         if (dX + DX > NX) {
    1009             // too wide for the rest of this row
    1010             if (dX == 0) {
    1011                 // alone on this row
    1012                 NY += DY;
    1013                 dX = 0;
    1014                 dY = 0;
    1015             } else {
    1016                 // start the next row
    1017                 NY += dY;
    1018                 dX = DX;
    1019                 dY = DY;
    1020             }
    1021         } else {
    1022             // extend this row
    1023             dX += DX;
    1024             dY = PS_MAX (dY, DY);
    1025         }
     1002        if (!keep) continue;
     1003
     1004        // how does this subimage get placed into the output image?
     1005        // DX = source->pixels->numCols
     1006        // DY = source->pixels->numRows
     1007
     1008        if (dX + DX > NX) {
     1009            // too wide for the rest of this row
     1010            if (dX == 0) {
     1011                // alone on this row
     1012                NY += DY;
     1013                dX = 0;
     1014                dY = 0;
     1015            } else {
     1016                // start the next row
     1017                NY += dY;
     1018                dX = DX;
     1019                dY = DY;
     1020            }
     1021        } else {
     1022            // extend this row
     1023            dX += DX;
     1024            dY = PS_MAX (dY, DY);
     1025        }
    10261026    }
    10271027    NY += DY;
    … …  
    10311031    psImageInit (outsat, 0.0);
    10321032
    1033     int Xo = 0;                         // starting corner of next box
    1034     int Yo = 0;                         // starting corner of next box
    1035     dY = 0;                             // height of row so far
     1033    int Xo = 0;                         // starting corner of next box
     1034    int Yo = 0;                         // starting corner of next box
     1035    dY = 0;                             // height of row so far
    10361036
    10371037    int nSAT = 0;
    … …  
    10431043        pmSource *source = sources->data[i];
    10441044
    1045         bool keep = false;
     1045        bool keep = false;
    10461046        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    1047             nSAT ++;
    1048             keep = true;
    1049         }           
    1050         if (!keep) continue;
    1051 
    1052         if (Xo + DX > NX) {
    1053             // too wide for the rest of this row
    1054             if (Xo == 0) {
    1055                 // place source alone on this row
    1056                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    1057                 psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1058                 psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded
    1059 
    1060                 Yo += DY;
    1061                 Xo = 0;
    1062                 dY = 0;
    1063             } else {
    1064                 // start the next row
    1065                 Yo += dY;
    1066                 Xo = 0;
    1067                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    1068                 psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1069                 psphotSetState (source, true, maskVal); // replace source (has been subtracted)
    1070 
    1071                 Xo = DX;
    1072                 dY = DY;
    1073             }
    1074         } else {
    1075             // extend this row
    1076             psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    1077             psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1078             psphotSetState (source, true, maskVal); // replace source (has been subtracted)
    1079 
    1080             Xo += DX;
    1081             dY = PS_MAX (dY, DY);
    1082         }
     1047            nSAT ++;
     1048            keep = true;
     1049        }
     1050        if (!keep) continue;
     1051
     1052        if (Xo + DX > NX) {
     1053            // too wide for the rest of this row
     1054            if (Xo == 0) {
     1055                // place source alone on this row
     1056                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1057                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
     1058                psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded
     1059
     1060                Yo += DY;
     1061                Xo = 0;
     1062                dY = 0;
     1063            } else {
     1064                // start the next row
     1065                Yo += dY;
     1066                Xo = 0;
     1067                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1068                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
     1069                psphotSetState (source, true, maskVal); // replace source (has been subtracted)
     1070
     1071                Xo = DX;
     1072                dY = DY;
     1073            }
     1074        } else {
     1075            // extend this row
     1076            psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1077            psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
     1078            psphotSetState (source, true, maskVal); // replace source (has been subtracted)
     1079
     1080            Xo += DX;
     1081            dY = PS_MAX (dY, DY);
     1082        }
    10831083    }
    10841084
    … …  
    10901090    fprintf (stdout, "[c]ontinue? ");
    10911091    if (!fgets(key, 8, stdin)) {
    1092         psWarning("Unable to read option");
     1092        psWarning("Unable to read option");
    10931093    }
    10941094
    … …  
    11171117    float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
    11181118    for (int iy = 0; iy < source->pixels->numRows; iy++) {
    1119         for (int ix = 0; ix < source->pixels->numCols; ix++) {
    1120             if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
    1121                 // rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
    1122                 rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
    1123                 Rb->data.F32[nb] = log10(rb->data.F32[nb]);
    1124                 fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]);
    1125                 nb++;
    1126             } else {
    1127                 // rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
    1128                 rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
    1129                 Rg->data.F32[ng] = log10(rg->data.F32[ng]);
    1130                 fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]);
    1131                 ng++;
    1132             }
    1133         }
    1134     }
    1135  
     1119        for (int ix = 0; ix < source->pixels->numCols; ix++) {
     1120            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
     1121                // rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
     1122                rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
     1123                Rb->data.F32[nb] = log10(rb->data.F32[nb]);
     1124                fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]);
     1125                nb++;
     1126            } else {
     1127                // rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
     1128                rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
     1129                Rg->data.F32[ng] = log10(rg->data.F32[ng]);
     1130                fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]);
     1131                ng++;
     1132            }
     1133        }
     1134    }
     1135
    11361136    // reset source Add/Sub state to recorded
    1137     psphotSetState (source, state, maskVal); 
     1137    psphotSetState (source, state, maskVal);
    11381138
    11391139    KapaInitGraph (&graphdata);
    … …  
    11481148    graphdata.ymax = +5.05;
    11491149    KapaSetLimits (myKapa, &graphdata);
    1150  
     1150
    11511151    KapaSetFont (myKapa, "helvetica", 14);
    11521152    KapaBox (myKapa, &graphdata);
    11531153    KapaSendLabel (myKapa, "radius (pixels)", KAPA_LABEL_XM);
    11541154    KapaSendLabel (myKapa, "log flux (counts)", KAPA_LABEL_YM);
    1155                
     1155
    11561156    graphdata.color = KapaColorByName ("black");
    11571157    graphdata.ptype = 2;
    … …  
    11691169    KapaPlotVector (myKapa, nb, rb->data.F32, "x");
    11701170    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
    1171  
     1171
    11721172    // ** loglog **
    11731173    KapaSelectSection (myKapa, "loglog");
    … …  
    11791179    graphdata.ymax = +5.05;
    11801180    KapaSetLimits (myKapa, &graphdata);
    1181  
     1181
    11821182    KapaSetFont (myKapa, "helvetica", 14);
    11831183    KapaBox (myKapa, &graphdata);
    11841184    KapaSendLabel (myKapa, "log radius (pixels)", KAPA_LABEL_XM);
    11851185    KapaSendLabel (myKapa, "log flux (counts)", KAPA_LABEL_YM);
    1186                
     1186
    11871187    graphdata.color = KapaColorByName ("black");
    11881188    graphdata.ptype = 2;
    … …  
    12121212bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources) {
    12131213
    1214     KapaSection section;  // put the positive profile in one and the residuals in another? 
     1214    KapaSection section;  // put the positive profile in one and the residuals in another?
    12151215
    12161216    if (!isVisual) return true;
    … …  
    12181218    if (kapa3 == -1) {
    12191219        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    1220         if (kapa3 == -1) {
    1221             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1222             isVisual = false;
    1223             return false;
    1224         }
    1225     } 
     1220        if (kapa3 == -1) {
     1221            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     1222            isVisual = false;
     1223            return false;
     1224        }
     1225    }
    12261226
    12271227    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    12571257        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
    12581258
    1259         psphotVisualPlotRadialProfile (kapa3, source, maskVal);
    1260 
    1261         // pause and wait for user input:
    1262         // continue, save (provide name), ??
    1263         char key[10];
    1264         fprintf (stdout, "[e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : ");
    1265         if (!fgets(key, 8, stdin)) {
    1266             psWarning("Unable to read option");
    1267         }
    1268         if (key[0] == 'e') {
    1269             KapaClearPlots (kapa3);
    1270         }
    1271         if (key[0] == 's') {
    1272             break;
    1273         }
    1274     }       
     1259        psphotVisualPlotRadialProfile (kapa3, source, maskVal);
     1260
     1261        // pause and wait for user input:
     1262        // continue, save (provide name), ??
     1263        char key[10];
     1264        fprintf (stdout, "[e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : ");
     1265        if (!fgets(key, 8, stdin)) {
     1266            psWarning("Unable to read option");
     1267        }
     1268        if (key[0] == 'e') {
     1269            KapaClearPlots (kapa3);
     1270        }
     1271        if (key[0] == 's') {
     1272            break;
     1273        }
     1274    }
    12751275
    12761276    return true;
    … …  
    12821282    int NoverlayO, NOVERLAYO;
    12831283    KiiOverlay *overlayE, *overlayO;
    1284  
     1284
    12851285    psEllipseMoments emoments;
    12861286    psEllipseAxes axes;
    … …  
    12891289
    12901290    if (kapa == -1) {
    1291         fprintf (stderr, "kapa not opened, skipping\n");
    1292         return false;
     1291        fprintf (stderr, "kapa not opened, skipping\n");
     1292        return false;
    12931293    }
    12941294
    … …  
    13041304    for (int i = 0; i < sources->n; i++) {
    13051305
    1306         float Xo, Yo, Rmaj, Rmin, cs, sn;
    1307 
    1308         pmSource *source = sources->data[i];
    1309         if (source == NULL) continue;
    1310 
    1311         pmMoments *moments = source->moments;
    1312         if (0) {
    1313             emoments.x2 = moments->Mxx;
    1314             emoments.y2 = moments->Myy;
    1315             emoments.xy = moments->Mxy;
    1316             Xo = moments->Mx;
    1317             Yo = moments->My;
    1318 
    1319             axes = psEllipseMomentsToAxes (emoments, 20.0);
    1320             Rmaj = 2.0*axes.major;
    1321             Rmin = 2.0*axes.minor;
    1322             cs = cos(axes.theta);
    1323             sn = sin(axes.theta);
    1324         } else {
    1325             Rmaj = Rmin = 5.0;
    1326             cs = 1.0;
    1327             sn = 0.0;
    1328             Xo = source->peak->xf;
    1329             Yo = source->peak->yf;
    1330         }
    1331 
    1332         unsigned short int flagMask = 0x01;
    1333         for (int j = 0; j < 8; j++) {
    1334             if (source->mode & flagMask) {
    1335                 overlayE[NoverlayE].type = KII_OVERLAY_LINE;
    1336                 overlayE[NoverlayE].x = Xo;
    1337                 overlayE[NoverlayE].y = Yo;
    1338 
    1339                 float phi = j*M_PI/4.0;
    1340                 overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
    1341                 overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
    1342                 overlayE[NoverlayE].angle = 0;
    1343                 overlayE[NoverlayE].text = NULL;
    1344                 NoverlayE ++;
    1345                 CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100);
    1346             }
    1347             flagMask <<= 1;
    1348 
    1349             if (source->mode & flagMask) {
    1350                 overlayO[NoverlayO].type = KII_OVERLAY_LINE;
    1351                 overlayO[NoverlayO].x = Xo + 1;
    1352                 overlayO[NoverlayO].y = Yo;
    1353 
    1354                 float phi = j*M_PI/4.0;
    1355                 overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
    1356                 overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
    1357                 overlayO[NoverlayO].angle = 0;
    1358                 overlayO[NoverlayO].text = NULL;
    1359                 NoverlayO ++;
    1360                 CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100);
    1361             }
    1362             flagMask <<= 1;
    1363         }
     1306        float Xo, Yo, Rmaj, Rmin, cs, sn;
     1307
     1308        pmSource *source = sources->data[i];
     1309        if (source == NULL) continue;
     1310
     1311        pmMoments *moments = source->moments;
     1312        if (0) {
     1313            emoments.x2 = moments->Mxx;
     1314            emoments.y2 = moments->Myy;
     1315            emoments.xy = moments->Mxy;
     1316            Xo = moments->Mx;
     1317            Yo = moments->My;
     1318
     1319            axes = psEllipseMomentsToAxes (emoments, 20.0);
     1320            Rmaj = 2.0*axes.major;
     1321            Rmin = 2.0*axes.minor;
     1322            cs = cos(axes.theta);
     1323            sn = sin(axes.theta);
     1324        } else {
     1325            Rmaj = Rmin = 5.0;
     1326            cs = 1.0;
     1327            sn = 0.0;
     1328            Xo = source->peak->xf;
     1329            Yo = source->peak->yf;
     1330        }
     1331
     1332        unsigned short int flagMask = 0x01;
     1333        for (int j = 0; j < 8; j++) {
     1334            if (source->mode & flagMask) {
     1335                overlayE[NoverlayE].type = KII_OVERLAY_LINE;
     1336                overlayE[NoverlayE].x = Xo;
     1337                overlayE[NoverlayE].y = Yo;
     1338
     1339                float phi = j*M_PI/4.0;
     1340                overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
     1341                overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
     1342                overlayE[NoverlayE].angle = 0;
     1343                overlayE[NoverlayE].text = NULL;
     1344                NoverlayE ++;
     1345                CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100);
     1346            }
     1347            flagMask <<= 1;
     1348
     1349            if (source->mode & flagMask) {
     1350                overlayO[NoverlayO].type = KII_OVERLAY_LINE;
     1351                overlayO[NoverlayO].x = Xo + 1;
     1352                overlayO[NoverlayO].y = Yo;
     1353
     1354                float phi = j*M_PI/4.0;
     1355                overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
     1356                overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
     1357                overlayO[NoverlayO].angle = 0;
     1358                overlayO[NoverlayO].text = NULL;
     1359                NoverlayO ++;
     1360                CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100);
     1361            }
     1362            flagMask <<= 1;
     1363        }
    13641364    }
    13651365
    … …  
    13761376    fprintf (stdout, "[c]ontinue? ");
    13771377    if (!fgets(key, 8, stdin)) {
    1378         psWarning("Unable to read option");
     1378        psWarning("Unable to read option");
    13791379    }
    13801380
    … …  
    13901390
    13911391    if (kapa == -1) {
    1392         fprintf (stderr, "kapa not opened, skipping\n");
    1393         return false;
     1392        fprintf (stderr, "kapa not opened, skipping\n");
     1393        return false;
    13941394    }
    13951395
    … …  
    14021402    for (int i = 0; i < sources->n; i++) {
    14031403
    1404         pmSource *source = sources->data[i];
    1405         if (source == NULL) continue;
    1406 
    1407         if (!(source->mode & PM_SOURCE_MODE_CR_LIMIT)) continue;
    1408 
    1409         overlay[Noverlay].type = KII_OVERLAY_BOX;
    1410         overlay[Noverlay].x = source->peak->xf;
    1411         overlay[Noverlay].y = source->peak->yf;
    1412 
    1413         overlay[Noverlay].dx = 4;
    1414         overlay[Noverlay].dy = 4;
    1415         overlay[Noverlay].angle = 0;
    1416         overlay[Noverlay].text = NULL;
    1417         Noverlay ++;
    1418         CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     1404        pmSource *source = sources->data[i];
     1405        if (source == NULL) continue;
     1406
     1407        if (!(source->mode & PM_SOURCE_MODE_CR_LIMIT)) continue;
     1408
     1409        overlay[Noverlay].type = KII_OVERLAY_BOX;
     1410        overlay[Noverlay].x = source->peak->xf;
     1411        overlay[Noverlay].y = source->peak->yf;
     1412
     1413        overlay[Noverlay].dx = 4;
     1414        overlay[Noverlay].dy = 4;
     1415        overlay[Noverlay].angle = 0;
     1416        overlay[Noverlay].text = NULL;
     1417        Noverlay ++;
     1418        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    14191419    }
    14201420    KiiLoadOverlay (kapa, overlay, Noverlay, "red");
    … …  
    14241424    for (int i = 0; i < sources->n; i++) {
    14251425
    1426         pmSource *source = sources->data[i];
    1427         if (source == NULL) continue;
    1428 
    1429         // mark EXTs with yellow circles
    1430         if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
    1431 
    1432         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    1433         overlay[Noverlay].x = source->peak->xf;
    1434         overlay[Noverlay].y = source->peak->yf;
    1435 
    1436         overlay[Noverlay].dx = 10;
    1437         overlay[Noverlay].dy = 10;
    1438         overlay[Noverlay].angle = 0;
    1439         overlay[Noverlay].text = NULL;
    1440         Noverlay ++;
    1441         CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     1426        pmSource *source = sources->data[i];
     1427        if (source == NULL) continue;
     1428
     1429        // mark EXTs with yellow circles
     1430        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
     1431
     1432        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     1433        overlay[Noverlay].x = source->peak->xf;
     1434        overlay[Noverlay].y = source->peak->yf;
     1435
     1436        overlay[Noverlay].dx = 10;
     1437        overlay[Noverlay].dy = 10;
     1438        overlay[Noverlay].angle = 0;
     1439        overlay[Noverlay].text = NULL;
     1440        Noverlay ++;
     1441        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    14421442    }
    14431443
    … …  
    14531453    fprintf (stdout, "[c]ontinue? ");
    14541454    if (!fgets(key, 8, stdin)) {
    1455         psWarning("Unable to read option");
     1455        psWarning("Unable to read option");
    14561456    }
    14571457
    … …  
    14681468    if (kapa3 == -1) {
    14691469        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    1470         if (kapa3 == -1) {
    1471             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1472             isVisual = false;
    1473             return false;
    1474         }
    1475     } 
     1470        if (kapa3 == -1) {
     1471            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     1472            isVisual = false;
     1473            return false;
     1474        }
     1475    }
    14761476
    14771477    KapaClearPlots (kapa3);
    … …  
    15001500    for (int i = 0; i < sources->n; i++) {
    15011501        pmSource *source = sources->data[i];
    1502         if (!source) continue;
     1502        if (!source) continue;
    15031503        if (source->type != PM_SOURCE_TYPE_STAR) continue;
    1504         if (!isfinite (source->crNsigma)) continue;
     1504        if (!isfinite (source->crNsigma)) continue;
    15051505
    15061506        x->data.F32[n] = -2.5*log10(source->peak->flux);
    … …  
    15551555    for (int i = 0; i < sources->n; i++) {
    15561556        pmSource *source = sources->data[i];
    1557         if (!source) continue;
     1557        if (!source) continue;
    15581558        if (source->type != PM_SOURCE_TYPE_STAR) continue;
    1559         if (!isfinite (source->extNsigma)) continue;
     1559        if (!isfinite (source->extNsigma)) continue;
    15601560
    15611561        x->data.F32[n] = -2.5*log10(source->peak->flux);
    … …  
    15971597    fprintf (stdout, "[c]ontinue? ");
    15981598    if (!fgets(key, 8, stdin)) {
    1599         psWarning("Unable to read option");
     1599        psWarning("Unable to read option");
    16001600    }
    16011601    return true;
    … …  
    16081608    if (kapa == -1) {
    16091609        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
    1610         if (kapa == -1) {
    1611             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1612             isVisual = false;
    1613             return false;
    1614         }
    1615     } 
     1610        if (kapa == -1) {
     1611            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     1612            isVisual = false;
     1613            return false;
     1614        }
     1615    }
    16161616
    16171617    psphotVisualScaleImage (kapa, readout->image, "resid", 1);
    … …  
    16221622    fprintf (stdout, "[c]ontinue? ");
    16231623    if (!fgets(key, 8, stdin)) {
    1624         psWarning("Unable to read option");
     1624        psWarning("Unable to read option");
    16251625    }
    16261626    return true;
    … …  
    16351635    if (kapa3 == -1) {
    16361636        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    1637         if (kapa3 == -1) {
    1638             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1639             isVisual = false;
    1640             return false;
    1641         }
    1642     } 
     1637        if (kapa3 == -1) {
     1638            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     1639            isVisual = false;
     1640            return false;
     1641        }
     1642    }
    16431643
    16441644    KapaClearPlots (kapa3);
    … …  
    16571657    for (int i = 0; i < sources->n; i++) {
    16581658        pmSource *source = sources->data[i];
    1659         if (!source) continue;
     1659        if (!source) continue;
    16601660        if (source->type != PM_SOURCE_TYPE_STAR) continue;
    1661         if (!isfinite (source->apMag)) continue;
    1662         if (!isfinite (source->psfMag)) continue;
     1661        if (!isfinite (source->apMag)) continue;
     1662        if (!isfinite (source->psfMag)) continue;
    16631663
    16641664        x->data.F32[n] = source->psfMag;
    … …  
    17051705    fprintf (stdout, "[c]ontinue? ");
    17061706    if (!fgets(key, 8, stdin)) {
    1707         psWarning("Unable to read option");
     1707        psWarning("Unable to read option");
    17081708    }
    17091709    return true;
Note: See TracChangeset for help on using the changeset viewer.