IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 10, 2013, 2:55:11 PM (12 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-20130904

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/psModules

  • trunk/psModules/src/objects/pmPCMdata.c

    r36089 r36375  
    4343
    4444# define USE_DELTA_PSF 0
    45 # define USE_1D_GAUSS 1
    4645
    4746static void pmPCMdataFree (pmPCMdata *pcm) {
    … …  
    5857    psFree (pcm->psfFFT);
    5958    psFree (pcm->constraint);
     59
    6060    psFree (pcm->smdata); // pre-allocated data for psImageSmooth_PreAlloc
     61    psFree (pcm->smdata2d); // pre-allocated data for psImageSmooth_PreAlloc
    6162    return;
    6263}
    … …  
    8889    }
    8990
     91    pcm->smdata = NULL;
     92    pcm->smdata2d = NULL;
     93
    9094    pcm->modelConv = NULL;
    9195    pcm->psf = NULL;
    … …  
    9498    pcm->nDOF = 0;
    9599
     100    pcm->poissonErrors = true;
     101
    96102    // full convolution with the PSF is expensive.  if we have to save time, we can do a 1D
    97103    // convolution with a Gaussian approximation to the kernel
    98104    pcm->use1Dgauss = false;
    99     pcm->nsigma = 3.0;
     105    pcm->nsigma = NAN; // this is set to something defined by the user
    100106    pcm->sigma = 1.0; // this should be set to something sensible when the psf is known
    101107
    … …  
    248254}
    249255
     256static int modelType_GAUSS = -1;
     257static int modelType_PS1_V1 = -1;
     258
     259// generate a Gaussian smoothing kernel for supplied sigma.  sigma here does not need to match
     260// that used to allocate the structure, but it is recommended
     261bool psImageSmoothCacheKernel_PS1_V1 (psImageSmoothCacheData *smdata, float sigma, float kappa) {
     262    // check for NULL structure elements?
     263
     264    int size = smdata->Nrange;
     265
     266    psFree (smdata->kernel);
     267    smdata->kernel = psVectorAlloc(2 * smdata->Nrange + 1, PS_TYPE_F32);
     268
     269    double sum = 0.0;                   // Sum of Gaussian, for normalization
     270    double factor = 1.0 / (sigma * M_SQRT2);    // Multiplier for i -> z
     271
     272    // PS1_V1 is a power-law with fitted linear term:
     273    // 1 / (1 + kappa z + z^1.666)  where z = (r/sigma)^2
     274
     275    // generate the kernel (not normalized)
     276    for (int i = -size, j = 0; i <= size; i++, j++) {
     277        float z = PS_SQR(i * factor);
     278        sum += smdata->kernel->data.F32[j] = 1.0 / (1 + kappa * z + pow(z,1.666));
     279    }
     280
     281    // renormalize kernel to integral of 1.0
     282    for (int i = 0; i < 2 * size + 1; i++) {
     283        smdata->kernel->data.F32[i] /= sum;
     284    }
     285
     286    return true;
     287}
     288
     289psImageSmoothCacheData *psImageSmoothCacheSetKernel (float *sigma, float *kappa, float nsigma, psImage *flux, pmModel *modelPSF) {
     290
     291    psAssert (modelPSF, "psf model must be defined");
     292   
     293    psEllipseAxes axes;
     294    bool useReff = pmModelUseReff (modelPSF->type);
     295    psF32 *PAR = modelPSF->params->data.F32;
     296    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
     297   
     298    *sigma = NAN;
     299    *kappa = NAN;
     300
     301    // XXX need to do this more carefully
     302    if (modelPSF->type == modelType_GAUSS) {
     303        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
     304        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
     305        *sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
     306    }
     307    if (modelPSF->type == modelType_PS1_V1) {
     308        *sigma = 0.5 * (axes.major + axes.minor);
     309        *kappa = PAR[PM_PAR_7];
     310    }
     311    psAssert (isfinite(*sigma), "invalid model type");
     312
     313    // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
     314    psImageSmoothCacheData *smdata = psImageSmoothCacheAlloc (flux, *sigma, nsigma);
     315
     316    if (modelPSF->type == modelType_GAUSS) {
     317        psImageSmoothCacheKernel_Gauss (smdata, *sigma);
     318    }
     319    if (modelPSF->type == modelType_PS1_V1) {
     320        psImageSmoothCacheKernel_PS1_V1 (smdata, *sigma, *kappa);
     321    }
     322
     323    return smdata;
     324}
     325
     326psImageSmooth2dCacheData *psImageSmooth2dCacheSetKernel (float *sigma, float *kappa, float nsigma, psImage *flux, pmModel *modelPSF) {
     327
     328    psAssert (modelPSF, "psf model must be defined");
     329   
     330    psEllipseAxes axes;
     331    bool useReff = pmModelUseReff (modelPSF->type);
     332    psF32 *PAR = modelPSF->params->data.F32;
     333    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
     334   
     335    *sigma = NAN;
     336    *kappa = NAN;
     337
     338    // XXX need to do this more carefully
     339    if (modelPSF->type == modelType_GAUSS) {
     340        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
     341        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
     342        *sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
     343    }
     344    if (modelPSF->type == modelType_PS1_V1) {
     345        *sigma = 0.5 * (axes.major + axes.minor);
     346        *kappa = PAR[PM_PAR_7];
     347    }
     348    psAssert (isfinite(*sigma), "invalid model type");
     349
     350    // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
     351    psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc (nsigma);
     352
     353    if (modelPSF->type == modelType_GAUSS) {
     354        psImageSmooth2dCacheKernel_Gauss (smdata, *sigma);
     355    }
     356    if (modelPSF->type == modelType_PS1_V1) {
     357        psImageSmooth2dCacheKernel_PS1_V1 (smdata, *sigma, *kappa);
     358    }
     359
     360    return smdata;
     361}
     362
    250363pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize) {
    251364
    252     // make sure we save a cached copy of the psf flux
    253     pmSourceCachePSF (source, maskVal);
    254 
    255     // convert the cached cached psf model for this source to a psKernel
    256     psKernel *psf = pmPCMkernelFromPSF (source, psfSize);
    257     if (!psf) {
    258         // NOTE: this only happens if the source is too close to an edge
    259         model->flags |= PM_MODEL_STATUS_BADARGS;
    260         return NULL;
    261     }
    262 
    263 # if (USE_DELTA_PSF)
    264     psImageInit (psf->image, 0.0);
    265     psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0;
    266 # endif
     365    modelType_GAUSS = pmModelClassGetType ("PS_MODEL_GAUSS");
     366    modelType_PS1_V1 = pmModelClassGetType ("PS_MODEL_PS1_V1");
    267367
    268368    // count the number of unmasked pixels:
    … …  
    298398    if (nPix <  nParams + 1) {
    299399        psTrace ("psModules.objects", 4, "insufficient valid pixels\n");
    300         psFree (psf);
    301400        psFree (constraint);
    302401        model->flags |= PM_MODEL_STATUS_BADARGS;
    … …  
    306405    // generate PCM data storage structure
    307406    pmPCMdata *pcm = pmPCMdataAlloc (params, constraint->paramMask, source);
    308 
    309     pcm->psf = psf;
    310407    pcm->modelConv = psMemIncrRefCounter(model);
    311408    pcm->constraint = constraint;
     409
     410    pcm->poissonErrors = fitOptions->poissonErrors;
     411    pcm->nsigma = fitOptions->nsigma;
    312412
    313413    pcm->nPix = nPix;
    … …  
    316416
    317417# if (USE_1D_GAUSS)
    318     pmModel *modelPSF = source->modelPSF;
    319     psAssert (modelPSF, "psf model must be defined");
    320    
    321     psEllipseAxes axes;
    322     bool useReff = pmModelUseReff (modelPSF->type);
    323     psF32 *PAR = modelPSF->params->data.F32;
    324     pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
    325    
    326     float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
    327     float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
    328418
    329419    pcm->use1Dgauss = true;
    330     pcm->sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
    331     pcm->nsigma = 2.0;
    332 
    333     pcm->smdata = psImageSmooth_PreAlloc_DataAlloc (source->pixels, pcm->sigma, pcm->nsigma);
     420    if (USE_1D_CACHE) {
     421        pcm->smdata = psImageSmoothCacheSetKernel (&pcm->sigma, &pcm->kappa, pcm->nsigma, source->pixels, source->modelPSF);
     422    } else {
     423        pcm->smdata2d = psImageSmooth2dCacheSetKernel (&pcm->sigma, &pcm->kappa, pcm->nsigma, source->pixels, source->modelPSF);
     424    }
     425
    334426# else
     427    // make sure we save a cached copy of the psf flux
     428    pmSourceCachePSF (source, maskVal);
     429
     430    // convert the cached cached psf model for this source to a psKernel
     431    psKernel *psf = pmPCMkernelFromPSF (source, psfSize);
     432    if (!psf) {
     433        // NOTE: this only happens if the source is too close to an edge
     434        model->flags |= PM_MODEL_STATUS_BADARGS;
     435        return NULL;
     436    }
     437
     438# if (USE_DELTA_PSF)
     439    psImageInit (psf->image, 0.0);
     440    psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0;
     441# endif
     442    pcm->psf = psf;
    335443    pcm->smdata = NULL;
    336444# endif
    … …  
    395503            pcm->dmodelsConvFlux->data[n] = psImageCopy (pcm->dmodelsConvFlux->data[n], source->pixels, PS_TYPE_F32);
    396504        }
    397         psFree(pcm->smdata);
    398         pcm->smdata = psImageSmooth_PreAlloc_DataAlloc (source->pixels, pcm->sigma, pcm->nsigma);
     505
     506        // If we have changed the window, we need to redefine the smoothing target vectors (but pcm->sigma,kappa,nsigma remain)
     507        if (USE_1D_CACHE) {
     508            psFree(pcm->smdata);
     509            pcm->smdata = psImageSmoothCacheAlloc (source->pixels, pcm->sigma, pcm->nsigma);
     510
     511            pmModel *modelPSF = source->modelPSF;
     512            if (modelPSF->type == modelType_GAUSS) {
     513                psImageSmoothCacheKernel_Gauss (pcm->smdata, pcm->sigma);
     514            }
     515            if (modelPSF->type == modelType_PS1_V1) {
     516                psImageSmoothCacheKernel_PS1_V1 (pcm->smdata, pcm->sigma, pcm->kappa);
     517            }
     518        } else {
     519            psFree(pcm->smdata2d);
     520            pcm->smdata2d = psImageSmooth2dCacheAlloc (pcm->nsigma);
     521
     522            pmModel *modelPSF = source->modelPSF;
     523            if (modelPSF->type == modelType_GAUSS) {
     524                // psImageSmooth2dCacheKernel_Gauss (pcm->smdata2d, pcm->sigma);
     525            }
     526            if (modelPSF->type == modelType_PS1_V1) {
     527                psImageSmooth2dCacheKernel_PS1_V1 (pcm->smdata2d, pcm->sigma, pcm->kappa);
     528            }
     529        }
    399530    }
    400531
    … …  
    403534
    404535// construct a realization of the source model
    405 bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize) {
     536bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize, float nsigma) {
    406537
    407538    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    420551    // convolve the model image with the PSF
    421552    if (USE_1D_GAUSS) {
    422         // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
    423         // * the model flux is not masked
    424         // * threading takes place above this level
    425553       
    426         // define the Gauss parameters from the psf
    427         pmModel *modelPSF = source->modelPSF;
    428         psAssert (modelPSF, "psf model must be defined");
    429    
    430         psEllipseAxes axes;
    431         bool useReff = pmModelUseReff (modelPSF->type);
    432         psF32 *PAR = modelPSF->params->data.F32;
    433         pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
    434    
    435         float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
    436         float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
    437 
    438         float sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
    439         float nsigma = 2.0;
    440 
    441         psImageSmooth (source->modelFlux, sigma, nsigma);
     554        float sigma = NAN;
     555        float kappa = NAN;
     556
     557        if (USE_1D_CACHE) {
     558            psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, nsigma, source->modelFlux, source->modelPSF);
     559            psImageSmoothCache_F32 (source->modelFlux, smdata);
     560            psFree (smdata);
     561        } else {
     562            psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheSetKernel (&sigma, &kappa, nsigma, source->modelFlux, source->modelPSF);
     563            psImageSmooth2dCache_F32 (source->modelFlux, smdata);
     564            psFree (smdata);
     565        }
     566        // old call: psImageSmooth (source->modelFlux, sigma, nsigma);
    442567    } else {
    443568        // make sure we save a cached copy of the psf flux
    … …  
    459584
    460585// construct a realization of the source model
    461 bool pmPCMMakeModel (pmSource *source, pmModel *model, psImageMaskType maskVal, int psfSize) {
     586bool pmPCMMakeModel (pmSource *source, pmModel *model, float Nsigma, psImageMaskType maskVal, int psfSize) {
    462587
    463588    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    468593
    469594    // modelFlux always has unity normalization (I0 = 1.0)
    470     pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
     595    // pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
     596    pmModelAdd (source->modelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_SKY | PM_MODEL_OP_NORM, maskVal);
    471597
    472598    // convolve the model image with the PSF
    473599    if (USE_1D_GAUSS) {
    474         // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
    475         // * the model flux is not masked
    476         // * threading takes place above this level
    477        
    478         // define the Gauss parameters from the psf
    479         pmModel *modelPSF = source->modelPSF;
    480         psAssert (modelPSF, "psf model must be defined");
    481    
    482         psEllipseAxes axes;
    483         bool useReff = pmModelUseReff (modelPSF->type);
    484         psF32 *PAR = modelPSF->params->data.F32;
    485         pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
    486    
    487         float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
    488         float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
    489 
    490         float sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
    491         float nsigma = 2.0;
    492 
    493         psImageSmooth (source->modelFlux, sigma, nsigma);
     600
     601        float sigma = NAN;
     602        float kappa = NAN;
     603
     604        if (USE_1D_CACHE) {
     605            psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, Nsigma, source->modelFlux, source->modelPSF);
     606            psImageSmoothCache_F32 (source->modelFlux, smdata);
     607            psFree (smdata);
     608        } else {
     609            psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheSetKernel (&sigma, &kappa, Nsigma, source->modelFlux, source->modelPSF);
     610            psImageSmooth2dCache_F32 (source->modelFlux, smdata);
     611            psFree (smdata);
     612        }
     613        // old call: psImageSmooth (source->modelFlux, sigma, nsigma);
    494614    } else {
    495615        // make sure we save a cached copy of the psf flux
    … …  
    509629    return true;
    510630}
    511 
Note: See TracChangeset for help on using the changeset viewer.