IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1863


Ignore:
Timestamp:
Sep 23, 2004, 8:30:57 AM (22 years ago)
Author:
desonia
Message:

added initial/untested version of psImageConvolution.

Location:
trunk/psLib/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r1840 r1863  
    55 *  @author Robert DeSonia, MHPCC
    66 *
    7  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2004-09-21 22:30:19 $
     7 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2004-09-23 18:30:57 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1111 */
    1212
     13#include <string.h>
     14
    1315#include "psImageConvolve.h"
    1416#include "psImageFFT.h"
     17#include "psMatrixVectorArithmetic.h"
    1518#include "psMemory.h"
    1619#include "psLogMsg.h"
    … …  
    2730    int numCols;
    2831
     32    // following is explicitly spelled out in the SDRS as a requirement
    2933    if (yMin > yMax) {
    3034        int temp = yMin;
    … …  
    3741    }
    3842
     43    // following is explicitly spelled out in the SDRS as a requirement
    3944    if (xMin > xMax) {
    4045        int temp = xMin;
    … …  
    5560    result->yMin = yMin;
    5661    result->yMax = yMax;
    57     result->p_data = psImageAlloc(numCols,numRows,PS_TYPE_F32);
    58     result->p_kernelRows = psAlloc(sizeof(psF32*)*numRows);
    59 
    60     psF32** kernelRows = result->p_kernelRows;
    61     psF32** imageRows = result->p_data->data.F32;
     62    result->image = psImageAlloc(numCols,numRows,PS_TYPE_KERNEL);
     63    result->p_kernelRows = psAlloc(sizeof(psKernelType*)*numRows);
     64
     65    psKernelType** kernelRows = result->p_kernelRows;
     66    psKernelType** imageRows = result->image->data.PS_TYPE_KERNEL_DATA;
    6267    for (int i = 0; i < numRows; i++) {
    6368        kernelRows[i] = imageRows[i] - xMin;
    … …  
    7378{
    7479    if (ptr != NULL) {
    75         psFree(ptr->p_data);
     80        psFree(ptr->image);
    7681        psFree(ptr->p_kernelRows);
    7782    }
    7883}
    7984
    80 psKernel* psKernelGenerate(const psVector* xShifts, const psVector* yShifts)
     85psKernel* psKernelGenerate(const psVector* tShifts,
     86                           const psVector* xShifts,
     87                           const psVector* yShifts,
     88                           bool relative)
    8189{
    8290    int x = 0;
    8391    int y = 0;
     92    int t = 0;
     93    int deltaT = 0;
    8494    int xMin = 0;
    8595    int xMax = 0;
    … …  
    8797    int yMax = 0;
    8898    int length = 0;
    89     psF32 relativeTime = 1.0f;            // fraction of total time for each shift clock
     99    psKernelType normalizeTime = 1.0;  // fraction of total time for each shift clock
    90100    psKernel* result = NULL;
    91     psF32** kernel = NULL;
     101    psKernelType** kernel = NULL;
    92102
    93103    // got non-NULL vectors?
    … …  
    100110
    101111    // types match?
    102     if (xShifts->type.type != yShifts->type.type) {
     112    if (xShifts->type.type != yShifts->type.type ||
     113            tShifts->type.type != xShifts->type.type) {
    103114        char* typeXStr;
    104115        char* typeYStr;
     116        char* typeTStr;
    105117        PS_TYPE_NAME(typeXStr,xShifts->type.type);
    106118        PS_TYPE_NAME(typeYStr,yShifts->type.type);
     119        PS_TYPE_NAME(typeTStr,tShifts->type.type);
    107120        psErrorMsg(PS_ERRORNAME_DOMAIN "psKernelGenerate",
    108121                   PS_ERR_BAD_PARAMETER_TYPE, true,
    109122                   PS_ERRORTEXT_psImageConvolve_SHIFT_TYPE_MISMATCH,
    110                    typeXStr, typeYStr);
     123                   typeTStr, typeXStr, typeYStr);
    111124        return NULL;
    112125    }
    113126
    114127    // determine the usable length of the shift vector
    115     if (xShifts->n != yShifts->n) {
     128    length = xShifts->n;
     129    if (length != yShifts->n ||
     130            length != tShifts->n) {
    116131        psLogMsg(__func__,PS_LOG_WARN,"Shift vectors found to be different sizes.");
    117         length = xShifts->n;
    118132        if (yShifts->n < length) {
    119133            length = yShifts->n;
    120134        }
    121     } else {
    122         length = xShifts->n;
    123     }
    124     relativeTime = 1 / (1+length); // 1+length 'cause we count origin before shift
    125     // e.g., 1 time unit at origin + a time unit for each shift clock
     135        if (tShifts->n < length) {
     136            length = tShifts->n;
     137        }
     138    }
    126139
    127140    #define KERNEL_GENERATE_CASE(TYPE) \
    128141case PS_TYPE_##TYPE: { \
     142        ps##TYPE *tShiftData = tShifts->data.TYPE; \
    129143        ps##TYPE *xShiftData = xShifts->data.TYPE; \
    130144        ps##TYPE *yShiftData = yShifts->data.TYPE; \
    131145        for (int lcv = 0; lcv < length; lcv++) { \
    132             x += xShiftData[lcv]; \
    133             y += yShiftData[lcv]; \
     146            if (relative) { \
     147                x += xShiftData[lcv]; \
     148                y += yShiftData[lcv]; \
     149            } else { \
     150                x = xShiftData[lcv]; \
     151                y = yShiftData[lcv]; \
     152            } \
    134153            if (x < xMin) { \
    135154                xMin = x; \
    … …  
    146165        kernel = result->kernel; \
    147166        \
    148         kernel[0][0]+=relativeTime;    /* mark the original point as one exposure. */ \
    149         \
    150         x = 0; \
    151         y = 0; \
    152         for (int i = 0; i < length; i++) { \
    153             x += xShiftData[i]; \
    154             y += yShiftData[i]; \
    155             kernel[y][x]+=relativeTime; \
     167        x = xShiftData[0]; \
     168        y = yShiftData[0]; \
     169        t = tShiftData[0]; \
     170        normalizeTime = 1.0 / (psKernelType)(tShiftData[length-1]-t); \
     171        for (int i = 1; i < length; i++) { \
     172            deltaT = t; \
     173            if (relative) { \
     174                t += tShiftData[i]; \
     175                x += xShiftData[i]; \
     176                y += yShiftData[i]; \
     177            } else { \
     178                t = tShiftData[i]; \
     179                x = xShiftData[i]; \
     180                y = yShiftData[i]; \
     181            } \
     182            deltaT = t - deltaT; \
     183            \
     184            kernel[y][x] += deltaT * normalizeTime; \
    156185        } \
    157186    }
    … …  
    186215psImage* psImageConvolve(psImage* out, const psImage* in, const psKernel* kernel, bool direct)
    187216{
    188 
    189 
    190     return NULL;
     217    if (in == NULL) {
     218        psFree(out);
     219        return NULL;
     220    }
     221
     222    if (kernel == NULL) {
     223        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
     224                   PS_ERR_BAD_PARAMETER_NULL, true,
     225                   PS_ERRORTEXT_psImageConvolve_KERNEL_NULL);
     226        psFree(out);
     227        return NULL;
     228    }
     229    int xMin = kernel->xMin;
     230    int xMax = kernel->xMax;
     231    int yMin = kernel->yMin;
     232    int yMax = kernel->yMax;
     233    psKernelType** kData = kernel->kernel;
     234
     235    // make the output image to the proper size and type
     236    int numRows = in->numRows;
     237    int numCols = in->numCols;
     238
     239
     240
     241    if (direct) {
     242        // spatial convolution
     243
     244        #define SPATIAL_CONVOLVE_CASE(TYPE) \
     245    case PS_TYPE_##TYPE: { \
     246            ps##TYPE** inData = in->data.TYPE; \
     247            out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \
     248            for (int row=0;row<numRows;row++) { \
     249                ps##TYPE* outRow = out->data.TYPE[row]; \
     250                for (int col=0;col<numCols;col++) { \
     251                    ps##TYPE pixel = 0.0; \
     252                    for (int kRow = yMin; kRow < yMax; kRow++) { \
     253                        if (row-kRow >= 0 && row-kRow < numRows) { \
     254                            for (int kCol = xMin; kCol < xMax; kCol++) { \
     255                                if (col-kCol >= 0 && col-kCol < numCols) { \
     256                                    pixel += kData[kRow][kCol] * inData[row-kRow][col-kCol]; \
     257                                } \
     258                            } \
     259                        } \
     260                    } \
     261                    outRow[col] = pixel; \
     262                } \
     263            } \
     264        } \
     265        break;
     266
     267        switch (in->type.type) {
     268            SPATIAL_CONVOLVE_CASE(U8)
     269            SPATIAL_CONVOLVE_CASE(U16)
     270            SPATIAL_CONVOLVE_CASE(U32)
     271            SPATIAL_CONVOLVE_CASE(U64)
     272            SPATIAL_CONVOLVE_CASE(S8)
     273            SPATIAL_CONVOLVE_CASE(S16)
     274            SPATIAL_CONVOLVE_CASE(S32)
     275            SPATIAL_CONVOLVE_CASE(S64)
     276            SPATIAL_CONVOLVE_CASE(F32)
     277            SPATIAL_CONVOLVE_CASE(F64)
     278            SPATIAL_CONVOLVE_CASE(C32)
     279            SPATIAL_CONVOLVE_CASE(C64)
     280
     281        default: {
     282                char* typeStr;
     283                PS_TYPE_NAME(typeStr,in->type.type);
     284                psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
     285                           PS_ERR_BAD_PARAMETER_TYPE, true,
     286                           PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED,
     287                           typeStr);
     288                psFree(out);
     289                return NULL;
     290
     291            }
     292        }
     293
     294
     295    } else {
     296        // fourier convolution
     297        int kRows = kernel->image->numRows;
     298        int kCols = kernel->image->numCols;
     299        int x0 = (numCols - kCols) / 2;
     300        int y0 = (numRows - kRows) / 2;
     301        psImage* paddedKernel;
     302
     303        // check to see if kernel is smaller, otherwise padding it up will fail.
     304        if (x0 < 0 || y0 < 0) {
     305            psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
     306                       PS_ERR_BAD_PARAMETER_SIZE, true,
     307                       PS_ERRORTEXT_psImageConvolve_KERNEL_TOO_LARGE,
     308                       kernel->image->numCols,kernel->image->numRows,
     309                       numCols, numRows);
     310            psFree(out);
     311            return NULL;
     312        }
     313
     314        // pad the kernel to the same size of image
     315        paddedKernel = psImageAlloc(numCols,numRows,PS_TYPE_KERNEL);
     316        memset(paddedKernel->data.V,0,sizeof(psKernelType)*numCols*numRows); // zeroize image
     317        for (int row=0;row<kRows;row++) {
     318            psKernelType* padRow = paddedKernel->data.PS_TYPE_KERNEL_DATA[row+y0];
     319            psKernelType* kernelRow = kernel->image->data.PS_TYPE_KERNEL_DATA[row];
     320            for (int col=0;col<kCols;col++) {
     321                padRow[col+x0] = kernelRow[col];
     322            }
     323        }
     324
     325        psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD);
     326        if (kernelFourier == NULL) {
     327            psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
     328                       PS_ERR_UNKNOWN, false,
     329                       PS_ERRORTEXT_psImageConvolve_KERNEL_FFT_FAILED);
     330            psFree(out);
     331            return NULL;
     332        }
     333
     334
     335        psImage* inFourier = psImageFFT(NULL, in, PS_FFT_FORWARD);
     336        if (inFourier == NULL) {
     337            psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
     338                       PS_ERR_UNKNOWN, false,
     339                       PS_ERRORTEXT_psImageConvolve_FFT_FAILED);
     340            psFree(out);
     341            return NULL;
     342        }
     343
     344        // convolution in fourier domain is just a pixel-wise multiplication
     345        psImage* outFourier = psImageAlloc(numCols,numRows,inFourier->type.type);
     346        psBinaryOp(outFourier,inFourier,"*",kernelFourier);
     347
     348        psImage* complexOut = psImageFFT(NULL, outFourier, PS_FFT_REVERSE);
     349        if (complexOut == NULL) {
     350            psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
     351                       PS_ERR_UNKNOWN, false,
     352                       PS_ERRORTEXT_psImageConvolve_FFT_FAILED);
     353            psFree(out);
     354            return NULL;
     355        }
     356
     357        // since our image was real to start, the result is real
     358        out = psImageReal(out,complexOut);
     359
     360        psFree(complexOut);
     361        psFree(kernelFourier);
     362        psFree(inFourier);
     363        psFree(paddedKernel);
     364
     365    }
     366
     367    return out;
    191368}

  • trunk/psLib/src/image/psImageConvolve.h

    r1653 r1863  
    77 *  @author Robert DeSonia, MHPCC
    88 *
    9  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-08-28 01:18:28 $
     9 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-23 18:30:57 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1313 */
    1414
    15 #ifndef PS_IMAGE_FFT_H
    16 #define PS_IMAGE_FFT_H
     15#ifndef PS_IMAGE_CONVOLVE_H
     16#define PS_IMAGE_CONVOLVE_H
    1717
    1818#include<stdbool.h>
    … …  
    2020#include "psImage.h"
    2121#include "psVector.h"
     22#include "psType.h"
     23
     24#define PS_TYPE_KERNEL PS_TYPE_F32     /**< the data member to use for kernel image */
     25#define PS_TYPE_KERNEL_DATA F32        /**< the data member to use for kernel image */
     26#define PS_TYPE_KERNEL_NAME "psF32"    /**< the data type for kernel as a string */
     27
     28typedef psF32 psKernelType;
    2229
    2330/** A convolution kernel */
    2431typedef struct
    2532{
     33    psImage* image;                    ///< Kernel data, in the form of an image
    2634    int xMin;                          ///< Most negative x index
    2735    int yMin;                          ///< Most negative y index
    2836    int xMax;                          ///< Most positive x index
    2937    int yMax;                          ///< Most positive y index
    30     psF32** kernel;                    ///< Pointer to the kernel data
    31     psF32** p_kernelRows;              ///< Pointer to the rows of the kernel data
    32     psImage* p_data;                   ///< Kernel data, in the form of an image
     38    psKernelType** kernel;             ///< Pointer to the kernel data
     39    psKernelType** p_kernelRows;       ///< Pointer to the rows of the kernel data; not intended for user use.
    3340}
    3441psKernel;
    … …  
    6269psKernel* psKernelAlloc(
    6370    int xMin,                          ///< Most negative x index
    64     int xMax,                          ///< Most negative y index
    65     int yMin,                          ///< Most positive x index
     71    int xMax,                          ///< Most positive x index
     72    int yMin,                          ///< Most negative y index
    6673    int yMax                           ///< Most positive y index
    6774);
    … …  
    8491 */
    8592psKernel* psKernelGenerate(
    86     const psVector* xShifts,           ///< list of shifts relative to a reference point
    87     const psVector* yShifts            ///< list of shifts relative to a reference point
     93    const psVector* tShifts,           ///< list of time shifts
     94    const psVector* xShifts,           ///< list of x-axis shifts
     95    const psVector* yShifts,           ///< list of y-axis shifts
     96    bool relative
    8897);
    8998

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

    r1840 r1863  
    55 *  @author Robert DeSonia, MHPCC
    66 *
    7  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2004-09-21 22:30:19 $
     7 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2004-09-23 18:30:57 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1111 */
    1212
     13#include <string.h>
     14
    1315#include "psImageConvolve.h"
    1416#include "psImageFFT.h"
     17#include "psMatrixVectorArithmetic.h"
    1518#include "psMemory.h"
    1619#include "psLogMsg.h"
    … …  
    2730    int numCols;
    2831
     32    // following is explicitly spelled out in the SDRS as a requirement
    2933    if (yMin > yMax) {
    3034        int temp = yMin;
    … …  
    3741    }
    3842
     43    // following is explicitly spelled out in the SDRS as a requirement
    3944    if (xMin > xMax) {
    4045        int temp = xMin;
    … …  
    5560    result->yMin = yMin;
    5661    result->yMax = yMax;
    57     result->p_data = psImageAlloc(numCols,numRows,PS_TYPE_F32);
    58     result->p_kernelRows = psAlloc(sizeof(psF32*)*numRows);
    59 
    60     psF32** kernelRows = result->p_kernelRows;
    61     psF32** imageRows = result->p_data->data.F32;
     62    result->image = psImageAlloc(numCols,numRows,PS_TYPE_KERNEL);
     63    result->p_kernelRows = psAlloc(sizeof(psKernelType*)*numRows);
     64
     65    psKernelType** kernelRows = result->p_kernelRows;
     66    psKernelType** imageRows = result->image->data.PS_TYPE_KERNEL_DATA;
    6267    for (int i = 0; i < numRows; i++) {
    6368        kernelRows[i] = imageRows[i] - xMin;
    … …  
    7378{
    7479    if (ptr != NULL) {
    75         psFree(ptr->p_data);
     80        psFree(ptr->image);
    7681        psFree(ptr->p_kernelRows);
    7782    }
    7883}
    7984
    80 psKernel* psKernelGenerate(const psVector* xShifts, const psVector* yShifts)
     85psKernel* psKernelGenerate(const psVector* tShifts,
     86                           const psVector* xShifts,
     87                           const psVector* yShifts,
     88                           bool relative)
    8189{
    8290    int x = 0;
    8391    int y = 0;
     92    int t = 0;
     93    int deltaT = 0;
    8494    int xMin = 0;
    8595    int xMax = 0;
    … …  
    8797    int yMax = 0;
    8898    int length = 0;
    89     psF32 relativeTime = 1.0f;            // fraction of total time for each shift clock
     99    psKernelType normalizeTime = 1.0;  // fraction of total time for each shift clock
    90100    psKernel* result = NULL;
    91     psF32** kernel = NULL;
     101    psKernelType** kernel = NULL;
    92102
    93103    // got non-NULL vectors?
    … …  
    100110
    101111    // types match?
    102     if (xShifts->type.type != yShifts->type.type) {
     112    if (xShifts->type.type != yShifts->type.type ||
     113            tShifts->type.type != xShifts->type.type) {
    103114        char* typeXStr;
    104115        char* typeYStr;
     116        char* typeTStr;
    105117        PS_TYPE_NAME(typeXStr,xShifts->type.type);
    106118        PS_TYPE_NAME(typeYStr,yShifts->type.type);
     119        PS_TYPE_NAME(typeTStr,tShifts->type.type);
    107120        psErrorMsg(PS_ERRORNAME_DOMAIN "psKernelGenerate",
    108121                   PS_ERR_BAD_PARAMETER_TYPE, true,
    109122                   PS_ERRORTEXT_psImageConvolve_SHIFT_TYPE_MISMATCH,
    110                    typeXStr, typeYStr);
     123                   typeTStr, typeXStr, typeYStr);
    111124        return NULL;
    112125    }
    113126
    114127    // determine the usable length of the shift vector
    115     if (xShifts->n != yShifts->n) {
     128    length = xShifts->n;
     129    if (length != yShifts->n ||
     130            length != tShifts->n) {
    116131        psLogMsg(__func__,PS_LOG_WARN,"Shift vectors found to be different sizes.");
    117         length = xShifts->n;
    118132        if (yShifts->n < length) {
    119133            length = yShifts->n;
    120134        }
    121     } else {
    122         length = xShifts->n;
    123     }
    124     relativeTime = 1 / (1+length); // 1+length 'cause we count origin before shift
    125     // e.g., 1 time unit at origin + a time unit for each shift clock
     135        if (tShifts->n < length) {
     136            length = tShifts->n;
     137        }
     138    }
    126139
    127140    #define KERNEL_GENERATE_CASE(TYPE) \
    128141case PS_TYPE_##TYPE: { \
     142        ps##TYPE *tShiftData = tShifts->data.TYPE; \
    129143        ps##TYPE *xShiftData = xShifts->data.TYPE; \
    130144        ps##TYPE *yShiftData = yShifts->data.TYPE; \
    131145        for (int lcv = 0; lcv < length; lcv++) { \
    132             x += xShiftData[lcv]; \
    133             y += yShiftData[lcv]; \
     146            if (relative) { \
     147                x += xShiftData[lcv]; \
     148                y += yShiftData[lcv]; \
     149            } else { \
     150                x = xShiftData[lcv]; \
     151                y = yShiftData[lcv]; \
     152            } \
    134153            if (x < xMin) { \
    135154                xMin = x; \
    … …  
    146165        kernel = result->kernel; \
    147166        \
    148         kernel[0][0]+=relativeTime;    /* mark the original point as one exposure. */ \
    149         \
    150         x = 0; \
    151         y = 0; \
    152         for (int i = 0; i < length; i++) { \
    153             x += xShiftData[i]; \
    154             y += yShiftData[i]; \
    155             kernel[y][x]+=relativeTime; \
     167        x = xShiftData[0]; \
     168        y = yShiftData[0]; \
     169        t = tShiftData[0]; \
     170        normalizeTime = 1.0 / (psKernelType)(tShiftData[length-1]-t); \
     171        for (int i = 1; i < length; i++) { \
     172            deltaT = t; \
     173            if (relative) { \
     174                t += tShiftData[i]; \
     175                x += xShiftData[i]; \
     176                y += yShiftData[i]; \
     177            } else { \
     178                t = tShiftData[i]; \
     179                x = xShiftData[i]; \
     180                y = yShiftData[i]; \
     181            } \
     182            deltaT = t - deltaT; \
     183            \
     184            kernel[y][x] += deltaT * normalizeTime; \
    156185        } \
    157186    }
    … …  
    186215psImage* psImageConvolve(psImage* out, const psImage* in, const psKernel* kernel, bool direct)
    187216{
    188 
    189 
    190     return NULL;
     217    if (in == NULL) {
     218        psFree(out);
     219        return NULL;
     220    }
     221
     222    if (kernel == NULL) {
     223        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
     224                   PS_ERR_BAD_PARAMETER_NULL, true,
     225                   PS_ERRORTEXT_psImageConvolve_KERNEL_NULL);
     226        psFree(out);
     227        return NULL;
     228    }
     229    int xMin = kernel->xMin;
     230    int xMax = kernel->xMax;
     231    int yMin = kernel->yMin;
     232    int yMax = kernel->yMax;
     233    psKernelType** kData = kernel->kernel;
     234
     235    // make the output image to the proper size and type
     236    int numRows = in->numRows;
     237    int numCols = in->numCols;
     238
     239
     240
     241    if (direct) {
     242        // spatial convolution
     243
     244        #define SPATIAL_CONVOLVE_CASE(TYPE) \
     245    case PS_TYPE_##TYPE: { \
     246            ps##TYPE** inData = in->data.TYPE; \
     247            out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \
     248            for (int row=0;row<numRows;row++) { \
     249                ps##TYPE* outRow = out->data.TYPE[row]; \
     250                for (int col=0;col<numCols;col++) { \
     251                    ps##TYPE pixel = 0.0; \
     252                    for (int kRow = yMin; kRow < yMax; kRow++) { \
     253                        if (row-kRow >= 0 && row-kRow < numRows) { \
     254                            for (int kCol = xMin; kCol < xMax; kCol++) { \
     255                                if (col-kCol >= 0 && col-kCol < numCols) { \
     256                                    pixel += kData[kRow][kCol] * inData[row-kRow][col-kCol]; \
     257                                } \
     258                            } \
     259                        } \
     260                    } \
     261                    outRow[col] = pixel; \
     262                } \
     263            } \
     264        } \
     265        break;
     266
     267        switch (in->type.type) {
     268            SPATIAL_CONVOLVE_CASE(U8)
     269            SPATIAL_CONVOLVE_CASE(U16)
     270            SPATIAL_CONVOLVE_CASE(U32)
     271            SPATIAL_CONVOLVE_CASE(U64)
     272            SPATIAL_CONVOLVE_CASE(S8)
     273            SPATIAL_CONVOLVE_CASE(S16)
     274            SPATIAL_CONVOLVE_CASE(S32)
     275            SPATIAL_CONVOLVE_CASE(S64)
     276            SPATIAL_CONVOLVE_CASE(F32)
     277            SPATIAL_CONVOLVE_CASE(F64)
     278            SPATIAL_CONVOLVE_CASE(C32)
     279            SPATIAL_CONVOLVE_CASE(C64)
     280
     281        default: {
     282                char* typeStr;
     283                PS_TYPE_NAME(typeStr,in->type.type);
     284                psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
     285                           PS_ERR_BAD_PARAMETER_TYPE, true,
     286                           PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED,
     287                           typeStr);
     288                psFree(out);
     289                return NULL;
     290
     291            }
     292        }
     293
     294
     295    } else {
     296        // fourier convolution
     297        int kRows = kernel->image->numRows;
     298        int kCols = kernel->image->numCols;
     299        int x0 = (numCols - kCols) / 2;
     300        int y0 = (numRows - kRows) / 2;
     301        psImage* paddedKernel;
     302
     303        // check to see if kernel is smaller, otherwise padding it up will fail.
     304        if (x0 < 0 || y0 < 0) {
     305            psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
     306                       PS_ERR_BAD_PARAMETER_SIZE, true,
     307                       PS_ERRORTEXT_psImageConvolve_KERNEL_TOO_LARGE,
     308                       kernel->image->numCols,kernel->image->numRows,
     309                       numCols, numRows);
     310            psFree(out);
     311            return NULL;
     312        }
     313
     314        // pad the kernel to the same size of image
     315        paddedKernel = psImageAlloc(numCols,numRows,PS_TYPE_KERNEL);
     316        memset(paddedKernel->data.V,0,sizeof(psKernelType)*numCols*numRows); // zeroize image
     317        for (int row=0;row<kRows;row++) {
     318            psKernelType* padRow = paddedKernel->data.PS_TYPE_KERNEL_DATA[row+y0];
     319            psKernelType* kernelRow = kernel->image->data.PS_TYPE_KERNEL_DATA[row];
     320            for (int col=0;col<kCols;col++) {
     321                padRow[col+x0] = kernelRow[col];
     322            }
     323        }
     324
     325        psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD);
     326        if (kernelFourier == NULL) {
     327            psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
     328                       PS_ERR_UNKNOWN, false,
     329                       PS_ERRORTEXT_psImageConvolve_KERNEL_FFT_FAILED);
     330            psFree(out);
     331            return NULL;
     332        }
     333
     334
     335        psImage* inFourier = psImageFFT(NULL, in, PS_FFT_FORWARD);
     336        if (inFourier == NULL) {
     337            psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
     338                       PS_ERR_UNKNOWN, false,
     339                       PS_ERRORTEXT_psImageConvolve_FFT_FAILED);
     340            psFree(out);
     341            return NULL;
     342        }
     343
     344        // convolution in fourier domain is just a pixel-wise multiplication
     345        psImage* outFourier = psImageAlloc(numCols,numRows,inFourier->type.type);
     346        psBinaryOp(outFourier,inFourier,"*",kernelFourier);
     347
     348        psImage* complexOut = psImageFFT(NULL, outFourier, PS_FFT_REVERSE);
     349        if (complexOut == NULL) {
     350            psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
     351                       PS_ERR_UNKNOWN, false,
     352                       PS_ERRORTEXT_psImageConvolve_FFT_FAILED);
     353            psFree(out);
     354            return NULL;
     355        }
     356
     357        // since our image was real to start, the result is real
     358        out = psImageReal(out,complexOut);
     359
     360        psFree(complexOut);
     361        psFree(kernelFourier);
     362        psFree(inFourier);
     363        psFree(paddedKernel);
     364
     365    }
     366
     367    return out;
    191368}

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

    r1653 r1863  
    77 *  @author Robert DeSonia, MHPCC
    88 *
    9  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-08-28 01:18:28 $
     9 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-23 18:30:57 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1313 */
    1414
    15 #ifndef PS_IMAGE_FFT_H
    16 #define PS_IMAGE_FFT_H
     15#ifndef PS_IMAGE_CONVOLVE_H
     16#define PS_IMAGE_CONVOLVE_H
    1717
    1818#include<stdbool.h>
    … …  
    2020#include "psImage.h"
    2121#include "psVector.h"
     22#include "psType.h"
     23
     24#define PS_TYPE_KERNEL PS_TYPE_F32     /**< the data member to use for kernel image */
     25#define PS_TYPE_KERNEL_DATA F32        /**< the data member to use for kernel image */
     26#define PS_TYPE_KERNEL_NAME "psF32"    /**< the data type for kernel as a string */
     27
     28typedef psF32 psKernelType;
    2229
    2330/** A convolution kernel */
    2431typedef struct
    2532{
     33    psImage* image;                    ///< Kernel data, in the form of an image
    2634    int xMin;                          ///< Most negative x index
    2735    int yMin;                          ///< Most negative y index
    2836    int xMax;                          ///< Most positive x index
    2937    int yMax;                          ///< Most positive y index
    30     psF32** kernel;                    ///< Pointer to the kernel data
    31     psF32** p_kernelRows;              ///< Pointer to the rows of the kernel data
    32     psImage* p_data;                   ///< Kernel data, in the form of an image
     38    psKernelType** kernel;             ///< Pointer to the kernel data
     39    psKernelType** p_kernelRows;       ///< Pointer to the rows of the kernel data; not intended for user use.
    3340}
    3441psKernel;
    … …  
    6269psKernel* psKernelAlloc(
    6370    int xMin,                          ///< Most negative x index
    64     int xMax,                          ///< Most negative y index
    65     int yMin,                          ///< Most positive x index
     71    int xMax,                          ///< Most positive x index
     72    int yMin,                          ///< Most negative y index
    6673    int yMax                           ///< Most positive y index
    6774);
    … …  
    8491 */
    8592psKernel* psKernelGenerate(
    86     const psVector* xShifts,           ///< list of shifts relative to a reference point
    87     const psVector* yShifts            ///< list of shifts relative to a reference point
     93    const psVector* tShifts,           ///< list of time shifts
     94    const psVector* xShifts,           ///< list of x-axis shifts
     95    const psVector* yShifts,           ///< list of y-axis shifts
     96    bool relative
    8897);
    8998
Note: See TracChangeset for help on using the changeset viewer.