研究光度立體法階段性小結和最佳化（可20ms獲取4個2500*2000灰度圖的Normal Map）。

Imageshop發表於2022-11-25

這個東西是我接觸的第一個非2D方面的演算法，到目前為止其實也沒有完全搞定，不過可能短時間內也無法突破。先把能搞定的搞定吧。

這個東西也有一大堆參考資料，不過呢，搜來搜去其實也就那些同樣的東西，個人覺得就屬這個文章最經典，既有說明，也有圖片，還有程式碼：

Photometric Stereo　　　　Chaman Singh Verma and Mon-Ju Wu

　 https://pages.cs.wisc.edu/~csverma/CS766_09/Stereo/stereo.html

　　另外，github上也應該有一些參考的資料吧，我主要參考的是 https://github.com/chaochaojnu這個中國小哥的部落格。

　　目前為止，我只實現了提取Albedo、Normal Map和Normal Vectors三個結果。

　　從硬體上講，這演算法應該需要一個固定位置的相機（應該是要和目標垂直吧），以及至少3個以上的平行光源，一般實際上可能需要至少4個以上的光源吧，然後每個光源單獨打光，單獨拍一張圖片，共得到N個不同的圖片，然後根據這個N個圖片，合成一個結果圖以及得到額外的梯度和高度資訊。

　　在Halcon中，有對應的photometric_stereo運算元實現該功能，該運算元除了要提供N個圖片，還需要提供 Slants和Tilts兩個引數，你去看他們的英語翻譯，其實都是傾斜角，個人理解Tilts就是光源在XY平面投影時和X軸的夾角，而Slants就是光源和XY平面的夾角。

　　在我剛剛提供的兩個連結裡，他們都不是直接提供 Slants和Tilts，而是直接利用標準物體在對應光源下拍照，得到幾幅標準影像，然後由標準影像的畫素值推算出對應的歸一化光源向量，這個方法也是不錯了，省去了相機和光源位置的標定。

　　有了這些引數，就可以進行演算法的執行了，對於Normal Map的獲取，在Photometric Stereo這個文章裡有一大堆推導，開始看不懂，慢慢的又覺得懂了，然後又有點懵逼，接著折騰又似乎清晰了。

　　其實不用管那麼多，我們看看Photometric Stereo給出的NormalMap.m程式碼裡的細節吧：

   for i = 1:nrows
        for j = 1:ncols
            if( maskImage(i,j) ) 
            for im = 1:numImages
                I(im) = double(grayimages(i,j,im));
                end
            [NP,R,fail] = PixelNormal(I, lightMatrix);
            surfNormals(i,j,1) =  NP(1);
            surfNormals(i,j,2) =  NP(2);
            surfNormals(i,j,3) =  NP(3);
            albedo(i,j)        = R;
            end
        end
   end

　　這裡的surfNormals就是對應上圖中的Normal Vectors，albedo就是反射率圖。

　　lightMatrix是多個光源的向量，I是多個圖對應的畫素值，這裡的關鍵在於PixelNormal函式。

unction [N,R, fail] = PixelNormal(I, L)
   fail = 0;
   I  = I';
   LT = L';
   A  = LT*L;
   b  = LT*I;
   g  = inv(A)*b;
   R  = norm(g);
   N  = g/R;
   if( norm(I) < 1.0E-06) 
       fprintf( ' Warning: Pixel intensity is zero \n' );
       N(1) = 0.0;
       N(2) = 0.0;
       N(3) = 0.0;
       R    = 0.0;
       fail = 1;
   end
end

　　仔細看這個函式，其實就是上面工時最後一部分的直接實現，什麼轉置、乘積、求逆、歸一化等等。

這個M程式碼要稍微修改才可以執行，我嘗試了下只是執行獲取Normal Map這一塊，4個500*500的大小的灰度圖，需要大概2s的時間，這個時間其實對於工程專案本身來說是沒有任何意義的。所以也驗證了一句話，matlab只是實驗裡的工具。

當然，M程式碼本身也常常只是用來作為驗證一個演算法的結果是否正確的第一步而已。

　　要真正讓他變得有意義，像這麼大的圖，一個合理的處理時間是不大於2ms，這個最佳化其實也不是很困難，只要你仔細的看看PixelNormal函式里的資料和程式碼。

　　裡面最耗時的其實是inv(A)，矩陣求逆需要用到LU分解，很少麻煩，但是實際上，我們注意到這裡的A值其實是由L一個值決定的，L是什麼，是光源的向量矩陣，這意味著什麼，就是L是不變的。沒有必要每次在迴圈裡計算矩陣的逆的。這是最重要的一個問題和耗時點所在。

　　瞭解到這個問題後，我們其他的最佳化手段就是程式碼層次上的了，比如用C++寫演算法、把PixelNormal這個小函式直接整合到迴圈內部、使用SIMD指令加速等等。

　　在幾個matlab程式碼裡，還要求提供一個mask圖，有這個的原因是很多立體光度法拍攝的圖片其實有很多黑色的部分或者說可以確定不是目標的部分，這部分如果處理，是會拖累演算法的速度的，因此用個標定好的MASK去刪除他，這也無可厚非，不過在halcon裡似乎沒有這個引數。

　　我目前也就只研究到這裡，至於後面的深度圖或者說是高度圖的實現，文章裡提供的都是解一個很大的稀疏矩陣，這個已經超出了我所能自行程式設計的範圍。暫時沒有能力去解決了。

在Halcon中，利用光度立體法去實現一些檢測目標的一個重要應用是透過photometric_stereo運算元獲取對應的gradient，然後在利用derivate_vector_field 獲得梯度的平均曲率場，我目前還不明白這個gradient到底代表了什麼值，是上面的M程式碼裡的surfNormals向量嗎？有沒有哪位朋友知道呢。

　　透過和halcon比較，目前獲取的反射率圖，基本還是差不多正確的，比如下面幾個halcon的測試圖：

研究光度立體法階段性小結和最佳化（可20ms獲取4個2500*2000灰度圖的Normal Map）。