IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5058


Ignore:
Timestamp:
Sep 15, 2005, 12:12:26 PM (21 years ago)
Author:
eugene
Message:

moved ApResid, etc into pmPSF from pmPSFtry

Location:
trunk/psphot
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/Makefile

    r4977 r5058  
    33        @echo "USAGE: make psphot"
    44
    5 # CC      =       gcc -g -Wall -O2
    6 CC      =       gcc -g -Wall
     5CC      =       gcc
    76SRC     =       src
    87BIN     =       bin
    … …  
    1615# LIBS  =       -lpsmodule $(LPSLIB)
    1716LIBS    =       $(LPSLIB)
    18 CFLAGS  =       $(INCS) -std=c99 -Wall -Werror
     17CFLAGS  =       $(INCS) -std=c99 -Wall -Werror -g
    1918LFLAGS  =       $(LIBS)
    2019

  • trunk/psphot/src/pmPSF.c

    r4977 r5058  
    2424    pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF));
    2525
    26     psf->type = type;
    27     psf->chisq = 0.0;
     26    psf->type     = type;
     27    psf->chisq    = 0.0;
     28    psf->ApResid  = 0.0;
     29    psf->dApResid = 0.0;
     30    psf->skyBias  = 0.0;
    2831
    2932    Nparams = pmModelParameterCount (type);

  • trunk/psphot/src/pmPSF.h

    r4977 r5058  
    1313    psArray *params;                    ///< Model parameters (psPolynomial2D)
    1414    float chisq;                        ///< PSF goodness statistic
     15    float ApResid;
     16    float dApResid;
     17    float skyBias;
    1518    int nPSFstars;                      ///< number of stars used to measure PSF
    1619}

  • trunk/psphot/src/pmPSFtry.c

    r4977 r5058  
    3131pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName) {
    3232
     33    pmModelType type;
     34
    3335    pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry));
    3436
    3537    // XXX probably need to increment ref counter
    36     test->modelType   = pmModelSetType (modelName);
    37     test->psf         = pmPSFAlloc (test->modelType);
    38     test->sources     = psMemCopy(sources);
    39     test->modelFLT    = psArrayAlloc (sources->n);
    40     test->modelPSF    = psArrayAlloc (sources->n);
    41     test->metric      = psVectorAlloc (sources->n, PS_TYPE_F64);
    42     test->fitMag      = psVectorAlloc (sources->n, PS_TYPE_F64);
    43     test->mask        = psVectorAlloc (sources->n, PS_TYPE_U8);
    44     test->ApResid     = 0;
    45     test->dApResid    = 0;
    46     test->skyBias     = 0;
     38    type           = pmModelSetType (modelName);
     39    test->psf      = pmPSFAlloc (type);
     40    test->sources  = psMemCopy(sources);
     41    test->modelFLT = psArrayAlloc (sources->n);
     42    test->modelPSF = psArrayAlloc (sources->n);
     43    test->metric   = psVectorAlloc (sources->n, PS_TYPE_F64);
     44    test->fitMag   = psVectorAlloc (sources->n, PS_TYPE_F64);
     45    test->mask     = psVectorAlloc (sources->n, PS_TYPE_U8);
    4746
    4847    for (int i = 0; i < test->modelFLT->n; i++) {
    … …  
    8483
    8584        pmSource *source = try->sources->data[i];
    86         pmModel  *model  = pmSourceModelGuess (source, try->modelType);
     85        pmModel  *model  = pmSourceModelGuess (source, try->psf->type);
    8786        x = source->peak->x;
    8887        y = source->peak->y;
    … …  
    164163    pmPSFtryMetric (try, RADIUS);
    165164    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
    166               modelName, try->ApResid, try->dApResid, try->skyBias);
     165              modelName, try->psf->ApResid, try->psf->dApResid, try->psf->skyBias);
    167166
    168167    return (try);
    … …  
    279278  stats = psVectorStats (stats, daResid, NULL, maskB, 1);
    280279
    281   try->ApResid = poly->coeff[0];
    282   try->dApResid = stats->clippedStdev;
    283   try->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
     280  try->psf->ApResid = poly->coeff[0];
     281  try->psf->dApResid = stats->clippedStdev;
     282  try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    284283
    285284  psFree (rflux);

  • trunk/psphot/src/pmPSFtry.h

    r4977 r5058  
    55// data to try a given PSF model type
    66typedef struct {
    7     pmModelType modelType;
    87    pmPSF      *psf;
    98    psArray    *sources;      // pointers to the original sources
    … …  
    1312    psVector   *metric;
    1413    psVector   *fitMag;
    15     float       ApResid;
    16     float       dApResid;
    17     float       skyBias;
    1814} pmPSFtry;
    1915

  • trunk/psphot/src/psLibUtils.h

    r5049 r5058  
    1515# define psMemCopy(A)(psMemIncrRefCounter((A)))
    1616
    17 // XXX EAM : bug : missing from psLib/*.h
    18 bool psMinimizeGaussNewtonDelta (psVector *delta,
    19                                  const psVector *params,
    20                                  const psVector *paramMask,
    21                                  const psArray  *x,
    22                                  const psVector *y,
    23                                  const psVector *yErr,
    24                                  psMinimizeLMChi2Func func);
    25 
     17// XXX EAM : my version using varience instead of stdev
    2618bool psMinimizeGaussNewtonDelta_EAM (psVector *delta,
    2719                                     const psVector *params,
    … …  
    3224                                     psMinimizeLMChi2Func func);
    3325
    34 // minimize
     26// minimize : using varience vs sigma and parameter limits
    3527psBool       p_psMinLM_GuessABP_EAM (psImage  *Alpha, psVector *Beta, psVector *Params, const psImage  *alpha, const psVector *beta, const psVector *params, const psVector *paramMask, const psVector *beta_lim, const psVector *params_min, const psVector *params_max, psF64 lambda);
    3628psBool       psMinimizeLMChi2_EAM(psMinimization *min, psImage *covar, psVector *params, const psVector *paramMask, const psArray *x, const psVector *y, const psVector *yErr, psMinimizeLMChi2Func func);

  • trunk/psphot/src/psphot.c

    r5049 r5058  
    4747
    4848    // write out data in appropriate format
    49     psphotOutput (imdata, config, sources);
     49    psphotOutput (imdata, config, sources, psf, sky);
    5050    exit (0);
    5151}

  • trunk/psphot/src/psphot.h

    r5049 r5058  
    2626bool         psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
    2727bool         psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats);
    28 void         psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources);
     28void         psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
    2929
    3030bool         psphotMarkPSF (pmSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE);
    … …  
    3838
    3939// output functions
    40 bool         pmSourcesWriteText (eamReadout *imdata, char *filename, psArray *sources);
    41 bool         pmSourcesWriteOBJ (eamReadout *imdata, char *filename, psArray *sources);
    42 bool         pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources);
    43 bool         pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources);
    44 bool         pmSourcesWriteSX (eamReadout *imdata, char *filename, psArray *sources);
     40bool         pmSourcesWriteText (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
     41bool         pmSourcesWriteOBJ (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
     42bool         pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
     43bool         pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
     44bool         pmSourcesWriteSX (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
    4545int          pmSourcesDophotType (pmSource *source);
    4646bool         pmPeaksWriteText (psArray *sources, char *filename);

  • trunk/psphot/src/psphotApplyPSF.c

    r5049 r5058  
    2828    // set the object surface-brightness limit for fitted pixels
    2929    float FLUX_LIMIT  = PSF_FIT_NSIGMA * sky->sampleStdev;
    30     psLogMsg ("psphot.apply_psf_model", 3, "fitting pixels with at least %f object counts\n", FLUX_LIMIT);
     30    psLogMsg ("psphot.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT);
    3131
    3232    // this function specifies the radius at this the model hits the given flux
    … …  
    7272        psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
    7373        if (!status || (model->params->data.F32[1] < 0)) {
    74           psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter);
     74          psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter);
    7575          source->type = PM_SOURCE_FAIL_FIT_PSF;  // better choice?
    7676          psFree (model);
    … …  
    9191
    9292    psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
    93     psLogMsg ("psphot", 3, "subtracted %d PSF objects\n", Nsub);
     93    psLogMsg ("psphot", 4, "subtracted %d PSF objects\n", Nsub);
    9494    return (true);
    9595}

  • trunk/psphot/src/psphotChoosePSF.c

    r5048 r5058  
    1212
    1313    // array to store candidate PSF stars
     14    int NSTARS = psMetadataLookupS32 (&status, config, "PSF_MAX_NSTARS");
     15    if (!status) NSTARS = sources->n;
    1416    stars = psArrayAlloc (sources->n);
    1517    stars->n = 0;
    1618
    1719    // select the candidate PSF stars (pointers to original sources)
    18     for (int i = 0; i < sources->n; i++) {
     20    for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) {
    1921        pmSource *source = sources->data[i];
    2022        if (source->type != PM_SOURCE_PSFSTAR) continue;
    … …  
    5557    try = models->data[0];
    5658    int bestN = 0;
    57     float bestM = try->dApResid;
     59    float bestM = try->psf->dApResid;
    5860    for (int i = 1; i < models->n; i++) {
    5961        try = models->data[i];
    60         float M = try->dApResid;
     62        float M = try->psf->dApResid;
    6163        if (M < bestM) {
    6264            bestM = M;
    … …  
    6769    // keep only the selected model:
    6870    try = models->data[bestN];
    69     modelName = pmModelGetType (try->modelType);
    70     psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, try->ApResid, try->dApResid);
     71    pmPSF *psf = psMemCopy(try->psf);
     72    psFree (models);                             // keep only the pmPSF resulting from this analysis
    7173
    72     // keep only the pmPSF resulting from this analysis
    73     pmPSF *psf = psMemCopy(try->psf);
    74     psFree (models);
     74    modelName = pmModelGetType (psf->type);
     75    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
     76
    7577    return (psf);
    7678}

  • trunk/psphot/src/psphotFitGalaxies.c

    r5049 r5058  
    8181        if (!status) {
    8282          // if the fit fails, we need to change the classification
    83           psLogMsg ("psphot", 3, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
     83          psLogMsg ("psphot", 5, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
    8484          source->type = PM_SOURCE_FAIL_FIT_GAL;  // better choice?
    8585          source->modelFLT = model;
    … …  
    9292        if (!goodfit) {
    9393          // if the fit fails, we need to change the classification
    94           psLogMsg ("psphot", 3, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
     94          psLogMsg ("psphot", 5, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
    9595          source->type = PM_SOURCE_POOR_FIT_GAL;  // better choice?
    9696          source->modelFLT = model;

  • trunk/psphot/src/psphotImageStats.c

    r5049 r5058  
    4949    sky->sampleStdev  = sqrt(sky->sampleMean/GAIN + PS_SQR(NOISE));
    5050
    51     psLogMsg ("psphot", 3, "stats: %f sec\n", psTimerMark ("psphot"));
     51    psLogMsg ("psphot", 4, "stats: %f sec\n", psTimerMark ("psphot"));
    5252    psLogMsg ("psphot", 3, "background: %f +/- %f\n", sky->sampleMean, sky->sampleStdev);
    5353

  • trunk/psphot/src/psphotMarkPSF.c

    r5048 r5058  
    3636    if (source->modelPSF->params->data.F32[1] >= SATURATE) {
    3737        if (source->type == PM_SOURCE_PSFSTAR) {
    38             psLogMsg ("psphot", 3, "PSFSTAR marked SATSTAR\n");
     38            psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n");
    3939        }
    4040        source->type = PM_SOURCE_SATSTAR;
    … …  
    4242    }
    4343    if (source->type == PM_SOURCE_SATSTAR) {
    44         psLogMsg ("psphot", 4, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n");
     44        psLogMsg ("psphot", 5, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n");
    4545        source->type = PM_SOURCE_GOODSTAR;
    4646    }
    … …  
    7070   
    7171    if (source->type == PM_SOURCE_PSFSTAR) {
    72         psLogMsg ("psphot", 3, "PSFSTAR demoted based on fit quality\n");
     72        psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality\n");
    7373    }
    7474

  • trunk/psphot/src/psphotOutput.c

    r5049 r5058  
    22
    33// output functions: we have several fixed modes (sx, obj, cmp)
    4 void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources) {
     4void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) {
    55
    66    bool status;
    … …  
    2121    }
    2222    if (!strcasecmp (outputMode, "TEXT")) {
    23         pmSourcesWriteText (imdata, outputFile, sources);
     23        pmSourcesWriteText (imdata, outputFile, sources, psf, sky);
    2424        return;
    2525    }
    2626
    2727    if (!strcasecmp (outputMode, "OBJ")) {
    28         pmSourcesWriteOBJ (imdata, outputFile, sources);
     28        pmSourcesWriteOBJ (imdata, outputFile, sources, psf, sky);
    2929        return;
    3030    }
    3131 
    3232    if (!strcasecmp (outputMode, "SX")) {
    33         pmSourcesWriteSX (imdata, outputFile, sources);
     33        pmSourcesWriteSX (imdata, outputFile, sources, psf, sky);
    3434        return;
    3535    }
    3636 
    3737    if (!strcasecmp (outputMode, "CMP")) {
    38         pmSourcesWriteCMP (imdata, outputFile, sources);
     38        pmSourcesWriteCMP (imdata, outputFile, sources, psf, sky);
    3939        return;
    4040    }
    4141 
    4242    if (!strcasecmp (outputMode, "CMF")) {
    43         pmSourcesWriteCMF (imdata, outputFile, sources);
     43        pmSourcesWriteCMF (imdata, outputFile, sources, psf, sky);
    4444        return;
    4545    }
    … …  
    4747}
    4848
    49 bool pmSourcesWriteText (eamReadout *imdata, char *filename, psArray *sources) {
     49bool pmSourcesWriteText (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky) {
    5050
    5151    char *name = (char *) psAlloc (strlen(filename) + 10);
    … …  
    6868
    6969// dophot-style output list with fixed line width
    70 bool pmSourcesWriteOBJ (eamReadout *imdata, char *filename, psArray *sources) {
     70bool pmSourcesWriteOBJ (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
    7171
    7272    int i, type, status;
    7373    pmModel *model;
    7474    psF32 *PAR, *dPAR;
    75     float sky, dmag, apMag, fitMag;
     75    float skyFlux, dmag, apMag, fitMag;
    7676
    7777    psLine *line = psLineAlloc (104);  // 104 is dophot-defined line length
    … …  
    118118
    119119        // save local sky for later
    120         sky = model->params->data.F32[0];
     120        skyFlux = model->params->data.F32[0];
    121121
    122122        // replace source flux
    … …  
    124124        model->params->data.F32[0] = 0;
    125125        pmSourceAddModel (source->pixels, source->mask, model, false);
    126         model->params->data.F32[0] = sky;
     126        model->params->data.F32[0] = skyFlux;
    127127
    128128        // measure object photometry
    129129        status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask);
     130        fitMag += psf->ApResid;
    130131
    131132        // subtract object, leave local sky
    … …  
    133134        model->params->data.F32[0] = 0;
    134135        pmSourceSubModel (source->pixels, source->mask, model, false);
    135         model->params->data.F32[0] = sky;
     136        model->params->data.F32[0] = skyFlux;
    136137
    137138        if (status == FALSE) continue;
    … …  
    159160
    160161// elixir/sextractor-style output list with fixed line width
    161 bool pmSourcesWriteSX (eamReadout *imdata, char *filename, psArray *sources) {
     162bool pmSourcesWriteSX (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
    162163
    163164    int i, status, flags;
    … …  
    224225        // measure object photometry
    225226        status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask);
     227        fitMag += psf->ApResid;
    226228
    227229        // subtract object, leave local sky
    … …  
    257259
    258260// elixir-style pseudo FITS table (header + ascii list)
    259 bool pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources) {
     261bool pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
    260262
    261263    int i, type, status;
    … …  
    328330        // measure object photometry
    329331        status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask);
     332        fitMag += psf->ApResid;
    330333
    331334        // subtract object, leave local sky
    … …  
    360363
    361364// elixir-style FITS table output (header + table in 1st extension)
    362 bool pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources) {
     365bool pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
    363366
    364367    // write imdata->header to file

  • trunk/psphot/src/psphotSetup.c

    r5049 r5058  
    2828    psFits *file = psFitsAlloc (input);
    2929    header = psFitsReadHeader (header, file);
    30     image = psFitsReadImage (NULL, file, region, 0);
     30    psImage *tmpimage = psFitsReadImage (NULL, file, region, 0);
     31    image = psImageCopy (NULL, tmpimage, PS_TYPE_F32);
    3132    psFree (file);
     33   
    3234    // psFree (input);
    3335    // XXX EAM - does the current API expect to get a copy or a view?
    … …  
    9092    }
    9193
    92     psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot"));
     94    psLogMsg ("psphot", 4, "load data: %f sec\n", psTimerMark ("psphot"));
    9395
    9496    // save the image data & return it

  • trunk/psphot/src/psphotSourceStats.c

    r5049 r5058  
    4747    }
    4848
    49     psLogMsg ("psphot", 3, "moments: %f sec\n", psTimerMark ("psphot"));
    50     psLogMsg ("psphot", 3, "%d sources\n", sources->n);
    51 
     49    psLogMsg ("psphot", 3, "%d moments: %f sec\n", sources->n, psTimerMark ("psphot"));
    5250    return (sources);
    5351}
Note: See TracChangeset for help on using the changeset viewer.