2019-09-27 21:48:56 Jammy_Frazer 阅读数 61
'''
Implementation of the secend class's image process methods, which contains:
1. furior transform
2. walsh transform
3. DCT transform
'''
import cv2
import numpy as np
import matplotlib.pyplot as plt

# load test img and transform it to grey scale
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
image = cv2.resize(image, (512, 512))
print(np.shape(image))

plt.subplot(221), plt.imshow(image, cmap='gray')
plt.title('origin'), plt.xticks([]), plt.yticks([])

# 1. furior transform
# do dft and show the magnitude spectrum
dft = cv2.dft(np.float32(image), flags=cv2.DFT_COMPLEX_OUTPUT)
magnitude_spectrum1 = 20*np.log(cv2.magnitude(dft[:, :, 0], dft[:, :, 1]))

dft_shift = np.fft.fftshift(dft)
magnitude_spectrum2 = 20*np.log(cv2.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))

#plt.subplot(121), plt.imshow(magnitude_spectrum1, cmap='gray')
#plt.title('dft'), plt.xticks([]), plt.yticks([])
plt.subplot(222), plt.imshow(magnitude_spectrum2, cmap='gray')
plt.title('dft_shift'), plt.xticks([]), plt.yticks([])

'''
# simple averaging filter without scaling parameter
mean_filter = np.ones((3,3))

# creating a guassian filter
x = cv2.getGaussianKernel(5,10)
gaussian = x*x.T

# different edge detecting filters
# scharr in x-direction
scharr = np.array([[-3, 0, 3],
[-10,0,10],
[-3, 0, 3]])
# sobel in x direction
sobel_x= np.array([[-1, 0, 1],
[-2, 0, 2],
[-1, 0, 1]])
# sobel in y direction
sobel_y= np.array([[-1,-2,-1],
[0, 0, 0],
[1, 2, 1]])
# laplacian ( high pass )
laplacian=np.array([[0, 1, 0],
[1,-4, 1],
[0, 1, 0]])

filters = [mean_filter, gaussian, laplacian, sobel_x, sobel_y, scharr]
filter_name = ['mean_filter', 'gaussian','laplacian', 'sobel_x', \
'sobel_y', 'scharr_x']
fft_filters = [np.fft.fft2(x) for x in filters]
fft_shift = [np.fft.fftshift(y) for y in fft_filters]
mag_spectrum = [np.log(np.abs(z)+1) for z in fft_shift]

for i in range(6):
plt.subplot(2, 3, i+1), plt.imshow(mag_spectrum[i], cmap = 'gray')
plt.title(filter_name[i]), plt.xticks([]), plt.yticks([])
plt.show()
'''

# 2. walsh transform
# something wrong with this part
def joint(h):
up = np.hstack((h, h))
down = np.hstack((h, -h))
H = np.vstack((up, down))
return H

def gen_walsh_matrix(size=512):
H = np.array([[1,1],[1,-1]])
iter = int(np.log2(size))-1
for i in range(iter):
H = joint(H)
return H
size = np.shape(image)[0]
H = gen_walsh_matrix(512)

print(image)
image_wt = H.T.dot(image).dot(H)/ 9
print(image_wt)

plt.subplot(223), plt.imshow(image_wt, cmap='gray')
plt.title('walsh transform'), plt.xticks([]), plt.yticks([])

# 3. dct with opencv
img_dct = cv2.dct(image.astype(np.float32))
img_dct_log = np.log(abs(img_dct))

plt.subplot(224), plt.imshow(img_dct_log, cmap='gray')
plt.title('dct transform'), plt.xticks([]), plt.yticks([])
plt.show()



2018-04-24 20:36:38 weixin_39569242 阅读数 1974

1）了解图像变换的意义和手段。

2）熟悉傅立叶变换的基本性质。

3）通过实验了解二维频谱的分布特点。

任意选择几幅图像，对其进行傅立叶变换，观察图像的频谱特征。

1）对图像进行傅里叶变换（包括移位和未移位），观察频谱信号

2）观察频谱的三维图形

①　移位前：

v 代码：

f=zeros(64,64);

f(15:50,15:50)=1;%输入64*64的黑色图像矩阵

F=fft2(f);%平移前傅立叶变换

F2=fftshift(abs(F));%频谱中心化

x=1:64;y=1:64;

%绘制平移前

figure(1)

subplot(1,3,1),imshow(f);title('平移前原始图像');

subplot(1,3,2),imshow(abs(F));title('图像傅里叶变换');

subplot(1,3,3),imshow(abs(F2));title('图像中心化傅里叶变换');

figure(2)

subplot(2,1,1),title('三维图像');

mesh(abs(real(F)),'FaceColor','white');

subplot(2,1,2),title('FFT变换');

mesh(x,y,F2(x,y),'FaceColor','red');

图1

2：彩色三维图像                          图3FFT变换

②　移位后：

ff=circshift(f,[-10 -5]);%将原始图像向左向上移，及X轴，Y轴负方向

FF=fft2(ff);%平移后傅立叶变换

FF2=fftshift(abs(FF));%频谱中心化

%绘制平移后

figure(3)

subplot(2,3,1),imshow(ff);title('平移后原始图像');

subplot(2,3,2),imshow(abs(FF));title('图像傅里叶变换');

subplot(2,3,3),imshow(abs(FF2));title('图像中心化傅里叶变换');

