一、基本原理：

    正弦條紋打在三維物體上，CCD記錄到的條紋由於受到三維物體高度的調製而發生扭曲，扭曲的條紋（deformed fringe）實質上為原始條紋在物體具有高度存在的位置有了附加相位，各點的相位表現為由CCD影象採集獲得的被調製的條紋數字影象的灰度值。通過扭曲的條紋和原始條紋比對計算得出相位變化值。又已知投影儀、CCD和物體的具體位置和之間的距離，利用數學關係可求出對應點的高度值，實現3D重建。公式推導如下：（下圖來自[1]）

    獲得正確的相位值並計算出相位差是得到物體真實高度資訊並進行3D重建的關鍵。

    在利用CCD進行數字採集過程中，由於環境等各種因素的影響，以及物體在空間上對光的反射程度和與CCD距離不同，不能夠得到理想的均勻而又高對比度的條紋，也就不能得到最真實的相位資訊。運用四步相移法處理此問題（此方法參考[1]）,過程分析如下：

式中 R(x,y) 表示所測物體表面的不均勻反射率，A(x,y) 是背景強度，B(x,y)/A(x,y) 表示的是光柵條紋的對比度。∅(x,y) 是相位值。

    考慮到以上引入的多個外部影響因子，利用四步相移法分別記錄四次初相不同的正弦條紋和被調製的條紋（每個相差π/2），聯立即可得到消去外部因子的真實相位值。表示式如下：

分別為初相等於0、π/2、π 和 3π/2 的條紋表示式，由三角函式變換可分別化成上式。

聯立解得  

    可以看到A、B、R三個外部因子被消去，留下相位的正切值。利用反正切函式即可求出相位分佈。

三、包裹相位處理（相位展開/相位解卷）

1、 相位包裹的含義：

    此處有另一個重點，數學中常用的反正切函式是二象限的反正切函式，值域僅是(-π/2,π/2)，而此處我們運用的反正切函式是四象限的反正切函式，值域是(-π,π)。區別如下：
二象限反正切函式（matlab中函式名為atan）輸入引數為正切值，返回正切值僅有兩個第一第四兩個象限的角度值，即(-π/2,π/2)。而四象限反正切函式（matlab中函式名為atan2）輸入引數為y和x，返回正切值包括四個象限，即(-π,π)。
    四象限反正切考慮到正切值的由來是兩個量的商，多考慮了x為負時的另一種情況，即可以將值域擴充套件一倍，延長到(-π,π)。表示式如下：（表示式來自Wikipedia）

    上圖為本人用matlab粗略匯出的四象限反正切函式曲線，各條曲線顏色對應下圖圓中角所在的區域。

2、解卷/相位展開/相位解包裹操作（phase unwrapping）

    上面講述的相位包裹表現為相位在主值區間內增加到超過主值區間的上限的時候，即超過π時，將發生跳變，也就是說繼續增加相位值會跳變到-π。在整體提取的影象中沿著相位增加的方向表現為鋸齒波的形狀。如下圖為沒有被物體調製過的條紋沿著條紋的方向橫截，被包裹的相點陣圖：

    相位展開最根本的思想就是找出相位發生跳變的地方，補上被減去的n個2π，使跳變出相位值和周圍的值變得連續、平滑。方法很多，請參考[1]，文中詳細講解質量法。
    本次使用matlab模擬直接使用了matlab中相位解卷函式unwrap(P,tol,1)。其中P為執行相位解卷的向量或矩陣；tol為相位跳變限度（jump tolerance），即差值大於tol被判定為發生跳變，預設為pi，此處我使用的時候填pi；第三個引數表示進行相位解卷操作的維度，向量只有一個維度，矩陣有兩個維度，則1表示行解卷，2表示列解卷。上一張圖相位展開完成後相位分佈如下圖，可以看到不存在跳變（鋸齒）。

四、傅立葉變換法還原相位

    運用此方法不進行相位解卷便可以獲得真實的相位資訊。原理如下：（此方法和原理公式來自[2]）
被包裹的相位表示為

其中φ(x,y)為真實的相位資訊，fx和fx標誌條紋頻率。
將被包裹的相位表示成複數形式

二維傅立葉變換

頻移後進行逆傅立葉變換

求複數表示的相位的相角即可得到真實相位值

兩次傅立葉變換亦可以應用傅立葉變換的頻移性質實現：

    此方法就是利用傅立葉變換頻移實現濾波操作，濾掉條紋部分的資訊。[2]對該方法提高精度的一些辦法進行了闡述，具體參考原文。
    此方法雖然可以直接得到相位資訊而不需要進行解卷，且精度較高，但是因為沒有解卷，若物體高度換算大於2π的相位差，則或出現重建變形，較高部分凹陷下去的情況。

五、matlab模擬結果

    模擬基本可以重建出半球的形狀。半球最高點高度理論上應該是半徑的長度，由於計算所用的正弦波週期是通過離散傅立葉變換頻譜峰值提取出的頻率，提取出的頻率存在誤差。

參考文獻：

