1. 写在前面

在前面的博客中，介绍了采用同心双圆环靶标下的相机标定算法与代码实践。在这篇博客中，博主将介绍结构光相机模型与逆相机法投影仪标定。完成本篇博客的学习内容后，你将收获投影仪逆相机法标定和结构光相机模型算法与代码实践经验。

本系列博客的完整项目代码皆位于博主的Github项目SLMaster中：

https://github.com/Practice3DVision/SLMaster

动动你的小指头给个Star并follow博主吧！你的支持是博主不懈的动力！

关注公众号，每天推送不迷路！更多技术博客在公众号与你分享！

2. 结构光系统模型

结构光系统模型也就是以数学形式表达条纹投影系统中的关系，一般将结构光系统模型分为相位高度模型和三角测量模型。在相位高度模型中，物体的高度是相对于参考平面计算的，高度是相位的函数。在三角测量模型中，将投影仪视为相机，从而构建立体视觉并通过三角测量原理恢复物体的深度。

博主对上述两种系统模型都标定实验过，个人认为三角测量模型比较方便好用，具体原因有以下几点：

• 相位高度模型标定略微麻烦点儿，具体体现在部分方法需要高精度位移机构
• 相位高度模型计算所得出的参数不具有直观性，无法通过标定参数判断标定是否准确
• 相位高度模型受景深影响较高点儿

所以博主在这里推荐三角测量模型作为结构光系统模型。在前述文章从0到1搭建一套属于你自己的高精度实时结构光3D相机（1）：硬件搭建中已经给出了关于三角测量模型的参考文献，现简单说一下其原理，如图1所示[1]：

图1. 三角测量原理

从图中可见，投影仪平面上的一点$x^p$与其光心$O_p$形成的射线与相机平面上的一点$x$与其光心$O_c$形成的射线相交于空间点$X$。因此，如果我们能确定相机平面上的点与投影仪平面上的点之间的对应关系，那么，我们就能够恢复出物体的深度信息。

所以关键问题就回到如何确定相机平面上的点与投影仪平面上的点之间的对应关系。

我们知道，投影仪平面上的点都是已知的，因为是我们通过烧录条纹图案确定的。相机通过接收投影仪投射的光线从而形成图像，如果我们能在相机成像平面上将条纹编码信息重新恢复出来，那么我们也就知道相机平面上的点与投影仪平面上的点之间的对应关系了！

是不是有点子熟悉呢，对的，博主前两篇博客从0到1搭建一套属于你自己的高精度实时结构光3D相机（6）：相移互补格雷码解码所提及到的相移互补格雷码正好能完成这个任务！绝对相位说白了其实就是投影仪的$u$或$v$坐标！

垂直方向的相移互补格雷码条纹能够确定相机平面与投影仪平面之间的有关$x$坐标的映射关系，而水平方向的相移互补格雷码条纹则能够确定相机平面与投影仪平面之间的有关$y$坐标的映射关系。因此，我们使用两组方向的相移互补格雷码就能确定出$(x_p,y_p)$与$(x_c,y_c)$之间的对应关系了，这种思路可以用图2来表示：

图2. 通过两组方向的条纹确定相机与投影仪之间的映射关系

问题来了！如果我们已经能确定相机平面上的点与投影仪平面上的点之间的对应关系了，那么我如何通过数学关系恢复深度呢？

让我们来推导一下吧！

根据相机成像模型，我们知道任意空间上的一点$P\left[X_w{Y}_w{Z}_w\right]$可以通过下式投影到相机平面上的点$\left[u_c{v}_c\right]$：

$$Z_c\left[u_c{v}_{c}1\right]=K_c\left[R_c|T_c\right]\left[X_w{Y}_w{Z}_{w}1\right]$$

式中，$Z_c$是相机坐标系下点$P$的深度，$K_c$是相机内参，$R_c$是世界坐标系到相机坐标系的旋转矩阵，$T_c$是世界坐标系到相机坐标系的平移矩阵。

