C语言图像处理卷积,图像处理之基础---卷积及其快速算法的C++实现

头文件：

/*

* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), [email protected]

*

* This program is free software; you can redistribute it and/or modify it

* under the terms of the GNU General Public License as published by the

* Free Software Foundation, either version 2 or any later version.

*

* Redistribution and use in source and binary forms, with or without

* modification, are permitted provided that the following conditions are met:

*

* 1. Redistributions of source code must retain the above copyright notice,

* this list of conditions and the following disclaimer.

*

* 2. Redistributions in binary form must reproduce the above copyright

* notice, this list of conditions and the following disclaimer in the

* documentation and/or other materials provided with the distribution.

*

* This program is distributed in the hope that it will be useful, but WITHOUT

* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for

* more details. A copy of the GNU General Public License is available at:

* http://www.fsf.org/licensing/licenses

*/

/*****************************************************************************

* convolution.h

*

* Linear convolution and polynomial multiplication.

*

* The convolution routine "conv" is implemented by it‘s definition in time

* domain. If the sequence to be convoluted are long, you should use the

* fast convolution algorithm "fastConv", which is implemented in frequency

* domain by usin FFT.

*

* Zhang Ming, 2010-01, Xi‘an Jiaotong University.

*****************************************************************************/

#ifndef CONVOLUTION_H

#define CONVOLUTION_H

#include

#include

#include

namespace splab

{

template Vector conv( const Vector&,

const Vector& );

template Vector convolution( const Vector&,

const Vector& );

template Vector fastConv( const Vector&,

const Vector& );

#include

}

// namespace splab

#endif

// CONVOLUTION_H

实现文件：

/*

* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), [email protected]

*

* This program is free software; you can redistribute it and/or modify it

* under the terms of the GNU General Public License as published by the

* Free Software Foundation, either version 2 or any later version.

*

* Redistribution and use in source and binary forms, with or without

* modification, are permitted provided that the following conditions are met:

*

* 1. Redistributions of source code must retain the above copyright notice,

* this list of conditions and the following disclaimer.

*

* 2. Redistributions in binary form must reproduce the above copyright

* notice, this list of conditions and the following disclaimer in the

* documentation and/or other materials provided with the distribution.

*

* This program is distributed in the hope that it will be useful, but WITHOUT

* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for

* more details. A copy of the GNU General Public License is available at:

* http://www.fsf.org/licensing/licenses

*/

/*****************************************************************************

* convolution-impl.h

*

* Implementation for linear convolution.

*

* Zhang Ming, 2010-01, Xi‘an Jiaotong University.

*****************************************************************************/

/**

* convolution and ploynonal multiplication.

*/

template

Vector conv( const Vector &signal, const Vector &filter )

{

if( signal.dim() < filter.dim() )

return convolution( filter, signal );

else

return convolution( signal, filter );

}

template

Vector convolution( const Vector &signal, const Vector &filter )

{

int sigLength = signal.dim();

int filLength = filter.dim();

assert( sigLength >= filLength );

int length = sigLength + filLength - 1;

Vector x(length);

for( int i=1; i<=length; ++i )

{

x(i) = 0;

if( i < filLength )

for( int j=1; j<=i; ++j )

x(i) += filter(j) * signal(i-j+1);

else if( i <= sigLength )

for( int j=1; j<=filLength; ++j )

x(i) += filter(j) * signal(i-j+1);

else

for( int j=i-sigLength+1; j<=filLength; ++j )

x(i) += filter(j) * signal(i-j+1);

}

return x;

}

/**

* Fast convolution by FFT.

*/

template

Vector fastConv( const Vector &xn, const Vector &yn )

{

int M = xn.dim(),

N = yn.dim();

Vector xnPadded = wextend( xn, N-1, "right", "zpd" ),

ynPadded = wextend( yn, M-1, "right", "zpd" );

return ifftc2r( fft(xnPadded) * fft(ynPadded) );

// Vector< complex > Zk = fft(xnPadded) * fft(ynPadded);

// return ifftc2r(Zk);

// return ifftc2r( fft(wextend(xn,N-1,"right","zpd")) * fft(wextend(yn,M-1,"right","zpd")) );

}

测试代码：

/*****************************************************************************

* convolution.cpp

*

* Convolution testing.

*

* Zhang Ming, 2010-01, Xi‘an Jiaotong University.

*****************************************************************************/

#define BOUNDS_CHECK

#include

#include

using namespace std;

using namespace splab;

typedef double Type;

const int M = 3;

const int N = 5;

int main()

{

Vector xn( M ), yn( N );

Vector zn;

for( int i=0; i

xn[i] = i;

for( int i=0; i

yn[i] = i-N/2;

// convolution

zn = conv( xn, yn );

cout << "xn: " << xn << endl << "yn: " << yn << endl;

cout << "convolution of xn and yn: " << zn << endl;

zn = fastConv( xn, yn );

cout << "fast convolution of xn and yn: " << zn << endl;

return 0;

}

运行结果：

xn: size: 3 by 1

0

1

2

yn: size: 5 by 1

-2

-1

0

1

2

convolution of xn and yn: size: 7 by 1

0

-2

-5

-2

1

4

4

fast convolution of xn and yn: size: 7 by 1

-2.53765e-016

-2

-5

-2

1

4

4

Process returned 0 (0x0) execution time : 0.078 s

Press any key to continue.

http://my.oschina.net/zmjerry/blog/3671

原文：http://www.cnblogs.com/pengkunfan/p/3946880.html