www.pudn.com > Criminisi-original-gray.zip > inpaint.m, change:2012-05-18,size:7868b


function [inpaintedImg,origImg,fillImg,C,D,fillMovie] = inpaint(imgFilename,fillFilename)
%INPAINT  Exemplar-based inpainting.
%
% Usage:   [inpaintedImg,origImg,fillImg,C,D,fillMovie] ...
%                = inpaint(imgFilename,fillFilename,fillColor)
% Inputs: 
%   imgFilename    原始图像.   
%   fillFilename   指定了修复区域的图像. 
%   fillColor      为了指定修复区域而由RGB三分量表示的一个颜色
% Outputs:
%   inpaintedImg   双精度度的M*N*3的修复好了的图像. 
%   origImg        双精度度的M*N*3的原始图像.
%   fillImg        双精度度的M*N*3的填充区域图像.
%   C              各次迭代中,M*N的置信度值的矩阵.
%   D              各次迭代中,M*N的数据项值的矩阵.
%   fillMovie      展示每次迭代后的填充区域.. 
%
% Example:
%   [i1,i2,i3,c,d,mov] = inpaint('bungee0.png','bungee1.png',[0 255 0]);
%  plotall;           % quick and dirty plotting script
%  close; movie(mov); % grab some popcorn 
%  author: Sooraj Bhat

warning off MATLAB:divideByZero %关闭警告
tic;
[img,fillImg_org,fillRegion,number_of_inpainted] = loadimgs(imgFilename,fillFilename);
imshow(fillRegion);
img = double(img);
origImg = img;  %双精度的原图像
ind = img2ind(img);  %调用img2ind()函数
sz = [size(img,1) size(img,2)];%sz=[img行数 img列数]
sourceRegion = ~fillRegion;
fillImg=fillImg_org;

% Initialize isophote values
[Ix(:,:) Iy(:,:)] = gradient(img(:,:));
Ix= Ix/255; Iy=Iy/255;
temp = Ix; Ix = -Iy; Iy = temp;  % Rotate gradient 90 degrees

% Initialize confidence and data terms
C = double(sourceRegion);   
D = repmat(-.1,sz);   %重复矩阵

% Visualization stuff
if nargout==6    %nargout输出参数的个数
  fillMovie(1).cdata=uint8(img); 
  fillMovie(1).colormap=[];
%   origImg(1,1,:) = fillColor;
  iter = 2;
end

% Seed 'rand' for reproducible results (good for testing)
rand('state',0);   %Resets the generator to its initial state.

img=fillImg_org;
% Loop until entire fill region has been covered
fill=find(fillRegion==1);

while any(fillRegion(:)) %any()非零时为真
     % Find contour & normalized gradients of fill region
     dR = find(conv2(single(fillRegion),[1,1,1;1,-8,1;1,1,1],'same')>0);  %single()变成单精度的,                               
                                                                          %conv2(a,b,'same')二维卷积,返回和a一样大小的矩阵;
                                                                          %find()返回满足条件的矩阵中元素的序列号列向量,按列
                                                                          %fillRegion原图大小 
                                                                     
                                                                       
     [Nx,Ny] = gradient(1-fillRegion);
     N = [Nx(dR(:)) Ny(dR(:))];
     % imshow(N);
     N = normr(N);     %按行归一化 
     N(~isfinite(N))=0; % 处理 NaN and Inf ,isfinite()返回向量,元素为有限则对应1,元素为无限则对应0
  
     % Compute confidences along the fill front
     for k=dR'
       Hp = getpatch(sz,k);      %调用getpatch()函数
       q = Hp(~(fillRegion(Hp))); %fillRegion(Hp)是9*9的块,里面待填充区域是1,源区域是0
                                 %q是列向量,由源区域像素点构成
       C(k) = sum(C(q))/numel(Hp);%numel()面积即像素数
     end
  
     % Compute patch priorities = confidence term * data term
     D(dR) = abs(Ix(dR).*N(:,1)+Iy(dR).*N(:,2)) + 0.00001;
     priorities = C(dR).* D(dR);
  
     % Find patch with maximum priority, Hp
     [unused,ndx] = max(priorities(:));% max()找到每列中的最大值放入unused,第几个放入ndx
     
     p = dR(ndx(1));
     [Hp,rows,cols] = getpatch(sz,p);
     toFill = fillRegion(Hp);%9*9二值块
  
     % Find exemplar that minimizes error, Hq
%      img1(:,:,1)=img;
%      img1(:,:,2)=img;
%      img1(:,:,3)=img;
     [X,map]=gray2ind(img);
     img1=ind2rgb(X,map);