figure(4)

subplot(2,1,1),title('三维图像');

mesh(abs(real(FF)),'FaceColor','white');

subplot(2,1,2),title('FFT变换');

mesh(x,y,FF2(x,y),'FaceColor','red');

图4

图5：彩色三维图像                            图6FFT变换

3）进行傅里叶逆变换，重建图像

4）根据幅值谱重建图像信息

5）根据相位谱重建图像信息

%平移之前的频谱、相位谱、频谱逆变换、相谱逆变换、傅里叶逆变换

xf1=log(abs(F));%频谱

xf2=angle(F);%相位谱

xf12=ifft2(xf1);

xf13=ifftshift(abs(xf12));%频谱逆变换

xf22=ifft2(xf2);

xf23=ifftshift(abs(xf22));%相位逆变换

xr1=ixff1=log(abs(FF));%平移后频谱

%平移之后频谱、相位谱

xff2=angle(FF);%平移后相位谱

fft2(F);%进行傅里叶逆变换

%平移前平移后相位谱、频谱对比

figure(5)

subplot(2,2,1),imshow(abs(xf2));title('傅里叶相位谱');

subplot(2,2,2),imshow(abs(xff2));title('平移后傅里叶相位谱');

subplot(2,2,3),imshow(abs(xf1));title('傅里叶幅度谱');

subplot(2,2,4),imshow(abs(xff1));title('平移后傅里叶幅度谱');

%平移前逆变换、相位谱逆变换、幅度谱逆变换

figure(6)

subplot(2,2,1),imshow(abs(xf23));title('傅里叶相位谱逆变换');

subplot(2,2,2),imshow(abs(xf13));title('傅里叶幅度谱逆变换');

subplot(2,2,3);imshow(abs(xr1));title('逆变换');

7：平移前后相位谱、幅度谱

8：平移前逆变换、相位谱逆变换、幅度谱逆变换

1、图像平移之后的傅里叶幅度谱不会发生变化，而仅仅是相位谱产生了一定的相移特性，如图7

2、图像进行傅里叶变换后，进行傅里叶逆变换会还原为原来的图像，如图8

3、相位谱逆变换之后会大致还原原图像的轮廓、幅度谱逆变换之后还原不出原图像，只会确定图像的精度。如图8

f=zeros(64,64);

f(15:50,15:50)=1;%输入64*64的黑色图像矩阵

F=fft2(f);%平移前傅立叶变换

F2=fftshift(abs(F));%频谱中心化

x=1:64;

y=1:64;

%平移之前的频谱、相位谱、频谱逆变换、相谱逆变换、傅里叶逆变换

xf1=log(abs(F));%频谱

xf2=angle(F);%相位谱

xf12=ifft2(xf1);

xf13=ifftshift(abs(xf12));%频谱逆变换

xf22=ifft2(xf2);

xf23=ifftshift(abs(xf22));%相位逆变换

xr1=ifft2(F);%进行傅里叶逆变换

%平移之后的傅里叶变换频谱、相位谱

ff=circshift(f,[-10 -5]);%将原始图像向左向上移，及X轴，Y轴负方向

FF=fft2(ff);%平移后傅立叶变换

FF2=fftshift(abs(FF));%频谱中心化

xff1=log(abs(FF));%平移后频谱

xff2=angle(FF);%平移后相位谱

%绘制平移前

figure(1)

subplot(1,3,1),imshow(f);title('平移前原始图像');

subplot(1,3,2),imshow(abs(F));title('图像傅里叶变换');

subplot(1,3,3),imshow(abs(F2));title('图像中心化傅里叶变换');

figure(2)

subplot(2,1,1),title('三维图像');

mesh(abs(real(F)),'FaceColor','white');

subplot(2,1,2),title('FFT变换');

mesh(x,y,F2(x,y),'FaceColor','red');

%绘制平移后

figure(3)

subplot(2,3,1),imshow(ff);title('平移后原始图像');

subplot(2,3,2),imshow(abs(FF));title('图像傅里叶变换');

subplot(2,3,3),imshow(abs(FF2));title('图像中心化傅里叶变换');

figure(4)

subplot(2,1,1),title('三维图像');

mesh(abs(real(FF)),'FaceColor','white');

subplot(2,1,2),title('FFT变换');

mesh(x,y,FF2(x,y),'FaceColor','red');

%平移前平移后相位谱、频谱对比

figure(5)

subplot(2,2,1),imshow(abs(xf2));title('傅里叶相位谱');

subplot(2,2,2),imshow(abs(xff2));title('平移后傅里叶相位谱');

subplot(2,2,3),imshow(abs(xf1));title('傅里叶幅度谱');

subplot(2,2,4),imshow(abs(xff1));title('平移后傅里叶幅度谱');

%平移前逆变换、相位谱逆变换、幅度谱逆变换

figure(6)

subplot(2,2,1),imshow(abs(xf23));title('傅里叶相位谱逆变换');

subplot(2,2,2),imshow(abs(xf13));title('傅里叶幅度谱逆变换');

subplot(2,2,3);imshow(abs(xr1));title('逆变换');

2014-10-14 09:52:54 u012596983 阅读数 14068

H.264中使用了4阶和8阶的阿达马变换来计算SATD，其变换矩阵为：

## SATD计算方法

