IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5980


Ignore:
Timestamp:
Jan 13, 2006, 8:24:10 PM (20 years ago)
Author:
eugene
Message:

substantial work on the aperture residuals

Location:
trunk/psphot
Files:
1 added
15 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/Makefile

    r5840 r5980  
    5858$(SRC)/psphotBlendFit.$(ARCH).o        \
    5959$(SRC)/psphotApResid.$(ARCH).o \
     60$(SRC)/psphotAssessPSF.$(ARCH).o \
    6061$(SRC)/psBicube.$(ARCH).o
    6162
    … …  
    8283$(SRC)/psphotTest.$(ARCH).o \
    8384$(SRC)/psphotTestArguments.$(ARCH).o \
    84 $(SRC)/pmSourceContour.$(ARCH).o \
    85 $(SRC)/psImageData.$(ARCH).o        \
    86 $(SRC)/psphotSetup.$(ARCH).o        \
    87 $(SRC)/psModulesUtils.$(ARCH).o     \
    88 $(SRC)/psSparse.$(ARCH).o
     85$(SRC)/psphotApResid.$(ARCH).o
    8986
    9087psphot: $(BIN)/psphot.$(ARCH)
    … …  
    10299$(TEST): $(SRC)/psphot.h
    103100
    104 INSTALL = psphot psphotTest psphotModelTest
     101INSTALL = psphot
    105102
    106103# dependancy rules for binary code #########################

  • trunk/psphot/doc/notes.txt

    r5840 r5980  
    1 
    2 Notes on psphot
    31
    422005.12.23

  • trunk/psphot/src/psLine.c

    r4977 r5980  
    4141    Nchar = vsnprintf (&line->line[line->Nline], nMax, format, ap);
    4242    line->Nline += PS_MIN (nMax - 1, Nchar);
     43    va_end (ap);
    4344
    4445    if (Nchar >= nMax) return (false);

  • trunk/psphot/src/psphot.c

    r5828 r5980  
    7878        psphotFullFit (imdata, config, sources, psf, sky);
    7979        psphotReplaceUnfit (sources);
    80         psphotApResid (sources, config, psf);
     80        psphotApResid (imdata, sources, config, psf);
    8181        break;
    8282
    … …  
    8585        psphotBlendFit (imdata, config, sources, psf, sky);
    8686        psphotReplaceUnfit (sources);
    87         psphotApResid (sources, config, psf);
     87        psphotApResid (imdata, sources, config, psf);
    8888        break;
    8989

  • trunk/psphot/src/psphot.h

    r5828 r5980  
    44# include <pslib.h>
    55# include <pmObjects.h>
     6# include <pmGrowthCurve.h>
    67# include <pmPSF.h>
    78# include <pmPSFtry.h>
    … …  
    3738int             psphotSaveImage (psMetadata *header, psImage *image, char *filename);
    3839bool            psphotDefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius);
     40bool            psphotRedefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius);
    3941void            psphotModelGroupInit (void);
    4042
    … …  
    8486bool psphotReplaceUnfit (psArray *sources);
    8587bool psphotDumpMoments (psMetadata *config, psArray *sources);
    86 bool psphotApResid (psArray *sources, psMetadata *config, pmPSF *psf);
     88bool psphotApResid (eamReadout *imdata, psArray *sources, psMetadata *config, pmPSF *psf);
     89bool psphotAssessPSF (eamReadout *imdata, psMetadata *config, pmPSF *psf);
     90
    8791bool psphotWritePSF (pmPSF *psf, char *filename);
    8892pmPSF *psphotReadPSF (char *filename);
    … …  
    9094bool psPolynomial2DtoMD (psMetadata *md, psPolynomial2D *poly, char *format, ...);
    9195bool psPolynomial3DtoMD (psMetadata *md, psPolynomial3D *poly, char *format, ...);
     96bool psPolynomial4DtoMD (psMetadata *md, psPolynomial4D *poly, char *format, ...);
    9297psPolynomial2D *psPolynomial2DfromMD (psMetadata *folder);
    9398psPolynomial3D *psPolynomial3DfromMD (psMetadata *folder);
     99psPolynomial4D *psPolynomial4DfromMD (psMetadata *folder);
    94100bool psphotModelTest (eamReadout *imdata, psMetadata *config);
    95101
    … …  
    112118psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y);
    113119psPlane psImageBicubeMin (psPolynomial2D *poly);
     120void psphotTestArguments (int *argc, char **argv);
     121
     122bool psphotGrowthCurve (eamReadout *imdata, pmPSF *psf);
     123
     124psPolynomial4D *psVectorChiClipFitPolynomial4D(
     125    psPolynomial4D *poly,
     126    psStats *stats,
     127    const psVector *mask,
     128    psMaskType maskValue,
     129    const psVector *f,
     130    const psVector *fErr,
     131    const psVector *x,
     132    const psVector *y,
     133    const psVector *z,
     134    const psVector *t);
     135
     136// XXX deprecated
     137// bool psphotApResidFullFit (psVector *x, psVector *y, psVector *radius, psVector *rflux, psVector *apresid, psVector *dmag);
     138// bool psphotApResidPolyFit (psVector *x, psVector *y, psVector *radius, psVector *rflux, psVector *apresid, psVector *dmag);

  • trunk/psphot/src/psphotApResid.c

    r5952 r5980  
    11# include "psphot.h"
    22
    3 bool psphotApResid (psArray *sources, psMetadata *config, pmPSF *psf) {
     3psPolynomial4D *psVectorChiClipFitPolynomial4D(
     4    psPolynomial4D *poly,
     5    psStats *stats,
     6    const psVector *mask,
     7    psMaskType maskValue,
     8    const psVector *f,
     9    const psVector *fErr,
     10    const psVector *x,
     11    const psVector *y,
     12    const psVector *z,
     13    const psVector *t)
     14{
     15    PS_ASSERT_POLY_NON_NULL(poly, NULL);
     16    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
     17    PS_ASSERT_PTR_NON_NULL(stats, NULL);
     18    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
     19    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
     20    if (mask != NULL) {
     21        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
     22        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
     23    }
     24    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
     25    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
     26    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
     27    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     28    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
     29    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
     30    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
     31    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
     32    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(z, NULL);
     33    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
     34    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
     35    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(t, NULL);
     36    PS_ASSERT_VECTOR_NON_NULL(fErr, NULL);
     37    PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
     38    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
     39
     40    // clipping range defined by min and max and/or clipSigma
     41    float minClipSigma;
     42    float maxClipSigma;
     43    if (isfinite(stats->max)) {
     44        maxClipSigma = +fabs(stats->max);
     45    } else {
     46        maxClipSigma = +fabs(stats->clipSigma);
     47    }
     48    if (isfinite(stats->min)) {
     49        minClipSigma = -fabs(stats->min);
     50    } else {
     51        minClipSigma = -fabs(stats->clipSigma);
     52    }
     53    psVector *fit   = NULL;
     54    psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64);
     55
     56    // eventual expansion: user supplies one of various stats option pairs,
     57    // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to
     58    // evaluate the clipping sigma
     59    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
     60    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     61
     62    for (int N = 0; N < stats->clipIter; N++) {
     63        int Nkeep = 0;
     64
     65        poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t);
     66        fit = psPolynomial4DEvalVector (poly, x, y, z, t);
     67        resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit);
     68
     69        stats  = psVectorStats (stats, resid, NULL, mask, maskValue);
     70        psTrace (".psphot.VectorClipFit", 5, "resid stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
     71
     72        // set mask if pts are not valid
     73        // we are masking out any point which is out of range
     74        // recovery is not allowed with this scheme
     75        for (int i = 0; i < resid->n; i++) {
     76            if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
     77                continue;
     78            }
     79            float sigma = hypot (psVectorGet (fErr, i), stats->sampleStdev);
     80            if (resid->data.F64[i] - stats->sampleMedian > sigma*maxClipSigma) {
     81                if (mask != NULL) {
     82                    mask->data.U8[i] |= 0x01;
     83                }
     84                continue;
     85            }
     86            if (resid->data.F64[i] - stats->sampleMedian < sigma*minClipSigma) {
     87                if (mask != NULL) {
     88                    mask->data.U8[i] |= 0x01;
     89                }
     90                continue;
     91            }
     92            Nkeep ++;
     93        }
     94
     95        psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n",
     96                 Nkeep, x->n);
     97
     98        stats->clippedNvalues = Nkeep;
     99        psFree (fit);
     100    }
     101    // Free local temporary variables
     102    psFree (resid);
     103
     104    if (poly == NULL) {
     105        psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data.  Returning NULL.\n");
     106        return(NULL);
     107    }
     108    return(poly);
     109}
     110
     111// measure the aperture residual statistics
     112bool psphotApResid (eamReadout *imdata, psArray *sources, psMetadata *config, pmPSF *psf) {
    4113
    5114    int Npsf;
    … …  
    8117    pmSource *source;
    9118
    10     // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
    11     float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
    12 
    13119    psTimerStart ("psphot");
    14120
    15     // XXX EAM : drop this if we know the list is sorted
    16     sources = psArraySort (sources, psphotSortBySN);
    17 
     121    // measure the aperture loss as a function of radius for PSF
     122    float REF_RADIUS = psMetadataLookupF32 (&status, config, "PSF_REF_RADIUS");
     123    psf->growth = pmGrowthCurveAlloc (3.0, REF_RADIUS, 0.1);
     124    psphotGrowthCurve (imdata, psf);
     125   
    18126    psVector *mask    = psVectorAlloc (300, PS_TYPE_U8);
    19127    psVector *xPos    = psVectorAlloc (300, PS_TYPE_F64);
    20128    psVector *yPos    = psVectorAlloc (300, PS_TYPE_F64);
    21     psVector *rflux   = psVectorAlloc (300, PS_TYPE_F64);
     129    psVector *flux    = psVectorAlloc (300, PS_TYPE_F64);
     130    psVector *r2rflux = psVectorAlloc (300, PS_TYPE_F64);
    22131    psVector *apResid = psVectorAlloc (300, PS_TYPE_F64);
    23     mask->n = xPos->n = yPos->n = rflux->n = apResid->n = 0;
     132    psVector *dMag    = psVectorAlloc (300, PS_TYPE_F64);
     133    mask->n = xPos->n = yPos->n = flux->n = r2rflux->n = apResid->n = dMag->n = 0;
    24134    Npsf = 0;
    25135
    26     // select the NNN brightest, non-saturated sources, or just select PSFSTARs?
    27     for (int i = 0; (i < sources->n) && (Npsf < 300); i++) {
     136    // select all good PM_SOURCE_STAR entries
     137    for (int i = 0; i < sources->n; i++) {
    28138        source = sources->data[i];
    29139
    … …  
    34144        if (source->mode &  PM_SOURCE_POOR) continue;
    35145
    36         // get magnitudes, uncorrected for (x, y, rflux)
    37         model = pmSourceMagnitudes (source, NULL, RADIUS);
     146        // get uncorrected magnitudes in scaled apertures
     147        model = pmSourceMagnitudes (source, NULL, 0);
    38148        if (model == NULL) continue;
    39149
    … …  
    41151        xPos->data.F64[Npsf] = model->params->data.F32[2];
    42152        yPos->data.F64[Npsf] = model->params->data.F32[3];
    43         rflux->data.F64[Npsf] = pow(10.0, 0.4*source->fitMag);
    44         apResid->data.F64[Npsf] = source->apMag - source->fitMag;
    45 
    46         psVectorExtend (mask, 100, 1);
    47         psVectorExtend (xPos, 100, 1);
    48         psVectorExtend (yPos, 100, 1);
    49         psVectorExtend (rflux, 100, 1);
     153
     154        flux->data.F64[Npsf] = pow(10.0, -0.4*source->fitMag);
     155        r2rflux->data.F64[Npsf] = PS_SQR(model->radius) / flux->data.F64[Npsf];
     156       
     157        apResid->data.F64[Npsf] = source->apMag + pmGrowthCurveCorrect (psf->growth, model->radius) - source->fitMag ;
     158
     159        // XXX sanity clip?
     160        // XXX need to see if all data were tossed?
     161        // XXX need to subtract median?
     162        if (fabs(apResid->data.F64[Npsf]) > 0.2) continue;
     163
     164        dMag->data.F64[Npsf] = model->dparams->data.F32[1] / model->params->data.F32[1];
     165
     166        psVectorExtend (mask,    100, 1);
     167        psVectorExtend (xPos,    100, 1);
     168        psVectorExtend (yPos,    100, 1);
     169        psVectorExtend (flux,    100, 1);
     170        psVectorExtend (r2rflux, 100, 1);
     171        psVectorExtend (dMag,    100, 1);
    50172        psVectorExtend (apResid, 100, 1);
    51173        Npsf ++;
    52174    }
    53     psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec\n", psTimerMark ("psphot"));
     175    psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec for %d objects\n", psTimerMark ("psphot"), Npsf);
     176
     177    // APTREND options : NONE SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL
     178    char *ApTrendOption = psMetadataLookupPtr (&status, config, "APTREND");
     179    if (!status) ApTrendOption = psStringCopy ("SKYBIAS");
    54180
    55181    // 3hi/1lo sigma clipping on the rflux vs metric fit
    56182    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    57     stats->min = 1.0;
     183    stats->min = 3.0;
    58184    stats->max = 3.0;
    59     stats->clipIter = 3;
    60 
    61     // first clip out objects which are too far from the median
    62     stats->clipIter = 1;
    63     maskToConstant (psf->ApTrend);
    64     psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux);
    65 
    66     // next, fit just SkyBias and clip out objects which are too far from the median
    67     stats->clipIter = 2;
    68     maskToSkyBias (psf->ApTrend);
    69     psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux);
    70 
    71     // finally, fit x, y, SkyBias and clip out objects which are too far from the median
    72     stats->clipIter = 2;
    73     maskToDefault (psf->ApTrend);
    74     psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux);
    75 
    76     // linear clipped fit of apResid to rflux, xPos, yPos
    77     # if (0)
    78     psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux);
    79     psf->skyBias  = psf->ApTrend->coeff[0][0][1] / (M_PI * PS_SQR(RADIUS));
    80     psf->ApResid  = psf->ApTrend->coeff[0][0][0];
    81     psf->dApResid = stats->sampleStdev;
    82     psf->ApTrend->coeff[0][0][1] = 0;
    83     # endif
     185
     186    // constant only
     187    if (!strcasecmp (ApTrendOption, "CONSTANT")) {
     188        stats->clipIter = 2;
     189        pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
     190        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     191
     192        stats->clipIter = 3;
     193        pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
     194        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     195    }
     196
     197    // constant and skybias only
     198    if (!strcasecmp (ApTrendOption, "SKYBIAS")) {
     199        // first clip out objects which are too far from the median
     200        stats->clipIter = 2;
     201        pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
     202        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     203
     204        // apply the fit
     205        stats->clipIter = 3;
     206        pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS);
     207        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     208    }
     209
     210    if (!strcasecmp (ApTrendOption, "SKYSAT")) {
     211        // first clip out objects which are too far from the median
     212        stats->clipIter = 2;
     213        pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
     214        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     215
     216        // apply the fit
     217        stats->clipIter = 2;
     218        pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS);
     219        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     220
     221        // apply the fit
     222        stats->clipIter = 3;
     223        pmPSF_MaskApTrend (psf, PM_PSF_SKYSAT);
     224        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     225    }
     226
     227    // constant and linear X,Y only
     228    if (!strcasecmp (ApTrendOption, "XY_LIN")) {
     229        // first clip out objects which are too far from the median
     230        stats->clipIter = 2;
     231        pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
     232        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     233
     234        // apply the fit
     235        stats->clipIter = 3;
     236        pmPSF_MaskApTrend (psf, PM_PSF_XY_LIN);
     237        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     238    }
     239
     240    // constant and quadratic X,Y only
     241    if (!strcasecmp (ApTrendOption, "XY_QUAD")) {
     242        // first clip out objects which are too far from the median
     243        stats->clipIter = 2;
     244        pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
     245        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     246
     247        // apply the fit
     248        stats->clipIter = 3;
     249        pmPSF_MaskApTrend (psf, PM_PSF_XY_QUAD);
     250        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     251    }
     252
     253     // constant and sky, linear X,Y only
     254    if (!strcasecmp (ApTrendOption, "SKY_XY_LIN")) {
     255        // first clip out objects which are too far from the median
     256        stats->clipIter = 2;
     257        pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
     258        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     259
     260        // apply the fit
     261        stats->clipIter = 3;
     262        pmPSF_MaskApTrend (psf, PM_PSF_SKY_XY_LIN);
     263        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     264    }
     265
     266     // constant and sky, linear X,Y only
     267    if (!strcasecmp (ApTrendOption, "SKYSAT_XY_LIN")) {
     268        // first clip out objects which are too far from the median
     269        stats->clipIter = 2;
     270        pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
     271        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     272
     273        // apply the fit
     274        stats->clipIter = 3;
     275        pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS);
     276        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     277
     278        // apply the fit
     279        stats->clipIter = 3;
     280        pmPSF_MaskApTrend (psf, PM_PSF_SKYSAT_XY_LIN);
     281        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     282    }
     283
     284    if (!strcasecmp (ApTrendOption, "ALL")) {
     285        // first clip out objects which are too far from the median
     286        stats->clipIter = 2;
     287        pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
     288        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     289
     290        // fit just SkyBias and clip out objects which are too far from the median
     291        stats->clipIter = 2;
     292        pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS);
     293        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     294
     295        // finally, fit x, y, SkyBias and clip out objects which are too far from the median
     296        stats->clipIter = 3;
     297        pmPSF_MaskApTrend (psf, PM_PSF_ALL);
     298        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     299    }
     300
     301# if (1)
     302    psPolynomial4D *poly = psf->ApTrend;
     303    for (int nt = 0; nt <= poly->nT; nt++) {
     304        for (int nz = 0; nz <= poly->nZ; nz++) {
     305            for (int ny = 0; ny <= poly->nY; ny++) {
     306                for (int nx = 0; nx <= poly->nX; nx++) {
     307                    if (poly->mask[nx][ny][nz][nt]) continue;
     308                    fprintf (stderr, "%d %d %d %d : %22.15g\n", nx, ny, nz, nt, poly->coeff[nx][ny][nz][nt]);
     309                }
     310            }
     311        }
     312    }
     313# endif
     314
     315    // construct the fitted values and the residuals
     316    psVector *fit   = psPolynomial4DEvalVector (psf->ApTrend, xPos, yPos, r2rflux, flux);
     317    psVector *resid = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) fit);
     318
     319# if (0)
     320    FILE *fout = fopen ("resid.dat", "w");
     321    for (int i = 0; i < resid->n; i++) {
     322        fprintf (fout, "%d %f %f %f %f  %f %f %f %f %d\n",
     323                 i,
     324                 (float) psVectorGet(xPos, i), (float) psVectorGet(yPos, i), (float) psVectorGet(r2rflux, i), (float) psVectorGet(flux, i),
     325                 (float) psVectorGet(apResid, i), (float) psVectorGet(fit, i), (float) psVectorGet(resid, i), (float) psVectorGet(dMag, i),
     326                 mask->data.U8[i]);
     327    }
     328    fclose (fout);
     329# endif
     330
     331    // measure scatter for sources with dMag < 0.01 (S/N = 100)
     332    psStats *residStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
     333    for (int i = 0; i < dMag->n; i++) {
     334        if (dMag->data.F64[i] > 0.01) {
     335            mask->data.U8[i] |= 0x02;
     336        }
     337    }
     338    residStats  = psVectorStats (residStats, resid, NULL, mask, 0x03);
     339
     340    // apply ApTrend results
     341    psf->skyBias  = psf->ApTrend->coeff[0][0][1][0];
     342    psf->skySat   = psf->ApTrend->coeff[0][0][0][1];
     343    psf->ApResid  = psf->ApTrend->coeff[0][0][0][0];
     344    psf->dApResid = residStats->sampleStdev;
     345    psf->ApTrend->coeff[0][0][1][0] = 0;
     346    psf->ApTrend->coeff[0][0][0][1] = 0;
     347    psf->nApResid = residStats->clippedNvalues;
    84348
    85349    /*
    … …  
    89353    */
    90354
    91     # if (0)
    92     psPolynomial3D *poly = psf->ApTrend;
    93     for (int nz = 0; nz <= poly->nZ; nz++) {
    94         for (int ny = 0; ny <= poly->nY; ny++) {
    95             for (int nx = 0; nx <= poly->nX; nx++) {
    96                 fprintf (stderr, "%d %d %d : %22.15g\n", nx, ny, nz, poly->coeff[nx][ny][nz]);
    97             }
    98         }
    99     }
    100     # endif
    101 
    102355    psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot"));
    103356    psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias\n", psf->ApResid, psf->dApResid, psf->skyBias);
    104 
    105     psFree (stats);
    106     psFree (mask);
    107     psFree (rflux);
    108     psFree (apResid);
     357    psLogMsg ("psphot.apresid", 4, "apresid trends: %f %f %f %f %f\n",
     358              1e3*psf->ApTrend->coeff[1][0][0][0],
     359              1e6*psf->ApTrend->coeff[2][0][0][0],
     360              1e6*psf->ApTrend->coeff[1][1][0][0],
     361              1e3*psf->ApTrend->coeff[0][1][0][0],
     362              1e6*psf->ApTrend->coeff[0][2][0][0]);
     363
     364    // psFree (stats);
     365    // psFree (mask);
     366    // psFree (rflux);
     367    // psFree (apResid);
    109368
    110369    return true;
    111370}
     371

  • trunk/psphot/src/psphotArguments.c

    r5837 r5980  
    77    int N;
    88    unsigned int Nfail;
    9     int mode = PS_DATA_STRING | PS_META_REPLACE;
     9
     10    psMetadata *options = psMetadataAlloc ();
    1011
    1112    // basic pslib options
    … …  
    1314    psArgumentVerbosity (argc, argv);
    1415
    15     // optional mask image - add to config
    16     char *mask = NULL;
     16    // command-line options -- place on options MD
     17    // mask image
    1718    if ((N = psArgumentGet (*argc, argv, "-mask"))) {
    1819        psArgumentRemove (N, argc, argv);
    19         mask = psStringCopy (argv[N]);
     20        psMetadataAddStr (options, PS_LIST_TAIL, "MASK_IMAGE", 0, "", psStringCopy(argv[N]));
    2021        psArgumentRemove (N, argc, argv);
    2122    }
    2223
    23     // optional weight image - add to config
    24     char *weight = NULL;
     24    // weight image
    2525    if ((N = psArgumentGet (*argc, argv, "-weight"))) {
    2626        psArgumentRemove (N, argc, argv);
    27         weight = psStringCopy (argv[N]);
     27        psMetadataAddStr (options, PS_LIST_TAIL, "WEIGHT_IMAGE", 0, "", psStringCopy(argv[N]));
    2828        psArgumentRemove (N, argc, argv);
    2929    }
    3030
    31     // optional output residual image - add to config
    32     char *resid = NULL;
     31    // output residual image
    3332    if ((N = psArgumentGet (*argc, argv, "-resid"))) {
    3433        psArgumentRemove (N, argc, argv);
    35         resid = psStringCopy (argv[N]);
     34        psMetadataAddStr (options, PS_LIST_TAIL, "RESID_IMAGE", 0, "", psStringCopy(argv[N]));
    3635        psArgumentRemove (N, argc, argv);
    3736    }
    3837
    39     // optional analysis region - add to config
    40     char *region = NULL;
     38    // analysis region
    4139    if ((N = psArgumentGet (*argc, argv, "-region"))) {
    4240        psArgumentRemove (N, argc, argv);
    43         region = psStringCopy (argv[N]);
     41        psMetadataAddStr (options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", psStringCopy(argv[N]));
    4442        psArgumentRemove (N, argc, argv);
    4543    }
    4644
    47     // optional output residual image - add to config
    48     char *photcode = NULL;
     45    // output residual image
     46    psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", 0, "", psStringCopy("NONE"));
    4947    if ((N = psArgumentGet (*argc, argv, "-photcode"))) {
    5048        psArgumentRemove (N, argc, argv);
    51         photcode = psStringCopy (argv[N]);
     49        psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", psStringCopy(argv[N]));
    5250        psArgumentRemove (N, argc, argv);
    5351    }
    5452
    55     char *psffile = NULL;
     53    // input PSF file
    5654    if ((N = psArgumentGet (*argc, argv, "-psf"))) {
    5755        psArgumentRemove (N, argc, argv);
    58         psffile = psStringCopy (argv[N]);
     56        psMetadataAddStr (options, PS_LIST_TAIL, "PSF_INPUT_FILE", 0, "", psStringCopy (argv[N]));
    5957        psArgumentRemove (N, argc, argv);
    6058    }
    6159
    62     bool ModelTest = false;
    63     float ModelTest_X = 0;
    64     float ModelTest_Y = 0;
    65     char *model = NULL;
    66     char *fitset = NULL;
    67     char *fitmode = NULL;
     60    // run the 0ltest model (requires X,Y coordinate)
    6861    if ((N = psArgumentGet (*argc, argv, "-modeltest"))) {
    69         ModelTest = true;
     62        psMetadataAddBool (options, PS_LIST_TAIL, "TEST_FIT",   0, "", true);
     63        psMetadataAddF32  (options, PS_LIST_TAIL, "TEST_FIT_X", 0, "", atof(argv[N+1]));
     64        psMetadataAddF32  (options, PS_LIST_TAIL, "TEST_FIT_Y", 0, "", atof(argv[N+2]));
     65
    7066        psArgumentRemove (N, argc, argv);
    71         ModelTest_X = atof (argv[N]);
    7267        psArgumentRemove (N, argc, argv);
    73         ModelTest_Y = atof (argv[N]);
    7468        psArgumentRemove (N, argc, argv);
    7569
     70        // specify the modeltest model
    7671        if ((N = psArgumentGet (*argc, argv, "-model"))) {
    7772            psArgumentRemove (N, argc, argv);
    78             model = psStringCopy (argv[N]);
     73            psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODEL", 0, "", psStringCopy (argv[N]));
    7974            psArgumentRemove (N, argc, argv);
    8075        }
    8176
     77        // specify the test fit mode
    8278        if ((N = psArgumentGet (*argc, argv, "-fitmode"))) {
    8379            psArgumentRemove (N, argc, argv);
    84             fitmode = psStringCopy (argv[N]);
     80            psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODE", 0, "", psStringCopy (argv[N]));
    8581            psArgumentRemove (N, argc, argv);
    8682        }
    8783        if ((N = psArgumentGet (*argc, argv, "-fitset"))) {
    8884            psArgumentRemove (N, argc, argv);
    89             fitset = psStringCopy (argv[N]);
     85            psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_SET", 0, "", psStringCopy (argv[N]));
    9086            psArgumentRemove (N, argc, argv);
    9187        }
     88    }
     89
     90    // other arbitrary options: -D key value (all added as string)
     91    while ((N = psArgumentGet (*argc, argv, "-D"))) {
     92        psArgumentRemove (N, argc, argv);
     93        psMetadataAddStr (options, PS_LIST_TAIL, argv[N], 0, "", argv[N+1]);
     94        psArgumentRemove (N, argc, argv);
     95        psArgumentRemove (N, argc, argv);
     96    }
     97
     98    // other arbitrary options: -Df key value (all added as float)
     99    while ((N = psArgumentGet (*argc, argv, "-Df"))) {
     100        psArgumentRemove (N, argc, argv);
     101        psMetadataAddF32 (options, PS_LIST_TAIL, argv[N], 0, "", atof(argv[N+1]));
     102        psArgumentRemove (N, argc, argv);
     103        psArgumentRemove (N, argc, argv);
     104    }
     105
     106    // other arbitrary options: -Di key value (all added as int)
     107    while ((N = psArgumentGet (*argc, argv, "-Di"))) {
     108        psArgumentRemove (N, argc, argv);
     109        psMetadataAddS32 (options, PS_LIST_TAIL, argv[N], 0, "", atoi(argv[N+1]));
     110        psArgumentRemove (N, argc, argv);
     111        psArgumentRemove (N, argc, argv);
    92112    }
    93113
    … …  
    96116    // load config information
    97117    psMetadata *config = psMetadataAlloc ();
    98     psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL);
     118
     119    psMetadataAdd (config, PS_LIST_TAIL, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL);
    99120    config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE);
     121
     122    psMetadataItem *item = NULL;
     123    psMetadataIterator *iter = psMetadataIteratorAlloc (options, PS_LIST_HEAD, NULL);
     124    while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
     125        psMetadataAddItem (config, item, PS_LIST_TAIL, PS_META_REPLACE);
     126        // psFree (item);
     127    }
     128    // psFree (iter);
    100129
    101130    // identify input image & optional weight & mask images
    102131    // command-line entries override config-file entries
    103     psMetadataAdd (config, PS_LIST_HEAD, "IMAGE",       mode, "", argv[1]);
    104     psMetadataAdd (config, PS_LIST_HEAD, "OUTPUT_FILE", mode, "", argv[2]);
     132    psMetadataAddStr (config, PS_LIST_TAIL, "IMAGE",       PS_META_REPLACE, "", argv[1]);
     133    psMetadataAddStr (config, PS_LIST_TAIL, "OUTPUT_FILE", PS_META_REPLACE, "", argv[2]);
    105134
    106     if (mask != NULL) {
    107         psMetadataAdd (config, PS_LIST_HEAD, "MASK_IMAGE", mode, "", mask);
    108     }
    109     if (weight != NULL) {
    110         psMetadataAdd (config, PS_LIST_HEAD, "WEIGHT_IMAGE", mode, "", weight);
    111     }
    112     if (resid != NULL) {
    113         psMetadataAdd (config, PS_LIST_HEAD, "RESID_IMAGE", mode, "", resid);
    114     }
    115     if (region != NULL) {
    116         psMetadataAdd (config, PS_LIST_HEAD, "ANALYSIS_REGION", mode, "", region);
    117     }
    118     if (photcode != NULL) {
    119         psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", photcode);
    120     } else {
    121         psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", "NONE");
    122     }   
    123     if (psffile != NULL) {
    124         psMetadataAdd (config, PS_LIST_HEAD, "PSF_INPUT_FILE", mode, "", psffile);
    125     }
    126 
    127     // model related options
    128     if (ModelTest) {
    129         psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT", PS_DATA_BOOL, "", true);
    130 
    131         int fmode = PS_DATA_F32 | PS_META_REPLACE;
    132         psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_X", fmode, "", ModelTest_X);
    133         psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_Y", fmode, "", ModelTest_Y);
    134 
    135         if (model != NULL) {
    136             psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_MODEL", mode, "", model);
    137         }
    138 
    139         if (fitmode != NULL) {
    140             psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_MODE", mode, "", fitmode);
    141         }
    142         if (fitset != NULL) {
    143             psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_SET", mode, "", fitset);
    144         }
    145     }
    146135    return (config);
    147136}

  • trunk/psphot/src/psphotDefinePixels.c

    r5134 r5980  
    1616    mySource->weight = psImageSubset(imdata->weight, srcRegion);
    1717    mySource->mask   = psImageSubset(imdata->mask,  srcRegion);
     18    mySource->region = srcRegion;
    1819
    19     return(mySource);
     20    return true;
    2021}
    2122
     23bool psphotRedefinePixels(pmSource *mySource,
     24                          const eamReadout *imdata,
     25                          psF32 x,
     26                          psF32 y,
     27                          psF32 Radius)
     28{
     29    bool extend;
     30    psRegion newRegion;
     31
     32    if (Radius == 0) return false;
     33
     34    // check to see if new region is completely contained within old region
     35    newRegion = psRegionForSquare (x, y, Radius);
     36    newRegion = psRegionForImage (imdata->image, newRegion);
     37
     38    extend = false;
     39    extend |= (int)(newRegion.x0) < (int)(mySource->region.x0);
     40    extend |= (int)(newRegion.x1) > (int)(mySource->region.x1);
     41    extend |= (int)(newRegion.y0) < (int)(mySource->region.y0);
     42    extend |= (int)(newRegion.y1) > (int)(mySource->region.y1);
     43
     44    extend |= (mySource->pixels == NULL);
     45    extend |= (mySource->weight == NULL);
     46    extend |= (mySource->mask == NULL);
     47
     48    // extend = true;
     49    if (extend) {
     50        // Grab a new subimage
     51        // psFree (mySource->pixels);
     52        // psFree (mySource->weight);
     53        // psFree (mySource->mask);
     54
     55        // fprintf (stderr, "extend %f,%f\n", x, y);
     56        mySource->pixels = psImageSubset(imdata->image,  newRegion);
     57        mySource->weight = psImageSubset(imdata->weight, newRegion);
     58        mySource->mask   = psImageSubset(imdata->mask,   newRegion);
     59        mySource->region = newRegion;
     60    }
     61
     62    return extend;
     63}
     64
     65
     66//**** this function was ALWAYS restricting to model->radius
     67//     should have left it larger (PSF_FIT_RAD) for the faint sources

  • trunk/psphot/src/psphotEnsemblePSF.c

    r5837 r5980  
    5858    sources = psArraySort (sources, psphotSortByY);
    5959
    60     float OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
    6160    float INNER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS");
    6261    float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
    … …  
    104103        // really saturated stars should be re-measured for a better centroid
    105104        // XXX EAM : move this to a 'clear satstar function'
    106         // XXX EAM : extend size of fit box around SATSTAR
    107105        if (inSource->mode &  PM_SOURCE_SATSTAR) {
    108106            status = pmSourceMoments (inSource, INNER_RADIUS);
    … …  
    143141        y = model->params->data.F32[3];
    144142
    145         // get function which specifies the radius at this the model hits the given flux
     143        // get function which specifies the radius at which the model hits the given flux
    146144        pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type);
    147145
    … …  
    151149
    152150        // if needed, ask for more object pixels
    153         if (model->radius > OUTER_RADIUS) {
    154             psphotDefinePixels (inSource, imdata, x, y, model->radius);
    155         }
     151        psphotRedefinePixels (inSource, imdata, x, y, model->radius);
    156152
    157153        // make temporary copies of the image pixels and mask

  • trunk/psphot/src/psphotFixedPSF.c

    r5828 r5980  
    1313
    1414    // we may set this differently here from the value used to mark likely saturated stars
    15     float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
    16 
    1715    float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
    1816    float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
    … …  
    5250        }
    5351       
    54         // check if we need to redefine the pixels
    55         // XXX EAM : a better test would examine the source pixels
    56         if (model->radius > OUTER_RADIUS) {
    57           // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER)
    58             psphotDefinePixels (source, imdata, x, y, model->radius);
    59         }
     52        psphotRedefinePixels (source, imdata, x, y, model->radius);
    6053
    6154        // fit PSF model (set/unset the pixel mask)

  • trunk/psphot/src/psphotMagnitudes.c

    r5837 r5980  
    1717        model = source->modelPSF;
    1818        if (model == NULL) return NULL;
    19         radius = apRadius;
     19        radius = (apRadius > 0) ? apRadius : model->radius;
    2020        isPSF = true;
    2121        break;
    … …  
    4646    // for PSFs, correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS
    4747    if (isPSF && (psf != NULL)) {
     48      if (psf->growth != NULL) {
     49          source->apMag += pmGrowthCurveCorrect (psf->growth, model->radius);
     50      }
     51
    4852      rflux   = pow (10.0, 0.4*source->fitMag);
    49       source->apMag  -= rflux * psf->skyBias * (M_PI * PS_SQR(apRadius));
    50       source->fitMag += psPolynomial3DEval (psf->ApTrend, x, y, 0.0);
     53      source->apMag  -= PS_SQR(model->radius)*rflux * psf->skyBias + psf->skySat / rflux;
     54      source->fitMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0);
    5155    }
     56
     57    // unmask aperture
     58    psImageKeepCircle (source->mask, x, y, radius, "AND", ~PSPHOT_MASK_MARKED);
    5259
    5360    // subtract object, leave local sky
    5461    pmSourceSubModel (source->pixels, source->mask, model, false, false);
    55 
    56     // unmask aperture
    57     psImageKeepCircle (source->mask, x, y, radius, "AND", ~PSPHOT_MASK_MARKED);
    5862
    5963    if (!status) return NULL;
    … …  
    6266
    6367/*
    64 aprMag' - fitMag = rflux*skyBias + ApTrend(x,y)
    65 (aprMag - rflux*skyBias) - fitMAg = ApTrend(x,y)
    66 (aprMag - rflux*skyBias) = fitMAg + ApTrend(x,y)
     68aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y)
     69(aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y)
     70(aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y)
     71
    6772*/

  • trunk/psphot/src/psphotOutput.c

    r5837 r5980  
    6565    }
    6666
    67     float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
     67    float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS");
    6868
    6969    // write sources with models
    … …  
    116116    }
    117117
    118     float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
     118    float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS");
    119119
    120120    // write sources with models
    … …  
    162162
    163163    // find config information for output header
    164     float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
     164    float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS");
    165165    float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT");
    166166    if (!status) ZERO_POINT = 25.0;
    … …  
    171171    psMetadataAdd (imdata->header, PS_LIST_TAIL, "ZERO_PT",  PS_DATA_F32 | PS_META_REPLACE, "zero point",          ZERO_POINT);
    172172    psMetadataAdd (imdata->header, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
    173     psMetadataAdd (imdata->header, PS_LIST_TAIL, "dAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
     173    psMetadataAdd (imdata->header, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
     174    psMetadataAdd (imdata->header, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid);
     175    psMetadataAdd (imdata->header, PS_LIST_TAIL, "NPSFSTAR", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nPSFstars);
    174176    psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skyBias);
     177    psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skySat);
    175178    psMetadataAdd (imdata->header, PS_LIST_TAIL, "PHOTCODE", PS_DATA_STRING | PS_META_REPLACE, "photometry code",     PHOTCODE);
    176179    psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_X",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM X",          0.0);
    … …  
    250253
    251254    // find config information for output header
    252     float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
     255    float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS");
    253256    float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT");
    254257    char *PHOTCODE = psMetadataLookupPtr (&status, config, "PHOTCODE");
    … …  
    367370    bool status;
    368371
    369     float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
     372    float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS");
    370373
    371374    f = fopen (filename, "w");
    … …  
    400403            fprintf (f, "%9.6f ", dPAR[j]);
    401404        }
    402         fprintf (f, ": %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",
    403                  source[0].type, source[0].mode,
     405        fprintf (f, ": %7.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",
     406                 source[0].apMag, source[0].type, source[0].mode,
    404407                 log10(model[0].chisq/model[0].nDOF),
    405408                 source[0].moments->SN,
    … …  
    453456            fprintf (f, "%9.6f ", dPAR[j]);
    454457        }
    455         fprintf (f, ": %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",
     458        fprintf (f, ": %7.4f  %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",
     459                 source->apMag,
    456460                 source[0].type, source[0].mode,
    457461                 log10(model[0].chisq/model[0].nDOF),
    … …  
    694698}
    695699
     700psPolynomial4D *psPolynomial4DfromMD (psMetadata *folder) {
     701
     702    bool status;
     703    char keyword[80];
     704
     705    // get polynomial orders
     706    // XXX add status failures tests
     707    int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
     708    int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
     709    int nZorder = psMetadataLookupS32 (&status, folder, "NORDER_Z");
     710    int nTorder = psMetadataLookupS32 (&status, folder, "NORDER_T");
     711
     712    psPolynomial4D *poly = psPolynomial4DAlloc (nXorder, nYorder, nZorder, nTorder, PS_POLYNOMIAL_ORD);
     713
     714    for (int nx = 0; nx < poly->nX + 1; nx++) {
     715        for (int ny = 0; ny < poly->nY + 1; ny++) {
     716            for (int nz = 0; nz < poly->nZ + 1; nz++) {
     717                for (int nt = 0; nt < poly->nT + 1; nt++) {
     718                    sprintf (keyword, "VAL_X%02d_Y%02d_Z%02d_T%02d", nx, ny, nz, nt);
     719                    poly->coeff[nx][ny][nz][nt] = psMetadataLookupF32 (&status, folder, keyword);
     720                    if (!status) poly->mask[nx][ny][nz][nt] = 1;
     721                   
     722                    sprintf (keyword, "ERR_X%02d_Y%02d_Z%02d_T%02d", nx, ny, nz, nt);
     723                    poly->coeffErr[nx][ny][nz][nt] = psMetadataLookupF32 (&status, folder, keyword);
     724                }
     725            }
     726        }
     727    }
     728    return (poly);
     729}
     730
    696731// XXX : these may need F64, or %g format for output
    697732bool psPolynomial2DtoMD (psMetadata *md, psPolynomial2D *poly, char *format, ...) {
    … …  
    774809}
    775810
     811bool psPolynomial4DtoMD (psMetadata *md, psPolynomial4D *poly, char *format, ...) {
     812
     813    int Nbyte;
     814    char tmp;
     815    char *root;
     816    va_list argp; 
     817
     818    va_start (argp, format);
     819    Nbyte = vsnprintf (&tmp, 0, format, argp);
     820    va_end (argp);
     821
     822    if (!Nbyte) return false;
     823
     824    va_start (argp, format);
     825    root = (char *) psAlloc (Nbyte + 1);
     826    memset (root, 0, Nbyte + 1);
     827    vsnprintf (root, Nbyte + 1, format, argp);
     828    va_end (argp);
     829
     830    psMetadata *folder = psMetadataAlloc ();
     831    psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 4D polynomial", folder);
     832    psFree (root);
     833
     834    // specify the polynomial orders
     835    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
     836    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
     837    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Z", PS_DATA_S32, "number of z orders", poly->nZ);
     838    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_T", PS_DATA_S32, "number of z orders", poly->nT);
     839
     840    // place polynomial entries on folder
     841    for (int nx = 0; nx < poly->nX + 1; nx++) {
     842        for (int ny = 0; ny < poly->nY + 1; ny++) {
     843            for (int nz = 0; nz < poly->nZ + 1; nz++) {
     844                for (int nt = 0; nt < poly->nT + 1; nt++) {
     845                    if (poly->mask[nx][ny][nz][nt]) continue;
     846                    psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d_Z%02d_T%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny][nz][nt], nx, ny, nz, nt);
     847                    psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d_Z%02d_T%02d", PS_DATA_F32, "polynomial coeffficient error", poly->coeffErr[nx][ny][nz][nt], nx, ny, nz, nt);
     848                }
     849            }
     850        }
     851    }
     852    return true;
     853}
     854
    776855bool psphotWritePSF (pmPSF *psf, char *filename) {
    777856
    … …  
    788867        psPolynomial2DtoMD (psfdata, poly, "PSF_PAR%02d", i);
    789868    }
    790     psPolynomial3DtoMD (psfdata, psf->ApTrend, "APTREND");
     869    psPolynomial4DtoMD (psfdata, psf->ApTrend, "APTREND");
    791870
    792871    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid);
    … …  
    828907    sprintf (keyword, "APTREND");
    829908    psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword);
    830     psPolynomial3D *poly = psPolynomial3DfromMD (folder);
     909    psPolynomial4D *poly = psPolynomial4DfromMD (folder);
    831910    psFree (psf->ApTrend);
    832911    psf->ApTrend = poly;

  • trunk/psphot/src/psphotRadiusChecks.c

    r5828 r5980  
    3131    }
    3232
    33     // check if we need to redefine the pixels
    34     // XXX EAM : a better test would examine the source pixels
     33# if (1)
     34    bool status = psphotRedefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius);
     35    return status;
     36
     37# else
     38
    3539    if (model->radius > PSF_OUTER_RADIUS) {
    3640        // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER)
    3741        psphotDefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius);
    3842        return true;
     43    } else {
     44        fprintf (stderr, "skipping %f,%f\n", model->params->data.F32[2], model->params->data.F32[3]);
    3945    }
    40 
    4146    return false;
     47# endif
    4248}
    4349
    … …  
    6874    if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius");
    6975
     76# if (1)
     77    // redefine the pixels if needed
     78    bool status = psphotRedefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius);
     79    return status;
     80
     81# else
     82
    7083    // check if we need to redefine the pixels
    7184    if (model->radius > GAL_OUTER_RADIUS) {
    … …  
    7689    }
    7790    return false;
     91# endif
    7892}

  • trunk/psphot/src/psphotTest.c

    r5837 r5980  
    11# include "psphot.h"
     2bool apResidSetPower (float value);
     3
     4void psExit (int status, char *process, char *format, ...) {
     5
     6    va_list ap;
     7
     8    va_start (ap, format);
     9    fprintf (stderr, "exiting %s\n", process);
     10    vfprintf (stderr, format, ap);
     11    va_end (ap);
     12
     13    exit (status);
     14}
    215
    316int main (int argc, char **argv) {
    417
    5     char line[80];
    6     psImage *image;
     18    // read in a file consisting of:
     19    // x y radius fitmag apmag dmag
    720
    8     image = psImageAlloc (512, 512, PS_DATA_F32);
    9     for (int j = 0; j < 512; j++) {
    10         for (int i = 0; i < 512; i++) {
    11             image->data.F32[j][i] = i;
    12         }
     21    double fit;
     22    double ap;
     23
     24    psphotTestArguments (&argc, argv);
     25
     26    FILE *f = fopen (argv[1], "r");
     27    if (f == NULL) psExit (2, "test", "no file\n");
     28
     29    float power = atof (argv[2]);
     30    apResidSetPower (power);
     31
     32    psVector *x     = psVectorAlloc (100, PS_TYPE_F64);
     33    psVector *y     = psVectorAlloc (100, PS_TYPE_F64);
     34    psVector *rad   = psVectorAlloc (100, PS_TYPE_F64);
     35    psVector *rflux = psVectorAlloc (100, PS_TYPE_F64);
     36    psVector *apres = psVectorAlloc (100, PS_TYPE_F64);
     37    psVector *dmag  = psVectorAlloc (100, PS_TYPE_F64);
     38
     39    int Npt = 0;
     40    while (fscanf (f, "%lf %lf %lf %lf %lf %lf",
     41                   &x->data.F64[Npt],
     42                   &y->data.F64[Npt],
     43                   &rad->data.F64[Npt],
     44                   &fit, &ap,
     45                   &dmag->data.F64[Npt]) != EOF) {
     46
     47        if (fabs(ap-fit) > 0.1) continue;
     48        if (fit > -10) continue;
     49
     50        apres->data.F64[Npt] = ap-fit;
     51        rflux->data.F64[Npt] = pow(10.0, 0.4*fit);
     52
     53        psVectorExtend (x,     100, 1);
     54        psVectorExtend (y,     100, 1);
     55        psVectorExtend (rad,   100, 1);
     56        psVectorExtend (rflux, 100, 1);
     57        psVectorExtend (apres, 100, 1);
     58        psVectorExtend (dmag,  100, 1);
     59        Npt ++;
    1360    }
     61    x->n = y->n = rad->n = rflux->n = apres->n = dmag->n = Npt;
    1462
    15     psMetadata *row = NULL;
    16 
    17     psArray *table = psArrayAlloc (100);
    18     table->n = 0;
    19    
    20     for (int i = 0; i < 100; i++) {
    21         row = psMetadataAlloc ();
    22         sprintf (line, "line %03d", i);
    23         psMetadataAdd (row, PS_LIST_HEAD, "FLOAT_DATA", PS_DATA_F32,    "", 0.2*i);
    24         psMetadataAdd (row, PS_LIST_HEAD, "INT_DATA",   PS_DATA_S32,    "", i);
    25         psMetadataAdd (row, PS_LIST_HEAD, "CHAR_DATA",  PS_DATA_STRING, "", line);
    26         psArrayAdd (table, 100, row);
    27         // psFree (row);
    28     }
    29 
    30     psFits *fits = psFitsOpen ("test.fits", "w");
    31     psFitsWriteImage (fits, NULL, image, 0);
    32     psFitsWriteTable (fits, NULL, table);
    33     psFitsClose (fits);
     63    psphotApResidPolyFit (x, y, rad, rflux, apres, dmag);
    3464
    3565    exit (0);

  • trunk/psphot/src/psphotTestArguments.c

    r5672 r5980  
    22static int usage ();
    33
    4 psMetadata *psphotTestArguments (int *argc, char **argv) {
    5 
    6   int N;
    7   unsigned int Nfail;
    8   int mode = PS_DATA_STRING | PS_META_REPLACE;
    9   psF32 UseCoords_X, UseCoords_Y;
     4void psphotTestArguments (int *argc, char **argv) {
    105
    116  // basic pslib options
    … …  
    138  psArgumentVerbosity (argc, argv);
    149
    15   // identify options in args list - these go on config
    16   bool UseCoords = false;
    17   if ((N = psArgumentGet (*argc, argv, "-coords"))) {
    18     UseCoords = true;
    19     psArgumentRemove (N, argc, argv);
    20     UseCoords_X = atof (argv[N]);
    21     psArgumentRemove (N, argc, argv);
    22     UseCoords_Y = atof (argv[N]);
    23     psArgumentRemove (N, argc, argv);
    24   }
    25 
    26   float threshold = -1;
    27   if ((N = psArgumentGet (*argc, argv, "-threshold"))) {
    28     psArgumentRemove (N, argc, argv);
    29     threshold = atof (argv[N]);
    30     psArgumentRemove (N, argc, argv);
    31   }
    32 
    33   float fraction = 0.5;
    34   if ((N = psArgumentGet (*argc, argv, "-frac"))) {
    35     psArgumentRemove (N, argc, argv);
    36     fraction = atof (argv[N]);
    37     psArgumentRemove (N, argc, argv);
    38   }
    39 
    4010  if (*argc != 3) usage ();
    4111
    42   // load config information
    43   psMetadata *config = psMetadataAlloc ();
    44   psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL);
    45   config = psMetadataConfigParse (config, &Nfail, argv[2], FALSE);
    46 
    47   // identify input image
    48   psMetadataAdd (config, PS_LIST_HEAD, "IMAGE", mode, "", argv[1]);
    49 
    50   mode = PS_DATA_F32 | PS_META_REPLACE;
    51   psMetadataAdd (config, PS_LIST_HEAD, "TEST_THRESHOLD", mode, "", threshold);
    52   psMetadataAdd (config, PS_LIST_HEAD, "TEST_FRACTION",  mode, "", fraction);
    53 
    54   if (UseCoords) {
    55     psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_X", mode, "", UseCoords_X);
    56     psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_Y", mode, "", UseCoords_Y);
    57   }
    58 
    59   return (config);
     12  return;
    6013}
    6114
    6215static int usage () {
    6316
    64     fprintf (stderr, "USAGE: psphotModelTest (image.fits) (config)\n");
    65     fprintf (stderr, "options: \n");
    66     fprintf (stderr, "  -coords x y   : specify object center\n");
    67     fprintf (stderr, "  -threshold (threshold) : (defaults to FRAC of peak\n");
    68     fprintf (stderr, "  -frac (FRAC) : specify how to determine threshold\n");
     17    fprintf (stderr, "USAGE: psphotModelTest (apresid.dat) (power)\n");
    6918    exit (2);
    7019}
Note: See TracChangeset for help on using the changeset viewer.