IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25496


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

Location:
branches/eam_branches/20090715/psModules
Files:
19 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psModules/src/extras/pmVisual.c

    r25476 r25496  
    2424#include "pmPSF.h"
    2525#include "pmPSFtry.h"
     26#include "pmSource.h"
    2627#include "pmFPAExtent.h"
    2728

  • branches/eam_branches/20090715/psModules/src/imcombine/pmPSFEnvelope.c

    r25407 r25496  
    320320    options->poissonErrorsParams = true;
    321321    options->stats = psStatsAlloc(PSF_STATS);
    322     options->radius = maxRadius;
     322    options->fitRadius = maxRadius;
     323    options->apRadius = maxRadius; // XXX need to decide if aperture mags need a different radius
    323324    options->psfTrendMode = PM_TREND_MAP;
    324325    options->psfTrendNx = xOrder;

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

    r25355 r25496  
    6262    tmp->nPix  = 0;
    6363    tmp->nIter = 0;
    64     tmp->radiusFit = 0;
     64    tmp->fitRadius = 0;
    6565    tmp->flags = PM_MODEL_STATUS_NONE;
    6666    tmp->residuals = NULL;              // XXX should the model own this memory?
    … …  
    109109    new->nIter     = model->nIter;
    110110    new->flags     = model->flags;
    111     new->radiusFit = model->radiusFit;
     111    new->fitRadius = model->fitRadius;
    112112
    113113    for (int i = 0; i < new->params->n; i++) {

  • branches/eam_branches/20090715/psModules/src/objects/pmModel.h

    r21516 r25496  
    9696    int nIter;                          ///< number of iterations to reach min
    9797    pmModelStatus flags;                ///< model status flags
    98     float radiusFit;                    ///< fit radius actually used
     98    float fitRadius;                    ///< fit radius actually used
    9999    pmResiduals *residuals;             ///< normalized PSF residuals
    100100

  • branches/eam_branches/20090715/psModules/src/objects/pmPSF.h

    r25476 r25496  
    7979    bool          poissonErrorsPhotLin; ///< use poission errors for linear model fitting
    8080    bool          poissonErrorsParams; ///< use poission errors for model parameter fitting
    81     float         radius;
     81    float         fitRadius;
     82    float         apRadius;
    8283    bool          chiFluxTrend;         // Fit a trend in Chi2 as a function of flux?
    8384} pmPSFOptions;

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

    r25476 r25496  
    7070        // set object mask to define valid pixels
    7171        // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)?
    72         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);
     72        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
    7373
    7474        // fit model as EXT, not PSF
    … …  
    7676
    7777        // clear object mask to define valid pixels
    78         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
     78        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal));
    7979
    8080        // exclude the poor fits

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

    r25476 r25496  
    6767            continue;
    6868        }
    69         source->modelPSF->radiusFit = options->radius;
    70         // XXXX use a different radius for the aperture magnitude than for the PSF fit?
     69        // PSF fit and aperture mags use different radii
     70        source->modelPSF->fitRadius = options->fitRadius;
     71        source->apRadius = options->apRadius;
    7172
    72         // set object mask to define valid pixels
    73         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);
     73        // set object mask to define valid pixels for PSF model fit
     74        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
    7475
    7576        // fit the PSF model to the source
    … …  
    7879        // skip poor fits
    7980        if (!status) {
    80             psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
     81            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal));
    8182            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
    8283            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
    … …  
    8485        }
    8586
    86         // XXX this function calculates the aperture magnitude, but we now use the moments->Sum as the flux
     87        // set object mask to define valid pixels for APERTURE magnitude
     88        if (options->fitRadius < options->apRadius) {
     89            // this is not really a recommended situation...
     90            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal));
     91            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal);
     92        }
     93        if (options->fitRadius > options->apRadius) {
     94            // ensure the aperture is defined
     95            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal);
     96        }
     97
     98        // This function calculates the psf and aperture magnitudes
    8799        status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal);
    88100        if (!status || isnan(source->apMag)) {
    89             psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
     101            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "AND", PS_NOT_IMAGE_MASK(markVal));
    90102            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
    91103            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
    … …  
    94106
    95107        // clear object mask to define valid pixels
    96         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
     108        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "AND", PS_NOT_IMAGE_MASK(markVal));
    97109
    98110        psfTry->fitMag->data.F32[i] = source->psfMag;
    99         psfTry->metric->data.F32[i] = -2.5*log10(source->moments->Sum) - source->psfMag;
     111        // psfTry->metric->data.F32[i] = -2.5*log10(source->moments->Sum) - source->psfMag;
     112        psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
    100113        psfTry->metricErr->data.F32[i] = source->errMag;
    101114
    … …  
    104117    }
    105118    psfTry->psf->nPSFstars = Npsf;
     119
     120    // XXX test code:
     121    pmSourceVisualShowModelFits (psfTry->psf, psfTry->sources, maskVal);
    106122
    107123    psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, psfTry->sources->n);

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

    r25476 r25496  
    184184    // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    185185    // This way, the parameters masked by one of the fits will be applied to the others
    186     for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
     186    // NOTE : trend->stats (below) points to the same data as psfTrendStats; set the value to 1
     187    // to limit the iteration within the loop
     188    int nIter = psf->psfTrendStats->clipIter;
     189    psf->psfTrendStats->clipIter = 1;
     190    for (int i = 0; i < nIter; i++) {
    187191        // XXX we are using the same stats structure on each pass: do we need to re-init it?
    188192        // XXX we hardwire this to SAMPLE stats above (psphotChoosePSF.c), hardwire here instead?
    … …  
    203207        psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n);
    204208        if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map);
    205         // pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask);
     209        pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask);
    206210
    207211        trend = psf->params->data[PM_PAR_E1];
    … …  
    212216        psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n);
    213217        if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map);
    214         // pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask);
     218        pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask);
    215219
    216220        trend = psf->params->data[PM_PAR_E2];
    … …  
    221225        psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n);
    222226        if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map);
    223         // pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask);
     227        pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask);
    224228
    225229        if (!status) {
    … …  
    228232        }
    229233    }
     234    psf->psfTrendStats->clipIter = nIter;
    230235
    231236    // test dump of the psf parameters

  • branches/eam_branches/20090715/psModules/src/objects/pmSource.h

    r25385 r25496  
    8282    float extNsigma;                    ///< Nsigma deviation from PSF to EXT
    8383    float sky, skyErr;                  ///< The sky and its error at the center of the object
     84    float apRadius;
    8485    psRegion region;                    ///< area on image covered by selected pixels
    8586    pmSourceExtendedPars *extpars;      ///< extended source parameters

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

    r25451 r25496  
    128128            nDOF = model->nDOF;
    129129            nPix = model->nPix;
    130             apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius?
     130            apRadius = source->apRadius;
    131131            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    132132        } else {
    … …  
    308308        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    309309        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
     310        source->apRadius  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");
    310311
    311312        // note that some older versions used PSF_PROBABILITY: this was not well defined.
    … …  
    313314        model->nDOF       = psMetadataLookupS32 (&status, row, "PSF_NDOF");
    314315        model->nPix       = psMetadataLookupS32 (&status, row, "PSF_NPIX");
    315         model->radiusFit  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");
    316316
    317317        source->moments = pmMomentsAlloc ();

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

    r25452 r25496  
    128128            nDOF = model->nDOF;
    129129            nPix = model->nPix;
    130             apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius?
     130            apRadius = source->apRadius;
    131131            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    132132        } else {
    … …  
    308308        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    309309        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
     310        source->apRadius  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");
    310311
    311312        // note that some older versions used PSF_PROBABILITY: this was not well defined.
    … …  
    313314        model->nDOF       = psMetadataLookupS32 (&status, row, "PSF_NDOF");
    314315        model->nPix       = psMetadataLookupS32 (&status, row, "PSF_NPIX");
    315         model->radiusFit  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");
    316316
    317317        source->moments = pmMomentsAlloc ();

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

    r24581 r25496  
    119119                 logChi, logChiNorm,
    120120                 source[0].peak->SN,
    121                  model[0].radiusFit,
     121                 source[0].apRadius,
    122122                 source[0].pixWeight,
    123123                 model[0].nDOF,
    … …  
    181181                 log10(model[0].chisqNorm/model[0].nDOF),
    182182                 source[0].peak->SN,
    183                  model[0].radiusFit,
     183                 source[0].apRadius,
    184184                 source[0].pixWeight,
    185185                 model[0].nDOF,

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

    r21511 r25496  
    201201    if (isfinite (source->apMag) && isPSF && psf) {
    202202        if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
    203             source->apMag += pmGrowthCurveCorrect (psf->growth, model->radiusFit);
     203            source->apMag += pmGrowthCurveCorrect (psf->growth, source->apRadius);
    204204        }
    205205        if (mode & PM_SOURCE_PHOT_APCORR) {
     206            // XXX this should be removed -- we no longer fit for the 'sky bias'
    206207            rflux   = pow (10.0, 0.4*source->psfMag);
    207             source->apMag -= PS_SQR(model->radiusFit)*rflux * psf->skyBias + psf->skySat / rflux;
     208            source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux;
    208209        }
    209210    }

  • 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:

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

    r25476 r25496  
    1919bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask);
    2020bool pmSourceVisualPlotPSFMetric (pmPSFtry *try);
     21bool pmSourceVisualShowModelFit (pmSource *source);
     22bool pmSourceVisualShowModelFits (pmPSF *psf, psArray *sources, psImageMaskType maskVal);
    2123
    2224/// @}

  • branches/eam_branches/20090715/psModules/test/objects/tap_pmGrowthCurve.c

    r25022 r25496  
    131131        source->mode = PM_SOURCE_MODE_PSFSTAR;
    132132
    133         source->modelPSF->radiusFit = 15.0;
     133        source->apRadius = 15.0;
    134134
    135135        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    136136        double refMag = source->apMag;
    137137
    138         source->modelPSF->radiusFit = 10.0;
    139         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    140         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    141 
    142         source->modelPSF->radiusFit = 8.0;
    143         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    144         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    145 
    146         source->modelPSF->radiusFit = 6.0;
     138        source->apRadius = 10.0;
     139        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     140        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     141
     142        source->apRadius = 8.0;
     143        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     144        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     145
     146        source->apRadius = 6.0;
    147147        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    148148        ok_float_tol(refMag - source->apMag, +0.0003, 0.0001, "growth offset is is %f", refMag - source->apMag);
    149149
    150         source->modelPSF->radiusFit = 4.0;
     150        source->apRadius = 4.0;
    151151        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    152152        ok_float_tol(refMag - source->apMag, +0.0020, 0.0001, "growth offset is is %f", refMag - source->apMag);
    153153
    154         source->modelPSF->radiusFit = 3.0;
     154        source->apRadius = 3.0;
    155155        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    156156        ok_float_tol(refMag - source->apMag, -0.0001, 0.0001, "growth offset is is %f", refMag - source->apMag);
    157157
    158         source->modelPSF->radiusFit = 2.0;
     158        source->apRadius = 2.0;
    159159        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    160160        ok_float_tol(refMag - source->apMag, -0.0075, 0.0001, "growth offset is is %f", refMag - source->apMag);
    … …  
    234234        source->mode = PM_SOURCE_MODE_PSFSTAR;
    235235
    236         source->modelPSF->radiusFit = 15.0;
     236        source->apRadius = 15.0;
    237237        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    238238        double refMag = source->apMag;
    239239
    240         source->modelPSF->radiusFit = 10.0;
    241         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    242         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    243 
    244         source->modelPSF->radiusFit = 8.0;
    245         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    246         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    247 
    248         source->modelPSF->radiusFit = 6.0;
     240        source->apRadius = 10.0;
     241        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     242        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     243
     244        source->apRadius = 8.0;
     245        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     246        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     247
     248        source->apRadius = 6.0;
    249249        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    250250        ok_float_tol(refMag - source->apMag, +0.0004, 0.0001, "growth offset is is %f", refMag - source->apMag);
    251251
    252         source->modelPSF->radiusFit = 4.0;
     252        source->apRadius = 4.0;
    253253        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    254254        ok_float_tol(refMag - source->apMag, +0.0026, 0.0001, "growth offset is is %f", refMag - source->apMag);
    255255
    256         source->modelPSF->radiusFit = 3.0;
     256        source->apRadius = 3.0;
    257257        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    258258        ok_float_tol(refMag - source->apMag, -0.0001, 0.0001, "growth offset is is %f", refMag - source->apMag);
    259259
    260         source->modelPSF->radiusFit = 2.0;
     260        source->apRadius = 2.0;
    261261        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    262262        ok_float_tol(refMag - source->apMag, -0.0103, 0.0001, "growth offset is is %f", refMag - source->apMag);
    … …  
    336336        source->mode = PM_SOURCE_MODE_PSFSTAR;
    337337
    338         source->modelPSF->radiusFit = 15.0;
     338        source->apRadius = 15.0;
    339339        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    340340        double refMag = source->apMag;
    341341
    342         source->modelPSF->radiusFit = 10.0;
    343         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    344         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    345 
    346         source->modelPSF->radiusFit = 8.0;
    347         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    348         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    349 
    350         source->modelPSF->radiusFit = 6.0;
     342        source->apRadius = 10.0;
     343        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     344        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     345
     346        source->apRadius = 8.0;
     347        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     348        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     349
     350        source->apRadius = 6.0;
    351351        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    352352        ok_float_tol(refMag - source->apMag, +0.0006, 0.0001, "growth offset is is %f", refMag - source->apMag);
    353353
    354         source->modelPSF->radiusFit = 4.0;
     354        source->apRadius = 4.0;
    355355        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    356356        ok_float_tol(refMag - source->apMag, +0.0038, 0.0001, "growth offset is is %f", refMag - source->apMag);
    357357
    358         source->modelPSF->radiusFit = 3.0;
    359         pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    360         ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
    361 
    362         source->modelPSF->radiusFit = 2.0;
     358        source->apRadius = 3.0;
     359        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
     360        ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag);
     361
     362        source->apRadius = 2.0;
    363363        pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0);
    364364        ok_float_tol(refMag - source->apMag, -0.0164, 0.0001, "growth offset is is %f", refMag - source->apMag);

  • branches/eam_branches/20090715/psModules/test/objects/tap_pmModel.c

    r21471 r25496  
    7777        ok(model->nDOF == 0, "pmModelAlloc() set pmModel->nDOF correctly");
    7878        ok(model->nIter == 0, "pmModelAlloc() set pmModel->nIter correctly");
    79         ok(model->radiusFit == 0, "pmModelAlloc() set pmModel->radiusFit correctly");
     79        ok(model->fitRadius == 0, "pmModelAlloc() set pmModel->fitRadius correctly");
    8080        ok(model->flags == PM_MODEL_STATUS_NONE, "pmModelAlloc() set pmModel->flags correctly");
    8181        ok(model->residuals == NULL, "pmModelAlloc() set pmModel->residuals correctly");
    … …  
    132132        modelSrc->nIter = 3;
    133133        modelSrc->flags = PM_MODEL_STATUS_NONE;
    134         modelSrc->radiusFit = 4;
     134        modelSrc->fitRadius = 4;
    135135        pmModel *modelDst = pmModelCopy(modelSrc);
    136136        ok(modelDst != NULL && psMemCheckModel(modelDst), "pmModelCopy() returned a non-NULL pmModel");
    … …  
    139139        ok(modelDst->nIter == 3, "pmModelCopy() set the pmModel->nIter member correctly");
    140140        ok(modelDst->flags == PM_MODEL_STATUS_NONE, "pmModelCopy() set the pmModel->flags member correctly");
    141         ok(modelDst->radiusFit == 4, "pmModelCopy() set the pmModel->radiusFit member correctly");
     141        ok(modelDst->fitRadius == 4, "pmModelCopy() set the pmModel->fitRadius member correctly");
    142142
    143143        psFree(modelSrc);

  • branches/eam_branches/20090715/psModules/test/objects/tap_pmModelUtils.c

    r15985 r25496  
    8181        ok(tmpModel->nIter == testModelPSF->nIter, "pmModelFromPSF() set the model->nIter correctly");
    8282        ok(tmpModel->flags == testModelPSF->flags, "pmModelFromPSF() set the model->flags correctly");
    83         ok(tmpModel->radiusFit == testModelPSF->radiusFit, "pmModelFromPSF() set the model->radiusFit correctly");
     83        ok(tmpModel->fitRadius == testModelPSF->fitRadius, "pmModelFromPSF() set the model->fitRadius correctly");
    8484        ok(tmpModel->modelFunc == testModelPSF->modelFunc, "pmModelFromPSF() set the model->modelFunc correctly");
    8585        ok(tmpModel->modelFlux == testModelPSF->modelFlux, "pmModelFromPSF() set the model->modelFlux correctly");
    … …  
    140140        ok(tmpModel->nIter == testModelPSF->nIter, "pmModelFromPSF() set the model->nIter correctly");
    141141        ok(tmpModel->flags == testModelPSF->flags, "pmModelFromPSF() set the model->flags correctly");
    142         ok(tmpModel->radiusFit == testModelPSF->radiusFit, "pmModelFromPSF() set the model->radiusFit correctly");
     142        ok(tmpModel->fitRadius == testModelPSF->fitRadius, "pmModelFromPSF() set the model->fitRadius correctly");
    143143        ok(tmpModel->modelFunc == testModelPSF->modelFunc, "pmModelFromPSF() set the model->modelFunc correctly");
    144144        ok(tmpModel->modelFlux == testModelPSF->modelFlux, "pmModelFromPSF() set the model->modelFlux correctly");

  • branches/eam_branches/20090715/psModules/test/objects/tap_pmSourcePhotometry.c

    r9922 r25496  
    9696    source->modelPSF = pmModelFromPSF (modelRef, psf);
    9797    source->modelPSF->dparams->data.F32[PM_PAR_I0] = 1;
    98     source->modelPSF->radiusFit = radius;
     98    source->apRadius = radius;
    9999
    100100    // measure photometry for centered source (fractional pix : 0.5,0.5)
Note: See TracChangeset for help on using the changeset viewer.