function xx=fwhtdya2d(data2)
xx=data2;
N=length(xx);
for i=1:N
xx1(i,:)=fwhtdyald(xx(i,:));
end
xx=zeros(N);
for j=1:N
xx(:,j)=fwhtdyald(xx1(:,j));
end

function x=fwhtdyald(data1)
N=length(data1);
x=bitrevorder(data1);
L=log2(N);
k1=N;k2=1;k3=N/2;
for i1=1:L
L1=1;
for i2=1:k2
for i3=1:k3
i=i3+L1-1;j=i+k3;
temp1=x(i);temp2=x(j);
x(i)=temp1+temp2;
x(j)=temp1-temp2;
end
L1=L1+k1;
end
k1=k1/2;k2=k2*2;k3=k3/2;
end
x=inv(N)*x;

walsh变换在图像压缩领域的应用：

I=imread('lena.png');
sig=rgb2gray(I);
sig=double(sig)/255;
%figure,imshow(sig);
[m_sig,n_sig]=size(sig);
sizi=16;
Snum=128;
hdcoe=blkproc(sig,[sizi sizi],'P1*x*P2',T,T);
coe=im2col(hdcoe,[sizi sizi],'distinct');
coe_temp=coe;
[Y,Ind]=sort(coe);
%舍去较小方差的系数
[m,n]=size(coe);
Snum=m-Snum;
for i=1:n
coe_temp(Ind(1:Snum),i)=0;
end
re_hdcoe=col2im(coe_temp,[sizi sizi],[m_sig n_sig],'distinct');
re_sig=blkproc(re_hdcoe,[sizi sizi],'P1*x*P2',T,T);
%figure,imshow(uint8(re_sig));
error=sig.^2-re_sig.^2;
MSE=sum(error(:)/prod(size(re_sig)))
subplot(1,2,1),imshow(sig);
subplot(1,2,2),imshow(uint8(re_sig));

