光度立体算法实现（opencv_c++）

一. 理论讲解部分

这边主要我也不太会，就只能给几个参考博客链接了
我只能根据得出的结论来讲，去讲解和实现算法

• https://pages.cs.wisc.edu/~csverma/CS766_09/Stereo/stereo.html – 一篇很好的数学原理讲解的英文文章
• https://zhuanlan.zhihu.com/p/205778332 – 知乎文章

另外：相关的数学知识链接

• 线性代数基础知识——常见矩阵的概念及其关系
• 线性代数基础–矩阵
• 一文让你彻底搞懂最小二乘法（超详细推导）
• 矩阵乘以其矩阵转置

二. 依赖求解结果手推实现

1. 首先我们可以知道，一幅图像的可见度有两个主要原因
                  I = p * L * l
    I: 为亮度值
    p: 为材料本身的反照率
    L: 光源方向向量（单位向量）
    l: 光源强度（一般为1就行）
        
2. 根据上述公式和之前的理论基础，进一步的公式可以写为：
                 I[(x,y),j] = p(x,y) * L(j) * l
    以上公式意为: 当前图像某一点(x,y)的亮度值为当前材质(x,y)点的光照率与当前光源方向向量的乘积
        
3. 但上述公式并没有考虑方向性，这个点(x,y)表面是朝哪倾斜的？光照打在这个点上，是正着照？斜着照？背着照？
   考虑上述后，最终公式变为：
    			I[(x,y),j] = p(x,y) * n * L(j) * l  (注：这里的n是点的单位法向量，L为当前光源方向向量)
        n: 法向量控制了当前点对光照的“响应角度”
        n * L(j)的乘积为：表面朝向与光照方向的夹角余弦值

$$I=\rho \cdot (n\cdot L)$$

我们现在要求的就是反照率p:物体表面的“固有亮度”

如何去求取呢？
因为n为表面的单位向量，p为常数，我们可以将np看作整体：
				令np = g， 进而：I = L*g
    因为L不是方矩阵，所以我们在其左侧同时乘以L矩阵的转置，这样就在g的左侧得到一个方矩阵。
转换下面截图的字母：自己理解哈，最后我们也可以相应的得到g向量的结果
最终是求取p（反照率），g又是n向量与p的乘积,且n又是单位向量，即我们求取模值就是对应的反照率值

三. 求解

1. 为什么至少需要三张图（不同光源方向）

光源方向L: 可以分解为3个方向，xyz
同理，我们的g向量，可以分为3个方向：xyz
				I[(x,y),j] = p(x,y) * n * L(j) * l
    将（x,y）去除，看成图像整体：
				I[j] = p * n * L(j) * l  =  g * L(j)
    由此可得，不同光源方向所拍的照片，g是不变的，只与本身的光源方向有关，我们将方向分解出来：

如上述图：
	I: 可以理解为像素值（这个我们拍照下来就有）
    l: 当前图像的光源方向向量（需要分解为xyz三个方向）  这个我们未知，但能求取知道（通过Slant和Tilt）
    g: 就是我们需要求的量
三个未知数（gx,gy,gz）,那当然我们需要三张图了

2. 求取光源方向的不同分量

我们知道：在拍取不同光源方向的图，需要记录两个角度，一个是Slant, 一个是Tilt
	Slant: 光束方向与相机中轴线的角度
	Tilt: 俯视看，光束投影与起点的夹角

• 求取lx,ly,lz分量

  

• 有了光照图，光源方向的不同分量，三张图足够求取g了
  

3. 代码实现求解反照率图

#include <iostream>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

#define Tilts(x) Tilts
#define Slants(x) Slants

