干涉图的延拓、criminism修复方法 的学习

目的：基于阶梯块的标定，需要对上下的条纹进行延拓

工具：matlab
参考文献：单小琴的硕士论文《单幅载波条纹图相位提取方法的研究》
和丁一飞的硕士论文《数字光栅投影测量关键技术研究》

重点：有效利用投影光栅条纹图像的条纹特征，即可信度值和等照度线

1、确定蒙版/掩膜，并将参数初始化

%% creatmask子函数
function [mask,counter1]=CreatMask(row,col,a,b,c,k)
counter=0;
mask=ones(row,col);
%掩膜的大小和原图像一般大
a=round(a);
b=round(b);
c=round(c);
k=round(k);
%a,b,c,k是对原图像使用imcrop后截取矩形的相关参数，详见imcrop的使用
for i=1:row
    for j=1:col
        if( a>j)||(j>(a+c))||( b>i)||(i>(b+k))
            mask(i,j)=0;
        else
            counter=counter+1;
        end
    end
end
%对矩形区域以外地方mask置为0，counter为确定像素点的数量
end

%% 主函数的第一部分：读图，确定掩膜图
w=double(imread('*.bmp'))/255;%读入图像，可以是JPG，BMP格式


figure;imshow(w,[]); 
title('原始图像');
[row,col]=size(w);
c=round(row/2);
bGrey = w(c,:);%c行的灰度值，可以看光栅的正弦性如何
figure;
plot(bGrey);title('原图x轴方向的相位分布');
%w=I(57:960,54:955);%取有效数据点矩阵


m1=255;
m2=255;
img= imresize(w, [m1,m2],'bicubic');
%适当将图像放缩，否则处理起来很耗时间，最后再将其大小还原
figure;imshow(img,[]);
title('放缩图');
 
rect3=[1 1 40 60];
a=rect3(1);
b=rect3(2);
c=rect3(3);
k=rect3(4);
%自己创建的参数，也可以直接在原图上截取优质的区域进行延拓
mask1 = CreatMask(m1,m2,a,b,c, k);%创建有效数据区域mask，去除边缘条纹数据

resourceInMask =img.* mask1;%对原图像进行掩膜
figure;imshow(resourceInMask,[]);
title('mask掩膜图');

%% 主函数的第二部分：参数初始化
fillImg=resourceInMask; 
fillRegion = fillImg == 0;%fillRegion为待修复点为1，已存在点为0
fillImg = double(fillImg);
origImg = fillImg;%都是待修复图像
s = size(fillImg);
ind = reshape(1:s(1)*s(2),s(1),s(2));%找出各点的在从列中得到的坐标重新调整矩阵的行数、列数、维数。
sz = [size(fillImg,1) size(fillImg,2)];
sourceRegion = ~fillRegion;%待修复点为0，已存在点为1；
% 初始化等照度线值
[Ix, Iy] = gradient(fillImg);%Gradient(F)函数求的是数值上的梯度，假设F为矩阵.
Ix = Ix ./ (255*sqrt(Ix.^2 + Iy.^2));
Ix(~isfinite(Ix))=0;%isFinite() 函数用于检查其参数是否是无穷大
Iy = Iy ./ (255*sqrt(Ix.^2 + Iy.^2));
Iy(~isfinite(Iy))=0; 
temp = Ix; Ix = -Iy; Iy = temp; % 等照度线方向，旋转梯度90度
% 初始化置信项和数据项
C = double(sourceRegion);%待修复点为0，已存在点为1；
D = repmat(-.1,sz);%repmat函数主要是用于快速的产生一个大的矩阵
% 'rand' 用于产生可重复的随机数，便于测试
rand('state',0);%以后产生的随机数都与第一次运行产生的相同

2、循环进行延拓