2015-05-02 11:16:53 fztfztfzt 阅读数 3083
/*************************************************************************	*	* 函数名称：	*   WALSH()	*	* 参数:	*   double* f			- 输入的时域序列	*   double* F			- 输出的频域序列	*   int r				- 2的幂数			*	* 返回值:	*   BOOL               - 成功返回TRUE，否则返回FALSE。	*	* 说明:	*   该函数进行一维快速沃尔什——哈达马变换。	*	************************************************************************/void WALSH(double*f,double*F,int r);/*************************************************************************	*	* 函数名称：	*   IWALSH()	*	* 参数:	*   double* f			- 输入的时域序列	*   double* F			- 输出的频域序列	*   int r				- 2的幂数			*	* 返回值:	*   BOOL               - 成功返回TRUE，否则返回FALSE。	*	* 说明:	*   该函数进行一维快速沃尔什——哈达马逆变换。	*	************************************************************************/void IWALSH(double*F,double*f,int r);/*************************************************************************	*	* 函数名称：	*   FreqWALSH()	*	* 参数:	*   double* f			- 输入的时域序列	*   double* F			- 输出的频域序列	*	 LONG width			- 图象宽度	*	 LONG height		- 图象高度	*	* 返回值:	*   BOOL               - 成功返回TRUE，否则返回FALSE。	*	* 说明:	*   该函数进行二维快速沃尔什——哈达玛变换。	*	************************************************************************/	BOOL FreqWALSH(double*f,double*F, LONG width, LONG height);/*************************************************************************	*	* 函数名称：	*   IFreqWALSH()	*	* 参数:	*   double* f			- 输入的时域序列	*   double* F			- 输出的频域序列	*	 LONG width		- 图象宽度	*	 LONG height		- 图象高度	*	* 返回值:	*   BOOL               - 成功返回TRUE，否则返回FALSE。	*	* 说明:	*   该函数进行二维快速沃尔什——哈达玛逆变换。	*	************************************************************************/	BOOL IFreqWALSH(double*f,double*F, LONG width, LONG height);/*************************************************************************	*	* 函数名称：	*   BmpWalsh()	*	* 参数:	*   BYTE* bmp,LONG width,LONG height	*	* 返回值:	*   BOOL               - 成功返回TRUE，否则返回FALSE。	*	* 说明:	*   该函数对图象进行二维快速沃尔什——哈达马变换。	*	************************************************************************/	BOOL BmpWalsh(BYTE* bmp,LONG width,LONG height);
voidMyProcess::WALSH(double*f,double*F,int r){// 循环变量	LONG	i;	LONG	j;	LONG	k;// 中间变量int		p;double* X;// 计算快速沃尔什变换点数	LONG N =1<< r;// 分配运算所需的数组double* X1 =newdouble[N];double* X2 =newdouble[N];// 将时域点写入数组X1	memcpy(X1, f,sizeof(double)* N);// 蝶形运算for(k =0; k < r; k++){for(j =0; j <1<<k; j++){for(i =0; i <1<<(r - k -1); i++){				p = j *(1<<(r - k));				X2[i + p]= X1[i + p]+ X1[i + p +(int)(1<<(r - k -1))];				X2[i + p +(int)(1<<(r - k -1))]= X1[i + p]- X1[i + p +(int)(1<<(r - k -1))];}}// 互换X1和X2  		X = X1;		X1 = X2;		X2 = X;}// 调整系数for(j =0; j < N; j++){		p =0;for(i =0; i < r; i++){if(j &(1<<i)){				p +=1<<(r - i -1);}}		F[j]= X1[p]/ N;}// 释放内存delete X1;delete X2;}voidMyProcess::IWALSH(double*F,double*f,int r){// 循环变量int		i;// 计算变换点数	LONG N =1<< r;// 调用快速沃尔什－哈达玛变换进行反变换	WALSH(F, f, r);// 调整系数for(i =0; i < N; i++)		f[i]*= N;}BOOL MyProcess::FreqWALSH(double*f,double*F, LONG width, LONG height){// 循环变量	LONG	i;	LONG	j;	LONG    k;// 进行离散余弦变换的宽度和高度（2的整数次方）	LONG w =1;	LONG h =1;int wp =0;int hp =0;// 计算进行离散余弦变换的宽度和高度（2的整数次方）while(w < width/3){		w *=2;		wp++;}while(h < height){		h *=2;		hp++;}// 分配内存double*TempIn=newdouble[h];double*TempOut=newdouble[h];// 对y方向进行离散余弦变换for(i =0; i < w *3; i++){// 抽取数据for(j =0; j < h; j++)TempIn[j]= f[j * w *3+ i];// 一维快速离散余弦变换		WALSH(TempIn,TempOut, hp);// 保存变换结果for(j =0; j < h; j++)			f[j * w *3+ i]=TempOut[j];}// 释放内存deleteTempIn;deleteTempOut;// 分配内存TempIn=newdouble[w];TempOut=newdouble[w];// 对x方向进行快速离散余弦变换for(i =0; i < h; i++){for(k =0; k <3; k++){// 抽取数据for(j =0; j < w; j++)TempIn[j]= f[i * w *3+ j *3+ k];// 一维快速离散余弦变换			WALSH(TempIn,TempOut, wp);// 保存变换结果for(j =0; j < w; j++)				F[i * w *3+ j *3+ k]=TempOut[j];}}// 释放内存deleteTempIn;deleteTempOut;return TRUE;}BOOL MyProcess::IFreqWALSH(double*f,double*F, LONG width, LONG height){// 循环变量	LONG	i;	LONG	j;	LONG    k;// 赋初值	LONG w =1;	LONG h =1;int wp =0;int hp =0;// 计算进行付立叶变换的宽度和高度（2的整数次方）while(w < width/3){		w *=2;		wp++;}while(h < height){		h *=2;		hp++;}// 分配内存double*TempIn=newdouble[w];double*TempOut=newdouble[w];// 对x方向进行快速付立叶变换for(i =0; i < h; i++){for(k =0; k <3; k++){// 抽取数据for(j =0; j < w; j++)TempIn[j]= F[i * w *3+ j *3+ k];// 一维快速傅立叶变换			IWALSH(TempIn,TempOut, wp);// 保存变换结果for(j =0; j < w; j++)				F[i * w *3+ j *3+ k]=TempOut[j];}}// 释放内存deleteTempIn;deleteTempOut;TempIn=newdouble[h];TempOut=newdouble[h];// 对y方向进行快速付立叶变换for(i =0; i < w *3; i++){// 抽取数据for(j =0; j < h; j++)TempIn[j]= F[j * w *3+ i];// 一维快速傅立叶变换		IWALSH(TempIn,TempOut, hp);// 保存变换结果for(j =0; j < h; j++)			F[j * w *3+ i]=TempOut[j];}// 释放内存deleteTempIn;deleteTempOut;for(i =0; i < h; i++){for(j =0; j < w *3; j++){if(i < height && j < width)*(f + i * width + j)= F[i * w *3+ j];}}return TRUE;}BOOL MyProcess::BmpWalsh(BYTE* bmp,LONG width,LONG height){// 循环变量	LONG	i;	LONG	j;// 进行沃尔什——哈达玛变换的宽度和高度（2的整数次方）	LONG w =1;	LONG h =1;int wp =0;int hp =0;// 计算进行离散余弦变换的宽度和高度（2的整数次方）while(w < width/3){		w *=2;		wp++;}while(h < height){		h *=2;		hp++;}// 分配内存double*f =newdouble[w * h *3];double*F =newdouble[w * h *3];// 向时域赋值并补零for(i =0; i < h; i++){for(j =0; j < w *3; j++){if(i < height && j < width)				f[i * w *3+ j]= bmp[width * i + j];else				f[w * i *3+ j]=0.0f;}}// 进行频谱分析if(FreqWALSH(f, F,width, height)== FALSE)return FALSE;// 更新所有象素for(i =0; i < height; i++){for(j =0; j < width; j++){// 判断是否超过255if(fabs(F[i * w *3+ j]*1000)>255){// 对于超过的，直接设置为255				bmp[width *(height -1- i)+ j]=255;}else{// 如果没有超过，则按实际计算结果赋值				bmp[width *(height -1- i)+ j]= fabs(F[i * w *3+ j]*1000);}}}//释放内存delete[] f;delete[] F;// 返回return TRUE;}

2016-10-24 17:18:48 liushuiwen101423 阅读数 6458

数字图像处理领域的二十四个典型算法及vc实现、第一章

--------------------------------------------------

像处理，是对图像进行分析、加工、和处理，使其满足视觉、心理以及其他要求的技术。图像处理是信号处理在图像域上的一个应用。目前大多数的图像是以数字形式存储，因而图像处理很多情况下指数字图像处理。

本文接下来，简单粗略介绍下数字图像处理领域中的24个经典算法，然后全部算法用vc实现。由于篇幅所限，只给出某一算法的主体代码。

ok，请细看。

算法介绍（

什么叫灰度图？任何颜色都有红、绿、蓝三原色组成，假如原来某点的颜色为RGB(R，G，B)，那么，我们可以通过下面几种方法，将其转换为灰度：
1.浮点算法：Gray=R*0.3+G*0.59+B*0.11
2.整数方法：Gray=(R*30+G*59+B*11)/100
3.移位方法：Gray =(R*28+G*151+B*77)>>8;
4.平均值法：Gray=（R+G+B）/3;
5.仅取绿色：Gray=G；
通过上述任一种方法求得Gray后，将原来的RGB(R,G,B)中的R,G,B统一用Gray替换，形成新的颜色RGB(Gray,Gray,Gray)，用它替换原来的RGB(R,G,B)就是灰度图了。

程序实现：
ok，知道了什么叫灰度图，下面，咱们就来实现此256色灰度图。

void Convert256toGray(HDIB hDIB)
{
LPSTR lpDIB;
// 由DIB句柄得到DIB指针并锁定DIB
lpDIB = (LPSTR) ::GlobalLock((HGLOBAL)hDIB);
// 指向DIB象素数据区的指针
LPSTR   lpDIBBits;
// 指向DIB象素的指针
BYTE * lpSrc;
// 图像宽度
LONG lWidth;
// 图像高度
LONG   lHeight;
// 图像每行的字节数
LONG lLineBytes;
// 指向BITMAPINFO结构的指针（Win3.0）
LPBITMAPINFO lpbmi;
// 指向BITMAPCOREINFO结构的指针
LPBITMAPCOREINFO lpbmc;
// 获取指向BITMAPINFO结构的指针（Win3.0）
lpbmi = (LPBITMAPINFO)lpDIB;
// 获取指向BITMAPCOREINFO结构的指针
lpbmc = (LPBITMAPCOREINFO)lpDIB;
// 灰度映射表
BYTE bMap[256];

// 计算灰度映射表（保存各个颜色的灰度值），并更新DIB调色板
int i,j;
for (i = 0; i < 256; i ++)
{
// 计算该颜色对应的灰度值
bMap[i] = (BYTE)(0.299 * lpbmi->bmiColors[i].rgbRed +
0.587 * lpbmi->bmiColors[i].rgbGreen +
0.114 * lpbmi->bmiColors[i].rgbBlue + 0.5);
// 更新DIB调色板红色分量
lpbmi->bmiColors[i].rgbRed = i;
// 更新DIB调色板绿色分量
lpbmi->bmiColors[i].rgbGreen = i;
// 更新DIB调色板蓝色分量
lpbmi->bmiColors[i].rgbBlue = i;
// 更新DIB调色板保留位
lpbmi->bmiColors[i].rgbReserved = 0;
}
// 找到DIB图像象素起始位置
lpDIBBits = ::FindDIBBits(lpDIB);
// 获取图像宽度
lWidth = ::DIBWidth(lpDIB);
// 获取图像高度
lHeight = ::DIBHeight(lpDIB);
// 计算图像每行的字节数
lLineBytes = WIDTHBYTES(lWidth * 8);
// 更换每个象素的颜色索引（即按照灰度映射表换成灰度值）

//逐行扫描
for(i = 0; i < lHeight; i++)
{
//逐列扫描
for(j = 0; j < lWidth; j++)
{
// 指向DIB第i行，第j个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
// 变换
*lpSrc = bMap[*lpSrc];
}
}
//解除锁定
::GlobalUnlock ((HGLOBAL)hDIB);

算法介绍：
有关Walsh变换的深入介绍，请看此论文：http://www.informatics.org.cn/doc/ucit200510/ucit20051005.pdf

程序实现:

double * f    - 指向时域值的指针
double * F    - 指向频域值的指针
r      －2的幂数

VOID WINAPI WALSH(double *f, double *F, int r)
{
// 沃尔什-哈达玛变换点数
LONG count;
// 循环变量
int  i,j,k;
// 中间变量
int  bfsize,p;
double *X1,*X2,*X;
// 计算快速沃尔什变换点数
count = 1 << r;
// 分配运算所需的数组
X1 = new double[count];
X2 = new double[count];
// 将时域点写入数组X1
memcpy(X1, f, sizeof(double) * count);

// 蝶形运算
for(k = 0; k < r; k++)
{
for(j = 0; j < 1<<k; j++)
{
bfsize = 1 << (r-k);
for(i = 0; i < bfsize / 2; i++)
{
p = j * bfsize;
X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
X2[i + p + bfsize / 2] = X1[i + p] - X1[i + p + bfsize / 2];
}
}
// 互换X1和X2
X = X1;
X1 = X2;
X2 = X;
}

// 调整系数
for(j = 0; j < count; j++)
{
p = 0;
for(i = 0; i < r; i++)
{
if (j & (1<<i))
{
p += 1 << (r-i-1);
}
}

F[j] = X1[p] / count;
}

// 释放内存
delete X1;
delete X2;
}

LPSTR lpDIBBits    - 指向源DIB图像指针
LONG  lWidth       - 源图像宽度（象素数）
LONG  lHeight      - 源图像高度（象素数）

BOOL WINAPI DIBWalsh1(LPSTR lpDIBBits, LONG lWidth, LONG lHeight)
{
// 指向源图像的指针
unsigned char* lpSrc;
// 循环变量
LONG i;
LONG j;
// 进行付立叶变换的宽度和高度（2的整数次方）
LONG w;
LONG h;
// 中间变量
double dTemp;
int  wp;
int  hp;
// 图像每行的字节数
LONG lLineBytes;
// 计算图像每行的字节数
lLineBytes = WIDTHBYTES(lWidth * 8);

// 赋初值
w = 1;
h = 1;
wp = 0;
hp = 0;

// 计算进行离散余弦变换的宽度和高度（2的整数次方）
while(w * 2 <= lWidth)
{
w *= 2;
wp++;
}

while(h * 2 <= lHeight)
{
h *= 2;
hp++;
}

// 分配内存
double *f = new double[w * h];
double *F = new double[w * h];

// 列
for(i = 0; i < w; i++)
{
// 行
for(j = 0; j < h; j++)
{
// 指向DIB第j行，第i个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - j) + i;

// 给时域赋值
f[j + i * w] = *(lpSrc);
}
}

// 调用快速沃尔什－哈达玛变换
WALSH(f, F, wp + hp);
// 列
for(i = 0; i < w; i++)
{
// 行
for(j = 0; j < h; j++)
{
// 计算频谱
dTemp = fabs(F[i * w + j] * 1000);

// 判断是否超过255
if (dTemp > 255)
{
// 对于超过的，直接设置为255
dTemp = 255;
}
// 指向DIB第j行，第i个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - j) + i;

// 更新源图像
* (lpSrc) = (BYTE)(dTemp);
}
}

//释放内存
delete f;
delete F;

// 返回
return TRUE;
}

算法描述：
二值化是图像分割的一种方法。在二值化图象的时候把大于某个临界灰度值的像素灰度设为灰度極大值，把小于这个值的像素灰度设为灰度極小值，从而实现二值化。
根据阈值选取的不同，二值化的算法分为固定阈值和自适应阈值。 比较常用的二值化方法则有：双峰法、P参数法、迭代法和OTSU法等。

程序实现：

void CMyDIPView::OnDraw(CDC* pDC)
{
CMyDIPDoc* pDoc = GetDocument();
ASSERT_VALID(pDoc);
if(pDoc->m_hDIB == NULL)
return ;
// TODO: add draw code for native data here
int i,j;
unsigned char *lpSrc;
LPSTR lpDIB = (LPSTR) ::GlobalLock((HGLOBAL) pDoc->m_hDIB);
int cxDIB = (int) ::DIBWidth(lpDIB);         // Size of DIB - x
int cyDIB = (int) ::DIBHeight(lpDIB);        // Size of DIB - y
LPSTR lpDIBBits=::FindDIBBits (lpDIB);
// 计算图像每行的字节数
long lLineBytes = WIDTHBYTES(cxDIB * 8);
// 每行
for(i = 0; i < cyDIB; i++)
{
// 每列
for(j = 0; j < cxDIB; j++)
{
// 指向DIB第i行，第j个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (cyDIB - 1 - i) + j;
// 计算新的灰度值
//*(lpSrc) = BYTE(255-*lpSrc);
}
}
::GlobalUnlock((HGLOBAL) pDoc->m_hDIB);
CRect rect(0,0,cxDIB,cyDIB), rcDIB(0,0,cxDIB,cyDIB);
::PaintDIB(pDC->m_hDC, &rect, pDoc->m_hDIB, &rcDIB, pDoc->m_palDIB);
}

{
int i,j;
unsigned char *lpSrc;
CMyDIPDoc* pDoc = GetDocument();
ASSERT_VALID(pDoc);
if(pDoc->m_hDIB == NULL)
return ;
LPSTR lpDIB = (LPSTR) ::GlobalLock((HGLOBAL) pDoc->m_hDIB);
LPSTR lpDIBBits=::FindDIBBits (lpDIB);
int cxDIB = (int) ::DIBWidth(lpDIB);         // Size of DIB - x
int cyDIB = (int) ::DIBHeight(lpDIB);        // Size of DIB - y
long lLineBytes = WIDTHBYTES(cxDIB * 8);     // 计算图像每行的字节数
const float c1=150,c2=2.5;
// 每行
for(i = 0; i < cyDIB; i++)
{
// 每列
for(j = 0; j < cxDIB; j++)
{
// 指向DIB第i行，第j个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (cyDIB - 1 - i) + j;

// 计算新的灰度值
if(*lpSrc<122) *lpSrc=BYTE(0);
else *lpSrc = BYTE(255);
}
}
::GlobalUnlock((HGLOBAL) pDoc->m_hDIB);
Invalidate(TRUE);
}

算法描述：
输入图像像元密度值(灰度、亮度值)按对数函数关系变换为输出图像。

程序实现：

//参数说明:
//LPSTR lpDIBBits:指向源DIB图像指针
//LONG  lWidth:源图像宽度（象素数）
//LONG  lHeight:源图像高度（象素数）
//BYTE  bThre:阈值
//程序说明:
//该函数用来对图像进行阈值变换。对于灰度值小于阈值的象素直接设置

BOOL WINAPI ThresholdTrans(LPSTR lpDIBBits, LONG lWidth, LONG lHeight, BYTE bThre)
{
// 指向源图像的指针
unsigned char* lpSrc;
// 循环变量
LONG i;
LONG j;
// 图像每行的字节数
LONG lLineBytes;
// 计算图像每行的字节数
lLineBytes = WIDTHBYTES(lWidth * 8);
// 每行
for(i = 0; i < lHeight; i++)
{
// 每列
for(j = 0; j < lWidth; j++)
{
// 指向DIB第i行，第j个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
// 判断是否小于阈值
if ((*lpSrc) < bThre)
{
// 直接赋值为0
*lpSrc = 0;
}
else
{
// 直接赋值为255
*lpSrc = 255;
}
}
}
// 返回
return TRUE;
}

算法描述：
关于此傅里叶变换算法的具体介绍，请参考本BLOG文章：十、从头到尾彻底理解傅里叶变换算法、上

程序实现：

complex<double> * TD - 指向时域数组的指针
complex<double> * FD - 指向频域数组的指针
r      －2的幂数，即迭代次数

VOID WINAPI FFT(complex<double> * TD, complex<double> * FD, int r)
{
// 付立叶变换点数
LONG count;
// 循环变量
int  i,j,k;
// 中间变量
int  bfsize,p;
// 角度
double angle;
complex<double> *W,*X1,*X2,*X;
// 计算付立叶变换点数
count = 1 << r;

// 分配运算所需存储器
W  = new complex<double>[count / 2];
X1 = new complex<double>[count];
X2 = new complex<double>[count];

// 计算加权系数
for(i = 0; i < count / 2; i++)
{
angle = -i * PI * 2 / count;
W[i] = complex<double> (cos(angle), sin(angle));
}

// 将时域点写入X1
memcpy(X1, TD, sizeof(complex<double>) * count);

// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (r-k);
for(i = 0; i < bfsize / 2; i++)
{
p = j * bfsize;
X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)];
}
}
X  = X1;
X1 = X2;
X2 = X;
}

