IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36281


Ignore:
Timestamp:
Nov 11, 2013, 11:20:15 AM (13 years ago)
Author:
eugene
Message:

add 2d cached convolution to the convolution concepts (clear this out when I am happy)

Location:
branches/eam_branches/ipp-20130904/psModules/src/objects
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130904/psModules/src/objects/pmPCM_MinimizeChisq.c

    r36276 r36281  
    473473# else
    474474    if (pcm->use1Dgauss) {
    475         // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
    476         // * the model flux is not masked
    477         // * threading takes place above this level
    478         pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type);
    479        
    480         // hard-wire this for now?
    481         psImageSmoothCache_F32 (pcm->modelConvFlux, pcm->smdata1);
    482         psImageSmoothCache_F32 (pcm->modelConvFlux, pcm->smdata2);
     475
     476        if (USE_1D_CACHE) {
     477            // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
     478            // * the model flux is not masked
     479            // * threading takes place above this level
     480            pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type);
     481            psImageSmoothCache_F32 (pcm->modelConvFlux, pcm->smdata);
     482        } else {
     483            pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type);
     484            psImageSmooth2dCache_F32 (pcm->modelConvFlux, pcm->smdata2d);
     485        }
    483486    } else {
    484487        psImageConvolveKernel (pcm->modelConvFlux, pcm->modelFlux, NULL, 0, pcm->psfFFT);
    … …  
    495498# else
    496499        if (pcm->use1Dgauss) {
    497             // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
    498             // * the model flux is not masked
    499             // * threading takes place above this level
    500             dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
    501             psImageSmoothCache_F32 (dmodelConv, pcm->smdata1);
    502             psImageSmoothCache_F32 (dmodelConv, pcm->smdata2);
     500            if (USE_1D_GAUSS) {
     501                // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
     502                // * the model flux is not masked
     503                // * threading takes place above this level
     504                dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
     505                psImageSmoothCache_F32 (dmodelConv, pcm->smdata);
     506            } else {
     507                dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
     508                psImageSmooth2dCache_F32 (dmodelConv, pcm->smdata2d);
     509            }
    503510        } else {
    504511            psImageConvolveKernel (dmodelConv, dmodel, NULL, 0, pcm->psfFFT);
    … …  
    516523
    517524        if (pcm->use1Dgauss) {
    518             // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
    519             // * the model flux is not masked
    520             // * threading takes place above this level
    521             dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
    522             psImageSmoothCache_F32 (dmodelConv, pcm->smdata1);
    523             psImageSmoothCache_F32 (dmodelConv, pcm->smdata2);
     525            if (USE_1D_GAUSS) {
     526                // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
     527                // * the model flux is not masked
     528                // * threading takes place above this level
     529                dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
     530                psImageSmoothCache_F32 (dmodelConv, pcm->smdata);
     531            } else {
     532                dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
     533                psImageSmooth2dCache_F32 (dmodelConv, pcm->smdata2d);
     534            }
    524535        } else {
    525536            psImageConvolveFFT (dmodelConv, dmodel, NULL, 0, pcm->psf);

  • branches/eam_branches/ipp-20130904/psModules/src/objects/pmPCMdata.c

    r36277 r36281  
    4343
    4444# define USE_DELTA_PSF 0
    45 # define USE_1D_GAUSS 1
     45# define KERNEL_NSIGMA 3.0
    4646
    4747static void pmPCMdataFree (pmPCMdata *pcm) {
    … …  
    104104    // convolution with a Gaussian approximation to the kernel
    105105    pcm->use1Dgauss = false;
    106     pcm->nsigma = 5.0;
     106    pcm->nsigma = KERNEL_NSIGMA;
    107107    pcm->sigma = 1.0; // this should be set to something sensible when the psf is known
    108108
    … …  
    288288}
    289289
    290 // set 2 gaussians (extend to more if interesting)
    291 psImageSmoothCacheData *psImageSmoothCacheSetKernel (float *sigma1, float *sigma2, psImage *flux, pmModel *modelPSF) {
     290psImageSmoothCacheData *psImageSmoothCacheSetKernel (float *sigma, float *kappa, float nsigma, psImage *flux, pmModel *modelPSF) {
    292291
    293292    psAssert (modelPSF, "psf model must be defined");
    … …  
    298297    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
    299298   
    300     float nsigma = 5.0;
    301     *sigma1 = NAN;
    302     *sigma2 = NAN;
     299    *sigma = NAN;
     300    *kappa = NAN;
    303301
    304302    // XXX need to do this more carefully
    … …  
    306304        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
    307305        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
    308         *sigma1 = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
    309         *sigma2 = 0.25 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
     306        *sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
    310307    }
    311308    if (modelPSF->type == modelType_PS1_V1) {
    … …  
    323320    if (modelPSF->type == modelType_PS1_V1) {
    324321        psImageSmoothCacheKernel_PS1_V1 (smdata, *sigma, *kappa);
     322    }
     323
     324    return smdata;
     325}
     326
     327psImageSmooth2dCacheData *psImageSmooth2dCacheSetKernel (float *sigma, float *kappa, float nsigma, psImage *flux, pmModel *modelPSF) {
     328
     329    psAssert (modelPSF, "psf model must be defined");
     330   
     331    psEllipseAxes axes;
     332    bool useReff = pmModelUseReff (modelPSF->type);
     333    psF32 *PAR = modelPSF->params->data.F32;
     334    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
     335   
     336    *sigma = NAN;
     337    *kappa = NAN;
     338
     339    // XXX need to do this more carefully
     340    if (modelPSF->type == modelType_GAUSS) {
     341        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
     342        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
     343        *sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
     344    }
     345    if (modelPSF->type == modelType_PS1_V1) {
     346        *sigma = 0.5 * (axes.major + axes.minor);
     347        *kappa = PAR[PM_PAR_7];
     348    }
     349    psAssert (isfinite(*sigma), "invalid model type");
     350
     351    // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
     352    psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc (nsigma);
     353
     354    if (modelPSF->type == modelType_GAUSS) {
     355        // psImageSmooth2dCacheKernel_Gauss (smdata, *sigma);
     356    }
     357    if (modelPSF->type == modelType_PS1_V1) {
     358        psImageSmooth2dCacheKernel_PS1_V1 (smdata, *sigma, *kappa);
    325359    }
    326360
    … …  
    384418
    385419    pcm->use1Dgauss = true;
    386     pcm->smdata = psImageSmoothCacheSetKernel (&pcm->sigma, &pcm->kappa, source->pixels, source->modelPSF);
     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    }
    387425
    388426# else
    … …  
    467505
    468506        // If we have changed the window, we need to redefine the smoothing target vectors (but pcm->sigma,kappa,nsigma remain)
    469         psFree(pcm->smdata);
    470         pcm->smdata = psImageSmoothCacheAlloc (source->pixels, pcm->sigma, pcm->nsigma);
    471 
    472         pmModel *modelPSF = source->modelPSF;
    473         if (modelPSF->type == modelType_GAUSS) {
    474             psImageSmoothCacheKernel_Gauss (pcm->smdata, pcm->sigma);
    475         }
    476         if (modelPSF->type == modelType_PS1_V1) {
    477             psImageSmoothCacheKernel_PS1_V1 (pcm->smdata, pcm->sigma, pcm->kappa);
     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            }
    478529        }
    479530    }
    … …  
    503554        float sigma = NAN;
    504555        float kappa = NAN;
    505         psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, source->modelFlux, source->modelPSF);
    506 
    507         psImageSmoothCache_F32 (source->modelFlux, smdata);
    508         psFree (smdata);
    509 
     556
     557        if (USE_1D_CACHE) {
     558            psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, KERNEL_NSIGMA, source->modelFlux, source->modelPSF);
     559            psImageSmoothCache_F32 (source->modelFlux, smdata);
     560            psFree (smdata);
     561        } else {
     562            psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheSetKernel (&sigma, &kappa, KERNEL_NSIGMA, source->modelFlux, source->modelPSF);
     563            psImageSmooth2dCache_F32 (source->modelFlux, smdata);
     564            psFree (smdata);
     565        }
    510566        // old call: psImageSmooth (source->modelFlux, sigma, nsigma);
    511567    } else {
    … …  
    544600        float sigma = NAN;
    545601        float kappa = NAN;
    546         psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, source->modelFlux, source->modelPSF);
    547 
    548         psImageSmoothCache_F32 (source->modelFlux, smdata);
    549         psFree (smdata);
    550 
     602
     603        if (USE_1D_CACHE) {
     604            psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, KERNEL_NSIGMA, source->modelFlux, source->modelPSF);
     605            psImageSmoothCache_F32 (source->modelFlux, smdata);
     606            psFree (smdata);
     607        } else {
     608            psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheSetKernel (&sigma, &kappa, KERNEL_NSIGMA, source->modelFlux, source->modelPSF);
     609            psImageSmooth2dCache_F32 (source->modelFlux, smdata);
     610            psFree (smdata);
     611        }
    551612        // old call: psImageSmooth (source->modelFlux, sigma, nsigma);
    552613    } else {

  • branches/eam_branches/ipp-20130904/psModules/src/objects/pmPCMdata.h

    r36276 r36281  
    1414/// @addtogroup Objects Object Detection / Analysis Functions
    1515/// @{
     16
     17// XXX this is basically for testing -- when I am happy with the convolution process, I'll strip this out
     18# define USE_1D_CACHE 0
     19# define USE_1D_GAUSS 1
    1620
    1721/** pmPCMdata : PSF Convolved Model data storage structure
    … …  
    4044    bool use1Dgauss;
    4145    float kappa;
     46    float sigma;
    4247    float nsigma;
    4348

  • branches/eam_branches/ipp-20130904/psModules/src/objects/pmSourceIO.c

    r36270 r36281  
    646646
    647647      case PM_FPA_FILE_CFF: {
     648        // determine the output table format
     649        psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSPHOT");
     650        if (!status) {
     651            psError(PS_ERR_UNKNOWN, true, "missing recipe PSPHOT in config data");
     652            return false;
     653        }
     654
    648655        hdu = pmFPAviewThisHDU (view, fpa);
    649656        pmConfigConformHeader(hdu->header, file->format);
Note: See TracChangeset for help on using the changeset viewer.