IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 5, 2009, 12:36:19 PM (17 years ago)
Author:
Paul Price
Message:

Adding function to return the kernel used for smoothing (so we can use it to calculate the covariance matrix). Consolidated code to generate a normalised Gaussian vector.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/imageops/psImageConvolve.c

    r21183 r21335  
    77/// @author Eugene Magnier, IfA
    88///
    9 /// @version $Revision: 1.82 $ $Name: not supported by cvs2svn $
    10 /// @date $Date: 2009-01-27 06:39:37 $
     9/// @version $Revision: 1.83 $ $Name: not supported by cvs2svn $
     10/// @date $Date: 2009-02-05 22:36:19 $
    1111///
    1212/// Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    131131
    132132
    133 psKernel *psKernelGenerate(const psVector *tShifts,
    134                            const psVector *xShifts,
    135                            const psVector *yShifts,
    136                            float totalTime,
    137                            bool xyRelative)
     133psKernel *psKernelGenerate(const psVector *tShifts, const psVector *xShifts, const psVector *yShifts,
     134                           float totalTime, bool xyRelative)
    138135{
    139136    PS_ASSERT_VECTOR_NON_NULL(tShifts, NULL);
    … …  
    232229}
    233230
    234 psImage *psImageConvolveDirect(psImage *out,
    235                                const psImage *in,
    236                                const psKernel *kernel)
     231psImage *psImageConvolveDirect(psImage *out, const psImage *in, const psKernel *kernel)
    237232{
    238233    PS_ASSERT_IMAGE_NON_NULL(in, NULL);
    … …  
    466461}
    467462
    468 bool psImageSmooth(psImage *image,
    469                    double  sigma,
    470                    double  Nsigma)
     463
     464// Generate normalised Gaussian vector
     465#define IMAGE_SMOOTH_GAUSS(TARGET, SIZE, SIGMA, TYPE) \
     466    psVector *TARGET = psVectorAlloc(2 * SIZE + 1, PS_TYPE_##TYPE); /* Gaussian */ \
     467    double sum = 0.0; /* Sum of Gaussian, for normalisation */ \
     468    double factor = -0.5/PS_SQR(SIGMA); /* Multiplier for exponential */ \
     469    for (int i = -SIZE, j = 0; i <= SIZE; i++, j++) { \
     470        sum += TARGET->data.TYPE[j] = exp(factor * PS_SQR(i)); \
     471    } \
     472    for (int i = 0; i <= 2 * SIZE + 1; i++) { \
     473        TARGET->data.TYPE[i] /= sum; \
     474    }
     475
     476psKernel *psImageSmoothKernel(float sigma, float nSigma)
     477{
     478    PS_ASSERT_FLOAT_LARGER_THAN(sigma, 0.0, NULL);
     479    PS_ASSERT_FLOAT_LARGER_THAN(nSigma, 0.0, NULL);
     480
     481    int size = sigma * nSigma + 0.5; // Half-size of kernel
     482    psKernel *kernel = psKernelAlloc(-size, size, -size, size); // Kernel to return
     483
     484    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
     485    psF32 *gauss = &gaussNorm->data.F32[size]; // Convenient kernel-like indexing for Gaussian
     486    for (int y = -size; y <= size; y++) {
     487        for (int x = -size; x <= size; x++) {
     488            kernel->kernel[y][x] = gauss[y] * gauss[x];
     489        }
     490    }
     491    psFree(gaussNorm);
     492
     493    return kernel;
     494}
     495
     496
     497bool psImageSmooth(psImage *image, double  sigma, double  Nsigma)
    471498{
    472499    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    … …  
    474501    // relevant terms
    475502    int Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side for convolution kernel
    476     int Npixel = 2*Nrange + 1;          // Total number of pixels in convolution kernel
    477503    int Nx = image->numCols;            // Number of columns
    478504    int Ny = image->numRows;            // Number of rows
    479505
    480     #define IMAGESMOOTH_CASE(TYPE) \
    481 case PS_TYPE_##TYPE: { \
     506#define IMAGESMOOTH_CASE(TYPE) \
     507  case PS_TYPE_##TYPE: { \
    482508        /* generate normalized gaussian */ \
    483         psVector *gaussnorm = psVectorAlloc(Npixel, PS_TYPE_##TYPE); \
    484         { \
    485             double sum = 0.0; \
    486             double factor = -0.5/(sigma*sigma); \
    487             for (int i = -Nrange; i < Nrange + 1; i++) { \
    488                 gaussnorm->data.TYPE[i+Nrange] = exp(factor*i*i); \
    489                 sum += gaussnorm->data.TYPE[i+Nrange]; \
    490             } \
    491             for (int i = -Nrange; i < Nrange + 1; i++) { \
    492                 gaussnorm->data.TYPE[i+Nrange] /= sum; \
    493             } \
    494         } \
     509        IMAGE_SMOOTH_GAUSS(gaussnorm, Nrange, sigma, TYPE) \
    495510        ps##TYPE *gauss = &gaussnorm->data.TYPE[Nrange]; \
    496511        \
    … …  
    711726}
    712727
    713 psImage *psImageSmoothMask(psImage *output,
    714                            const psImage *image,
    715                            const psImage *mask,
    716                            psImageMaskType maskVal,
    717                            float sigma,
    718                            float numSigma,
    719                            float minGauss)
     728psImage *psImageSmoothMask(psImage *output, const psImage *image, const psImage *mask,
     729                           psImageMaskType maskVal, float sigma, float numSigma, float minGauss)
    720730{
    721731    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    … …  
    726736    int numCols = image->numCols, numRows = image->numRows; // Size of image
    727737    int xLast = numCols - 1, yLast = numRows - 1; // Last index
    728 
    729     // relevant terms
    730738    int size = sigma * numSigma + 0.5;  // Half-size of kernel
    731     int fullSize = 2*size + 1;          // Total number of pixels in convolution kernel
    732739
    733740    // Generate normalized gaussian
    734     psVector *gaussNorm = psVectorAlloc(fullSize + 1, PS_TYPE_F32); // Normalised Gaussian
    735     {
    736         double sum = 0.0;
    737         double norm = -0.5/(sigma*sigma);
    738         for (int i = 0, u = -size; i <= fullSize; i++, u++) {
    739             float value = expf(norm*PS_SQR(u));
    740             gaussNorm->data.F32[i] = value;
    741             sum += value;
    742         }
    743         sum = 1.0 / sum;
    744         for (int i = 0; i <= fullSize; i++) {
    745             gaussNorm->data.F32[i] *= sum;
    746         }
    747     }
     741    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
    748742    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
    749743
    … …  
    828822}
    829823
    830 bool psImageSmoothMask_ScanRows_F32 (psImage *calculation,
    831                                      psImage *calcMask,
    832                                      const psImage *image,
    833                                      const psImage *mask,
    834                                      psImageMaskType maskVal,
    835                                      psVector *gaussNorm,
    836                                      float minGauss,
    837                                      int size,
    838                                      int rowStart,
    839                                      int rowStop)
     824bool psImageSmoothMask_ScanRows_F32(psImage *calculation, psImage *calcMask, const psImage *image,
     825                                    const psImage *mask, psImageMaskType maskVal, psVector *gaussNorm,
     826                                    float minGauss, int size, int rowStart, int rowStop)
    840827{
    841828    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
    … …  
    873860}
    874861
    875 bool psImageSmoothMask_ScanCols_F32 (psImage *output,
    876                                      psImage *calculation,
    877                                      psImage *calcMask,
    878                                      psImageMaskType maskVal,
    879                                      psVector *gaussNorm,
    880                                      float minGauss,
    881                                      int size,
    882                                      int colStart,
    883                                      int colStop)
     862bool psImageSmoothMask_ScanCols_F32(psImage *output, psImage *calculation, psImage *calcMask,
     863                                    psImageMaskType maskVal, psVector *gaussNorm, float minGauss,
     864                                    int size, int colStart, int colStop)
    884865{
    885866    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
    … …  
    953934}
    954935
    955 psImage *psImageSmoothMask_Threaded(psImage *output,
    956                                     const psImage *image,
    957                                     const psImage *mask,
    958                                     psImageMaskType maskVal,
    959                                     float sigma,
    960                                     float numSigma,
    961                                     float minGauss)
     936psImage *psImageSmoothMask_Threaded(psImage *output, const psImage *image, const psImage *mask,
     937                                    psImageMaskType maskVal, float sigma, float numSigma, float minGauss)
    962938{
    963939    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    … …  
    971947    int scanCols = nThreads ? (numCols / 4 / nThreads) : numCols;
    972948    int scanRows = nThreads ? (numRows / 4 / nThreads) : numRows;
    973 
    974     // relevant terms
    975949    int size = sigma * numSigma + 0.5;  // Half-size of kernel
    976     int fullSize = 2*size + 1;          // Total number of pixels in convolution kernel
    977950
    978951    // Generate normalized gaussian
    979     psVector *gaussNorm = psVectorAlloc(fullSize + 1, PS_TYPE_F32); // Normalised Gaussian
    980     {
    981         double sum = 0.0;
    982         double norm = -0.5/(sigma*sigma);
    983         for (int i = 0, u = -size; i <= fullSize; i++, u++) {
    984             float value = expf(norm*PS_SQR(u));
    985             gaussNorm->data.F32[i] = value;
    986             sum += value;
    987         }
    988         sum = 1.0 / sum;
    989         for (int i = 0; i <= fullSize; i++) {
    990             gaussNorm->data.F32[i] *= sum;
    991         }
    992     }
     952    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
    993953
    994954    switch (image->type.type) {
    … …  
    10891049
    10901050
    1091 bool psImageSmoothMaskF32 (psImage *image,
    1092                            psImage *mask,
    1093                            psImageMaskType maskVal,
    1094                            double  sigma,
    1095                            double  Nsigma)
     1051bool psImageSmoothMaskF32(psImage *image, psImage *mask, psImageMaskType maskVal,
     1052                          double  sigma, double  Nsigma)
    10961053{
    10971054    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    … …  
    11001057    // relevant terms
    11011058    int Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side for convolution kernel
    1102     int Npixel = 2*Nrange + 1;          // Total number of pixels in convolution kernel
    11031059    int Nx = image->numCols;            // Number of columns
    11041060    int Ny = image->numRows;            // Number of rows
    11051061
    11061062    /* generate normalized gaussian */
    1107     psVector *gaussnorm = psVectorAlloc(Npixel, PS_TYPE_F32);
    1108     {
    1109         double sum = 0.0;
    1110         double factor = -0.5/(sigma*sigma);
    1111         for (int i = -Nrange; i < Nrange + 1; i++) {
    1112             gaussnorm->data.F32[i+Nrange] = exp(factor*i*i);
    1113             sum += gaussnorm->data.F32[i+Nrange];
    1114         }
    1115         for (int i = -Nrange; i < Nrange + 1; i++) {
    1116             gaussnorm->data.F32[i+Nrange] /= sum;
    1117         }
    1118     }
     1063    IMAGE_SMOOTH_GAUSS(gaussnorm, Nrange, sigma, F32);
    11191064    psF32 *gauss = &gaussnorm->data.F32[Nrange];
    11201065
    … …  
    11481093    psFree(calculation);
    11491094
    1150     // XXX test
    1151     if (0) {
    1152         psFree(gaussnorm);
    1153         return true;
    1154     }
    1155 
    11561095    /** Smooth in Y direction  **/
    11571096    // allocate and save Nrange extra row vectors for storage. these will be
Note: See TracChangeset for help on using the changeset viewer.