基於結構光的相移法三維重建

1、基本原理：
    正弦條紋打在三維物體上，CCD記錄到的條紋因爲受到三維物體高度的調製而發生扭曲，扭曲的條紋（deformed fringe）實質上爲原始條紋在物體具備高度存在的位置有了附加相位，各點的相位表現爲由CCD圖像採集得到的被調製的條紋數字圖像的灰度值。經過扭曲的條紋和原始條紋比對計算得出相位變化值。

又已知投影儀、CCD和物體的具體位置和之間的距離，利用數學關係可求出對應點的高度值，實現3D重建。公式推導以下：（下圖來自[1]）

 

2、四步相移法
    

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

    在利用CCD進行數字採集過程當中，因爲環境等各類因素的影響，以及物體在空間上對光的反射程度和與CCD距離不一樣，不可以獲得理想的均勻而又高對比度的條紋，也就不能獲得最真實的相位信息。

運用四步相移法處理此問題（此方法參考[1]）,過程分析以下：

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

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

表達式以下：

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

聯立解得  

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

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

一、 相位包裹的含義：

如圖所示（圖片來自[1]），下圖爲真實相位值，上圖爲被包裹的相位。即不在反正切函數(-π,π)的主值區間內的相位值被加上或減去n個2π（n爲整數），使其值被調整了主值區間(-π,π)內，稱爲相位包裹。被包裹的相位（wrapped phase）並非真正的相位信息，須要進行相位展開，還原出其原始的相位值。
    此處有另外一個重點，數學中經常使用的反正切函數是二象限的反正切函數，值域僅是(-π/2,π/2)，而此處咱們運用的反正切函數是四象限的反正切函數，值域是(-π,π)。區別以下：
二象限反正切函數（matlab中函數名爲atan）輸入參數爲正切值，返回正切值僅有兩個第一第四兩個象限的角度值，即(-π/2,π/2)。而四象限反正切函數（matlab中函數名爲atan2）輸入參數爲y和x，返回正切值包括四個象限，即(-π,π)。

 四象限反正切考慮到正切值的由來是兩個量的商，多考慮了x爲負時的另外一種狀況，便可以將值域擴展一倍，延長到(-π,π)。

表達式以下：（表達式來自Wikipedia）

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

二、解卷/相位展開/相位解包裹操做（phase unwrapping）
    上面講述的相位包裹表現爲相位在主值區間內增長到超過主值區間的上限的時候，即超過π時，將發生跳變，也就是說繼續增長相位值會跳變到-π。在總體提取的圖像中沿着相位增長的方向表現爲鋸齒波的形狀。以下圖爲沒有被物體調製過的條紋沿着條紋的方向橫截，

被包裹的相位圖：

相位展開最根本的思想就是找出相位發生跳變的地方，補上被減去的n個2π，使跳變出相位值和周圍的值變得連續、平滑。方法不少，請參考[1]，文中詳細講解質量法。

本次使用matlab仿真直接使用了matlab中相位解卷函數unwrap(P,tol,1)。其中P爲執行相位解卷的向量或矩陣；tol爲相位跳變限度（jump tolerance），即差值大於tol被斷定爲發生跳變，默認爲pi，此處我使用的時候填pi；第三個參數表示進行相位解卷操做的維度，向量只有一個維度，矩陣有兩個維度，則1表示行解卷，2表示列解卷。上一張圖相位展開完成後相位分佈以下圖，

能夠看到不存在跳變（鋸齒）。

4、傅里葉變換法還原相位

    運用此方法不進行相位解卷即可以得到真實的相位信息。原理以下：（此方法和原理公式來自[2]）
被包裹的相位表示爲

其中φ(x,y)爲真實的相位信息，fx和fx標誌條紋頻率。
將被包裹的相位表示成複數形式
二維傅里葉變換

頻移後進行逆傅里葉變換

求複數表示的相位的相角便可獲得真實相位值

兩次傅里葉變換亦能夠應用傅里葉變換的頻移性質實現：

 此方法就是利用傅里葉變換頻移實現濾波操做，濾掉條紋部分的信息。[2]對該方法提升精度的一些辦法進行了闡述，具體參考原文。
    此方法雖然能夠直接獲得相位信息而不須要進行解卷，且精度較高，可是由於沒有解卷，若物體高度換算大於2π的相位差，則或出現重建變形，較高部分凹陷下去的狀況。

5、matlab仿真結果

仿真基本能夠重建出半球的形狀。半球最高點高度理論上應該是半徑的長度，因爲計算所用的正弦波週期是經過離散傅里葉變換頻譜峯值提取出的頻率，提取出的頻率存在偏差。

 

參考文獻：

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

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

6、仿真matlab代碼

<span style="font-size:14px;">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
</span>
  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

end