简介：面向摄影测量、三维重建与无人机航测方向的学习者，这份VC++工程实现了直接线性变换（DLT）算法，用于把二维像点坐标与物方三维坐标关联起来，完成投影中心、内参等参数的线性求解，适合已具备一定摄影测量与C++基础、希望理解算法落地过程的中高级读者。压缩包共38个文件、约323KB，以cpp、h源码与dsw/dsp工程文件为主体，另含txt控制点坐标与说明文本、rc资源脚本及obj、pdb等编译中间文件，可直接用VC6打开编译调试；源码按对话框、主程序与标准头文件分模块组织，关键步骤配有注释，便于对照公式梳理数据输入、方程构建、参数估计与结果输出的完整链路。目前已有622人学习。读者可从中获得一套可运行的DLT实现范例、控制点样例数据与工程配置参考，用于课程实验、算法复现和航测项目中的坐标解算练习，也能借助调试产物排查编译与运行问题。

1. 摄影测量 DLT 程序：六个控制点反算相机参数的一套实用方案

上周帮某实验室核对一组近景摄影测量的标定结果，拿到的资料很朴素：一份写着六个控制点物方坐标的 txt、一份对应的像点像素坐标，外加一个 VC++ 工程。要求是把相机的主点、焦距、旋转和平移全解出来，重投影残差压在 0.3 像素以内。核心就是 DLT——直接线性变换，把共线条件方程展开成线性最小二乘，不用初值、不用迭代、不用碰非线性优化。这份 VC++ 程序代码做的就是这件事：读控制点、组设计矩阵、SVD 求解、反解内参、输出残差表。适合手里有控制点坐标、想快速验一遍相机位姿的人，也适合把 DLT 当成光束法平差初值来源的人。

2. 工程结构与编译环境：VC++ 里 MFC 外壳与控制台内核怎么分层

拿到一个 VC++ 摄影测量工程，第一件事不是点运行，而是看它分了几个工程、算法代码有没有和界面绑死。绑死的那种最难改：想批量跑一百张影像，得在对话框里点一百次。合理的三层划分是界面层、数据层、算法层，算法层不引用任何 MFC 头文件，能单独编成一个静态库，这样控制台入口和对话框入口共用同一份求解代码。

常见做法是把算法抽成DltCore静态库，MFC 只负责选文件、填控制点、把残差表画出来。这样二次开发时换界面不用动算法，做批量实验时直接调库函数就行。

2.1 目录结构与各文件职责

一个能跑的三维 DLT 工程，目录大致长这样：

dlt_solver/ ├─ DLTApp/ MFC 对话框工程，只负责选文件、填参数、显示残差表 │ ├─ DLTAppDlg.cpp/.h 界面事件，调用 DltCore 的接口 │ ├─ res/ 图标与资源 │ └─ DLTApp.vcxproj ├─ DLTConsole/ 控制台入口，适合脚本化批量跑 │ └─ main.cpp ├─ DltCore/ 纯 C++ 算法静态库，不引用 MFC │ ├─ include/dlt.h 对外接口与 ControlPoint 结构体 │ ├─ src/dlt.cpp 构建设计矩阵 + SVD 求解 │ ├─ src/matrix.cpp 3x3/3x4 小矩阵工具：RQ 分解、求逆、行列式 │ ├─ src/evaluate.cpp 重投影残差、RMS、最大偏差统计 │ └─ DltCore.vcxproj ├─ data/ │ ├─ control_points.txt 物方 X Y Z + 像方 x y，一行一个点 │ └─ result_template.txt 输出格式模板 └─ DLT.sln

dlt.h里最关键的是一个结构体和两个函数声明，结构体字段顺序决定了 txt 文件的列顺序，改格式的时候两边一起改，别只改一边：

// dlt.h —— 对外暴露的最小接口 struct ControlPoint { double X, Y, Z; // 物方坐标，单位统一为米（或统一为毫米） double x, y; // 像点坐标，单位统一为像素 }; // 求解 3x4 投影矩阵，rms 返回重投影中误差（像素） bool solveDlt(const std::vector<ControlPoint>& pts, cv::Mat& P, double& rms); // 从 P 中拆出内参 K、旋转 R、平移 t 和摄影中心 C bool decomposeP(const cv::Mat& P, cv::Mat& K, cv::Mat& R, cv::Mat& t, cv::Point3d& C);

