1. 参考文献
2. TV图像修复模型
2.1 TV模型
% demo_TV.m
% Author: HSW
% Date: 2015/3/25
% HARBIN INSTITUTE OF TECHNOLOGY
%
% set matlab
close all;
clear all;
clc;options.null = 0;
% read image
Img = imread('Image\butterfly.bmp');
% Img = imread('Image\peppers.bmp');
%Img = rgb2gray(Img);
Img = im2double(Img);
if max(Img(:)) < 2Img = Img*255;
endFlagColor = (size(Img,3) == 3);% set mask
SetMask = 2;
if SetMask == 1% read maskMask = imread('Mask\peppers_mask512.bmp');Mask = Mask > 5;MFlagColor = (size(Mask,3) == 3);if MFlagColor ~= FlagColor && FlagColor == 1Mask = repmat(Mask,[1,1,3]);elseif MFlagColor ~= FlagColor && FlagColor == 0Mask = Mask(:,:,1);end
elseif SetMask == 2% Interactively set maskif not(exist('grab_mode'))grab_mode = 'line';endoptions.grab_mode = grab_mode;if not(exist('grab_radius'))grab_radius = 4;endswitch grab_modecase 'points'options.r = grab_radius;U = grab_inpainting_mask(Img,options);case 'line'options.r = grab_radius;[U,options.point_list] = grab_inpainting_mask(Img,options);end %switchIin = find(U(:,:,1) == Inf);Iout = find(U(:,:,1) ~= Inf);m1 = length(Iin); % 缺损点的总数% product the maskMask = zeros(size(Img));if FlagColor == 1tmpMask = zeros([size(Img,1),size(Img,2)]);for channel = 1:3tmpMask(Iin) = 1 ;Mask(:,:,channel) = tmpMask;endelseMask(Iin) = 1; % 缺损区域为1end % if FlagColor
end% if SetMask
FlagNoise = 0; % default noiseless
if FlagNoise == 1sigma = 25; nImg = Img + sigma*randn(size(Img));
elsenImg = Img;
end
nImg = (1-Mask).*nImg;
PSNRin = 10*log10(255^2/mean((Img(:)-nImg(:)).^2));
InImg = nImg;
% Initial Image
if FlagColor == 1Positions = find(Mask(:,:,1) == 1);for channel = 1:3tmpnImg = nImg(:,:,channel);tmpnImg(Positions) = floor(255*rand(1,length(Positions))) + 1;InImg(:,:,channel) = tmpnImg;end
elsePositions = find(Mask == 1);randValue = floor(255*rand(1,length(Positions))) + 1;InImg(Positions) = randValue;
end% Main Total Variable Model
IterNum = 300;
lambda = 0.001;
if FlagNoise == 1lambda = 1;
end
I = InImg;
for iter = 1:IterNumI = TV_inpaint_gray_color(I,Mask,lambda,FlagNoise);if mod(iter,100) == 0figure;imshow(I/255,[]);title(['Results of iterNum = ',num2str(iter)]);end
end %for iterI = max(0,min(I,255));
PSNRout = 10*log10(255^2/mean((I(:) - Img(:)).^2));
figure;
subplot(1,3,1);
imshow(Img/255,[]);
title('Original Image');
subplot(1,3,2);
imshow(nImg/255,[]);
title(['Masked Image PSNR = ',num2str(PSNRin), ' dB']);
subplot(1,3,3);
imshow(I/255,[]);
title(['Inpainting Image PSNR = ', num2str(PSNRout), ' dB']);
2.2 全变分模型
function I = TV_inpaint_gray_color(I,M,lambda,FlagNoise)
% input:
% I: 待修复图像
% lambda: 保真项系数,在修复实验中不需要处理
% M：缺损区域mask, 缺损区域取值为1
% FlagNoise: 标示是否含噪声
% output;
% I: 修复结果
% Author: HSW
% Date: 2015/3/25
% HARBIN INSTITUTE OF TECHNOLOGY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dims = ndims(I);
epsilon = 0.001;
switch dimscase 3 %彩色图像if FlagNoise == 0for channel = 1:dimsGradientE_m = sqrt((I([2:end,end],:,channel) - I(:,:,channel)).^2 + ((I([2:end,end],[2:end,end],channel) + I(:,[2:end,end],channel) - I(:,[1,1:end-1],channel) - I([2:end,end],[1,1:end-1],channel))./4).^2);GradientS_m = sqrt( (I(:,:,channel) - I(:,[1,1:end-1],channel) ).^2 + ((I([2:end,end],:,channel) - I([1,1:end-1],:,channel) + I([2:end,end],[1,1:end-1],channel) - I([1,1:end-1],[1,1:end-1],channel) )./4).^2 );GradientW_m = sqrt( ( I(:,:,channel) - I([1,1:end-1],:,channel) ).^2 + ((I(:,[2:end,end],channel) - I(:,[1,1:end-1],channel) + I([1,1:end-1],[2:end,end],channel) - I([1,1:end-1],[1,1:end-1],channel))./4).^2 );GradientN_m = sqrt( ( ( I([2:end,end],:,channel) - I([1,1:end-1],:,channel) + I([2:end,end],[2:end,end],channel) - I([1,1:end-1],[2:end,end],channel))./4).^2 + ( I(:,[2:end,end],channel) - I(:,:,channel) ).^2 );GradientE = 1./(GradientE_m + epsilon);GradientS = 1./(GradientS_m + epsilon);GradientW = 1./(GradientW_m + epsilon);GradientN = 1./(GradientN_m + epsilon);Temp = GradientE + GradientS + GradientW + GradientN;Temp = Temp + (Temp == 0).*epsilon;I(:,:,channel) = (GradientE.*I([2:end,end],:,channel) + GradientW.*I([1,1:end-1],:,channel) + GradientN.*I(:,[2:end,end],channel) + GradientS.*I(:,[1,1:end-1],channel)).*M(:,:,channel)./Temp + (1-M(:,:,channel)).*I(:,:,channel);endelsefor channel = 1:dimsGradientE_m = sqrt((I([2:end,end],:,channel) - I(:,:,channel)).^2 + ((I([2:end,end],[2:end,end],channel) + I(:,[2:end,end],channel) - I(:,[1,1:end-1],channel) - I([2:end,end],[1,1:end-1],channel))./4).^2);GradientS_m = sqrt( (I(:,:,channel) - I(:,[1,1:end-1],channel) ).^2 + ((I([2:end,end],:,channel) - I([1,1:end-1],:,channel) + I([2:end,end],[1,1:end-1],channel) - I([1,1:end-1],[1,1:end-1],channel) )./4).^2 );GradientW_m = sqrt( ( I(:,:,channel) - I([1,1:end-1],:,channel) ).^2 + ((I(:,[2:end,end],channel) - I(:,[1,1:end-1],channel) + I([1,1:end-1],[2:end,end],channel) - I([1,1:end-1],[1,1:end-1],channel))./4).^2 );GradientN_m = sqrt( ( ( I([2:end,end],:,channel) - I([1,1:end-1],:,channel) + I([2:end,end],[2:end,end],channel) - I([1,1:end-1],[2:end,end],channel))./4).^2 + ( I(:,[2:end,end],channel) - I(:,:,channel) ).^2 );GradientE = 1./(GradientE_m + epsilon);GradientS = 1./(GradientS_m + epsilon);GradientW = 1./(GradientW_m + epsilon);GradientN = 1./(GradientN_m + epsilon);Temp = GradientE + GradientS + GradientW + GradientN + lambda;Temp = Temp + (Temp == 0).*epsilon;I(:,:,channel) = (GradientE.*I([2:end,end],:,channel) + GradientW.*I([1,1:end-1],:,channel) + GradientN.*I(:,[2:end,end],channel) + GradientS.*I(:,[1,1:end-1],channel) + lambda*I(:,:,channel)).*M(:,:,channel)./Temp + (1-M(:,:,channel)).*I(:,:,channel);endend% if FlagNoiseotherwise %灰度图像if FlagNoise == 0GradientE_m = sqrt((I([2:end,end],:) - I).^2 + ((I([2:end,end],[2:end,end]) + I(:,[2:end,end]) - I(:,[1,1:end-1]) - I([2:end,end],[1,1:end-1]))./4).^2);GradientS_m = sqrt( (I - I(:,[1,1:end-1]) ).^2 + ((I([2:end,end],:) - I([1,1:end-1],:) + I([2:end,end],[1,1:end-1]) - I([1,1:end-1],[1,1:end-1]) )./4).^2 );GradientW_m = sqrt( ( I - I([1,1:end-1],:) ).^2 + ((I(:,[2:end,end]) - I(:,[1,1:end-1]) + I([1,1:end-1],[2:end,end]) - I([1,1:end-1],[1,1:end-1]))./4).^2 );GradientN_m = sqrt( ( ( I([2:end,end],:) - I([1,1:end-1],:) + I([2:end,end],[2:end,end]) - I([1,1:end-1],[2:end,end]))./4).^2 + ( I(:,[2:end,end]) - I ).^2 );GradientE = 1./(GradientE_m + epsilon);GradientS = 1./(GradientS_m + epsilon);GradientW = 1./(GradientW_m + epsilon);GradientN = 1./(GradientN_m + epsilon);Temp = GradientE + GradientS + GradientW + GradientN;Temp = Temp + (Temp == 0).*epsilon;I = (GradientE.*I([2:end,end],:) + GradientW.*I([1,1:end-1],:) + GradientN.*I(:,[2:end,end]) + GradientS.*I(:,[1,1:end-1])).*M./Temp + (1-M).*I;elseGradientE_m = sqrt((I([2:end,end],:) - I).^2 + ((I([2:end,end],[2:end,end]) + I(:,[2:end,end]) - I(:,[1,1:end-1]) - I([2:end,end],[1,1:end-1]))./4).^2);GradientS_m = sqrt( (I - I(:,[1,1:end-1]) ).^2 + ((I([2:end,end],:) - I([1,1:end-1],:) + I([2:end,end],[1,1:end-1]) - I([1,1:end-1],[1,1:end-1]) )./4).^2 );GradientW_m = sqrt( ( I - I([1,1:end-1],:) ).^2 + ((I(:,[2:end,end]) - I(:,[1,1:end-1]) + I([1,1:end-1],[2:end,end]) - I([1,1:end-1],[1,1:end-1]))./4).^2 );GradientN_m = sqrt( ( ( I([2:end,end],:) - I([1,1:end-1],:) + I([2:end,end],[2:end,end]) - I([1,1:end-1],[2:end,end]))./4).^2 + ( I(:,[2:end,end]) - I ).^2 );GradientE = 1./(GradientE_m + epsilon);GradientS = 1./(GradientS_m + epsilon);GradientW = 1./(GradientW_m + epsilon);GradientN = 1./(GradientN_m + epsilon);Temp = GradientE + GradientS + GradientW + GradientN + lambda;Temp = Temp + (Temp == 0).*epsilon;I = (GradientE.*I([2:end,end],:) + GradientW.*I([1,1:end-1],:) + GradientN.*I(:,[2:end,end]) + GradientS.*I(:,[1,1:end-1]) + lambda*I).*M./Temp + (1-M).*I;end %if FlagNoise
end% switch dims
end% functionfunction v = getoptions(options, name, v, mendatory)% getoptions - retrieve options parameter
%
% v = getoptions(options, 'entry', v0);
% is equivalent to the code:
% if isfield(options, 'entry')
% v = options.entry;
% else
% v = v0;
% end
%
% Copyright (c) 2007 Gabriel Peyreif nargin<4mendatory = 0;
endif isfield(options, name)v = eval(['options.' name ';']);
elseif mendatoryerror(['You have to provide options.' name '.']);
end function [U,point_list] = grab_inpainting_mask(M, options)% grab_inpainting_mask - create a mask from user input
%
% U = grab_inpainting_mask(M, options);
%
% Select set of point in an image (useful to select a region for
% inpainting). The set of point is U==Inf.
%
% options.r is the radius for selection (default r=5).
%
% Selection stops with right click.
%
% Set options.mode='points' to gather disconnected points.
% Set options.mode='line' to gather connected lines.
%
% Copyright (c) 2006 Gabriel Peyreif nargin==3 && method==1U = grab_inpainting_mask_old(M, options);return;
endoptions.null = 0;
r = getoptions(options, 'r', 5);
method = getoptions(options, 'mode', 'points');if strcmp(method, 'line')if not(isfield(options, 'point_list'))[V,point_list] = pick_polygons(rescale(sum(M,3)),r);elsepoint_list = options.point_list;V = draw_polygons(rescale(sum(M,3)),r,point_list);end U = M; U(V==1) = Inf;return;
endm = size(M,1);
n = size(M,2);
s = size(M,3);U = sum(M,3)/3;
b = 1;
[Y,X] = meshgrid(1:n,1:m);
point_list = [];
while b==1clf;hold on;imagesc(rescale(M)); axis image; axis off; colormap gray(256);[y,x,b] = ginput(1);point_list(:,end+1) = [x;y];I = find((X-x).^2 + (Y-y).^2 <= r^2 );U(I) = Inf;for k=1:sMa = M(:,:,k);Ma(I) = 0;M(:,:,k) = Ma;end
end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sk = draw_polygons(mask,r,point_list)sk = mask*0;
for i=1:length(point_list)pl = point_list{i};for k=2:length(pl)sk = draw_line(sk,pl(1,k-1),pl(2,k-1),pl(1,k),pl(2,k),r);end
endfunction [sk,point_list] = pick_polygons(mask,r)% pick_polygons - ask for the user to build a set of curves
%
% sk = pick_polygons(mask,r);
%
% mask is a background image (should be in [0,1] approx).
%
% The user right-click on a set of point which create a curve.
% Left click stop a curve.
% Another left click stop the process.
%
% Copyright (c) 2007 Gabriel Peyren = size(mask,1);sk = zeros(n);
point_list = {};
b = 1;
while b(end)==1% draw a lineclf;imagesc(mask+sk); axis image; axis off;colormap gray(256);[y1,x1,b] = ginput(1);pl = [x1;y1];while b==1clf;imagesc(mask+sk); axis image; axis off;[y2,x2,c] = ginput(1);if c~=1if length(pl)>1point_list{end+1} = pl;endbreak;endpl(:,end+1) = [x2;y2];sk = draw_line(sk,x1,y1,x2,y2,r);x1 = x2; y1 = y2;end
end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sk = draw_line(sk,x1,y1,x2,y2,r)n = size(sk,1);
[Y,X] = meshgrid(1:n,1:n);
q = 80;
t = linspace(0,1,q);
x = x1*t+x2*(1-t); y = y1*t+y2*(1-t);
if r==0x = round( x ); y = round( y );sk( x+(y-1)*n ) = 1;
elsefor k=1:qI = find((X-x(k)).^2 + (Y-y(k)).^2 <= r^2 );sk(I) = 1;end
endfunction U = grab_inpainting_mask_old(M, r)% grab_inpainting_mask - create a mask from user input
%
% U = grab_inpainting_mask(M, r);
%
% r is the radius for selection (default r=5).
%
% Selection stops with right click.
%
% Copyright (c) 2006 Gabriel Peyr?if nargin<2r = 5;
endm = size(M,1);
n = size(M,2);
s = size(M,3);U = sum(M,3)/3;
b = 1;
[Y,X] = meshgrid(1:n,1:m);
while b==1clf;hold on;imagesc(rescale(M)); axis image; axis off; colormap gray(256);[y,x,b] = ginput(1);I = find((X-x).^2 + (Y-y).^2 <= r^2 );U(I) = Inf;for k=1:sMa = M(:,:,k);Ma(I) = 0;M(:,:,k) = Ma;end
end
function y = rescale(x,a,b)% rescale - rescale data in [a,b]
%
% y = rescale(x,a,b);
%
% Copyright (c) 2004 Gabriel Peyr?if nargin<2a = 0;
end
if nargin<3b = 1;
endm = min(x(:));
M = max(x(:));if M-m<epsy = x;
elsey = (b-a) * (x-m)/(M-m) + a;
end3. 模型效果
