二元回归解决图像恢复问题（图像去噪）

1 概述

       对于本次实验，进行了多次方法的更新迭代，最终使用多变量固定窗口二元高斯回归方法，自我测试效果达到：0.8噪音比下，原受损图像距离774.9727，处理后距离33.2362；0.6噪音比下，原受损图像距离671.1107，处理后距离20.1649；与一元线性回归相比提升三分之一以上。效果良好。

Figure 0

      

从过程来说，依次进行了如下尝试与测试：

              按行一元高斯线性回归

              整幅图片二元高斯线性回归

              增加Phi数量进行整幅图片二元高斯线性回归

              将图片均分为10 x10 按区域进行二元高斯线性回归

              将图片均分为100x 100 按区域进行二元高斯线性回归

              将图片均分为100x 100 按区域进行二元高斯线性回归 利用双边滤波器平滑处理

              将图片均分为100x 100 按区域进行二元高斯线性回归 增加sigma维度

              将图片分为固定大小5x5区域按区域进行二元高斯线性回归 增加sigma维度

              将图片分为固定大小5x5区域按区进行二元高斯线性回归 加sigma维度 并行化

       效果如Figure 0 所示。

接下来将依次介绍详细过程及测试效果。

2 算法依次介绍

2.1 按行一元高斯线性回归

       这是助教提供的算法。思路是对每一行，根据非噪音进行回归分析，算出权重，然后计算出信息丢失点。其中，根据该点通道值是否为零判断是否是信息丢失点。Phi函数为高斯函数，sigma为0.01，Phi函数数量为50.

       Figure 1为测试原图。640x640

       在0.6的噪音比下，受损图片与原图的距离为671.1107。在利用按行线性回归后，效果如Figure 3，距离为31.9975。

Figure 1 测试原图

 

Figure 2 0.6噪音比图

 

Figure 3

 

       在0.8的噪音比下，图片如Figure 4，距离达到774.9727。利用一元线性回归处理后，效果如Figure 5，距离达到54.2619.

Figure 4

 

Figure 5

2.2 整幅图片二元高斯线性回归

       直接对整幅图片，以行、列为二维变量，其他不变进行二元高斯线性回归。高斯函数参数与按行一元相似。在0.6噪音比下，效果如Figure 6，距离变化为671.1107 - 201.1378。在0.8噪音比下，效果如Figure 7，距离变化为774.9727 - 231.9823。

       效果非常差。我认为是在整幅图像下，像素点之间的相关性较弱，难以根据50个参数建立起回归函数。即使将参数扩展到1000时，效果也并未发生较大变化。0.6噪音比下，距离变化为671.1107 - 201.1402；0.8噪音比下，距离变化为774.9727 - 231.9911.

Figure 6

Figure 7

      

代码：

      

basisNum = 50; % define the number of basis functions.

sigma = 0.01; % define the standard deviation.

Phi_mu_x = linspace(1, cols, basisNum)/cols; % set the mean x value of each basis function

Phi_mu_y = linspace(1, rows, basisNum)/rows; % set the mean y value of each basis function

Phi_mu = [Phi_mu_x(:) Phi_mu_y(:)];

Phi_sigma = sigma * ones(2, basisNum); % here we set the standard deviation to the same value for brevity.

 

% use pixel index as the independent variable in the regression function

x = 1:cols;

x = (x - min(x)) / (max(x)-min(x));

y = 1:rows;

y = (y - min(y)) / (max(y)-min(y));

[X, Y] = meshgrid(x, y);

r = [X(:) Y(:)];

 

% select the missing pixels randomly

resImg = corrImg;

 

% for each channel

