锡膏高度判断（使用五步相位移法）--matlab

线路问题进行到一半的时候，被分配着去做了锡膏的3D检测算法。趁着两者都有了些许进展，终于是可以来记录一下相关的情况了~

由于解相位算法不是我做的，直接拿了同事的成果，这里就不讲了。主要是这算法也存在着问题，并且有个部分不知道怎么得到的。。。汗

    计算高度的方法总结就是：根据解相位的图，拟合出锡膏周围基准面，再进行计算得到高度。

遇到的问题

1：解出来的相位在边界处，沿x方向在左边界和右边界存在着翘曲的情况，这导致我在后续处理中遇到了一些困难

    解决方法：将处理的图片沿x和y方向在最早需要拓展的像素值上再多拓展了20个点，后续进行计算高度前将之裁掉。

2：图片锡膏周围存在阴影处，解出来的相位基本是近乎水平但是乱序波动的情况。

    解决方法：可以根据灰度判断出阴影的具体位置，根据解相位出来的规律，在左边界存在的阴影所有点都取来作为拟合点、右边界所有点都不取，最后得到一个拟合的基准面。

function  X = Ls1(~)
 

 

% Ref1 = imread('E:\Desktop\3DSPI\pmp\2\4\140_20_0.bmp');
% Ref2 = imread('E:\Desktop\3DSPI\pmp\2\4\140_20_1.bmp');
% Ref3 = imread('E:\Desktop\3DSPI\pmp\2\4\140_20_2.bmp');
% Ref4 = imread('E:\Desktop\3DSPI\pmp\2\4\140_20_3.bmp');
% Ref5 = imread('E:\Desktop\3DSPI\pmp\2\4\140_20_4.bmp');

%边缘全阴影读取
% Ref1 = imread('E:\Desktop\3DSPI\pmp\2\3\130_20_0.bmp');
% Ref2 = imread('E:\Desktop\3DSPI\pmp\2\3\130_20_1.bmp');
% Ref3 = imread('E:\Desktop\3DSPI\pmp\2\3\130_20_2.bmp');
% Ref4 = imread('E:\Desktop\3DSPI\pmp\2\3\130_20_3.bmp');
% Ref5 = imread('E:\Desktop\3DSPI\pmp\2\3\130_20_4.bmp');

Ref1 = imread('E:\Desktop\3DSPI\pmp\3\右边界2\140_20_0.bmp');
Ref2 = imread('E:\Desktop\3DSPI\pmp\3\右边界2\140_20_1.bmp');
Ref3 = imread('E:\Desktop\3DSPI\pmp\3\右边界2\140_20_2.bmp');
Ref4 = imread('E:\Desktop\3DSPI\pmp\3\右边界2\140_20_3.bmp');
Ref5 = imread('E:\Desktop\3DSPI\pmp\3\右边界2\140_20_4.bmp');

% Ref1 = imread('E:\Desktop\3DSPI\metch\3\120_10.bmp');
% Ref2 = imread('E:\Desktop\3DSPI\metch\3\120_11.bmp');
% Ref3 = imread('E:\Desktop\3DSPI\metch\3\120_12.bmp');
% Ref4 = imread('E:\Desktop\3DSPI\metch\3\120_13.bmp');
% Ref5 = imread('E:\Desktop\3DSPI\metch\3\120_14.bmp');

 Test1 = Ref4();
 
 Ref1 = im2double(Ref1);
 Ref2 = im2double(Ref2);
 Ref3 = im2double(Ref3);
 Ref4 = im2double(Ref4);
 Ref5 = im2double(Ref5);

 
%  temp1 = Ref4-Ref2;
%  temp2 = Ref1-Ref3;
 
temp1 = 2*Ref2-2*Ref4;
temp2 = 2*Ref3-Ref1-Ref5;

thetaREF1 = atan2(temp1,temp2); %M自带计算反三角函数

b = thetaREF1;  

lamda=10;%调节参数
[M,N]=size(b);
dx=zeros(M,N);
dy=zeros(M,N);
 

%求包裹相位的相位差
for n=1:N
  for m=1:M-1
         if((b(m+1,n)-b(m,n))>pi)
           dx(m,n)=b(m+1,n)-b(m,n)-2*pi;
         elseif((b(m+1,n)-b(m,n))<-pi)
           dx(m,n)=b(m+1,n)-b(m,n)+2*pi;
         else
           dx(m,n)=b(m+1,n)-b(m,n);
         end
      end
 end
 for n=1:N
     dx(M,n)=0;
 end
 
 for m=1:M
   for n=1:N-1
         if((b(m,n+1)-b(m,n))>pi)
         dy(m,n)=b(m,n+1)-b(m,n)-2*pi;
         elseif((b(m,n+1)-b(m,n))<-pi)
         dy(m,n)=b(m,n+1)-b(m,n)+2*pi;
         else
          dy(m,n)=b(m,n+1)-b(m,n);
         end
     end
 end
 for m=1:M
     dy(m,N)=0;
 end
 
