[转]多频外差法(三频四相)理论及代码

在提取相位时,如果图像内的周期数仅为一个,则相对相位就是绝对相位,则求得的相位主值就是它的相位值。

实际中,我们投射的光栅并不是一个周期,那么在整个测量的空间中就会求得多个相同的相位主值,这时就需要对包裹的相位值进行展开。

相位包裹又是什么意思呢?

相位包裹的主要原因是相移法解相时使用了反正切函数,用atan2函数,得到四个象限的反正切,所以计算的相位都是在(-pi,pi]之间,也就是被包裹。

图片

因而真实的相位值和被包裹的相位主值之间是相差 2npi的,其中 n 为整数,即为

图片

多频外差,字面上意思就是多个频率,利用外差原理。

论文里看到的解释是这样的:

(1)多频外差编码法利用拍频原理,将两个或两个以上周期相近但不相同的相移编码光组合成一组编码光,把解得的不同周期的多个相位作差,将包裹相位的小周期放大为相位差的大周期,直到相位差信号的周期包含整个测量视场。

图片

图片

这种方法中是用周期表示,也就是条纹的宽度,而且(T2>T1)。

(2) 李中伟博士的博士论文里是这样的:

图片

图片

图片

如果图片对应的是频率,这里的值应该是小于1的值,即正弦波的个数除以整个图像方向上的像素值,它的周期才会是条纹的宽度,但在图中图片所对应的并非这种含义,所以我认为这里解释和图示是有矛盾的。

正确图示可能是这样的:

图片

从双频外差可以看出是将小的周期扩展到了一个大的周期内,如果频率差为1可以达到全场范围内只有一个周期的相位的目的,展开仅与频率相位有关,不受颜色的影响。

代码如下:

% 程序开始
clc;
close all;
clear;
 
% 图片的初始化
width = 1024;  
heigth = 768;
 
% 三频率
% 这个可以参见李中伟的博士论文
freq = [70 64 59];  %像素单位为个数,可以看做频率;正弦函数为周期含义
 
% 利用分块矩阵C存储3组共计12张图
% 三种频率,四组相位
C = cell(3,4);  
for i=1:3
    for j=1:4
        C{i,j} = zeros(heigth,width);
    end
end
 
% 利用余弦函数计算12张图的灰度值
% 图像的生成
% 三种频率,四组相位
for i = 1:3 % 对应三种不同的频率
    for  j = 0:3 % 对应四种相位
        for k = 1:width
            C{i,j+1}(:,k) = 128+127*sin(2*pi*k*freq(i)/width+j*pi/2);
        end
    end
end
 
% 对灰度值进行归一化处理
for i = 1:3
    for j = 1:4
        C{i,j} = mat2gray(C{i,j});
    end
end
 
% 显示12张图
% for i = 1:3
%     for j = 1:4
%         n = 4*(i-1)+j;
%         h = figure(n);
%         imshow(C{i,j});
%     end
% end
 
% 初始化三组处理后的图片灰度矩阵
% phi也是分块矩阵
% 存储相位主值图像
phi = cell(3,1);
for i = 1:3
    phi{i,1} = zeros(heigth,width);
end
 
% 求取相位差
% 计算每种频率对应的相位主值
% 输出三种频率的相位主值,用于相差计算
for i = 1:3 % 对于3组中的每一组图片,每一组相同频率的有四张图片
     I1 = C{i,1};
     I2 = C{i,2};
     I3 = C{i,3};
     I4 = C{i,4};
      for g = 1:heigth
          for k = 1:width          
            if     I4(g,k)==I2(g,k)&&I1(g,k)>I3(g,k) %四个特殊位置
                    phi{i,1}(g,k)=0;
            elseif I4(g,k)==I2(g,k)&&I1(g,k)<I3(g,k) %四个特殊位置
                    phi{i,1}(g,k)=pi; 
            elseif I1(g,k)==I3(g,k)&&I4(g,k)>I2(g,k) %四个特殊位置
                    phi{i,1}(g,k)=pi/2;
            elseif I1(g,k)==I3(g,k)&&I4(g,k)<I2(g,k) %四个特殊位置
                    phi{i,1}(g,k)=3*pi/2;
            elseif I1(g,k)<I3(g,k) %二三象限
                    phi{i,1}(g,k)=atan((I4(g,k)-I2(g,k))./(I1(g,k)-I3(g,k)))+pi;
            elseif I1(g,k)>I3(g,k)&&I4(g,k)>I2(g,k) %第一象限
                    phi{i,1}(g,k)=atan((I4(g,k)-I2(g,k))./(I1(g,k)-I3(g,k)));
            elseif I1(g,k)>I3(g,k)&&I4(g,k)<I2(g,k) %第四象限
                    phi{i,1}(g,k)=atan((I4(g,k)-I2(g,k))./(I1(g,k)-I3(g,k)))+2*pi;  
            end
          end             
     end
end
  
% 计算相差
% 保存矩阵,用于多频相差的计算
PH1 = phi{1,1};   %频率1
PH2 = phi{2,1};   %频率2
PH3 = phi{3,1};   %频率3
 
% 初始化相差变量
% 多频相差
PH12 = zeros(heigth,width);
PH23 = zeros(heigth,width);
PH123 = zeros(heigth,width);
 
% 计算相差
% 相差计算
% 解相
for g = 1:heigth
  for k = 1:width
    % 计算第一组和第二组的相差
    if PH1(g,k)>PH2(g,k)
      PH12(g,k) = PH1(g,k)-PH2(g,k);
    else
      PH12(g,k) = PH1(g,k)+2*pi-PH2(g,k);
    end
    % 计算第二组和第三组的相差
    if PH2(g,k)>PH3(g,k)
      PH23(g,k) = PH2(g,k)-PH3(g,k);
    else
      PH23(g,k) = PH2(g,k)+2*pi-PH3(g,k);
    end
     % plot(1,k);
  end
end
 
% 计算最终相差
% 相差图案
% 相位解包裹  相位展开
for g = 1:heigth
  for k = 1:width 
    if PH12(g,k)>PH23(g,k)
      PH123(g,k) = PH12(g,k)-PH23(g,k);
    else
      PH123(g,k) = PH12(g,k)+2*pi-PH23(g,k);
    end
  end
end

% 显示
figure,imshow(mat2gray(PH12));title('1,2外差');   imwrite(mat2gray(PH12),'12外差.bmp');
figure,imshow(mat2gray(PH23));title('2,3外差');   imwrite(mat2gray(PH23),'23外差.bmp');
figure,imshow(mat2gray(PH123));title('1,2,3外差');imwrite(mat2gray(PH123),'123外差.bmp');
figure,imshow(mat2gray(PH1));title('1相位主值');  imwrite(mat2gray(PH1),'1相位主值.bmp');  
figure,imshow(mat2gray(PH2));title('2相位主值');  imwrite(mat2gray(PH2),'2相位主值.bmp');  
figure,imshow(mat2gray(PH3));title('3相位主值');  imwrite(mat2gray(PH3),'3相位主值.bmp');

结果如下:

相位主值1

图片

相位主值2

图片

相位主值3

图片

12外差计算

图片

23外差计算

图片

123外差计算

图片

文章来源:
https://blog.csdn.net/qq_27353621/article/details/124235712


---------------------
作者:马少爷
来源:CSDN
原文:https://blog.csdn.net/qq_27353621/article/details/124235712
版权声明:本文为作者原创文章,转载请附上博文链接!
内容解析By:CSDN,CNBLOG博客文章一键转载插件

  • 1
    点赞
  • 8
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值