数据层单独拎出来是有原因的。控制点文件经常被 Excel 编辑过，制表符、空格、全角逗号混着来，读文件的代码比算法代码更容易出问题。data/control_points.txt里我一般要求第一行写点数，后面每行五列，行尾不留空格。读入时先数一下实际行数对不对得上，对不上直接报错，不要静默截断——静默截断最坑人，明明给了十二个点，只读了六个，残差看着还挺小。

2.2 依赖配置：OpenCV 与纯手写矩阵两条路线

矩阵求解有两条路，选哪条取决于你要不要把这个工程拷到没有第三方库的机器上跑。

路线 A 的编译配置踩坑集中在运行时库。OpenCV 官方包默认/MD，如果你的工程在“代码生成 → 运行时库”里选了/MT，链接时会报一堆LNK2005 已经在 xxx.obj 中定义。翻车基本都是这个原因，把整个解决方案的运行时库调成一致的/MD（Debug 用/MDd）就好。

Debug 和 Release 的 dll 名字不一样，opencv_world4xxd.dll带 d 的是 Debug 版，Release 目录里放了带 d 的 dll 会直接起不来，报“无法定位程序输入点”。我一般把两个 dll 都拷到输出目录，省得来回换。

路线 B 也不是完全从零写。11 个未知数、至少 12 个方程，把A^T A求出来做 Jacobi 特征分解取最小特征值对应的特征向量，也能出结果，代码量大概一百多行。代价是条件数差的时候精度掉得比 SVD 快，控制点分布不理想时差距明显。

2.3 第一次运行的检查清单

工程能编译通过不等于能算对，按顺序过一遍这几项，能省掉大半排查时间。

1. 平台选x64。x86 下默认的 double 对齐和内存布局一般没问题，但如果链接了 x64 的 OpenCV 库，编译期就会报错，别在这上面浪费时间。
2. 字符集选“使用 Unicode 字符集”。MFC 工程改成多字节字符集后，CString到std::string的隐式转换会失效，读文件路径那行会编译不过。
3. 控制点文件用绝对路径，或者把工作目录设成$(ProjectDir)..\data。VS 默认工作目录是工程目录，不是输出目录，相对路径找不到文件是新人最常撞的一堵墙。
4. 至少六个控制点，且不能在同一个平面上，不能在一条直线上。少于六个，方程数不够，SVD 出来的解没有物理意义。
5. 先跑控制台入口，看到残差表输出正常，再去点 MFC 界面的按钮。界面出错的信息量远小于控制台。

3. DLT 数学模型拆解：一个点两行方程，六个点凑够 11 个未知数

DLT 的出发点特别朴素：共线条件方程说，物点、像点、摄影中心三点共线。把这句话写成代数式，会得到一个分式方程；把分母乘过去，就变成线性的了。线性之后就能用最小二乘，不用给初值，不用迭代，这也是 DLT 在工程现场一直有人用的原因——它不挑输入，给点就能算。

3.1 共线条件方程怎么线性化

设投影矩阵P是 3x4，元素记为L1到L11（第十二个元素固定为 1，作为尺度归一）：

x = (L1·X + L2·Y + L3·Z + L4) / (L9·X + L10·Y + L11·Z + 1) y = (L5·X + L6·Y + L7·Z + L8) / (L9·X + L10·Y + L11·Z + 1)

交叉相乘，把未知数挪到左边、观测值挪到右边，一个控制点拆出两行：

L1·X + L2·Y + L3·Z + L4 - x·L9·X - x·L10·Y - x·L11·Z = x L5·X + L6·Y + L7·Z + L8 - y·L9·X - y·L10·Y - y·L11·Z = y

每个点两行，六个点十二行，十一个未知数，方程数刚好压过未知数，这就是“最少六个控制点”的来处。点多当然更好，最小二乘会自动把矛盾摊平。

需要注意的是，上面这两行里X、Y、Z和x、y的量级差得非常远——物方坐标可能是几十米，像点坐标是几千像素。直接把这些数塞进矩阵，A^T A的条件数会大到离谱，解出来的 L9、L10、L11 全是噪声。所以真正能用的实现里，构建设计矩阵之前必须先归一化，这一步放到第 4 章展开。

3.2 设计矩阵 A 的构建代码

把上面两行照着敲进循环里就行，注意列的对应关系别错位。第 0～3 列对应 L1～L4，第 4～7 列对应 L5～L8，第 8～10 列对应 L9～L11，第 11 列是常数项的系数：

