MATLAB 拉普拉斯残差金字塔

过去研究过高斯金字塔，不过今天看来，当时似乎搞错了，金字塔上下层应该是采样的关系，而不是缩放的关系，而且不同层的sigma应该也有所不同。

拉普拉斯金字塔在图像融合中有所应用，方法是首先对两个待融合图像求拉普拉斯残差金字塔，然后用模板对每一级残差图像进行融合得到融合后图像的残差金字塔，然后对这个金字塔进行重构就能得到最终的融合图像，图像各尺度细节得到保留。（注：融合时模板一般会先用高斯函数模糊一下）

不过这里不实现融合，只实现拉普拉斯金字塔的建立。

建立金字塔的方框图如下（摘自《数字图像处理（第二版）》）：

通过j-1级近似图像和j级残差图像就能恢复j级输入图像。

近似滤波器就是拉普拉斯滤波模板。

插入滤波器是对拉普拉斯滤波模板采样形成的新滤波模板。

使用的拉普拉斯滤波模板如下：

w=1/256*[1  4   6   4  1; 
                4 16 24 16 4;
                6 24 36 24 6;
                4 16 24 16 4;
                1  4   6   4 1];

产生的残差金字塔：

前面四张图像都是残差，最后一张图像是原图的4级近似，通过最后一张图和残差图，就能重构回原图。

matlab代码如下：

main.m

clear all; close all; clc;

img=double(imread('lena.jpg'));
[m n]=size(img);

w=1/256*[1  4  6  4 1;      %拉普拉斯滤波器
         4 16 24 16 4;
         6 24 36 24 6;
         4 16 24 16 4;
         1  4  6  4 1];

imgn{1}=img;
for i=2:5                   %滤波，下采样
   imgn{i}=imfilter(imgn{i-1},w,'replicate');
   imgn{i}=imgn{i}(1:2:size(imgn{i},1)-1,1:2:size(imgn{i},2)-1); %i-1级近似
end
       
for i=5:-1:2        %调整图像大小
   imgn{i-1}=imgn{i-1}(1:2*size(imgn{i},1),1:2*size(imgn{i},2)); 
end

 for i=1:4          %获得残差图像，i级预测残差
    imgn{i}=imgn{i}-expand(imgn{i+1},w);     
 end
 
for i=4:-1:1        %残差图像重构原图像
    imgn{i}=imgn{i}+expand(imgn{i+1},w);
end

imshow(uint8(imgn{1}));

expand.m

function re=expand(img,w)

    img=double(img);
    w=w*4;
    
    [m n]=size(img);
    [M N]=size(w);
    
    %插入滤波器                                   
    w_up_left=w(1:2:M,1:2:N); 
    w_up_right=w(1:2:M,2:2:N); 
    w_down_left=w(2:2:M,1:2:N); 
    w_down_right=w(2:2:M,2:2:N); 

    img_up_left=imfilter(img,w_up_left,'replicate','same');
    img_up_right=imfilter(img,w_up_right,'replicate','same');   
    img_down_left=imfilter(img,w_down_left,'replicate','same');  
    img_down_right=imfilter(img,w_down_right,'replicate','same');  
  
    re= zeros(m*2,n*2);             %上采样
    re(1:2:m*2,1:2:n*2)=img_up_left;
    re(2:2:m*2,1:2:n*2)=img_up_right;
    re(1:2:m*2,2:2:n*2)=img_down_left;
    re(2:2:m*2,2:2:n*2)=img_down_right;

end