[1] 潘玉玲. 相移法三維測量系統[D]. 重慶大學, 2013.

[2] Du G, Wang M, Zhou C, et al. Improved method for phasewraps reduction in profilometry[J]. 2016.

六、模擬matlab程式碼

function my_3D_reconstruction
u = pi/25;v = pi/25;  % 給定初始條紋角頻率
fringe1 = sin_wave(u,v,0);
fringe2 = sin_wave(u,v,pi/2);
fringe3 = sin_wave(u,v,pi);
fringe4 = sin_wave(u,v,pi*3/2);
figure(1);set(gcf,'Name','四步相移原始正弦條紋','NumberTitle','off');
subplot(221);imshow(fringe1);title('相移：0');
subplot(222);imshow(fringe2);title('相移：pi/2');
subplot(223);imshow(fringe3);title('相移：pi');
subplot(224);imshow(fringe4);title('相移：pi*3/2');

h = halfball(65,256,256); % 半徑150mm  球心座標（256,256）模擬假設一個畫素為1mm
figure(2);surf(h);colormap(jet);set(gcf,'Name','原半球物體','NumberTitle','off');shading interp;

D = 155; L = 660; T = 50/sqrt(2); % 引數單位 mm
delta_phi = 2*pi*D*h./(L-h)/T; % T為正弦條紋的空間週期
deformed_fringe1 = fringe_modulation(u,v,0,delta_phi); % 生成四步相移的被調製條紋
deformed_fringe2 = fringe_modulation(u,v,pi/2,delta_phi);
deformed_fringe3 = fringe_modulation(u,v,pi,delta_phi);
deformed_fringe4 = fringe_modulation(u,v,pi*3/2,delta_phi);
figure(4);set(gcf,'Name','四步相移被物體調製的條紋','NumberTitle','off');
subplot(221);imshow(deformed_fringe1);title('相移：0');
subplot(222);imshow(deformed_fringe2);title('相移：pi/2');
subplot(223);imshow(deformed_fringe3);title('相移：pi');
subplot(224);imshow(deformed_fringe4);title('相移：pi*3/2');

F = fft2(fringe1-0.5); % 傅立葉變換檢測峰值提取正弦條紋的空間頻率fx、fy
figure(3);mesh(abs(F));set(gcf,'Name','正弦波頻譜','NumberTitle','off');
[row,col] = find(abs(F) == max(max(abs(F))));
fx = (row(1)-1)/512;  
fy = (col(1)-1)/512;
 Tx = fx^(-1);
 Ty = fy^(-1);    % 提取出正弦條紋的頻率,換算成空間週期
 T1 = Tx*Ty/sqrt(Tx^2+Ty^2);

wrapped_phase = atan2((deformed_fringe4-deformed_fringe2),(deformed_fringe1-deformed_fringe3));
original_phase = atan2((fringe4-fringe2),(fringe1-fringe3));  
%相位包裹的被調製的條紋和原始條紋計算

%----------傅立葉變換頻移法還原相位------------------------------------------
k = exp(1i*wrapped_phase);
[x,y] = meshgrid(1:512,1:512);
k = k.*exp(-2*1i*pi*(fx*x+fy*y));
k = atan2(imag(k),real(k));

f = exp(1i*original_phase);
[x,y] = meshgrid(1:512,1:512);
f = f.*exp(-2*1i*pi*(fx*x+fy*y));
f = atan2(imag(f),real(f));

p = k-f;
high = L*T1*p./(2*pi*D+T*p);
figure(6);set(gcf,'Name','頻移相位還原3D重建模型','NumberTitle','off');surf(high);colormap(jet);shading interp;
%--------------------------------------------------------------------------
%-----------相位展開還原相位------------------------------------------------
unwrapped_org= unwrap(original_phase,pi,2); % 襯底列解卷
unwrapped_org= unwrap(unwrapped_org,pi);    % 襯底行解卷

unwrapped_phase= unwrap(wrapped_phase,pi,2);    % 物體列解卷
unwrapped_phase= unwrap(unwrapped_phase,pi);    % 物體行解卷

delta = unwrapped_phase-unwrapped_org; % 相位差
H = L*T1*delta./(2*pi*D+T*delta);    % 計算高度
figure(5);set(gcf,'Name','相位解卷3D重建模型','NumberTitle','off');surf(H);colormap(jet);shading interp;
%--------------------------------------------------------------------------


    function output = halfball(r,x,y)  % 半球繪製函式
        temp = zeros(512);
        for m = 1:512
            for n = 1:512
                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
    function output = fringe_modulation(u,v,phi,delta_phi) % 物體產生相移條紋函式(x頻率，y頻率，原條紋初相，相移量)
        for a = 1:512
            for b = 1:512
                img(a,b) = (1+sin(u*a+v*b+delta_phi(a,b)+phi))/2;
            end 
        end
        output = img;
    end

  function output = sin_wave(u,v,phi)   % u、v為角頻率,phi為初相
     for m = 1:512
        for n = 1:512
          temp(m,n) = (1+sin(u*m+v*n+phi))/2;
        end
     end
    output = temp;
  end