// 构建设计矩阵：A 是 2n x 12，最后求出的是齐次解，再除以第 12 个分量 static void buildDesignMatrix(const std::vector<ControlPoint>& pts, const std::vector<cv::Point2d>& nImg, const std::vector<cv::Point3d>& nObj, cv::Mat& A) { const int n = static_cast<int>(pts.size()); A = cv::Mat::zeros(2 * n, 12, CV_64F); // 必须是 double，float 精度不够 for (int i = 0; i < n; ++i) { const double X = nObj[i].x, Y = nObj[i].y, Z = nObj[i].z; const double x = nImg[i].x, y = nImg[i].y; double* r1 = A.ptr<double>(2 * i); // 第一个方程：x 方向 r1[0] = X; r1[1] = Y; r1[2] = Z; r1[3] = 1.0; r1[8] = -x * X; r1[9] = -x * Y; r1[10] = -x * Z; r1[11] = -x; double* r2 = A.ptr<double>(2 * i + 1); // 第二个方程：y 方向 r2[4] = X; r2[5] = Y; r2[6] = Z; r2[7] = 1.0; r2[8] = -y * X; r2[9] = -y * Y; r2[10] = -y * Z; r2[11] = -y; } }

参数说明几句。nObj和nImg是归一化之后的坐标，不是原始坐标，这一点和很多教科书直接写原始坐标不同，但工程上必须这么干。r1[11] = -x这一项是把等号右边的观测值挪到左边，最后求齐次解A·L = 0的时候，解向量的第 12 个分量理论上等于 L 的尺度因子，除一下就把L12 = 1的约束还原了。

矩阵用CV_64F，不要图省内存用CV_32F。float 只有 7 位有效数字，而设计矩阵里x·X这一项可能是 1e7 量级，误差直接吃掉小数点后好几位，残差表看着还行，内参反解出来的焦距能差出几十像素。

3.3 为什么用 SVD 而不是正规方程求逆

教科书写 DLT 一般写L = (A^T A)^{-1} A^T b，看着干净，实际上手就翻车。A^T A的条件数是A的条件数的平方，A本身条件数 1e4，平方之后 1e8，double 也扛不住。更麻烦的是控制点共面的时候，A的秩直接掉到 8 或 9，A^T A奇异，求逆那行要么抛异常，要么返回一堆 inf。

cv::SVD::solveZ干的事情是：对A做奇异值分解，取最小奇异值对应的右奇异向量作为解。这个解就是min ‖A·L‖在‖L‖=1约束下的答案，秩亏的时候它依然能给出一个确定的结果，而且奇异值排序本身就是最好的诊断工具——第二小和最小奇异值的比值一旦小到 1e-3 量级，说明有约束方向没被控制点约束住，得回去看点位分布。

我一般在求解之后加一段打印，把最小的三个奇异值打出来。正常数据下三者差距明显，出现两个几乎相等的小奇异值，基本可以断定点共面或者近似共线，这时候不要急着信结果。

4. 关键实现：归一化、RQ 内参反解与重投影残差

这一章是真正能抄作业的部分。前面把方程列清楚了，落到代码上还有三件事：坐标归一化、反归一化、从投影矩阵里拆内参。三件事漏一件，结果就不对，而且不一定报错，属于那种「数出来了但是没法用」的错误。

4.1 坐标归一化：数值稳定的前提

归一化的思路是 Hartley 提的那套：把像点平移到质心、缩放到平均距离等于sqrt(2)；把物点平移到质心、缩放到平均距离等于sqrt(3)。缩放之后所有坐标都落在 1 附近，设计矩阵各列量级一致，条件数立刻降下来。

// 归一化：像点平移到质心、平均距离拉到 sqrt(2)；物点拉到 sqrt(3) static void normalizePoints(const std::vector<ControlPoint>& pts, cv::Mat& T, cv::Mat& U, std::vector<cv::Point2d>& nImg, std::vector<cv::Point3d>& nObj) { const double n = static_cast<double>(pts.size()); double mx = 0, my = 0, MX = 0, MY = 0, MZ = 0; for (const auto& p : pts) { mx += p.x; my += p.y; MX += p.X; MY += p.Y; MZ += p.Z; } mx /= n; my /= n; MX /= n; MY /= n; MZ /= n; double dImg = 0, dObj = 0; for (const auto& p : pts) { dImg += std::hypot(p.x - mx, p.y - my); dObj += std::sqrt((p.X-MX)*(p.X-MX) + (p.Y-MY)*(p.Y-MY) + (p.Z-MZ)*(p.Z-MZ)); } const double sImg = std::sqrt(2.0) / (dImg / n); // 平均距离 -> sqrt(2) const double sObj = std::sqrt(3.0) / (dObj / n); // 平均距离 -> sqrt(3) T = (cv::Mat_<double>(3,3) << sImg, 0, -sImg*mx, 0, sImg, -sImg*my, 0, 0, 1); U = (cv::Mat_<double>(4,4) << sObj, 0, 0, -sObj*MX, 0, sObj, 0, -sObj*MY, 0, 0, sObj, -sObj*MZ, 0, 0, 0, 1); nImg.clear(); nObj.clear(); for (const auto& p : pts) { nImg.emplace_back(sImg * (p.x - mx), sImg * (p.y - my)); nObj.emplace_back(sObj * (p.X - MX), sObj * (p.Y - MY), sObj * (p.Z - MZ)); } }