// 重新排序
for(j = 0; j < count; j++)
{
p = 0;
for(i = 0; i < r; i++)
{
if (j&(1<<i))
{
p+=1<<(r-i-1);
}
}
FD[j]=X1[p];
}

// 释放内存
delete W;
delete X1;
delete X2;
}

LPSTR lpDIBBits    - 指向源DIB图像指针
LONG  lWidth       - 源图像宽度（象素数）
LONG  lHeight      - 源图像高度（象素数）

BOOL WINAPI Fourier(LPSTR lpDIBBits, LONG lWidth, LONG lHeight)
{
// 指向源图像的指针
unsigned char* lpSrc;
// 中间变量
double dTemp;
// 循环变量
LONG i;
LONG j;

// 进行付立叶变换的宽度和高度（2的整数次方）
LONG w;
LONG h;
int  wp;
int  hp;

// 图像每行的字节数
LONG lLineBytes;

// 计算图像每行的字节数
lLineBytes = WIDTHBYTES(lWidth * 8);

// 赋初值
w = 1;
h = 1;
wp = 0;
hp = 0;

// 计算进行付立叶变换的宽度和高度（2的整数次方）
while(w * 2 <= lWidth)
{
w *= 2;
wp++;
}

while(h * 2 <= lHeight)
{
h *= 2;
hp++;
}

// 分配内存
complex<double> *TD = new complex<double>[w * h];
complex<double> *FD = new complex<double>[w * h];

// 行
for(i = 0; i < h; i++)
{
// 列
for(j = 0; j < w; j++)
{
// 指向DIB第i行，第j个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;

// 给时域赋值
TD[j + w * i] = complex<double>(*(lpSrc), 0);
}
}

for(i = 0; i < h; i++)
{
// 对y方向进行快速付立叶变换
FFT(&TD[w * i], &FD[w * i], wp);
}

// 保存变换结果
for(i = 0; i < h; i++)
{
for(j = 0; j < w; j++)
{
TD[i + h * j] = FD[j + w * i];
}
}

for(i = 0; i < w; i++)
{
// 对x方向进行快速付立叶变换
FFT(&TD[i * h], &FD[i * h], hp);
}

// 行
for(i = 0; i < h; i++)
{
// 列
for(j = 0; j < w; j++)
{
// 计算频谱
dTemp = sqrt(FD[j * h + i].real() * FD[j * h + i].real() +
FD[j * h + i].imag() * FD[j * h + i].imag()) / 100;
// 判断是否超过255
if (dTemp > 255)
{
// 对于超过的，直接设置为255
dTemp = 255;
}
// 指向DIB第(i<h/2 ? i+h/2 : i-h/2)行，第(j<w/2 ? j+w/2 : j-w/2)个象素的指针
// 此处不直接取i和j，是为了将变换后的原点移到中心
//lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
lpSrc = (unsigned char*)lpDIBBits + lLineBytes *
(lHeight - 1 - (i<h/2 ? i+h/2 : i-h/2)) + (j<w/2 ? j+w/2 : j-w/2);

// 更新源图像
* (lpSrc) = (BYTE)(dTemp);
}
}

// 删除临时变量
delete TD;
delete FD;

// 返回
return TRUE;
}

