1701010226
我和数学有约：趣味数学及算法解析 [:1701004270]
1701010227
我和数学有约：趣味数学及算法解析 9.3 还原真实场景
1701010228
1701010229
生活中常常遇到的图像，如图9-2和图9-3所示。
1701010230
1701010231
1701010232
1701010233
图9-2 部分清晰 图9-3 部分模糊 如图9-2所示，该图形中一部分是清晰的，一部分是模糊的。如图9-3所示，也是一部分模糊，一部分清晰，然而图9-2和图9-3两幅图恰好是互补的，两幅图恰好构成一副清晰的彩色图，然而对于这样的一种工况，在我们的相机拍照过程中经常出现，最常见的就是模糊了背景，突出了照相视觉中的人脸部分，然而面对这两幅图像，怎么完好的再现真实的场景呢？
1701010234
1701010235
【问题】有时候我们采集的图像常常带有模糊的成分，如何还原这些图像，重构出我们视觉看到的图形环境，显得至关重要。
1701010236
1701010237
【分析】
1701010238
1701010239
两幅互补的图像合成为一副清晰的彩色图，这是一种图像融合过程。
1701010240
1701010241
图像融合是指将来自不同探测器的多特征图像进行合并，去除图像间的冗余信息，保留不同图像的有用信息，从而得到一个包含被测物体信息更完整的图像。定义中的多特征图像是指各种探测器根据不同的成像机理来测量待测物体的不同物理特性而得到的图像，并以此作为融合处理的信息源。
1701010242
1701010243
图像融合在图像处理领域中属于图像增强的研究范畴。图像增强是对图像质量在一般意义上的改善，其目的是通过一系列技术去改善图像的视觉效果或将图像转换成一种更适合于人眼观察或机器自动分析的形式。图像增强技术是改善图像质量所采用的一种较为通用的方法。例如，在遥感应用中，常用两种传感器获得同一地域的高分辨图像和多光谱图像，融合这两种图像，可以获得细节和轮廓都得到改善的图像网。又例如，将融合技术用于同一数码相机在同一场景拍摄的聚焦点不同的两幅图像，可以获得比原图像清晰度高的新图像。但图像增强方法还缺乏统一的理论，因而一些增强方法往往带有针对性，以至于对某类图像效果较好的增强方法未必一定适用于另一类图像。
1701010244
1701010245
在实际应用中，要找出一种有效的方法常常必须经过广泛的实验和分析。图像增强技术也经常用于改善图像中某种特征的可探测性，从而使人或计算机更容易识别和分析。例如，减少图像模糊会使观察某种感兴趣的特征变得更容易些。有时可能需要彻底改变图像的视觉效果，以便突出重要特征的可观察性。在这种情况下，可以把图像增强理解为增强感兴趣特征的可识别性，而非改善视觉感官质量。
1701010246
1701010247
图像融合是20世纪70年代后期提出的新概念，是多传感器信息中可视信息部分的融合，是将多源信道所采集的关于同一目标图像经过一定的图像处理，提取各自信道的信息，最后综合成统一图像或综合图像特性以供观察或进一步处理。它是一门综合了传感器、图像处理、信号处理、显示、计算机和人工智能等技术的现代高新技术。由于图像融合系统具有突出的探测优越性（时空覆盖宽、目标分辨力与测量维数高、重构能力好、兀余性、互补性、时间优越性及相对低成本性等），在技术先进国家受到高度重视并已取得相当的进展。
1701010248
1701010249
近年来，随着小波理论的发展和应用，小波变换在时、频两域突出信号局部特征的能力己被越来越多的人所认可，各种类型的小波变换也相应地产生并得到了极为广泛的应用和发展。应用正交小波融合多传感器图像，在融合过程中，小波变换显示了它在压缩、方向信息上的正交性和有效性。
1701010250
1701010251
若对二维图像进行N层的小波分解，最终将有（3N＋1）个不同频带，其中包含3N个高频带和一个低频带，这样对不同分解层、不同频带就可以采用不同的融合规则分别进行融合处理。融合规则是图像融合的另一个重要方面，对融合图像起着同图像分解一样的作用，人们对融合规则的研究也在不断的升温，并已取得了很大的进展。
1701010252
1701010253
以前，人们提出在变换域的基本融合规则是绝对值最大选取，即保留带通中较大的绝对值用于重构，因为较大值与图像中的边缘、线性和区域边界相对应。Burt和Kolcz0ski用匹配和突出测量方法改进了这个规则，合成操作、选择或平均都是在匹配测量方法基础上实现的。
1701010254
1701010255
目前来看，融合规则有两大类，即基于像素的融合规则和基于区域（窗口）的融合规则，其中后一种融合规则对要融合的图像的要求较之前一种融合规则宽松，这对于图像融合的应用是非常有利的。
1701010256
1701010257
小波分析与Fourier变换相比，小波变换是空间（时间）和频率的局部变换，因而能有效地从信号中提取信息。通过伸缩和平移等运算功能可对函数或信号进行多尺度的细化分析，解决了Fourier变换不能解决的许多困难问题。
1701010258
1701010259
小波分析是一个新的数学分支，它是泛函分析、Fourier分析、样调分析和数值分析的完美结晶；小波分析是时间—尺度分析和多分辨分析的一种新技术，它在信号分析、语音合成、图像识别、计算机视觉、数据压缩、地震勘探、大气与海洋波分析等方面的研究都取得了有科学意义和应用价值的成果。
1701010260
1701010261
波分析的应用领域十分广泛。例如，在数学方面，它已用于数值分析、构造快速数值方法、曲线曲面构造、微分方程求解和控制论等。在信号分析方面的滤波、去噪声、压缩和传递等。在图像处理方面的图像压缩、分类、识别与诊断和去污等。在医学成像方面的减少B超、CT、核磁共振成像的时间和提高分辨率等。
1701010262
1701010263
对于如图9-2和图9-3所示的模糊图像融合，首先对彩色图像变成YIQ空间，对所有图像的各个颜色分量进行小波分解，然后通过小波反变换和YIQ反变换即得结果图像。
1701010264
1701010265
编写MATLAB程序如下：
1701010266
1701010267
clc,clear,close all %清屏和清除变量 warning off %消除警告 count=2; %融合图像数 level=2; %小波分解的层数 colors=3; %彩色图像的通道数(黑白图像为1，彩色图像为3) location = ‘test/’; %做颜色空间变换，变成亮度Y，色度I，饱和度Q for i=1:count I_origin = im2double(imread([location num2str(i) ‘.jpg’])); I_origin_yiq(i,:,:,:) = rgb2ntsc(I_origin); %subplot(2,count+1,i); figure; imshow(I_origin); title([‘原图像’ num2str(i)]); end %对所有图像的各个颜色分量进行小波分解 for i=1:count for d=1:colors I_cur(:,:)=I_origin_yiq(i,:,:,d); %进行小波分解 [C_cur,S]=wavedec2(I_cur,level,‘db4’); %记录第i个图像的第d个颜色分量的分解系数 C(i,d,:)=C_cur; end end %分别在分解的各层处理系数 %初始化为第一个图像的系数 Ca_len = S(1,1)*S(1,2); C_len = length(C_cur); C_fusion = zeros(3,C_len); %选择系数在3x3区域内方差比较大的 for l=1:level %该层的大小 row = S(level+2-l,1); col = S(level+2-l,2); l_select = ones(row,col); %记录每一点选择的哪幅图片的系数,映射表 l_ch_f = zeros(3,row,col); %初始化 l_cv_f = zeros(3,row,col); l_cd_f = zeros(3,row,col); l_ch = zeros(count,3,row,col); l_cv = zeros(count,3,row,col); l_cd = zeros(count,3,row,col); %提取第一张图像的高频系数 for d=1:colors l_ch(1,d,:,:) = detcoef2(‘h’,C(1,d,:),S,l); l_cv(1,d,:,:) = detcoef2(‘v’,C(1,d,:),S,l); l_cd(1,d,:,:) = detcoef2(‘d’,C(1,d,:),S,l); end %最后一层才提取低频系数 if l==level l_ca_f = zeros(3,row,col); l_ca = zeros(count,3,row,col); for d=1:colors l_ca(1,d,:,:) = appcoef2(C(1,d,:),S,‘db4’,l); end end %到当前图像为止细节分量的最大的方差和，只计算亮度的差别，初始化为第一幅图像 % l_f_var = kvar(l_ch_f) + kvar(l_cv_f) + kvar(l_cd_f); l_f_var = kgrads(l_ch(1,1,:,:),‘h’) + kgrads(l_cv(1,1,:,:),‘v’); %提取所有图像在第l层的高频系数，并计算方差，构造映射表 for i=2:count for d=1:colors l_ch(i,d,:,:) = detcoef2(‘h’,C(i,d,:),S,l); l_cv(i,d,:,:) = detcoef2(‘v’,C(i,d,:),S,l); l_cd(i,d,:,:) = detcoef2(‘d’,C(i,d,:),S,l); end l_ch_var = kgrads(l_ch(i,1,:,:),‘h’); l_cv_var = kgrads(l_cv(i,1,:,:),‘v’); if l==level for d=1:colors l_ca(i,d,:,:) = appcoef2(C(i,d,:),S,‘db4’,l); end end %在同一分辨率进行分组，使得各个朝向的选择是一样的。这里只选取了水平和垂直的 两个方向 for r=1:row for c=1:col if l_ch_var(r,c) + l_cv_var(r,c) > l_f_var(r,c) l_select(r,c)=i; l_f_var(r,c)=l_ch_var(r,c) + l_cv_var(r,c); end end end end %一致性检验：如果某点周围的绝大部分是从i图像中获得的，则该点也从i图像获得 l_select_con = kconsist(l_select); %选取系数，由于是在小波变换域上选取系数，虽然只是用亮度来衡量，但也必须对色度和 饱和度进行小波分解，选取系数 for r=1:row for c=1:col for d=1:colors l_ch_f(d,r,c) = l_ch(l_select_con(r,c),d,r,c); l_cv_f(d,r,c) = l_cv(l_select_con(r,c),d,r,c); l_cd_f(d,r,c) = l_cd(l_select_con(r,c),d,r,c); end end end if l==level for r=1:row for c=1:col for d=1:colors l_ca_f(d,r,c) = l_ca(l_select_con(r,c),d,r,c); end end end end %定位第l层系数在C中的起始位置 S_prod = prod(S,2)*3; beg = Ca_len + sum(S_prod(2:(level+1-l))); beg1 = beg+row*col; beg2 = beg1+row*col; %保存第l层的融合系数 for d=1:colors C_fusion(d,(beg+1):(beg+row*col)) = l_ch_f(d,:); C_fusion(d,(beg1+1):(beg1+row*col)) = l_cv_f(d,:); C_fusion(d,(beg2+1):(beg2+row*col)) = l_cd_f(d,:); end end for d=1:colors C_fusion(d,1:Ca_len) = l_ca_f(d,:); end %重构原信号 for d=1:colors I_result(:,:,d) = waverec2(C_fusion(d,:),S,‘db4’); end if colors==1 I_result(:,:,2) = I_origin_yiq(1,:,:,2); I_result(:,:,3) = I_origin_yiq(1,:,:,3); end I_result = ntsc2rgb(I_result); imwrite(I_result,[location ‘f’ num2str(count) num2str(level) num2str (colors) ‘.jpg’]); imwrite(l_select_con/count,[location ‘mgrads’ num2str(count) num2str (level) ‘con.jpg’]); figure; imshow(I_result,[]); title([‘融合图像’]);
1701010268
1701010269
其中求解图像梯度的程序如下：
1701010270
1701010271
function xresult = kgrads(x,d) %计算梯度， d=‘h’，检测行与行之间，即列方向上的差异(梯度)，即水平的边缘 if ndims(x)==3 %判断x的维数 xresult = zeros(size(x,2),size(x,3)); xresult(:) = x(:); %赋值 x = xresult; elseif ndims(x)==4 xresult = zeros(size(x,3),size(x,4)); xresult(:) = x(:); x = xresult; else xresult = x; end [row,col] = size(x); blocks = 3; %3x3大小窗口差分的平方和 %计算差分 if d==‘v’ for r=1:row for c=1:(col-1) xresult(r,c) = (xresult(r,c) - xresult(r,c+1))^2; %列与列之间，即行方向上的差异 end xresult(r,col) = xresult(r,col-1); end else for c=1:col for r=1:(row-1) xresult(r,c) = (xresult(r,c) - xresult(r+1,c))^2; %行与行之间，即列方向上的差异 end xresult(row,c) = xresult(row-1,c); end end patch = (blocks-1)/2; xx = zeros(row+patch*2,col+patch*2); %初始化 xx((patch+1):(row+patch),(patch+1):(col+patch)) = xresult(:,:); %赋值 for r=1:row for c=1:col s_grads = 0; for rr=r:(r+patch*2) for cc=c:(c+patch*2) s_grads = s_grads + xx(rr,cc); %计算梯度 end end xresult(r,c) = s_grads; end end
1701010272
1701010273
运行程序输出图形如图9-4所示。
1701010274
1701010275
