IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5672


Ignore:
Timestamp:
Dec 4, 2005, 5:33:32 AM (20 years ago)
Author:
eugene
Message:

simultaneous PSF fitting, better background, basic deblending

Location:
trunk/psphot
Files:
7 added
14 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/Makefile

    r5653 r5672  
    2929$(SRC)/psphotApplyPSF.$(ARCH).o    \
    3030$(SRC)/psphotFixedPSF.$(ARCH).o    \
     31$(SRC)/psphotEnsemblePSF.$(ARCH).o  \
     32$(SRC)/psphotReapplyPSF.$(ARCH).o  \
    3133$(SRC)/psphotFitGalaxies.$(ARCH).o \
    3234$(SRC)/psphotOutput.$(ARCH).o      \
    … …  
    3638$(SRC)/psphotDefinePixels.$(ARCH).o\
    3739$(SRC)/psphotMagnitudes.$(ARCH).o  \
     40$(SRC)/psphotImageBackground.$(ARCH).o  \
    3841$(SRC)/psLine.$(ARCH).o            \
    3942$(SRC)/psModulesUtils.$(ARCH).o    \
    4043$(SRC)/pmPeaksSigmaLimit.$(ARCH).o \
    4144$(SRC)/pmSourceFitFixed.$(ARCH).o  \
     45$(SRC)/psSparse.$(ARCH).o            \
    4246$(SRC)/psImageData.$(ARCH).o
    4347
    … …  
    6973TEST = \
    7074$(SRC)/psphotTest.$(ARCH).o \
     75$(SRC)/psphotTestArguments.$(ARCH).o \
     76$(SRC)/pmSourceContour.$(ARCH).o \
     77$(SRC)/psImageData.$(ARCH).o        \
     78$(SRC)/psphotSetup.$(ARCH).o        \
     79$(SRC)/psModulesUtils.$(ARCH).o     \
    7180$(SRC)/psSparse.$(ARCH).o
    7281
    … …  
    8594$(TEST): $(SRC)/psphot.h
    8695
    87 INSTALL = psphot psphotModelTest
     96INSTALL = psphot psphotTest psphotModelTest
    8897
    8998# dependancy rules for binary code #########################

  • trunk/psphot/src/modelTestFitSource.c

    r5622 r5672  
    77    bool status;
    88    int modelType;
    9     float sky, obsMag, fitMag, value;
     9    float obsMag, fitMag, value;
    1010    char name[64];
    1111
    … …  
    9696
    9797    // subtract object, leave local sky
    98     sky = model->params->data.F32[0];
    99     model->params->data.F32[0] = 0;
    100     pmSourceSubModel (source->pixels, source->mask, model, false);
    101     model->params->data.F32[0] = sky;
     98    pmSourceSubModel (source->pixels, source->mask, model, false, false);
    10299   
    103100    // write out

  • trunk/psphot/src/pmPeaksSigmaLimit.c

    r5617 r5672  
    2323    // set peak threshold
    2424    NSIGMA    = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT");
    25     threshold = NSIGMA*sky->sampleStdev + sky->sampleMean;
     25   
     26    // threshold = NSIGMA*sky->sampleStdev + sky->sampleMean;
     27    threshold = NSIGMA*sky->sampleStdev;
    2628    psLogMsg ("psphot", 3, "threshold: %f DN\n", threshold);
    2729

  • trunk/psphot/src/psphot.c

    r5643 r5672  
    1111    pmPSFClump   psfClump;
    1212    bool         status;
     13
     14    psphotModelGroupInit ();
    1315
    1416    config = psphotArguments (&argc, argv);
    … …  
    2527    // measure image stats for initial guess
    2628    sky = psphotImageStats (imdata, config);
     29
     30    // psPolynomial2D *skyModel = psphotImageBackground (imdata, config, sky);
     31    psphotImageBackground (imdata, config, sky);
    2732
    2833    // find the peaks in the image
    … …  
    5358
    5459    switch (FITMODE) {
     60      case -2:
     61        psphotEnsemblePSF (imdata, config, sources, psf, sky);
     62        psphotReapplyPSF (imdata, config, sources, psf, sky);
     63        break;
     64
     65      case -1:
     66        psphotEnsemblePSF (imdata, config, sources, psf, sky);
     67        break;
     68
    5569      case 0:
    5670        psphotFixedPSF (imdata, config, sources, psf, sky);
    … …  
    6377
    6478      case 2:
    65       // fit extended objects with galaxy models
     79        // fit extended objects with galaxy models
    6680        psphotApplyPSF (imdata, config, sources, psf, sky);
    6781        psphotFitGalaxies (imdata, config, sources, sky);
    … …  
    6983
    7084      case 3:
    71       // fit extended objects with galaxy models
     85        // fit extended objects with galaxy models
    7286        psphotFixedPSF (imdata, config, sources, psf, sky);
    7387        psphotFitGalaxies (imdata, config, sources, sky);

  • trunk/psphot/src/psphot.h

    r5654 r5672  
    33# include <unistd.h>   // for unlink
    44# include <pslib.h>
     5# include <pmObjects.h>
     6# include <pmPSF.h>
     7# include <pmPSFtry.h>
     8# include <pmModelGroup.h>
    59# include "psLibUtils.h"
    6 # include "pmObjects_EAM.h"
    710# include "psModulesUtils.h"
    8 # include "pmPSF.h"
    9 # include "pmPSFtry.h"
    10 # include "pmModelGroup.h"
    1111# include "psSparse.h"
     12
     13# define psMemCopy(A)(psMemIncrRefCounter((A)))
    1214
    1315typedef struct {
    … …  
    3234bool         psphotSubtractPSF (pmSource *source);
    3335int          psphotSortBySN (const void **a, const void **b);
     36int          psphotSortByY (const void **a, const void **b);
    3437int          psphotSaveImage (psMetadata *header, psImage *image, char *filename);
    3538bool         psphotDefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius);
    … …  
    6164psMetadata *modelTestArguments (int *argc, char **argv);
    6265bool modelTestFitSource (eamReadout *imdata, psMetadata *config);
     66bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
     67float psphotCrossProduct (pmSource *Mi, pmSource *Mj);
     68psPolynomial2D *psphotImageBackground (eamReadout *imdata, psMetadata *config, psStats *sky);
     69bool psphotReapplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
     70
     71pmModel *pmModelCopy (pmModel *model);
     72psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold);
     73psMetadata *psphotTestArguments (int *argc, char **argv);

  • trunk/psphot/src/psphotFitGalaxies.c

    r5593 r5672  
    66    float x;
    77    float y;
    8     float sky;
    98    int   Nfit = 0;
    109    int   Nfail = 0;
    … …  
    104103
    105104        // subtract object, leave local sky
    106         sky = model->params->data.F32[0];
    107         model->params->data.F32[0] = 0;
    108         pmSourceSubModel (source->pixels, source->mask, model, false);
    109         model->params->data.F32[0] = sky;
    110 
     105        pmSourceSubModel (source->pixels, source->mask, model, false, false);
    111106    }
    112107    psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d failures, %d total iterations)\n", psTimerMark ("psphot"), Nfit, Nfail, Niter);

  • trunk/psphot/src/psphotFixedPSF.c

    r5643 r5672  
    88    float x;
    99    float y;
    10     float Sky;
    1110    int   Nfit = 0;
    1211
    … …  
    7776
    7877        // subtract object, leave local sky
    79         Sky = model->params->data.F32[0];
    80         model->params->data.F32[0] = 0;
    81         pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false);
    82         model->params->data.F32[0] = Sky;
     78        pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
    8379    }
    8480

  • trunk/psphot/src/psphotMagnitudes.c

    r5124 r5672  
    2424
    2525    // replace source flux
    26     // XXX EAM : need to add 'sky' boolean option here
    27     model->params->data.F32[0] = 0;
    28     pmSourceAddModel (source->pixels, source->mask, model, false);
    29     model->params->data.F32[0] = sky;
     26    pmSourceAddModel (source->pixels, source->mask, model, false, false);
    3027
    3128    // measure object photometry
    … …  
    3835
    3936    // subtract object, leave local sky
    40     // XXX EAM : need to add 'sky' boolean option here
    41     model->params->data.F32[0] = 0;
    42     pmSourceSubModel (source->pixels, source->mask, model, false);
    43     model->params->data.F32[0] = sky;
     37    pmSourceSubModel (source->pixels, source->mask, model, false, false);
    4438
    4539    // unmask aperture

  • trunk/psphot/src/psphotMarkPSF.c

    r5058 r5672  
    3030    float nSx, nSy, Chi;
    3131
    32     if (source->modelPSF == NULL) return (false);
     32    pmModel *model = source->modelPSF;
    3333
    34     // if object has fitted peak above saturation, label as SATSTAR
     34    // do we actually have a valid PSF model?
     35    if (model == NULL) return (false);
     36
     37    // did the model fit fail for one or another reason?
     38    switch (model->status) {
     39      case PM_MODEL_BADARGS:
     40      case PM_MODEL_UNTRIED:
     41        source->type = PM_SOURCE_FAIL_FIT_PSF;  // better choice?
     42        return false;
     43      case PM_MODEL_SUCCESS:
     44        break;
     45      case PM_MODEL_NONCONVERGE:
     46      case PM_MODEL_OFFIMAGE:
     47        psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", model->params->data.F32[2], model->params->data.F32[3], model->nIter);
     48        source->type = PM_SOURCE_FAIL_FIT_PSF;  // better choice?
     49        return false;
     50      default:
     51        return false;
     52    }
     53
     54    // if the object has fitted peak above saturation, label as SATSTAR
     55    // this is a valid PSF object, but ignore the other quality tests
    3556    // remember: fit does not use saturated pixels (masked)
    36     if (source->modelPSF->params->data.F32[1] >= SATURATE) {
     57    if (model->params->data.F32[1] >= SATURATE) {
    3758        if (source->type == PM_SOURCE_PSFSTAR) {
    3859            psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n");
    … …  
    4162        return (true);
    4263    }
     64
     65    // if the object has a fitted peak below 0, the fit did not converge cleanly
     66    if (model->params->data.F32[1] < 0) {
     67        source->type = PM_SOURCE_FAIL_FIT_PSF;
     68        return (false);
     69    }
     70
     71    // if the source was predicted to be a SATSTAR, but it fitted below saturation,
     72    // make a note to the user
    4373    if (source->type == PM_SOURCE_SATSTAR) {
    4474        psLogMsg ("psphot", 5, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n");
    … …  
    4676    }
    4777
    48     SN  = source->modelPSF->params->data.F32[1]/source->modelPSF->dparams->data.F32[1];
    49     SX  = source->modelPSF->params->data.F32[4];
    50     SY  = source->modelPSF->params->data.F32[5];
    51     dSX = source->modelPSF->dparams->data.F32[4];
    52     dSY = source->modelPSF->dparams->data.F32[5];
    53     Chi = source->modelPSF->chisq / source->modelPSF->nDOF;
     78    SN  = model->params->data.F32[1]/model->dparams->data.F32[1];
     79    SX  = model->params->data.F32[4];
     80    SY  = model->params->data.F32[5];
     81    dSX = model->dparams->data.F32[4];
     82    dSY = model->dparams->data.F32[5];
     83    Chi = model->chisq / model->nDOF;
    5484
    5585    // swing of sigma_x,y in sigmas

  • trunk/psphot/src/psphotOutput.c

    r5628 r5672  
    583583    unlink (filename);
    584584    psFits *fits = psFitsAlloc (filename);
    585     psFitsWriteImage (fits, header, image, 0);
     585    psFitsWriteImage (fits, NULL, image, 0);
    586586    psFree (fits);
    587587    return (TRUE);

  • trunk/psphot/src/psphotSetup.c

    r5622 r5672  
    8484    // XXX EAM does the mask need to grow?
    8585    float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION");
    86     for (int i = 0; i < image->numRows; i++) {
     86    for (int i = 0; status && (i < image->numRows); i++) {
    8787        for (int j = 0; j < image->numCols; j++) {
    8888            if (image->data.F32[i][j] >= SATURATION) {
    8989                mask->data.U8[i][j] |= PSPHOT_MASK_SATURATED;
     90            }
     91        }
     92    }
     93
     94    // mask the pixels below min threshold
     95    float MIN_VALID = psMetadataLookupF32 (&status, config, "MIN_VALID_PIXEL");
     96    for (int i = 0; status && (i < image->numRows); i++) {
     97        for (int j = 0; j < image->numCols; j++) {
     98            if (image->data.F32[i][j] < MIN_VALID) {
     99                mask->data.U8[i][j] |= PSPHOT_MASK_INVALID;
    90100            }
    91101        }

  • trunk/psphot/src/psphotSortBySN.c

    r4949 r5672  
    1818}
    1919
     20// sort by Y (ascending)
     21int psphotSortByY (const void **a, const void **b)
     22{
     23    pmSource *A = *(pmSource **)a;
     24    pmSource *B = *(pmSource **)b;
     25
     26    psF32 fA = (A->moments == NULL) ? 0 : A->moments->y;
     27    psF32 fB = (B->moments == NULL) ? 0 : B->moments->y;
     28
     29    psF32 diff = fA - fB;
     30    if (diff > FLT_EPSILON) return (+1);
     31    if (diff < FLT_EPSILON) return (-1);
     32    return (0);
     33}
     34

  • trunk/psphot/src/psphotSubtractPSF.c

    r5049 r5672  
    66bool psphotSubtractPSF (pmSource *source)
    77{
    8   float sky;
    98  pmModel *model;
    10   psImage *pixels;
    119
    1210  // only subtract successful fits
    … …  
    2321  }         
    2422
    25   pixels = source->pixels;
    26 
    2723  // subtract object, leave local sky
    28   sky = model->params->data.F32[0];
    29   model->params->data.F32[0] = 0;
    30   pmSourceSubModel (pixels, source->mask, model, false);
    31   model->params->data.F32[0] = sky;
     24  pmSourceSubModel (source->pixels, source->mask, model, false, false);
    3225
    3326  // XXX EAM : amplify the weight matrix a la dophot?

  • trunk/psphot/src/psphotTest.c

    r5654 r5672  
    33int main (int argc, char **argv) {
    44
    5   psSparseMatrixTest ();
     5    bool status;
    66
     7    psMetadata *config = psphotTestArguments (&argc, argv);
     8    eamReadout *imdata = psphotSetup (config);
     9
     10    int x = psMetadataLookupF32 (&status, config, "TEST_FIT_X");
     11    int y = psMetadataLookupF32 (&status, config, "TEST_FIT_Y");
     12    float fraction = psMetadataLookupF32 (&status, config, "TEST_FRACTION");
     13    float threshold = psMetadataLookupF32 (&status, config, "TEST_THRESHOLD");
     14
     15    if (threshold == -1) {
     16        threshold = fraction * imdata->image->data.F32[y][x];
     17    }
     18
     19    psArray *contour = pmSourceContour_EAM (imdata->image, x, y, threshold);
     20    if (contour == NULL) {
     21        fprintf (stderr, "error: invalid contour\n");
     22        exit (1);
     23    }
     24
     25    psVector *xV = contour->data[0];
     26    psVector *yV = contour->data[1];
     27
     28    for (int i = 0; i < xV->n; i++) {
     29        fprintf (stdout, "%d %f %f\n", i, xV->data.F32[i], yV->data.F32[i]);
     30    }
     31
     32    // psSparseMatrixTest ();
     33    exit (0);
    734}
Note: See TracChangeset for help on using the changeset viewer.