July附注：此傅里叶变换算法，在本BLOG内有深入具体的介绍，请参考本BLOG内其它文章。

算法描述：
离散余弦变换（DCT for Discrete Cosine Transform）是与傅里叶变换相关的一种变换，它类似于离散傅里叶变换（DFT for Discrete Fourier Transform），但是只使用实数。
离散余弦变换相当于一个长度大概是它两倍的离散傅里叶变换，这个离散傅里叶变换是对一个实偶函数进行的（因为一个实偶函数的傅里叶变换仍然是一个实偶函数），在有些变形里面需要将输入或者输出的位置移动半个单位（DCT有8种标准类型，其中4种是常见的）。

程序实现：

complex<double> * TD - 指向时域数组的指针
complex<double> * FD - 指向频域数组的指针
r         －2的幂数，即迭代次数

VOID WINAPI FFT(complex<double> * TD, complex<double> * FD, int r)
{
// 付立叶变换点数
LONG count;
// 循环变量
int  i,j,k;
// 中间变量
int  bfsize,p;
// 角度
double angle;

complex<double> *W,*X1,*X2,*X;
// 计算付立叶变换点数
count = 1 << r;

// 分配运算所需存储器
W  = new complex<double>[count / 2];
X1 = new complex<double>[count];
X2 = new complex<double>[count];

// 计算加权系数
for(i = 0; i < count / 2; i++)
{
angle = -i * PI * 2 / count;
W[i] = complex<double> (cos(angle), sin(angle));
}

// 将时域点写入X1
memcpy(X1, TD, sizeof(complex<double>) * count);

// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (r-k);
for(i = 0; i < bfsize / 2; i++)
{
p = j * bfsize;
X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)];
}
}
X  = X1;
X1 = X2;
X2 = X;
}

