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.imread('test.jpeg')
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)通过实验了解二维频谱的分布特点。

4)了解余弦变换或WalshHadamard变换

二、实验内容

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

   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

简介:

阿达马变换Hadamard transform),或称沃尔什-阿达玛转换,是一种广义傅立叶变换(Fourier transforms),作为变换编码的一种在视频编码当中使用有很久的历史。在近来的视频编码标准中,阿达马变换多被用来计算SATD(一种视频残差信号大小的衡量)。

数位信号处理大型集成电路算法的领域中,阿达马变换是一种简单且重要的算法之一,主要能针对频谱做快速的分析。

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



SATD计算方法

当计算4x4块\begin{bmatrix}L_4\end{bmatrix}的SATD时,先使用下面的方法进行二维的阿达马变换:


然后计算\begin{bmatrix}L_4'\end{bmatrix}所有系数绝对值之和并归一化

类似的,当计算8x8块\begin{bmatrix}L_8\end{bmatrix}的SATD时,先使用下面的方法进行二维的Hadamard变换:


然后计算\begin{bmatrix}L_8'\end{bmatrix}所有系数绝对值之和并归一化

下面是实现二维快速Walsh-Hadamard 变换matlab代码:

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变换在图像压缩领域的应用:
Walsh-Hadamard变换是实时的、对称的正交变换。因为该变换只包括加、减,没有任何乘、除运算,所以有快速算法。Walsh-Hadamard变换在图像处理中

的主要 用途是压缩编码。下面用一个图像压缩的例子讲述matlab如何实现Walsh-Hadamard变换。

将图像分割成为16*16的子图像后,每个子图像经过变换将有很多系数。按照每个系数的方差来排序。保留方差较大的系数,舍弃方差较小的系数。程序如下:

