%?測試多頻外差原理
%?模擬程序
%?圖像尺寸?1024×768
%?共4張/組×3組=12張圖
%?3組頻率依次為1/701/641/59
%?本程序基于李中偉的博士論文來實現
%?《基于數字光柵投影的結構光三維測量技術與系統的研究》

%?程序開始
clc;
close?all;
clear;

%?圖片的初始化
width?=?1024;??
heigth?=?768;

%?三頻率
%?這個可以參見李中偉的博士論文
freq?=?[70?64?32?16];

%?利用分塊矩陣C存儲3組共計12張圖
%?三種頻率，四組相位
C?=?cell（44）;??
for?i=1:4
????for?j=1:4
????????C{ij}?=?zeros（heigthwidth）;
????end
end

%?利用余弦函數計算12張圖的灰度值
%?圖像的生成
%?三種頻率，四組相位
for?i?=?1:4?%?對應三種不同的頻率
????for??j?=?0:3?%?對應四種相位
????????for?k?=?1:width
????????????C{ij+1}（:k）?=?128+127*sin（2*pi*k*freq（i）/width+j*pi/2）;
????????end
????end
end

%?對灰度值進行歸一化處理
for?i?=?1:4
????for?j?=?1:4
????????C{ij}?=?mat2gray（C{ij}）;
????end
end

%?顯示12張圖
for?i?=?1:4
????for?j?=?1:4
?????????n?=?4*（i-1）+j;
?????????h?=?figure（n）;
?????????imshow（C{ij}）;
?????end
?end

%?初始化三組處理后的圖片灰度矩陣
%?phi也是分塊矩陣
%?存儲相位主值圖像
phi?=?cell（41）;
for?i?=?1:4
????phi{i1}?=?zeros（heigthwidth）;
end

%?求取相位差
%?計算每種頻率對應的相位主值
%?輸出三種頻率的相位主值，用于相差計算
for?i?=?1:4?%?對于3組中的每一組圖片
?????I1?=?C{i1};
?????I2?=?C{i2};
?????I3?=?C{i3};
?????I4?=?C{i4};
??????for?g?=?1:heigth
??????????for?k?=?1:width??????????
????????????if?I4（gk）==I2（gk）&&I1（gk）>I3（gk）?%四個特殊位置
????????????????????phi{i1}（gk）=0;
????????????elseif?I4（gk）==I2（gk）&&I1（gk）????????????????????phi{i1}（gk）=pi;?
????????????elseif?I1（gk）==I3（gk）&&I4（gk）>I2（gk）?%四個特殊位置
????????????????????phi{i1}（gk）=pi/2;
????????????elseif?I1（gk）==I3（gk）&&I4（gk）????????????????????phi{i1}（gk）=3*pi/2;
????????????elseif?I1（gk）????????????????????phi{i1}（gk）=atan（（I4（gk）-I2（gk））./（I1（gk）-I3（gk）））+pi;
????????????elseif?I1（gk）>I3（gk）&&I4（gk）>I2（gk）?%第一象限
????????????????????phi{i1}（gk）=atan（（I4（gk）-I2（gk））./（I1（gk）-I3（gk）））;
????????????elseif?I1（gk）>I3（gk）&&I4（gk）????????????????????phi{i1}（gk）=atan（（I4（gk）-I2（gk））./（I1（gk）-I3（gk）））+2*pi;??
????????????end
??????????end?????????????
?????end
end
??
%?計算相差