【图像处理中的数学修炼】一书之代码(1/3)

数字图像处理对数学的要求颇高，这不禁令很多学习者望而却步。在阅读图像处理方面的论文时，面对梯度、散度、黑塞矩阵、傅里叶变换等这些本该在微积分中早已耳熟能详的概念时，很多人仍然感觉一筹莫展。为了弭平图像处理道路上的数学险阻，帮助更多人学好数字图像处理，并更快地具备深入研究的能力。笔者特别撰写了这本《图像处理中的数学修炼》（该书现已由清华大学出版社正式出版）。欲了解《图像处理中的数学修炼》的更多详细内容，你可以参考本书的目录。

通常，我不喜欢翻开一本技术书，里面满篇满篇的都是代码。我也不希望用代码去凑页数或者挤占宝贵的篇幅。就《图像处理中的数学修炼》一书而言，更是如此。本书的重点在于全面系统地对图像处理中可能用到的各种数学基础加以总结归纳，所以数学原理才是我们所关注的重点！

但是如果把数学原理同图像处理实践割裂开来，又可能令某些读者觉得“英雄无用武之地”！为了真正帮读者建立起从数学原理到图像处理实践的桥梁，本书特意在后半部分安排了一些比较实用且非常经典的图像处理算法介绍，而书籍的前半部分则是纯数学部分，后半部分中的经典算法大量使用了前半部分里的数学知识，我们希望以这种方式来帮助读者夯实基础、巩固所学。

在后半部分的图像处理算法介绍部分，偶尔为了便于理解，本书配备了必要的MATLAB代码。但所用到的代码也非常有限的。因为这并不是一本教你怎么使用MATLAB的书！全书的代码不过六百多行，主要是为了对一些比较经典但又比较复杂的算法加以解释。这些“比较经典但又比较复杂的算法”主要包括（但不限于）：

• 基于暗通道的图像去雾实现
• 对比度有限的自适应直方图均衡（CLAHE）
• 自适应直方图均衡（AHE）
• 基于小波变换的图像融合（失焦图像修复）
• 基于泊松方程的泊松融合算法（基于迭代）
• 基于泊松方程的泊松融合算法（有限差分方法解偏微分方程）
• 基于主成分提取（PCA）的图像压缩

特别感谢前期试读本书的读者给予笔者的意见反馈，以及CSND博客上【图像处理中的数学原理详解】专栏的读者提出的宝贵见解。下面所列之代码已经结合之前收到的意见反馈进行了调整和优化，并修正了此前存在的一些小问题。非常感谢这部分读者为提高本书质量所做之努力。此书文字部分的勘误可以参见《 图像处理中的数学修炼》一书之勘误表 。

P270

i=double(imread(‘vase.tif‘));
[C,S]=wavedec2(i,2,‘db1‘);
a2=appcoef2(C,S,‘db1‘,2);
dh1=detcoef2(‘h‘,C,S,1);
dv1=detcoef2(‘v‘,C,S,1);
dd1=detcoef2(‘d‘,C,S,1);
dh2=detcoef2(‘h‘,C,S,2);
dv2=detcoef2(‘v‘,C,S,2);
dd2=detcoef2(‘d‘,C,S,2);
[x,y]=size(i);
img = zeros(x,y);
img(1:x/4,1:y/4) =im2uint8(mat2gray(a2));
img(((x/4)+1):y/2,1:y/4) = im2uint8(mat2gray(dv2));
img(((x/4)+1):x/2,1:y/4) = im2uint8(mat2gray(dv2));
img(1:x/4,((y/4)+1):y/2) = im2uint8(mat2gray(dh2));
img(((x/4)+1):x/2,((y/4)+1):y/2) = im2uint8(mat2gray(dd2));
img(((x/2)+1):x,1:y/2) = im2uint8(mat2gray(dv1));
img(1:x/2,((y/2)+1):y) = im2uint8(mat2gray(dh1));
img(((x/2)+1):x,((y/2)+1):y) = im2uint8(mat2gray(dd1));
imshow(img,[]);

P272

X1 = imread(‘cathe1.bmp‘);
X2 = imread(‘cathe2.bmp‘);
XFUS = wfusimg(X1,X2,‘sym4‘,5,‘mean‘,‘max‘);
imshow(XFUS,[]);

P273