I=imread('lena.png');
sig=rgb2gray(I);
sig=double(sig)/255;
%figure,imshow(sig);
[m_sig,n_sig]=size(sig);
sizi=16;
Snum=128;
T=hadamard(sizi);
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
  1. /*************************************************************************
  2. *
  3. * 函数名称:
  4. * WALSH()
  5. *
  6. * 参数:
  7. * double* f - 输入的时域序列
  8. * double* F - 输出的频域序列
  9. * int r - 2的幂数
  10. *
  11. * 返回值:
  12. * BOOL - 成功返回TRUE,否则返回FALSE。
  13. *
  14. * 说明:
  15. * 该函数进行一维快速沃尔什——哈达马变换。
  16. *
  17. ************************************************************************/
  18. void WALSH(double*f,double*F,int r);
  19. /*************************************************************************
  20. *
  21. * 函数名称:
  22. * IWALSH()
  23. *
  24. * 参数:
  25. * double* f - 输入的时域序列
  26. * double* F - 输出的频域序列
  27. * int r - 2的幂数
  28. *
  29. * 返回值:
  30. * BOOL - 成功返回TRUE,否则返回FALSE。
  31. *
  32. * 说明:
  33. * 该函数进行一维快速沃尔什——哈达马逆变换。
  34. *
  35. ************************************************************************/
  36. void IWALSH(double*F,double*f,int r);
  37. /*************************************************************************
  38. *
  39. * 函数名称:
  40. * FreqWALSH()
  41. *
  42. * 参数:
  43. * double* f - 输入的时域序列
  44. * double* F - 输出的频域序列
  45. * LONG width - 图象宽度
  46. * LONG height - 图象高度
  47. *
  48. * 返回值:
  49. * BOOL - 成功返回TRUE,否则返回FALSE。
  50. *
  51. * 说明:
  52. * 该函数进行二维快速沃尔什——哈达玛变换。
  53. *
  54. ************************************************************************/
  55. BOOL FreqWALSH(double*f,double*F, LONG width, LONG height);
  56. /*************************************************************************
  57. *
  58. * 函数名称:
  59. * IFreqWALSH()
  60. *
  61. * 参数:
  62. * double* f - 输入的时域序列
  63. * double* F - 输出的频域序列
  64. * LONG width - 图象宽度
  65. * LONG height - 图象高度
  66. *
  67. * 返回值:
  68. * BOOL - 成功返回TRUE,否则返回FALSE。
  69. *
  70. * 说明:
  71. * 该函数进行二维快速沃尔什——哈达玛逆变换。
  72. *
  73. ************************************************************************/
  74. BOOL IFreqWALSH(double*f,double*F, LONG width, LONG height);
  75. /*************************************************************************
  76. *
  77. * 函数名称:
  78. * BmpWalsh()
  79. *
  80. * 参数:
  81. * BYTE* bmp,LONG width,LONG height
  82. *
  83. * 返回值:
  84. * BOOL - 成功返回TRUE,否则返回FALSE。
  85. *
  86. * 说明:
  87. * 该函数对图象进行二维快速沃尔什——哈达马变换。
  88. *
  89. ************************************************************************/
  90. BOOL BmpWalsh(BYTE* bmp,LONG width,LONG height);
  1. voidMyProcess::WALSH(double*f,double*F,int r)
  2. {
  3. // 循环变量
  4. LONG i;
  5. LONG j;
  6. LONG k;
  7. // 中间变量
  8. int p;
  9. double* X;
  10. // 计算快速沃尔什变换点数
  11. LONG N =1<< r;
  12. // 分配运算所需的数组
  13. double* X1 =newdouble[N];
  14. double* X2 =newdouble[N];
  15. // 将时域点写入数组X1
  16. memcpy(X1, f,sizeof(double)* N);
  17. // 蝶形运算
  18. for(k =0; k < r; k++)
  19. {
  20. for(j =0; j <1<<k; j++)
  21. {
  22. for(i =0; i <1<<(r - k -1); i++)
  23. {
  24. p = j *(1<<(r - k));
  25. X2[i + p]= X1[i + p]+ X1[i + p +(int)(1<<(r - k -1))];
  26. X2[i + p +(int)(1<<(r - k -1))]= X1[i + p]- X1[i + p +(int)(1<<(r - k -1))];
  27. }
  28. }
  29. // 互换X1和X2
  30. X = X1;
  31. X1 = X2;
  32. X2 = X;
  33. }
  34. // 调整系数
  35. for(j =0; j < N; j++)
  36. {
  37. p =0;
  38. for(i =0; i < r; i++)
  39. {
  40. if(j &(1<<i))
  41. {
  42. p +=1<<(r - i -1);
  43. }
  44. }
  45. F[j]= X1[p]/ N;
  46. }
  47. // 释放内存
  48. delete X1;
  49. delete X2;
  50. }
  51. voidMyProcess::IWALSH(double*F,double*f,int r)
  52. {
  53. // 循环变量
  54. int i;
  55. // 计算变换点数
  56. LONG N =1<< r;
  57. // 调用快速沃尔什-哈达玛变换进行反变换
  58. WALSH(F, f, r);
  59. // 调整系数
  60. for(i =0; i < N; i++)
  61. f[i]*= N;
  62. }
  63. BOOL MyProcess::FreqWALSH(double*f,double*F, LONG width, LONG height)
  64. {
  65. // 循环变量
  66. LONG i;
  67. LONG j;
  68. LONG k;
  69. // 进行离散余弦变换的宽度和高度(2的整数次方)
  70. LONG w =1;
  71. LONG h =1;
  72. int wp =0;
  73. int hp =0;
  74. // 计算进行离散余弦变换的宽度和高度(2的整数次方)
  75. while(w < width/3)
  76. {
  77. w *=2;
  78. wp++;
  79. }
  80. while(h < height)
  81. {
  82. h *=2;
  83. hp++;
  84. }
  85. // 分配内存
  86. double*TempIn=newdouble[h];
  87. double*TempOut=newdouble[h];
  88. // 对y方向进行离散余弦变换
  89. for(i =0; i < w *3; i++)
  90. {
  91. // 抽取数据
  92. for(j =0; j < h; j++)
  93. TempIn[j]= f[j * w *3+ i];
  94. // 一维快速离散余弦变换
  95. WALSH(TempIn,TempOut, hp);
  96. // 保存变换结果
  97. for(j =0; j < h; j++)
  98. f[j * w *3+ i]=TempOut[j];
  99. }
  100. // 释放内存
  101. deleteTempIn;
  102. deleteTempOut;
  103. // 分配内存
  104. TempIn=newdouble[w];
  105. TempOut=newdouble[w];
  106. // 对x方向进行快速离散余弦变换
  107. for(i =0; i < h; i++)
  108. {
  109. for(k =0; k <3; k++)
  110. {
  111. // 抽取数据
  112. for(j =0; j < w; j++)
  113. TempIn[j]= f[i * w *3+ j *3+ k];
  114. // 一维快速离散余弦变换
  115. WALSH(TempIn,TempOut, wp);
  116. // 保存变换结果
  117. for(j =0; j < w; j++)
  118. F[i * w *3+ j *3+ k]=TempOut[j];
  119. }
  120. }
  121. // 释放内存
  122. deleteTempIn;
  123. deleteTempOut;
  124. return TRUE;
  125. }
  126. BOOL MyProcess::IFreqWALSH(double*f,double*F, LONG width, LONG height)
  127. {
  128. // 循环变量
  129. LONG i;
  130. LONG j;
  131. LONG k;
  132. // 赋初值
  133. LONG w =1;
  134. LONG h =1;
  135. int wp =0;
  136. int hp =0;
  137. // 计算进行付立叶变换的宽度和高度(2的整数次方)
  138. while(w < width/3)
  139. {
  140. w *=2;
  141. wp++;
  142. }
  143. while(h < height)
  144. {
  145. h *=2;
  146. hp++;
  147. }
  148. // 分配内存
  149. double*TempIn=newdouble[w];
  150. double*TempOut=newdouble[w];
  151. // 对x方向进行快速付立叶变换
  152. for(i =0; i < h; i++)
  153. {
  154. for(k =0; k <3; k++)
  155. {
  156. // 抽取数据
  157. for(j =0; j < w; j++)
  158. TempIn[j]= F[i * w *3+ j *3+ k];
  159. // 一维快速傅立叶变换
  160. IWALSH(TempIn,TempOut, wp);
  161. // 保存变换结果
  162. for(j =0; j < w; j++)
  163. F[i * w *3+ j *3+ k]=TempOut[j];
  164. }
  165. }
  166. // 释放内存
  167. deleteTempIn;
  168. deleteTempOut;
  169. TempIn=newdouble[h];
  170. TempOut=newdouble[h];
  171. // 对y方向进行快速付立叶变换
  172. for(i =0; i < w *3; i++)
  173. {
  174. // 抽取数据
  175. for(j =0; j < h; j++)
  176. TempIn[j]= F[j * w *3+ i];
  177. // 一维快速傅立叶变换
  178. IWALSH(TempIn,TempOut, hp);
  179. // 保存变换结果
  180. for(j =0; j < h; j++)
  181. F[j * w *3+ i]=TempOut[j];
  182. }
  183. // 释放内存
  184. deleteTempIn;
  185. deleteTempOut;
  186. for(i =0; i < h; i++)
  187. {
  188. for(j =0; j < w *3; j++)
  189. {
  190. if(i < height && j < width)
  191. *(f + i * width + j)= F[i * w *3+ j];
  192. }
  193. }
  194. return TRUE;
  195. }
  196. BOOL MyProcess::BmpWalsh(BYTE* bmp,LONG width,LONG height)
  197. {
  198. // 循环变量
  199. LONG i;
  200. LONG j;
  201. // 进行沃尔什——哈达玛变换的宽度和高度(2的整数次方)
  202. LONG w =1;
  203. LONG h =1;
  204. int wp =0;
  205. int hp =0;
  206. // 计算进行离散余弦变换的宽度和高度(2的整数次方)
  207. while(w < width/3)
  208. {
  209. w *=2;
  210. wp++;
  211. }
  212. while(h < height)
  213. {
  214. h *=2;
  215. hp++;
  216. }
  217. // 分配内存
  218. double*f =newdouble[w * h *3];
  219. double*F =newdouble[w * h *3];
  220. // 向时域赋值并补零
  221. for(i =0; i < h; i++)
  222. {
  223. for(j =0; j < w *3; j++)
  224. {
  225. if(i < height && j < width)
  226. f[i * w *3+ j]= bmp[width * i + j];
  227. else
  228. f[w * i *3+ j]=0.0f;
  229. }
  230. }
  231. // 进行频谱分析
  232. if(FreqWALSH(f, F,width, height)== FALSE)
  233. return FALSE;
  234. // 更新所有象素
  235. for(i =0; i < height; i++)
  236. {
  237. for(j =0; j < width; j++)
  238. {
  239. // 判断是否超过255
  240. if(fabs(F[i * w *3+ j]*1000)>255)
  241. {
  242. // 对于超过的,直接设置为255
  243. bmp[width *(height -1- i)+ j]=255;
  244. }
  245. else
  246. {
  247. // 如果没有超过,则按实际计算结果赋值
  248. bmp[width *(height -1- i)+ j]= fabs(F[i * w *3+ j]*1000);
  249. }
  250. }
  251. }
  252. //释放内存
  253. delete[] f;
  254. delete[] F;
  255. // 返回
  256. return TRUE;
  257. }

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

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

