研究光度立体法阶段性小结和优化（可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的测试图：