Changeset 34404 for trunk/psphot/src/psphotVisual.c

Timestamp:

Sep 5, 2012, 4:19:30 PM (14 years ago)

Author:

eugene

Message:

adding SATSTAR profile code (but keeping it disabled for now)

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphotVisual.c (modified) (14 diffs)

trunk/psphot/src/psphotVisual.c

@@ -262 +262 @@
 }
 
+static psImage *posImage = NULL;
+static psImage *delImage = NULL;
+
 bool psphotVisualShowImage (pmReadout *readout) {
 
@@ -277 +280 @@
     psphotVisualScaleImage (kapa, readout->variance, readout->mask, "variance", 1.0, 1);
     psphotVisualScaleImage (kapa, readout->image, readout->mask, "image", sqrt(factor), 0);
+
+    if (posImage == NULL) {
+        posImage = psImageCopy (NULL, readout->image, PS_TYPE_F32);
+    }
 
     pmVisualAskUser(NULL);
@@ -1139 +1146 @@
     // after displaying (as an image) the psf stars, we cycle throught them and display their
     // radial profiles:
-    psphotVisualPlotRadialProfiles (recipe, sources);
+    psphotVisualPlotRadialProfiles (recipe, sources, PM_SOURCE_MODE_PSFSTAR);
 
     return true;
@@ -1273 +1280 @@
 }
 
-static void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1) 
+void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1) 
 {
     float x[2], y[2];
@@ -1285 +1292 @@
 }
 
-bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal) {
+bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal, pmSourceMode showmode) {
 
     Graphdata graphdata;
+
+    float Rmax = (showmode & PM_SOURCE_MODE_SATSTAR) ? 100.0 : 30.0;
 
     bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
@@ -1302 +1311 @@
     int ng = 0;
     int nb = 0;
-    float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
-    float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
+
+    float Xo = NAN;
+    float Yo = NAN;
+
+    if (source->modelPSF) {
+        Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
+        Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
+    } else {
+        Xo = source->moments->Mx - source->pixels->col0;
+        Yo = source->moments->My - source->pixels->row0;
+    }
+
     for (int iy = 0; iy < source->pixels->numRows; iy++) {
         for (int ix = 0; ix < source->pixels->numCols; ix++) {
-            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
+            if ((source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal)) {
                 rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
                 // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
@@ -1322 +1341 @@
     }
 
-    // generate model profiles (major and minor axis):
-    // create a model with theta = 0.0 so major and minor axes are equiv to x and y:
-    psEllipseShape rawShape, rotShape;
-
-    rawShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
-    rawShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
-    rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0);
-
-    axes.theta = 0.0;
-
-    rotShape = psEllipseAxesToShape (axes);
-
-    psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32);
-    for (int i = 0; i < source->modelPSF->params->n; i++) {
-        params->data.F32[i] = source->modelPSF->params->data.F32[i];
-    }
-    params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2;
-    params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2;
-    params->data.F32[PM_PAR_SXY] = rotShape.sxy;
-    params->data.F32[PM_PAR_XPOS] = 0.0;
-    params->data.F32[PM_PAR_YPOS] = 0.0;
-
-    psVector *rmod = psVectorAlloc(300, PS_TYPE_F32);
-    psVector *fmaj = psVectorAlloc(300, PS_TYPE_F32);
-    psVector *fmin = psVectorAlloc(300, PS_TYPE_F32);
-
-    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
-
-    float r = 0.0;
-    for (int i = 0; i < rmod->n; i++) {
-        r = i*0.1;
-        rmod->data.F32[i] = r;
-
-        coord->data.F32[1] = r;
-        coord->data.F32[0] = 0.0;
-        fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
-
-        coord->data.F32[0] = r;
-        coord->data.F32[1] = 0.0;
-        fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
-    }
-    psFree (coord);
-    psFree (params);
-
-    float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
-    float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
-    if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR); 
-
-    psEllipseMoments emoments;
-    emoments.x2 = source->moments->Mxx;
-    emoments.xy = source->moments->Mxy;
-    emoments.y2 = source->moments->Myy;
-    axes = psEllipseMomentsToAxes (emoments, 20.0);
-    float MOMENTS_MAJOR = 2.355*axes.major;
-    float MOMENTS_MINOR = 2.355*axes.minor;
-
-    float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
-
-    // reset source Add/Sub state to recorded
-    if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-
     KapaInitGraph (&graphdata);
 
