IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 23, 2009, 11:16:52 AM (17 years ago)
Author:
eugene
Message:

add some additional psf model visualization, adjustments to existing visuals; allow independent aperture and fit radii; return to psfMag - apMag for psf fit metric; use psf fit metric to test 2D order

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.c

    r25476 r25496  
    77#include "pmPSF.h"
    88#include "pmPSFtry.h"
     9#include "pmSource.h"
    910#include "pmSourceVisual.h"
    1011
    … …  
    1718
    1819static int kapa1 = -1;
     20static int kapa2 = -1;
    1921static bool plotPSF = true;
    20 // static int kapa2 = -1;
    2122// static int kapa3 = -1;
    2223
    … …  
    3334    Graphdata graphdata;
    3435
    35     if (!pmVisualIsVisual()) return true;
     36    // if (!pmVisualIsVisual()) return true;
    3637
    3738    if (kapa1 == -1) {
    … …  
    4445    }
    4546
    46     KapaClearPlots (kapa1);
     47    KapaClearSections (kapa1);
    4748    KapaInitGraph (&graphdata);
    4849
    … …  
    112113}
    113114
     115// to see the structure of the psf model, place the sources in a fake image 1/10th the size
     116// at their appropriate relative location. later sources stomp on earlier sources
     117bool pmSourceVisualShowModelFits (pmPSF *psf, psArray *sources, psImageMaskType maskVal) {
     118
     119    // if (!pmVisualIsVisual()) return true;
     120
     121    if (kapa2 == -1) {
     122        kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images");
     123        if (kapa2 == -1) {
     124            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     125            pmVisualSetVisual(false);
     126            return false;
     127        }
     128    }
     129
     130    // create images 1/10 scale:
     131    psImage *image = psImageAlloc (0.1*psf->fieldNx, 0.1*psf->fieldNy, PS_TYPE_F32);
     132    psImage *model = psImageAlloc (0.1*psf->fieldNx, 0.1*psf->fieldNy, PS_TYPE_F32);
     133    psImage *resid = psImageAlloc (0.1*psf->fieldNx, 0.1*psf->fieldNy, PS_TYPE_F32);
     134    psImageInit (image, 0.0);
     135    psImageInit (model, 0.0);
     136    psImageInit (resid, 0.0);
     137
     138    FILE *f = fopen ("stats.dat", "w");
     139    for (int i = 0; i < sources->n; i++) {
     140        pmSource *source = sources->data[i];
     141        if (!source) continue;
     142        if (!source->pixels) continue;
     143
     144        pmSourceCacheModel (source, maskVal);
     145        if (!source->modelFlux) continue;
     146
     147        pmModel *srcModel = pmSourceGetModel (NULL, source);
     148        if (!model) continue;
     149
     150        float norm = srcModel->params->data.F32[PM_PAR_I0];
     151
     152        int Xo = 0.1*srcModel->params->data.F32[PM_PAR_XPOS];
     153        int Yo = 0.1*srcModel->params->data.F32[PM_PAR_YPOS];
     154
     155        // insert source pixels in the image at 1/10th offset (XXX ignore centering for now)
     156        // XXX calculate the pixel scatter values, chisq
     157        for (int iy = 0; iy < source->pixels->numRows; iy++) {
     158            int jy = iy + Yo;
     159            if (jy >= image->numRows) continue;
     160            for (int ix = 0; ix < source->pixels->numCols; ix++) {
     161                int jx = ix + Xo;
     162                if (jx >= image->numCols) continue;
     163                if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) continue;
     164                if (source->modelFlux->data.F32[iy][ix] < 0.001) continue;
     165                image->data.F32[jy][jx] = source->pixels->data.F32[iy][ix];
     166                model->data.F32[jy][jx] = source->modelFlux->data.F32[iy][ix];
     167                resid->data.F32[jy][jx] = source->pixels->data.F32[iy][ix] - norm*source->modelFlux->data.F32[iy][ix];
     168            }
     169        }
     170
     171        fprintf (f, "%d, %d -> %f  %f  %f  %f\n", Xo, Yo, norm, source->psfMag, source->apMag, -2.5*log10(source->moments->Sum));
     172    }
     173    fclose (f);
     174
     175    // KapaClearSections (kapa2);
     176    pmVisualScaleImage (kapa2, image, "image", 0, true);
     177    pmVisualScaleImage (kapa2, model, "model", 1, true);
     178    pmVisualScaleImage (kapa2, resid, "resid", 2, true);
     179
     180    {
     181        psFits *fits = psFitsOpen ("image.fits", "w");
     182        psFitsWriteImage (fits, NULL, image, 0, NULL);
     183        psFitsClose (fits);
     184        fits = psFitsOpen ("model.fits", "w");
     185        psFitsWriteImage (fits, NULL, model, 0, NULL);
     186        psFitsClose (fits);
     187        fits = psFitsOpen ("resid.fits", "w");
     188        psFitsWriteImage (fits, NULL, resid, 0, NULL);
     189        psFitsClose (fits);
     190    }
     191
     192    psFree (image);
     193    psFree (model);
     194    psFree (resid);
     195
     196    // pause and wait for user input:
     197    // continue, save (provide name), ??
     198    char key[10];
     199    fprintf (stdout, "[c]ontinue? ");
     200    if (!fgets(key, 8, stdin)) {
     201        psWarning("Unable to read option");
     202    }
     203    return true;
     204}
     205
     206bool pmSourceVisualShowModelFit (pmSource *source) {
     207
     208    // if (!pmVisualIsVisual()) return true;
     209    if (!source->pixels) return false;
     210    if (!source->modelFlux) return false;
     211
     212    if (kapa2 == -1) {
     213        kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images");
     214        if (kapa2 == -1) {
     215            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     216            pmVisualSetVisual(false);
     217            return false;
     218        }
     219    }
     220
     221    // KapaClearSections (kapa2);
     222    pmVisualScaleImage (kapa2, source->pixels, "source", 0, false);
     223    pmVisualScaleImage (kapa2, source->modelFlux, "model", 1, false);
     224
     225    pmModel *model = pmSourceGetModel (NULL, source);
     226    float norm = model->params->data.F32[PM_PAR_I0];
     227
     228    psImage *resid = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
     229    for (int iy = 0; iy < source->pixels->numRows; iy++) {
     230        for (int ix = 0; ix < source->pixels->numCols; ix++) {
     231            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
     232                resid->data.F32[iy][ix] = NAN;
     233                continue;
     234            }
     235            resid->data.F32[iy][ix] = source->pixels->data.F32[iy][ix] - norm*source->modelFlux->data.F32[iy][ix];
     236        }
     237    }
     238    pmVisualScaleImage (kapa2, resid, "resid", 2, false);
     239
     240    psFree (resid);
     241
     242    // pause and wait for user input:
     243    // continue, save (provide name), ??
     244    char key[10];
     245    fprintf (stdout, "[c]ontinue? ");
     246    if (!fgets(key, 8, stdin)) {
     247        psWarning("Unable to read option");
     248    }
     249    return true;
     250}
     251
    114252bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) {
    115253
    … …  
    121259
    122260    // XXX skip for now
    123     return true;
     261    // return true;
    124262
    125263    if (kapa1 == -1) {
    … …  
    135273    KapaInitGraph (&graphdata);
    136274
    137     float min = +1e32;
    138     float max = -1e32;
    139     float Min = +1e32;
    140     float Max = -1e32;
     275    // float min = +1e32;
     276    // float max = -1e32;
     277    // float Min = +1e32;
     278    // float Max = -1e32;
     279
     280    float Xmin = +1e32;
     281    float Xmax = -1e32;
     282    float Ymin = +1e32;
     283    float Ymax = -1e32;
     284    float Fmin = +1e32;
     285    float Fmax = -1e32;
    141286
    142287    psVector *resid = psVectorAlloc (x->n, PS_TYPE_F32);
    143288    psVector *model = psVectorAlloc (x->n, PS_TYPE_F32);
    144289
     290    psVector *xm = psVectorAlloc (x->n, PS_TYPE_F32);
     291    psVector *ym = psVectorAlloc (x->n, PS_TYPE_F32);
     292    psVector *Fm = psVectorAlloc (x->n, PS_TYPE_F32);
     293
     294    int n = 0;
    145295    for (int i = 0; i < x->n; i++) {
    146296        model->data.F32[i] = pmTrend2DEval (trend, x->data.F32[i], y->data.F32[i]);
    147297        resid->data.F32[i] = param->data.F32[i] - model->data.F32[i];
    148298        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
    149         min = PS_MIN (min, resid->data.F32[i]);
    150         max = PS_MAX (max, resid->data.F32[i]);
    151         Min = PS_MIN (min, param->data.F32[i]);
    152         Max = PS_MAX (max, param->data.F32[i]);
    153     }
    154 
    155     psVector *xn = psVectorAlloc (x->n, PS_TYPE_F32);
    156     psVector *yn = psVectorAlloc (x->n, PS_TYPE_F32);
    157     psVector *zn = psVectorAlloc (x->n, PS_TYPE_F32);
    158     psVector *Zn = psVectorAlloc (x->n, PS_TYPE_F32);
    159     psVector *Fn = psVectorAlloc (x->n, PS_TYPE_F32);
    160     for (int i = 0; i < x->n; i++) {
    161         xn->data.F32[i] = x->data.F32[i] / 5000.0;
    162         yn->data.F32[i] = y->data.F32[i] / 5000.0;
    163         zn->data.F32[i] = (resid->data.F32[i] - min) / (max - min);
    164         Zn->data.F32[i] = (param->data.F32[i] - Min) / (Max - Min);
    165         Fn->data.F32[i] = (model->data.F32[i] - Min) / (Max - Min);
    166     }
     299        Xmin = PS_MIN (Xmin, x->data.F32[i]);
     300        Xmax = PS_MAX (Xmax, x->data.F32[i]);
     301        Ymin = PS_MIN (Ymin, y->data.F32[i]);
     302        Ymax = PS_MAX (Ymax, y->data.F32[i]);
     303        Fmin = PS_MIN (Fmin, param->data.F32[i]);
     304        Fmax = PS_MAX (Fmax, param->data.F32[i]);
     305        xm->data.F32[n] = x->data.F32[i];
     306        ym->data.F32[n] = y->data.F32[i];
     307        Fm->data.F32[n] = param->data.F32[i];
     308        n++;
     309    }
     310    xm->n = ym->n = Fm->n = n;
     311
     312    // psVector *xn = psVectorAlloc (x->n, PS_TYPE_F32);
     313    // psVector *yn = psVectorAlloc (x->n, PS_TYPE_F32);
     314    // psVector *zn = psVectorAlloc (x->n, PS_TYPE_F32);
     315    // psVector *Zn = psVectorAlloc (x->n, PS_TYPE_F32);
     316    // psVector *Fn = psVectorAlloc (x->n, PS_TYPE_F32);
     317    // for (int i = 0; i < x->n; i++) {
     318    //     xn->data.F32[i] = x->data.F32[i] / 5000.0;
     319    //     yn->data.F32[i] = y->data.F32[i] / 5000.0;
     320    //     zn->data.F32[i] = (resid->data.F32[i] - min) / (max - min);
     321    //     Zn->data.F32[i] = (param->data.F32[i] - Min) / (Max - Min);
     322    //     Fn->data.F32[i] = (model->data.F32[i] - Min) / (Max - Min);
     323    // }
    167324
    168325    // view 1 on resid
    169     section.dx = 0.5;
    170     section.dy = 0.33;
     326    section.dx = 1.0;
     327    section.dy = 0.5;
    171328    section.x = 0.0;
    172329    section.y = 0.0;
    … …  
    175332    KapaSetSection (kapa1, &section);
    176333    psFree (section.name);
    177     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
     334
     335    graphdata.color = KapaColorByName ("black");
     336    graphdata.xmin = Xmin;
     337    graphdata.xmax = Xmax;
     338    graphdata.ymin = Fmin;
     339    graphdata.ymax = Fmax;
     340
     341    {
     342        float range;
     343        range = graphdata.xmax - graphdata.xmin;
     344        graphdata.xmax += 0.05*range;
     345        graphdata.xmin -= 0.05*range;
     346        range = graphdata.ymax - graphdata.ymin;
     347        graphdata.ymax += 0.05*range;
     348        graphdata.ymin -= 0.05*range;
     349    }
     350
     351    KapaSetLimits (kapa1, &graphdata);
     352    KapaSetFont (kapa1, "helvetica", 14);
     353    KapaBox (kapa1, &graphdata);
     354    KapaSendLabel (kapa1, "X (pixels)", KAPA_LABEL_XM);
     355    KapaSendLabel (kapa1, "Model Param", KAPA_LABEL_YM);
     356
     357    graphdata.ptype = 2;
     358    graphdata.size = 1.0;
     359    graphdata.style = 2;
     360    KapaPrepPlot (kapa1,   x->n, &graphdata);
     361    KapaPlotVector (kapa1, x->n, x->data.F32, "x");
     362    KapaPlotVector (kapa1, x->n, param->data.F32, "y");
     363
     364    graphdata.color = KapaColorByName ("red");
     365    graphdata.ptype = 1;
     366    KapaPrepPlot (kapa1,   xm->n, &graphdata);
     367    KapaPlotVector (kapa1, xm->n, xm->data.F32, "x");
     368    KapaPlotVector (kapa1, xm->n, Fm->data.F32, "y");
     369
     370    graphdata.color = KapaColorByName ("blue");
     371    graphdata.ptype = 1;
     372    KapaPrepPlot (kapa1,   x->n, &graphdata);
     373    KapaPlotVector (kapa1, x->n, x->data.F32, "x");
     374    KapaPlotVector (kapa1, x->n, model->data.F32, "y");
    178375
    179376    // view 2 on resid
    180     section.dx = 0.5;
    181     section.dy = 0.33;
    182     section.x = 0.5;
    183     section.y = 0.0;
     377    section.dx = 1.0;
     378    section.dy = 0.5;
     379    section.x = 0.0;
     380    section.y = 0.5;
    184381    section.name = NULL;
    185382    psStringAppend (&section.name, "a2");
    186383    KapaSetSection (kapa1, &section);
    187384    psFree (section.name);
    188     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
    189 
     385
     386    graphdata.color = KapaColorByName ("black");
     387    graphdata.xmin = Ymin;
     388    graphdata.xmax = Ymax;
     389    graphdata.ymin = Fmin;
     390    graphdata.ymax = Fmax;
     391    {
     392        float range;
     393        range = graphdata.xmax - graphdata.xmin;
     394        graphdata.xmax += 0.05*range;
     395        graphdata.xmin -= 0.05*range;
     396        range = graphdata.ymax - graphdata.ymin;
     397        graphdata.ymax += 0.05*range;
     398        graphdata.ymin -= 0.05*range;
     399    }
     400
     401    KapaSetLimits (kapa1, &graphdata);
     402    KapaSetFont (kapa1, "helvetica", 14);
     403    KapaBox (kapa1, &graphdata);
     404    KapaSendLabel (kapa1, "Y (pixels)", KAPA_LABEL_XM);
     405    KapaSendLabel (kapa1, "Model Param", KAPA_LABEL_YM);
     406
     407    graphdata.ptype = 2;
     408    graphdata.size = 1.0;
     409    graphdata.style = 2;
     410    KapaPrepPlot (kapa1,   y->n, &graphdata);
     411    KapaPlotVector (kapa1, y->n, y->data.F32, "x");
     412    KapaPlotVector (kapa1, y->n, param->data.F32, "y");
     413
     414    graphdata.color = KapaColorByName ("red");
     415    graphdata.ptype = 1;
     416    KapaPrepPlot (kapa1,   xm->n, &graphdata);
     417    KapaPlotVector (kapa1, xm->n, ym->data.F32, "x");
     418    KapaPlotVector (kapa1, xm->n, Fm->data.F32, "y");
     419
     420    graphdata.color = KapaColorByName ("blue");
     421    graphdata.ptype = 1;
     422    KapaPrepPlot (kapa1,   y->n, &graphdata);
     423    KapaPlotVector (kapa1, y->n, y->data.F32, "x");
     424    KapaPlotVector (kapa1, y->n, model->data.F32, "y");
     425
     426# if (0)
    190427    // view 3 on resid
    191428    section.dx = 0.5;
    … …  
    231468    psFree (section.name);
    232469    pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Fn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
     470# endif
    233471
    234472    psFree (resid);
    235 
    236     psFree (xn);
    237     psFree (yn);
    238     psFree (zn);
    239     psFree (Zn);
     473    psFree (model);
    240474
    241475    // pause and wait for user input:
Note: See TracChangeset for help on using the changeset viewer.