数字图像处理学习笔记(1.2)---位图的读写、几何变换、傅里叶变换、直方图均衡

#include"bmp.h"
#include<cmath>
#include<cstring>
#include<cstdio>

#define PI 3.1415926

//说明：对输入图像进行快速傅立叶变换，要求输入图像的宽和高必须是2的幂次方

void Bitmap::fourier()
{
int lineByte = (width_p*bitCount / 8 + 3) / 4 * 4;

//申请输出图像缓冲区，并初始化为0
unsigned char* m_pImgDataOut = new unsigned char[lineByte*height_p];
memset(m_pImgDataOut, 0, width_p*height_p);

//申请傅立叶缓冲区，并初始化为0
ComplexNumber* m_pFFTBuf = new ComplexNumber[width_p*height_p];
memset(m_pFFTBuf, 0, sizeof(ComplexNumber)*width_p*height_p);

//输入图像数据进行二维傅立叶变换
ImgFFT2D(dataBuf, m_pImgDataOut, m_pFFTBuf);
delete[] m_pFFTBuf;
delete dataBuf;
dataBuf = m_pImgDataOut;
}

//说明：对复数结构体数组arrayBuf进行一维快速傅立叶变换，变换后的结果仍存回arrayBuf中

void Bitmap::FFT1D(ComplexNumber *arrayBuf, int n)
{
//循环变量
int i, k, r;

//申请临时复数缓冲区buf1，长度为n
ComplexNumber *buf1 = new ComplexNumber
;

//将arrayBuf拷贝进buf1
memcpy(buf1, arrayBuf, sizeof(ComplexNumber)*n);

//申请临时复数缓冲区buf2，长度为n
ComplexNumber *buf2 = new ComplexNumber
;

//将arrayBuf数组元素基2抽取并重新排列
//若0、1、2、3、4、5、6、7八点序列对调后变作0、4、2、6、1、5、3、7
int t1, t2;
for (r = 1; pow(2, r)<n; r++){
t1 = pow(2, r);
t2 = pow(2, r - 1);
for (k = 0; k<t1; k++){
for (i = 0; i<n / t1; i++){
buf2[k*n / t1 + i].real = buf1[k / 2 * n / t2 + i * 2 + k % 2].real;
buf2[k*n / t1 + i].imag = buf1[k / 2 * n / t2 + i * 2 + k % 2].imag;
}
}
memcpy(buf1, buf2, sizeof(ComplexNumber)*n);
}

//采用蝶型算法进行快速傅立叶变换
//buf1是第r级的输入，buf2存放第r级的输出
float c, s;
for (r = 1; pow(2, r) <= n; r++){
t1 = pow(2, r);
for (k = 0; k<n / t1; k++){
for (i = t1 / 2; i<t1; i++){
//加权因子
c = cos(-2 * PI*(i - t1 / 2) / t1);
s = sin(-2 * PI*(i - t1 / 2) / t1);
buf1[k*t1 + i].real = buf2[k*t1 + i].real*c - buf2[k*t1 + i].imag*s;
buf1[k*t1 + i].imag = buf2[k*t1 + i].imag*c + buf2[k*t1 + i].real*s;
}
}
for (k = 0; k<n / t1; k++){
for (i = 0; i<t1 / 2; i++){
buf2[k*t1 + i].real = buf1[k*t1 + i].real + buf1[k*t1 + i + t1 / 2].real;
buf2[k*t1 + i].imag = buf1[k*t1 + i].imag + buf1[k*t1 + i + t1 / 2].imag;
}
for (i = t1 / 2; i<t1; i++){
buf2[k*t1 + i].real = buf1[k*t1 + i - t1 / 2].real - buf1[k*t1 + i].real;
buf2[k*t1 + i].imag = buf1[k*t1 + i - t1 / 2].imag - buf1[k*t1 + i].imag;
}
}

//第r级的输出存入buf1,作为下一级的输入数据
memcpy(buf1, buf2, sizeof(ComplexNumber)*n);
}

//傅立叶变换的结果存入arrayBuf
memcpy(arrayBuf, buf2, sizeof(ComplexNumber)*n);

//释放缓冲区
delete[]buf2;
delete[]buf1;

}

//说明：图像数据二维快速傅立叶变换将图像数据变成复数形式,进行两个一维快速
//   傅立叶变换,变换后的频谱以图像形式存入imgBufOut，此处要求图像宽和高
//   都为2的幂次方

void Bitmap::ImgFFT2D(unsigned char* imgBuf,unsigned char *imgBufOut,ComplexNumber* m_pFFTBuf)
{
//循环变量
int i, j, u, v;

//图像数据变成复数类型存入m_pFFTBuf
for (i = 0; i<width_p*height_p; i++){
m_pFFTBuf[i].real = imgBuf[i];
m_pFFTBuf[i].imag = 0;
}

//申请ComplexNumber结构体数组,长度为height_p
ComplexNumber *array = new ComplexNumber[height_p];

//先纵向一维快速傅立叶变换
for (u = 0; u<width_p; u++){
for (v = 0; v<height_p; v++){
array[v].real = m_pFFTBuf[v*width_p + u].real;
array[v].imag = m_pFFTBuf[v*width_p + u].imag;
}
FFT1D(array, height_p);
for (v = 0; v<height_p; v++){
m_pFFTBuf[v*width_p + u].real = array[v].real;
m_pFFTBuf[v*width_p + u].imag = array[v].imag;
}
}
delete[]array;

//再横向一维快速傅立叶变换
for (v = 0; v<height_p; v++){
FFT1D(m_pFFTBuf + v*width_p, width_p);
}

//将频谱图以图像形式存入imgBufOut
float t;
int i0, j0;
for (i = 0; i<height_p; i++){
//i0 = i;
//j0 = j;
for (j = 0; j<width_p; j++){
if (i<height_p / 2)
i0 = i + height_p / 2;
else
i0 = i - height_p / 2;
if (j<width_p / 2)
j0 = j + width_p / 2;
else
j0 = j - width_p / 2;

t = sqrt(m_pFFTBuf[i0*width_p + j0].real*m_pFFTBuf[i0*width_p + j0].real
+ m_pFFTBuf[i0*width_p + j0].imag*m_pFFTBuf[i0*width_p + j0].imag);
t = t / 500;
if (t>255)
imgBufOut[i*width_p + j] = 255;
else
imgBufOut[i*width_p + j] = t;
}
}

}

#include"bmp.h"
#include<iostream>

using namespace std;

int main()
{
char* fileName = "qianxun.bmp";
Bitmap* bmp = new Bitmap();
bmp->read(fileName);
//bmp->translation(10, 10);
//bmp->zoom(2, 2, 'n');//'n','l','c';默认为'n'邻近插值缩放
//bmp->rotate(90);
bmp->fourier();
bmp->write("fourier.bmp");
delete bmp;
return 1;
}