X1 = imread(‘cathe1.bmp‘);
X2 = imread(‘cathe2.bmp‘);
M1 = double(X1) / 256;
M2 = double(X2) / 256;
N = 4;
wtype = ‘sym4‘;
[c0,s0] = wavedec2(M1, N, wtype);
[c1,s1] = wavedec2(M2, N, wtype);
length = size(c1);
Coef_Fusion = zeros(1,length(2));
%低频系数的处理，取平均值
Coef_Fusion(1:s1(1,1)) = (c0(1:s1(1,1))+c1(1:s1(1,1)))/2;
%处理高频系数，取绝对值大者，这里用到了矩阵乘法
MM1 = c0(s1(1,1)+1:length(2));
MM2 = c1(s1(1,1)+1:length(2));
mm = (abs(MM1)) > (abs(MM2));
Y  = (mm.*MM1) + ((~mm).*MM2);
Coef_Fusion(s1(1,1)+1:length(2)) = Y;
%重构
Y = waverec2(Coef_Fusion,s0,wtype);
imshow(Y,[]);

P274

I = imread(‘noise_lena.bmp‘);
[thr,sorh,keepapp] = ddencmp(‘den‘,‘wv‘,I);
de_I = wdencmp(‘gbl‘,I,‘sym4‘,2,thr,sorh,keepapp);
imwrite(im2uint8(mat2gray(de_I)), ‘denoise_lena.bmp‘);

P298

I = imread(‘baboon.bmp‘);  
I1 = double(I);  
T = hadamard(8);  
myFun1 = @(block_struct)T*block_struct.data*T/64;  
H = blockproc(I1, [8 8], myFun1);  
H(abs(H)<3.5)=0;  
myFun2 = @(block_struct)T*block_struct.data*T;  
I2 = blockproc(H, [8 8], myFun2);  
subplot(121), imshow(I1,[]), title(‘original image‘);  
subplot(122), imshow(I2,[]), title(‘zipped image‘);  

P299

I = imread(‘baboon.bmp‘);  
I1 = double(I);  
[m n] =size(I);  
sizi = 8;  
num = 16;  
%分块进行离散沃尔什变换  
T = hadamard(sizi);  
myFun1 = @(block_struct)T*block_struct.data*T/(sizi.^2);  
hdcoe = blockproc(I1, [sizi, sizi], myFun1);  
%重新排列系数  
coe = im2col(hdcoe,  [sizi, sizi], ‘distinct‘);  
coe_t = abs(coe);  
[Y, ind] = sort(coe_t);  
%舍去绝对值较小的系数  
[m_c, n_c] = size(coe);  
for i = 1:n_c  
coe(ind(1:num, i), i)=0;  
end  
%重建图像  
re_hdcoe = col2im(coe, [sizi, sizi], [m, n], ‘distinct‘);  
myFun2 = @(block_struct)T*block_struct.data*T;  
re_s = blockproc(re_hdcoe, [sizi, sizi], myFun2);  
subplot(121), imshow(I1,[]), title(‘original image‘);  
subplot(122), imshow(re_s,[]), title(‘compressed image‘);  

P307

I = imread(‘baboon.bmp‘);  
x = double(I)/255;  
[m,n]=size(x);  
y =[];  
%拆解图像  
for i = 1:m/8;  
    for j = 1:n/8;  
        ii = (i-1)*8+1;  
        jj = (j-1)*8+1;  
        y_app = reshape(x(ii:ii+7,jj:jj+7),1,64);  
        y=[y;y_app];  
    end  
end  

%KL变换  
[COEFF,SCORE,latent] = princomp(y);  
kl = y * COEFF;  

kl1 = kl;  
kl2 = kl;  
kl3 = kl;  

%置零压缩过程  
kl1(:, 33:64)=0;  
kl2(:, 17:64)=0;  
kl3(:, 9:64)=0;  

%KL逆变换  
kl_i = kl*COEFF‘;  
kl1_i = kl1*COEFF‘;  
kl2_i = kl2*COEFF‘;  
kl3_i = kl3*COEFF‘;  

image = ones(256,256);  
image1 = ones(256,256);  
image2 = ones(256,256);  
image3 = ones(256,256);  

k=1;  
%重组图像  
for i = 1:m/8;  
    for j = 1:n/8;  

        y = reshape(kl_i(k, 1:64),8,8);  
        y1 = reshape(kl1_i(k, 1:64),8,8);  
        y2 = reshape(kl2_i(k, 1:64),8,8);  
        y3 = reshape(kl3_i(k, 1:64),8,8);  

        ii = (i-1)*8+1;  
        jj = (j-1)*8+1;  

        image(ii:ii+7,jj:jj+7) = y;  
        image1(ii:ii+7,jj:jj+7) = y1;  
        image2(ii:ii+7,jj:jj+7) = y2;  
        image3(ii:ii+7,jj:jj+7) = y3;  

        k=k+1;  
    end  