@@ -1392 +1348 @@
     // examine sources to set data range
     graphdata.xmin =  -0.05;
-    graphdata.xmax = +30.05;
+    graphdata.xmax = Rmax + 0.05;
     graphdata.ymin = -0.05;
     graphdata.ymax = +8.05;
@@ -1418 +1374 @@
     KapaPlotVector (myKapa, nb, fb->data.F32, "y");
 
-    graphdata.color = KapaColorByName ("blue");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
-    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
-    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
-    plotline (myKapa, &graphdata, 0.0, logHM, 30.0, logHM);
-    plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, 0.0, 0.5*FWHM_MINOR, 5.0);
-    graphdata.ltype = 1;
-    plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, 0.0, 0.5*MOMENTS_MINOR, 5.0);
-    graphdata.ltype = 0;
-        
-    graphdata.color = KapaColorByName ("green");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
-    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
-    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
-    plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, 0.0, 0.5*FWHM_MAJOR, 5.0);
-    graphdata.ltype = 1;
-    plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, 0.0, 0.5*MOMENTS_MAJOR, 5.0);
-    graphdata.ltype = 0;
-        
-    for (int i = 0; i < rmod->n; i++) {
-        rmod->data.F32[i] = log10(rmod->data.F32[i]);
-    }
-
     // ** loglog **
     KapaSelectSection (myKapa, "loglog");
@@ -1452 +1379 @@
     // examine sources to set data range
     graphdata.xmin = -1.51;
-    graphdata.xmax = +1.51;
+    graphdata.xmax = log10(Rmax) + 0.02;
     graphdata.ymin = -0.05;
     graphdata.ymax = +8.05;
@@ -1479 +1406 @@
     KapaPlotVector (myKapa, nb, fb->data.F32, "y");
 
-    graphdata.color = KapaColorByName ("blue");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
-    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
-    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
-
-    graphdata.color = KapaColorByName ("green");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
-    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
-    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
-
-    psFree (rmod);
-    psFree (fmin);
-    psFree (fmaj);
+    if (source->modelPSF) {
+        // generate model profiles (major and minor axis):
+        // create a model with theta = 0.0 so major and minor axes are equiv to x and y:
+        psEllipseShape rawShape, rotShape;
+
+        rawShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
+        rawShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
+        rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
+
+        psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0);
+
+        axes.theta = 0.0;
+
+        rotShape = psEllipseAxesToShape (axes);
+
+        psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32);
+        for (int i = 0; i < source->modelPSF->params->n; i++) {
+            params->data.F32[i] = source->modelPSF->params->data.F32[i];
+        }
+        params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2;
+        params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2;
+        params->data.F32[PM_PAR_SXY] = rotShape.sxy;
+        params->data.F32[PM_PAR_XPOS] = 0.0;
+        params->data.F32[PM_PAR_YPOS] = 0.0;
+
+        psVector *rmod = psVectorAlloc(Rmax*10, PS_TYPE_F32);
+        psVector *fmaj = psVectorAlloc(Rmax*10, PS_TYPE_F32);
+        psVector *fmin = psVectorAlloc(Rmax*10, PS_TYPE_F32);
+
+        psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+
+        float r = 0.0;
+        for (int i = 0; i < rmod->n; i++) {
+            r = i*0.1;
+            rmod->data.F32[i] = r;
+
+            coord->data.F32[1] = r;
+            coord->data.F32[0] = 0.0;
+            fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
+
+            coord->data.F32[0] = r;
+            coord->data.F32[1] = 0.0;
+            fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
+        }
+        psFree (coord);
+        psFree (params);
+
+        float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
+        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+        if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR); 
+
+        psEllipseMoments emoments;
+        emoments.x2 = source->moments->Mxx;
+        emoments.xy = source->moments->Mxy;
+        emoments.y2 = source->moments->Myy;
+        axes = psEllipseMomentsToAxes (emoments, 20.0);
+        float MOMENTS_MAJOR = 2.355*axes.major;
+        float MOMENTS_MINOR = 2.355*axes.minor;
+
+        float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
+
+        // reset source Add/Sub state to recorded
+        if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+
+        // ** linlog **
+        KapaSelectSection (myKapa, "linlog");
+        KapaGetGraphData (myKapa, &graphdata);
+        graphdata.color = KapaColorByName ("blue");
+        graphdata.ptype = 0;
+        graphdata.size = 0.0;
+        graphdata.style = 0;
+        KapaPrepPlot   (myKapa, rmod->n, &graphdata);
+        KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
+        KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
+        plotline (myKapa, &graphdata, graphdata.xmin, logHM, graphdata.xmax, logHM);
+        plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, graphdata.ymin, 0.5*FWHM_MINOR, graphdata.ymax);
+        graphdata.ltype = 1;
+        plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, graphdata.ymin, 0.5*MOMENTS_MINOR, graphdata.ymax);
+        graphdata.ltype = 0;
+        
+        graphdata.color = KapaColorByName ("green");
+        graphdata.ptype = 0;
+        graphdata.size = 0.0;
+        graphdata.style = 0;
+        KapaPrepPlot   (myKapa, rmod->n, &graphdata);
+        KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
+        KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
+        plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, graphdata.ymin, 0.5*FWHM_MAJOR, graphdata.ymax);
+        graphdata.ltype = 1;
+        plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, graphdata.ymin, 0.5*MOMENTS_MAJOR, graphdata.ymax);
+        graphdata.ltype = 0;
+        
+        for (int i = 0; i < rmod->n; i++) {
+            rmod->data.F32[i] = log10(rmod->data.F32[i]);
+        }
+
+        // ** loglog **
+        KapaSelectSection (myKapa, "loglog");
+        KapaGetGraphData (myKapa, &graphdata);
+        graphdata.color = KapaColorByName ("blue");
+        graphdata.ptype = 0;
+        graphdata.size = 0.0;
+        graphdata.style = 0;
+        KapaPrepPlot   (myKapa, rmod->n, &graphdata);
+        KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
+        KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
+
+        graphdata.color = KapaColorByName ("green");
+        graphdata.ptype = 0;
+        graphdata.size = 0.0;
+        graphdata.style = 0;
+        KapaPrepPlot   (myKapa, rmod->n, &graphdata);
+        KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
+        KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
+
+        psFree (rmod);
+        psFree (fmin);
+        psFree (fmaj);
+    }
 
     psFree (rg);