for k = 1:channels

    % select the missing pixels

    msk = noiseMask(:, :, k);

    msk = msk(:,:);

    misId = find(msk<1);

    misNum = length(misId);

    ddId = find(msk>=1);

    ddNum = length(ddId);

   

    % compute the coefficients

    Phi = [ones(ddNum, 1) zeros(ddNum, basisNum-1)];

    for j = 2: basisNum

        Phi(:, j) = mvnpdf( r(ddId',:), Phi_mu(j-1), Phi_sigma(:,j-1)') * sqrt(2*pi) * Phi_sigma(j-1);

    end

    corrImg_channal = corrImg(:,:,k);

    w = pinv(Phi' * Phi) * Phi' * corrImg_channal(ddId);

 

    % restore the missing values

    Phi1 = [ones(misNum, 1) zeros(misNum, basisNum-1)];

    for j = 2: basisNum

        Phi1(:, j) = mvnpdf( r(misId,:), Phi_mu(j-1), Phi_sigma(:,j-1)') * sqrt(2*pi) * Phi_sigma(j-1);

    end

    resImg_channal = resImg(:,:,k);

    resImg_channal(misId) = w' * Phi1';

    resImg(:,:,k) = resImg_channal;

end

resImg = min(resImg, 1);

resImg = max(resImg, 0);

 

2.3 将图片均分为10 x 10 按区域进行二元高斯线性回归

       很自然的想到，将整幅图片划分为几个部分，在每个部分内进行二元高斯线性回归。实验发现，果然要比整幅图进行回归效果好很多。在0.6噪音比下，距离变化为671.1107 - 67.6776；在0.8噪音比下，距离变化为774.9727 - 78.4643。效果如Figure 8, 9所示

Figure 8

 

Figure 9

 

      Code:

      

%% ==================learn the coefficents in row and coloumn linear regression function area by area=================

% In this section, we use gaussian kernels as the basis functions. And we

% do regression analysis row and colomun at a time.

% but we first divde the image into some square areas

 

smooth = false; % if smooth

areaNum = [100, 100];

basisNum = 50; % define the number of basis functions.

sigma = 10; % define the standard deviation.

Phi_mu_x = linspace(1, cols, basisNum)/cols; % set the mean x value of each basis function

Phi_mu_y = linspace(1, rows, basisNum)/rows; % set the mean y value of each basis function

Phi_mu = [Phi_mu_x(:) Phi_mu_y(:)];

Phi_sigma = sigma * ones(2, basisNum); % here we set the standard deviation to the same value for brevity.

 

% use pixel index as the independent variable in the regression function

x = 1:cols;

x = (x - min(x)) / (max(x)-min(x));

y = 1:rows;

y = (y - min(y)) / (max(y)-min(y));

[X, Y] = meshgrid(x, y);

r = [X(:) Y(:)];

 

% select the missing pixels randomly

resImg = corrImg;

 

% divide into areas

delta = floor([rows, cols] ./ areaNum);

areaBound_x = [1:delta(1):rows rows];

areaBound_y = [1:delta(2):cols cols];

 

% for each channel

for k = 1:channels

    for m = 1:length(areaBound_x)-1% row of area

        for n = 1:length(areaBound_y)-1% column of area

            % select the missing pixels in this area

            msk = noiseMask(areaBound_x(m):areaBound_x(m+1), areaBound_y(n):areaBound_y(n+1), k);

            msk = msk(:,:);

            misId = find(msk<1);

            misNum = length(misId);

            ddId = find(msk>=1);

            ddNum = length(ddId);

 

            % compute the coefficients

            Phi = [ones(ddNum, 1) zeros(ddNum, basisNum-1)];

            for j = 2: basisNum

                Phi(:, j) = mvnpdf( r(ddId',:), Phi_mu(j-1), Phi_sigma(:,j-1)');

            end

            corrImg_sub = corrImg(areaBound_x(m):areaBound_x(m+1), areaBound_y(n):areaBound_y(n+1), k);

            w = pinv(Phi' * Phi) * Phi' * corrImg_sub(ddId);

 

            % restore the missing values

            Phi1 = [ones(misNum, 1) zeros(misNum, basisNum-1)];

            for j = 2: basisNum

                Phi1(:, j) = mvnpdf( r(misId,:), Phi_mu(j-1), Phi_sigma(:,j-1)');

            end

            resImg_sub = resImg(areaBound_x(m):areaBound_x(m+1), areaBound_y(n):areaBound_y(n+1), k);

            resImg_sub(misId) = w' * Phi1';

            resImg(areaBound_x(m):areaBound_x(m+1), areaBound_y(n):areaBound_y(n+1), k) = resImg_sub;

        end

    end

end

resImg = min(resImg, 1);

resImg = max(resImg, 0);

 

 

2.4 将图片均分为100 x 100 按区域进行二元高斯线性回归

       可以看到由于分块过大，图片割裂感明显，将之改为100x100，效果显著提升。

       0.6噪音比下，671.1107 - 20.1649；0.8噪音比下，774.9727 - 74.5682。效果如Figure 10，11.

Figure 10

 

Figure 11

 

 

2.5 将图片均分为100 x 100 按区域进行二元高斯线性回归并利用双边滤波器进行平滑处理

       由于分块计算，图片看起来不够圆滑，于是选择滤波器进行平滑处理。比较理想的带有去噪功能的滤波器为高斯滤波器与双边滤波器。鉴于双边滤波器不仅能去噪和平滑，而且能不影响原图边界，因而选择了高斯滤波器。但高斯滤波器在实验过程中，虽然在图片观感上更为流畅和自然，但在大部分情况下距离上升。在少部分情况下，距离下降。经过分析发现，对回归后仍有噪点的图像进行滤波，效果会提升，其他反而会下降。

       在0.6与0.8下的效果如Figure12 13

Figure 12

Figure 13

 

2.6 将图片均分为100 x 100 按区域进行二元高斯线性回归增加sigma维度

       在0.8噪音比下不断测试，思考认为，增大sigma会加强其对高噪音下的效果表现。但是较大的sigma会降低低噪音比环境下的效果。因而想增加sigma的维度，增加Phi函数，利用最大似然估计算出最适合的sigma。

       因而在原先的基础上，增加了sigma的维度。

       以[100 1 0.01]为sigma参数进行实验，实验结果为0.6噪音比下，671.1107 - 20.1649；0.8噪音比下，774.9727 - 29.441。效果如Figure14 15

Figure 14

Figure 15

       Code:

smooth = false; % if smooth

delta = [5 5];

sub_basisNum = 50; % define the number of basis functions of each sigma

sigma = [100 1 0.01]; % define the standard deviation.

basisNum = sub_basisNum * length(sigma);

Phi_mu_x = linspace(1, cols, sub_basisNum)/cols; % set the mean x value of each basis function

Phi_mu_y = linspace(1, rows, sub_basisNum)/rows; % set the mean y value of each basis function

Phi_mu = [];

Phi_sigma = [];

 

for i=1:length(sigma)

    Phi_mu_new = [Phi_mu;Phi_mu_x(:) Phi_mu_y(:)];

    Phi_sigma_new = [Phi_sigma sigma(i) * ones(2, basisNum)]; % here we set the standard deviation to the same value for brevity.

    Phi_mu = Phi_mu_new;

    Phi_sigma = Phi_sigma_new;

end

 

% use pixel index as the independent variable in the regression function

x = 1:cols;

x = (x - min(x)) / (max(x)-min(x));

y = 1:rows;

y = (y - min(y)) / (max(y)-min(y));

[X, Y] = meshgrid(x, y);

r = [X(:) Y(:)];

 

% select the missing pixels randomly

resImg = corrImg;

 

% divide into areas

areaBound_x = [1:delta(1):rows rows];

areaBound_y = [1:delta(2):cols cols];

 

% for each channel

for k = 1:channels

    for m = 1:length(areaBound_x)-1% row of area

        for n = 1:length(areaBound_y)-1% column of area

            % select the missing pixels in this area

            msk = noiseMask(areaBound_x(m):areaBound_x(m+1), areaBound_y(n):areaBound_y(n+1), k);

            msk = msk(:,:);

            misId = find(msk<1);

            misNum = length(misId);

            ddId = find(msk>=1);

            ddNum = length(ddId);

 

            % compute the coefficients

            Phi = [ones(ddNum, 1) zeros(ddNum, basisNum-1)];

            for j = 2: basisNum

                Phi(:, j) = mvnpdf( r(ddId',:), Phi_mu(j-1), Phi_sigma(:,j-1)');

            end

            corrImg_sub = corrImg(areaBound_x(m):areaBound_x(m+1), areaBound_y(n):areaBound_y(n+1), k);

            w = pinv(Phi' * Phi) * Phi' * corrImg_sub(ddId);

 

            % restore the missing values

            Phi1 = [ones(misNum, 1) zeros(misNum, basisNum-1)];

            for j = 2: basisNum

                Phi1(:, j) = mvnpdf( r(misId,:), Phi_mu(j-1), Phi_sigma(:,j-1)');

            end

            resImg_sub = resImg(areaBound_x(m):areaBound_x(m+1), areaBound_y(n):areaBound_y(n+1), k);

            resImg_sub(misId) = w' * Phi1';

            resImg(areaBound_x(m):areaBound_x(m+1), areaBound_y(n):areaBound_y(n+1), k) = resImg_sub;

        end

    end

end

resImg = min(resImg, 1);

resImg = max(resImg, 0);

 

%smooth

if smooth

    smooth_d = 6; 

    smooth_sigma = [3 0.1]; 

    resImg = BilateralFilt2(double(resImg), smooth_d, smooth_sigma); 

end

 

 

2.7 将图片分为固定大小5x5区域 按区域进行二元高斯线性回归增加sigma维度

但将图片均分并不是合理的做法，回归效果应和具体窗口大小有关，如果对于不同大小图片进行相同份数均分，其效果实际并不相同。所以讲原先均分的区域划分法改成固定窗口大小的区域划分法。效果变化不大，但对于各种分辨率的图片效果更为稳定。

       最后，为加快速度，进行了并行化。

smooth = false; % if smooth

delta = [5 5];

sub_basisNum = 50; % define the number of basis functions of each sigma

sigma = [100 1 0.01]; % define the standard deviation.

basisNum = sub_basisNum * length(sigma);

Phi_mu_x = linspace(1, cols, sub_basisNum)/cols; % set the mean x value of each basis function

Phi_mu_y = linspace(1, rows, sub_basisNum)/rows; % set the mean y value of each basis function

Phi_mu = [];

Phi_sigma = [];

 

for i=1:length(sigma)

    Phi_mu_new = [Phi_mu;Phi_mu_x(:) Phi_mu_y(:)];

    Phi_sigma_new = [Phi_sigma sigma(i) * ones(2, basisNum)]; % here we set the standard deviation to the same value for brevity.

    Phi_mu = Phi_mu_new;

    Phi_sigma = Phi_sigma_new;

end

 

% use pixel index as the independent variable in the regression function

x = 1:cols;

x = (x - min(x)) / (max(x)-min(x));

y = 1:rows;

y = (y - min(y)) / (max(y)-min(y));

[X, Y] = meshgrid(x, y);

r = [X(:) Y(:)];

 

% select the missing pixels randomly

resImg = corrImg;

 

% divide into areas

areaBound_x = [1:delta(1):rows rows];

areaBound_y = [1:delta(2):cols cols];

 

parpool('local', 2);

% for each channel

parfor k = 1:channels

    for m = 1:length(areaBound_x)-1% row of area

        for n = 1:length(areaBound_y)-1% column of area

            % select the missing pixels in this area

            msk = noiseMask(areaBound_x(m):areaBound_x(m+1), areaBound_y(n):areaBound_y(n+1), k);

            msk = msk(:,:);

            misId = find(msk<1);

            misNum = length(misId);

            ddId = find(msk>=1);

            ddNum = length(ddId);

 

            % compute the coefficients

            Phi = [ones(ddNum, 1) zeros(ddNum, basisNum-1)];

            for j = 2: basisNum

                Phi(:, j) = mvnpdf( r(ddId',:), Phi_mu(j-1), Phi_sigma(:,j-1)');

            end

            corrImg_sub = corrImg(areaBound_x(m):areaBound_x(m+1), areaBound_y(n):areaBound_y(n+1), k);

            w = pinv(Phi' * Phi) * Phi' * corrImg_sub(ddId);

 

            % restore the missing values

            Phi1 = [ones(misNum, 1) zeros(misNum, basisNum-1)];

            for j = 2: basisNum

                Phi1(:, j) = mvnpdf( r(misId,:), Phi_mu(j-1), Phi_sigma(:,j-1)');

            end

            x0 = areaBound_x(m);

            x1 = areaBound_x(m+1);

            y0 = areaBound_y(n);

            y1 = areaBound_y(n+1);

            resImg_sub = resImg(x0:x1, y0:y1, k);

            resImg_sub(misId) = w' * Phi1';

            resImg(x0:x1, y0:y1, k) = resImg_sub;

        end

    end

end

parpool close;

resImg = min(resImg, 1);

resImg = max(resImg, 0);

 

%smooth

if smooth

    smooth_d = 6; 

    smooth_sigma = [3 0.1]; 

    resImg = BilateralFilt2(double(resImg), smooth_d, smooth_sigma); 

end

 

 

 

 

2.8 小结

       最终采用了5x5固定大小窗口内做二元高斯回归的方法进行。看情况进行双边滤波器滤波，默认不进行滤波。并利用多组sigma作为Phi的参数，以此适应多噪音比环境下的回归效果。

       此后又对程序中给的山水画png进行了测试，0.6噪音比下, 349.2799 - 84.1094,效果如图Figure 16。 0.8噪音比下，403.3304- 108.6377，效果如Figure 17。而用按行一元高斯分布回归则是403.3304 - 186.1008，效果如Figure 18。

       分析可知，本算法对于细节丰富、复杂的图像，效果不如趋于统一图像好，在这样的情形下，效果与一元高斯分布回归效果相差不大。但在高噪音比下，效果要明显好于一元高斯分布回归。

 

Figure 16

Figure 17

Figure 18

 

3 实验结果

A:

                            

B:

                           

C: