1. 参考文献
2. TV图像修复模型
2.1 TV模型
% demo_TV.m% Author: HSW% Date: /3/25% HARBIN INSTITUTE OF TECHNOLOGY%% set matlabclose all;clear all;clc;options.null = 0;% read imageImg = 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 maskSetMask = 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);endelseif 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 FlagColorend% if SetMaskFlagNoise = 0; % default noiselessif 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 Imageif FlagColor == 1Positions = find(Mask(:,:,1) == 1);for channel = 1:3tmpnImg = nImg(:,:,channel);tmpnImg(Positions) = floor(255*rand(1,length(Positions))) + 1;InImg(:,:,channel) = tmpnImg;endelsePositions = find(Mask == 1);randValue = floor(255*rand(1,length(Positions))) + 1;InImg(Positions) = randValue;end% Main Total Variable ModelIterNum = 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)]);endend %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: /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 FlagNoiseend% switch dimsend% function
function 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) 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) 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;endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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);endendfunction [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) 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;endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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;endendfunction 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) 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;endend
function y = rescale(x,a,b)% rescale - rescale data in [a,b]%% y = rescale(x,a,b);%% Copyright (c) Gabriel Peyr?if nargin<2a = 0;endif nargin<3b = 1;endm = min(x(:));M = max(x(:));if M-m<epsy = x;elsey = (b-a) * (x-m)/(M-m) + a;end
3. 模型效果
