IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 24, 2013, 12:09:27 PM (13 years ago)
Author:
eugene
Message:

first version of psphotFullForce which builds

Location:
branches/eam_branches/ipp-20130904/psphot/src
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130904/psphot/src

    • Property svn:ignore

      •  

        old new  
        2424psphotModelTest
        2525psphotMinimal
         26psphotFullForce

  • branches/eam_branches/ipp-20130904/psphot/src/Makefile.am

    r34317 r36148  
    2525libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    2626
    27 bin_PROGRAMS = psphot psphotForced psphotMinimal psphotMakePSF psphotStack psphotModelTest
     27bin_PROGRAMS = psphot psphotForced psphotFullForce psphotMinimal psphotMakePSF psphotStack psphotModelTest
    2828# bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy
    2929
    … …  
    3535psphotForced_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    3636psphotForced_LDADD = libpsphot.la
     37
     38psphotFullForce_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     39psphotFullForce_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     40psphotFullForce_LDADD = libpsphot.la
    3741
    3842psphotMinimal_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    … …  
    7781        psphotForced.c             \
    7882        psphotForcedArguments.c    \
     83        psphotParseCamera.c        \
     84        psphotImageLoop.c          \
     85        psphotMosaicChip.c         \
     86        psphotCleanup.c
     87
     88# generalized forced photometry (including Kron, Petro, and Models) of specified positions given a specified psf
     89psphotFullForce_SOURCES = \
     90        psphotFullForce.c             \
     91        psphotFullForceArguments.c         \
    7992        psphotParseCamera.c        \
    8093        psphotImageLoop.c          \
    … …  
    155168        psphotReadoutMinimal.c         \
    156169        psphotForcedReadout.c          \
     170        psphotFullForceReadout.c       \
    157171        psphotMakePSFReadout.c         \
    158172        psphotModelBackground.c        \
    … …  
    198212        psphotRadialPlot.c             \
    199213        psphotKronIterate.c            \
     214        psphotKronFlux.c            \
     215        psphotPetroFlux.c            \
     216        psphotGalaxyShape.c            \
    200217        psphotRadialProfileWings.c     \
    201218        psphotDeblendSatstars.c        \

  • branches/eam_branches/ipp-20130904/psphot/src/psphot.h

    r36146 r36148  
    342342bool psphotForcedReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
    343343
     344pmConfig *psphotFullForceArguments(int argc, char **argv);
     345bool psphotFullForceReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
     346
    344347pmConfig *psphotMinimalArguments(int argc, char **argv);
    345348bool psphotReadoutMinimal(pmConfig *config, const pmFPAview *view, const char *filerule);
    … …  
    498501bool psphotKronIterate_Threaded (psThreadJob *job);
    499502
     503bool psphotKronFlux (pmConfig *config, const pmFPAview *view, const char *filerule);
     504bool psphotKronFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char *filerule, pmReadout *readout, psArray *sources);
     505bool psphotKronFlux_Threaded (psThreadJob *job);
     506bool psphotKronFluxSource (pmSource *source, psImageMaskType maskVal);
     507
     508bool psphotPetroFlux (pmConfig *config, const pmFPAview *view, const char *filerule);
     509bool psphotPetroFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources);
     510bool psphotPetroFlux_Threaded (psThreadJob *job);
     511bool psphotPetroFluxSource (pmSource *source, psImageMaskType maskVal);
     512
     513bool psphotGalaxyShape (pmConfig *config, const pmFPAview *view, const char *filerule);
     514bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf);
     515bool psphotGalaxyShape_Threaded (psThreadJob *job);
     516bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, int psfSize);
     517bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize);
     518
    500519bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule);
    501520bool psphotRadialProfileWingsReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources);

  • branches/eam_branches/ipp-20130904/psphot/src/psphotBlendFit.c

    r36086 r36148  
    273273        pmSource *source = sources->data[i];
    274274
    275 # if (0)
    276         int TEST_ON = false;
     275# if (PS_TRACE_ON)
    277276# define TEST_X 653
    278277# define TEST_Y 466
    279         if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     278# define TESTING 1
     279        int TEST_ON = false;
     280        if (TESTING && (fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    280281            fprintf (stderr, "test object\n");
    281282            psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);

  • branches/eam_branches/ipp-20130904/psphot/src/psphotFullForceReadout.c

    r36146 r36148  
    5757
    5858    // generate a psf model for any readouts which need one
    59     psphotFullForcePSF (config, view, filerule);
     59    // psphotFullForcePSF (config, view, filerule);
     60    // XXX I'm not sure we need a different algorithm here or not : we are supplying the
     61    // sources; we should mark with a flag bit the ones we actually want to use as PSF
     62    // stars (this means we need to supply this info in the load).
     63    psphotChoosePSF (config, view, filerule, true);
    6064
    6165    // Construct an initial model for each object, set the radius to fitRadius, set circular
    … …  
    7377
    7478    // measure kron fluxes
    75     psphotKronFlux (config, view, filerule, false);
     79    psphotKronFlux (config, view, filerule);
    7680
    7781    // measure petro fluxes

  • branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c

    r36146 r36148  
    11# include "psphotInternal.h"
    22
    3 bool psphotGalaxyShape (pmConfig *config, const pmFPAview *view, const char *filerule, int pass)
     3bool psphotGalaxyShape (pmConfig *config, const pmFPAview *view, const char *filerule)
    44{
    55    bool status = true;
    … …  
    3838        // psAssert (psf, "missing psf?");
    3939
    40         if (!psphotGalaxyShapeReadout (config, recipe, view, filerule, readout, sources, psf, i, pass)) {
     40        if (!psphotGalaxyShapeReadout (config, recipe, view, filerule, readout, sources, psf)) {
    4141            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
    4242            return false;
    … …  
    4646}
    4747
    48 bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) {
     48bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf) {
    4949
    5050    bool status = false;
    … …  
    5555    }
    5656
    57     psTimerStart ("psphot.kron");
     57    psTimerStart ("psphot.galaxy");
    5858
    5959    // determine the number of allowed threads
    … …  
    6363    }
    6464
    65     bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
    66     if (!status) {
    67         KRON_APPLY_WEIGHT = true;
    68     }
    69 
    70     bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
    71     if (!status) {
    72         KRON_APPLY_WINDOW = false;
    73     }
    74 
    7565    // bit-masks to test for good/bad pixels
    7666    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    … …  
    8373    // maskVal is used to test for rejected pixels, and must include markVal
    8474    maskVal |= markVal;
     75
     76    // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
     77    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
     78    assert (status);
     79
     80    // source fitting parameters for extended source fits
     81    int fitIter = psMetadataLookupS32(&status, recipe, "EXT_FIT_ITER"); // Max number of fit iterations
     82    assert (status && fitIter > 0);
     83
     84    float fitMinTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MIN_TOL"); // Fit tolerance
     85    if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
     86        fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance
     87        if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
     88            psAbort("PSF_FIT_MIN_TOL (and PSF_FIT_TOL) not defined or positive");
     89        }
     90    }
     91
     92    float fitMaxTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_TOL"); // Fit tolerance
     93    if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) {
     94        fitMaxTol = 1.0;
     95    }
     96
     97    bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
     98    if (!status) {
     99        // default to the old method (chisqConvergence)
     100        chisqConvergence = true;
     101    }
     102
     103    int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance
     104    if (!status) {
     105        // default to the old method (chisqConvergence)
     106        gainFactorMode = 0;
     107    }
     108
     109    // perform full extended source non-linear fits?
     110    bool isInteractive = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_INTERACTIVE");
     111    if (!status) isInteractive = false;
     112
     113    // Define source fitting parameters for extended source fits
     114    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
     115    fitOptions->mode           = PM_SOURCE_FIT_EXT_AND_SKY;
     116    fitOptions->saveCovariance = true;  // XXX make this a user option?
     117    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
     118    fitOptions->nIter          = fitIter;
     119    fitOptions->minTol         = fitMinTol;
     120    fitOptions->maxTol         = fitMaxTol;
     121
     122    fitOptions->gainFactorMode   = gainFactorMode;
     123    fitOptions->chisqConvergence = chisqConvergence;
     124    fitOptions->isInteractive    = isInteractive;
    85125
    86126    // source analysis is done in S/N order (brightest first)
    … …  
    105145
    106146            // allocate a job -- if threads are not defined, this just runs the job
    107             psThreadJob *job = psThreadJobAlloc ("PSPHOT_KRON_ITERATE");
     147            psThreadJob *job = psThreadJobAlloc ("PSPHOT_GALAXY_SHAPES");
    108148
    109149            psArrayAdd(job->args, 1, readout);
    110150            psArrayAdd(job->args, 1, cells->data[j]); // sources
    111             PS_ARRAY_ADD_SCALAR(job->args, markVal,            PS_TYPE_IMAGE_MASK);
    112             PS_ARRAY_ADD_SCALAR(job->args, maskVal,            PS_TYPE_IMAGE_MASK);
    113             PS_ARRAY_ADD_SCALAR(job->args, pass,               PS_TYPE_S32);
     151            psArrayAdd(job->args, 1, fitOptions);
     152            PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK);
     153            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
     154            PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32);
    114155
    115156// set this to 0 to run without threading
    … …  
    146187    psFree (cellGroups);
    147188
    148     psLogMsg ("psphot.kron", PS_LOG_WARN, "measure kron fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
     189    psLogMsg ("psphot.galaxy", PS_LOG_WARN, "measure galaxy shapes : %f sec for %ld objects\n", psTimerMark ("psphot.galaxy"), sources->n);
    149190    return true;
    150191}
    … …  
    152193bool psphotGalaxyShape_Threaded (psThreadJob *job) {
    153194
    154     pmReadout *readout              = job->args->data[0];
    155     psArray *sources                = job->args->data[2];
    156     psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    157     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
     195    pmReadout *readout             = job->args->data[0];
     196    psArray *sources               = job->args->data[1];
     197    pmSourceFitOptions *fitOptions = job->args->data[2];
     198    psImageMaskType markVal        = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
     199    psImageMaskType maskVal        = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
     200    int psfSize                    = PS_SCALAR_VALUE(job->args->data[5],S32);
    158201
    159202    for (int i = 0; i < sources->n; i++) {
    160 
    161         // XXX what do i need for the kron flux (given a supplied kron radius?)
    162203
    163204        pmSource *source = sources->data[i];
    … …  
    184225
    185226        // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    186         bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
     227        pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    187228        psAssert (source->pixels, "WTF?");
    188229
    189230        // this function populates moments->Mrf,GalaxyShape,GalaxyShapeErr
    190231        // do the following for a set of shapes (Ex,Ey)
    191         psphotGalaxyShapeSource (source, maskVal);
     232        psphotGalaxyShapeGrid (source, fitOptions, maskVal, psfSize);
    192233        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    193234
    … …  
    198239}
    199240
    200 bool psphotGalaxyShapeGrid (pmSource *source, psImageMaskType maskVal, float psfSize) {
    201 
    202     pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
    203 
    204   // I have some source guess (e0, e1, e2)
    205   psEllipsePol guessPol;
    206 
    207   // I am going to retain the value of e0, and grid search around e1,e2
    208   // I need to study a valid range (and distribution?) for e1,e2)
    209 
    210   for (float de1 = -0.2; de1 <= +0.2; de1 += 0.1) {
    211     for (float de2 = -0.2; de2 <= +0.2; de2 += 0.1) {
     241bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, int psfSize) {
     242
     243    // allocate the model
     244    // XXX need to figure out how to get the model into the program
     245    pmModelType modelType = pmModelClassGetType ("PS_MODEL_EXP");
     246
     247    pmModel *model = pmModelAlloc(modelType);
     248    if (!model) {
     249        return NULL;
     250    }
     251
     252    pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
     253
     254    psF32 *PAR = pcm->modelConv->params->data.F32;
     255
     256    // I have some source guess (e0, e1, e2)
     257    psEllipsePol guessPol;
     258
     259    // I am going to retain the value of e0, and grid search around e1,e2
     260    // I need to study a valid range (and distribution?) for e1,e2)
     261
     262    for (float de1 = -0.2; de1 <= +0.2; de1 += 0.1) {
     263        for (float de2 = -0.2; de2 <= +0.2; de2 += 0.1) {
    212264 
    213       psEllipsePol testPol = guessPol;
    214       testPol.e1 += de1;
    215       testPol.e2 += de2;
    216 
    217       psEllipseShape shape;
    218       psEllipsePolToShape (&shape, testPol);
    219      
    220       PAR[PM_PAR_SXX] = shape.sxx;
    221       PAR[PM_PAR_SXY] = shape.sxy;
    222       PAR[PM_PAR_SYY] = shape.syy;
    223 
    224       // where do the results go??
    225       psphotGalaxyShapeSource (pcm, source, maskVal);
    226     }
    227   }
    228 
    229   return true;
     265            psEllipsePol testPol = guessPol;
     266            testPol.e1 += de1;
     267            testPol.e2 += de2;
     268
     269            psEllipseAxes axes = psEllipsePolToAxes (testPol, 0.1);
     270
     271            pmPSF_AxesToModel (PAR, axes, modelType);
     272
     273            // where do the results go??
     274            psphotGalaxyShapeSource (pcm, source, maskVal, psfSize);
     275        }
     276    }
     277
     278    return true;
    230279}
    231280
    232281// fit the given model to the source and find chisq & normalization
    233282// XXX is this a single-component model? sersic with a supplied index, Reff, axis ratio, and theta?
    234 bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal) {
     283bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize) {
    235284
    236285    PS_ASSERT_PTR_NON_NULL(source, false);
    237286    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    238287    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    239 
    240     psF32 *PAR = pcm->modelConv->params->data.F32;
    241288
    242289    // generate the modelFlux
    … …  
    249296
    250297    for (int iy = 0; iy < source->pixels->numRows; iy++) {
    251       for (int ix = 0; ix < source->pixels->numCols; ix++) {
    252         // skip masked points
    253         if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
    254           continue;
     298        for (int ix = 0; ix < source->pixels->numCols; ix++) {
     299            // skip masked points
     300            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
     301                continue;
     302            }
     303            // skip zero-variance points
     304            if (source->variance->data.F32[iy][ix] == 0) {
     305                continue;
     306            }
     307            // skip nan value points
     308            if (!isfinite(source->pixels->data.F32[iy][ix])) {
     309                continue;
     310            }
     311
     312            float fy = source->pixels->data.F32[iy][ix];
     313            float fm = source->modelFlux->data.F32[iy][ix];
     314            float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
     315
     316            YY += PS_SQR(fy) * wt;
     317            YM += fm * fy * wt;
     318            MM += PS_SQR(fm) * wt;
    255319        }
    256         // skip zero-variance points
    257         if (source->variance->data.F32[iy][ix] == 0) {
    258           continue;
    259         }
    260         // skip nan value points
    261         if (!isfinite(source->pixels->data.F32[iy][ix])) {
    262           continue;
    263         }
    264 
    265         float fy = source->pixels->data.F32[iy][ix];
    266         float fm = source->modelFlux->data.F32[iy][ix];
    267         float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
    268 
    269         YY += PS_SQR(fy) * wt;
    270         YM += fm * fy * wt;
    271         MM += PS_SQR(fm) * wt;
    272       }
    273320    }
    274321
    … …  
    276323    float Chisq = YY - 2 * Io * YM + Io * Io * MM;
    277324
    278     return true;
    279 }
     325    fprintf (stderr, "Io: %f, Chisq: %f\n", Io, Chisq);
     326
     327    return true;
     328}

  • branches/eam_branches/ipp-20130904/psphot/src/psphotKronFlux.c

    r36146 r36148  
    11# include "psphotInternal.h"
    22
    3 bool psphotKronFlux (pmConfig *config, const pmFPAview *view, const char *filerule, int pass)
     3bool psphotKronFlux (pmConfig *config, const pmFPAview *view, const char *filerule)
    44{
    55    bool status = true;
    … …  
    3535        psAssert (sources, "missing sources?");
    3636
    37         pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    38         // psAssert (psf, "missing psf?");
    39 
    40         if (!psphotKronFluxReadout (config, recipe, view, filerule, readout, sources, psf, i, pass)) {
     37        if (!psphotKronFluxReadout (config, recipe, view, filerule, readout, sources)) {
    4138            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
    4239            return false;
    … …  
    4643}
    4744
    48 bool psphotKronFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) {
     45bool psphotKronFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char *filerule, pmReadout *readout, psArray *sources) {
    4946
    5047    bool status = false;
    … …  
    6158    if (!status) {
    6259        nThreads = 0;
    63     }
    64 
    65     bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
    66     if (!status) {
    67         KRON_APPLY_WEIGHT = true;
    68     }
    69 
    70     bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
    71     if (!status) {
    72         KRON_APPLY_WINDOW = false;
    7360    }
    7461
    … …  
    10592
    10693            // allocate a job -- if threads are not defined, this just runs the job
    107             psThreadJob *job = psThreadJobAlloc ("PSPHOT_KRON_ITERATE");
     94            psThreadJob *job = psThreadJobAlloc ("PSPHOT_KRON_FLUX");
    10895
    10996            psArrayAdd(job->args, 1, readout);
    … …  
    11198            PS_ARRAY_ADD_SCALAR(job->args, markVal,            PS_TYPE_IMAGE_MASK);
    11299            PS_ARRAY_ADD_SCALAR(job->args, maskVal,            PS_TYPE_IMAGE_MASK);
    113             PS_ARRAY_ADD_SCALAR(job->args, pass,               PS_TYPE_S32);
    114100
    115101// set this to 0 to run without threading
    … …  
    153139
    154140    pmReadout *readout              = job->args->data[0];
    155     psArray *sources                = job->args->data[2];
    156     psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    157     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
     141    psArray *sources                = job->args->data[1];
     142    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA);
     143    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
    158144
    159145    for (int i = 0; i < sources->n; i++) {
    160 
    161         // XXX what do i need for the kron flux (given a supplied kron radius?)
    162146
    163147        pmSource *source = sources->data[i];
    … …  
    180164        }
    181165
    182         // is the right??
    183         float windowRadius = 6.0*source->moments->Mrf ;
     166        // is this right??
     167        float windowRadius = 2.5*source->moments->Mrf ;
    184168
    185169        // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    186         bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
     170        pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    187171        psAssert (source->pixels, "WTF?");
    188172
    … …  
    205189
    206190    // the peak position is less accurate but less subject to extreme deviations
     191    // XXX do this or just use the peak position?
    207192    float dX = source->moments->Mx - source->peak->xf;
    208193    float dY = source->moments->My - source->peak->yf;
    … …  
    211196    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
    212197
     198    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     199    // xCM, yCM from pixel coords to pixel index here.
     200    psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
     201    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
     202
    213203    // Calculate the Kron magnitude (make this block optional?)
    214204    // XXX set the aperture here
    215     float radKron  = 2.5*Mrf;
     205    float radKron  = 2.5*source->moments->Mrf;
    216206    float radKron2 = radKron*radKron;
    217207
    … …  
    222212    // set vPix to the source pixels (it may have been set to the
    223213    // smoothed image above)
    224     vPix = source->pixels->data.F32;
    225     vWgt = source->variance->data.F32;
    226     vMsk = source->mask->data.F32;
     214    psF32 **vPix = source->pixels->data.F32;
     215    psF32 **vWgt = source->variance->data.F32;
     216    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    227217
    228218    for (psS32 row = 0; row < source->pixels->numRows ; row++) {

  • branches/eam_branches/ipp-20130904/psphot/src/psphotPetroFlux.c

    r36146 r36148  
    11# include "psphotInternal.h"
    22
    3 bool psphotPetroFlux (pmConfig *config, const pmFPAview *view, const char *filerule, int pass)
     3bool psphotPetroFlux (pmConfig *config, const pmFPAview *view, const char *filerule)
    44{
    55    bool status = true;
    … …  
    3535        psAssert (sources, "missing sources?");
    3636
    37         pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    38         // psAssert (psf, "missing psf?");
    39 
    40         if (!psphotPetroFluxReadout (config, recipe, view, filerule, readout, sources, psf, i, pass)) {
     37        if (!psphotPetroFluxReadout (config, recipe, view, filerule, readout, sources)) {
    4138            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
    4239            return false;
    … …  
    4643}
    4744
    48 bool psphotPetroFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) {
     45bool psphotPetroFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources) {
    4946
    5047    bool status = false;
    5148
    5249    if (!sources->n) {
    53         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping kron fluxes");
     50        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping petrosian fluxes");
    5451        return true;
    5552    }
    5653
    57     psTimerStart ("psphot.kron");
     54    psTimerStart ("psphot.petro");
    5855
    5956    // determine the number of allowed threads
    … …  
    6158    if (!status) {
    6259        nThreads = 0;
    63     }
    64 
    65     bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
    66     if (!status) {
    67         KRON_APPLY_WEIGHT = true;
    68     }
    69 
    70     bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
    71     if (!status) {
    72         KRON_APPLY_WINDOW = false;
    7360    }
    7461
    … …  
    10592
    10693            // allocate a job -- if threads are not defined, this just runs the job
    107             psThreadJob *job = psThreadJobAlloc ("PSPHOT_KRON_ITERATE");
     94            psThreadJob *job = psThreadJobAlloc ("PSPHOT_PETRO_FLUX");
    10895
    10996            psArrayAdd(job->args, 1, readout);
    … …  
    11198            PS_ARRAY_ADD_SCALAR(job->args, markVal,            PS_TYPE_IMAGE_MASK);
    11299            PS_ARRAY_ADD_SCALAR(job->args, maskVal,            PS_TYPE_IMAGE_MASK);
    113             PS_ARRAY_ADD_SCALAR(job->args, pass,               PS_TYPE_S32);
    114100
    115101// set this to 0 to run without threading
    … …  
    146132    psFree (cellGroups);
    147133
    148     psLogMsg ("psphot.kron", PS_LOG_WARN, "measure kron fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
     134    psLogMsg ("psphot.petro", PS_LOG_WARN, "measure petro fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.petro"), sources->n);
    149135    return true;
    150136}
    … …  
    153139
    154140    pmReadout *readout              = job->args->data[0];
    155     psArray *sources                = job->args->data[2];
    156     psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    157     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
     141    psArray *sources                = job->args->data[1];
     142    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA);
     143    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
    158144
    159145    for (int i = 0; i < sources->n; i++) {
    160 
    161         // XXX what do i need for the kron flux (given a supplied kron radius?)
    162146
    163147        pmSource *source = sources->data[i];
    … …  
    169153        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
    170154
    171         // XXX where are we storing the supplied kron radius?
     155        // XXX where are we storing the supplied petro radius?
    172156        if (!isfinite(source->moments->Mrf)) continue;
    173157
    … …  
    180164        }
    181165
    182         // is the right??
    183         float windowRadius = 6.0*source->moments->Mrf ;
     166        // is the right??  XXX need the correct location for Rpet
     167        float windowRadius = 2.0*source->moments->Mrf ;
    184168
    185169        // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    186         bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
     170        pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    187171        psAssert (source->pixels, "WTF?");
    188172
    … …  
    211195    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
    212196
     197    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     198    // xCM, yCM from pixel coords to pixel index here.
     199    psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
     200    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
     201
    213202    // Calculate the Petro magnitude (make this block optional?)
    214203    // XXX set the aperture here
    215     float radPetro  = 2.5*Mrf;
     204    float radPetro  = 2.0*source->moments->Mrf;
    216205    float radPetro2 = radPetro*radPetro;
    217206
    … …  
    222211    // set vPix to the source pixels (it may have been set to the
    223212    // smoothed image above)
    224     vPix = source->pixels->data.F32;
    225     vWgt = source->variance->data.F32;
    226     vMsk = source->mask->data.F32;
     213    psF32 **vPix = source->pixels->data.F32;
     214    psF32 **vWgt = source->variance->data.F32;
     215    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    227216
    228217    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    253242
    254243    // return the flux and error to parameters?
    255     source->moments->PetroFlux    = Sum;
    256     source->moments->PetroFluxErr = sqrt(Var);
    257 
    258     return true;
    259 }
     244    // source->moments->PetroFlux    = Sum;
     245    // source->moments->PetroFluxErr = sqrt(Var);
     246
     247    fprintf (stderr, "petro flux: %f +/- %f\n", Sum, sqrt(Var));
     248
     249    return true;
     250}

  • branches/eam_branches/ipp-20130904/psphot/src/psphotSetThreads.c

    r35769 r36148  
    4747    psFree(task);
    4848
     49    task = psThreadTaskAlloc("PSPHOT_KRON_FLUX", 4);
     50    task->function = &psphotKronFlux_Threaded;
     51    psThreadTaskAdd(task);
     52    psFree(task);
     53
     54    task = psThreadTaskAlloc("PSPHOT_PETRO_FLUX", 4);
     55    task->function = &psphotPetroFlux_Threaded;
     56    psThreadTaskAdd(task);
     57    psFree(task);
     58
     59    task = psThreadTaskAlloc("PSPHOT_GALAXY_SHAPE", 6);
     60    task->function = &psphotGalaxyShape_Threaded;
     61    psThreadTaskAdd(task);
     62    psFree(task);
     63
    4964    task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 10);
    5065    task->function = &psphotBlendFit_Threaded;
Note: See TracChangeset for help on using the changeset viewer.