首页 > 代码库 > 《数字图像处理原理与实践（MATLAB版）》一书之代码Part8

《数字图像处理原理与实践（MATLAB版）》一书之代码Part8

2024-08-09 14:52:15   211人阅读

本文系《数字图像处理原理与实践（MATLAB版）》一书之代码系列的Part8，辑录该书第375至第415页之代码，供有需要读者下载研究使用。至此全书代码发布已经接近尾声，希望这些源码能够对有需要的读者有所帮助。代码执行结果请参见原书配图，建议下载代码前阅读下文：

关于《数字图像处理原理与实践（MATLAB版）》一书代码发布的说明

http://blog.csdn.net/baimafujinji/article/details/40987807

P385-1

function y = tv(X)
    [M,N] = size(X);
    Dh = diff(X,[],1);
    Dh = [Dh;zeros(1,N)];
    Dv = diff(X,[],2);
    Dv = [Dv zeros(M,1)];
    y = sum(sum(sqrt(Dh.^2+Dv.^2)));
end

P385-2

function y = atv(X)
    [M,N] = size(X);
    Dh = -diff(X,[],1);
    Dh = [Dh;zeros(1,N)];
    Dv = -diff(X,[],2);
    Dv = [Dv zeros(M,1)];
    y = sum(sum(abs(Dh)+abs(Dv)));
end

P387

I = double(rgb2gray(imread(‘Lena.bmp‘)));
I0 = I;
ep=1; dt=0.25; lam=0;
ep2=ep^2; [ny,nx]=size(I);
iter = 80;

for i=1:iter,
    % 中心差法计算梯度和微分
    % WN  N  EN
    % W   O  E
    % WS  S  ES
    I_x = (I(:,[2:nx nx])-I(:,[1 1:nx-1]))/2; % Ix = (E-W)/2
    I_y = (I([2:ny ny],:)-I([1 1:ny-1],:))/2; % Iy = (S-N)/2
    I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I; % Ixx = E+W-2*O
    I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I; % Iyy = S+N-2*O
    Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1]);
    Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);
    I_xy = (Dp-Dm)/4; % Ixy = Iyx = ((ES+WN)-(EN+WS))/4
    
    Num = I_xx.*(ep2+I_y.^2)-2*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^2);
    Den = (ep2+I_x.^2+I_y.^2).^(3/2);
    I_t = Num./Den + lam.*(I0-I);
    I=I+dt*I_t;  %% evolve image by dt
end
imshow(I,[]);

P400

RGB = imread(‘moon.jpg‘);
I = rgb2gray(RGB);
J = imnoise(I,‘gaussian‘,0,0.025);
imshow(J)
K = wiener2(J,[5 5]);
figure, imshow(K)

P401

I = im2double(imread(‘cameraman.tif‘));
imshow(I), title(‘Original Image‘);

%模拟运动模糊，生成一个点扩散函数PSF，相应的线性运动长度为21个像素，角度为11
LEN = 21;
THETA = 11;
PSF = fspecial(‘motion‘, LEN, THETA);
blurred = imfilter(I, PSF, ‘conv‘, ‘circular‘);
figure, imshow(blurred), title(‘Blurred Image‘);

%对运动模糊图像进行复原
result_w1= deconvwnr(blurred, PSF);
figure, imshow(result_w1), title(‘Restoration of Blurred Image‘)

%模拟加性噪声
noise_mean = 0;
noise_var = 0.0001;
blurred_noisy = imnoise(blurred, ‘gaussian‘, noise_mean, noise_var);
figure, imshow(blurred_noisy), title(‘Simulate Blur and Noise‘)

%假设没有噪声的情况下复原图像
estimated_nsr = 0;
wnr2 = deconvwnr(blurred_noisy, PSF, estimated_nsr);
figure, imshow(wnr2)
title(‘Restoration of Blurred, Noisy Image Using NSR = 0‘)

%使用一个更好的信噪功率比评估值来复原图像
estimated_nsr = noise_var / var(I(:));
wnr3 = deconvwnr(blurred_noisy, PSF, estimated_nsr);
figure, imshow(wnr3)
title(‘Restoration of Blurred, Noisy Image Using NSR‘);

P404

