IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5828


Ignore:
Timestamp:
Dec 22, 2005, 4:20:27 PM (20 years ago)
Author:
eugene
Message:

various cleanups, fixes, additions...

Location:
trunk/psphot
Files:
5 added
1 deleted
24 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/Makefile

    r5802 r5828  
    2121
    2222PSPHOT = \
    23 $(SRC)/psphot.$(ARCH).o            \
    24 $(SRC)/psphotArguments.$(ARCH).o   \
    25 $(SRC)/psphotSetup.$(ARCH).o       \
    26 $(SRC)/psphotImageStats.$(ARCH).o  \
    27 $(SRC)/psphotSourceStats.$(ARCH).o \
    28 $(SRC)/psphotChoosePSF.$(ARCH).o   \
    29 $(SRC)/psphotApplyPSF.$(ARCH).o    \
    30 $(SRC)/psphotFixedPSF.$(ARCH).o    \
     23$(SRC)/psphot.$(ARCH).o             \
     24$(SRC)/psphotArguments.$(ARCH).o    \
     25$(SRC)/psphotSetup.$(ARCH).o        \
     26$(SRC)/psphotImageStats.$(ARCH).o   \
     27$(SRC)/psphotSourceStats.$(ARCH).o  \
     28$(SRC)/psphotChoosePSF.$(ARCH).o    \
     29$(SRC)/psphotApplyPSF.$(ARCH).o     \
     30$(SRC)/psphotFixedPSF.$(ARCH).o     \
    3131$(SRC)/psphotEnsemblePSF.$(ARCH).o  \
    32 $(SRC)/psphotReapplyPSF.$(ARCH).o  \
    33 $(SRC)/psphotFitGalaxies.$(ARCH).o \
    34 $(SRC)/psphotOutput.$(ARCH).o      \
    35 $(SRC)/psphotMarkPSF.$(ARCH).o     \
    36 $(SRC)/psphotSortBySN.$(ARCH).o    \
    37 $(SRC)/psphotDefinePixels.$(ARCH).o\
    38 $(SRC)/psphotMagnitudes.$(ARCH).o  \
    39 $(SRC)/psphotImageBackground.$(ARCH).o  \
    40 $(SRC)/psphotBasicDeblend.$(ARCH).o  \
     32$(SRC)/psphotReapplyPSF.$(ARCH).o   \
     33$(SRC)/psphotFitGalaxies.$(ARCH).o  \
     34$(SRC)/psphotOutput.$(ARCH).o       \
     35$(SRC)/psphotMarkPSF.$(ARCH).o      \
     36$(SRC)/psphotSortBySN.$(ARCH).o     \
     37$(SRC)/psphotDefinePixels.$(ARCH).o \
     38$(SRC)/psphotMagnitudes.$(ARCH).o   \
     39$(SRC)/psphotRadiusChecks.$(ARCH).o \
     40$(SRC)/psphotReplaceUnfit.$(ARCH).o \
     41$(SRC)/psphotEvalPSF.$(ARCH).o      \
     42$(SRC)/psphotEvalFLT.$(ARCH).o      \
     43$(SRC)/psphotFullFit.$(ARCH).o      \
     44$(SRC)/pmSourceContour.$(ARCH).o    \
     45$(SRC)/psLine.$(ARCH).o             \
     46$(SRC)/psModulesUtils.$(ARCH).o     \
     47$(SRC)/pmPeaksSigmaLimit.$(ARCH).o  \
     48$(SRC)/pmSourceFitFixed.$(ARCH).o   \
     49$(SRC)/psSparse.$(ARCH).o           \
     50$(SRC)/psImageData.$(ARCH).o        \
     51$(SRC)/psphotModelTest.$(ARCH).o    \
     52$(SRC)/psphotImageBackground.$(ARCH).o \
     53$(SRC)/psphotBasicDeblend.$(ARCH).o    \
    4154$(SRC)/psphotModelGroupInit.$(ARCH).o  \
    42 $(SRC)/psphotRadiusChecks.$(ARCH).o  \
    43 $(SRC)/psphotReplaceUnfit.$(ARCH).o  \
    44 $(SRC)/psphotEvalPSF.$(ARCH).o  \
    45 $(SRC)/psphotEvalFLT.$(ARCH).o  \
    46 $(SRC)/psphotFullFit.$(ARCH).o  \
    47 $(SRC)/pmSourceContour.$(ARCH).o \
    48 $(SRC)/psLine.$(ARCH).o            \
    49 $(SRC)/psModulesUtils.$(ARCH).o    \
    50 $(SRC)/pmPeaksSigmaLimit.$(ARCH).o \
    51 $(SRC)/pmSourceFitFixed.$(ARCH).o  \
    52 $(SRC)/psSparse.$(ARCH).o            \
    53 $(SRC)/psImageData.$(ARCH).o        \
    54 $(SRC)/psphotApResid.$(ARCH).o
    55 
    56 PSMODULES = \
    57 $(SRC)/psEllipse.$(ARCH).o         \
    58 $(SRC)/pmPSF.$(ARCH).o             \
    59 $(SRC)/pmPSFtry.$(ARCH).o          \
    60 $(SRC)/pmModelGroup.$(ARCH).o      \
    61 $(SRC)/pmObjects_EAM.$(ARCH).o
     55$(SRC)/pmModelFitSet.$(ARCH).o  \
     56$(SRC)/pmSourceFitSet.$(ARCH).o        \
     57$(SRC)/psphotSourceFits.$(ARCH).o        \
     58$(SRC)/psphotBlendFit.$(ARCH).o        \
     59$(SRC)/psphotApResid.$(ARCH).o \
     60$(SRC)/psBicube.$(ARCH).o
    6261
    6362MODELS = \
    … …  
    9190psphot: $(BIN)/psphot.$(ARCH)
    9291$(BIN)/psphot.$(ARCH) : $(PSPHOT)
    93 $(PSPHOT) $(PSMODULES): $(SRC)/psphot.h
     92$(PSPHOT) : $(SRC)/psphot.h
    9493
    9594modeltest.install: psphotModelTest.install

  • trunk/psphot/src/pmPeaksSigmaLimit.c

    r5772 r5828  
    2626    threshold = NSIGMA*sky->sampleStdev + sky->sampleMean;
    2727    // threshold = NSIGMA*sky->sampleStdev;
    28     psLogMsg ("psphot", 3, "threshold: %f DN\n", threshold);
     28    psLogMsg ("psphot", 4, "threshold: %f DN\n", threshold);
    2929
    3030    // find the peaks in the smoothed image

  • trunk/psphot/src/psModulesUtils.c

    r5802 r5828  
    110110    pmModel *new = pmModelAlloc (model->type);
    111111   
    112     new->chisq = model->chisq;
    113     new->nIter = model->nIter;
     112    new->chisq  = model->chisq;
     113    new->nDOF   = model->nDOF;
     114    new->nIter  = model->nIter;
     115    new->status = model->status;
     116    new->radius = model->radius;
    114117
    115118    for (int i = 0; i < new->params->n; i++) {

  • trunk/psphot/src/psphot.c

    r5802 r5828  
    1212    bool         status;
    1313
     14    psTimerStart ("complete");
     15
    1416    psphotModelGroupInit ();
    1517
    1618    config = psphotArguments (&argc, argv);
    1719
    18     // load input data (image and config)
    19     // create or load mask and weight images
    20     // we have memory leaks here -- may be from psMetadata
     20    // load input data (config and images (signal, noise, mask)
     21    // XXX we have memory leaks here -- may be from psMetadata
    2122    imdata = psphotSetup (config);
     23
    2224    char *breakPt = psMetadataLookupPtr (&status, config, "BREAK_POINT");
    23     if (!status) {
    24       breakPt = psStringCopy ("NONE");
     25    if (!status) breakPt = psStringCopy ("NONE");
     26
     27    int FITMODE = psMetadataLookupS32 (&status, config, "FIT_MODE");
     28    if (!status) FITMODE = 2;
     29
     30    // run model fitting tests on a single source
     31    if (psMetadataLookupBool (&status, config, "TEST_FIT")) {
     32        psphotModelTest (imdata, config);
     33        exit (0);
    2534    }
    2635
    … …  
    2837    sky = psphotImageStats (imdata, config);
    2938
    30     // psPolynomial2D *skyModel = psphotImageBackground (imdata, config, sky);
     39    // generate a background model (currently, 2D polynomial)
     40    // XXX this should be available to be re-added to the original image
    3141    psphotImageBackground (imdata, config, sky);
    3242
    … …  
    5262    // source analysis is done in S/N order (brightest first)
    5363    sources = psArraySort (sources, psphotSortBySN);
    54 
    5564    psphotDumpMoments (config, sources);
    5665    if (!strcasecmp (breakPt, "CLASS")) exit (0);
    … …  
    5968    psf = psphotChoosePSF (config, sources, sky);
    6069    if (!strcasecmp (breakPt, "PSFFIT")) exit (0);
    61 
    62     // XXX add this as conditional output
    63     psphotSamplePSFs (psf, imdata->image);
    64 
    65     int FITMODE = psMetadataLookupS32 (&status, config, "FIT_MODE");
    66     if (!status) FITMODE = 2;
    6770
    6871    switch (FITMODE) {
    … …  
    7982
    8083      case 2:
    81         // fit extended objects with galaxy models
    82         psphotApplyPSF (imdata, config, sources, psf, sky);
    83         psphotFitGalaxies (imdata, config, sources, sky);
     84        psphotEnsemblePSF (imdata, config, sources, psf, sky);
     85        psphotBlendFit (imdata, config, sources, psf, sky);
     86        psphotReplaceUnfit (sources);
     87        psphotApResid (sources, config, psf);
    8488        break;
    8589
    … …  
    9195      case 4:
    9296        // fit extended objects with galaxy models
    93         psphotFixedPSF (imdata, config, sources, psf, sky);
     97        psphotApplyPSF (imdata, config, sources, psf, sky);
    9498        psphotFitGalaxies (imdata, config, sources, sky);
    9599        break;
    … …  
    97101
    98102    // write out data in appropriate format
    99     // measure aperture correction
    100     // psphotApertureCorrection (sources, config, sky);
     103    psphotSamplePSFs (config, psf, imdata->image);
    101104    psphotOutput (imdata, config, sources, psf, sky);
     105    psLogMsg ("psphot", 3, "wrote output: %f sec\n", psTimerMark ("psphot"));
     106    psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete"));
    102107    exit (0);
    103108}

  • trunk/psphot/src/psphot.h

    r5802 r5828  
    6363
    6464// psphotModelTest functions
    65 psMetadata     *modelTestArguments (int *argc, char **argv);
    66 bool            modelTestFitSource (eamReadout *imdata, psMetadata *config);
    6765bool            psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
    6866float           psphotCrossProduct (pmSource *Mi, pmSource *Mj);
    … …  
    7270pmModel        *pmModelCopy (pmModel *model);
    7371psArray        *pmSourceContour_EAM (psImage *image, int x, int y, float threshold);
    74 psMetadata     *psphotTestArguments (int *argc, char **argv);
    7572bool            psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky);
    7673
    7774bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
    7875bool psphotInitLimitsPSF (psMetadata *config);
    79 bool psphotEvalPSF (pmSource *source);
    80 bool psphotEvalFLT (pmSource *source);
     76bool psphotEvalPSF (pmSource *source, pmModel *model);
     77bool psphotEvalDBL (pmSource *source, pmModel *model);
     78bool psphotEvalFLT (pmSource *source, pmModel *model);
    8179bool psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type);
    82 bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source);
     80bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source, pmModel *model);
    8381bool psphotInitRadiusFLT (psMetadata *config, psStats *sky, pmModelType type);
    84 bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source);
    85 bool psphotSamplePSFs (pmPSF *psf, psImage *image);
     82bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source, pmModel *model);
     83bool psphotSamplePSFs (psMetadata *config, pmPSF *psf, psImage *image);
    8684bool psphotReplaceUnfit (psArray *sources);
    8785bool psphotDumpMoments (psMetadata *config, psArray *sources);
    8886bool psphotApResid (psArray *sources, psMetadata *config, pmPSF *psf);
     87bool psphotWritePSF (pmPSF *psf, char *filename);
     88pmPSF *psphotReadPSF (char *filename);
     89
     90bool psPolynomial2DtoMD (psMetadata *md, psPolynomial2D *poly, char *format, ...);
     91bool psPolynomial3DtoMD (psMetadata *md, psPolynomial3D *poly, char *format, ...);
     92psPolynomial2D *psPolynomial2DfromMD (psMetadata *folder);
     93psPolynomial3D *psPolynomial3DfromMD (psMetadata *folder);
     94bool psphotModelTest (eamReadout *imdata, psMetadata *config);
     95
     96bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF);
     97bool pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF);
     98psF32 pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x);
     99bool pmModelFitSetInit (pmModelType type);
     100
     101bool psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
     102
     103bool psphotInitLimitsFLT (psMetadata *config, psStats *sky);
     104
     105bool psphotFitPSF (eamReadout *imdata, pmSource *source);
     106bool psphotFitBlend (eamReadout *imdata, pmSource *source);
     107bool psphotFitBlob (eamReadout *imdata, pmSource *source, psArray *sources);
     108
     109pmModel *psphotFitFLT (eamReadout *imdata, pmSource *source);
     110psArray *psphotFitDBL (eamReadout *imdata, pmSource *source);
     111
     112psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y);
     113psPlane psImageBicubeMin (psPolynomial2D *poly);

  • trunk/psphot/src/psphotApResid.c

    r5802 r5828  
    2121    psVector *rflux   = psVectorAlloc (300, PS_TYPE_F64);
    2222    psVector *apResid = psVectorAlloc (300, PS_TYPE_F64);
    23     mask->n = rflux->n = apResid->n = 0;
     23    mask->n = xPos->n = yPos->n = rflux->n = apResid->n = 0;
    2424    Npsf = 0;
    2525
    … …  
    4343
    4444        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
    45         // XXX EAM : add in source flux
    4645        status = pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask);
    4746        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
     47
     48        pmSourceSubModel (source->pixels, source->mask, model, false, false);
    4849        if (!status) continue;
    4950
    … …  
    5556
    5657        psVectorExtend (mask, 100, 1);
     58        psVectorExtend (xPos, 100, 1);
     59        psVectorExtend (yPos, 100, 1);
    5760        psVectorExtend (rflux, 100, 1);
    5861        psVectorExtend (apResid, 100, 1);
    5962        Npsf ++;
    60 
    61         pmSourceSubModel (source->pixels, source->mask, model, false, false);
    62         psMemCheckCorruption (true);
    6363    }
     64    psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec\n", psTimerMark ("psphot"));
    6465
    6566    // 3hi/1lo sigma clipping on the rflux vs metric fit
    … …  
    7071
    7172    // linear clipped fit of apResid to rflux, xPos, yPos
    72     psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, rflux, xPos, yPos);
     73    psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux);
    7374    psf->skyBias  = psf->ApTrend->coeff[0][0][1] / (M_PI * PS_SQR(RADIUS));
    7475    psf->ApResid  = psf->ApTrend->coeff[0][0][0];
    7576    psf->dApResid = stats->sampleStdev;
     77    psf->ApTrend->coeff[0][0][1] = 0;
    7678
    77     psLogMsg ("ApResid", 4, "aperture residual: %f +/- %f : %f bias\n", psf->ApResid, psf->dApResid, psf->skyBias);
     79    /*
     80      (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
     81      (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
     82      (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
     83    */
     84
     85    # if (0)
     86    psPolynomial3D *poly = psf->ApTrend;
     87    for (int nz = 0; nz <= poly->nZ; nz++) {
     88        for (int ny = 0; ny <= poly->nY; ny++) {
     89            for (int nx = 0; nx <= poly->nX; nx++) {
     90                fprintf (stderr, "%d %d %d : %22.15g\n", nx, ny, nz, poly->coeff[nx][ny][nz]);
     91            }
     92        }
     93    }
     94    # endif
     95
     96    psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot"));
     97    psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias\n", psf->ApResid, psf->dApResid, psf->skyBias);
    7898
    7999    psFree (stats);

  • trunk/psphot/src/psphotApplyPSF.c

    r5772 r5828  
    3434
    3535        // sets the model radius (via source->model) and adjusts pixels as needed
    36         psphotCheckRadiusPSF (imdata, source);
     36        psphotCheckRadiusPSF (imdata, source, model);
    3737
    3838        x = model->params->data.F32[2];
    … …  
    4848
    4949        // check if model fit is acceptable
    50         if (psphotEvalPSF (source)) {
     50        if (psphotEvalPSF (source, source->modelPSF)) {
    5151            pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
    52             source->mode |= PM_SOURCE_SUBTRACTED;
     52            source->mode |=  PM_SOURCE_SUBTRACTED;
     53            source->mode &= ~PM_SOURCE_TEMPSUB;
    5354            Nsub ++;
    5455        }

  • trunk/psphot/src/psphotArguments.c

    r5607 r5828  
    55psMetadata *psphotArguments (int *argc, char **argv) {
    66
    7   int N;
    8   unsigned int Nfail;
    9   int mode = PS_DATA_STRING | PS_META_REPLACE;
     7    int N;
     8    unsigned int Nfail;
     9    int mode = PS_DATA_STRING | PS_META_REPLACE;
    1010
    11   // basic pslib options
    12   fprintf (stderr, "starting... %s\n", psLibVersion());
    13   psLogSetFormat ("M");
    14   psArgumentVerbosity (argc, argv);
     11    // basic pslib options
     12    psLogSetFormat ("M");
     13    psArgumentVerbosity (argc, argv);
    1514
    16   // optional mask image - add to config
    17   char *mask = NULL;
    18   if ((N = psArgumentGet (*argc, argv, "-mask"))) {
    19     psArgumentRemove (N, argc, argv);
    20     mask = psStringCopy (argv[N]);
    21     psArgumentRemove (N, argc, argv);
    22   }
     15    // optional mask image - add to config
     16    char *mask = NULL;
     17    if ((N = psArgumentGet (*argc, argv, "-mask"))) {
     18        psArgumentRemove (N, argc, argv);
     19        mask = psStringCopy (argv[N]);
     20        psArgumentRemove (N, argc, argv);
     21    }
    2322
    24   // optional weight image - add to config
    25   char *weight = NULL;
    26   if ((N = psArgumentGet (*argc, argv, "-weight"))) {
    27     psArgumentRemove (N, argc, argv);
    28     weight = psStringCopy (argv[N]);
    29     psArgumentRemove (N, argc, argv);
    30   }
     23    // optional weight image - add to config
     24    char *weight = NULL;
     25    if ((N = psArgumentGet (*argc, argv, "-weight"))) {
     26        psArgumentRemove (N, argc, argv);
     27        weight = psStringCopy (argv[N]);
     28        psArgumentRemove (N, argc, argv);
     29    }
    3130
    32   // optional output residual image - add to config
    33   char *resid = NULL;
    34   if ((N = psArgumentGet (*argc, argv, "-resid"))) {
    35     psArgumentRemove (N, argc, argv);
    36     resid = psStringCopy (argv[N]);
    37     psArgumentRemove (N, argc, argv);
    38   }
     31    // optional output residual image - add to config
     32    char *resid = NULL;
     33    if ((N = psArgumentGet (*argc, argv, "-resid"))) {
     34        psArgumentRemove (N, argc, argv);
     35        resid = psStringCopy (argv[N]);
     36        psArgumentRemove (N, argc, argv);
     37    }
    3938
    40   // optional output residual image - add to config
    41   char *photcode = NULL;
    42   if ((N = psArgumentGet (*argc, argv, "-photcode"))) {
    43     psArgumentRemove (N, argc, argv);
    44     photcode = psStringCopy (argv[N]);
    45     psArgumentRemove (N, argc, argv);
    46   }
     39    // optional analysis region - add to config
     40    char *region = NULL;
     41    if ((N = psArgumentGet (*argc, argv, "-region"))) {
     42        psArgumentRemove (N, argc, argv);
     43        region = psStringCopy (argv[N]);
     44        psArgumentRemove (N, argc, argv);
     45    }
    4746
    48   if (*argc != 4) usage ();
     47    // optional output residual image - add to config
     48    char *photcode = NULL;
     49    if ((N = psArgumentGet (*argc, argv, "-photcode"))) {
     50        psArgumentRemove (N, argc, argv);
     51        photcode = psStringCopy (argv[N]);
     52        psArgumentRemove (N, argc, argv);
     53    }
    4954
    50   // load config information
    51   psMetadata *config = psMetadataAlloc ();
    52   psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL);
    53   config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE);
     55    char *psffile = NULL;
     56    if ((N = psArgumentGet (*argc, argv, "-psf"))) {
     57        psArgumentRemove (N, argc, argv);
     58        psffile = psStringCopy (argv[N]);
     59        psArgumentRemove (N, argc, argv);
     60    }
    5461
    55   // identify input image & optional weight & mask images
    56   // command-line entries override config-file entries
    57   psMetadataAdd (config, PS_LIST_HEAD, "IMAGE",       mode, "", argv[1]);
    58   psMetadataAdd (config, PS_LIST_HEAD, "OUTPUT_FILE", mode, "", argv[2]);
     62    bool ModelTest = false;
     63    float ModelTest_X, ModelTest_Y;
     64    char *model = NULL;
     65    char *fitset = NULL;
     66    char *fitmode = NULL;
     67    if ((N = psArgumentGet (*argc, argv, "-modeltest"))) {
     68        ModelTest = true;
     69        psArgumentRemove (N, argc, argv);
     70        ModelTest_X = atof (argv[N]);
     71        psArgumentRemove (N, argc, argv);
     72        ModelTest_Y = atof (argv[N]);
     73        psArgumentRemove (N, argc, argv);
    5974
    60   if (mask != NULL) {
    61     psMetadataAdd (config, PS_LIST_HEAD, "MASK_IMAGE", mode, "", mask);
    62   }
    63   if (weight != NULL) {
    64     psMetadataAdd (config, PS_LIST_HEAD, "WEIGHT_IMAGE", mode, "", weight);
    65   }
    66   if (resid != NULL) {
    67     psMetadataAdd (config, PS_LIST_HEAD, "RESID_IMAGE", mode, "", resid);
    68   }
    69   if (photcode != NULL) {
    70     psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", photcode);
    71   } else {
    72     psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", "NONE");
    73   }   
    74   return (config);
     75        if ((N = psArgumentGet (*argc, argv, "-model"))) {
     76            psArgumentRemove (N, argc, argv);
     77            model = psStringCopy (argv[N]);
     78            psArgumentRemove (N, argc, argv);
     79        }
     80
     81        if ((N = psArgumentGet (*argc, argv, "-fitmode"))) {
     82            psArgumentRemove (N, argc, argv);
     83            fitmode = psStringCopy (argv[N]);
     84            psArgumentRemove (N, argc, argv);
     85        }
     86        if ((N = psArgumentGet (*argc, argv, "-fitset"))) {
     87            psArgumentRemove (N, argc, argv);
     88            fitset = psStringCopy (argv[N]);
     89            psArgumentRemove (N, argc, argv);
     90        }
     91    }
     92
     93    if (*argc != 4) usage ();
     94
     95    // load config information
     96    psMetadata *config = psMetadataAlloc ();
     97    psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL);
     98    config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE);
     99
     100    // identify input image & optional weight & mask images
     101    // command-line entries override config-file entries
     102    psMetadataAdd (config, PS_LIST_HEAD, "IMAGE",       mode, "", argv[1]);
     103    psMetadataAdd (config, PS_LIST_HEAD, "OUTPUT_FILE", mode, "", argv[2]);
     104
     105    if (mask != NULL) {
     106        psMetadataAdd (config, PS_LIST_HEAD, "MASK_IMAGE", mode, "", mask);
     107    }
     108    if (weight != NULL) {
     109        psMetadataAdd (config, PS_LIST_HEAD, "WEIGHT_IMAGE", mode, "", weight);
     110    }
     111    if (resid != NULL) {
     112        psMetadataAdd (config, PS_LIST_HEAD, "RESID_IMAGE", mode, "", resid);
     113    }
     114    if (region != NULL) {
     115        psMetadataAdd (config, PS_LIST_HEAD, "ANALYSIS_REGION", mode, "", region);
     116    }
     117    if (photcode != NULL) {
     118        psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", photcode);
     119    } else {
     120        psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", "NONE");
     121    }   
     122    if (psffile != NULL) {
     123        psMetadataAdd (config, PS_LIST_HEAD, "PSF_INPUT_FILE", mode, "", psffile);
     124    }
     125
     126    // model related options
     127    if (ModelTest) {
     128        psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT", PS_DATA_BOOL, "", true);
     129
     130        int fmode = PS_DATA_F32 | PS_META_REPLACE;
     131        psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_X", fmode, "", ModelTest_X);
     132        psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_Y", fmode, "", ModelTest_Y);
     133
     134        if (model != NULL) {
     135            psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_MODEL", mode, "", model);
     136        }
     137
     138        if (fitmode != NULL) {
     139            psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_MODE", mode, "", fitmode);
     140        }
     141        if (fitset != NULL) {
     142            psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_SET", mode, "", fitset);
     143        }
     144    }
     145    return (config);
    75146}
    76147

  • trunk/psphot/src/psphotBasicDeblend.c

    r5772 r5828  
    88    pmSource *source, *testSource;
    99
    10     FILE *f = fopen ("deblend.dat", "w");
    11     if (f == NULL) psAbort ("psphot", "can't open deblend.dat output file");
    1210    int Nblend = 0;
    1311
    … …  
    8280            // XXX EAM : should the contour input coordinate be in parent or subimage coords? parent, for now
    8381            psArray *contour = pmSourceContour_EAM (source->pixels, source->peak->x, source->peak->y, threshold);
    84             if (contour == NULL) {
    85                 fprintf (stderr, "below threshold? invalid peak?\n");
    86                 continue;
    87             }
     82            if (contour == NULL) continue;
    8883
    8984            // the source contour consists of two vectors, xv and yv.  the contour is
    … …  
    10297                    if (xv->data.F32[j+1] < testSource->peak->x) break;
    10398
     99                    # if (0)
    104100                    int xp0 = source->moments->x - source->pixels->col0;
    105101                    int xp1 = source->peak->x - source->pixels->col0;
    … …  
    111107                    int yp2 = testSource->moments->y - testSource->pixels->row0;
    112108                    int yp3 = testSource->peak->y - testSource->pixels->row0;
     109                   
     110                    fprintf (f, "%d %d (%f, %f) :  %d %d (%f, %f)  vs %f\n",
     111                             source->peak->x, source->peak->y,
     112                             source->pixels->data.F32[yp0][xp0], source->pixels->data.F32[yp1][xp1],
     113                             testSource->peak->x, testSource->peak->y,
     114                             testSource->pixels->data.F32[yp2][xp2], testSource->pixels->data.F32[yp3][xp3], threshold
     115                       );
     116                    # endif
     117                   
     118                    testSource->mode |= PM_SOURCE_BLEND;
    113119
    114                     fprintf (f, "%d %d (%f, %f) :  %d %d (%f, %f)  vs %f\n",
    115                              source->peak->x, source->peak->y,
    116                              source->pixels->data.F32[yp0][xp0], source->pixels->data.F32[yp1][xp1],
    117                              testSource->peak->x, testSource->peak->y,
    118                              testSource->pixels->data.F32[yp2][xp2], testSource->pixels->data.F32[yp3][xp3], threshold
    119                              );
     120                    // add this to the list of source->blends
     121                    if (source->blends == NULL) {
     122                        source->blends = psArrayAlloc (16);
     123                        source->blends->n = 0;
     124                    }
     125                    psArrayAdd (source->blends, 16, testSource);
    120126
    121                     testSource->mode |= PM_SOURCE_BLEND;
    122127                    Nblend ++;
    123128                    j = xv->n;
    … …  
    126131        }
    127132    }
    128     fclose (f);
    129     psLogMsg ("psphot.deblend", 3, "identified %d blended objects (%f)\n", Nblend, psTimerMark ("psphot"));
     133    psLogMsg ("psphot.deblend", 3, "identified %d blended objects (%f sec)\n", Nblend, psTimerMark ("psphot"));
    130134    return true;
    131135}
    132 
    133 
    134 

  • trunk/psphot/src/psphotChoosePSF.c

    r5802 r5828  
    1010    psArray        *stars = NULL;
    1111    psMetadataItem *item  = NULL;
     12
     13    psTimerStart ("psphot");
    1214
    1315    // array to store candidate PSF stars
    … …  
    2325        if (source->mode & PM_SOURCE_PSFSTAR) psArrayAdd (stars, 200, source);
    2426    }
    25     psLogMsg ("psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n);
     27    psLogMsg ("psphot.pspsf", 4, "selected candidate %d PSF objects\n", stars->n);
    2628
    2729    // get the fixed PSF fit radius
    … …  
    8385
    8486    modelName = pmModelGetType (psf->type);
     87    psLogMsg ("psphot.pspsf", 3, "select psf model: %f sec\n", psTimerMark ("psphot"));
    8588    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
    8689

  • trunk/psphot/src/psphotEnsemblePSF.c

    r5802 r5828  
    11# include "psphot.h"
     2
     3psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y) {
     4
     5    int ix = x - image->col0;
     6    int iy = y - image->row0;
     7
     8    psF32 *Fm = &image->data.F32[iy - 1][ix];
     9    psF32 *Fo = &image->data.F32[iy + 0][ix];
     10    psF32 *Fp = &image->data.F32[iy + 1][ix];
     11
     12    double Fxm = Fm[-1] + Fo[-1] + Fp[-1];
     13    double Fxp = Fm[+1] + Fo[+1] + Fp[+1];
     14    double Fym = Fm[-1] + Fm[+0] + Fm[+1];
     15    double Fyp = Fp[-1] + Fp[+0] + Fp[+1];
     16    double Foo = Fym + Fyp + Fo[-1] + Fo[+0] + Fo[+1];
     17
     18    psPolynomial2D *poly = psPolynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD);
     19    poly->mask[2][2] = 1;
     20    poly->mask[1][2] = 1;
     21    poly->mask[2][1] = 1;
     22
     23    poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ;
     24
     25    poly->coeff[1][0] = (Fxp - Fxm)/6.0;
     26    poly->coeff[0][1] = (Fyp - Fym)/6.0;
     27   
     28    poly->coeff[2][0] = (Fxp + Fxm)/2.0 - Foo/3.0;
     29    poly->coeff[0][2] = (Fyp + Fym)/2.0 - Foo/3.0;
     30   
     31    poly->coeff[1][1] = (Fp[+1] + Fm[-1] - Fm[+1] - Fp[-1])/4.0;
     32   
     33    return (poly);
     34}
     35
     36psPlane psImageBicubeMin (psPolynomial2D *poly) {
     37
     38    psPlane min;
     39
     40    min.xErr = min.yErr = 0;
     41
     42    double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]);
     43
     44    min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det;
     45    min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det;
     46    return (min);
     47}
    248
    349bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) {
    … …  
    2773    index->n = 0;
    2874
     75    bool UseAnalysisRegion = false;
     76    psRegion AnalysisRegion;
     77    char *region = psMetadataLookupPtr (&status, config, "ANALYSIS_REGION");
     78    if (status) {
     79        UseAnalysisRegion = true;
     80        AnalysisRegion = psRegionFromString (region);
     81        psLogMsg ("psphotEnsemblePSF", 4, "using region %f,%f - %f,%f\n",
     82                  AnalysisRegion.x0, AnalysisRegion.y0,
     83                  AnalysisRegion.x1, AnalysisRegion.y1);
     84    }
     85
    2986    for (int i = 0; i < sources->n; i++) {
    3087        pmSource *inSource = sources->data[i];
    … …  
    3592        if (inSource->type == PM_SOURCE_DEFECT) continue;
    3693        if (inSource->type == PM_SOURCE_SATURATED) continue;
     94
     95        if (UseAnalysisRegion) {
     96            if (inSource->moments->x < AnalysisRegion.x0) continue;
     97            if (inSource->moments->y < AnalysisRegion.y0) continue;
     98            if (inSource->moments->x > AnalysisRegion.x1) continue;
     99            if (inSource->moments->y > AnalysisRegion.y1) continue;
     100        }
    37101
    38102        pmSource *otSource = pmSourceAlloc ();
    … …  
    51115            modelFLT->params->data.F32[2] = inSource->moments->x;
    52116            modelFLT->params->data.F32[3] = inSource->moments->y;
    53         }
    54         // XXX EAM : add option to peak-up on peak (for non-sat objects)
     117        } else {
     118            // peak-up on peak (for non-sat objects)
     119
     120            int ix = inSource->peak->x;
     121            int iy = inSource->peak->y;
     122
     123            psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
     124            psPlane min = psImageBicubeMin (bicube);
     125
     126            psTrace ("psphotEnsemblePSF", 5, "peak coord: %f %f -> %f %f\n",
     127                     modelFLT->params->data.F32[2], modelFLT->params->data.F32[3], min.x + ix, min.y + iy);
     128           
     129            // if min point is too deviant, keep the old value
     130            if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
     131                modelFLT->params->data.F32[2] = min.x + ix;
     132                modelFLT->params->data.F32[3] = min.y + iy;
     133            }
     134            psFree (bicube);
     135        }
    55136
    56137        // set PSF parameters for this model
    … …  
    149230        // subtract object
    150231        pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
     232        Fi->mode |= PM_SOURCE_SUBTRACTED;
     233        Fi->mode |= PM_SOURCE_TEMPSUB;
    151234    }
    152235
    153236    // XXX EAM : need to free up many things here
    154237
    155     psLogMsg ("psphot.emsemble", 4, "apply models: %f\n", psTimerMark ("psphot"));
     238    psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f\n", psTimerMark ("psphot"));
    156239    return true;
    157240}

  • trunk/psphot/src/psphotEvalFLT.c

    r5802 r5828  
    11# include "psphot.h"
    22
    3 bool psphotEvalFLT (pmSource *source)
     3bool psphotEvalFLT (pmSource *source, pmModel *model)
    44{
    55    int keep;
    6 
    7     pmModel *model = source->modelFLT;
    86
    97    // do we actually have a valid FLT model?

  • trunk/psphot/src/psphotEvalPSF.c

    r5802 r5828  
    4646// examine the model->status, fit parameters, etc and decide if the model succeeded
    4747// set the source->type and source->mode appropriately
    48 bool psphotEvalPSF (pmSource *source) {
     48bool psphotEvalPSF (pmSource *source, pmModel *model) {
    4949
    5050    int keep;
    5151    float dSX, dSY, SX, SY, SN;
    5252    float nSx, nSy, Chi;
    53 
    54     pmModel *model = source->modelPSF;
    5553
    5654    // do we actually have a valid PSF model?
    … …  
    127125    // this source is not a star, warn if it was a PSFSTAR
    128126    if (source->mode & PM_SOURCE_PSFSTAR) {
    129         psphotSaveImage (NULL, source->pixels, "failpx.fits");
    130         psphotSaveImage (NULL, source->mask, "failmk.fits");
    131         psphotSaveImage (NULL, source->weight, "failwt.fits");
    132127        psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality   (%f, %f  :  %f %f %f %f)\n",
     128                  model->params->data.F32[2], model->params->data.F32[3], nSx, nSy, SN, Chi);
     129    } else {
     130        psLogMsg ("psphot", 5, "fails PSF fit (%f, %f  :  %f %f %f %f)\n",
    133131                  model->params->data.F32[2], model->params->data.F32[3], nSx, nSy, SN, Chi);
    134132    }
    … …  
    150148    return false;
    151149}       
     150
     151// examine the model->status, fit parameters, etc and decide if the model succeeded
     152// set the source->type and source->mode appropriately
     153bool psphotEvalDBL (pmSource *source, pmModel *model) {
     154
     155    // do we actually have a valid PSF model?
     156    if (model == NULL) {
     157        source->mode &= ~PM_SOURCE_FITTED;
     158        return false;
     159    }
     160
     161    // did the model fit fail for one or another reason?
     162    switch (model->status) {
     163      case PM_MODEL_SUCCESS:
     164        break;
     165      case PM_MODEL_UNTRIED:
     166        source->mode &= ~PM_SOURCE_FITTED;
     167        return false;
     168      case PM_MODEL_BADARGS:
     169      case PM_MODEL_NONCONVERGE:
     170      case PM_MODEL_OFFIMAGE:
     171      default:
     172        source->mode |= PM_SOURCE_FAIL;
     173        return false;
     174    }
     175
     176    // unless we prove otherwise, this object is a star.
     177    source->type = PM_SOURCE_STAR;
     178
     179    // the following source->mode information pertains to modelPSF:
     180    source->mode |= PM_SOURCE_PSFMODEL;
     181
     182    // if the object has fitted peak above saturation, label as SATSTAR
     183    // this is a valid PSF object, but ignore the other quality tests
     184    // remember: fit does not use saturated pixels (masked)
     185    // XXX no extended object can saturate and stay extended...
     186    if (model->params->data.F32[1] >= SATURATION) {
     187        if (source->mode & PM_SOURCE_PSFSTAR) {
     188            psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n");
     189        }
     190        source->mode |=  PM_SOURCE_SATSTAR;
     191        return true;
     192    }
     193
     194    // if the object has a fitted peak below 0, the fit did not converge cleanly
     195    if (model->params->data.F32[1] <= 0) {
     196        source->mode |= PM_SOURCE_FAIL;
     197        return false;
     198    }
     199
     200    // if the source was predicted to be a SATSTAR, but it fitted below saturation,
     201    // make a note to the user
     202    if (source->mode & PM_SOURCE_SATSTAR) {
     203        psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n");
     204        source->mode &= ~PM_SOURCE_SATSTAR;
     205    }
     206    return true;
     207}       

  • trunk/psphot/src/psphotFitGalaxies.c

    r5772 r5828  
    4343
    4444        // sets the model radius (via source->model) and adjusts pixels as needed
    45         psphotCheckRadiusFLT (imdata, source);
     45        psphotCheckRadiusFLT (imdata, source, model);
    4646
    4747        // fit FLT (not PSF) model (set/unset the pixel mask)
    … …  
    5656        if (pmModelFitStatus (model)) {
    5757            pmSourceSubModel (source->pixels, source->mask, model, false, false);
    58             source->mode |= PM_SOURCE_SUBTRACTED;
     58            source->mode |=  PM_SOURCE_SUBTRACTED;
     59            source->mode &= ~PM_SOURCE_TEMPSUB;
    5960            Nsub ++;
    6061        }

  • trunk/psphot/src/psphotFixedPSF.c

    r5772 r5828  
    7777        // subtract object, leave local sky
    7878        pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
     79        source->mode |=  PM_SOURCE_SUBTRACTED;
     80        source->mode &= ~PM_SOURCE_TEMPSUB;
    7981    }
    8082

  • trunk/psphot/src/psphotFullFit.c

    r5802 r5828  
    11# include "psphot.h"
    2 
    3 # define CLEAR_TRACE \
    4   if (TEST_SRC_ON) { \
    5     psTraceSetLevel (".pmObjects.pmSourceFitModel", 0); \
    6     psTraceSetLevel ("psMinimizeLMChi2", 0); \
    7     TEST_SRC_ON = false; } \
    82
    93bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky)
    … …  
    3226    psphotInitRadiusFLT (config, sky, modelTypeFLT);
    3327
    34     float TEST_SRC_X = psMetadataLookupF32 (&status, config, "TEST_SRC_X");
    35     float TEST_SRC_Y = psMetadataLookupF32 (&status, config, "TEST_SRC_Y");
    36     bool TEST_SRC_ON = false;
    37     bool TEST_SRC    = status;
    38     fprintf (stderr, "test src: %f %f\n", TEST_SRC_X, TEST_SRC_Y);
    39 
    4028    for (int i = 0; i < sources->n; i++) {
    4129
    … …  
    5644        // replace object in image
    5745        pmSourceAddModel (source->pixels, source->mask, modelPSF, false, false);
     46        source->mode &= ~PM_SOURCE_SUBTRACTED;
    5847
    59         psphotCheckRadiusPSF (imdata, source);
     48        psphotCheckRadiusPSF (imdata, source, modelPSF);
    6049
    6150        x = modelPSF->params->data.F32[2];
    6251        y = modelPSF->params->data.F32[3];
    63 
    64         if (TEST_SRC && (fabs(TEST_SRC_X - x) < 5) && (fabs(TEST_SRC_Y - y) < 5)) {
    65             psTraceSetLevel (".pmObjects.pmSourceFitModel", 6);
    66             psTraceSetLevel ("psMinimizeLMChi2", 6);
    67             TEST_SRC_ON = true;
    68         }
    6952
    7053        // fit PSF model (set/unset the pixel mask)
    … …  
    7861
    7962        // does the PSF model succeed?
    80         if (psphotEvalPSF (source)) {
     63        if (psphotEvalPSF (source, source->modelPSF)) {
    8164            pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
    82             source->mode |= PM_SOURCE_SUBTRACTED;
     65            source->mode |=  PM_SOURCE_SUBTRACTED;
     66            source->mode &= ~PM_SOURCE_TEMPSUB;
    8367            psFree (LIN);
    8468            Nsub ++;
    85             CLEAR_TRACE;
    8669            continue;
    8770        }
     71
     72        // use the moments to get the extended source angle
     73        // displace on either side of the main peak by 1sigma or so
     74        // cut the peak in half?  use the main peak?
     75        // also try fit set on (unsaturated?) objects with a blend within 2-3 sigma
     76        // allow all possible nearby blends to be fit?
     77        // EllipseMoments moments;
     78        // moments.x2 = source->moments->Sx;
     79        // moments.y2 = source->moments->Sy;
     80        // moments.xy = source->moments->Sxy;
     81        // EllipseAxes axes = EllipseMomentsToAxes (moments);
    8882
    8983        // skip the source if we don't think it is extended
    … …  
    10195        pmModel *modelFLT = source->modelFLT;
    10296
    103         psphotCheckRadiusFLT (imdata, source);
     97        psphotCheckRadiusFLT (imdata, source, modelFLT);
    10498
    10599        x = modelFLT->params->data.F32[2];
    … …  
    116110
    117111        // does the FLT model succeed?
    118         if (psphotEvalFLT (source)) {
     112        if (psphotEvalFLT (source, source->modelFLT)) {
    119113            pmSourceSubModel (source->pixels, source->mask, source->modelFLT, false, false);
    120             source->mode |= PM_SOURCE_SUBTRACTED;
     114            source->mode |=  PM_SOURCE_SUBTRACTED;
     115            source->mode &= ~PM_SOURCE_TEMPSUB;
    121116            psFree (LIN);
    122117            Nsub ++;
    123             CLEAR_TRACE;
    124118            continue;
    125119        }
    126120
    127         // subtract PSF for object, leave local sky
     121        // re-subtract PSF for object, leave local sky
    128122    subLINEAR:
    129123        psFree (source->modelPSF);
    … …  
    131125        pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
    132126        source->mode |= PM_SOURCE_SUBTRACTED;
    133         source->mode |= PM_SOURCE_LINEAR;
     127        source->mode |= PM_SOURCE_TEMPSUB;
    134128        Nsub ++;
    135         CLEAR_TRACE;
    136129    }
    137130

  • trunk/psphot/src/psphotImageBackground.c

    r5772 r5828  
    3838    }
    3939    x->n = y->n = z->n = Nout;
    40     psLogMsg ("psphot", 5, "back: %f sec (select %d points)\n", psTimerMark ("psphot"), Nout);
    4140
    4241    psPolynomial2D *skyModel = psPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
    … …  
    4948
    5049    skyModel = psVectorFitPolynomial2D (skyModel, NULL, 0, z, NULL, x, y);
    51     psLogMsg ("psphot", 5, "back: %f sec (fit model)\n", psTimerMark ("psphot"));
    52     fprintf (stderr, "model: %f %f %f %f\n",
    53              skyModel->coeff[0][0], skyModel->coeff[1][0],
    54              skyModel->coeff[0][1], skyModel->coeff[1][1]);
    5550             
    5651    // this is a very inefficient way to evaluate the function..
    … …  
    6661    sky->sampleMean = 0.0;
    6762
    68     psLogMsg ("psphot", 5, "back: %f sec (fit model)\n", psTimerMark ("psphot"));
    69 
    70     psphotSaveImage (imdata->header, imdata->image, "backsub.fits");
     63    psLogMsg ("psphot", 3, "fit background model: %f sec\n", psTimerMark ("psphot"));
    7164    return (skyModel);
    7265}

  • trunk/psphot/src/psphotImageStats.c

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

  • trunk/psphot/src/psphotMagnitudes.c

    r5802 r5828  
    3737      rflux   = pow (10.0, 0.4*source->fitMag);
    3838      source->apMag  -= rflux * psf->skyBias * (M_PI * PS_SQR(apRadius));
    39       source->fitMag += psf->ApResid;
     39      source->fitMag += psPolynomial3DEval (psf->ApTrend, x, y, 0.0);
    4040    }
    4141
    … …  
    7070    return (model);
    7171}
     72
     73/*
     74aprMag' - fitMag = rflux*skyBias + ApTrend(x,y)
     75(aprMag - rflux*skyBias) - fitMAg = ApTrend(x,y)
     76(aprMag - rflux*skyBias) = fitMAg + ApTrend(x,y)
     77*/

  • trunk/psphot/src/psphotModelTest.c

    r5802 r5828  
    11# include "psphot.h"
     2# include "psEllipse.h"
    23
    3 int main (int argc, char **argv) {
     4bool psphotModelTest (eamReadout *imdata, psMetadata *config) {
    45
    5     psphotModelGroupInit ();
     6    bool status;
     7    int modelType;
     8    float obsMag, fitMag, value;
     9    char name[64];
     10    pmPSF *psf;
    611
    7     psMetadata *config = modelTestArguments (&argc, argv);
    8     eamReadout *imdata = psphotSetup (config);
     12    psMetadataItem *item  = NULL;
    913
    10     modelTestFitSource (imdata, config);
     14    // what fitting mode to use?
     15    char *psfModeWord = psMetadataLookupPtr (&status, config, "TEST_FIT_MODE");
     16    if (!status) {
     17        psfModeWord = psStringCopy ("FLT");
     18    }
     19    bool psfMode = !strcasecmp (psfModeWord, "PSF");
    1120
    12     exit (0);
     21    // in psfMode, psf sets the model type
     22    if (psfMode) {
     23        char *psfFile = psMetadataLookupPtr (&status, config, "PSF_INPUT_FILE");
     24        if (!status) psAbort ("psphotModelTest", "PSF_INPUT_FILE not supplied");
     25        psf = psphotReadPSF (psfFile);
     26        modelType = psf->type;
     27    } else {
     28        // find the model: supplied by user or first in the PSF_MODEL list
     29        char *modelName  = psMetadataLookupPtr (&status, config, "TEST_FIT_MODEL");
     30        if (modelName == NULL) {
     31            // get the list pointers for the PSF_MODEL entries
     32            psMetadataItem *mdi = psMetadataLookup (config, "PSF_MODEL");
     33            if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection");
     34
     35            // take the first list element
     36            psList *list = (psList *) mdi->data.list;
     37            item = psListGet (list, PS_LIST_HEAD);
     38            modelName = item->data.V;
     39        }
     40        modelType = pmModelSetType (modelName);
     41        if (modelType < 0) psAbort ("fitsource", "unknown model %s", modelName);
     42    }
     43
     44    // find the fitting parameters (try test values first)
     45    float INNER = psMetadataLookupF32 (&status, config, "TEST_FIT_INNER_RADIUS");
     46    if (!status) {
     47        INNER = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS");
     48    }
     49
     50    float OUTER = psMetadataLookupF32 (&status, config, "TEST_FIT_OUTER_RADIUS");
     51    if (!status) {
     52        OUTER = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
     53    }
     54
     55    float RADIUS = psMetadataLookupF32 (&status, config, "TEST_FIT_RADIUS");
     56    if (!status) {
     57        RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
     58    }
     59
     60    float mRADIUS = psMetadataLookupF32 (&status, config, "TEST_MOMENTS_RADIUS");
     61    if (!status) {
     62        mRADIUS = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS");
     63    }
     64
     65    // define the source of interest
     66    float xObj     = psMetadataLookupF32 (&status, config, "TEST_FIT_X");
     67    float yObj     = psMetadataLookupF32 (&status, config, "TEST_FIT_Y");
     68
     69    // construct the source structures
     70    pmSource *source = pmSourceAlloc();
     71    source->peak = pmPeakAlloc (xObj, yObj, 0, 0);
     72    psphotDefinePixels (source, imdata, xObj, yObj, OUTER);
     73
     74    // find the local sky
     75    status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER);
     76    if (!status) psAbort ("fitsource", "pmSourceLocalSky error");
     77
     78    // get the source moments
     79    status = pmSourceMoments (source, mRADIUS);
     80    if (!status) psAbort ("fitsource", "pmSourceLocalSky error");
     81    source->peak->counts = source->moments->Peak;
     82
     83    fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->x, source->moments->y);
     84
     85    EllipseMoments moments;
     86    moments.x2 = source->moments->Sx;
     87    moments.y2 = source->moments->Sy;
     88    moments.xy = source->moments->Sxy;
     89    EllipseAxes axes = EllipseMomentsToAxes (moments);
     90
     91    fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
     92
     93    // get the initial model parameter guess
     94    pmModel *model = pmSourceModelGuess (source, modelType);
     95    // if any parameters are defined, use those values
     96    int nParams = pmModelParameterCount (modelType);
     97    psF32 *params = model->params->data.F32;
     98    for (int i = 0; i < nParams; i++) {
     99        if (i == 2) {
     100            params[i] = xObj;
     101            continue;
     102        }
     103        if (i == 3) {
     104            params[i] = yObj;
     105            continue;
     106        }
     107        sprintf (name, "TEST_FIT_PAR%d", i);
     108        value = psMetadataLookupF32 (&status, config, name);
     109        if (status) {
     110            params[i] = value;
     111        }
     112    }
     113
     114    float area = params[4]*params[5];
     115    fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
     116
     117    if (psfMode) {
     118        pmModel *modelPSF = pmModelFromPSF (model, psf);
     119        psFree (model);
     120        model = modelPSF;
     121        params = model->params->data.F32;
     122    }
     123
     124    // list model input shape
     125    EllipseShape shape;
     126    shape.sx  = 1.4 / model->params->data.F32[4];
     127    shape.sy  = 1.4 / model->params->data.F32[5];
     128    shape.sxy = model->params->data.F32[6];
     129    axes = EllipseShapeToAxes (shape);
     130
     131    fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
     132
     133
     134    fprintf (stderr, "input parameters: \n");
     135    for (int i = 0; i < nParams; i++) {
     136        fprintf (stderr, "%d : %f\n", i, params[i]);
     137    }
     138
     139    // define the pixels used for the fit
     140    psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PSPHOT_MASK_MARKED);
     141
     142    char *fitset = psMetadataLookupPtr (&status, config, "TEST_FIT_SET");
     143    if (status) {
     144        status = psphotFitSet (source, model, fitset, psfMode);
     145        return status;
     146    }
     147
     148    status = pmSourceFitModel (source, model, psfMode);
     149
     150    // measure the source mags
     151    pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask);
     152    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag);
     153
     154    // write out positive object
     155    psphotSaveImage (NULL, source->pixels, "object.fits");
     156
     157    // subtract object, leave local sky
     158    pmSourceSubModel (source->pixels, source->mask, model, false, false);
     159   
     160    fprintf (stderr, "output parameters: \n");
     161    for (int i = 0; i < nParams; i++) {
     162        fprintf (stderr, "%d : %f\n", i, params[i]);
     163    }
     164
     165    // write out
     166    psphotSaveImage (NULL, source->pixels, "resid.fits");
     167    psphotSaveImage (NULL, source->mask, "mask.fits");
     168    return true;
    13169}

  • trunk/psphot/src/psphotOutput.c

    r5802 r5828  
    55
    66    bool status;
     7
     8    psTimerStart ("psphot");
    79
    810    char *outputMode = psMetadataLookupPtr (&status, config, "OUTPUT_MODE");
    911    char *outputFile = psMetadataLookupPtr (&status, config, "OUTPUT_FILE");
    1012    char *residImage = psMetadataLookupPtr (&status, config, "RESID_IMAGE");
     13    char *psfFile    = psMetadataLookupPtr (&status, config, "PSF_OUTPUT_FILE");
    1114
    1215    if (residImage != NULL) psphotSaveImage (imdata->header, imdata->image, residImage);
     16    if (psfFile != NULL) psphotWritePSF (psf, psfFile);
    1317
    1418    if (outputFile == NULL) {
    … …  
    4448        return;
    4549    }
     50
    4651    psAbort ("psphot", "unknown output mode %s", outputMode);
    4752}
    … …  
    287292
    288293    // set NAXIS to 0
    289     fprintf (stderr, "setting naxis\n");
    290294    mdi = psMetadataLookup (imdata->header, "NAXIS");
    291295    mdi->data.S32 = 0;
    … …  
    540544    }
    541545
    542     fprintf (stderr, "writing out moments\n");
    543546    for (i = 0; i < sources->n; i++) {
    544547        source = sources->data[i];
    … …  
    604607
    605608
    606 bool psphotSamplePSFs (pmPSF *psf, psImage *image) {
    607 
    608   // make sample PSFs for 4 corners and the center
     609bool psphotSamplePSFs (psMetadata *config, pmPSF *psf, psImage *image) {
     610
     611    bool status;
     612
     613    // make sample PSFs for 4 corners and the center
    609614    psImage *sample;
     615
     616    // optional dump of all rough source data
     617    char *output = psMetadataLookupPtr (&status, config, "PSF_SAMPLE_FILE");
     618    if (!status) return false;
     619    if (output == NULL) return false;
     620    if (output[0] == 0) return false;
    610621
    611622    pmModel *modelFLT = pmModelAlloc (psf->type);
    … …  
    613624    modelFLT->params->data.F32[1] = 1;
    614625
    615     psFits *fits = psFitsOpen ("sample.psf.fits", "w");
     626    psFits *fits = psFitsOpen (output, "w");
    616627
    617628    sample = pmModelPSFatXY (image, modelFLT, psf, 25, 25, 25, 25);
    … …  
    630641}
    631642
     643psPolynomial2D *psPolynomial2DfromMD (psMetadata *folder) {
     644
     645    bool status;
     646    char keyword[80];
     647
     648    // get polynomial orders
     649    // XXX add status failures tests
     650    int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
     651    int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
     652
     653    psPolynomial2D *poly = psPolynomial2DAlloc (nXorder, nYorder, PS_POLYNOMIAL_ORD);
     654
     655    for (int nx = 0; nx < poly->nX + 1; nx++) {
     656        for (int ny = 0; ny < poly->nY + 1; ny++) {
     657            sprintf (keyword, "VAL_X%02d_Y%02d", nx, ny);
     658            poly->coeff[nx][ny] = psMetadataLookupF32 (&status, folder, keyword);
     659            if (!status) poly->mask[nx][ny] = 1;
     660
     661            sprintf (keyword, "ERR_X%02d_Y%02d", nx, ny);
     662            poly->coeffErr[nx][ny] = psMetadataLookupF32 (&status, folder, keyword);
     663        }
     664    }
     665    return (poly);
     666}
     667
     668psPolynomial3D *psPolynomial3DfromMD (psMetadata *folder) {
     669
     670    bool status;
     671    char keyword[80];
     672
     673    // get polynomial orders
     674    // XXX add status failures tests
     675    int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
     676    int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
     677    int nZorder = psMetadataLookupS32 (&status, folder, "NORDER_Z");
     678
     679    psPolynomial3D *poly = psPolynomial3DAlloc (nXorder, nYorder, nZorder, PS_POLYNOMIAL_ORD);
     680
     681    for (int nx = 0; nx < poly->nX + 1; nx++) {
     682        for (int ny = 0; ny < poly->nY + 1; ny++) {
     683            for (int nz = 0; nz < poly->nZ + 1; nz++) {
     684                sprintf (keyword, "VAL_X%02d_Y%02d_Z%02d", nx, ny, nz);
     685                poly->coeff[nx][ny][nz] = psMetadataLookupF32 (&status, folder, keyword);
     686                if (!status) poly->mask[nx][ny][nz] = 1;
     687
     688                sprintf (keyword, "ERR_X%02d_Y%02d_Z%02d", nx, ny, nz);
     689                poly->coeffErr[nx][ny][nz] = psMetadataLookupF32 (&status, folder, keyword);
     690            }
     691        }
     692    }
     693    return (poly);
     694}
     695
     696// XXX : these may need F64, or %g format for output
     697bool psPolynomial2DtoMD (psMetadata *md, psPolynomial2D *poly, char *format, ...) {
     698
     699    int Nbyte;
     700    char tmp;
     701    char *root;
     702    va_list argp; 
     703
     704    va_start (argp, format);
     705    Nbyte = vsnprintf (&tmp, 0, format, argp);
     706    va_end (argp);
     707
     708    if (!Nbyte) return false;
     709
     710    va_start (argp, format);
     711    root = (char *) psAlloc (Nbyte + 1);
     712    memset (root, 0, Nbyte + 1);
     713    vsnprintf (root, Nbyte + 1, format, argp);
     714    va_end (argp);
     715
     716    psMetadata *folder = psMetadataAlloc ();
     717    psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 2D polynomial", folder);
     718    psFree (root);
     719
     720    // specify the polynomial orders
     721    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
     722    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
     723
     724    // place polynomial entries on folder
     725    for (int nx = 0; nx < poly->nX + 1; nx++) {
     726        for (int ny = 0; ny < poly->nY + 1; ny++) {
     727            if (poly->mask[nx][ny]) continue;
     728            psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny], nx, ny);
     729            psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d", PS_DATA_F32, "polynomial coefficient error", poly->coeffErr[nx][ny], nx, ny);
     730        }
     731    }
     732    return true;
     733}
     734
     735bool psPolynomial3DtoMD (psMetadata *md, psPolynomial3D *poly, char *format, ...) {
     736
     737    int Nbyte;
     738    char tmp;
     739    char *root;
     740    va_list argp; 
     741
     742    va_start (argp, format);
     743    Nbyte = vsnprintf (&tmp, 0, format, argp);
     744    va_end (argp);
     745
     746    if (!Nbyte) return false;
     747
     748    va_start (argp, format);
     749    root = (char *) psAlloc (Nbyte + 1);
     750    memset (root, 0, Nbyte + 1);
     751    vsnprintf (root, Nbyte + 1, format, argp);
     752    va_end (argp);
     753
     754    psMetadata *folder = psMetadataAlloc ();
     755    psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 3D polynomial", folder);
     756    psFree (root);
     757
     758    // specify the polynomial orders
     759    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
     760    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
     761    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Z", PS_DATA_S32, "number of z orders", poly->nZ);
     762
     763    // place polynomial entries on folder
     764    for (int nx = 0; nx < poly->nX + 1; nx++) {
     765        for (int ny = 0; ny < poly->nY + 1; ny++) {
     766            for (int nz = 0; nz < poly->nZ + 1; nz++) {
     767                if (poly->mask[nx][ny][nz]) continue;
     768                psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d_Z%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny][nz], nx, ny, nz);
     769                psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d_Z%02d", PS_DATA_F32, "polynomial coeffficient error", poly->coeffErr[nx][ny][nz], nx, ny, nz);
     770            }
     771        }
     772    }
     773    return true;
     774}
     775
     776bool psphotWritePSF (pmPSF *psf, char *filename) {
     777
     778    psMetadata *psfdata = psMetadataAlloc ();
     779
     780    char *modelName = pmModelGetType (psf->type);
     781    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName);
     782
     783    int nPar = pmModelParameterCount (psf->type)    ;
     784    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
     785
     786    for (int i = 0; i < nPar - 4; i++) {
     787        psPolynomial2D *poly = psf->params->data[i];
     788        psPolynomial2DtoMD (psfdata, poly, "PSF_PAR%02d", i);
     789    }
     790    psPolynomial3DtoMD (psfdata, psf->ApTrend, "APTREND");
     791
     792    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid);
     793    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid);
     794    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
     795
     796    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq);
     797    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars);
     798   
     799    psMetadataConfigWrite (psfdata, filename);
     800    psFree (psfdata);
     801
     802    return true;
     803}
     804
     805pmPSF *psphotReadPSF (char *filename) {
     806
     807    bool status;
     808    unsigned int Nfail;
     809    char keyword[80];
     810
     811    psMetadata *psfdata = psMetadataConfigParse (NULL, &Nfail, filename, FALSE);
     812
     813    char *modelName = psMetadataLookupPtr (&status, psfdata, "PSF_MODEL_NAME");
     814    pmModelType type = pmModelSetType (modelName);
     815
     816    pmPSF *psf = pmPSFAlloc (type);
     817
     818    int nPar = psMetadataLookupS32 (&status, psfdata, "PSF_MODEL_NPAR");
     819    if (nPar != pmModelParameterCount (psf->type)) psAbort ("read PSF" , "mismatch model par count");
     820
     821    for (int i = 0; i < nPar - 4; i++) {
     822        sprintf (keyword, "PSF_PAR%02d", i);
     823        psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword);
     824        psPolynomial2D *poly = psPolynomial2DfromMD (folder);
     825        psFree (psf->params->data[i]);
     826        psf->params->data[i] = poly;
     827    }
     828    sprintf (keyword, "APTREND");
     829    psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword);
     830    psPolynomial3D *poly = psPolynomial3DfromMD (folder);
     831    psFree (psf->ApTrend);
     832    psf->ApTrend = poly;
     833
     834    psf->ApResid = psMetadataLookupF32 (&status, psfdata, "PSF_AP_RESID");
     835    psf->dApResid = psMetadataLookupF32 (&status, psfdata, "PSF_dAP_RESID");
     836    psf->skyBias = psMetadataLookupF32 (&status, psfdata, "PSF_SKY_BIAS");
     837
     838    psf->chisq = psMetadataLookupF32 (&status, psfdata, "PSF_CHISQ");
     839    psf->nPSFstars = psMetadataLookupS32 (&status, psfdata, "PSF_NSTARS");
     840
     841    psFree (psfdata);
     842    return (psf);
     843}
     844
    632845bool psphotDumpMoments (psMetadata *config, psArray *sources) {
    633846

  • trunk/psphot/src/psphotRadiusChecks.c

    r5802 r5828  
    2020}
    2121
    22 bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source) {
    23 
    24     pmModel *model = source->modelPSF;
     22bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source, pmModel *model) {
    2523
    2624    // set the fit radius based on the object flux limit and the model
    … …  
    6462}
    6563
    66 bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source) {
    67 
    68     pmModel *model = source->modelFLT;
     64bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source, pmModel *model) {
    6965
    7066    // set the fit radius based on the object flux limit and the model

  • trunk/psphot/src/psphotReplaceUnfit.c

    r5802 r5828  
    33bool psphotReplaceUnfit (psArray *sources) {
    44
    5     pmModel *model;
    65    pmSource *source;
    76
    … …  
    1110      source = sources->data[i];
    1211
    13         // skip non-astronomical objects (very likely defects)
    14         // these were never fitted and subtracted
    15         if (source->mode &  PM_SOURCE_BLEND) continue;
    16         if (source->type == PM_SOURCE_SATURATED) continue;
     12        if (!(source->mode & PM_SOURCE_TEMPSUB)) continue;
     13        if (source->modelPSF == NULL) continue;
    1714
    18         if (source->type == PM_SOURCE_DEFECT) {
    19             if (source->mode &  PM_SOURCE_SUBTRACTED) goto addPSF;
    20             continue;
    21         }
    22 
    23         if (source->type == PM_SOURCE_GALAXY) {
    24             if (source->mode & PM_SOURCE_FAIL) goto addPSF;
    25             if (source->mode & PM_SOURCE_POOR) goto addPSF;
    26             continue;
    27         }
    28 
    29         // need to skip the successful fits
    30         if (source->type == PM_SOURCE_STAR) {
    31             if (source->mode & PM_SOURCE_FAIL) goto addPSF;
    32             if (source->mode & PM_SOURCE_POOR) goto addPSF;
    33             continue;
    34         }
    35 
    36     addPSF:
    37         model = source->modelPSF;
    38         if (model == NULL) continue;
    39 
    40         pmSourceAddModel (source->pixels, source->mask, model, false, false);
     15        pmSourceAddModel (source->pixels, source->mask, source->modelPSF, false, false);
    4116        source->mode &= ~PM_SOURCE_SUBTRACTED;
     17        source->mode &= ~PM_SOURCE_TEMPSUB;
    4218    }
    43     psLogMsg ("psphot.replace", 4, "replace models: %f (%d objects)\n", psTimerMark ("psphot"), sources->n);
     19    psLogMsg ("psphot.replace", 3, "replace unfitted models: %f sec (%d objects)\n", psTimerMark ("psphot"), sources->n);
    4420    return true;
    4521}
     22
     23# if (0)
     24if (source->mode &  PM_SOURCE_BLEND) continue;
     25if (source->type == PM_SOURCE_SATURATED) continue;
     26
     27if (source->type == PM_SOURCE_DEFECT) {
     28    if (source->mode &  PM_SOURCE_SUBTRACTED) goto addPSF;
     29    continue;
     30}
     31
     32if (source->type == PM_SOURCE_GALAXY) {
     33    if (source->mode & PM_SOURCE_FAIL) goto addPSF;
     34    if (source->mode & PM_SOURCE_POOR) goto addPSF;
     35    continue;
     36}
     37
     38// need to skip the successful fits
     39if (source->type == PM_SOURCE_STAR) {
     40    if (source->mode & PM_SOURCE_FAIL) goto addPSF;
     41    if (source->mode & PM_SOURCE_POOR) goto addPSF;
     42    continue;
     43}
     44# endif

  • trunk/psphot/src/psphotSetup.c

    r5802 r5828  
    102102    }
    103103
    104     psLogMsg ("psphot", 4, "load data: %f sec\n", psTimerMark ("psphot"));
     104    psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot"));
    105105
    106     psphotSaveImage (NULL, weight, "weight.fits");
    107     psphotSaveImage (NULL, mask, "mask.fits");
    108 
    109     // save the image data & return it
     106    // return image data
    110107    eamReadout *imdata = eamReadoutAlloc(image, weight, mask, header);
    111108    return (imdata);
Note: See TracChangeset for help on using the changeset viewer.