T是 3x3 的像方变换，U是 4x4 的物方变换。之所以物方用 4x4，是因为它作用在齐次坐标(X, Y, Z, 1)上，平移要占第四列。这两个矩阵后面还要用，必须原样留着，不能算完就丢。

有一点容易忽略：物方缩放系数取sqrt(3)而不是sqrt(2)，因为物方是三维的，平均距离的量纲和像方的二维不一样。这两个常数写死的目的是让归一化后的点坐标尺度统一，也不是什么玄学，换个常数结果一样能用，只是数值稳定性差一点。

4.2 SVD 求解与反归一化

归一化之后的设计矩阵交给 SVD，拿到解向量之后还要乘回去。反归一化的公式是P = T⁻¹ · P' · U，其中P'是归一化坐标系下的解：

bool solveDlt(const std::vector<ControlPoint>& pts, cv::Mat& P, double& rms) { if (pts.size() < 6) return false; // 三维 DLT 的硬性下限 cv::Mat T, U; std::vector<cv::Point2d> nImg; std::vector<cv::Point3d> nObj; normalizePoints(pts, T, U, nImg, nObj); cv::Mat A; buildDesignMatrix(pts, nImg, nObj, A); cv::SVD svd(A, cv::SVD::FULL_UV); cv::Mat L = svd.vt.row(11).t(); // 最小奇异值对应的右奇异向量 cv::Mat Pn = L.reshape(1, 3); // 排成 3x4，此时还在归一化坐标系 P = T.inv() * Pn * U; // 反归一化：3x3 * 3x4 * 4x4 = 3x4 P /= cv::norm(P.row(2)); // 尺度归一，第三行范数取 1 rms = computeRms(P, pts); // 用原始坐标算残差，别用归一化的 return true; }

几处要说明。svd.vt.row(11)取的是最后一行，对于 2n x 12 的矩阵，vt是 12x12，最后一行对应最小奇异值，这就是solveZ内部做的事情，写开更直观，也方便中间打断点看奇异值。L.reshape(1, 3)是按行重排成 3 行，因为 OpenCV 的 reshape 是行优先连续内存，12 个元素刚好排成 3x4。

尺度归一那行P /= cv::norm(P.row(2))不能省。SVD 给出的解满足‖L‖ = 1，但投影矩阵是齐次量，只确定到一个比例因子。归一化到第三行范数为 1，后面拆内参时K(3,3) = 1的假设才成立。

残差一定要用原始坐标算，用归一化坐标算出来的数字是缩放过的，看着漂亮但没有物理意义。这一点我在别人的代码里见过好几次，残差表显示 0.001，实际重投影差了快两个像素。

4.3 RQ 分解反解内参与外参

有了P，内参就是它左边 3x3 子矩阵的 RQ 分解。P = K·[R|t]，K是上三角，R是正交阵。RQ 分解用 Givens 旋转从下往上把下三角元素消成零，代码不长，但下标容易写错：

// 3x3 的 RQ 分解：M = K * R，K 上三角（内参），R 正交（旋转） static void rq3(const double* M, double* K, double* R) { for (int i = 0; i < 9; ++i) { K[i] = M[i]; R[i] = 0.0; } R[0] = R[4] = R[8] = 1.0; for (int i = 2; i >= 1; --i) { // 从最后一列往前消 const double a = K[i*3 + i]; // 对角元素 const double b = K[i*3 + (i-1)]; // 待消的下三角元素 const double r = std::hypot(a, b); if (r < 1e-15) continue; // 该列已退化，跳过 const double c = a / r, s = -b / r; for (int j = 0; j < 3; ++j) { // K 右乘 Givens，混合第 i-1、i 列 const double u = K[j*3 + (i-1)], v = K[j*3 + i]; K[j*3 + (i-1)] = c * u + s * v; K[j*3 + i] = -s * u + c * v; } for (int j = 0; j < 3; ++j) { // R 同步左乘转置，保持 M = K*R const double u = R[(i-1)*3 + j], v = R[i*3 + j]; R[(i-1)*3 + j] = c * u - s * v; R[i*3 + j] = s * u + c * v; } } // 尺度归一到 K(3,3)=1，再把前两行对角元调成正 const double scale = K[8]; if (std::abs(scale) > 1e-15) for (int i = 0; i < 9; ++i) K[i] /= scale; for (int i = 0; i < 2; ++i) { if (K[i*3 + i] < 0) { for (int j = 0; j < 3; ++j) { K[i*3 + j] = -K[i*3 + j]; R[i*3 + j] = -R[i*3 + j]; } } } }

拆完之后，K(0,0)是 x 方向等效焦距，K(1,1)是 y 方向等效焦距，K(0,2)、K(1,2)是主点，K(0,1)是倾斜因子。工程上倾斜因子应该接近 0，如果解出来是几百，说明像点坐标里有系统性误差——比如扫描件没校正、镜头畸变没处理，或者像点量测的时候用了不同的原点。

外参的求法是把P左乘K⁻¹，得到[R|t]，第四列就是t，摄影中心C = -Rᵀ·t。有个细节：算出来的R行列式如果是 -1，说明整体符号反了，把P乘 -1 重来一遍就行，别硬凑。

bool decomposeP(const cv::Mat& P, cv::Mat& K, cv::Mat& R, cv::Mat& t, cv::Point3d& C) { cv::Mat M = P(cv::Rect(0, 0, 3, 3)).clone(); cv::Mat p4 = P.col(3).clone(); double k[9], r[9]; rq3(M.ptr<double>(), k, r); K = cv::Mat(3, 3, CV_64F, k).clone(); R = cv::Mat(3, 3, CV_64F, r).clone(); if (cv::determinant(R) < 0) { // 反射解，翻转符号重来 K = -K; R = -R; p4 = -p4; } t = K.inv() * p4; // [R|t] = K^-1 * P C = -R.t() * t; // 摄影中心在世界坐标系下 return true; }

4.4 残差评定与输出格式

残差表是唯一能让你判断「这次算得对不对」的东西，一定要逐点输出，不能只给一个 RMS。某个点残差特别大，往往是那个点的物方坐标抄错了或者像点戳歪了，只看总量根本发现不了。

// 重投影残差：把物方点投回像平面，与观测值做差 double computeRms(const cv::Mat& P, const std::vector<ControlPoint>& pts) { double sum = 0.0; for (const auto& cp : pts) { const double X[4] = { cp.X, cp.Y, cp.Z, 1.0 }; double u = 0, v = 0, w = 0; for (int j = 0; j < 4; ++j) { u += P.at<double>(0, j) * X[j]; v += P.at<double>(1, j) * X[j]; w += P.at<double>(2, j) * X[j]; // 齐次分母，接近 0 说明点在相机背后 } if (std::abs(w) < 1e-12) continue; const double dx = u / w - cp.x; const double dy = v / w - cp.y; sum += dx * dx + dy * dy; } return std::sqrt(sum / std::max<size_t>(1, pts.size())); }

输出我习惯用固定宽度对齐，一行一个点，带上dx、dy和点位误差，最后一行给 RMS 和最大值。格式稳定了，才好拿 Excel 或者 Python 画图看残差分布有没有方向性——残差如果都朝同一个方向偏，那是模型问题（比如镜头畸变），不是量测问题。

5. 避坑与排查：控制点共面、单位打架和条件数爆表的现场

DLT 的代码量不大，翻车基本都翻在数据和数值上。下面几条是按「现象 → 原因 → 解决」整理的现场记录，遇到问题时对着看，比漫无目的地改代码快得多。

5.1 数据层面的三个坑

现象：六个控制点算出来残差几乎为零，但拿去投影一个不在控制点范围内的点，偏了几十个像素。

原因：六个控制点全部布在同一面墙上，Z坐标完全相同。此时设计矩阵第 2 列和第 3 列线性相关，秩掉到 10 甚至更低，SVD 给出的解在L9、L10、L11上完全是噪声，控制点范围内因为方程刚好被满足所以残差漂亮，一旦外推立刻露馅。

解决：控制点必须在深度方向上有变化，理想情况是构成一个立方体的角点，退一步也要保证至少两个不同深度的平面。实在只能在平面上做，那就别用三维 DLT，改用平面单应矩阵 `

本文还有配套的精品资源，点击获取