// 重新排序
for(j = 0; j < count; j++)
{
p = 0;
for(i = 0; i < r; i++)
{
if (j&(1<<i))
{
p+=1<<(r-i-1);
}
}
FD[j]=X1[p];
}

// 释放内存
delete W;
delete X1;
delete X2;
}

double * f    - 指向时域值的指针
double * F    - 指向频域值的指针
r      －2的幂数

VOID WINAPI DCT(double *f, double *F, int r)
{
// 离散余弦变换点数
LONG count;
// 循环变量
int  i;
// 中间变量
double dTemp;

complex<double> *X;
// 计算离散余弦变换点数
count = 1<<r;

// 分配内存
X = new complex<double>[count*2];
// 赋初值为0
memset(X, 0, sizeof(complex<double>) * count * 2);
// 将时域点写入数组X
for(i=0;i<count;i++)
{
X[i] = complex<double> (f[i], 0);
}

// 调用快速付立叶变换
FFT(X,X,r+1);
// 调整系数
dTemp = 1/sqrt(count);
// 求F[0]
F[0] = X[0].real() * dTemp;
dTemp *= sqrt(2);
// 求F[u]
for(i = 1; i < count; i++)
{
F[i]=(X[i].real() * cos(i*PI/(count*2)) + X[i].imag() * sin(i*PI/(count*2))) * dTemp;
}

// 释放内存
delete X;
}

