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

這個東西是我接觸的第一個非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的測驗圖：