IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 26, 2012, 1:05:27 PM (14 years ago)
Author:
eugene
Message:

visual adjustment to examine saturated stars

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20120805/psphot/src/psphotVisual.c

    r31452 r34360  
    11391139    // after displaying (as an image) the psf stars, we cycle throught them and display their
    11401140    // radial profiles:
    1141     psphotVisualPlotRadialProfiles (recipe, sources);
     1141    psphotVisualPlotRadialProfiles (recipe, sources, PM_SOURCE_MODE_PSFSTAR);
    11421142
    11431143    return true;
    … …  
    13021302    int ng = 0;
    13031303    int nb = 0;
    1304     float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
    1305     float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
     1304
     1305    float Xo = NAN;
     1306    float Yo = NAN;
     1307
     1308    if (source->modelPSF) {
     1309        Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
     1310        Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
     1311    } else {
     1312        Xo = source->moments->Mx;
     1313        Yo = source->moments->My;
     1314    }
     1315
    13061316    for (int iy = 0; iy < source->pixels->numRows; iy++) {
    13071317        for (int ix = 0; ix < source->pixels->numCols; ix++) {
    … …  
    13221332    }
    13231333
    1324     // generate model profiles (major and minor axis):
    1325     // create a model with theta = 0.0 so major and minor axes are equiv to x and y:
    1326     psEllipseShape rawShape, rotShape;
    1327 
    1328     rawShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
    1329     rawShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
    1330     rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
    1331 
    1332     psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0);
    1333 
    1334     axes.theta = 0.0;
    1335 
    1336     rotShape = psEllipseAxesToShape (axes);
    1337 
    1338     psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32);
    1339     for (int i = 0; i < source->modelPSF->params->n; i++) {
    1340         params->data.F32[i] = source->modelPSF->params->data.F32[i];
    1341     }
    1342     params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2;
    1343     params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2;
    1344     params->data.F32[PM_PAR_SXY] = rotShape.sxy;
    1345     params->data.F32[PM_PAR_XPOS] = 0.0;
    1346     params->data.F32[PM_PAR_YPOS] = 0.0;
    1347 
    1348     psVector *rmod = psVectorAlloc(300, PS_TYPE_F32);
    1349     psVector *fmaj = psVectorAlloc(300, PS_TYPE_F32);
    1350     psVector *fmin = psVectorAlloc(300, PS_TYPE_F32);
    1351 
    1352     psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    1353 
    1354     float r = 0.0;
    1355     for (int i = 0; i < rmod->n; i++) {
    1356         r = i*0.1;
    1357         rmod->data.F32[i] = r;
    1358 
    1359         coord->data.F32[1] = r;
    1360         coord->data.F32[0] = 0.0;
    1361         fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
    1362 
    1363         coord->data.F32[0] = r;
    1364         coord->data.F32[1] = 0.0;
    1365         fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
    1366     }
    1367     psFree (coord);
    1368     psFree (params);
    1369 
    1370     float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
    1371     float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
    1372     if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR);
    1373 
    1374     psEllipseMoments emoments;
    1375     emoments.x2 = source->moments->Mxx;
    1376     emoments.xy = source->moments->Mxy;
    1377     emoments.y2 = source->moments->Myy;
    1378     axes = psEllipseMomentsToAxes (emoments, 20.0);
    1379     float MOMENTS_MAJOR = 2.355*axes.major;
    1380     float MOMENTS_MINOR = 2.355*axes.minor;
    1381 
    1382     float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
    1383 
    1384     // reset source Add/Sub state to recorded
    1385     if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    1386 
    13871334    KapaInitGraph (&graphdata);
    13881335
    … …  
    14171364    KapaPlotVector (myKapa, nb, rb->data.F32, "x");
    14181365    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
    1419 
    1420     graphdata.color = KapaColorByName ("blue");
    1421     graphdata.ptype = 0;
    1422     graphdata.size = 0.0;
    1423     graphdata.style = 0;
    1424     KapaPrepPlot   (myKapa, rmod->n, &graphdata);
    1425     KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
    1426     KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
    1427     plotline (myKapa, &graphdata, 0.0, logHM, 30.0, logHM);
    1428     plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, 0.0, 0.5*FWHM_MINOR, 5.0);
    1429     graphdata.ltype = 1;
    1430     plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, 0.0, 0.5*MOMENTS_MINOR, 5.0);
    1431     graphdata.ltype = 0;
    1432        
    1433     graphdata.color = KapaColorByName ("green");
    1434     graphdata.ptype = 0;
    1435     graphdata.size = 0.0;
    1436     graphdata.style = 0;
    1437     KapaPrepPlot   (myKapa, rmod->n, &graphdata);
    1438     KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
    1439     KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
    1440     plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, 0.0, 0.5*FWHM_MAJOR, 5.0);
    1441     graphdata.ltype = 1;
    1442     plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, 0.0, 0.5*MOMENTS_MAJOR, 5.0);
    1443     graphdata.ltype = 0;
    1444        
    1445     for (int i = 0; i < rmod->n; i++) {
    1446         rmod->data.F32[i] = log10(rmod->data.F32[i]);
    1447     }
    14481366
    14491367    // ** loglog **
    … …  
    14791397    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
    14801398
    1481     graphdata.color = KapaColorByName ("blue");
    1482     graphdata.ptype = 0;
    1483     graphdata.size = 0.0;
    1484     graphdata.style = 0;
    1485     KapaPrepPlot   (myKapa, rmod->n, &graphdata);
    1486     KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
    1487     KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
    1488 
    1489     graphdata.color = KapaColorByName ("green");
    1490     graphdata.ptype = 0;
    1491     graphdata.size = 0.0;
    1492     graphdata.style = 0;
    1493     KapaPrepPlot   (myKapa, rmod->n, &graphdata);
    1494     KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
    1495     KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
    1496 
    1497     psFree (rmod);
    1498     psFree (fmin);
    1499     psFree (fmaj);
     1399    if (source->modelPSF) {
     1400        // generate model profiles (major and minor axis):
     1401        // create a model with theta = 0.0 so major and minor axes are equiv to x and y:
     1402        psEllipseShape rawShape, rotShape;
     1403
     1404        rawShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
     1405        rawShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
     1406        rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
     1407
     1408        psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0);
     1409
     1410        axes.theta = 0.0;
     1411
     1412        rotShape = psEllipseAxesToShape (axes);
     1413
     1414        psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32);
     1415        for (int i = 0; i < source->modelPSF->params->n; i++) {
     1416            params->data.F32[i] = source->modelPSF->params->data.F32[i];
     1417        }
     1418        params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2;
     1419        params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2;
     1420        params->data.F32[PM_PAR_SXY] = rotShape.sxy;
     1421        params->data.F32[PM_PAR_XPOS] = 0.0;
     1422        params->data.F32[PM_PAR_YPOS] = 0.0;
     1423
     1424        psVector *rmod = psVectorAlloc(300, PS_TYPE_F32);
     1425        psVector *fmaj = psVectorAlloc(300, PS_TYPE_F32);
     1426        psVector *fmin = psVectorAlloc(300, PS_TYPE_F32);
     1427
     1428        psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     1429
     1430        float r = 0.0;
     1431        for (int i = 0; i < rmod->n; i++) {
     1432            r = i*0.1;
     1433            rmod->data.F32[i] = r;
     1434
     1435            coord->data.F32[1] = r;
     1436            coord->data.F32[0] = 0.0;
     1437            fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
     1438
     1439            coord->data.F32[0] = r;
     1440            coord->data.F32[1] = 0.0;
     1441            fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
     1442        }
     1443        psFree (coord);
     1444        psFree (params);
     1445
     1446        float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
     1447        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
     1448        if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR);
     1449
     1450        psEllipseMoments emoments;
     1451        emoments.x2 = source->moments->Mxx;
     1452        emoments.xy = source->moments->Mxy;
     1453        emoments.y2 = source->moments->Myy;
     1454        axes = psEllipseMomentsToAxes (emoments, 20.0);
     1455        float MOMENTS_MAJOR = 2.355*axes.major;
     1456        float MOMENTS_MINOR = 2.355*axes.minor;
     1457
     1458        float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
     1459
     1460        // reset source Add/Sub state to recorded
     1461        if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     1462
     1463        // ** linlog **
     1464        KapaSelectSection (myKapa, "linlog");
     1465        KapaGetGraphData (myKapa, &graphdata);
     1466        graphdata.color = KapaColorByName ("blue");
     1467        graphdata.ptype = 0;
     1468        graphdata.size = 0.0;
     1469        graphdata.style = 0;
     1470        KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1471        KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1472        KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
     1473        plotline (myKapa, &graphdata, graphdata.xmin, logHM, graphdata.xmax, logHM);
     1474        plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, graphdata.ymin, 0.5*FWHM_MINOR, graphdata.ymax);
     1475        graphdata.ltype = 1;
     1476        plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, graphdata.ymin, 0.5*MOMENTS_MINOR, graphdata.ymax);
     1477        graphdata.ltype = 0;
     1478       
     1479        graphdata.color = KapaColorByName ("green");
     1480        graphdata.ptype = 0;
     1481        graphdata.size = 0.0;
     1482        graphdata.style = 0;
     1483        KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1484        KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1485        KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
     1486        plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, graphdata.ymin, 0.5*FWHM_MAJOR, graphdata.ymax);
     1487        graphdata.ltype = 1;
     1488        plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, graphdata.ymin, 0.5*MOMENTS_MAJOR, graphdata.ymax);
     1489        graphdata.ltype = 0;
     1490       
     1491        for (int i = 0; i < rmod->n; i++) {
     1492            rmod->data.F32[i] = log10(rmod->data.F32[i]);
     1493        }
     1494
     1495        // ** loglog **
     1496        KapaSelectSection (myKapa, "loglog");
     1497        KapaGetGraphData (myKapa, &graphdata);
     1498        graphdata.color = KapaColorByName ("blue");
     1499        graphdata.ptype = 0;
     1500        graphdata.size = 0.0;
     1501        graphdata.style = 0;
     1502        KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1503        KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1504        KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
     1505
     1506        graphdata.color = KapaColorByName ("green");
     1507        graphdata.ptype = 0;
     1508        graphdata.size = 0.0;
     1509        graphdata.style = 0;
     1510        KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1511        KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1512        KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
     1513
     1514        psFree (rmod);
     1515        psFree (fmin);
     1516        psFree (fmaj);
     1517    }
    15001518
    15011519    psFree (rg);
    … …  
    15081526}
    15091527
    1510 bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources) {
     1528bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources, pmSourceMode showmode) {
    15111529
    15121530    KapaSection section;  // put the positive profile in one and the residuals in another?
    … …  
    15491567
    15501568        pmSource *source = sources->data[i];
    1551         if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
     1569        if (!(source->mode & showmode)) continue;
    15521570
    15531571        psphotVisualPlotRadialProfile (myKapa, source, maskVal);
Note: See TracChangeset for help on using the changeset viewer.