LPSTR lpDIBBits    - 指向源DIB图像指针
LONG  lWidth       - 源图像宽度（象素数）
LONG  lHeight      - 源图像高度（象素数）

BOOL WINAPI DIBDct(LPSTR lpDIBBits, LONG lWidth, LONG lHeight)
{
// 指向源图像的指针
unsigned char* lpSrc;
// 循环变量
LONG i;
LONG j;
// 进行付立叶变换的宽度和高度（2的整数次方）
LONG w;
LONG h;
// 中间变量
double dTemp;
int  wp;
int  hp;

// 图像每行的字节数
LONG lLineBytes;
// 计算图像每行的字节数
lLineBytes = WIDTHBYTES(lWidth * 8);

// 赋初值
w = 1;
h = 1;
wp = 0;
hp = 0;

// 计算进行离散余弦变换的宽度和高度（2的整数次方）
while(w * 2 <= lWidth)
{
w *= 2;
wp++;
}

while(h * 2 <= lHeight)
{
h *= 2;
hp++;
}
// 分配内存
double *f = new double[w * h];
double *F = new double[w * h];

// 行
for(i = 0; i < h; i++)
{
// 列
for(j = 0; j < w; j++)
{
// 指向DIB第i行，第j个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;

// 给时域赋值
f[j + i * w] = *(lpSrc);
}
}

for(i = 0; i < h; i++)
{
// 对y方向进行离散余弦变换
DCT(&f[w * i], &F[w * i], wp);
}

// 保存计算结果
for(i = 0; i < h; i++)
{
for(j = 0; j < w; j++)
{
f[j * h + i] = F[j + w * i];
}
}

for(j = 0; j < w; j++)
{
// 对x方向进行离散余弦变换
DCT(&f[j * h], &F[j * h], hp);
}
// 行
for(i = 0; i < h; i++)
{
// 列
for(j = 0; j < w; j++)
{
// 计算频谱
dTemp = fabs(F[j*h+i]);

// 判断是否超过255
if (dTemp > 255)
{
// 对于超过的，直接设置为255
dTemp = 255;
}

// 指向DIB第y行，第x个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;

// 更新源图像
* (lpSrc) = (BYTE)(dTemp);
}
}

// 释放内存
delete f;
delete F;

// 返回
return TRUE;
}

变化效果：