p2=zeros(M,N);
%求相位差
 for m=2:M
     for n=2:N
         p2(m,n)=dx(m,n)-dx(m-1,n)+dy(m,n)-dy(m,n-1);
     end
 end
 p2(1,1)=dx(1,1)+dy(1,1);
 for n=2:N
     p2(1,n)=dx(1,n)+dy(1,n)-dy(1,n-1);
 end
 for m=2:M
     p2(m,1)=dx(m,1)-dx(m-1,1)+dy(m,1);
 end
 

% %DCT变换
 pp=dct2(p2);        %M自带离散余弦变换函数
 
 fi=zeros(M,N);
 for m=1:M
     for n=1:N  
       fi(m,n)=pp(m,n)./((2+16*lamda)*(cos(pi*(m-1)/M)+cos(pi*(n-1)/N))-2*lamda*(cos(2*pi*(m-1)/M)+cos(2*pi*(n-1)/N))-8*lamda*cos(pi*(m-1)/M)*cos(pi*(n-1)/N)-20*lamda-4);
       %fi(m,n) = pp(m,n) ./ (2*(cos(pi*(m-1)/M)+cos(pi*(n-1)/N)-2));
    end
 end
 fi(1,1)=pp(1,1);
 
 y01=idct2(fi);
 
 

y01 = y01();

M = M-20;
N = N-20;
y01 = y01(1:M,1:N);
mesh(y01());title('解包相位');
[M,N]=size(y01);

后续进行曲面拟合的思路如下：在锡膏的周围往外扩散20个像素点的区域，沿x或y方向找取该方向上最直的那条线段或多条线段，判定为可选取拟合点的线，最后找取了很多点进行拟合，得到基准面，再与解相位的面进行相减得到相位差，最后根据标定的系数进行高度的计算。

关键步骤代码：

找取最直的直线，根据二次梯度进行判断：

ys = y01(i,:);
Test1s = Test1(i,:);
%所有点的间距是固定的，那么可以通过y值的差异来标识斜率的不同，当然若间距不固定，就用正常的斜率求法
diffYs = diff(ys);              %M自带函数：数组中（i+1）-i
diff2Ys = diff(diffYs);

根据阈值判断可取的直线起始点 

    k = 0.01;
    diff2Ys( find(abs(diff2Ys < k)) ) = 0;

    zeroIndexs = find(abs(diff2Ys < k));
    flags = ones(1,length(diff2Ys));
    flags( zeroIndexs ) = 0;
    
    %flags里的0值存储的就是斜率相近的点，即二阶导数接近0的点。其他置1.求出每个0值序列，选取要求的那个序列
    diff_flags = diff( flags );
    
    start_0 = find( diff_flags == -1 ) + 1; %1-0为起点，第一二个数为零是start_0会少第一个数1
    end_0 = find( diff_flags == 1 );        %此处会遗漏当最后一个点diff为零时的end点

    %line_num中最大的数的位置代表似直线最长的一段
    line_num = end_0 - start_0;
    
    %line_num最大数坐标映射到start_0和end_0代表曲线的起始和终点
    [m,index] = max(line_num);         %index为位置
    start_line = start_0(index);
    end_line = end_0(index);

%选取拟合点 
   if end_0(index)-start_0(index) > 15
       get_y = randperm(end_0(index)-start_0(index),5) + start_0(index);
       for get_num = 1:5
            standard_pro = [standard_pro;[i get_y(get_num)  ys(get_y(get_num))]];
       end
       
   else
       continue;
   end

 值得一提的是对于阴影区域的高度计算障碍排除：         ，像这样的。

不知道是解相位算法的问题，还是通用的问题：在此处解出来的相位全是混乱波动的。

方法：

    %找到阴影处边界，将梯度加速度置1
    Test1s(find(Test1s<2)) = 0;
    
    shadow_nf = find(Test1s ~=0,1,'first');
    shadow_n = find(Test1s~=0,1,'last');
    diff2Ys(shadow_n + 1: N) = 1;
    diff2Ys(1: shadow_nf - 1) = 1;

最后进行拟合：

%拟合曲面
pro_x=standard_pro(:,1);
pro_y=standard_pro(:,2);
pro_z=standard_pro(:,3);

[X,Y,Z]=griddata(pro_y,pro_x,pro_z,linspace(1,N,N)',linspace(1,M,M),'v4');  %M自带拟合曲面函数

figure;mesh(X,Y,Z);title('profit');


phase_diff = y01(1:M,1:N) - Z;
height = 315* phase_diff;

最后一些结果展示：

存在问题：在阴影与锡膏的交界处，总会产生翘曲或者凹陷导致高度的突增或消失，暂时还未能解决。。。

后续