由于投影仪的成像模型也可以用相机成像模型来表示，因此，对于投影仪成像平面上的一点$\left[u_p{v}_p\right]$，我们也有下式成立：

$$Z_p\left[u_p{v}_{p}1\right]=K_p\left[R_p|T_p\right]\left[X_w{Y}_w{Z}_{w}1\right]$$

式中，$Z_p$是投影仪坐标系下点$P$的深度，$K_p$是相机内参，$R_p$是世界坐标系到相机坐标系的旋转矩阵，$T_p$是世界坐标系到相机坐标系的平移矩阵。

对于矩阵$K\left[R|T\right]$，我们可以用投影矩阵$P=K\left[R|T\right]$来表示。因此，上述两式可以改写为：

$$Z_c\left[u_c{v}_{c}1\right]=P_c\left[X_w{Y}_w{Z}_{w}1\right]$$

$$Z_p\left[u_p{v}_{p}1\right]=P_p\left[X_w{Y}_w{Z}_{w}1\right]$$

通常，我们将世界坐标系定义为与相机坐标系相重合（见注释），则可通过变换并联立上述两式得：

$$A{\left[X_w{Y}_w{Z}_{w}1\right]}^T=0$$

此时$R_c$为单位矩阵，$T_c$为零矩阵，$R_p$为相机坐标系到投影仪坐标系的旋转矩阵，$T_p$为相机坐标系到投影仪坐标系的平移矩阵。

式中，
$$A=\left[\begin{array}{c}{u}_c{P}_3^c-P_1^c\\ v_c{P}_3^c-P_2^c\\ u_p{P}_3^p-P_1^p\\ u_p{P}_3^p-P_2^p\end{array}\right]$$
其中， $P_i^c$表示相机投影矩阵$P_c$的第$i$行，$P_i^p$表示投影仪投影矩阵$P_p$的第$i$行。

因此，对式上式解方程即可得到$\left[X_w{Y}_w{Z}_w\right]$（也就是相机坐标系下点的空间坐标）。仔细观察上式可以发现，它其实是一个冗余方程（见下注释），因此我们可以通过只投影一组方向的条纹即可恢复空间三维信息。

观察矩阵A的表达式，相机成像平面上的一点$(u_c，v_c)$仅与投影仪成像平面上的一点的单方向坐标$u_p$有关（水平方向条纹则为$v_p$）。

因此，对于相机上一点$(u_c，v_c)$，他的空间三维坐标可以通过下式进行计算：

$$\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\end{array}\right]=\left[\begin{array}{ccc}{p}_{11}^c-p_{31}^c{u}^c& p_{12}^c-p_{32}^c{u}^c& p_{13}^c-p_{33}^c{u}^c\\ p_{21}^c-p_{31}^c{v}^c& p_{22}^c-p_{32}^c{v}^c& p_{23}^c-p_{33}^c{v}^c\\ p_{11}^p-p_{31}^p{u}^p& p_{12}^p-p_{32}^p{u}^p& p_{13}^p-p_{33}^p{u}^p\end{array}\right]\left[\begin{array}{c}{p}_{34}^c{u}^c-p_{14}^c\\ p_{34}^c{v}^c-p_{24}^c\\ p_{34}^p{u}^p-p_{14}^p\end{array}\right]$$

3. 逆相机法标定

既然我们已经推导出能够恢复空间三维信息的结构光系统模型，那么，我们根据上述方程计算即可。但是有一个问题，投影矩阵$P$依赖于内参与外参，相机内参我们可以通过前篇博客相机标定算法得到，相机外参是已知的，那么，投影仪内参与投影仪外参呢？

还是那句话，投影仪成像模型可以简单的看作是相机成像模型！

那么我们就又回到了前述相机标定中的流程，如何查找已知信息靶标的特征点！可是，投影仪无法像相机一样进行拍照，那如何获取投影仪所拍摄的图片呢？