%载入原始图像
I = imread(‘board.tif‘);
I = I(50+[1:256],2+[1:256],:);
figure, imshow(I)
title(‘Original Image‘)

%模糊处理
PSF = fspecial(‘gaussian‘,5,5);
Blurred = imfilter(I,PSF,‘symmetric‘,‘conv‘);
figure, imshow(Blurred);
title(‘Blurred Image‘)

%添加噪声
V = 0.002;
BlurredNoisy = imnoise(Blurred,‘gaussian‘,0,V);
figure, imshow(BlurredNoisy)
title(‘Blurred and Noisy Image‘)

P405

%利用露茜-理查德森算法复原图像（5次迭代）
luc1 = deconvlucy(BlurredNoisy, PSF, 5);
figure, imshow(luc1)

%利用露茜-理查德森算法复原图像
luc1_cell = deconvlucy({BlurredNoisy}, PSF, 5);
luc2_cell = deconvlucy(luc1_cell, PSF);
luc2 = im2uint8(luc2_cell{2});
figure, imshow(luc2);

%控制阈值的复原效果
DAMPAR = im2uint8(3*sqrt(V));
luc3 = deconvlucy(BlurredNoisy, PSF, 15, DAMPAR);
figure, imshow(luc3);

P412

%求一幅图像的暗通道图像，窗口大小为15*15
imageRGB = imread(‘picture.bmp‘);
imageRGB = double(imageRGB);
imageRGB = imageRGB./255;
dark = darkChannel(imageRGB);

%保存中间结果
%imwrite(dark,‘stadium_darkchannel.png‘);

% 选取暗通道中最亮的0.1%像素，从而求得大气光
[m, n, ~] = size(imageRGB);
imsize = m * n;
numpx = floor(imsize/1000);
JDarkVec = reshape(dark,imsize,1);
ImVec = reshape(imageRGB,imsize,3);

[JDarkVec, indices] = sort(JDarkVec);
indices = indices(imsize-numpx+1:end);

atmSum = zeros(1,3);
for ind = 1:numpx
    atmSum = atmSum + ImVec(indices(ind),:);
end

atmospheric = atmSum / numpx;

%求解透射率，并通过omega参数来选择保留一定程度雾霾，以免损坏真实感

omega = 0.95;

im = zeros(size(imageRGB));

for ind = 1:3
    im(:,:,ind) = imageRGB(:,:,ind)./atmospheric(ind);
end

dark_2 = darkChannel(im);

t = 1-omega*dark_2;

%用于保存中间结果
%imwrite(t, ‘house_t.png‘);

%通过导向滤波来获得更为精细的透射率图
r = 60;
eps = 10^-6;

refined_t = guidedfilter_color(imageRGB, t, r, eps);

%用于保存中间结果
%imwrite(refined_t , ‘house_r_t.png‘);

refinedRadiance = getRadiance(atmospheric, imageRGB, refined_t);

%用于保存去雾操作后的结果
imwrite(refinedRadiance , ‘house_dehaze.png‘);

P413

function dark = darkChannel(imRGB)

r=imRGB(:,:,1);
g=imRGB(:,:,2);
b=imRGB(:,:,3);

[m n] = size(r);
a = zeros(m,n);
for i = 1: m
    for j = 1: n
        a(i,j) = min(r(i,j), g(i,j));
        a(i,j)= min(a(i,j), b(i,j));
    end
end

d = ones(15,15);
fun = @(block_struct)min(min(block_struct.data))*d;
dark = blockproc(a, [15 15], fun);

dark = dark(1:m, 1:n);

补充函数：原书中限于篇幅未能全部列出之函数（主要是执行上述代码时需要调用的函数）如下，其中导向滤波的代码来自算法原作者的个人主页。

补充函数1——boxfilter

function imDst = boxfilter(imSrc, r)

%   BOXFILTER   O(1) time box filtering using cumulative sum
%
%   - Definition imDst(x, y)=sum(sum(imSrc(x-r:x+r,y-r:y+r)));
%   - Running time independent of r;
%   - Equivalent to the function: colfilt(imSrc, [2*r+1, 2*r+1], ‘sliding‘, @sum);
%   - But much faster.

[hei, wid] = size(imSrc);
imDst = zeros(size(imSrc));

