三维形貌测量----光学位相轮廓术（PMP）

光学位相测量轮廓术(PMP)是一种基于光栅投影和相位测量的光学三维面形测量技术, 可有效实现物

体表面三维轮廓的自动测量。参考文献：点击打开链接

PMP三维型貌测量MATLAB仿真：

% ************************************************************************
%           PMP[Phase measuring profilometry]仿真程序
% ************************************************************************
% ************************************************************************
%                           初始化
% ************************************************************************
clc
clear all
close all
% ************************************************************************
%                           构造待测物
% ************************************************************************
Row=300;
Column=300;
x=1:Row;%参考面横向像素范围
y=1:Column;%参考面纵向像素范围
[X,Y]=meshgrid(x,y);
 z=real(0.4*sqrt(100^2-(X-150).^2-(Y-150).^2))*0.01;
% M =300;
% z = peaks(M);

figure(1)
mesh(X,Y,z)
title('Profile to be measured','FontSize',20)
xlabel('x [pixel]','FontSize',20)
ylabel('y [pixel]','FontSize',20)
axis tight
zlabel('Hight [pixel]','FontSize',20)
set(gca,'FontSize',20,'linewidth',2.5,'color','w','FontWeight','bold','FontName','Times New Roman');
set(0,'defaultfigurecolor','w');
% ************************************************************************
%                           参考正弦光栅光源
% ************************************************************************
A=1;%光场强度系数,实验中为未知待求量
B=1;%光场强度系数,实验中为未知待求量
f=1/60;%正弦光栅周期设定为6个像素
phi=0;%光场本征相位分布值，实验中为未知待求量
Ir=A+B*cos(2*pi*f*X+phi);
figure(2)
imshow(Ir)
title('Reference grating intensity distribution','FontSize',20)
xlabel('x [pixel]','FontSize',20)
ylabel('y [pixel]','FontSize',20)
set(gca,'FontSize',20,'linewidth',2.5,'color','w','FontWeight','bold','FontName','Times New Roman');
set(0,'defaultfigurecolor','w');
% ************************************************************************
%                           变形正弦光栅光强分布
% ************************************************************************
d=20;%光栅光源光心与CCD光心连线距离，单位pixel
l=100;%CCD光心到参考平面的距离，单位pixel
I0=A+B*cos(2*pi*f*X+phi+d*z/l);
figure(3)
imshow(I0)
title('Deformed grating intensity distribution','FontSize',20)
xlabel('x [pixel]','FontSize',20)
ylabel('y [pixel]','FontSize',20)
set(gca,'FontSize',20,'linewidth',2.5,'color','w','FontWeight','bold','FontName','Times New Roman');
set(0,'defaultfigurecolor','w');
figure;mesh(I0);
% ************************************************************************
%    为求出变形光栅图样中各点的截断相位值，对变形图样进行N（N>=3）步相移操作
%                  并根据相移图像计算截断相位分布矩阵
% ************************************************************************
N=3;%以三步相移法为例
ITcos=zeros(length(x),length(y));
ITsin=zeros(length(x),length(y));
for j=1:N
   eval(['I0',num2str(j),'=A+B*cos(2*pi*f*X+phi+d*z/l+2*pi*(j-1)/N);'])
   figure(3+j)
   eval(['imshow(I0',num2str(j),')'])
   title(['Deformed grating with ',num2str(j),' phase shift'],'FontSize',20)
   xlabel('x [pixel]','FontSize',20)
   ylabel('y [pixel]','FontSize',20)
   set(gca,'FontSize',20,'linewidth',2.5,'color','w','FontWeight','bold','FontName','Times New Roman');
   set(0,'defaultfigurecolor','w');
   eval(['ITcos=ITcos+I0',num2str(j),'*cos(2*pi*(j-1)/N);'])
   eval(['ITsin=ITsin+I0',num2str(j),'*sin(2*pi*(j-1)/N);'])
end
figure(3+N+1)
imshow(ITcos)
title([num2str(N),'Step truncated phase diagram'],'FontSize',20)
xlabel('x [pixel]','FontSize',20)
ylabel('y [pixel]','FontSize',20)
set(gca,'FontSize',20,'linewidth',2.5,'color','w','FontWeight','bold','FontName','Times New Roman');
set(0,'defaultfigurecolor','w');
figure(3+N+2)
imshow(ITsin)
title([num2str(N),'Step truncated phase diagram'],'FontSize',20)
xlabel('x [pixel]','FontSize',20)
ylabel('y [pixel]','FontSize',20)
set(gca,'FontSize',20,'linewidth',2.5,'color','w','FontWeight','bold','FontName','Times New Roman');
set(0,'defaultfigurecolor','w');
% ************************************************************************
%             将截断相位分布用相位展开法恢复成原有的连续相位分布
% ************************************************************************
Phi=-atan2(ITsin,ITcos);
figure(3+N+3)
imshow(Phi)
title('Wrapped Phase','FontSize',20)
xlabel('x [pixel]','FontSize',20)
ylabel('y [pixel]','FontSize',20)
set(gca,'FontSize',20,'linewidth',2.5,'color','w','FontWeight','bold','FontName','Times New Roman');
set(0,'defaultfigurecolor','w');

h=Phi*l/d;
H=h-ones(Row,1)*h(1,:);
figure;
mesh(X,Y,H)
%***********方法1**************
Phi_unwrapped=Phi;
for i=1:Row
Phi_unwrapped(i,:) = unwrap(Phi_unwrapped(i,:));
end
for i=1:Column
Phi_unwrapped(:,i) = unwrap(Phi_unwrapped(:,i));
end

figure(3+N+4)
imagesc(Phi_unwrapped)
title('Unwrapped Phase','FontSize',20)
xlabel('x [pixel]','FontSize',20)
ylabel('y [pixel]','FontSize',20)
set(gca,'FontSize',20,'linewidth',2.5,'color','w','FontWeight','bold','FontName','Times New Roman');
set(0,'defaultfigurecolor','w');
% ************************************************************************
%                                 三维轮廓重构
% ************************************************************************
h=Phi_unwrapped*l/d;
H=h-ones(Row,1)*h(1,:);
figure(3+N+5)
mesh(X,Y,H)
title('Rebuilding','FontSize',20)
xlabel('x [pixels]','FontSize',20)
ylabel('y [pixels]','FontSize',20)
zlabel('Hight [pixels]','FontSize',20)
axis tight
set(gca,'FontSize',20,'linewidth',2.5,'color','w','FontWeight','bold','FontName','Times New Roman');
set(0,'defaultfigurecolor','w');