【图像修复】基于criminis算法实现图像修复附matlab代码

clear ;clc close all; %% 原始图像 srcImage= rgb2gray(imread(’bird12mask.bmp’)); % srcImage= imread(’lena.jpg’); srcImage= srcImage(150:375,50:298); % %% 破损图像 % damage_size_m=20; damage_size_n=20; %破损块大小 damageImage= srcImage; % damageImage(20:20+damage_size_m,125:125+damage_size_n)= 0; % damageImage= rgb2gray(imread(’pepper4mask.bmp’)); % damageImage= imread(’111..jpg’); % % figure % imshow(uint8(damageImage)); %破损图像 % imwrite(uint8(damageImage),’LLmask.jpg’) %% 等照度 damageImage=double(damageImage); [Gx, Gy] = gradient(damageImage); %图像梯度 Gx = Gx / 255; Gy = Gy / 255; temp = Gx; Gx = -Gy; Gy = temp; %梯度旋转90°为等照度 %% 初始化 % C 置信度项，为与图像大小相同的矩阵 % D 数据项，为与图像大小相同的矩阵 fillRegion = damageImage == 0; %fillRegion为待修复点为1，已存在点为0 sourceRegion = ~fillRegion; %sourceRegion待修复点为0，已存在点为1 sz = size(damageImage); ind = reshape(1:sz(1)*sz(2),sz(1),sz(2)); %找出各点的在从列中得到的坐标 C = double(sourceRegion); %待修复点为0，已存在点为1； D = repmat(-.1,sz); con=0;%填充次数初始化 %% 开始填充 tic %开始计时 while any(fillRegion(:)) %any函数：检测矩阵中是否有非零元素，如果有，则返回1，否则，返回0 %% 找外边界 dR = find(conv2(double(fillRegion),[1,1,1;1,-8,1;1,1,1],’same’)>0); %外 进行卷积以找到边界，dR为一组列向量，表示边界坐标 % dR = find(conv2(double(fillRegion),[1,1,1;1,-8,1;1,1,1],’same’)<0); %% 计算数据项，计算法向量，计算等照度 %更新等照度 [Gx, Gy] = gradient(double(damageImage)); %破损图梯度 Gx = Gx / 255; Gy = Gy / 255; temp = Gx; Gx = -Gy; Gy = temp; %计算法向量 sourceRegion = ~fillRegion; C1 = double(sourceRegion); [Nx, Ny] = gradient(C1); N = [Nx(dR(:)), Ny(dR(:))]; %找出边界点的梯度，在二值图像中梯度等于法线方向 N = normr(N); %求单位法向量 N(~isfinite(N))=0; % handle NaN and Inf（标准化梯度)求出边界点的法线方向 [Nm,Nn]=size(N); %将(1,1)归一化 for i=1:Nm for j=1:Nn if N(i,1)==1&N(i,2)==1 N(i,1)=sqrt(2)/2;N(i,2)=sqrt(2)/2; end end end %计算数据项D(P),只算边界上的 D(dR) = abs(Gx(dR).*N(:,1)+Gy(dR).*N(:,2)); %% 计算置信度项C(p) 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 %%计算 P = C(dR).*D(dR); %计算优先权P(p) %% 找到权值优先所对应的点p % 找到最大优先权值的块Hp，Hp为一9*9矩阵，矩阵值是待修复块坐标 [unused,ndx] = max(P); %返回P的最大值及P的位置 p = dR(ndx(1)); %p为待修复块的中心坐标点 [Hp,rows,cols] = getpatch(sz,p); %找到优先权最大的块，Hp为一9*9矩阵，矩阵值是待修复块坐标 toFill = fillRegion(Hp); %fillRegion为待修复点为1，已存在点为0，toFill为9*9的块 %% 找到最优匹配快 Hq为一9*9矩阵，矩阵值为最优匹配块的坐标 %找到最优匹配快，Hq为一9*9矩阵，矩阵值为最优匹配块的坐标 Hq = bestexemplar(damageImage, damageImage(rows,cols),fillRegion); %Hq = bestexemplar1(damageImage, damageImage(rows,cols), toFill’, sourceRegion); % Update fill region 更新破损区域 fillRegion(Hp(toFill)) = false; %Hp(toFill)表示只留下待修复点的坐标，使得待修复点的值为0 % Propagate confidence C(Hp(toFill)) = C(p); %被去掉那块的修复块的置信度被更新，C代表sourceRegion，待修复点为0，已存在点为1；用优先权最大的修复点的置信度来代替优先权最大的块的置信度 % Gx(Hp(toFill)) = Gx(Hq(toFill)); %用最优匹配块的向量梯度来更新优先权最大的块的梯度 % Gy(Hp(toFill)) = Gy(Hq(toFill)); % Copy image data from Hq to Hp ind(Hp(toFill)) = ind(Hq(toFill)); %最优匹配块的坐标来作为待修复块的坐标 damageImage(rows,cols) = ind2img(ind(rows,cols),damageImage); %ind(rows,cols)表示待修复块的坐标范围，该范围已被最优匹配块的值所替代 con=con+1; %% 保存图像 % imwrite(uint8(damageImage),strcat(’indoor2_’,num2str(con),’.bmp’)) end toc count=con %填充次数 OKImg= damageImage; figure, title(’原图’) subplot(121) imshow(uint8(srcImage)); %破损图像 title(’修复图’) subplot(122) imshow(OKImg, []) imwrite(uint8(OKImg),’bbGX.bmp’)