[College of mathematics] improved TV Model in 057th Image Restoration

$ \ BF abstract $: This article gives a detailed description of section 6.2 of "methods of partial differential equations in image processing" Edited by Wang dakai.

 

$ \ BF keywords $: image restoration; TV model; matlab programming

 

1. Preface

In the process of image formation, transmission, and storage, the image quality may degrade ). the degraded image can be described by the mathematical model $ \ bee \ label {1: d} u_0 = h_d * F + n \ EEE $.

(1) $ f (x, y) $ is the ideal image;

(2) $ h_d (x, y) $ is the point-Spread Function (PSF) of the imaging system. The ideal PSF is $ \ delta (\ cdot) $, the actual PSF is divided

(A) blur caused by a defocus (point size changed by electron beam scan): $ \ Bex h_d (x, y) = e ^ {-\ frac {x ^ 2 + y ^ 2} {2 \ Sigma ^ 2 }}, \ EEx $

(B) Effect of atmospheric turbulence (changes in large-scale motion state due to small scale disturbance): $ \ Bex \ calf (h_d) (u, v) \ equiv H (u, v) = E ^ {-K (U ^ 2 + V ^ 2) ^ \ frac {5} {6}, \ EEx $

(C) motion blur: $ \ Bex h_d (u, v) = T \ frac {\ sin [\ Pi (Au + bv)]} {\ Pi (Au + bv)} e ^ {-J \ Pi (Au + bv)}, \ EEx $ where $ T $ is the shutter (shutter, specifies the camera's exposure time. $ (a, B) $ indicates the object's speed;

(3) $ n (x, y) $ is a supplementary noise. In this paper, we use the Gaussian white noise (white Gaussian noise, WGN) where the mean is zero and the variance is $ \ Sigma $ ).

 

Image Restoration is used to restore the original image based on the recorded $ u_0 $. traditional restoration methods mainly involve Inverse Filtering or deconvolution (for example, the deconvolution filter theory of winner) and constrained least square method.

 

Now, we will briefly introduce the constrained Least Square Method Related to the partial encryption restoration. its starting point is: output image under the constraints of $ \ Bex \ iint _ \ Omega \ sev {h_d * u-u_0} ^ 2 \ RD x \ RD y = \ const \ EEx $ images are less blurred than images to be processed (What are the requirements of image processing ), so blur should be almost the same as the image to be processed) should be as smooth as possible.

 

