IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 17396 for trunk/psphot/src


Ignore:
Timestamp:
Apr 8, 2008, 8:36:06 AM (18 years ago)
Author:
eugene
Message:

merging from eam_branch_20080324 : psphot work on extended source fitting and related I/O functions

Location:
trunk/psphot/src
Files:
2 added
19 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/Makefile.am

    r16820 r17396  
    1818
    1919libpsphot_la_SOURCES = \
    20         psphotErrorCodes.c       \
    21         pmFootprint.c            \
    22         psphotVersion.c          \
    23         psphotModelGroupInit.c   \
    24         psphotMaskReadout.c      \
    25         psphotDefineFiles.c      \
    26         psphotReadout.c          \
    27         psphotModelBackground.c  \
    28         psphotSubtractBackground.c \
    29         psphotFindDetections.c   \
    30         psphotFindPeaks.c        \
    31         psphotSourceStats.c      \
    32         psphotRoughClass.c       \
    33         psphotBasicDeblend.c     \
    34         psphotChoosePSF.c        \
    35         psphotGuessModels.c      \
    36         psphotFitSourcesLinear.c \
    37         psphotBlendFit.c         \
    38         psphotReplaceUnfit.c     \
    39         psphotApResid.c          \
    40         psphotMagnitudes.c       \
    41         psphotSkyReplace.c       \
    42         psphotEvalPSF.c          \
    43         psphotEvalFLT.c          \
    44         psphotSourceFits.c       \
    45         psphotRadiusChecks.c     \
    46         psphotOutput.c           \
    47         psphotGrowthCurve.c      \
    48         psphotFakeSources.c      \
    49         psphotModelWithPSF.c     \
    50         psphotExtendedSources.c  \
    51         psphotRadialProfile.c    \
    52         psphotPetrosian.c        \
    53         psphotIsophotal.c        \
    54         psphotAnnuli.c           \
    55         psphotKron.c             \
    56         psphotKernelFromPSF.c    \
    57         psphotPSFConvModel.c     \
    58         psphotModelTest.c        \
    59         psphotFitSet.c           \
    60         psphotSourceFreePixels.c \
    61         psphotSummaryPlots.c     \
    62         psphotMergeSources.c     \
    63         psphotLoadPSF.c          \
    64         psphotReadoutCleanup.c   \
    65         psphotSourcePlots.c      \
    66         psphotRadialPlot.c       \
    67         psphotDeblendSatstars.c  \
    68         psphotMosaicSubimage.c   \
    69         psphotMakeResiduals.c    \
    70         psphotSourceSize.c       \
    71         psphotDiagnosticPlots.c  \
    72         psphotMakeFluxScale.c    \
    73         psphotCheckStarDistribution.c \
     20        psphotErrorCodes.c             \
     21        pmFootprint.c                  \
     22        psphotVersion.c                \
     23        psphotModelGroupInit.c         \
     24        psphotMaskReadout.c            \
     25        psphotDefineFiles.c            \
     26        psphotReadout.c                \
     27        psphotModelBackground.c        \
     28        psphotSubtractBackground.c     \
     29        psphotFindDetections.c         \
     30        psphotFindPeaks.c              \
     31        psphotSourceStats.c            \
     32        psphotRoughClass.c             \
     33        psphotBasicDeblend.c           \
     34        psphotChoosePSF.c              \
     35        psphotGuessModels.c            \
     36        psphotFitSourcesLinear.c       \
     37        psphotBlendFit.c               \
     38        psphotReplaceUnfit.c           \
     39        psphotApResid.c                \
     40        psphotMagnitudes.c             \
     41        psphotSkyReplace.c             \
     42        psphotEvalPSF.c                \
     43        psphotEvalFLT.c                \
     44        psphotSourceFits.c             \
     45        psphotRadiusChecks.c           \
     46        psphotOutput.c                 \
     47        psphotGrowthCurve.c            \
     48        psphotFakeSources.c            \
     49        psphotModelWithPSF.c           \
     50        psphotExtendedSourceAnalysis.c \
     51        psphotExtendedSourceFits.c     \
     52        psphotRadialProfile.c          \
     53        psphotPetrosian.c              \
     54        psphotIsophotal.c              \
     55        psphotAnnuli.c                 \
     56        psphotKron.c                   \
     57        psphotKernelFromPSF.c          \
     58        psphotPSFConvModel.c           \
     59        psphotModelTest.c              \
     60        psphotFitSet.c                 \
     61        psphotSourceFreePixels.c       \
     62        psphotSummaryPlots.c           \
     63        psphotMergeSources.c           \
     64        psphotLoadPSF.c                \
     65        psphotReadoutCleanup.c         \
     66        psphotSourcePlots.c            \
     67        psphotRadialPlot.c             \
     68        psphotDeblendSatstars.c        \
     69        psphotMosaicSubimage.c         \
     70        psphotMakeResiduals.c          \
     71        psphotSourceSize.c             \
     72        psphotDiagnosticPlots.c        \
     73        psphotMakeFluxScale.c          \
     74        psphotCheckStarDistribution.c  \
    7475        psphotAddNoise.c
    7576

  • trunk/psphot/src/psphot.h

    r16820 r17396  
    4545bool            psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background);
    4646bool            psphotSkyReplace (pmConfig *config, const pmFPAview *view);
    47 bool            psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe);
     47bool            psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
     48bool            psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe);
    4849
    4950// used by psphotFindDetections
    … …  
    8788bool            psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, psMaskType maskVal);
    8889bool            psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psMaskType maskVal);
    89 pmModel        *psphotFitEXT (pmReadout *readout, pmSource *source, psMaskType maskVal);
     90pmModel        *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal);
    9091psArray        *psphotFitDBL (pmReadout *readout, pmSource *source, psMaskType maskVal);
    9192
    … …  
    117118bool            psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal);
    118119
    119 bool            psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal);
     120pmModel        *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal, int psfSize);
     121
    120122psKernel       *psphotKernelFromPSF (pmSource *source, int nPix);
    121123

  • trunk/psphot/src/psphotAnnuli.c

    r15562 r17396  
    1 # include "psphot.h"
     1# include "psphotInternal.h"
    22
    33bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
    … …  
    5555  source->extpars->annuli->fluxVar = fluxSquare;
    5656
     57  psFree (pixelCount);
     58
    5759  return true;
    5860}

  • trunk/psphot/src/psphotBlendFit.c

    r16820 r17396  
    4040        // skip non-astronomical objects (very likely defects)
    4141        if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
     42        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) continue;
    4243        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    4344        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    … …  
    7980        // XXX re-consider conditions under which the source has EXT fit:
    8081        // I should try EXT if the source size measurement says it is large
    81         if (psphotFitBlend (readout, source, psf, maskVal)) {
    82             psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y);
    83             Npsf ++;
    84             continue;
    85         }
    86         if (psphotFitBlob (readout, source, sources, psf, maskVal)) {
    87             psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y);
    88             Next ++;
    89             continue;
    90         }
     82        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     83            if (psphotFitBlob (readout, source, sources, psf, maskVal)) {
     84                source->type = PM_SOURCE_TYPE_EXTENDED;
     85                psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y);
     86                Next ++;
     87                continue;
     88            }
     89        } else {
     90            if (psphotFitBlend (readout, source, psf, maskVal)) {
     91                source->type = PM_SOURCE_TYPE_STAR;
     92                psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y);
     93                Npsf ++;
     94                continue;
     95            }
     96        }
    9197
    9298        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->moments->x, source->moments->y);

  • trunk/psphot/src/psphotEvalFLT.c

    r15060 r17396  
    88    if (model == NULL) {
    99        source->mode &= ~PM_SOURCE_MODE_FITTED;
     10        psTrace ("psphot", 5, "no model fitted?\n");
    1011        return false;
    1112    }
    … …  
    1415    if (!(model->flags & PM_MODEL_STATUS_FITTED)) {
    1516        source->mode &= ~PM_SOURCE_MODE_FITTED;
     17        psTrace ("psphot", 5, "no model fitted?\n");
    1618        return false;
    1719    }
     20
    1821    // did the model fit fail for one or another reason?
    1922    if (model->flags & (PM_MODEL_STATUS_BADARGS |
    … …  
    2225        source->mode |= PM_SOURCE_MODE_FAIL;
    2326        psLogMsg ("psphot", 5, "EXT fail fit\n");
     27        psTrace ("psphot", 5, "EXT fail fit\n");
     28
     29        if (model->flags & PM_MODEL_STATUS_OFFIMAGE) {
     30          psTrace ("psphot", 5, "off image\n");
     31        }
     32        if (model->flags & PM_MODEL_STATUS_BADARGS) {
     33          psTrace ("psphot", 5, "bad args\n");
     34        }
     35        if (model->flags & PM_MODEL_STATUS_NONCONVERGE) {
     36          psTrace ("psphot", 5, "non converge\n");
     37        }
    2438        return false;
    2539    }
    … …  
    3549    if (model->params->data.F32[PM_PAR_I0] <= 0.02) {
    3650        source->mode |= PM_SOURCE_MODE_FAIL;
     51        psTrace ("psphot", 5, "model central intensity ~ zero\n");
    3752        return false;
    3853    }
    … …  
    4358    // poor-quality fit; only keep if nothing else works...
    4459    psLogMsg ("psphot", 5, "EXT poor fit\n");
     60    psTrace ("psphot", 5, "EXT poor fit\n");
     61
    4562    source->mode |= PM_SOURCE_MODE_POOR;
    4663    return false;

  • trunk/psphot/src/psphotFindDetections.c

    r16820 r17396  
    88    bool status;
    99    int pass;
     10
     11    psTimerStart ("psphot");
    1012
    1113    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    6769    psphotCullPeaks(readout->image, readout->weight, recipe, detections->footprints);
    6870    detections->peaks = pmFootprintArrayToPeaks(detections->footprints);
    69     psLogMsg ("psphot", PS_LOG_INFO, "peaks: %ld, total footprints: %ld\n", detections->peaks->n, detections->footprints->n);
     71    psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot"));
    7072
    7173    return detections;

  • trunk/psphot/src/psphotFindPeaks.c

    r16820 r17396  
    1111
    1212    // smooth the image and weight map
    13     psTimerStart ("psphot");
     13    psTimerStart ("peaks");
    1414
    1515    // XXX if we have been supplied a PSF, we can use that to set the FWHM_X,FWHM_Y values
    … …  
    3232    psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32);
    3333    psImageSmoothMaskF32 (smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH);
    34     psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("psphot"));
     34    psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("peaks"));
    3535
    3636    // smooth the weight, applying the mask as we go
    3737    psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32);
    3838    psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH);
    39     psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("psphot"));
     39    psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("peaks"));
    4040
    4141    psImage *mask = readout->mask;
    … …  
    6161        }
    6262    }
    63     psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("psphot"));
     63    psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("peaks"));
    6464
    6565    // optionally save example images under trace
    … …  
    7070    }
    7171
    72     psTimerStart ("psphot");
     72    psTimerStart ("peaks");
    7373    // set peak threshold
    7474
    … …  
    140140        pmPeaksWriteText (peaks, output);
    141141    }
    142     psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot"));
     142    psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("peaks"));
    143143
    144144    // If they asked us to return a set of pmFootprints, not a raw pmPeak list,
    … …  
    167167        peaks = footprints;             // well, you know what I mean
    168168
    169         psLogMsg ("psphot", PS_LOG_INFO, "%ld footprints: %f sec\n", peaks->n, psTimerMark ("psphot"));
     169        psLogMsg ("psphot", PS_LOG_INFO, "%ld footprints: %f sec\n", peaks->n, psTimerMark ("peaks"));
    170170    }
    171171

  • trunk/psphot/src/psphotIsophotal.c

    r15562 r17396  
    1 # include "psphot.h"
     1# include "psphotInternal.h"
     2
     3static float ISOPHOT_FLUX = NAN;
    24
    35bool psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
    … …  
    1416
    1517  // flux at which to measure isophotal parameters
    16   // XXX cache this?
    17   float ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX");
     18  if (!isfinite(ISOPHOT_FLUX)) {
     19    // XXX ISOPHOTAL_FLUX should be specified in mags, need the zero point to get counts/sec
     20    ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX");
     21    assert (status);
     22  }
    1823
    1924  // find the first bin below the flux level and the last above the level
    20   // XXX can this be done faster with disection?
     25  // XXX can this be done faster with bisection?
    2126  // XXX do I need to worry about crazy outliers? 
    2227  // XXX should i be smoothing or fitting the curve?
    … …  
    2732    if ((firstBelow < 0) && (flux->data.F32[i] < ISOPHOT_FLUX)) firstBelow = i;
    2833  }
    29 
     34  // if we don't go out far enough, we have a problem...
     35  if (lastAbove == flux->n - 1) {
     36    psTrace ("psphot", 5, "did not go out far enough to reach isophotal magnitude");
     37    // XXX raise a flag ?
     38    return false;
     39  }
     40  if (firstBelow < 0) {
     41    psTrace ("psphot", 5, "did not go out far enough to bound isophotal magnitude: error unmeasured");
     42    // XXX raise a flag ?
     43    lastAbove = firstBelow;
     44    return false;
     45  }
     46 
    3047  // need to examine pixels in this vicinity
    3148  float isophotalFluxFirst = 0;
    … …  
    5673  source->extpars->isophot->radErr = isophotalRadErr;
    5774
     75  psTrace ("psphot", 5, "Isophot flux:%f +/- %f @ %f +/- %f for %f, %f\n",
     76           source->extpars->isophot->mag, source->extpars->isophot->magErr,
     77           source->extpars->isophot->rad, source->extpars->isophot->radErr,
     78           source->peak->xf, source->peak->yf);
     79
    5880  return true;
    5981

  • trunk/psphot/src/psphotKernelFromPSF.c

    r14338 r17396  
    1 # include "psphot.h"
     1# include "psphotInternal.h"
    22
    33psKernel *psphotKernelFromPSF (pmSource *source, int nPix) {
    44
    5   assert (source);
    6   assert (source->psfFlux); // XXX build if needed?
     5    assert (source);
     6    assert (source->psfFlux); // XXX build if needed?
    77
    8   int x0 = source->peak->xf - source->psfFlux->col0;
    9   int y0 = source->peak->yf - source->psfFlux->row0;
     8    int x0 = source->peak->xf - source->psfFlux->col0;
     9    int y0 = source->peak->yf - source->psfFlux->row0;
    1010
    11   // need to decide on the size: dynamically? statically?
    12   psKernel *psf = psKernelAlloc (-nPix, +nPix, -nPix, +nPix);
     11    // need to decide on the size: dynamically? statically?
     12    psKernel *psf = psKernelAlloc (-nPix, +nPix, -nPix, +nPix);
    1313
    14   for (int j = psf->yMin; j <= psf->yMax; j++) {
    15     for (int i = psf->xMin; i <= psf->xMax; i++) {
    16       psf->kernel[j][i] = source->psfFlux->data.F32[y0 + j][x0 + i];
     14    // XXX we should just re-construct a PSF at this location
     15    // psModelAdd (psf->image, NULL, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM | PM_MODEL_OP_CENTER);
     16 
     17    // if the realized PSF for this object does not cover the full kernel, give up for now
     18    if (x0 + psf->xMin < 0) goto escape;
     19    if (x0 + psf->xMax >= source->psfFlux->numCols) goto escape;
     20    if (y0 + psf->yMin < 0) goto escape;
     21    if (y0 + psf->yMax >= source->psfFlux->numRows) goto escape;
     22
     23    double sum = 0.0;
     24    for (int j = psf->yMin; j <= psf->yMax; j++) {
     25        for (int i = psf->xMin; i <= psf->xMax; i++) {
     26            double value = source->psfFlux->data.F32[y0 + j][x0 + i];
     27            psf->kernel[j][i] = value;
     28            sum += value;
     29        }
    1730    }
    18   }
    19   return psf;
     31    assert (sum > 0.0);
     32
     33    // psf must be normalized (integral = 1.0)
     34    for (int i = 0; i < psf->image->numRows; i++) {
     35        for (int j = 0; j < psf->image->numCols; j++) {
     36            psf->image->data.F32[i][j] /= sum;
     37        }
     38    }
     39
     40    return psf;
     41
     42escape:
     43    psFree (psf);
     44    return NULL;
    2045}

  • trunk/psphot/src/psphotKron.c

    r13983 r17396  
    1 # include "psphot.h"
     1# include "psphotInternal.h"
    22
    33bool psphotKron (pmSource *source, psMetadata *recipe, psMaskType maskVal) {

  • trunk/psphot/src/psphotModelTest.c

    r16820 r17396  
    208208        source->modelPSF = pmModelFromPSF (model, psf);
    209209        source->modelEXT = model;
    210         status = psphotPSFConvModel (source, recipe, maskVal);
    211         model = source->modelConv;
     210
     211        // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
     212        int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
     213        assert (status);
     214
     215        model = psphotPSFConvModel (readout, source, modelType, maskVal, psfSize);
    212216        params = model->params->data.F32;
    213217    } else {

  • trunk/psphot/src/psphotModelWithPSF.c

    r14347 r17396  
    1 # include "psphot.h"
     1# include "psphotInternal.h"
     2# define SAVE_IMAGES 0
    23
    34bool psphotModelWithPSF_LMM (
    … …  
    3637
    3738    // allocate internal arrays (current vs Guess)
    38     psImage *alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    39     psImage *Alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    40     psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F32);
    41     psVector *Beta   = psVectorAlloc(params->n, PS_TYPE_F32);
     39    psImage *Alpha = NULL;
     40    psVector *Beta = NULL;
     41
     42    // Alpha & Beta only contain elements to represent the unmasked parameters
     43    if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) {
     44        psAbort ("programming error: no unmasked parameters to be fit\n");
     45    }
     46   
     47    // allocate internal arrays (current vs Guess)
     48    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
     49    psVector *beta   = psVectorAlloc(Beta->n, PS_TYPE_F32);
    4250    psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
    4351
    … …  
    7987        // dump some useful info if trace is defined
    8088        if (psTraceGetLevel("psphot") >= 6) {
    81             p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (1)");
    82             p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (1)");
     89            p_psImagePrint(psTraceGetDestination(), Alpha, "Alpha guess (1)");
     90            p_psVectorPrint(psTraceGetDestination(), Beta, "Beta guess (1)");
     91            p_psVectorPrint(psTraceGetDestination(), beta, "beta current (1)");
    8392        }
    8493        if (psTraceGetLevel("psphot") >= 5) {
    … …  
    116125            beta   = psVectorCopy(beta, Beta, PS_TYPE_F32);
    117126            params = psVectorCopy(params, Params, PS_TYPE_F32);
    118             lambda *= 0.1;
     127            lambda *= 0.25;
    119128
    120129            // save the new convolved model image
    … …  
    130139    // construct & return the covariance matrix (if requested)
    131140    if (covar != NULL) {
    132         if (!psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
     141        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
    133142            psTrace ("psphot", 5, "failure to calculate covariance matrix\n");
    134143        }
     144        // set covar values which are not masked
     145        psImageInit (covar, 0.0);
     146        for (int j = 0, J = 0; j < params->n; j++) {
     147            if (paramMask && (paramMask->data.U8[j])) {
     148                covar->data.F32[j][j] = 1.0;
     149                continue;
     150            }
     151            for (int k = 0, K = 0; k < params->n; k++) {
     152                if (paramMask && (paramMask->data.U8[k])) continue;
     153                covar->data.F32[j][k] = Alpha->data.F32[J][K];
     154                K++;
     155            }
     156            J++;
     157        }
    135158    }
    136159
    … …  
    177200    }
    178201
    179     // generate the model and derivative images for this parameter set
     202    // 1 *** generate the model and derivative images for this parameter set
    180203
    181204    // storage for model derivatives
    … …  
    223246
    224247            for (int n = 0; n < params->n; n++) {
    225               if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
    226               psImage *dmodel = pcm->dmodels->data[n];
    227               dmodel->data.F32[i][j] = deriv->data.F32[n];
     248                if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
     249                psImage *dmodel = pcm->dmodels->data[n];
     250                dmodel->data.F32[i][j] = deriv->data.F32[n];
    228251            }
    229252        }
    … …  
    231254    psFree(coord);
    232255    psFree(deriv);
    233 
    234256
    235257    // convolve model and dmodel arrays with PSF
    236258    psImageConvolveDirect (pcm->modelConv, pcm->model, psf);
    237259    for (int n = 0; n < pcm->dmodels->n; n++) {
    238       if (pcm->dmodels->data[n] == NULL) continue;
    239       psImage *dmodel = pcm->dmodels->data[n];
    240       psImage *dmodelConv = pcm->dmodelsConv->data[n];
    241       psImageConvolveDirect (dmodelConv, dmodel, psf);
     260        if (pcm->dmodels->data[n] == NULL) continue;
     261        psImage *dmodel = pcm->dmodels->data[n];
     262        psImage *dmodelConv = pcm->dmodelsConv->data[n];
     263        psImageConvolveDirect (dmodelConv, dmodel, psf);
    242264    }
    243265
    244266    // XXX TEST : SAVE IMAGES
    245 # if (1)
     267# if (SAVE_IMAGES)
    246268    psphotSaveImage (NULL, psf->image, "psf.fits");
    247269    psphotSaveImage (NULL, pcm->model, "model.fits");
    … …  
    252274# endif
    253275
     276    // 2 *** accumulate alpha & beta
     277
    254278    // zero alpha and beta for summing below
    255279    psImageInit (alpha, 0.0);
    … …  
    259283    for (psS32 i = 0; i < source->pixels->numRows; i++) {
    260284        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    261           // XXX are we doing the right thing with the mask?
     285            // XXX are we doing the right thing with the mask?
    262286            // skip masked points
    263287            if (source->maskObj->data.U8[i][j]) {
    … …  
    279303            chisq += PS_SQR(delta) * yweight;
    280304
    281             if (isnan(delta))
    282               psAbort("nan in delta");
    283             if (isnan(chisq))
    284               psAbort("nan in chisq");
    285 
    286             for (psS32 n1 = 0; n1 < params->n; n1++) {
    287               if ((paramMask != NULL) && (paramMask->data.U8[n1])) {
    288                 continue;
    289               }
    290               psImage *dmodel = pcm->dmodelsConv->data[n1];
    291               float weight = dmodel->data.F32[i][j] * yweight;
    292               for (psS32 n2 = 0; n2 <= n1; n2++) {
    293                 if ((paramMask != NULL) && (paramMask->data.U8[n2])) {
    294                   continue;
    295                 }
    296                 dmodel = pcm->dmodelsConv->data[n2];
    297                 alpha->data.F32[n1][n2] += weight * dmodel->data.F32[i][j];
    298               }
    299               beta->data.F32[n1] += weight * delta;
     305            if (isnan(delta)) psAbort("nan in delta");
     306            if (isnan(chisq)) psAbort("nan in chisq");
     307
     308            // alpha & beta only contain unmasked elements
     309            for (int n1 = 0, N1 = 0; n1 < params->n; n1++) {
     310                if ((paramMask != NULL) && (paramMask->data.U8[n1])) continue;
     311                psImage *dmodel = pcm->dmodelsConv->data[n1];
     312                float weight = dmodel->data.F32[i][j] * yweight;
     313                for (int n2 = 0, N2 = 0; n2 <= n1; n2++) {
     314                    if ((paramMask != NULL) && (paramMask->data.U8[n2])) continue;
     315                    dmodel = pcm->dmodelsConv->data[n2];
     316                    alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j];
     317                    N2++;
     318                }
     319                beta->data.F32[N1] += weight * delta;
     320                N1++;
    300321            }
    301322        }
    … …  
    303324
    304325    // calculate lower-left half of alpha
    305     for (psS32 j = 1; j < params->n; j++) {
     326    for (psS32 j = 1; j < alpha->numCols; j++) {
    306327        for (psS32 k = 0; k < j; k++) {
    307328            alpha->data.F32[k][j] = alpha->data.F32[j][k];
    308         }
    309     }
    310 
    311     // fill in pivots if we apply a mask
    312     if (paramMask != NULL) {
    313         for (psS32 j = 0; j < params->n; j++) {
    314             if (paramMask->data.U8[j]) {
    315                 alpha->data.F32[j][j] = 1;
    316                 beta->data.F32[j] = 1;
    317             }
    318329        }
    319330    }
    … …  
    345356    pcm->dmodels = psArrayAlloc (params->n);
    346357    for (psS32 n = 0; n < params->n; n++) {
    347       if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
    348       pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     358        pcm->dmodels->data[n] = NULL;
     359        if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
     360        pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
    349361    }
    350362
    … …  
    353365    pcm->dmodelsConv = psArrayAlloc (params->n);
    354366    for (psS32 n = 0; n < params->n; n++) {
    355       if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
    356       pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     367        pcm->dmodelsConv->data[n] = NULL;
     368        if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
     369        pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
    357370    }
    358371
    … …  
    378391 
    379392 while () {
    380    GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda)
    381    dLinear = dLinear(Beta, beta, lambda);
    382    chisq = SetABX (alpha, beta, params, paramMask, x, y, dy, func)
    383    convergence tests...
     393 GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda)
     394 dLinear = dLinear(Beta, beta, lambda);
     395 chisq = SetABX (alpha, beta, params, paramMask, x, y, dy, func)
     396 convergence tests...
    384397 }
    385398 
    … …  
    388401 ** GuessABP:
    389402
    390       f_c = sum_i (kern_i * func (x_i; p_o))
    391 
    392 df_c/dp_o = d/dp_o [sum_i (kern_i * func (x_i; p_o))]
    393 
    394 df_c/dp_o = sum_i (d/dp_o [kern_i * func (x_i; p_o)])
    395 
    396 df_c/dp_o = sum_i (kern_i * d/dp_o [func (x_i; p_o)])
    397 
    398 - generate image arrays for func, dfunc/dp_j (not masked)
    399 - convolve each with psf
    400 - measure delta = f_conv - data
    401 - etc
     403 f_c = sum_i (kern_i * func (x_i; p_o))
     404
     405 df_c/dp_o = d/dp_o [sum_i (kern_i * func (x_i; p_o))]
     406
     407 df_c/dp_o = sum_i (d/dp_o [kern_i * func (x_i; p_o)])
     408
     409 df_c/dp_o = sum_i (kern_i * d/dp_o [func (x_i; p_o)])
     410
     411 - generate image arrays for func, dfunc/dp_j (not masked)
     412 - convolve each with psf
     413 - measure delta = f_conv - data
     414 - etc
    402415*/
    403416

  • trunk/psphot/src/psphotPSFConvModel.c

    r14655 r17396  
    1 # include "psphot.h"
     1# include "psphotInternal.h"
     2# define USE_DELTA_PSF 0
    23
    34// save as static values so they may be set externally
    45static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
    56static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
    6 // static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0;
    7 // static bool  PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;
    87
    98// input source has both modelPSF and modelEXT.  on successful exit, we set the
    109// modelConv to contain the fitted parameters, and the modelFlux to contain the
    1110// convolved model image.
    12 bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
     11pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal, int psfSize) {
    1312   
    14     bool status;
    15 
    16     int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
    17     if (!status) {
    18         psfSize = 2;
    19     }
    20 
    2113    // make sure we save a cached copy of the psf flux
    2214    pmSourceCachePSF (source, maskVal);
    2315
    2416    // convert the cached cached psf model for this source to a psKernel
    25     // XXX for the moment, hard-wire the kernel to be 5x5 (2 pix radius)
    26     // XXX for the moment, hard-wire the kernel to be 9x9 (4 pix radius)
    2717    psKernel *psf = psphotKernelFromPSF (source, psfSize);
     18    if (!psf) return NULL;
    2819
    29     // psf must be normalized (integral = 1.0)
    30     double sum = 0.0;
    31     for (int i = 0; i < psf->image->numRows; i++) {
    32         for (int j = 0; j < psf->image->numCols; j++) {
    33             sum += psf->image->data.F32[i][j];
    34         }
    35     }
    36     assert (sum > 0.0);
    37     for (int i = 0; i < psf->image->numRows; i++) {
    38         for (int j = 0; j < psf->image->numCols; j++) {
    39             psf->image->data.F32[i][j] /= sum;
    40         }
    41     }
    42 
    43 # if (0)
    44     // XXX sanity check: convolve with delta function should behave like unconvolved version
    45     for (int i = 0; i < psf->image->numRows; i++) {
    46         for (int j = 0; j < psf->image->numCols; j++) {
    47             psf->image->data.F32[i][j] = 0.0;
    48         }
    49     }
    50     psf->image->data.F32[2][2] = 1.0;
     20# if (USE_DELTA_PSF)
     21    psImageInit (psf->image, 0.0);
     22    psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0;
    5123# endif
    5224
    … …  
    5426    // XXX we could modify the parameter values or even the model
    5527    // here based on the observed seeing (some lookup table...)
    56     pmModel *modelConv = pmModelCopy (source->modelEXT);
     28
     29    // use the source moments, etc to guess basic model parameters
     30    pmModel *modelConv = pmSourceModelGuess (source, modelType);
     31    if (!modelConv) {
     32        psFree (psf);
     33        return NULL;
     34    }
     35
     36    // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
     37    psEllipseShape psfShape;
     38    psfShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
     39    psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
     40    psfShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
     41    psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0);
     42
     43    // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
     44    psEllipseShape extShape;
     45    extShape.sx  = modelConv->params->data.F32[PM_PAR_SXX] / M_SQRT2;
     46    extShape.sxy = modelConv->params->data.F32[PM_PAR_SXY];
     47    extShape.sy  = modelConv->params->data.F32[PM_PAR_SYY] / M_SQRT2;
     48    psEllipseAxes extAxes = psEllipseShapeToAxes (extShape, 20.0);
     49
     50    // decrease the initial guess ellipse by psf_minor axis:
     51    psEllipseAxes extAxesMod;
     52    extAxesMod.major = sqrt (PS_MAX (1.0, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor)));
     53    extAxesMod.minor = sqrt (PS_MAX (1.0, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor)));
     54    extAxesMod.theta = extAxes.theta;
     55
     56    psEllipseShape extShapeMod = psEllipseAxesToShape (extAxesMod);
     57    modelConv->params->data.F32[PM_PAR_SXX] = extShapeMod.sx * M_SQRT2;
     58    modelConv->params->data.F32[PM_PAR_SXY] = extShapeMod.sxy;
     59    modelConv->params->data.F32[PM_PAR_SYY] = extShapeMod.sy * M_SQRT2;
     60
     61    // increase the initial guess central intensity by 2pi r^2:
     62    modelConv->params->data.F32[PM_PAR_I0] *= (1.0 + PS_SQR(psfAxes.minor) / PS_SQR(extAxesMod.minor));
     63
    5764    psVector *params  = modelConv->params;
    5865    psVector *dparams = modelConv->dparams;
     66
     67    psphotCheckRadiusEXT (readout, source, modelConv);
    5968
    6069    // create the minimization constraints
    … …  
    91100    psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
    92101
    93     // renormalize output model image
     102    // renormalize output model image (generated by fitting process)
    94103    float Io = params->data.F32[PM_PAR_I0];
    95104    for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
    … …  
    115124    onPic &= (params->data.F32[PM_PAR_YPOS] >= source->pixels->row0);
    116125    onPic &= (params->data.F32[PM_PAR_YPOS] <  source->pixels->row0 + source->pixels->numRows);
    117     if (!onPic) {
    118         modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE;
    119     }
     126    if (!onPic) modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE;
    120127
    121     source->mode |= PM_SOURCE_MODE_FITTED;
    122     source->modelConv = modelConv;
     128    source->mode |= PM_SOURCE_MODE_FITTED; // XXX is this needed?
    123129
     130    psFree(psf);
    124131    psFree(myMin);
    125132    psFree(covar);
    126133    psFree(constraint);
    127134
    128     bool retval = (onPic && fitStatus);
    129     psTrace("psphot", 5, "---- %s(%d) end ----\n", __func__, retval);
    130     return(retval);
     135    return modelConv;
    131136}

  • trunk/psphot/src/psphotPetrosian.c

    r15800 r17396  
    1 # include "psphot.h"
     1# include "psphotInternal.h"
     2
     3static float PETROSIAN_R0 = NAN;
     4static float PETROSIAN_RF = NAN;
    25
    36bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
    … …  
    1417
    1518  // flux at which to measure isophotal parameters
    16   // XXX cache this?
    17   float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0");
    18   float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO");
     19  if (!isfinite(PETROSIAN_R0)) {
     20    PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0");
     21    PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO");
     22    assert (status);
     23  }
    1924
    2025  // first find flux at R0
    21   int firstBelow = -1;
    22   int lastAbove = -1;
     26  int firstAbove = -1;
     27  int lastBelow = -1;
    2328  for (int i = 0; i < radius->n; i++) {
    24     if (radius->data.F32[i] > PETROSIAN_R0) lastAbove = i;
    25     if ((firstBelow < 0) && (flux->data.F32[i] < PETROSIAN_R0)) firstBelow = i;
     29    if (radius->data.F32[i] < PETROSIAN_R0) lastBelow = i;
     30    if ((firstAbove < 0) && (radius->data.F32[i] > PETROSIAN_R0)) firstAbove = i;
     31  }
     32  // if we don't go out far enough, we have a problem...
     33  if (lastBelow == radius->n - 1) {
     34    psTrace ("psphot", 5, "did not go out far enough to reach petrosian reference radius...");
     35    // XXX skip object? raise a flag ?
     36    return false;
     37  }
     38  if (firstAbove < 0) {
     39    psTrace ("psphot", 5, "did not go out far enough to bound petrosian reference radius");
     40    // XXX raise a flag ?
     41    return false;
    2642  }
    2743
    … …  
    2945  float fluxR0 = 0.0;
    3046  int fluxRn = 0;
    31   for (int i = PS_MIN(firstBelow, lastAbove); i <= PS_MAX(firstBelow, lastAbove); i++) {
     47  for (int i = PS_MIN(firstAbove, lastBelow); i <= PS_MAX(firstAbove, lastBelow); i++) {
    3248    fluxR0 += flux->data.F32[i];
    3349    fluxRn ++;
    … …  
    4258  // XXX do I need to worry about crazy outliers? 
    4359  // XXX should i be smoothing or fitting the curve?
    44   firstBelow = -1;
    45   lastAbove = -1;
     60  int firstBelow = -1;
     61  int lastAbove = -1;
    4662  for (int i = 0; i < flux->n; i++) {
    4763    if (flux->data.F32[i] > fluxRP) lastAbove = i;
    4864    if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i;
     65  }
     66  // if we don't go out far enough, we have a problem...
     67  if (lastAbove == radius->n - 1) {
     68    psTrace ("psphot", 5, "did not go out far enough to reach petrosian radius...");
     69    // XXX skip object? raise a flag ?
     70    return false;
     71  }
     72  if (firstBelow < 0) {
     73    psTrace ("psphot", 5, "did not go out far enough to bound petrosian radius");
     74    // XXX raise a flag ?
     75    return false;
    4976  }
    5077
    … …  
    78105  source->extpars->petrosian->radErr = radErr;
    79106
     107  psTrace ("psphot", 5, "Petrosian flux:%f +/- %f @ %f +/- %f for %f, %f\n",
     108           source->extpars->petrosian->mag, source->extpars->petrosian->magErr,
     109           source->extpars->petrosian->rad, source->extpars->petrosian->radErr,
     110           source->peak->xf, source->peak->yf);
     111
    80112  return true;
    81113

  • trunk/psphot/src/psphotRadialProfile.c

    r15819 r17396  
    4141
    4242    int n = 0;
    43     float Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
    44     float Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
     43   
     44    float Xo = 0.0;
     45    float Yo = 0.0;
     46
     47    if (source->modelEXT) {
     48      Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
     49      Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
     50    } else {
     51      Xo = source->peak->xf - source->pixels->col0;
     52      Yo = source->peak->yf - source->pixels->row0;
     53    }
    4554    for (int iy = 0; iy < source->pixels->numRows; iy++) {
    4655        for (int ix = 0; ix < source->pixels->numCols; ix++) {

  • trunk/psphot/src/psphotReadout.c

    r17112 r17396  
    117117    psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
    118118
    119     // XXX test the CR/EXT measurement
    120     // XXX we need to call this here and after the second pass: option for starting number?
     119    // identify CRs and extended sources
    121120    psphotSourceSize (readout, sources, recipe, 0);
    122121    if (!strcasecmp (breakPt, "ENSEMBLE")) {
    … …  
    127126    psphotBlendFit (readout, sources, recipe, psf);
    128127
    129     // replace all sources (make this part of psphotFitSourcesLinear?)
     128    // replace all sources
    130129    psphotReplaceAll (sources, recipe);
    131130
    132131    // linear fit to include all sources
    133132    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
     133
     134    // if we only do one pass, skip to extended source analysis
    134135    if (!strcasecmp (breakPt, "PASS1")) {
    135         goto finish;
     136        goto pass1finish;
    136137    }
    137138
    … …  
    145146    psphotAddNoise (readout, sources, recipe);
    146147
     148    // find fainter sources (pass 2)
    147149    detections = psphotFindDetections (detections, readout, recipe);
    148150
    149     // remove noise for subtracted objects
     151    // remove noise for subtracted objects (ie, return to normal noise level)
    150152    psphotSubNoise (readout, sources, recipe);
    151153
    … …  
    172174    psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
    173175
     176pass1finish:
     177
    174178    // measure source size for the remaining sources
    175179    psphotSourceSize (readout, sources, recipe, 0);
    176180
    177     psphotExtendedSources (readout, sources, recipe);
     181    psphotExtendedSourceAnalysis (readout, sources, recipe);
     182
     183    psphotExtendedSourceFits (readout, sources, recipe);
    178184
    179185finish:

  • trunk/psphot/src/psphotSourceFits.c

    r16820 r17396  
    220220    pmSource *tmpSrc = pmSourceAlloc ();
    221221
    222     pmModel *EXT = psphotFitEXT (readout, source, maskVal);
     222    pmModel *EXT = psphotFitEXT (readout, source, modelTypeEXT, maskVal);
    223223    okEXT = psphotEvalEXT (tmpSrc, EXT);
    224224    chiEXT = EXT->chisq / EXT->nDOF;
    … …  
    267267    // save new model
    268268    source->modelEXT = EXT;
     269    source->type = PM_SOURCE_TYPE_EXTENDED;
     270    source->mode |= PM_SOURCE_MODE_EXTMODEL;
    269271
    270272    // build cached model and subtract
    … …  
    349351}
    350352
    351 pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, psMaskType maskVal) {
     353pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal) {
    352354
    353355    NfitEXT ++;
    354356
    355357    // use the source moments, etc to guess basic model parameters
    356     pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT);
     358    pmModel *EXT = pmSourceModelGuess (source, modelType);
    357359    PS_ASSERT (EXT, NULL);
    358360

  • trunk/psphot/src/psphotSourceSize.c

    r17113 r17396  
    11# include "psphotInternal.h"
    22# include <gsl/gsl_sf_gamma.h>
     3
     4float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro);
    35
    46// we need to call this function after sources have been fitted to the PSF model and
    … …  
    1416    psTimerStart ("psphot");
    1517
    16     float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CRNSIGMA.LIMIT");
     18    float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT");
     19    assert (status);
     20
     21    float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT");
    1722    assert (status);
    1823
    … …  
    2429        if (isfinite(source->crNsigma)) continue;
    2530
    26         source->crNsigma  = -1.0;
    27         source->extNsigma = -1.0;
    28 
    2931        // source must have been subtracted
    30         source->crNsigma  = -3.0;
    3132        if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
    3233
    33         psF32 **resid = source->pixels->data.F32;
     34        psF32 **resid  = source->pixels->data.F32;
    3435        psF32 **weight = source->weight->data.F32;
    35         psU8 **mask = source->maskObj->data.U8;
     36        psU8 **mask    = source->maskObj->data.U8;
     37
     38        // check for extendedness: measure the delta flux significance at the 1 sigma contour
     39        source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, source->modelPSF, 1.0);
     40
     41        // XXX prevent a source from being both CR and EXT?
     42        if (source->extNsigma > EXT_NSIGMA_LIMIT) {
     43          source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
     44        }
    3645
    3746        int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
    3847        int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
    3948
    40         // skip sources which are too close to a boundary
    41         source->crNsigma  = -4.0;
     49        // XXX for now, skip sources which are too close to a boundary
     50        // XXX raise a flag?
    4251        if (xPeak < 1) continue;
    4352        if (xPeak > source->pixels->numCols - 2) continue;
    … …  
    4655
    4756        // XXX for now, just skip any sources with masked pixels
    48         source->crNsigma  = -5.0;
     57        // XXX raise a flag?
    4958        bool keep = true;
    5059        for (int iy = -1; (iy <= +1) && keep; iy++) {
    … …  
    113122        }
    114123        source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
    115         source->extNsigma = (nEXT > 0) ? fabs(fEXT) / nEXT : 0.0;
    116124        // NOTE: abs needed to make the Nsigma value positive
    117125
    … …  
    122130
    123131        if (source->crNsigma > CR_NSIGMA_LIMIT) {
    124           source->mode |= PM_SOURCE_MODE_CRLIMIT;
     132          source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    125133        }
    126134    }
    … …  
    169177 */
    170178
     179
     180// given the PSF ellipse parameters, navigate around the 1sigma contour, return the total
     181// deviation in sigmas.  This is measure on the residual image - should we ignore negative
     182// deviations?
     183float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro) {
     184
     185    psF32 *PAR = model->params->data.F32;
     186
     187    // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY;
     188    // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0;
     189    // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2];
     190    // y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A];
     191
     192    // min/max value of x is where T -> 0
     193    // solve this for x2:
     194    float Q = Ro * PS_SQR(PAR[PM_PAR_SXX]) / (1.0 - PS_SQR(PAR[PM_PAR_SXX]*PAR[PM_PAR_SYY]*PAR[PM_PAR_SXY]) / 4.0);
     195    if (Q < 0.0) return NAN; // ellipse is imaginary
     196
     197    int xMax = sqrt(Q);
     198    int xMin = -1.0*xMax;
     199
     200    int nPts = 0;
     201    float nSigma = 0.0;
     202
     203    for (int x = xMin; x <= xMax; x++) {
     204        float A = PS_SQR (1.0 / PAR[PM_PAR_SYY]);
     205        float B = x * PAR[PM_PAR_SXY];
     206        float C = PS_SQR (x / PAR[PM_PAR_SXX]) - Ro;
     207
     208        float T = PS_SQR(B) - 4*A*C;
     209        if (T < 0.0) continue;
     210   
     211        float yP = (-B + sqrt (T)) / (2.0 * A);
     212        float yM = (-B - sqrt (T)) / (2.0 * A);
     213
     214        int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
     215        int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
     216        int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
     217
     218        if (xPix < 0) continue;
     219        if (xPix >= image->numCols) continue;
     220
     221        if ((yPixM >= 0) && (yPixM < image->numRows)) {
     222            if (!mask || !mask->data.U8[yPixM][xPix]) {
     223                float dSigma = image->data.F32[yPixM][xPix] / sqrt (weight->data.F32[yPixM][xPix]);
     224                nSigma += dSigma;
     225                nPts ++;
     226            }
     227        }
     228       
     229        if (yPixM == yPixP) continue;
     230
     231        if ((yPixP >= 0) && (yPixP < image->numRows)) {
     232            if (!mask || !mask->data.U8[yPixP][xPix]) {
     233                float dSigma = image->data.F32[yPixP][xPix] / sqrt (weight->data.F32[yPixP][xPix]);
     234                nSigma += dSigma;
     235                nPts ++;
     236            }
     237        }
     238    }   
     239    nSigma /= nPts;
     240    return nSigma;
     241}

  • trunk/psphot/src/psphotSourceStats.c

    r16820 r17396  
    3737    for (int i = 0; i < peaks->n; i++) {
    3838
     39        pmPeak *peak = peaks->data[i];
     40        if (peak->assigned) continue;
     41
    3942        // create a new source, add peak
    4043        pmSource *source = pmSourceAlloc();
    41         source->peak = (pmPeak *)psMemIncrRefCounter(peaks->data[i]);
     44        source->peak = psMemIncrRefCounter(peak);
    4245
    4346        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
    4447        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
    4548        if (!strcasecmp (breakPt, "PEAKS")) {
     49            peak->assigned = true;
    4650            psArrayAdd (sources, 100, source);
    4751            psFree (source);
    … …  
    4953        }
    5054
    51         // skip faint sources
     55        // skip faint sources for moments measurement
    5256        if (source->peak->SN < MIN_SN) {
     57            peak->assigned = true;
    5358            psArrayAdd (sources, 100, source);
    5459            psFree (source);
    … …  
    5762
    5863        // measure a local sky value
    59         // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
     64        // the local sky is now ignored; kept here for reference only
    6065        status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, mark);
    6166        if (!status) {
    … …  
    8085        if (status) {
    8186            // add to the source array
     87            peak->assigned = true;
    8288            psArrayAdd (sources, 100, source);
    8389            psFree (source);
    … …  
    9399        if (status) {
    94100            // add to the source array
     101            peak->assigned = true;
    95102            psArrayAdd (sources, 100, source);
    96103            psFree (source);
Note: See TracChangeset for help on using the changeset viewer.