IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 2, 2009, 3:11:32 PM (17 years ago)
Author:
eugene
Message:

check in changes from genes development branch : extensive changes to moments calculation, psf model generation, aperture residuals

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmSourceVisual.c

    r23242 r25754  
    55#include <pslib.h>
    66#include "pmTrend2D.h"
     7#include "pmPSF.h"
     8#include "pmPSFtry.h"
     9#include "pmSource.h"
    710#include "pmSourceVisual.h"
    811
    … …  
    1518
    1619static int kapa1 = -1;
     20static int kapa2 = -1;
    1721static bool plotPSF = true;
    18 // static int kapa2 = -1;
    1922// static int kapa3 = -1;
    2023
    … …  
    2730bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi);
    2831
    29 
    30 bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) {
    31 
    32     KapaSection section;  // put the positive profile in one and the residuals in another?
     32bool pmSourceVisualPlotPSFMetric (pmPSFtry *psfTry) {
    3333
    3434    Graphdata graphdata;
    3535
    36     if (!pmVisualIsVisual() || !plotPSF) return true;
     36    if (!pmVisualIsVisual()) return true;
    3737
    3838    if (kapa1 == -1) {
    … …  
    4545    }
    4646
     47    KapaClearSections (kapa1);
     48    KapaInitGraph (&graphdata);
     49
     50    psVector *x = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32);
     51    psVector *y = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32);
     52    psVector *dy = psVectorAllocEmpty(psfTry->sources->n, PS_TYPE_F32);
     53
     54    graphdata.xmin = +32.0;
     55    graphdata.xmax = -32.0;
     56    graphdata.ymin = +32.0;
     57    graphdata.ymax = -32.0;
     58
     59    // construct the plot vectors
     60    int n = 0;
     61    for (int i = 0; i < psfTry->sources->n; i++) {
     62        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
     63        x->data.F32[n] = psfTry->fitMag->data.F32[i];
     64        y->data.F32[n] = psfTry->metric->data.F32[i];
     65        dy->data.F32[n] = psfTry->metricErr->data.F32[i];
     66        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
     67        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
     68        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
     69        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);
     70        n++;
     71    }
     72    x->n = y->n = dy->n = n;
     73
     74    float range;
     75    range = graphdata.xmax - graphdata.xmin;
     76    graphdata.xmax += 0.05*range;
     77    graphdata.xmin -= 0.05*range;
     78    range = graphdata.ymax - graphdata.ymin;
     79    graphdata.ymax += 0.05*range;
     80    graphdata.ymin -= 0.05*range;
     81
     82    // better choice for range?
     83    // graphdata.xmin = -17.0;
     84    // graphdata.xmax =  -9.0;
     85    graphdata.ymin = -0.51;
     86    graphdata.ymax = +0.51;
     87
     88    KapaSetLimits (kapa1, &graphdata);
     89
     90    KapaSetFont (kapa1, "helvetica", 14);
     91    KapaBox (kapa1, &graphdata);
     92    KapaSendLabel (kapa1, "PSF Mag", KAPA_LABEL_XM);
     93    KapaSendLabel (kapa1, "Ap Mag - PSF Mag", KAPA_LABEL_YM);
     94
     95    graphdata.color = KapaColorByName ("black");
     96    graphdata.ptype = 2;
     97    graphdata.size = 0.5;
     98    graphdata.style = 2;
     99    graphdata.etype |= 0x01;
     100
     101    KapaPrepPlot (kapa1, n, &graphdata);
     102    KapaPlotVector (kapa1, n, x->data.F32, "x");
     103    KapaPlotVector (kapa1, n, y->data.F32, "y");
     104    KapaPlotVector (kapa1, n, dy->data.F32, "dym");
     105    KapaPlotVector (kapa1, n, dy->data.F32, "dyp");
     106
     107    psFree (x);
     108    psFree (y);
     109    psFree (dy);
     110
     111    pmVisualAskUser(NULL);
     112    return true;
     113}
     114
     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    for (int i = sources->n - 1; i >= 0; i--) {
     139        pmSource *source = sources->data[i];
     140        if (!source) continue;
     141        if (!source->pixels) continue;
     142
     143        pmSourceCacheModel (source, maskVal);
     144        if (!source->modelFlux) continue;
     145
     146        pmModel *srcModel = pmSourceGetModel (NULL, source);
     147        if (!model) continue;
     148
     149        float norm = srcModel->params->data.F32[PM_PAR_I0];
     150
     151        int Xo = 0.1*srcModel->params->data.F32[PM_PAR_XPOS];
     152        int Yo = 0.1*srcModel->params->data.F32[PM_PAR_YPOS];
     153
     154        // insert source pixels in the image at 1/10th offset
     155        for (int iy = 0; iy < source->pixels->numRows; iy++) {
     156            int jy = iy + Yo;
     157            if (jy >= image->numRows) continue;
     158            for (int ix = 0; ix < source->pixels->numCols; ix++) {
     159                int jx = ix + Xo;
     160                if (jx >= image->numCols) continue;
     161                if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) continue;
     162                if (source->modelFlux->data.F32[iy][ix] < 0.001) continue;
     163                image->data.F32[jy][jx] = source->pixels->data.F32[iy][ix];
     164                model->data.F32[jy][jx] = source->modelFlux->data.F32[iy][ix];
     165                resid->data.F32[jy][jx] = source->pixels->data.F32[iy][ix] - norm*source->modelFlux->data.F32[iy][ix];
     166            }
     167        }
     168    }
     169
     170    // KapaClearSections (kapa2);
     171    pmVisualScaleImage (kapa2, image, "image", 0, true);
     172    pmVisualScaleImage (kapa2, model, "model", 1, true);
     173    pmVisualScaleImage (kapa2, resid, "resid", 2, true);
     174
     175# ifdef DEBUG
     176    {
     177        psFits *fits = psFitsOpen ("image.fits", "w");
     178        psFitsWriteImage (fits, NULL, image, 0, NULL);
     179        psFitsClose (fits);
     180        fits = psFitsOpen ("model.fits", "w");
     181        psFitsWriteImage (fits, NULL, model, 0, NULL);
     182        psFitsClose (fits);
     183        fits = psFitsOpen ("resid.fits", "w");
     184        psFitsWriteImage (fits, NULL, resid, 0, NULL);
     185        psFitsClose (fits);
     186    }
     187# endif
     188
     189    psFree (image);
     190    psFree (model);
     191    psFree (resid);
     192
     193    pmVisualAskUser(NULL);
     194    return true;
     195}
     196
     197bool pmSourceVisualShowModelFit (pmSource *source) {
     198
     199    if (!pmVisualIsVisual()) return true;
     200    if (!source->pixels) return false;
     201    if (!source->modelFlux) return false;
     202
     203    if (kapa2 == -1) {
     204        kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images");
     205        if (kapa2 == -1) {
     206            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     207            pmVisualSetVisual(false);
     208            return false;
     209        }
     210    }
     211
     212    // KapaClearSections (kapa2);
     213    pmVisualScaleImage (kapa2, source->pixels, "source", 0, false);
     214    pmVisualScaleImage (kapa2, source->modelFlux, "model", 1, false);
     215
     216    pmModel *model = pmSourceGetModel (NULL, source);
     217    float norm = model->params->data.F32[PM_PAR_I0];
     218
     219    psImage *resid = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
     220    for (int iy = 0; iy < source->pixels->numRows; iy++) {
     221        for (int ix = 0; ix < source->pixels->numCols; ix++) {
     222            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
     223                resid->data.F32[iy][ix] = NAN;
     224                continue;
     225            }
     226            resid->data.F32[iy][ix] = source->pixels->data.F32[iy][ix] - norm*source->modelFlux->data.F32[iy][ix];
     227        }
     228    }
     229    pmVisualScaleImage (kapa2, resid, "resid", 2, false);
     230
     231    psFree (resid);
     232
     233    pmVisualAskUser(NULL);
     234    return true;
     235}
     236
     237bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) {
     238
     239    KapaSection section;  // put the positive profile in one and the residuals in another?
     240
     241    Graphdata graphdata;
     242
     243    if (!pmVisualIsVisual() || !plotPSF) return true;
     244
     245    if (kapa1 == -1) {
     246        kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots");
     247        if (kapa1 == -1) {
     248            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     249            pmVisualSetVisual(false);
     250            return false;
     251        }
     252    }
     253
    47254    KapaClearPlots (kapa1);
    48255    KapaInitGraph (&graphdata);
    49256
    50     float min = +1e32;
    51     float max = -1e32;
    52     float Min = +1e32;
    53     float Max = -1e32;
     257    float Xmin = +1e32;
     258    float Xmax = -1e32;
     259    float Ymin = +1e32;
     260    float Ymax = -1e32;
     261    float Fmin = +1e32;
     262    float Fmax = -1e32;
    54263
    55264    psVector *resid = psVectorAlloc (x->n, PS_TYPE_F32);
    56265    psVector *model = psVectorAlloc (x->n, PS_TYPE_F32);
    57266
     267    psVector *xm = psVectorAlloc (x->n, PS_TYPE_F32);
     268    psVector *ym = psVectorAlloc (x->n, PS_TYPE_F32);
     269    psVector *Fm = psVectorAlloc (x->n, PS_TYPE_F32);
     270
     271    int n = 0;
    58272    for (int i = 0; i < x->n; i++) {
    59273        model->data.F32[i] = pmTrend2DEval (trend, x->data.F32[i], y->data.F32[i]);
    60274        resid->data.F32[i] = param->data.F32[i] - model->data.F32[i];
    61275        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
    62         min = PS_MIN (min, resid->data.F32[i]);
    63         max = PS_MAX (max, resid->data.F32[i]);
    64         Min = PS_MIN (min, param->data.F32[i]);
    65         Max = PS_MAX (max, param->data.F32[i]);
    66     }
    67 
    68     psVector *xn = psVectorAlloc (x->n, PS_TYPE_F32);
    69     psVector *yn = psVectorAlloc (x->n, PS_TYPE_F32);
    70     psVector *zn = psVectorAlloc (x->n, PS_TYPE_F32);
    71     psVector *Zn = psVectorAlloc (x->n, PS_TYPE_F32);
    72     psVector *Fn = psVectorAlloc (x->n, PS_TYPE_F32);
    73     for (int i = 0; i < x->n; i++) {
    74         xn->data.F32[i] = x->data.F32[i] / 5000.0;
    75         yn->data.F32[i] = y->data.F32[i] / 5000.0;
    76         zn->data.F32[i] = (resid->data.F32[i] - min) / (max - min);
    77         Zn->data.F32[i] = (param->data.F32[i] - Min) / (Max - Min);
    78         Fn->data.F32[i] = (model->data.F32[i] - Min) / (Max - Min);
    79     }
     276        Xmin = PS_MIN (Xmin, x->data.F32[i]);
     277        Xmax = PS_MAX (Xmax, x->data.F32[i]);
     278        Ymin = PS_MIN (Ymin, y->data.F32[i]);
     279        Ymax = PS_MAX (Ymax, y->data.F32[i]);
     280        Fmin = PS_MIN (Fmin, param->data.F32[i]);
     281        Fmax = PS_MAX (Fmax, param->data.F32[i]);
     282        xm->data.F32[n] = x->data.F32[i];
     283        ym->data.F32[n] = y->data.F32[i];
     284        Fm->data.F32[n] = param->data.F32[i];
     285        n++;
     286    }
     287    xm->n = ym->n = Fm->n = n;
    80288
    81289    // view 1 on resid
    82     section.dx = 0.5;
    83     section.dy = 0.33;
     290    section.dx = 1.0;
     291    section.dy = 0.5;
    84292    section.x = 0.0;
    85293    section.y = 0.0;
    … …  
    88296    KapaSetSection (kapa1, &section);
    89297    psFree (section.name);
    90     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
     298
     299    graphdata.color = KapaColorByName ("black");
     300    graphdata.xmin = Xmin;
     301    graphdata.xmax = Xmax;
     302    graphdata.ymin = Fmin;
     303    graphdata.ymax = Fmax;
     304
     305    {
     306        float range;
     307        range = graphdata.xmax - graphdata.xmin;
     308        graphdata.xmax += 0.05*range;
     309        graphdata.xmin -= 0.05*range;
     310        range = graphdata.ymax - graphdata.ymin;
     311        graphdata.ymax += 0.05*range;
     312        graphdata.ymin -= 0.05*range;
     313    }
     314
     315    KapaSetLimits (kapa1, &graphdata);
     316    KapaSetFont (kapa1, "helvetica", 14);
     317    KapaBox (kapa1, &graphdata);
     318    KapaSendLabel (kapa1, "X (pixels)", KAPA_LABEL_XM);
     319    KapaSendLabel (kapa1, "Model Param", KAPA_LABEL_YM);
     320
     321    graphdata.ptype = 2;
     322    graphdata.size = 1.0;
     323    graphdata.style = 2;
     324    KapaPrepPlot (kapa1,   x->n, &graphdata);
     325    KapaPlotVector (kapa1, x->n, x->data.F32, "x");
     326    KapaPlotVector (kapa1, x->n, param->data.F32, "y");
     327
     328    graphdata.color = KapaColorByName ("red");
     329    graphdata.ptype = 1;
     330    KapaPrepPlot (kapa1,   xm->n, &graphdata);
     331    KapaPlotVector (kapa1, xm->n, xm->data.F32, "x");
     332    KapaPlotVector (kapa1, xm->n, Fm->data.F32, "y");
     333
     334    graphdata.color = KapaColorByName ("blue");
     335    graphdata.ptype = 1;
     336    KapaPrepPlot (kapa1,   x->n, &graphdata);
     337    KapaPlotVector (kapa1, x->n, x->data.F32, "x");
     338    KapaPlotVector (kapa1, x->n, model->data.F32, "y");
    91339
    92340    // view 2 on resid
    93     section.dx = 0.5;
    94     section.dy = 0.33;
    95     section.x = 0.5;
    96     section.y = 0.0;
     341    section.dx = 1.0;
     342    section.dy = 0.5;
     343    section.x = 0.0;
     344    section.y = 0.5;
    97345    section.name = NULL;
    98346    psStringAppend (&section.name, "a2");
    99347    KapaSetSection (kapa1, &section);
    100348    psFree (section.name);
    101     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
    102 
    103     // view 3 on resid
    104     section.dx = 0.5;
    105     section.dy = 0.33;
    106     section.x = 0.0;
    107     section.y = 0.33;
    108     section.name = NULL;
    109     psStringAppend (&section.name, "a3");
    110     KapaSetSection (kapa1, &section);
    111     psFree (section.name);
    112     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Zn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
    113 
    114     // view 4 on resid
    115     section.dx = 0.5;
    116     section.dy = 0.33;
    117     section.x = 0.5;
    118     section.y = 0.33;
    119     section.name = NULL;
    120     psStringAppend (&section.name, "a4");
    121     KapaSetSection (kapa1, &section);
    122     psFree (section.name);
    123     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Zn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
    124 
    125     // view 5 on resid
    126     section.dx = 0.5;
    127     section.dy = 0.33;
    128     section.x = 0.0;
    129     section.y = 0.66;
    130     section.name = NULL;
    131     psStringAppend (&section.name, "a5");
    132     KapaSetSection (kapa1, &section);
    133     psFree (section.name);
    134     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Fn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
    135 
    136     // view 6 on resid
    137     section.dx = 0.5;
    138     section.dy = 0.33;
    139     section.x = 0.5;
    140     section.y = 0.66;
    141     section.name = NULL;
    142     psStringAppend (&section.name, "a6");
    143     KapaSetSection (kapa1, &section);
    144     psFree (section.name);
    145     pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Fn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);
     349
     350    graphdata.color = KapaColorByName ("black");
     351    graphdata.xmin = Ymin;
     352    graphdata.xmax = Ymax;
     353    graphdata.ymin = Fmin;
     354    graphdata.ymax = Fmax;
     355    {
     356        float range;
     357        range = graphdata.xmax - graphdata.xmin;
     358        graphdata.xmax += 0.05*range;
     359        graphdata.xmin -= 0.05*range;
     360        range = graphdata.ymax - graphdata.ymin;
     361        graphdata.ymax += 0.05*range;
     362        graphdata.ymin -= 0.05*range;
     363    }
     364
     365    KapaSetLimits (kapa1, &graphdata);
     366    KapaSetFont (kapa1, "helvetica", 14);
     367    KapaBox (kapa1, &graphdata);
     368    KapaSendLabel (kapa1, "Y (pixels)", KAPA_LABEL_XM);
     369    KapaSendLabel (kapa1, "Model Param", KAPA_LABEL_YM);
     370
     371    graphdata.ptype = 2;
     372    graphdata.size = 1.0;
     373    graphdata.style = 2;
     374    KapaPrepPlot (kapa1,   y->n, &graphdata);
     375    KapaPlotVector (kapa1, y->n, y->data.F32, "x");
     376    KapaPlotVector (kapa1, y->n, param->data.F32, "y");
     377
     378    graphdata.color = KapaColorByName ("red");
     379    graphdata.ptype = 1;
     380    KapaPrepPlot (kapa1,   xm->n, &graphdata);
     381    KapaPlotVector (kapa1, xm->n, ym->data.F32, "x");
     382    KapaPlotVector (kapa1, xm->n, Fm->data.F32, "y");
     383
     384    graphdata.color = KapaColorByName ("blue");
     385    graphdata.ptype = 1;
     386    KapaPrepPlot (kapa1,   y->n, &graphdata);
     387    KapaPlotVector (kapa1, y->n, y->data.F32, "x");
     388    KapaPlotVector (kapa1, y->n, model->data.F32, "y");
    146389
    147390    psFree (resid);
    148 
    149     psFree (xn);
    150     psFree (yn);
    151     psFree (zn);
    152     psFree (Zn);
     391    psFree (model);
    153392
    154393    // pause and wait for user input:
    … …  
    159398}
    160399
    161 // send in normalized points
     400// Somewhat broken 3D plotting function (was used by pmSourceVisualPSFModelResid, but not anymore)
    162401bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi) {
     402
     403    return true;
    163404
    164405    psVector *xv = psVectorAlloc (PS_MAX(6, 2*xn->n), PS_TYPE_F32);
    … …  
    192433    KapaSetLimits (myKapa, graphdata);
    193434
    194     // KapaSetFont (myKapa, "helvetica", 14);
    195     // KapaBox (myKapa, graphdata);
    196     // KapaSendLabel (myKapa, "&ss&h_x| (pixels)", KAPA_LABEL_XM);
    197     // KapaSendLabel (myKapa, "&ss&h_y| (pixels)", KAPA_LABEL_YM);
    198 
    199435    graphdata->color = KapaColorByName ("black");
    200436    graphdata->ptype = 100;
Note: See TracChangeset for help on using the changeset viewer.