@@ -1508 +1535 @@
 }
 
-bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources) {
+bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources, pmSourceMode showmode) {
 
     KapaSection section;  // put the positive profile in one and the residuals in another?
@@ -1549 +1576 @@
 
         pmSource *source = sources->data[i];
-        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
-
-        psphotVisualPlotRadialProfile (myKapa, source, maskVal);
+        if (!(source->mode & showmode)) continue;
+
+        psphotVisualPlotRadialProfile (myKapa, source, maskVal, showmode);
+
+        if (pmVisualTestLevel("psphot.image", 1)) {
+          int display = psphotKapaChannel (1);
+          KiiCenter (display, source->peak->xf, source->peak->yf, 1);
+          KiiOverlay overlay;
+          overlay.x = source->peak->xf;
+          overlay.y = source->peak->yf;
+          overlay.dx = 5;
+          overlay.dy = 5;
+          overlay.angle = 0.0;
+          overlay.type = KiiOverlayTypeByName ("circle");
+          KiiLoadOverlay (display, &overlay, 1, "red");
+          overlay.x = source->moments->Mx;
+          overlay.y = source->moments->My;
+          overlay.dx = 8;
+          overlay.dy = 8;
+          overlay.angle = 0.0;
+          overlay.type = KiiOverlayTypeByName ("circle");
+          KiiLoadOverlay (display, &overlay, 1, "blue");
+        }
 
         // pause and wait for user input:
@@ -2540 +2587 @@
     }
 
-    if (reshow) {
+    if (false && reshow) {
         psphotVisualShowMask (myKapa, readout->mask, "mask", 2);
         psphotVisualScaleImage (myKapa, readout->variance, readout->mask, "variance", 1.0, 1);
     }
-    psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", sqrt(factor), 0);
+
+    if (posImage) {
+        delImage = (psImage *) psBinaryOp(delImage, posImage, "-", readout->image);
+        psphotVisualScaleImage (myKapa, posImage, readout->mask, "posimage", sqrt(factor), 0);
+        psphotVisualScaleImage (myKapa, delImage, readout->mask, "delimage", sqrt(factor), 2);
+    }
+
+    psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", sqrt(factor), 1);
 
     pmVisualAskUser(NULL);