频域算法广泛用于图像处理。常用的去噪算法和边缘增强算法的实现都使用2D频域算法。这里介绍一下2D频域算法和快速频域算法的C语言实现。
频域定义
即
步骤：
1）滑动内核，使其中心位于输入图像g的(i,j)像素上；
2）用上式求和c语言实现图像处理c语言实现图像处理，得到输出图像的(i,j)像素值；
3） 完全操作，直到得到输出图像的所有像素值。
1、二维频域算法的实现
bool convolve2DSlow(unsigned char* in, unsigned char* out, int dataSizeX, int dataSizeY,
float* kernel, int kernelSizeX, int kernelSizeY)
{
int i, j, m, n, mm, nn;
int kCenterX, kCenterY; // center index of kernel
float sum; // temp accumulation buffer
int rowIndex, colIndex;
// check validity of params
if(!in || !out || !kernel) return false;
if(dataSizeX <= 0 || kernelSizeX <= 0) return false;
// find center position of kernel (half of kernel size)
kCenterX = kernelSizeX / 2;
kCenterY = kernelSizeY / 2;
for(i=0; i < dataSizeY; ++i) // rows
{
for(j=0; j < dataSizeX; ++j) // columns
{
sum = 0; // init to 0 before sum
for(m=0; m < kernelSizeY; ++m) // kernel rows
{
mm = kernelSizeY - 1 - m; // row index of flipped kernel
for(n=0; n = 0 && rowIndex = 0 && colIndex < dataSizeX)
sum += in[dataSizeX * rowIndex + colIndex] * kernel[kernelSizeX * mm + nn];
}
}
out[dataSizeX * i + j] = (unsigned char)((float)fabs(sum) + 0.5f);
}
}
return true;
}2、二维快速频域算法实现
快速频域算法需要将一个核 2D 核 h 拆分为一个行向量和一个列向量的乘积。实际实现是将二维频域拆分为两个一维频域。设求和的大小为n*n，算法复杂度从n*n增加到2n。
bool convolve2DFast(unsigned char* in, unsigned char* out, int dataSizeX, int dataSizeY,
float* kernel, int kernelSizeX, int kernelSizeY)
{
int i, j, m, n, x, y, t;
unsigned char **inPtr, *outPtr, *ptr;
int kCenterX, kCenterY;
int rowEnd, colEnd; // ending indice for section divider
float sum; // temp accumulation buffer
int k, kSize;
// check validity of params
if(!in || !out || !kernel) return false;
if(dataSizeX <= 0 || kernelSizeX > 1;
kCenterY = kernelSizeY >> 1;
kSize = kernelSizeX * kernelSizeY; // total kernel size
// allocate memeory for multi-cursor
inPtr = new unsigned char*[kSize];
if(!inPtr) return false; // allocation error
// set initial position of multi-cursor, NOTE: it is swapped instead of kernel
ptr = in + (dataSizeX * kCenterY + kCenterX); // the first cursor is shifted (kCenterX, kCenterY)
for(m=0, t=0; m < kernelSizeY; ++m)
{
for(n=0; n < kernelSizeX; ++n, ++t)
{
inPtr[t] = ptr - n;
}
ptr -= dataSizeX;
}
// init working pointers
outPtr = out;
rowEnd = dataSizeY - kCenterY; // bottom row partition divider
colEnd = dataSizeX - kCenterX; // right column partition divider
// convolve rows from index=0 to index=kCenterY-1
y = kCenterY;
for(i=0; i < kCenterY; ++i)
{
// partition #1 ***********************************
x = kCenterX;
for(j=0; j < kCenterX; ++j) // column from index=0 to index=kCenterX-1
{
sum = 0;
t = 0;
for(m=0; m <= y; ++m)
{
for(n=0; n <= x; ++n)
{
sum += *inPtr[t] * kernel[t];
++t;
}
t += (kernelSizeX - x - 1); // jump to next row
}
// store output
*outPtr = (unsigned char)((float)fabs(sum) + 0.5f);
++outPtr;
++x;
for(k=0; k < kSize; ++k) ++inPtr[k]; // move all cursors to next
}
// partition #2 ***********************************
for(j=kCenterX; j < colEnd; ++j) // column from index=kCenterX to index=(dataSizeX-kCenterX-1)
{
sum = 0;
t = 0;
for(m=0; m <= y; ++m)
{
for(n=0; n < kernelSizeX; ++n)
{
sum += *inPtr[t] * kernel[t];
++t;
}
}
// store output
*outPtr = (unsigned char)((float)fabs(sum) + 0.5f);
++outPtr;
++x;
for(k=0; k < kSize; ++k) ++inPtr[k]; // move all cursors to next
}
// partition #3 ***********************************
x = 1;
for(j=colEnd; j < dataSizeX; ++j) // column from index=(dataSizeX-kCenter) to index=(dataSizeX-1)
{
sum = 0;
t = x;
for(m=0; m <= y; ++m)
{
for(n=x; n < kernelSizeX; ++n)
{
sum += *inPtr[t] * kernel[t];
++t;
}
t += x; // jump to next row
}
// store output
*outPtr = (unsigned char)((float)fabs(sum) + 0.5f);
++outPtr;
++x;
for(k=0; k < kSize; ++k) ++inPtr[k]; // move all cursors to next
}
++y; // add one more row to convolve for next run
}
// convolve rows from index=kCenterY to index=(dataSizeY-kCenterY-1)
for(i= kCenterY; i < rowEnd; ++i) // number of rows
{
// partition #4 ***********************************
x = kCenterX;
for(j=0; j < kCenterX; ++j) // column from index=0 to index=kCenterX-1
{
sum = 0;
t = 0;
for(m=0; m < kernelSizeY; ++m)
{
for(n=0; n <= x; ++n)
{
sum += *inPtr[t] * kernel[t];
++t;
}
t += (kernelSizeX - x - 1);
}
// store output
*outPtr = (unsigned char)((float)fabs(sum) + 0.5f);
++outPtr;
++x;
for(k=0; k < kSize; ++k) ++inPtr[k]; // move all cursors to next
}
// partition #5 ***********************************
for(j = kCenterX; j < colEnd; ++j) // column from index=kCenterX to index=(dataSizeX-kCenterX-1)
{
sum = 0;
t = 0;
for(m=0; m < kernelSizeY; ++m)
{
for(n=0; n < kernelSizeX; ++n)
{
sum += *inPtr[t] * kernel[t];
++inPtr[t]; // in this partition, all cursors are used to convolve. moving cursors to next is safe here
++t;
}
}
// store output
*outPtr = (unsigned char)((float)fabs(sum) + 0.5f);
++outPtr;
++x;
}
// partition #6 ***********************************
x = 1;
for(j=colEnd; j < dataSizeX; ++j) // column from index=(dataSizeX-kCenter) to index=(dataSizeX-1)
{
sum = 0;
t = x;
for(m=0; m < kernelSizeY; ++m)
{
for(n=x; n < kernelSizeX; ++n)
{
sum += *inPtr[t] * kernel[t];
++t;
}
t += x;
}
// store output
*outPtr = (unsigned char)((float)fabs(sum) + 0.5f);
++outPtr;
++x;
for(k=0; k < kSize; ++k) ++inPtr[k]; // move all cursors to next
}
}
// convolve rows from index=(dataSizeY-kCenterY) to index=(dataSizeY-1)
y = 1;
for(i= rowEnd; i < dataSizeY; ++i) // number of rows
{
// partition #7 ***********************************
x = kCenterX;
for(j=0; j < kCenterX; ++j) // column from index=0 to index=kCenterX-1
{
sum = 0;
t = kernelSizeX * y;
for(m=y; m < kernelSizeY; ++m)
{
for(n=0; n <= x; ++n)
{
sum += *inPtr[t] * kernel[t];
++t;
}
t += (kernelSizeX - x - 1);
}
// store output
*outPtr = (unsigned char)((float)fabs(sum) + 0.5f);
++outPtr;
++x;
for(k=0; k < kSize; ++k) ++inPtr[k]; // move all cursors to next
}
// partition #8 ***********************************
for(j=kCenterX; j < colEnd; ++j) // column from index=kCenterX to index=(dataSizeX-kCenterX-1)
{
sum = 0;
t = kernelSizeX * y;
for(m=y; m < kernelSizeY; ++m)
{
for(n=0; n < kernelSizeX; ++n)
{
sum += *inPtr[t] * kernel[t];
++t;
}
}
// store output
*outPtr = (unsigned char)((float)fabs(sum) + 0.5f);
++outPtr;
++x;
for(k=0; k < kSize; ++k) ++inPtr[k];
}
// partition #9 ***********************************
x = 1;
for(j=colEnd; j < dataSizeX; ++j) // column from index=(dataSizeX-kCenter) to index=(dataSizeX-1)
{
sum = 0;
t = kernelSizeX * y + x;
for(m=y; m < kernelSizeY; ++m)
{
for(n=x; n < kernelSizeX; ++n)
{
sum += *inPtr[t] * kernel[t];
++t;
}
t += x;
}
// store output
*outPtr = (unsigned char)((float)fabs(sum) + 0.5f);
++outPtr;
++x;
for(k=0; k < kSize; ++k) ++inPtr[k]; // move all cursors to next
}
++y; // the starting row index is increased
}
return true;
}2017.10.29