%      imshow(uint8(img1),[]);title('aa');
     Hq = bestexemplar(img1,img1(rows,cols,:),toFill',sourceRegion);    %调用bestexemplar()
     %qq = Hp(~(fillRegion(Hp)));
     % Update fill region
     fillRegion(Hp(toFill)) = false;
  
     % Propagate confidence & isophote values
      
     C(Hp(toFill))  = C(p);

     Ix(Hp(toFill)) = Ix(Hq(toFill));
     Iy(Hp(toFill)) = Iy(Hq(toFill));
  
     % Copy image data from Hq to Hp
    ind(Hp(toFill)) = ind(Hq(toFill));
    img(rows,cols) = ind2img(ind(rows,cols),origImg);  
    figure(8);
    imshow(uint8(img));

    % Visualization stuff
%     if nargout==6
%        ind2 = ind;
%        ind2(fillRegion) = 1;
%        fillMovie(iter).cdata=uint8(ind2img(ind2,origImg)); 
%        fillMovie(iter).colormap=[];
%     end
    iter = iter+1;

inpaintedImg=double(img);
fill_mse=mse(origImg(fill(:))-inpaintedImg(fill(:)))
fill_all=mse(origImg-inpaintedImg)
psnr_img=10*log10(255^2/fill_all)


end
inpaintedImg=double(img);
fill_mse=mse(origImg(fill(:))-inpaintedImg(fill(:)))
fill_all=mse(origImg-inpaintedImg)
psnr_img=10*log10(255^2/fill_all)

toc;

%---------------------------------------------------------------------
% Scans over the entire image (with a sliding window)
% for the exemplar with the lowest error. Calls a MEX function.
%---------------------------------------------------------------------
function Hq = bestexemplar(img,Ip,toFill,sourceRegion)
m=size(Ip,1);n=size(Ip,2);mm=size(img,1); nn=size(img,2);
best = bestexemplarhelper(mm,nn,m,n,img,Ip,toFill,sourceRegion);
Hq = sub2ndx(best(1):best(2),(best(3):best(4))',mm);


%---------------------------------------------------------------------
% Returns the indices for a 9x9 patch centered at pixel p.
%---------------------------------------------------------------------
function [Hp,rows,cols] = getpatch(sz,p)
% [x,y] = ind2sub(sz,p);  % 2*w+1 == the patch size
w=3; p=p-1; y=floor(p/sz(1))+1; p=rem(p,sz(1)); x=floor(p)+1;  %floor(x)小于x的最大的整数,rem(x,y)=x - n.*y, n=fix(x./y)=floor(x./y)
                                                               %x是p的行数,y是p的列数

rows = max(x-w,1):min(x+w,sz(1));          %9*9的块 rows是行向量
cols = (max(y-w,1):min(y+w,sz(2)))';       %cols是列向量
Hp = sub2ndx(rows,cols,sz(1));     %调用sub2ndx()


%---------------------------------------------------------------------
% Converts the (rows,cols) subscript-style indices to Matlab index-style
% indices.  Unforunately, 'sub2ind' cannot be used for this.
%---------------------------------------------------------------------
function N = sub2ndx(rows,cols,nTotalRows)
X = rows(ones(length(cols),1),:);% ones()是length(cols)*1的全1列向量
%imshow(X);title('x');          %X是9*9矩阵,每一列的元素都一样,都是行数
%[yy,ee]=size(X)
Y = cols(:,ones(1,length(rows)));%Y是9*9矩阵,每一行的元素都一样,都是列数
N = X+(Y-1)*nTotalRows;  %索引矩阵


%---------------------------------------------------------------------
% Converts an indexed image into an RGB image, using 'img' as a colormap
%---------------------------------------------------------------------
function img2 = ind2img(ind,img)
img2=img(ind); 


%---------------------------------------------------------------------
% Converts an RGB image into a indexed image, using the image itself as
% the colormap.
%---------------------------------------------------------------------
function ind = img2ind(img)
s=size(img); ind=reshape(1:s(1)*s(2),s(1),s(2));%1到s(1)*s(2)按列给出s(1)行s(2)列的矩阵,按列赋值


%---------------------------------------------------------------------
% Loads the an image and it's fill region, using 'fillColor' as a marker
% value for knowing which pixels are to be filled.
%---------------------------------------------------------------------
function [img,fillImg,fillRegion,number_of_inpainted]= loadimgs(imgFilename,fillFilename)
img = imread(imgFilename);a= imread(fillFilename); fillImg=a(:,:,1);
fillRegion = fillImg(:,:)==255;
a=single(fillRegion);
number_of_inpainted=sum(sum(a))
mse_broken=mse(double(img)-double(fillImg))
psnr_broken=10*log10(255^2/mse_broken)