事实上，我们仅需要靶标特征点在投影仪成像平面上的位置信息即可。那么，我们就可以通过前述的两组相移互补格雷码方法，先从相机成像平面知道这些特征点的位置，然后，将这些特征点映射到投影仪成像平面即可！仔细想想，还真是那么回事，照着逻辑往后推就行了！

因此，我们可以得到以下投影仪标定算法流程，也将这种算法称为逆相机法：

• 通过相机拍摄靶标图片
• 在相机成像平面下通过靶标特征点提取算法提取特征点
• 通过绝对相位映射将相机成像平面下的特征点映射至投影仪成像平面
• 与相机标定类似，连续变换多个位姿并重复上述过程
• 使用张氏标定法进行内参、外参估计

哇塞！是不是瞬间明了了呢？具体效果如何呢？请看实验环节！

4. 代码实践与实验

1）代码实践
我们首先来看看代码如何编写。由于投影仪标定需要两组方向的相移互补格雷码条纹，如果一组一组的存入则需要较多的内存，并且也无法在标定过程中知道当前的摆放位姿在投影仪成像平面上是否合理，特征点提取精度情况等等。所以，博主将其该环节编写为在线式的投影仪标定，使用者可以在标定过程中实时查看当前位姿靶标情况。

投影仪标定函数位于SLMaster项目CalibrateEngine.cpp中的captureOnce函数，让我们来看看具体的代码：

bool CalibrateEngine::captureOnce() {
    //清除当前标定误差分布
    if (!projectorErrorDistributes_.empty()) {
        projectorErrorDistributes_.clear();
    }
	//获取结构光相机并准备
    auto slCamera = CameraEngine::instance()->getSLCamera();

    std::string camManufactor, leftCameraName, rightCameraName, colorCameraName, dlpEvmName;
    slCamera->getStringAttribute("2D Camera Manufactor", camManufactor);
    const device::camera::CameraFactory::CameraManufactor manufator = camManufactor == "Huaray" ? device::camera::CameraFactory::Huaray : device::camera::CameraFactory::Halcon;
    slCamera->getStringAttribute("Left Camera Name", leftCameraName);
    slCamera->getStringAttribute("DLP Evm", dlpEvmName);

    device::projector::Projector* projector = slCamera->getProjectorFactory()->getProjector(dlpEvmName);

    device::camera::Camera* leftCamera = nullptr, *rightCamera = nullptr, *colorCamera = nullptr;
    leftCamera = slCamera->getCameraFactory()->getCamera(leftCameraName, manufator);
    if(slCamera->getStringAttribute("Right Camera Name", rightCameraName)) {
        rightCamera = slCamera->getCameraFactory()->getCamera(rightCameraName, manufator);
    }
    if(slCamera->getStringAttribute("Color Camera Name", colorCameraName)) {
        colorCamera = slCamera->getCameraFactory()->getCamera(colorCameraName, manufator);
    }
    //投影一次
    projector->project(false);

    const int imgSizeWaitFor = CameraEngine::instance()->getNumberAttribute("Total Fringes");
    const int totalExposureTime = (CameraEngine::instance()->getNumberAttribute("Pre Exposure Time") + CameraEngine::instance()->getNumberAttribute("Exposure Time") + CameraEngine::instance()->getNumberAttribute("Aft Exposure Time")) * imgSizeWaitFor;
    auto endTime = std::chrono::steady_clock::now() + std::chrono::duration<int, std::ratio<1, 1000000>>(totalExposureTime + 1000000);
    while (leftCamera->getImgs().size() != imgSizeWaitFor) {
        if (std::chrono::steady_clock::now() > endTime) {
            leftCamera->clearImgs();
            if (rightCamera) {
                rightCamera->clearImgs();
            }
            if (colorCamera) {
                colorCamera->clearImgs();
            }
            return false;
        }
    }
    //软触发一次获取纹理图片（由于博主这套相机是离焦，如果聚焦状态下，可以用相移图案取平均作为纹理图案）
    leftCamera->setTrigMode(device::camera::trigSoftware);
    leftCamera->setNumberAttribute("ExposureTime", 100000);
    cv::Mat texture = leftCamera->capture();
    if(texture.type() == CV_8UC3) {
        cv::cvtColor(texture, texture, cv::COLOR_BGR2GRAY);
    }
    leftCamera->setNumberAttribute("ExposureTime", CameraEngine::instance()->getNumberAttribute("Exposure Time"));
    //获取当前投影仪条纹烧录记录表，用于确定水平条纹或垂直条纹以及条纹周期等情况
    auto orderTablesRecord = CameraEngine::instance()->getOrderTableRecord();
    static cv::Mat honrizonUnwrapMap, verticalUnwrapMap, normHUnwrapMap, normVUnwrapMap, textureMap;
    const int dlpHeight = CameraEngine::instance()->getStringAttribute("DLP Height").toInt();
    const int dlpWidth = CameraEngine::instance()->getStringAttribute("DLP Width").toInt();
    //依次按烧录的条纹序列进行解相
    for (int i = 0; i < orderTablesRecord.size(); ++i) {
        std::vector<cv::Mat> imgs;
        int index = 0;
        while(index++ < orderTablesRecord[i].patternsNum_) {
            cv::Mat img = leftCamera->popImg();

            if(img.type() == CV_8UC3) {
                cv::cvtColor(img, img, cv::COLOR_BGR2GRAY);
            }

            imgs.push_back(img);
        }
       
        cv::structured_light::SinusCompleGrayCodePattern::Params params;
        params.confidenceThreshold = CameraEngine::instance()->getNumberAttribute("Contrast Threshold");
        params.height = dlpHeight;
        params.width = dlpWidth;
        params.horizontal = !orderTablesRecord[i].isVertical_;
        params.shiftTime = orderTablesRecord[i].shiftTime_;
        params.nbrOfPeriods = std::pow(2, orderTablesRecord[i].patternsNum_ - orderTablesRecord[i].shiftTime_ - 1);
        auto pattern = cv::structured_light::SinusCompleGrayCodePattern::create(params);

        cv::Mat confidenceMap, wrappedPhaseMap, floorMap, unwrapMap, normalizeUnwrapMap;

        pattern->computeConfidenceMap(imgs, confidenceMap);
        pattern->computePhaseMap(imgs, wrappedPhaseMap);
        pattern->computeFloorMap(imgs, confidenceMap, wrappedPhaseMap, floorMap);
        pattern->unwrapPhaseMap(wrappedPhaseMap, floorMap, unwrapMap, confidenceMap > params.confidenceThreshold);

        cv::normalize(unwrapMap, normalizeUnwrapMap, 0, 255, cv::NORM_MINMAX);

        if(orderTablesRecord[i].isVertical_) {
            verticalUnwrapMap = unwrapMap;
            textureMap = texture;
            normVUnwrapMap = normalizeUnwrapMap;
            normVUnwrapMap.convertTo(normVUnwrapMap, CV_8UC1);
            onlineProjVertiPaintItem_->updateImage(QImage(normVUnwrapMap.data, normVUnwrapMap.cols, normVUnwrapMap.rows, normVUnwrapMap.step, QImage::Format_Grayscale8));
        }
        else {
            honrizonUnwrapMap = unwrapMap;
            normHUnwrapMap = normalizeUnwrapMap;
            normHUnwrapMap.convertTo(normHUnwrapMap, CV_8UC1);
            onlineProjHonriPaintItem_->updateImage(QImage(normHUnwrapMap.data, normHUnwrapMap.cols, normHUnwrapMap.rows, normHUnwrapMap.step, QImage::Format_Grayscale8));
        }
    }

    leftCamera->clearImgs();
    if (rightCamera) {
        rightCamera->clearImgs();
    }
    if (colorCamera) {
        colorCamera->clearImgs();
    }

    cv::Mat projectorImg = cv::Mat::zeros(dlpHeight, dlpWidth, CV_8UC3);

    auto calibrator = getCalibrator(projectorCaliParams_.targetType_);
    std::vector<cv::Point2f> featurePoints;

    //提取当前纹理图案下的靶标特征点
    if(!calibrator->findFeaturePoints(textureMap, cv::Size(projectorCaliParams_.rowNum_, projectorCaliParams_.colNum_), featurePoints)) {
        return false;
    }
    else {
        cv::cvtColor(textureMap, textureMap, cv::COLOR_GRAY2BGR);
        cv::drawChessboardCorners(textureMap, cv::Size(projectorCaliParams_.rowNum_, projectorCaliParams_.colNum_), featurePoints, true);
        onlineProjColorPaintItem_->updateImage(QImage(textureMap.data, textureMap.cols, textureMap.rows, textureMap.step, QImage::Format_BGR888));
        //将相机成像平面上的靶标特征点坐标映射至投影仪成像平面上
        std::vector<cv::Point2f> remapPoints;
        for (auto camPoint : featurePoints) {
            auto projPoint = remapProjectorPoint(honrizonUnwrapMap, verticalUnwrapMap, camPoint);
            remapPoints.emplace_back(projPoint);
        }

        projCameraPoints_.emplace_back(featurePoints);
        projectorPoints_.emplace_back(remapPoints);

        cv::drawChessboardCorners(projectorImg, cv::Size(projectorCaliParams_.rowNum_, projectorCaliParams_.colNum_), remapPoints, true);
        onlineProjPaintItem_->updateImage(QImage(projectorImg.data, projectorImg.cols, projectorImg.rows, projectorImg.step, QImage::Format_BGR888));
    }

    std::vector<cv::Mat> groupImg;
    groupImg.emplace_back(normHUnwrapMap);
    groupImg.emplace_back(normVUnwrapMap);
    groupImg.emplace_back(textureMap);
    groupImg.emplace_back(projectorImg);
    projectorCaliImgs_.emplace_back(groupImg);

    if(rightCamera) {
        rightCamera->clearImgs();
    }

    if(colorCamera) {
        colorCamera->clearImgs();
    }

    const QString index = projectorModel_->imgPaths().empty() ? "0" : QString::number(projectorModel_->imgPaths()[projectorModel_->imgPaths().size() - 1].toInt() + 1);
    projectorModel_->emplace_back(index);
    emit projectorModelChanged();

    return true;
}

紧接着，我们再来看看映射到投影仪成像平面下的具体函数remapProjectorPoint：

const cv::Point2f CalibrateEngine::remapProjectorPoint(const cv::Mat& honrizonPhaseMap, const cv::Mat& verticalPhaseMap, const cv::Point2f& camPoint) {
//双线性插值
#ifdef LinearInter
    int index_Y_Upper = std::ceil(camPoint.y);
    int index_Y_Lower = std::floor(camPoint.y);
    int index_X_Upper = std::ceil(camPoint.x);
    int index_X_Lower = std::floor(camPoint.x);
    const float vPLU = verticalPhaseMap.ptr<float>(index_Y_Upper)[index_X_Lower];
    const float vPRU = verticalPhaseMap.ptr<float>(index_Y_Upper)[index_X_Upper];
    const float vPLD = verticalPhaseMap.ptr<float>(index_Y_Lower)[index_X_Lower];
    const float vPRD = verticalPhaseMap.ptr<float>(index_Y_Lower)[index_X_Upper];
    const float tPLU = honrizonPhaseMap.ptr<float>(index_Y_Upper)[index_X_Lower];
    const float tPRU = honrizonPhaseMap.ptr<float>(index_Y_Upper)[index_X_Upper];
    const float tPLD = honrizonPhaseMap.ptr<float>(index_Y_Lower)[index_X_Lower];
    const float tPRD = honrizonPhaseMap.ptr<float>(index_Y_Lower)[index_X_Upper];
    const float vfR1 = (index_X_Upper - camPoint.x) / (index_X_Upper - index_X_Lower) * vPLD + (camPoint.x - index_X_Lower) / (index_X_Upper - index_X_Lower) * vPRD;
    const float vfR2 = (index_X_Upper - camPoint.x) / (index_X_Upper - index_X_Lower) * vPLU + (camPoint.x - index_X_Lower) / (index_X_Upper - index_X_Lower) * vPRU;
    const float tfR1 = (index_X_Upper - camPoint.x) / (index_X_Upper - index_X_Lower) * tPLD + (camPoint.x - index_X_Lower) / (index_X_Upper - index_X_Lower) * tPRD;
    const float tfR2 = (index_X_Upper - camPoint.x) / (index_X_Upper - index_X_Lower) * tPLU + (camPoint.x - index_X_Lower) / (index_X_Upper - index_X_Lower) * tPRU;
    const float verticalPhaseValue = (index_Y_Upper - camPoint.y) / (index_Y_Upper - index_Y_Lower) * vfR1 + (camPoint.y - index_Y_Lower) / (index_Y_Upper - index_Y_Lower) * vfR2;
    const float transversePhaseValue = (index_Y_Upper - camPoint.y) / (index_Y_Upper - index_Y_Lower) * tfR1 + (camPoint.y - index_Y_Lower) / (index_Y_Upper - index_Y_Lower) * tfR2;
//最近邻插值
#else 
    int index_Y = std::round(camPoint.y);
    int index_X = std::round(camPoint.x);
    const float verticalPhaseValue = verticalPhaseMap.ptr<float>(index_Y)[index_X];
    const float transversePhaseValue = honrizonPhaseMap.ptr<float>(index_Y)[index_X];
#endif
//映射至投影仪坐标系
    const float xLocation = (verticalPhaseValue) / CV_2PI * projectorCaliParams_.verticalPitch_;
    const float yLocation = (transversePhaseValue) / CV_2PI * projectorCaliParams_.honrizonPitch_;
    return cv::Point2f(xLocation, yLocation);
}

我们通过双线性插值进一步提高特征点提取精度，效果显著，不信？那你可以试试喔。

2）实验

首先我们需要一块标定板，由于黑色对光的反射率较低，因此棋盘格、圆环靶标容易受到干扰，所以我们使用黑底白色前景的圆点标定板进行投影仪标定。

打开我们的SLMaster软件，按如下步骤进行，你就能得到你想要的结果了！

图4. 软件标定流程

然后我们可以将标定文件进行保存。前面说过，我们的三角测量模型是可以通过标定参数确定标定的好不好的！则么实施的呢？我们打开保存好的标定文件，查看Tlr参数，如图5所示。然后实际测量相机与投影仪之间的相对位置，对比一下就知道了！

图5. 标定所得的相机到投影仪的平移矩阵

通过结构光系统模型进行三维重建的实验呢？

这一部分环节将在之后的多视几何约束相关算法进行。博主不单独使用结构光系统模型进行三维重建，如果你有需要，可以调用博主写的方法进行实验（见reverse_camera.cu）。

3. 总结

在这篇博客中，博主介绍了结构光系统模型与逆相机法投影仪标定算法 ，在下一小节中，博主将开始介绍最佳脉冲宽度调制方法的二值正弦条纹生成算法。

参考文献

[1] Feng, Shijie & Zuo, Chao & Zhang, Liang & Tao, Tianyang & Hu, Yan & Yin, Wei & Qian, Jiaming & Chen, Qian. (2021). Calibration of fringe projection profilometry: A comparative review. Optics and Lasers in Engineering. 143. 106622. 10.1016/j.optlaseng.2021.106622.

本系列文章将持续更新，如果有等不及想要获得代码的小伙伴，可访问博主Github中的SLMaster项目，动动你的小手指，follow and star！你们的关注是博主持续的动力！

Github：https://github.com/Practice3DVision

QQ群：229441078

公众号：实战3D视觉

扫一扫关注公众号