%cumulative sum over Y axis
imCum = cumsum(imSrc, 1);
%difference over Y axis
imDst(1:r+1, :) = imCum(1+r:2*r+1, :);
imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);
imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);

%cumulative sum over X axis
imCum = cumsum(imDst, 2);
%difference over Y axis
imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);
imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);
imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);
end

补充函数2——guidedfilter_color

function q = guidedfilter_color(I, p, r, eps)
%   GUIDEDFILTER_COLOR   O(1) time implementation of guided filter using a color image as the guidance.
%
%   - guidance image: I (should be a color (RGB) image)
%   - filtering input image: p (should be a gray-scale/single channel image)
%   - local window radius: r
%   - regularization parameter: eps

[hei, wid] = size(p);
N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.

mean_I_r = boxfilter(I(:, :, 1), r) ./ N;
mean_I_g = boxfilter(I(:, :, 2), r) ./ N;
mean_I_b = boxfilter(I(:, :, 3), r) ./ N;

mean_p = boxfilter(p, r) ./ N;

mean_Ip_r = boxfilter(I(:, :, 1).*p, r) ./ N;
mean_Ip_g = boxfilter(I(:, :, 2).*p, r) ./ N;
mean_Ip_b = boxfilter(I(:, :, 3).*p, r) ./ N;

% covariance of (I, p) in each local patch.
cov_Ip_r = mean_Ip_r - mean_I_r .* mean_p;
cov_Ip_g = mean_Ip_g - mean_I_g .* mean_p;
cov_Ip_b = mean_Ip_b - mean_I_b .* mean_p;

% variance of I in each local patch: the matrix Sigma in Eqn (14).
% Note the variance in each local patch is a 3x3 symmetric matrix:
%           rr, rg, rb
%   Sigma = rg, gg, gb
%           rb, gb, bb
var_I_rr = boxfilter(I(:, :, 1).*I(:, :, 1), r) ./ N - mean_I_r .*  mean_I_r;
var_I_rg = boxfilter(I(:, :, 1).*I(:, :, 2), r) ./ N - mean_I_r .*  mean_I_g;
var_I_rb = boxfilter(I(:, :, 1).*I(:, :, 3), r) ./ N - mean_I_r .*  mean_I_b;
var_I_gg = boxfilter(I(:, :, 2).*I(:, :, 2), r) ./ N - mean_I_g .*  mean_I_g;
var_I_gb = boxfilter(I(:, :, 2).*I(:, :, 3), r) ./ N - mean_I_g .*  mean_I_b;
var_I_bb = boxfilter(I(:, :, 3).*I(:, :, 3), r) ./ N - mean_I_b .*  mean_I_b;

a = zeros(hei, wid, 3);
for y=1:hei
    for x=1:wid        
        Sigma = [var_I_rr(y, x), var_I_rg(y, x), var_I_rb(y, x);
            var_I_rg(y, x), var_I_gg(y, x), var_I_gb(y, x);
            var_I_rb(y, x), var_I_gb(y, x), var_I_bb(y, x)];
        %Sigma = Sigma + eps * eye(3);
        
        cov_Ip = [cov_Ip_r(y, x), cov_Ip_g(y, x), cov_Ip_b(y, x)];        
        
        a(y, x, :) = cov_Ip * inv(Sigma + eps * eye(3)); % Eqn. (14) in the paper;
    end
end

b = mean_p - a(:, :, 1) .* mean_I_r - a(:, :, 2) .* mean_I_g - a(:, :, 3) .* mean_I_b; % Eqn. (15) in the paper;

q = (boxfilter(a(:, :, 1), r).* I(:, :, 1)...
+ boxfilter(a(:, :, 2), r).* I(:, :, 2)...
+ boxfilter(a(:, :, 3), r).* I(:, :, 3)...
+ boxfilter(b, r)) ./ N;  % Eqn. (16) in the paper;
end

补充函数3——getRadiance

function J = getRadiance(A,im,tMat)
t0 = 0.1;
J = zeros(size(im));
for ind = 1:3
   J(:,:,ind) = A(ind) + (im(:,:,ind) - A(ind))./max(tMat,t0);
end

J = J./(max(max(max(J))));


（代码发布未完，请待后续...）

《数字图像处理原理与实践（MATLAB版）》一书之代码Part8


联系
我们