Code Part8 and MATLAB Digital Image Processing in digital image processing principles and practices (matlab)

Code Part8 and MATLAB Digital Image Processing in digital image processing principles and practices (matlab)
This article is Part 8 of the Code series in the book "Digital Image Processing Principles and Practices (MATLAB)", which records the code from pages 375th to 415th for readers to download and use. So far, the code release of this book is coming to an end. I hope these source codes can help readers who need them. For the code execution result, see the original book layout. We recommend that you read the following before downloading the Code:

Introduction to the code release in the book "Digital Image Processing Principles and Practices (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 (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 (abs (Dh) + abs (Dv )));
End

Lifecycle 7

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,
% Center Difference Method for gradient and differentiation Calculation
% WN N EN
% W O E
% WS S ES
I _x = (I (:, [2: nx])-I (:, [1 1: nx-1])/2; % Ix = (E-W)/2
I _y = (I ([2: ny], :)-I ([1 1: ny-1], :))/2; % Iy = (S-N)/2
I _xx = I (:, [2: nx]) + I (:, [1 1: nx-1])-2 * I; % Ixx = E + W-2 * O
I _yy = I ([2: ny], :) + I ([1 1: ny-1], :)-2 * I; % Iyy = S + N-2 * O
Dp = I ([2: ny], [2: nx]) + I ([1: ny-1], [1: nx-1]);
Dm = I ([1: ny-1], [2: nx]) + I ([2: ny], [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 ');

% Simulate motion blur, generate a point spread function PSF, the corresponding linear motion length is 21 pixels, angle is 11
LEN = 21;
THETA = 11;
PSF = fspecial ('Demo', LEN, THETA );
Blurred = imfilter (I, PSF, 'conv', 'Circular ');
Figure, imshow (blurred), title ('blurred image ');

% Restore Motion Blur Images
Result_w1 = deconvwnr (blurred, PSF );
Figure, imshow (result_w1), title ('restoration of Blurred image ')

% Simulate addition Noise
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 ')

% Assume that the image is restored without noise
Estimated_nsr = 0;
Wnr2 = deconvwnr (blurred_noisy, PSF, estimated_nsr );
Figure, imshow (wnr2)
Title ('restoration of Blurred, Noisy Image Using NSR = 0 ')

% Use a better signal-noise power ratio rating to restore images
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

% Load Original Image
I = imread ('board. tif ');
I = I (50 + [], 2 + [], :);
Figure, imshow (I)
Title ('al inal image ')

% Fuzzy processing
PSF = fspecial ('gaussian ', 5, 5 );
Blurred = imfilter (I, PSF, 'shortric ', 'conv ');
Figure, imshow (Blurred );
Title ('blurred image ')

% Add Noise
V = 0.002;
BlurredNoisy = imnoise (Blurred, 'gaussian, 0, V );
Figure, imshow (BlurredNoisy)
Title ('blurred and Noisy image ')

P405

% Image restoration using Lucy-Richards algorithm (5 iterations)
Luc1 = deconvlucy (BlurredNoisy, PSF, 5 );
Figure, imshow (luc1)

% Use the Lucy-Richards algorithm to restore images
Luc1_cell = deconvlucy ({BlurredNoisy}, PSF, 5 );
Luc2_cell = deconvlucy (luc1_cell, PSF );
Luc2 = im2uint8 (luc2_cell {2 });
Figure, imshow (luc2 );

% Control threshold restoration Effect
DAMPAR = im2uint8 (3 * sqrt (V ));
Luc3 = deconvlucy (BlurredNoisy, PSF, 15, DAMPAR );
Figure, imshow (luc3 );

P412

% Calculate the dark channel image of an image. The window size is 15*15.
ImageRGB = imread('picture.bmp ');
ImageRGB = double (imageRGB );
ImageRGB = imageRGB./255;
Dark = darkChannel (imageRGB );

% Save intermediate results
Required imwrite(dark,'stadium_darkchannel.png ');

% Select the brightest 0.1% pixels in the dark channel to obtain the atmospheric light
[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;

% Solve the transmittance, and retain a certain degree of haze through the omega parameter to avoid damage to the sense of reality

ω = 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;

% Is used to save intermediate results
% Imwrite (t, 'house_t.png ');

% Get more precise transmittance chart through guided Filtering
R = 60;
Eps = 10 ^-6;

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

% Is used to save intermediate results
% Imwrite (refined_t, 'house_r_t.png ');

RefinedRadiance = getRadiance (atmospheric, imageRGB, refined_t );

% Is used to save the result after the defog operation
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 (block_struct.data) * d;
Dark = blockproc (a, [15 15], fun );

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

Supplementary function: In the original book, only functions that are not listed in full (mainly functions that need to be called when the above code is executed) are as follows. The code for guided filtering comes from the personal homepage of the original author of the algorithm.

Supplement function 1 -- boxfilter

Function imDst = boxfilter (imSrc, r)

% Boxfilter o (1) time box filtering using cumulative sum
%
%-Definition imDst (x, y) = 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

Supplement function 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 (shocould be a color (RGB) image)
%-Filtering input image: p (shocould 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 blocks t 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 linear Ric 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

Supplement function 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 (J ))));


(Code release is not complete. Please wait for the future ...)