In this way, smooth measurement is the key to the least square method. the traditional method is to use $ \ DPS {\ iint _ \ Omega \ sev {\ n u} ^ 2 \ RD x \ RD y} $ as a smooth measure, however, it is incompatible with the inherent features of the image-there is a mutation (edge. therefore, Rudin, Osher, and Fatime first proposed to use $ \ DPS {\ iint _ \ Omega \ sev {du} $ as a smooth measurement, and created a total variation (TV) the restoration method, or ROF method. the TV method is a bit because the BV function allows skip while a binary image (or a more general slice constant image) it is a member ($ \ DPS {\ chi_e \ In BV (\ Omega )} $ but $ \ DPS {\ chi_e \ not \ In w ^ {1, 2} (\ Omega)} $ ).

 

2. Model Creation

A TV model with $ \ DPS {\ iint _ \ Omega \ sev {du }}$ as a smooth measurement can be expressed as the following functional functions $ \ bee \ label {2: e} e (u) = \ iint _ \ Omega \ sev {du} + \ frac {\ Lambda} {2} \ iint _ \ Omega \ sev {h_d * u-u_0} ^ 2 \ RD x \ RD y, \ EEE $ where $ \ Lambda $ is a multiplier. use formula $ \ Bex h (x) = H (-x) \ Ra \ int f \ cdot (G * H) = \ int (F * H) \ cdot g \ EEx $ we know that \ eqref {2: e} has an Euler-Laplace equation $ \ bee \ label {2: s}-\ Div \ sex {\ frac {\ n u} {\ sev {\ n u }}+ \ Lambda h_d * \ sex {h_d * u-u_0} = 0, \ EEE $ the corresponding gradient descent flow (GDF) is $ \ bee \ label {2: e} \ frac {\ P u} {\ p t} = \ Div \ sex {\ frac {\ n u} {\ sev {\ n u}-\ Lambda h_d * \ sex {h_d * u-u_0 }. \ EEE $ restore the image by using 'pde \ eqref {2: e, although it has its advantages (the least square method for smoothing measurements is compared with $ \ DPS {\ iint _ \ Omega \ sev {\ n u} ^ 2 \ RD x \ RD y} $ ), however

(1) On the one hand, the time step must be small, while the CPU overhead is high;

(2) On the other hand, a step (staircase) effect ($ \ DPS {\ chi_e \ In BV (\ Omega)} $) is generated when the image processing morphological principles are not fully met ). to alleviate these two defects, marquina and Osher proposed the following improvements: $ \ bee \ label {2: e_im} \ frac {\ P u} {\ p t} = \ sev {\ n u} \ Div \ sex {\ frac {\ n u} {\ sev {\ n U }}- \ Lambda \ sev {\ n u} h_d * \ sex {h_d * u-u_0 }. \ EEE $ \ eqref {2: e_im} has the following advantages:

(1) $ \ DPS {\ sev {\ n u} \ Div \ sex {\ frac {\ n u} {\ sev {\ n u }}=\ Kappa \ sev {\ n u }}$ non-linear diffusion, it is combined with the second item-Hamilton-Markov (HJ), and the CFL conditions required for numerical calculation are relaxed (the diffusion term uses the central difference (central difference ), HJ adopts upwind scheme ));

(2) \ eqref {2: e_im} remains unchanged under monotonic gray-scale transformation, while the morphological principle is almost satisfied (\ eqref {2: e_im} each item contains $ \ sev {\ n u} $, which is an inherent feature of the image-edge), and the step effect can be mitigated.

 

3. Numerical Algorithm

\ Eqref {2: e_im} can be directly solved using the explicit solution (explicit scheme: $ \ Bex U ^ {n + 1 }_{ IJ} = u ^ N _ {IJ} + \ lap t \ sex {P ^ N _ {IJ} + q ^ N _ {IJ }}, \ EEx $ where

(1) $ P ^ N _ {IJ} $ is the result of the center difference of the diffusion term: $ \ Bex P ^ N _ {IJ} = \ frac {\ sex {d ^ {(0 )} _ Yu ^ N _ {IJ }}^ 2 \ cdot d ^ {(0) }_{ XX} U ^ N _ {IJ}-2D ^ {(0 )} _ Xu ^ N _ {IJ} \ cdot d ^ {(0)} _ Yu ^ N _ {IJ} \ cdot d ^ {(0 )} _ {XY} U ^ N _ {IJ} + \ sex {d ^ {(0 )} _ Xu ^ N _ {IJ }}^ 2 \ cdot d ^ {(0 )} _ {YY} U ^ N _ {IJ }}{\ sex {d ^ {(0 )} _ Xu ^ N _ {IJ }}^ 2 + \ sex {d ^ {(0) }_yu ^ N _ {IJ }}^ 2 + \ ve }; \ EEx $

(2) $ q ^ N _ {IJ} $ is the result of HJ entry Windward format: $ \ Bex Q ^ N _ {IJ} & =&-\ Lambda \ beta ^ N _ {IJ} \ sev {\ Delta U ^ N _ {IJ }}, \ EEx $ \ Bex \ beta ^ N _ {IJ} = \ sex {h_d * u ^ n-u_0 }}_{ IJ }, \ EEx $ \ Bex \ sev {\ Delta U ^ N _ {IJ} & = & \ MAX \ sex {\ beta ^ N _ {IJ }, 0 }\\& \ cdot \ SQRT {\ MAX \ sex {d ^ {(-)} _ Xu ^ N _ {IJ }, 0} ^ 2 + \ min \ sex {d ^ {(+)} _ Xu ^ N _ {IJ} ^ 2 \ MAX \ sex {d ^ {(-)} _ Yu ^ N _ {IJ}, 0} ^ 2 + \ min \ sex {d ^ {(+ )} _ Yu ^ N _ {IJ }}^ 2 }\\& & + \ min \ sex {\ beta ^ N _ {IJ }, 0 }\\& \ cdot \ SQRT {\ min \ sex {d ^ {(-)} _ Xu ^ N _ {IJ }, 0} ^ 2 + \ MAX \ sex {d ^ {(+)} _ Xu ^ N _ {IJ} ^ 2 \ min \ sex {d ^ {(-)} _ Yu ^ N _ {IJ}, 0} ^ 2 + \ MAX \ sex {d ^ {(+)} _ Yu ^ N _ {IJ} ^ 2 }. \ EEx $

 

4. Numerical Test

Use the following Matlab code:

Clear all;
Close all;
CLC;

N = 1000; % set the number of iterations
D = 200; % output the current image 200 times per Iteration
Nn = 3; % initialize the number of output images

I01_imread('lenna.bmp ');
I0 = double (I0 );

Figure (1); imshow (uint8 (I0 ));

[NY, nx] = size (I0 );

I0 = Gauss (I0, 7, 3); % Gaussian Noise
I0 = I0 + 10 * randn ([NY, nx]); % supplementary Noise
Figure (2); imshow (uint8 (I0 ));

Lambda = 0.2; % Laplace Multiplier
Dt = 0.01; % time step
I = I0; % image initialization, storing Windward solution results
J = I0; % image initialization, store Roe Windward solution results

% Iteration started
For n = 1: N
IX = 0.5 * (I (:, [2: NX, nx])-I (:, [: nx-1]);
Iy = 0.5 * (I ([2: NY, NY], :)-I ([: ny-1], :));
Gradient = IX. ^ 2 + Iy. ^ 2 + EPS;

Ix_back = I-I (:, [: nx-1]);
Ix_forward = I (:, [2: NX, nx])-I;
Iy_back = I-I ([: ny-1], :);
Iy_forward = I ([2: NY, NY], :)-I;

Ixx = (I (:, [2: NX, nx])-I)-(I-I (:, [: nx-1]);
Ixy = 0.25 * (I ([2: NY, NY], [2: NX, nx])-I ([2: NY, NY], [1, 1: nx-1])...
-(I ([: ny-1], [2: NX, nx])-I ([: ny-1], [: nx-1]);
Iyy = (I ([2: NY, NY], :)-I)-(I-I ([: ny-1], :));

% Central difference
Diffusion = (Iy. ^ 2. * Ixx-2 * IX. * Iy. * ixy + ix. ^ 2. * iyy)./gradient;

% Windward Solution
Beta = Gauss (I, 7,3)-I0, 7,3 );

Upwind = max (Beta, 0 ).*...
SQRT (...
Max (ix_back, 0). ^ 2 + min (ix_forward, 0). ^ 2...
+ Max (iy_back, 0). ^ 2 + min (iy_forward, 0). ^ 2...
)...
+ Min (Beta, 0 ).*...
SQRT (...
Min (ix_back, 0). ^ 2 + max (ix_forward, 0). ^ 2...
+ Min (iy_back, 0). ^ 2 + max (iy_forward, 0). ^ 2...
);

% Roe upwind format
Upwind_x = ix_back;
Upwind_x (Beta. * IX <0) = ix_forward (Beta. * IX <0 );

Upwind_y = iy_back;
Upwind_y (Beta. * Iy <0) = iy_forward (Beta. * Iy <0 );

Roe_upwind = beta. * SQRT (upwind_x. ^ 2 + upwind_y. ^ 2 );

% Update Image
I = I + dt * (diffusion-Lambda * upwind );
J = J + dt * (diffusion-Lambda * roe_upwind );

% Output image
If Mod (n, d) = 0
Figure (NN );
Subplot (1, 2); imshow (uint8 (I ));
Subplot (1, 2); imshow (uint8 (j ));
Nn = nn + 1;
End
End

You can add noise to the lenna image and restore it. For details, see:

 

Here we call the following Matlab function to add Gaussian noise to the original image:

Function ig = Gauss (I, window, SIGMA)

% Function guass () implements smooth Guassian Filter
% Parameter description
% I --- functions to be smoothed
% Window --- Gaussian Kernel size (ODD)
% Simga --- Gaussian variance
% Ig --- return the Guassian smooth image

[NY, nx] = size (I );

Half_window = (window-1)/2; % convenient to center

If ny X = (-half_window: half_window );
Filter = exp (-X. ^ 2/(2 * sigma); % One-dimensional Guassian Function
Filter = filter/sum (filter); % Normalization
% Image Extension
XL = mean (I (:, [1: half_window]);
Xr = mean (I (:, [nx-half_window + 1, nx]);
I = [XL * ones (NY, half_window) I XR * ones (NY, half_window)];
% Expanded to column NX + window-1
IG = Conv (I, filter );
% Form (nx + window-1) + (window)-1 = NX + 2 (window-1) Column
IG = ig (:, window + half_window + 1: NX + window + half_window );
% The first convolution must use the data of the original image.
% The minimum value of image data used for item l in Convolution = L-Windows
%> = (Windows-1)/2 + 1 = position of the original image in the extended image
Else
% Two-dimensional convolution
X = ones (window, 1) * (-half_window: half_window); % abscissa
Y = x'; % Y coordinate

Filter = exp (-(X. ^ 2 + Y. ^ 2/(2 * sigma); % Two-Dimensional Guassian Function
Filter = filter/(sum (filter); % Normalization

% Image Extension
If (half_window> 1)
Xleft = mean (I (:, [1: half_window]) ';
% MATLAB calculates the average value for the column and returns the row vector.
Xright = mean (I (: [nx-half_window + 1: NX]) ';
Else
Xleft = I (:, 1 );
Xright = I (:, nx );
End

I = [xleft * ones (1, half_window), I, xright * ones (1, half_window)];

If (half_window> 1)
Xup = mean (I ([1: half_window], :));
Xdown = mean (I ([ny-half_window + 1, NY], :));
Else
Xup = I (1 ,:);
Xdown = I (NY ,:);
End

I = [ones (half_window, 1) * xup; I; ones (half_window, 1) * xdown];

IG = conv2 (I, filter, 'valid ');

End