end  

P356-1

mountains = double(imread(‘./img/mountain.jpg‘));
moon = double(imread(‘./img/moon.png‘));
sizeSrc = size(mountains);
sizeDst = size(moon);

gradient_inner = moon(1:sizeDst(1)-2,2:sizeDst(2)-1,:)...
    + moon(3:sizeDst(1),2:sizeDst(2)-1,:)...
    + moon(2:sizeDst(1)-1,1:sizeDst(2)-2,:)...
    + moon(2:sizeDst(1)-1,3:sizeDst(2),:)...
    - 4*moon(2:sizeDst(1)-1,2:sizeDst(2)-1,:);

P356-2

Lap = [0, 1, 0;1, -4, 1;0, 1, 0];

I1 = conv2(double(moon(:,:,1)), double(Lap));
I2 = conv2(double(moon(:,:,2)), double(Lap));
I3 = conv2(double(moon(:,:,3)), double(Lap));
gradient_inner(:, :, 1) = I1(3:sizeDst(1),3:sizeDst(2));
gradient_inner(:, :, 2) = I2(3:sizeDst(1),3:sizeDst(2));
gradient_inner(:, :, 3) = I3(3:sizeDst(1),3:sizeDst(2));

P357

dstX = 350;dstY = 100;
rebuilt = mountains(dstY:dstY+sizeDst(1)-1,dstX:dstX+sizeDst(2)-1,:);

for n = [1:1000]
    rebuilt(2:2:sizeDst(1)-1,2:2:sizeDst(2)-1,:)= ...
        (rebuilt(1:2:sizeDst(1)-2 , 2:2:sizeDst(2)-1,:)...
        +rebuilt(3:2:sizeDst(1) , 2:2:sizeDst(2)-1,:)...
        +rebuilt(2:2:sizeDst(1)-1 , 1:2:sizeDst(2)-2,:)...
        +rebuilt(2:2:sizeDst(1)-1 , 3:2:sizeDst(2),:)...
        -gradient_inner(1:2:sizeDst(1)-2 , 1:2:sizeDst(2)-2,:))/4;
     rebuilt(3:2:sizeDst(1)-1,3:2:sizeDst(2)-1,:)= ...
        (rebuilt(2:2:sizeDst(1)-2 , 3:2:sizeDst(2)-1,:)...
        +rebuilt(4:2:sizeDst(1) , 3:2:sizeDst(2)-1,:)...
        +rebuilt(3:2:sizeDst(1)-1 , 2:2:sizeDst(2)-2,:)...
        +rebuilt(3:2:sizeDst(1)-1 , 4:2:sizeDst(2),:)...
        -gradient_inner(2:2:sizeDst(1)-2 , 2:2:sizeDst(2)-2,:))/4;
     rebuilt(3:2:sizeDst(1)-1,2:2:sizeDst(2)-1,:)= ...
        (rebuilt(2:2:sizeDst(1)-2 , 2:2:sizeDst(2)-1,:)...
        +rebuilt(4:2:sizeDst(1) , 2:2:sizeDst(2)-1,:)...
        +rebuilt(3:2:sizeDst(1)-1 , 1:2:sizeDst(2)-2,:)...
        +rebuilt(3:2:sizeDst(1)-1 , 3:2:sizeDst(2),:)...
        -gradient_inner(2:2:sizeDst(1)-2 , 1:2:sizeDst(2)-2,:))/4;
    rebuilt(2:2:sizeDst(1)-1 , 3:2:sizeDst(2)-1,:)= ...
        (rebuilt(1:2:sizeDst(1)-2 , 3:2:sizeDst(2)-1,:)...
        +rebuilt(3:2:sizeDst(1) , 3:2:sizeDst(2)-1,:)...
        +rebuilt(2:2:sizeDst(1)-1 , 2:2:sizeDst(2)-2,:)...
        +rebuilt(2:2:sizeDst(1)-1 , 4:2:sizeDst(2),:)...
        -gradient_inner(1:2:sizeDst(1)-2 , 2:2:sizeDst(2)-2,:))/4;
end

mountains(dstY:sizeDst(1)+dstY-1,dstX:sizeDst(2)+dstX-1,:) = rebuilt;
figure,imshow(uint8(mountains));

本书中用于实验的图片下载【百度网盘】

如果你是图像处理的同道中人，欢迎加入图像处理算法交流群（单击链接查看群号），本书中的补充代码亦可在此群中获取，如果你在阅读本书时有任何问题也可以在QQ群里同笔者讨论。