参考:百度百科、维基百科、vc数字图像处理。
--------------------------------------------------
数字图像处理领域的二十四个典型算法及vc实现、第一章
一、256色转灰度图
二、Walsh变换
三、二值化变换
四、阈值变换
五、傅立叶变换
六、离散余弦变换

数字图像处理领域的二十四个典型算法及vc实现、第二章
七、高斯平滑
八、图像平移
九、图像缩放
十、图像旋转
数字图像处理领域的二十四个典型算法及vc实现、第三章

 

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

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

      ok,请细看。

一、256色转灰度图
    算法介绍(
百度百科)
    什么叫灰度图?任何颜色都有红、绿、蓝三原色组成,假如原来某点的颜色为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)就是灰度图了。

灰度分为256阶。所以,用灰度表示的图像称作灰度图。

    程序实现:
    ok,知道了什么叫灰度图,下面,咱们就来实现此256色灰度图。
这个Convert256toGray(),即是将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变换
    算法介绍:
    有关Walsh变换的深入介绍,请看此论文:http://www.informatics.org.cn/doc/ucit200510/ucit20051005.pdf

    程序实现:

函数名称:WALSH()
参数:
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;
}

函数名称:DIBWalsh1()
参数:
LPSTR lpDIBBits    - 指向源DIB图像指针
LONG  lWidth       - 源图像宽度(象素数)
LONG  lHeight      - 源图像高度(象素数)
返回值:BOOL               - 成功返回TRUE,否则返回FALSE。
说明:该函数用来对图像进行沃尔什-哈达玛变换。于上面不同的是,此处是将二维
矩阵转换成一个列向量,然后对该列向量进行一次一维沃尔什-哈达玛变换。

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);
}