%% 主函数的第三部分
% d循环，直到整个填充区域被覆盖
while any(fillRegion(:))
    % 找到填充区域的轮廓和归一化梯度
    dR = find(conv2(double(fillRegion),[1,1,1;1,-8,1;1,1,1],'same')>0);%进行卷积以找到边界，dR为一组列向量，表示边界坐标
    
    [Nx, Ny] = gradient(C);
    N = [Nx(dR(:)), Ny(dR(:))];%找出边界点的梯度，在二值图像中梯度等于法线方向
    N = normr(N); %normr(X)接受单个矩阵或矩阵的单元数组，并返回行规范化为1的矩阵。
    N(~isfinite(N))=0; % h处理NaN和Inf（标准化梯度)求出边界点的法线方向
    % 计算沿着边界点的置信度
    for k = dR'
        Hp = getpatch(sz,k);%Hp为一9*9矩阵，矩阵值是待修复块坐标.返回以像素p为中心的9x9补丁的索引 ，getpatch是子函数，后面再看
        q = Hp(~(fillRegion(Hp)));%q为一个列向量矩阵,为该修复块已知像素点的坐标
        C(k) = sum(C(q))/numel(Hp);%计算置信度
    end
    % 计算补丁优先级=置信项*数据项
    D(dR) = abs(Ix(dR).*N(:,1)+Iy(dR).*N(:,2))+0.001;%计算的是边界点的梯度值，而不是待修复块的梯度最大值
    priorities = C(dR).* D(dR);
    % 找到最高优先级的补丁，Hp
    [unused,ndx] = max(priorities);
    p = dR(ndx(1));%p为待修复块的中心坐标点
    [Hp,rows,cols] = getpatch(sz,p);%找到优先权最大的块，Hp为一9*9矩阵，矩阵值是待修复块坐标
    toFill = fillRegion(Hp);%fillRegion为待修复点为1，已存在点为0，toFill为待修复块，也为待修复点为1，已存在点为0
    %找到最小化错误的范例，Hq
    Hq = bestexemplar(origImg,origImg(rows,cols),toFill',sourceRegion);%找到最优匹配快，Hq为一9*9矩阵，矩阵值为最优匹配块的坐标
    % 更新填充区域
    fillRegion(Hp(toFill)) = false; %Hp(toFill)表示只留下待修复点的坐标，使得待修复点的值为0
    
    % 传递置信值和等照度线值
    C(Hp(toFill))  = C(p); %被去掉那块的修复快的置信度被更新，C代表sourceRegion，待修复点为0，已存在点为1；用优先权最大的修复点的置信度来代替优先权最大的块的置信度
    Ix(Hp(toFill)) = Ix(Hq(toFill));%用最优匹配块的向量梯度来更新优先权最大的块的梯度
    Iy(Hp(toFill)) = Iy(Hq(toFill));
    % 从 Hq 给 Hp复制图像数据
    ind(Hp(toFill)) = ind(Hq(toFill));%最优匹配块的坐标来作为待修复块的坐标
    origImg(rows,cols) = ind2img(ind(rows,cols),origImg);%ind(rows,cols)表示待修复块的坐标范围，该范围已被最优匹配块的值所替代
end

3、子函数的补充

%  getpatch(sz,p)返回以像素p为中心的9x9补丁的索引。
function [Hp,rows,cols] = getpatch(sz,p)
% [x,y] = ind2sub(sz,p);  % 2*w+1 ==补丁大小
w=4; p=p-1; y=floor(p/sz(1))+1; p=rem(p,sz(1)); x=floor(p)+1;
rows = max(x-w,1):min(x+w,sz(1));
cols = (max(y-w,1):min(y+w,sz(2)))';
Hp = sub2ndx(rows,cols,sz(1));%将(rows,cols)订阅样式的索引转换为Matlab索引样式的索引。

function img2 = ind2img(ind,img)
    img2 = img(ind);%只是在该范围内的img值被赋予到img2

% 将(rows,cols)订阅样式的索引转换为Matlab索引样式的索引。不幸的是，“sub2ind”不能用于此。
function N = sub2ndx(rows,cols,nTotalRows)
X = rows(ones(length(cols),1),:);
Y = cols(:,ones(1,length(rows)));
N = X+(Y-1)*nTotalRows;

%扫描整个图像(使用滑动窗口)，寻找错误最小的范例。
%调用MEX函数bestexemplarhelper1（）。这个函数没办法展示出来。
function Hq = bestexemplar(img,Ip,toFill,sourceRegion)
m=size(Ip,1); mm=size(img,1); n=size(Ip,2); nn=size(img,2);
best = bestexemplarhelper1(mm,nn,m,n,img,Ip,toFill,sourceRegion);
Hq = sub2ndx(best(1):best(2),(best(3):best(4))',mm);

4、最后一部分，看现象

inpaintedImg=imresize(origImg, [row,col],'bicubic');
bGrey = inpaintedImg(c,:);%c行的灰度值
figure;
plot(bGrey);title('延拓图x轴方向的相位分布');
filename=['1.bmp']; 
imwrite(inpaintedImg,filename,'bmp');%对图像进行存储
figure, imshow(inpaintedImg, []);
title('条纹延拓图');