1. 参考文献
2. BSCB模型代码
2.1 BSCB模型demo
% demo_BSCB.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 = 1;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 = 1;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 SetMasknImg = (1-Mask).*Img;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 BSCB ModelIterNum = 3000;I = BSCB_Diffusion(InImg,FlagColor,Mask,0.1);for iter = 1:IterNumI = BSCB_Inpainting(I,FlagColor,Mask,0.2);if mod(iter,500) == 0 figure; imshow(I/255,[]); title(['Results of IterNum = ', num2str(iter)]); I = BSCB_Diffusion(I,FlagColor,Mask,0.2);end endI = 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 BSCB图像修复模型实现
function Img = BSCB_Inpainting(I,FlagColor,M,delta_t)% input:% I: 待修复图像% M：缺损区域mask,缺损区域取值为1 % FlagColor:确定是否为彩色图像% delta_t:时间差分% output:%Img: 修复图像% Author： HSW% Date: /3/25% HARBIN INSTITUTE OF TECHNOLOGY%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Img = zeros(size(I));if FlagColor == 1for channel = 1:3Ix = (I([2:end,end],:,channel) - I([1,1:end-1],:,channel))./2;Iy = (I(:,[2:end,end],channel) - I(:,[1,1:end-1],channel))./2;Ixx = I([2:end,end],:,channel) + I([1,1:end-1],:,channel) - 2*I(:,:,channel);Iyy = I(:,[2:end,end],channel) + I(:,[1,1:end-1],channel) - 2*I(:,:,channel);Laplace = Ixx + Iyy;Laplacex = (Laplace([2:end,end],:) - Laplace([1,1:end-1],:))./2;Laplacey = (Laplace(:,[2:end,end]) - Laplace(:,[1,1:end-1]))./2;Temp = sqrt(Ix.^2 + Iy.^2 + 0.00001);normNx = (-Iy)./Temp;normNy = Ix./Temp;Beta = Laplacex.*normNx + Laplacey.*normNy;dist1 = min(0,I(:,:,channel) - I([1,1:end-1],:,channel)).^2 + max(0, I([2:end,end],:,channel) - I(:,:,channel) ).^2 ...+ min(0, I(:,:,channel) - I(:,[1,1:end-1],channel)).^2 + max(0, I(:,[2:end,end],channel) - I(:,:,channel)).^2;dist2 = max(0, I(:,:,channel) - I([1,1:end-1],:,channel)).^2 + min(0, I([2:end,end],:,channel) - I(:,:,channel)).^2 ...+ max(0, I(:,:,channel) - I(:,[1,1:end-1],channel)).^2 + min(0,I(:,[2:end,end],channel) - I(:,:,channel)).^2;MuD1 = sqrt(dist1);MuD2 = sqrt(dist2);MuD = Beta.*((Beta > 0).* MuD1 + (1 - (Beta > 0)).* MuD2);signMuD = sign(MuD); Img(:,:,channel) = I(:,:,channel) + delta_t*M(:,:,channel).*signMuD.*sqrt(sqrt(signMuD.*MuD));end% for channel elseIx = (I([2:end,end],:) - I([1,1:end-1],:))./2;Iy = (I(:,[2:end,end]) - I(:,[1,1:end-1]))./2;Ixx = I([2:end,end],:) + I([1,1:end-1],:) - 2*I;Iyy = I(:,[2:end,end]) + I(:,[1,1:end-1]) - 2*I;Laplace = Ixx + Iyy;Laplacex = (Laplace([2:end,end],:) - Laplace([1,1:end-1],:))./2;Laplacey = (Laplace(:,[2:end,end]) - Laplace(:,[1,1:end-1]))./2;Temp = sqrt(Ix.^2 + Iy.^2 + 0.00001);normNx = -Iy./Temp;normNy = Ix./Temp;Beta = Laplacex.*normNx + Laplacey.*normNy;dist1 = min(0,I - I([1,1:end-1],:)).^2 + max(0, I([2:end,end],:) - I ).^2 ...+ min(0, I - I(:,[1,1:end-1])).^2 + max(0, I(:,[2:end,end]) - I).^2;dist2 = max(0, I - I([1,1:end-1],:)).^2 + min(0, I([2:end,end],:) - I).^2 ...+ max(0, I - I(:,[1,1:end-1])).^2 + min(0,I(:,[2:end,end]) - I).^2;MuD1 = sqrt(dist1);MuD2 = sqrt(dist2);MuD = Beta.*(( Beta > 0 ).* MuD1 + ( 1 - (Beta > 0) ).* MuD2);signMuD = sign(MuD); Img = I + delta_t*M.*signMuD.*sqrt(sqrt(signMuD.*MuD));%Img = I + delta_t*M.*MuD; %I = M.*I + (1-M).*Img; end% if FlagColor == 1end % function
2.3 BSCB扩散模型
function Img = BSCB_Diffusion(I,FlagColor,M,delta_t)% input:% I: 待修复图像% M： 缺损区域mask,缺损区域取值为1% FlagColor: 标记是否为彩色% delta_t: 时间差分% output：% Img: 扩散后的图像% Author： HSW% Date: /3/25% HARBIN INSTITUTE OF TECHNOLOGY%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Img = zeros(size(I));if FlagColor == 1dims = size(I,3);for channel = 1:dimsIx = ( I([2:end,end],:,channel) - I([1,1:end-1],:,channel) )./2;Iy = ( I(:,[2:end,end],channel) - I(:,[1,1:end-1],channel) )./2;Ixx = I([2:end,end],:,channel) + I([1,1:end-1],:,channel) - 2*I(:,:,channel);Iyy = I(:,[2:end,end],channel) + I(:,[1,1:end-1],channel) - 2*I(:,:,channel);Ixy = ((I([2:end,end],[2:end,end],channel) + I([1,1:end-1],[2:end,end],channel))./2 - (I([2:end,end],[1,1:end-1],channel) + I([1,1:end-1],[1,1:end-1],channel))./2)./2;TempNorm = Ix.^2 + Iy.^2 + 0.00001;DMu = (Ixx.*(Iy).^2 + Iyy.*(Ix).^2 - 2*Ix.*Iy.*Ixy)./TempNorm;Img(:,:,channel) = I(:,:,channel) + delta_t*M(:,:,channel).*DMu;endelseIx = ( I([2:end,end],:) - I([1,1:end-1],:) )./2;Iy = ( I(:,[2:end,end]) - I(:,[1,1:end-1]) )./2;Ixx = I([2:end,end],:) + I([1,1:end-1],:) - 2*I;Iyy = I(:,[2:end,end]) + I(:,[1,1:end-1]) - 2*I;Ixy = ((I([2:end,end],[2:end,end]) + I([1,1:end-1],[2:end,end]))./2 - (I([2:end,end],[1,1:end-1]) + I([1,1:end-1],[1,1:end-1]))./2)./2;TempNorm = Ix.^2 + Iy.^2 + 0.00001;DMu = (Ixx.*(Iy).^2 + Iyy.*(Ix).^2 - 2*Ix.*Iy.*Ixy)./TempNorm;Img = I + delta_t*M.*DMu;end %if FlagColorend %function
2.4 检索选项参数
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
2.5 创建掩模Mask
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
3. 模型效果