int main()
{
    const int NUM_IMGS = 4;
    // 定义图片路径
    const string MODEL = "./images/blister_";
    // 定义4张图的各个拍照角度
    vector<float> Tilts(blister_tilts) = { 6.1,95.0,-176.1,-86.8 };
    vector<float> Slants(blister_slants) = { 41.4,42.6,41.7,40.9 };
    // 存储图像
    vector<Mat> modelImages;
    // 反照率图
    cv::Mat AlbedoMap;
    for (int i = 0; i < NUM_IMGS; i++) 
    {
        std::string path = MODEL + to_string(i) + ".png";
        cv::Mat Model = imread(path, IMREAD_GRAYSCALE);
        if (Model.empty() == true)
        {
            std::cout << "Read Image " << path << " failed!" << std::endl;
            return -1;
        }
        modelImages.push_back(Model);
    }
    auto ret = PhotometricStereo(modelImages, AlbedoMap, NUM_IMGS, Slants, Tilts);
    if (0 != ret)
    {
        std::cout << "Algo Run failed!" << std::endl;
        return -1;
    }
    // 从原图中挑一张出来看对比
    cv::Mat SrcImage = imread(MODEL + "1.png", IMREAD_GRAYSCALE);
    cv::imshow("SrcImage", SrcImage);
    //反照率图
    cv::imshow("Albedomap", AlbedoMap);
    cv::waitKey(0);
    return 0;
}

uint32_t PhotometricStereo(const vector<Mat>& srcImages,
	Mat& Albedo,
	int ImageCount,
	vector<float> Slants,
	vector<float> Tilts
    )
{
    if (srcImages.size() < ImageCount ||
        Slants.size() < ImageCount ||
        Tilts.size() < ImageCount)
    {
        return -1;
    }

    //【1】准备数据
    const int NUM_IMGS = 4;
    vector<Mat> modelImages;
    // 四光源、xyz三个轴
    Mat Lights(NUM_IMGS, 3, CV_32F);//保存光源矩阵  （4*3）的矩阵
    // 4: 4张图      3: 3光源分量（xyz）
    //【2】通过给定角度计算Lights
    for (int i = 0; i < NUM_IMGS; i++) 
    {
        cv::Mat Model = srcImages[i];
        if (Model.empty())
        {
            return -1;
        }
        //计算光源各方向向量分量
        /* x,y are swapped here */
        float z = std::cos(CV_2PI * Slants[i] / 360.0);
        float xy = std::sin(CV_2PI * Slants[i] / 360.0);
        float x = std::sin(CV_2PI * Tilts[i] / 360.0) * xy;
        float y = std::cos(CV_2PI * Tilts[i] / 360.0) * xy;
	
        Vec3f L(x, y, z);
        Lights.at<float>(i, 0) = L[0];
        Lights.at<float>(i, 1) = L[1];
        Lights.at<float>(i, 2) = L[2];
        modelImages.push_back(Model);
    }

    const int height = modelImages[0].rows;
    const int width = modelImages[0].cols;
    //【3】计算反照率ρ和表面法向量N
    /* light directions, surface normals, p,q gradients */
    cv::Mat LightsInv;// 定义L的逆矩阵，L-1
    // 直接使用api方法求取逆矩阵
    cv::invert(Lights, LightsInv, cv::DECOMP_SVD);

    cv::Mat Normals(height, width, CV_32FC3, cv::Scalar::all(0));//法向量图
    cv::Mat AlbedoMap(height, width, CV_32F, cv::Scalar::all(0));//反照率图
    /* estimate surface normals and p,q gradients */
    for (int x = 0; x < width; x++) 
    {
        for (int y = 0; y < height; y++) 
         // 这里是循环每一点(x，y).因为每个点都有属于自己的反照率
        {
            Vec<float, NUM_IMGS> I;
            for (int i = 0; i < NUM_IMGS; i++) 
            {
                // 取出n张图对应位置的像素值：I: [img1(x,y),... imgi(x,y)]
                I[i] = modelImages[i].at<uchar>(Point(x, y));
            }
		   // 对应最下面截图，直接用逆矩阵与每亮度矩阵相乘得到 g
            cv::Mat n = LightsInv * cv::Mat(I);
            // 求取反照率模值，就是对应位置点的亮度值 
            float kd = sqrt(n.dot(n));
            if (kd > 0) { n = n / kd; }
            if (n.at<float>(2, 0) == 0) { n.at<float>(2, 0) = 1.0; }
            int legit = 1;
            // 下面除以255是归一化,将反照率赋值给（x,y）点
            AlbedoMap.at<float>(cv::Point(x, y)) = kd / 255;
        }
    }
    //【4】结果图像输出
    cv::normalize(AlbedoMap, AlbedoMap, 0, 255, cv::NORM_MINMAX, CV_8U);
    Albedo = AlbedoMap;
    return 0;
}

代码运行所需图片点此链接