void CMyDIPView::OnMenuitem32778() 
{
 // TODO: Add your command handler code here
 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:阈值
//程序说明:
//该函数用来对图像进行阈值变换。对于灰度值小于阈值的象素直接设置
灰度值为0;灰度值大于阈值的象素直接设置为255。
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文章:十、从头到尾彻底理解傅里叶变换算法、上

    程序实现:

函数名称:FFT()
参数:
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;
}

函数名称:Fourier()
参数:
LPSTR lpDIBBits    - 指向源DIB图像指针
LONG  lWidth       - 源图像宽度(象素数)
LONG  lHeight      - 源图像高度(象素数)
返回值:BOOL               - 成功返回TRUE,否则返回FALSE。
说明:该函数用来对图像进行付立叶变换。
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种是常见的)。

    程序实现:

函数名称:FFT()
参数:
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;
}

函数名称:DCT()
参数:
double * f    - 指向时域值的指针
double * F    - 指向频域值的指针
r      -2的幂数
返回值:无。
说明:该函数用来实现快速离散余弦变换,利用2N点的快速付立叶变换来实现离散余弦变换。
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;
}

函数名称:DIBDct()
参数:
LPSTR lpDIBBits    - 指向源DIB图像指针
LONG  lWidth       - 源图像宽度(象素数)
LONG  lHeight      - 源图像高度(象素数)
返回值:BOOL               - 成功返回TRUE,否则返回FALSE。
说明:该函数用来对图像进行离散余弦变换。
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;
}


    变化效果:

没有更多推荐了,返回首页