1. 参考文献
2. FFM模型实现
2.1 FFM模型代码
%%%%%%%%%%%%%%%%%%%%%%%%%%% Demo of Fast March Method %%%%%%%%%%%%%%%%%%%%%% Date: -10-23% author: HSW% HARBIN INSTITUTE OF TECHNOLOGY%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%close allclear allclc;% 读入图像BAND = 1;KNOWN = 2;INSIDE = 3;% 读入图像I = imread('C10.bmp');M = imread('mask2.png');% 彩色图像if size(M,3) == 3 M = rgb2gray(M);M = M>0;endfigure; imshow(I, []); title('原图像'); I = im2double(I);Temp(:,:,1) = M;Temp(:,:,2) = M;Temp(:,:,3) = M;I = (1-Temp).*I; %+ Temp.*0.5;% % 灰度图像% if ndims(M)%M = rgb2gray(M);%I = rgb2gray(I);% end% M = M<1;% I = double(I);% Temp = M;% I = (1-Temp).*I + Temp.*255;figure;imshow(I,[]);title('原图像');% 标记[m,n,dims] = size(I);M1 = (bwdist(M,'cityblock') <= 1);%初始化flag值flag = INSIDE*M + KNOWN*(1-M);flag =(1-M1+M).*flag+BAND*(M1-M);%初始化 TT = zeros([m,n]) + 1e6.*M;IterMax = 10;for iter = 1:IterMax[I] = fast_march_method(I,T,flag);endfigure;imshow(I,[]);title('修复结果');
2.2 FFM模型
function [I,T,flag] = fast_march_method(I,T,flag)% input:%I: 待修复图像%M: 缺损区域mask，缺损区域 = 1%T: 距离标记矩阵%flag:标记矩阵，BAND/KNOWN/INSIDE% output:% I: 修复图像%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%BAND = 1;KNOWN = 2;% INSIDE = 3;[m,n,dims] = size(I);iter = 0;while find(flag == BAND) %存在待修复点% 获取T最小的点的坐标xposition = 1;yposition = 1;Tmin = 1e6;for i = 1:mfor j = 1:nif (flag(i,j) == BAND) && (T(i,j) < Tmin )xposition = i;yposition = j;Tmin = T(i,j);endendendflag(xposition,yposition) = KNOWN;switch dimscase 3 %彩色图像[I,T,flag] = four_solve(I,T,flag,xposition,yposition,3); %更新四领域，注意三个通道要同时处理otherwise %灰度图像[I,T,flag] = four_solve(I,T,flag,xposition,yposition,1); %更新四领域end %switchiter = iter + 1;if mod(iter,2000) == 0figure;imshow(I,[]);endend %whileend %functionfunction [I,T,flag] = four_solve(I,T,flag,xposition,yposition,channel)% input:%xposition: x坐标%yposition: y坐标%T：距离标记矩阵%flag: 标记BAND,KNOWN,INSIDE%channel：图像通道% output:%T: 距离标记矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%pSize = 1;BAND = 1;KNOWN = 2;INSIDE = 3;FourNeighbor = [xposition-1,yposition; xposition,yposition-1; xposition+1,yposition; xposition,yposition+1];for i = 1:4x = FourNeighbor(i,1);y = FourNeighbor(i,2);if ( flag(x,y) ~= KNOWN)if (flag(x,y) == INSIDE)flag(x,y) = BAND;I = fast_march_method_inpaint(I,T,flag,x,y,channel,pSize); %进行修复endT(x,y) = min(T(x,y),min(four_solve_min(T,flag,x-1,y,x,y-1),min(four_solve_min(T,flag,x+1,y,x,y-1),min(four_solve_min(T,flag,x-1,y,x,y+1),four_solve_min(T,flag,x+1,y,x,y+1))))); %更新距离endend%for iend% functionfunction sol = four_solve_min(T,flag,x1,y1,x2,y2)% input:%坐标点（x1,y1),(x2,y2)% output:%sol: 距离%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%KNOWN = 2;sol = 100000;if flag(x1,y1) == KNOWNif flag(x2,y2) == KNOWNr = sqrt( abs(2 - ((T(x1,y1) - T(x2,y2))^2)));s = (T(x1,y1) + T(x2,y2) - r)/2;if (s >= T(x1,y1)) && (s >= T(x2,y2))sol = s;elses = s + r;if ( (s >= T(x1,y1)) && (s >= T(x2,y2)) )sol = s;endendelsesol = 1 + T(x1,y1);endelseif flag(x2,y2) == KNOWNsol = 1 + T(x1,y1);end %ifend %functionfunction I = fast_march_method_inpaint(I,T,flag,x,y,channel,pSize)% input:% I: 待修复图像% x: 待修复点x坐标% y: 待修复点y坐标% channel: 颜色通道% pSize: 利用领域块参数（2*pSize + 1) X (2pSize + 1)% output：% I：修复图像%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%KNOWN = 2;INSIDE = 3;sum1 = 0;sum2 = 0;sum_Ex0n=0;sum_Ex0d=0;sum_Ey0n=0;sum_Ey0d=0;dims = channel;switch dimscase 1 %灰度图像%求等照度线方向for k = x - pSize:x+pSizefor l = y-pSize:y+pSizeif flag(k,l) == KNOWNGradientIx = I(k,l) - I(k-1,l-1);GradientIy = I(k-1,l) - I(k,l-1);IsoIx = -GradientIy;IsoIy = GradientIx;d0=1/(1+sqrt((k-x)^2+(l-y)^2));sum_Ex0n=sum_Ex0n+d0*IsoIx;sum_Ex0d=sum_Ex0d+d0;sum_Ey0n=sum_Ey0n+d0*IsoIy;sum_Ey0d=sum_Ey0d+d0;Dst(k,l) = d0;%式（2）IsoIx = IsoIx/(1+sqrt(IsoIx^2+IsoIy^2));IsoIy = IsoIy/(1+sqrt(IsoIx^2+IsoIy^2));Dirctx(k,l) = IsoIx;Dircty(k,l) = IsoIy;endendend% (x,y)的等照度线方向。Ex0=sum_Ex0n/sum_Ex0d;Ey0= sum_Ey0n/sum_Ey0d;% (x,y)归一化Ex0=Ex0/(1+sqrt(Ex0^2+Ey0^2));Ey0=Ey0/(1+sqrt(Ex0^2+Ey0^2));%估值I(x,y)for k = x - pSize:x+pSizefor l = y-pSize:y+pSizeif flag(k,l) == KNOWNQ= exp(Ex0*Dirctx(k,l)+Ey0* Dircty(k,l));% 式（1）W(k,l)=Q*Dst(k,l)*(1/(1+T(k,l)));% if flag(k+1,l) == KNOWN && flag(k-1,l) == KNOWN && flag(k,l+1) == KNOWN && flag(k,l-1) == KNOWN%sum1 = sum1+(I(k,l) + (I(k+1,l) - I(k-1,l))*(k-x) + (I(k,l+1) - I(k,l-1))*(l-y))*W(k,l);% elseif flag(k+1,l) == KNOWN && flag(k-1,l) == KNOWN%sum1 = sum1 + ( I(k,l) + (I(k+1,l) - I(k-1,l))*(k-x) + (I(k-1,l) - I(k,l-1))*(l-y))*W(k,l);% elseif flag(k,l+1) == KNOWN && flag(k,l-1) == KNOWN%sum1 = sum1 + ( I(k,l) + (I(k,l) - I(k-1,l-1))*(k-x) + (I(k,l+1) - I(k,l-1))*(l-y))*W(k,l);% else%sum1 = sum1 + I(k,l)*W(k,l);% endsum1 = sum1 + I(k,l)*W(k,l);sum2=sum2+W(k,l);endendendI(x,y)=sum1/sum2;case 3for channel = 1:dims%求等照度线方向for k = (x - pSize):(x+pSize)for l = (y-pSize):(y+pSize)if flag(k,l) == KNOWNGradientIx = I(k,l,channel) - I(k-1,l-1,channel);GradientIy = I(k-1,l,channel) - I(k,l-1,channel);IsoIx = -GradientIy;IsoIy = GradientIx;d0=1/(1+sqrt((k-x)^2+(l-y)^2));sum_Ex0n = sum_Ex0n+d0*IsoIx;sum_Ex0d = sum_Ex0d+d0;sum_Ey0n = sum_Ey0n+d0*IsoIy;sum_Ey0d = sum_Ey0d+d0;Dst(k,l) = d0;IsoIx = IsoIx./(0.1+sqrt(IsoIx^2+IsoIy^2));IsoIy = IsoIy./(0.1+sqrt(IsoIx^2+IsoIy^2));Dirctx(k,l) = IsoIx;Dircty(k,l) = IsoIy;endendend%进行估值Ex0 = sum_Ex0n/sum_Ex0d;% (i,j)的等照度线方向。Ey0 = sum_Ey0n/sum_Ey0d;Ex0 = Ex0/(0.1+sqrt(Ex0^2+Ey0^2));Ey0 = Ey0/(0.1+sqrt(Ex0^2+Ey0^2));for k = (x - pSize):(x+pSize)for l = (y-pSize):(y+pSize)if flag(k,l) == KNOWNQ = exp(Ex0*Dirctx(k,l)+Ey0* Dircty(k,l));% 式（1）W = Q*Dst(k,l)*(1/(1+abs(T(k,l) - T(x,y)) ));%if flag(k+1,l) == KNOWN && flag(k-1,l) == KNOWN && flag(k,l+1) == KNOWN && flag(k,l-1) == KNOWN% sum1 = sum1+(I(k,l,channel) + (I(k+1,l,channel) - I(k-1,l,channel))*(k-x) + (I(k,l+1,channel) - I(k,l-1,channel))*(l-y))*W(k,l);%elseif flag(k+1,l) == KNOWN && flag(k-1,l) == KNOWN% sum1 = sum1 + ( I(k,l,channel) + (I(k+1,l,channel) - I(k-1,l,channel))*(k-x))*W(k,l);%elseif flag(k,l+1) == KNOWN && flag(k,l-1) == KNOWN% sum1 = sum1 + ( I(k,l,channel) + (I(k,l+1,channel) - I(k,l-1,channel))*(l-y))*W(k,l);%else% sum1 = sum1 + I(k,l,channel)*W(k,l);%endsum1 = sum1 + I(k,l,channel)*W;sum2=sum2+W;end%if endendI(x,y,channel) = sum1/sum2; end %for channelotherwisedisplay('ERROR !');end %switchend %function
3. 模型效果
