IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 3, 2010, 8:41:49 AM (16 years ago)
Author:
eugene
Message:

updates from trunk

Location:
branches/tap_branches
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/tap_branches

  • branches/tap_branches/psphot

  • branches/tap_branches/psphot/src

    • Property svn:ignore

      •  

        old new  
        1919psphotMomentsStudy
        2020psphotPetrosianStudy
         21psphotForced
         22psphotMakePSF
         23psphotStack

  • branches/tap_branches/psphot/src/psphotVisual.c

    r25852 r27838  
    6060}
    6161
     62bool psphotVisualEraseOverlays (int channel, char *overlay) {
     63
     64    int myKapa = psphotKapaChannel (channel);
     65    if (!(strcasecmp (overlay, "all"))) {
     66      KiiEraseOverlay (myKapa, "red");
     67      KiiEraseOverlay (myKapa, "green");
     68      KiiEraseOverlay (myKapa, "blue");
     69      KiiEraseOverlay (myKapa, "yellow");
     70      return true;
     71    }
     72    KiiEraseOverlay (myKapa, overlay);
     73    return true;
     74}
     75
    6276bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel) {
    6377
    … …  
    106120}
    107121
    108 bool psphotVisualScaleImage (int kapaFD, psImage *inImage, const char *name, int channel) {
     122bool psphotVisualScaleImage (int kapaFD, psImage *inImage, psImage *inMask, const char *name, int channel) {
    109123
    110124    KiiImage image;
    … …  
    116130    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    117131    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    118     if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) {
     132    if (!psImageBackground(stats, NULL, inImage, inMask, 0xffff, rng)) {
    119133        fprintf (stderr, "failed to get background values\n");
    120134        return false;
    … …  
    171185    if (kapa == -1) return false;
    172186
    173     // psphotVisualShowMask (kapa, readout->mask, "mask", 2);
    174     psphotVisualScaleImage (kapa, readout->variance, "variance", 1);
    175     psphotVisualScaleImage (kapa, readout->image, "image", 0);
     187    psphotVisualShowMask (kapa, readout->mask, "mask", 2);
     188    psphotVisualScaleImage (kapa, readout->variance, readout->mask, "variance", 1);
     189    psphotVisualScaleImage (kapa, readout->image, readout->mask, "image", 0);
    176190
    177191    pmVisualAskUser(NULL);
    … …  
    197211    }
    198212
    199     psphotVisualScaleImage (kapa, backgnd->image, "backgnd", 2);
    200     psphotVisualScaleImage (kapa, readout->image, "backsub", 0);
     213    psphotVisualScaleImage (kapa, backgnd->image, readout->mask, "backgnd", 2);
     214    psphotVisualScaleImage (kapa, readout->image, readout->mask, "backsub", 0);
    201215
    202216    pmVisualAskUser(NULL);
    … …  
    204218}
    205219
    206 bool psphotVisualShowSignificance (psImage *image) {
     220bool psphotVisualShowSignificance (psImage *image, float min, float max) {
    207221
    208222    if (!pmVisualIsVisual()) return true;
    … …  
    211225    if (kapa == -1) return false;
    212226
    213     psphotVisualRangeImage (kapa, image, "signif", 2, -1.0, 25.0*25.0);
     227    psphotVisualRangeImage (kapa, image, "signif", 2, min, max);
    214228
    215229    pmVisualAskUser(NULL);
    … …  
    247261        overlay[Noverlay].text = NULL;
    248262        Noverlay ++;
    249 
    250 # if (0)
    251         overlay[Noverlay].type = KII_OVERLAY_BOX;
    252         overlay[Noverlay].x = peak->x;
    253         overlay[Noverlay].y = peak->y;
    254         overlay[Noverlay].dx = 1.0;
    255         overlay[Noverlay].dy = 1.0;
    256         overlay[Noverlay].angle = 0.0;
    257         overlay[Noverlay].text = NULL;
    258         Noverlay ++;
    259 
    260         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    261         overlay[Noverlay].x = peak->xf;
    262         overlay[Noverlay].y = peak->yf;
    263         overlay[Noverlay].dx = 2.0;
    264         overlay[Noverlay].dy = 2.0;
    265         overlay[Noverlay].angle = 0.0;
    266         overlay[Noverlay].text = NULL;
    267         Noverlay ++;
    268 # endif
    269     }
    270 
    271 # if (1)
    272     overlay[Noverlay].type = KII_OVERLAY_BOX;
    273     overlay[Noverlay].x = 10.0;
    274     overlay[Noverlay].y = 10.0;
    275     overlay[Noverlay].dx = 1.0;
    276     overlay[Noverlay].dy = 1.0;
    277     overlay[Noverlay].angle = 0.0;
    278     overlay[Noverlay].text = NULL;
    279     Noverlay ++;
    280 # endif
     263    }
    281264
    282265    KiiLoadOverlay (kapa, overlay, Noverlay, "red");
    … …  
    463446    KapaSetFont (myKapa, "courier", 14);
    464447
     448    section.bg = KapaColorByName ("none"); // XXX probably should be 'none'
     449
    465450    float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
    466451
    … …  
    494479    Xmin = PS_MAX(Xmin, -0.1);
    495480    Ymin = PS_MAX(Ymin, -0.1);
     481
     482    // XXX test: hardwire plot limits
     483    // Xmin = -0.1; Ymin = -0.1;
     484    // Xmax = 20.1; Ymax = 20.1;
    496485
    497486    // storage vectors for data to be plotted
    … …  
    551540    KapaSetLimits (myKapa, &graphdata);
    552541
     542    graphdata.padXm = NAN;
     543    graphdata.padYm = NAN;
     544    graphdata.padXp = 0.5;
     545    graphdata.padYp = 0.5;
    553546    KapaBox (myKapa, &graphdata);
     547
    554548    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
    555549    KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM);
    … …  
    560554    graphdata.style = 2;
    561555    KapaPrepPlot (myKapa, nF, &graphdata);
     556
    562557    KapaPlotVector (myKapa, nF, xFaint->data.F32, "x");
    563558    KapaPlotVector (myKapa, nF, yFaint->data.F32, "y");
    … …  
    575570    section.dy = 0.25;
    576571    section.x  = 0.00;
    577     section.y  = 0.80;
     572    section.y  = 0.75;
    578573    section.name = psStringCopy ("MagMyy");
    579574    KapaSetSection (myKapa, &section);
    … …  
    587582    KapaSetLimits (myKapa, &graphdata);
    588583
     584    graphdata.padXm = 0.5;
     585    graphdata.padYm = NAN;
     586    graphdata.padXp = NAN;
     587    graphdata.padYp = 0.5;
    589588    strcpy (graphdata.labels, "0210");
    590589    KapaBox (myKapa, &graphdata);
    … …  
    611610    section.dx = 0.25;
    612611    section.dy = 0.75;
    613     section.x  = 0.80;
     612    section.x  = 0.75;
    614613    section.y  = 0.00;
    615614    section.name = psStringCopy ("MagMxx");
    … …  
    624623    KapaSetLimits (myKapa, &graphdata);
    625624
     625    graphdata.padXm = NAN;
     626    graphdata.padYm = 0.5;
     627    graphdata.padXp = 0.5;
     628    graphdata.padYp = NAN;
    626629    strcpy (graphdata.labels, "2001");
    627630    KapaBox (myKapa, &graphdata);
    … …  
    654657
    655658        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
    656         float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, analysis, "PSF_CLUMP_NSIGMA");
     659        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
    657660
    658661        graphdata.color = KapaColorByName ("blue");
    … …  
    11071110    }
    11081111
    1109     psphotVisualScaleImage (myKapa, outsat, "satstar", 2);
     1112    psphotVisualScaleImage (myKapa, outsat, NULL, "satstar", 2);
    11101113
    11111114    pmVisualAskUser(NULL);
    11121115    psFree (outsat);
    11131116    return true;
     1117}
     1118
     1119static void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1)
     1120{
     1121    float x[2], y[2];
     1122    x[0] = x0;
     1123    x[1] = x1;
     1124    y[0] = y0;
     1125    y[1] = y1;
     1126    KapaPrepPlot   (myKapa, 2, graphdata);
     1127    KapaPlotVector (myKapa, 2, x, "x");
     1128    KapaPlotVector (myKapa, 2, y, "y");
    11141129}
    11151130
    … …  
    11511166    }
    11521167
     1168    // generate model profiles (major and minor axis):
     1169    // create a model with theta = 0.0 so major and minor axes are equiv to x and y:
     1170    psEllipseShape rawShape, rotShape;
     1171
     1172    rawShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
     1173    rawShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
     1174    rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
     1175
     1176    psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0);
     1177
     1178    axes.theta = 0.0;
     1179
     1180    rotShape = psEllipseAxesToShape (axes);
     1181
     1182    psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32);
     1183    for (int i = 0; i < source->modelPSF->params->n; i++) {
     1184        params->data.F32[i] = source->modelPSF->params->data.F32[i];
     1185    }
     1186    params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2;
     1187    params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2;
     1188    params->data.F32[PM_PAR_SXY] = rotShape.sxy;
     1189    params->data.F32[PM_PAR_XPOS] = 0.0;
     1190    params->data.F32[PM_PAR_YPOS] = 0.0;
     1191
     1192    psVector *rmod = psVectorAlloc(300, PS_TYPE_F32);
     1193    psVector *fmaj = psVectorAlloc(300, PS_TYPE_F32);
     1194    psVector *fmin = psVectorAlloc(300, PS_TYPE_F32);
     1195
     1196    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     1197
     1198    float r = 0.0;
     1199    for (int i = 0; i < rmod->n; i++) {
     1200        r = i*0.1;
     1201        rmod->data.F32[i] = r;
     1202
     1203        coord->data.F32[1] = r;
     1204        coord->data.F32[0] = 0.0;
     1205        fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
     1206
     1207        coord->data.F32[0] = r;
     1208        coord->data.F32[1] = 0.0;
     1209        fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
     1210    }
     1211    psFree (coord);
     1212    psFree (params);
     1213
     1214    float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
     1215    float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
     1216    if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR);
     1217
     1218    psEllipseMoments emoments;
     1219    emoments.x2 = source->moments->Mxx;
     1220    emoments.xy = source->moments->Mxy;
     1221    emoments.y2 = source->moments->Myy;
     1222    axes = psEllipseMomentsToAxes (emoments, 20.0);
     1223    float MOMENTS_MAJOR = 2.355*axes.major;
     1224    float MOMENTS_MINOR = 2.355*axes.minor;
     1225
     1226    float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
     1227
    11531228    // reset source Add/Sub state to recorded
    11541229    if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    … …  
    11871262    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
    11881263
     1264    graphdata.color = KapaColorByName ("blue");
     1265    graphdata.ptype = 0;
     1266    graphdata.size = 0.0;
     1267    graphdata.style = 0;
     1268    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1269    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1270    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
     1271    plotline (myKapa, &graphdata, 0.0, logHM, 30.0, logHM);
     1272    plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, 0.0, 0.5*FWHM_MINOR, 5.0);
     1273    graphdata.ltype = 1;
     1274    plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, 0.0, 0.5*MOMENTS_MINOR, 5.0);
     1275    graphdata.ltype = 0;
     1276       
     1277    graphdata.color = KapaColorByName ("green");
     1278    graphdata.ptype = 0;
     1279    graphdata.size = 0.0;
     1280    graphdata.style = 0;
     1281    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1282    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1283    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
     1284    plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, 0.0, 0.5*FWHM_MAJOR, 5.0);
     1285    graphdata.ltype = 1;
     1286    plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, 0.0, 0.5*MOMENTS_MAJOR, 5.0);
     1287    graphdata.ltype = 0;
     1288       
     1289    for (int i = 0; i < rmod->n; i++) {
     1290        rmod->data.F32[i] = log10(rmod->data.F32[i]);
     1291    }
     1292
    11891293    // ** loglog **
    11901294    KapaSelectSection (myKapa, "loglog");
    … …  
    11951299    graphdata.ymin = -0.05;
    11961300    graphdata.ymax = +5.05;
     1301    graphdata.color = KapaColorByName ("black");
    11971302    KapaSetLimits (myKapa, &graphdata);
    11981303
    … …  
    12171322    KapaPlotVector (myKapa, nb, Rb->data.F32, "x");
    12181323    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
     1324
     1325    graphdata.color = KapaColorByName ("blue");
     1326    graphdata.ptype = 0;
     1327    graphdata.size = 0.0;
     1328    graphdata.style = 0;
     1329    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1330    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1331    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
     1332
     1333    graphdata.color = KapaColorByName ("green");
     1334    graphdata.ptype = 0;
     1335    graphdata.size = 0.0;
     1336    graphdata.style = 0;
     1337    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
     1338    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
     1339    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
     1340
     1341    psFree (rmod);
     1342    psFree (fmin);
     1343    psFree (fmaj);
    12191344
    12201345    psFree (rg);
    … …  
    12401365    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    12411366    assert (maskVal);
     1367
     1368    section.bg  = KapaColorByName ("none"); // XXX probably should be 'none'
    12421369
    12431370    KapaClearPlots (myKapa);
    … …  
    14051532        if (source == NULL) continue;
    14061533
    1407         // if (source->type != type) continue;
    14081534        if (mode) {
    14091535            if (keep) {
    … …  
    14781604    KapaSetFont (myKapa, "courier", 14);
    14791605
     1606    section.bg  = KapaColorByName ("none"); // XXX probably should be 'none'
     1607
    14801608    // select the max psfX,Y values for the plot limits
    14811609    float Xmin = 1000.0, Xmax = 0.0;
    … …  
    15291657    psVector *sDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    15301658
     1659    psVector *xLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1660    psVector *yLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1661    psVector *mLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1662    psVector *sLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     1663
    15311664    psVector *xCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    15321665    psVector *yCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    … …  
    15391672    int nPSF = 0;
    15401673    int nDEF = 0;
     1674    int nLOW = 0;
    15411675    int nCR  = 0;
    15421676    for (int i = 0; i < sources->n; i++) {
    15431677        pmSource *source = sources->data[i];
    15441678        if (source->moments == NULL) continue;
     1679
     1680        // only plot the measured sources...
     1681        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) continue;
    15451682
    15461683        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
    … …  
    15741711            continue;
    15751712        }
    1576         if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
     1713        if (source->errMag > 0.1) {
     1714            xLOW->data.F32[nLOW] = source->moments->Mxx;
     1715            yLOW->data.F32[nLOW] = source->moments->Myy;
     1716            mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum);
     1717            sLOW->data.F32[nLOW] = source->extNsigma;
     1718            nLOW++;
    15771719            continue;
    15781720        }
    … …  
    15831725        nPSF++;
    15841726    }
     1727
    15851728    xSAT->n = nSAT;
    15861729    ySAT->n = nSAT;
    … …  
    16071750    mDEF->n = nDEF;
    16081751    sDEF->n = nDEF;
     1752
     1753    xLOW->n = nLOW;
     1754    yLOW->n = nLOW;
     1755    mLOW->n = nLOW;
     1756    sLOW->n = nLOW;
    16091757
    16101758    // four sections: MxxMyy, MagMxx, MagMyy, MagSigma
    … …  
    16261774    KapaSetLimits (myKapa, &graphdata);
    16271775
     1776    graphdata.padXm = NAN;
     1777    graphdata.padYm = NAN;
     1778    graphdata.padXp = 0.5;
     1779    graphdata.padYp = 0.5;
    16281780    KapaBox (myKapa, &graphdata);
    16291781    KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM);
    … …  
    16691821    KapaPlotVector (myKapa, nSAT, xSAT->data.F32, "x");
    16701822    KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y");
     1823
     1824    graphdata.color = KapaColorByName ("black");
     1825    graphdata.ptype = 7;
     1826    graphdata.size = 1.0;
     1827    graphdata.style = 2;
     1828    KapaPrepPlot   (myKapa, nLOW, &graphdata);
     1829    KapaPlotVector (myKapa, nLOW, xLOW->data.F32, "x");
     1830    KapaPlotVector (myKapa, nLOW, yLOW->data.F32, "y");
    16711831
    16721832    // second section: MagMyy
    … …  
    16811841    graphdata.color = KapaColorByName ("black");
    16821842    graphdata.xmin = -17.1;
    1683     graphdata.xmax =  -7.9;
     1843    graphdata.xmax =  -6.9;
    16841844    graphdata.ymin = Ymin;
    16851845    graphdata.ymax = Ymax;
    16861846    KapaSetLimits (myKapa, &graphdata);
    16871847
     1848    graphdata.padXm = 0.5;
     1849    graphdata.padYm = NAN;
     1850    graphdata.padXp = NAN;
     1851    graphdata.padYp = 0.5;
    16881852    strcpy (graphdata.labels, "0210");
    16891853    KapaBox (myKapa, &graphdata);
    … …  
    17341898    section.dx = 0.25;
    17351899    section.dy = 0.60;
    1736     section.x  = 0.80;
     1900    section.x  = 0.75;
    17371901    section.y  = 0.00;
    17381902    section.name = psStringCopy ("MagMxx");
    … …  
    17431907    graphdata.xmin = Xmin;
    17441908    graphdata.xmax = Xmax;
    1745     graphdata.ymin =  -7.9;
     1909    graphdata.ymin =  -6.9;
    17461910    graphdata.ymax = -17.1;
    17471911    KapaSetLimits (myKapa, &graphdata);
    17481912
     1913    graphdata.padXm = NAN;
     1914    graphdata.padYm = 0.5;
     1915    graphdata.padXp = 0.5;
     1916    graphdata.padYp = NAN;
    17491917    strcpy (graphdata.labels, "2001");
    17501918    KapaBox (myKapa, &graphdata);
    … …  
    17941962    // fourth section: MagSigma
    17951963    section.dx = 0.75;
    1796     section.dy = 0.15;
     1964    section.dy = 0.20;
    17971965    section.x  = 0.00;
    1798     section.y  = 0.65;
     1966    section.y  = 0.60;
    17991967    section.name = psStringCopy ("MagSigma");
    18001968    KapaSetSection (myKapa, &section);
    … …  
    18021970
    18031971    graphdata.color = KapaColorByName ("black");
    1804     graphdata.xmax =  -7.9;
     1972    graphdata.xmax =  -6.9;
    18051973    graphdata.xmin = -17.1;
    18061974    graphdata.ymin = -20.1;
    … …  
    18081976    KapaSetLimits (myKapa, &graphdata);
    18091977
     1978    graphdata.padXm = 0.5;
     1979    graphdata.padYm = NAN;
     1980    graphdata.padXp = 0.5;
     1981    graphdata.padYp = 0.5;
    18101982    strcpy (graphdata.labels, "0100");
    18111983    KapaBox (myKapa, &graphdata);
    … …  
    18622034
    18632035        int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS");
    1864         float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, analysis, "PSF_CLUMP_NSIGMA");
     2036        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
    18652037
    18662038        graphdata.color = KapaColorByName ("blue");
    … …  
    19012073    psFree (ySAT);
    19022074    psFree (mSAT);
     2075    psFree (sSAT);
    19032076
    19042077    psFree (xEXT);
    19052078    psFree (yEXT);
    19062079    psFree (mEXT);
     2080    psFree (sEXT);
    19072081
    19082082    psFree (xPSF);
    19092083    psFree (yPSF);
    19102084    psFree (mPSF);
     2085    psFree (sPSF);
    19112086
    19122087    psFree (xDEF);
    19132088    psFree (yDEF);
    19142089    psFree (mDEF);
     2090    psFree (sDEF);
     2091
     2092    psFree (xLOW);
     2093    psFree (yLOW);
     2094    psFree (mLOW);
     2095    psFree (sLOW);
    19152096
    19162097    psFree (xCR);
    19172098    psFree (yCR);
    19182099    psFree (mCR);
     2100    psFree (sCR);
    19192101
    19202102    pmVisualAskUser(NULL);
    … …  
    19292111    if (myKapa == -1) return false;
    19302112
    1931     psphotVisualScaleImage (myKapa, readout->image, "resid", 1);
     2113    psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", 1);
    19322114
    19332115    pmVisualAskUser(NULL);
    … …  
    20542236}
    20552237
     2238bool psphotVisualPlotChisq (psArray *sources) {
     2239
     2240    Graphdata graphdata;
     2241
     2242    if (!pmVisualIsVisual()) return true;
     2243
     2244    int myKapa = psphotKapaChannel (2);
     2245    if (myKapa == -1) return false;
     2246
     2247    KapaClearPlots (myKapa);
     2248    KapaInitGraph (&graphdata);
     2249
     2250    psVector *x = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2251    psVector *y = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     2252
     2253    graphdata.xmin = +32.0;
     2254    graphdata.xmax = -32.0;
     2255    graphdata.ymin = +32.0;
     2256    graphdata.ymax = -32.0;
     2257
     2258    FILE *f = fopen ("chisq.dat", "w");
     2259
     2260    // construct the plot vectors
     2261    int n = 0;
     2262    for (int i = 0; i < sources->n; i++) {
     2263        pmSource *source = sources->data[i];
     2264        if (!source) continue;
     2265        if (source->type != PM_SOURCE_TYPE_STAR) continue;
     2266        if (!source->moments) continue;
     2267        if (!isfinite(source->moments->Sum)) continue;
     2268        if (!source->modelPSF) continue;
     2269        if (!isfinite(source->modelPSF->chisq)) continue;
     2270
     2271        x->data.F32[n] = -2.5*log10(source->moments->Sum);
     2272        y->data.F32[n] = source->modelPSF->chisq / source->modelPSF->nDOF;
     2273        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
     2274        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
     2275        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
     2276        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);
     2277
     2278        fprintf (f, "%d %d %f %f\n", i, n, x->data.F32[n], y->data.F32[n]);
     2279
     2280        n++;
     2281    }
     2282    x->n = y->n = n;
     2283    fclose (f);
     2284
     2285    float range;
     2286    range = graphdata.xmax - graphdata.xmin;
     2287    graphdata.xmax += 0.05*range;
     2288    graphdata.xmin -= 0.05*range;
     2289    range = graphdata.ymax - graphdata.ymin;
     2290    graphdata.ymax += 0.05*range;
     2291    graphdata.ymin -= 0.05*range;
     2292
     2293    // XXX test
     2294    graphdata.xmin = -17.0;
     2295    graphdata.xmax =  -3.0;
     2296    graphdata.ymin =  -0.1;
     2297    graphdata.ymax = +10.1;
     2298
     2299    KapaSetLimits (myKapa, &graphdata);
     2300
     2301    KapaSetFont (myKapa, "helvetica", 14);
     2302    KapaBox (myKapa, &graphdata);
     2303    KapaSendLabel (myKapa, "PSF Mag", KAPA_LABEL_XM);
     2304    KapaSendLabel (myKapa, "ChiSq", KAPA_LABEL_YM);
     2305
     2306    graphdata.color = KapaColorByName ("black");
     2307    graphdata.ptype = 2;
     2308    graphdata.size = 0.5;
     2309    graphdata.style = 2;
     2310    KapaPrepPlot (myKapa, n, &graphdata);
     2311    KapaPlotVector (myKapa, n, x->data.F32, "x");
     2312    KapaPlotVector (myKapa, n, y->data.F32, "y");
     2313
     2314    psFree (x);
     2315    psFree (y);
     2316
     2317    pmVisualAskUser(NULL);
     2318    return true;
     2319}
     2320
    20562321bool psphotVisualShowPetrosians (psArray *sources) {
    20572322
    … …  
    20732338        if (!source) continue;
    20742339        if (!source->extpars) continue;
    2075         if (!source->extpars->profile) continue;
    2076         if (!source->extpars->petrosian_80) continue;
    2077 
    2078         pmSourceRadialProfile *profile = source->extpars->profile;
    2079         pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80;
     2340        if (!source->extpars->petProfile) continue;
     2341
     2342        float petrosianRadius = source->extpars->petrosianRadius;
     2343        psEllipseAxes *axes = &source->extpars->axes;
    20802344
    20812345        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    20822346        overlay[Noverlay].x = source->peak->xf;
    20832347        overlay[Noverlay].y = source->peak->yf;
    2084         overlay[Noverlay].dx = 2.0*petrosian->radius;
    2085         overlay[Noverlay].dy = 2.0*petrosian->radius*profile->axes.minor/profile->axes.major;
    2086         overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
     2348        overlay[Noverlay].dx = 1.0*petrosianRadius;
     2349        overlay[Noverlay].dy = 1.0*petrosianRadius*axes->minor/axes->major;
     2350        overlay[Noverlay].angle = axes->theta * PS_DEG_RAD;
    20872351        overlay[Noverlay].text = NULL;
    20882352        Noverlay ++;
    20892353        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    20902354
    2091         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    2092         overlay[Noverlay].x = source->peak->xf;
    2093         overlay[Noverlay].y = source->peak->yf;
    2094         overlay[Noverlay].dx = 4.0*petrosian->radius;
    2095         overlay[Noverlay].dy = 4.0*petrosian->radius*profile->axes.minor/profile->axes.major;
    2096         overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD;
    2097         overlay[Noverlay].text = NULL;
    2098         Noverlay ++;
    2099         CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     2355        // overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     2356        // overlay[Noverlay].x = source->peak->xf;
     2357        // overlay[Noverlay].y = source->peak->yf;
     2358        // overlay[Noverlay].dx = 2.0*petrosianRadius;
     2359        // overlay[Noverlay].dy = 2.0*petrosianRadius*axes->minor/axes->major;
     2360        // overlay[Noverlay].angle = axes->theta * PS_DEG_RAD;
     2361        // overlay[Noverlay].text = NULL;
     2362        // Noverlay ++;
     2363        // CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    21002364    }
    21012365
Note: See TracChangeset for help on using the changeset viewer.