matlab实现单峰物体复原--光栅投影-多频外差

matlab仿真实现单峰物体复原–光栅投影-多频外差

多频外差

关于多频外差上一章已经讲到，其中关于处理步骤不太详细，这里会结合代码来进行讲解.代码是在理想环境下模拟调制后的图像，并用多频外差进行求解主项位值。

%  目标  条栅三维重建并用三频四相法解包裹相位
% 参数设置
% 画布大小  
W = 1024;
H = 1140;


% 三个不同的频率  保证最后计算出来的等效周期大于整幅图像
f = [70 64 59];
%T = [1/70 1/64 1/59];
T = [H/70 H/64 H/59];
A = 1; % 背景光强
B = 1; % 调制光强
% 模型参数
d = 20; % 投影仪与相机的距离
l = 100; % 相机与参考面之间的距离


% 生成条纹图像 
Ce1 = cell(3,4);   % 用来存储调制前的条纹图像
for i = 1:3
    for j = 1:4
            Ce1{i,j} = zeros(W,H);
    end
end
for i=1:3
    pha = linspace(0,f(i)*2*pi,H);
    pha = repmat(pha' ,1, W);
    for j = 0:3
    I0 = A + B * cos( pha' + j*pi/2 ); %光强
    Ce1{i,j+1}=I0;%记录光强
    end
end
% 归一化处理
for i = 1:3
    for j = 1:4
        Ce1{i,j} = mat2gray(Ce1{i,j});
    end
end
% 打印条纹图像
% for i = 1:3
%     for j = 1:4
%          tmp=Ce1{i,j};
%          if i==1
%            figure(j); %弹出窗口
%            imshow(Ce1{i,j},[]); %图像显示
%          elseif i>1
%            figure(2.^(i)+j); % 弹出窗口
%            imshow(Ce1{i,j},[]);  %图像显示
%            
%          end      
%     end
% end

% 生成物体
h = halfball(65,512,570); % m  球半径150m心坐标（256,256）仿真假设一个像素为1mm
figure(2);mesh(h);colormap(jet);set(gcf,'Name','原半球物体','NumberTitle','off');shading interp;
% 计算调制后的图像
% 因为是仿真，所以这里调制后的条纹图像表达式可以根据第一章的高度相位公式进行反推 将高度当成  %  已知量
Ce2 = cell(3,4);
for i = 1:3
    for j = 1:4
        Ce2{i,j} = zeros(W,H);
    end
end
de = zeros(W,H);
for i=1:3
    pha = linspace(0,f(i)*2*pi,H);
    pha = repmat(pha' ,1, W);    
    de = 2*pi*d*h ./((l-h)*T(i));
    for j = 0:3
    Ir = A + B * cos( pha'+j*pi/2  + de ); %光强
    Ce2{i,j+1}=Ir;
    end
end

% 打印条纹图像
% for i = 1:3
%     for j = 1:4
%          tmp=Ce2{i,j};
%          if i==1
%            figure(j); %弹出窗口
%            imshow(Ce2{i,j},[]); %图像显示
%          elseif i>1
%            figure(2.^(i)+j); % 弹出窗口
%            imshow(Ce2{i,j},[]);  %图像显示
%            
%          end      
%     end
% end

% 求解参考面的包裹相位
pha_C1 = atan2(Ce1{1,4}-Ce1{1,2}, Ce1{1,1}- Ce1{1,3});
pha_C2 = atan2(Ce1{2,4}-Ce1{2,2}, Ce1{2,1}- Ce1{2,3});
pha_C3 = atan2(Ce1{3,4}-Ce1{3,2}, Ce1{3,1}- Ce1{3,3});
% 求解物品面的包裹相位
pha_W1 = atan2(Ce2{1,4}-Ce2{1,2}, Ce2{1,1}- Ce2{1,3});
pha_W2 = atan2(Ce2{2,4}-Ce2{2,2}, Ce2{2,1}- Ce2{2,3});
pha_W3 = atan2(Ce2{3,4}-Ce2{3,2}, Ce2{3,1}- Ce2{3,3});
% figure(30);
% mesh(pha_W2);

% 求解等效周期
T12 = (T(1) * T(2)) / (T(2) - T(1));
T23 = (T(2) * T(3)) / (T(3) - T(2));
T123 = (T12 * T23) / (T23 - T12);

% 将参考面的截断相位映射到0-2pi
pha_C1(pha_C1 < 0) = pha_C1(pha_C1 < 0) + 2*pi ;
pha_C2(pha_C2 < 0) = pha_C2(pha_C2 < 0) + 2*pi ;
pha_C3(pha_C3 < 0) = pha_C3(pha_C3 < 0) + 2*pi ;

pha_W1(pha_W1 < 0) = pha_W1(pha_W1 < 0) + 2*pi ;
pha_W2(pha_W2 < 0) = pha_W2(pha_W2 < 0) + 2*pi ;
pha_W3(pha_W3 < 0) = pha_W3(pha_W3 < 0) + 2*pi ;


% 计算参考面外差   
pha_C12 = pha_C1 - pha_C2;
pha_C12(pha_C12 < 0) = pha_C12(pha_C12 < 0) + 2*pi;
pha_C23 = pha_C2 - pha_C3;
pha_C23(pha_C23 < 0) = pha_C23(pha_C23 < 0) + 2*pi;
pha_C123 = pha_C12 - pha_C23;
pha_C123(pha_C123 < 0) = pha_C123(pha_C123 < 0) + 2*pi;




% 计算物品面
pha_W12 = pha_W1 - pha_W2;
pha_W12(pha_W12 < 0) = pha_W12(pha_W12 < 0) + 2*pi;
pha_W23 = pha_W2 - pha_W3;
pha_W23(pha_W23 < 0) = pha_W23(pha_W23 < 0) + 2*pi;
pha_W123 = pha_W12 - pha_W23;
pha_W123(pha_W123 < 0) = pha_W123(pha_W123 < 0) + 2*pi;


% 求解参考面的绝对相位选择频率最高的那个
% 这里要注意求解顺序，要先用求的得pha_W123这个绝对相位对pha_W23或者pha_W12进行展开，
 %再利用已经展开的绝对相位展开pha_W1，pha_W2，pha_W3  对参考平面的展开也是如此
K_C23 = round((pha_C123*T123/T23 - pha_C23)/ (2*pi) );  
pha_TC23 = pha_C23 + 2*pi*K_C23;
% figure(18);
% mesh(pha_TC23);

K_C1= round((pha_TC23*T23/T(1) - pha_C1)/ (2*pi));
pha_TC1 = pha_C1 + 2*pi*K_C1;
figure(18);
mesh(pha_TC1);

K_W23 = round(((pha_W123*T123)/T23 - pha_W23)/ (2*pi));
pha_TW23 = pha_W23 + 2*pi*K_W23;
% figure(19);
% mesh(pha_TW23);

K_W1 = round(((pha_TW23*T23)/T(1)- pha_W1)/ (2*pi));
pha_TW1 = pha_W1 + 2*pi*K_W1;
figure(19);
mesh(pha_TW1);

% 求解实际物体高度
ph =   pha_TW1-pha_TC1;
Th = (ph*T(1)*l)./(ph*T(1)+(2*pi)*d);
figure(20);
mesh(Th);

bb = h - Th;
sum = 0;
for i = 1:W
    for j = 1:H
        sum = sum + bb(i,j);
    end
end
sum;



function output = halfball(r,x,y)  % 半球绘制函数
        temp = zeros(1024,1140);
        for m = 1:1024
            for n = 1:1140
                if((m-x)^2 + (n-y)^2 < r^2)
                    temp(m,n) = sqrt(r^2 - (x-m)^2 - (y-n)^2);
                end
            end
        end
        output = temp;
    end

代码和注释应该比较详细，相信对照原理可以顺着读下来。
学习注定是个漫长的过程，如有错误欢迎指正，留下博客更多的是对自己的学习过程进行梳理和总结。禁止转载哦！