简介面向摄影测量与计算机视觉学习者的DLT直接线性变换VC实现代码包重点解决由二维像点坐标反算三维物方坐标的参数求解问题适合已具备测量平差与C基础、希望动手复现算法的中高级读者。压缩包共38个文件、约323KB以cpp与h源文件为核心配合dsw、dsp工程配置与rc、ico界面资源构成一套可直接编译的MFC对话框程序另附4个txt控制点坐标与读数文件及ReadMe说明便于对照输入输出。该资源已有622人学习下载。代码围绕至少5对同名点构建线性方程组逐步估计投影中心位置与焦距、主点等内参并给出最小二乘线性化与结果输出的完整流程注释覆盖关键函数读者可借此理解DLT的方程组织方式、数据组织格式与调试思路并迁移到航测、三维重建等实际项目中。1. 六个控制点换一组相机参数DLT 在 VC 里到底解决什么问题给你一台普通相机、六个以上已知坐标的控制点再在照片上点出它们的像素位置能不能反推出相机在哪、朝哪、焦距多少顺带把未知点的三维坐标算出来直接线性变换DLT就是干这个的把非线性的共线条件方程乘开移项硬掰成一组线性方程不用给初值解 11 个 L 系数就完事。这套 VC 程序代码围绕这条链路组织读控制点文件、组误差方程、解法方程、还原内方位元素和摄站坐标、输出残差。它适合摄影测量课程设计、近景测量小工具的原型验证以及手里有 MFC 底子、想找一个能直接编译的完整 DLT 例子的人。2. 工程骨架拆解MFC 界面与解算内核怎么分家一套能长期改的 DLT 程序最怕把矩阵运算、文件读写、界面刷新全堆在按钮响应函数里。这套 VC 工程的思路是把「解算内核」抽成不依赖 MFC 的纯 C 类对话框只负责三件事选文件、调内核、把结果填到控件。这样后面想改成控制台批处理、或者封装成 DLL 给别的系统调用都不用重写算法。2.1 目录结构与文件职责工程按典型的 MFC 对话框程序组织源码目录大致是这样DLT/ ├─ DLT.sln ├─ DLT/ │ ├─ DLT.vcxproj │ ├─ DLT.h / DLT.cpp // 应用入口注册主对话框 │ ├─ DLTDlg.h / DLTDlg.cpp // 界面逻辑读文件、触发解算、显示结果 │ ├─ DltSolver.h / .cpp // 解算内核组方程、求解、还原相机参数 │ ├─ Matrix.h / .cpp // 矩阵工具转置、乘、高斯消元、SVD 备选 │ ├─ ControlPoint.h // 控制点与结果结构体 │ ├─ stdafx.h / .cpp // 预编译头 │ └─ res/ // 图标、对话框资源 └─ data/ └─ control_points.txt // 示例控制点文件职责二次开发时的改动频率DLTDlg.cpp界面事件、文件选择、结果展示中换展示方式才动DltSolver.cppDLT 解算主流程高换算法、加粗差检测都在这Matrix.cpp线性代数底座低除非换 EigenControlPoint.h数据结构中加字段要同步读写control_points.txt输入样例高换项目就换数据这张表的价值在于告诉你改哪里安全。界面和算法分离之后DltSolver 里不出现CWnd、CString这类 MFC 类型全部用std::vector、std::string和裸double。我第一次接手类似代码时解算函数里混着UpdateData(TRUE)想把算法搬到控制台测试光剥 MFC 依赖就花了一下午这种坑完全可以提前避开。2.2 控制点结构与输入文件格式DLT 的输入非常朴素每个控制点需要物方三维坐标和对应的像点二维坐标。结构体定义如下字段顺序和文本文件列顺序保持一致// ControlPoint.h struct ControlPoint { int id; // 点号仅用于输出定位不参与解算 double X, Y, Z; // 物方空间坐标单位与工程约定一致米或毫米 double x, y; // 像点坐标单位像素原点在图像左上角 }; struct DltResult { double L[11]; // DLT 的 11 个系数顺序 L1~L11 double x0, y0; // 像主点像素 double f, fx, fy; // 焦距像素fx/fy 为 x、y 方向分别反算的值 double R[3][3]; // 旋转矩阵行优先 double Xs, Ys, Zs; // 摄站坐标与物方坐标同单位 double sigma0; // 单位权中误差像素 std::vectordouble vx, vy; // 各点 x、y 方向残差 };物方坐标和像点坐标的单位必须各自统一这一点后面避坑章还会展开。文本文件采用空白分隔一行一个点1 0.0000 0.0000 0.0000 1024.5 768.2 2 2.0000 0.0000 0.0500 1180.3 770.1 3 2.0000 1.5000 0.0200 1176.8 512.4 4 0.0000 1.5000 0.0000 1020.7 509.9 5 0.0000 0.0000 1.2000 1026.1 402.3 6 2.0000 1.5000 1.2500 1182.0 398.7 7 1.0000 0.7500 0.6000 1101.2 601.5读文件的实现不要用while(!feof(fp))这是老 VC 代码里最常见的翻车点最后一行会被重复读一次。用fscanf的返回值判断更稳// 读取控制点文件返回 false 表示文件打不开或点数不足 bool LoadControlPoints(const char* path, std::vectorControlPoint pts) { FILE* fp fopen(path, r); if (!fp) return false; pts.clear(); ControlPoint p; // fscanf 返回成功匹配的字段数6 个字段全部读到才算一条有效记录 while (fscanf(fp, %d %lf %lf %lf %lf %lf, p.id, p.X, p.Y, p.Z, p.x, p.y) 6) { pts.push_back(p); } fclose(fp); return pts.size() 6; // 11 个未知数至少 6 个点才可能满秩 }fscanf的格式串里%d对应点号五个%lf依次是 X、Y、Z、x、y。返回值等于 6 说明整行读全了遇到空行或注释行会提前退出循环。这里没有做注释行过滤实际项目里如果要在文件里写#开头的说明得先fgets读整行再sscanf直接fscanf会在注释行卡住。2.3 界面层与解算层的调用关系对话框上一般放四个东西一个「打开文件」按钮、一个控制点列表控件、一个「解算」按钮、一个多行编辑框显示结果。点击解算时不做任何计算只把数据往内核递// DLTDlg.cpp 中的解算按钮响应节选 void CDLTDlg::OnBnClickedSolve() { std::vectorControlPoint pts; if (!LoadControlPoints(m_strPath.c_str(), pts)) { AfxMessageBox(_T(控制点读取失败或少于 6 个)); return; } DltSolver solver; DltResult result; if (!solver.Solve(pts, result)) { AfxMessageBox(_T(法方程奇异检查控制点是否共面)); return; } CString msg; msg.Format(_T(x0%.2f y0%.2f f%.2f\r\nXs%.4f Ys%.4f Zs%.4f\r\nm0%.3f pixel), result.x0, result.y0, result.f, result.Xs, result.Ys, result.Zs, result.sigma0); SetDlgItemText(IDC_EDIT_RESULT, msg); FillResidualList(pts, result); // 残差填到列表控件 }这段代码本身没有算法但它决定了工程的扩展性。DltSolver::Solve的签名只依赖标准容器返回值用bool表达成功失败错误信息由界面层负责翻译成人话。后面想加批量处理只要在控制台程序里复用LoadControlPoints和DltSolver即可界面代码一行不用碰。3. 算法主线从共线条件方程到 11 个 L 系数的完整推导与代码DLT 的核心不是「解方程」三个字而是搞清楚 11 个未知数从哪来、每个控制点贡献几行、解完之后怎么把 L 系数翻译回相机参数。这一章按数据流顺序拆先写误差方程再解法方程再做参数还原最后算残差。每一步都对应工程里的一个函数。3.1 共线条件方程的线性化形式摄影测量的共线条件方程描述的是「物点、摄站点、像点三点共线」x - x0 -f · [a1(X-Xs) b1(Y-Ys) c1(Z-Zs)] / [a3(X-Xs) b3(Y-Ys) c3(Z-Zs)] y - y0 -f · [a2(X-Xs) b2(Y-Ys) c2(Z-Zs)] / [a3(X-Xs) b3(Y-Ys) c3(Z-Zs)]这里 a、b、c 是旋转矩阵的九个元素Xs、Ys、Zs 是摄站坐标x0、y0、f 是内方位元素。直接解这组非线性方程需要初值初值给不好就发散。DLT 的做法是把分子分母整体除一下令L1 (x0·a3 - f·a1)/D L2 (x0·b3 - f·b1)/D L3 (x0·c3 - f·c1)/D L4 -(L1·Xs L2·Ys L3·Zs) ... 以此类推整理后得到只含 L 系数的线性分式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)分母最后那个 1 是把比例因子吸收进 L 系数的结果等价于强制 L12 1。这个归一化很关键不固定 L12方程有无数个成比例的解法方程永远秩亏。交叉相乘再移项每个控制点得到两个线性方程X·L1 Y·L2 Z·L3 L4 - x·X·L9 - x·Y·L10 - x·Z·L11 x X·L5 Y·L6 Z·L7 L8 - y·X·L9 - y·Y·L10 - y·Z·L11 y到这里非线性问题就变成了标准的A·L B最小二乘问题6 个控制点给 12 行方程、11 个未知数刚好超定一行实际用 8 到 12 个点更稳。3.2 误差方程矩阵的组装组装误差方程的函数是整条链路的第一环把每个控制点拆成两行写进二维数组。注意 x 方向方程只碰 L1~L4 和 L9~L11y 方向方程只碰 L5~L8 和 L9~L11这个稀疏规律别写错列号// 由 n 个控制点组误差方程 A·L BA 为 (2n)×11 bool BuildDesignMatrix(const std::vectorControlPoint pts, std::vectorstd::vectordouble A, std::vectordouble B) { const int n (int)pts.size(); if (n 6) return false; A.assign(2 * n, std::vectordouble(11, 0.0)); B.assign(2 * n, 0.0); for (int i 0; i n; i) { const ControlPoint p pts[i]; const int r 2 * i; // x 方程行号 A[r][0] p.X; A[r][1] p.Y; A[r][2] p.Z; A[r][3] 1.0; A[r][8] -p.x * p.X; A[r][9] -p.x * p.Y; A[r][10] -p.x * p.Z; B[r] p.x; A[r1][4] p.X; A[r1][5] p.Y; A[r1][6] p.Z; A[r1][7] 1.0; A[r1][8] -p.y * p.X; A[r1][9] -p.y * p.Y; A[r1][10] -p.y * p.Z; B[r1] p.y; } return true; }参数说明pts是读入的控制点数组A是输出设计矩阵行数 2n、列数 11B是观测向量。第 8、9、10 列下标从 0 起在 x、y 两个方程里都要填符号是负的这是交叉相乘移项的结果写反了焦距会算成负数。这一步纯属矩阵填充没有数值风险但列号错了后面全盘皆输我一般会拿第一个点手算一遍两行方程对一对。3.3 法方程与求解为什么优先用归一化加列主元有了 A 和 B最直接的做法是组法方程N AᵀA、U AᵀB解N·L U。法方程的好处是维度固定 11×11坏处是条件数平方放大。控制点坐标动不动几千像素直接组出来的 N 矩阵元素量级能差十几个数量级高斯消元选主元稍微马虎一点就丢精度。工程里先用列主元高斯-约当消元写完再加一步坐标归一化。归一化的思路是把所有控制点的像点坐标平移到质心、物方坐标缩放到 [-1, 1] 附近解完 L 系数再反变换回来。这一步在 DLT 里效果立竿见影因为它把法方程的条件数压下来一大截。// 解法方程 (AᵀA)·L AᵀB列主元高斯-约当消元 bool SolveNormalEquation(const std::vectorstd::vectordouble A, const std::vectordouble B, std::vectordouble L) { const int m 11; std::vectorstd::vectordouble N(m, std::vectordouble(m 1, 0.0)); // 组法方程增广矩阵最后一列放 AᵀB for (int i 0; i m; i) { for (int j 0; j m; j) { double s 0.0; for (size_t k 0; k A.size(); k) s A[k][i] * A[k][j]; N[i][j] s; } double u 0.0; for (size_t k 0; k A.size(); k) u A[k][i] * B[k]; N[i][m] u; } // 列主元消元 for (int c 0; c m; c) { int pivot c; for (int r c 1; r m; r) if (fabs(N[r][c]) fabs(N[pivot][c])) pivot r; // 主元过小控制点共面或分布退化的典型信号 if (fabs(N[pivot][c]) 1e-12) return false; std::swap(N[c], N[pivot]); double d N[c][c]; for (int j c; j m; j) N[c][j] / d; for (int r 0; r m; r) { if (r c) continue; double f N[r][c]; for (int j c; j m; j) N[r][j] - f * N[c][j]; } } L.resize(m); for (int i 0; i m; i) L[i] N[i][m]; return true; }参数说明N是 11×12 的增广矩阵最后一列是右端项1e-12这个阈值是经验值物方坐标单位是毫米时够用单位换成米要相应调小。如果这里返回false先别怀疑代码八成是控制点布设的问题第五章会专门讲怎么排查。想更稳可以把这部分换成奇异值分解A直接做 SVD取最小奇异值对应的右奇异向量连法方程都不用组代价是多几十行矩阵代码。3.4 从 L 系数还原内方位元素与旋转矩阵解出 11 个 L 系数只是中间结果工程上真正要的是 x0、y0、f、旋转矩阵和摄站坐标。还原公式推导起来不复杂但符号约定容易错尤其是图像坐标系 y 轴朝下这一条// 由 L 系数还原内方位元素和旋转矩阵 bool RecoverCamera(const std::vectordouble L, DltResult out) { const double L1L[0], L2L[1], L3L[2], L4L[3]; const double L5L[4], L6L[5], L7L[6], L8L[7]; const double L9L[8], L10L[9], L11L[10]; double w L9*L9 L10*L10 L11*L11; if (w 1e-15) return false; double rho 1.0 / sqrt(w); // 像主点分母二次型的投影 out.x0 -(L1*L9 L2*L10 L3*L11) / w; out.y0 -(L5*L9 L6*L10 L7*L11) / w; // 焦距两个方向分别反算正常情况应接近相等 double fx2 (L1*L1 L2*L2 L3*L3) / w - out.x0*out.x0; double fy2 (L5*L5 L6*L6 L7*L7) / w - out.y0*out.y0; if (fx2 0 || fy2 0) return false; out.fx sqrt(fx2); out.fy sqrt(fy2); out.f 0.5 * (out.fx out.fy); // 旋转矩阵第三行直接由 L9~L11 归一化得到 out.R[2][0] rho * L9; out.R[2][1] rho * L10; out.R[2][2] rho * L11; // 前两行由内参反算再做叉乘正交化 double r1[3], r2[3], r3[3] { out.R[2][0], out.R[2][1], out.R[2][2] }; r1[0] (L1 - out.x0*L9) / out.fx; r1[1] (L2 - out.x0*L10) / out.fx; r1[2] (L3 - out.x0*L11) / out.fx; r2[0] (L5 - out.y0*L9) / out.fy; r2[1] (L6 - out.y0*L10) / out.fy; r2[2] (L7 - out.y0*L11) / out.fy; // 叉乘重建第二行再重建第一行压制测量噪声 r2[0] r3[1]*r1[2] - r3[2]*r1[1]; r2[1] r3[2]*r1[0] - r3[0]*r1[2]; r2[2] r3[0]*r1[1] - r3[1]*r1[0]; double n2 sqrt(r2[0]*r2[0] r2[1]*r2[1] r2[2]*r2[2]); if (n2 1e-12) return false; for (int i 0; i 3; i) r2[i] / n2; r1[0] r2[1]*r3[2] - r2[2]*r3[1]; r1[1] r2[2]*r3[0] - r2[0]*r3[2]; r1[2] r2[0]*r3[1] - r2[1]*r3[0]; for (int i 0; i 3; i) { out.R[0][i] r1[i]; out.R[1][i] r2[i]; } return true; }参数说明L是上一步解出的 11 个系数out同时接收内参和旋转矩阵。rho是比例归一化因子作用是把第三行旋转向量变回单位向量。fx2、fy2理论上应该相等实际会有零点几像素的差异取平均得到f。符号约定方面如果输入的像点 y 轴朝上而不是朝下x0、y0 会整体反号写完之后务必用反投影残差验证一遍别死记公式。摄站坐标的还原更直接。把 L1~L3、L5~L7、L9~L11 看成 3×3 矩阵 ML4、L8、1 看成第四列 p4相机中心就是C -M⁻¹·p4// 解 3×3 线性方程组 M·C -p4得到摄站坐标 bool RecoverStation(const std::vectordouble L, DltResult out) { double M[3][3] { { L[0], L[1], L[2] }, { L[4], L[5], L[6] }, { L[8], L[9], L[10] } }; double b[3] { -L[3], -L[7], -1.0 }; double det M[0][0]*(M[1][1]*M[2][2] - M[1][2]*M[2][1]) - M[0][1]*(M[1][0]*M[2][2] - M[1][2]*M[2][0]) M[0][2]*(M[1][0]*M[2][1] - M[1][1]*M[2][0]); if (fabs(det) 1e-15) return false; // 克拉默法则三个分量分别替换一列求行列式 double C[3]; for (int k 0; k 3; k) { double T[3][3]; for (int i 0; i 3; i) for (int j 0; j 3; j) T[i][j] (j k) ? b[i] : M[i][j]; C[k] (T[0][0]*(T[1][1]*T[2][2] - T[1][2]*T[2][1]) - T[0][1]*(T[1][0]*T[2][2] - T[1][2]*T[2][0]) T[0][2]*(T[1][0]*T[2][1] - T[1][1]*T[2][0])) / det; } out.Xs C[0]; out.Ys C[1]; out.Zs C[2]; return true; }克拉默法则在这里够用3×3 求逆没必要上通用矩阵库。det接近 0 说明三个方向的 L 系数线性相关通常还是控制点布设问题。算出来的摄站坐标单位跟物方坐标一致如果物方用的是米结果就是米别自己再乘 1000。3.5 精度评定残差、单位权中误差与粗差定位解算不输出残差等于耍流氓。DLT 的多余观测数是2n - 11用残差平方和除以多余观测数开根号得到单位权中误差单位是像素。反投影代码同时把每个点的残差存下来列表控件里哪一行残差大哪一行就是可疑粗差// 反投影计算残差与单位权中误差 void ComputeResidual(const std::vectorControlPoint pts, const std::vectordouble L, DltResult out) { const int n (int)pts.size(); out.vx.assign(n, 0.0); out.vy.assign(n, 0.0); double sum 0.0; for (int i 0; i n; i) { const ControlPoint p pts[i]; double den L[8]*p.X L[9]*p.Y L[10]*p.Z 1.0; if (fabs(den) 1e-15) { out.sigma0 1e9; return; } double px (L[0]*p.X L[1]*p.Y L[2]*p.Z L[3]) / den; double py (L[4]*p.X L[5]*p.Y L[6]*p.Z L[7]) / den; out.vx[i] p.x - px; out.vy[i] p.y - py; sum out.vx[i]*out.vx[i] out.vy[i]*out.vy[i]; } // 多余观测数 2n - 11n6 时只有 1精度指标参考意义有限 int redundant 2 * n - 11; out.sigma0 (redundant 0) ? sqrt(sum / redundant) : 0.0; }参数说明vx、vy分别是 x、y 方向残差sigma0是单位权中误差。n6时多余观测只有 1sigma0数值上好看但统计意义很弱实际项目一般建议至少 8 个控制点。判断标准方面像点坐标人工量测误差通常在 0.5 到 2 个像素sigma0明显超过这个量级要么是粗差要么是控制点坐标本身有问题别急着调算法参数。4. 编译跑通环境配置、控制点文件与结果自检代码看完下一步是让它在这台机器上跑起来。VC 工程的坑集中在项目属性上尤其是 MFC 链接方式和运行库匹配配错一个就是一堆链接错误。这一章从新建配置讲到结果验证按顺序走一遍。4.1 环境准备与项目属性设置打开解决方案之前先确认开发环境装了 MFC 组件。VC 安装器里 MFC 是可选组件默认不一定勾选缺了会在编译时报「无法打开 afxwin.h」。装好之后按下面的属性表逐项核对属性项推荐值说明平台x64老工程可能只有 Win32新建一个 x64 配置MFC 使用方式在静态库中使用 MFC目标机不用装运行库拷贝 exe 就能跑运行库多线程 (/MT)与静态 MFC 配套选错会报 LNK2005 重定义字符集使用多字节字符集老代码里char路径处理不用改SDL 检查否否则fopen、fscanf会被当成不安全函数报错C 语言标准默认C14 即可代码只用了容器和基础语法不用开新标准「SDL 检查」这一项单独说一句它把一批 C 运行时函数标记为不安全fopen、sprintf都会触发 C4996 警告工程如果把警告当错误就直接编不过。关掉是最省事的做法想保留检查就把这些函数换成_s后缀版本但代码里会多出一堆参数。编译时如果报LNK2019: 无法解析的外部符号并且符号名里带Afx基本就是 MFC 链接方式配错了报LNK2005: 已经在 msvcrt.lib 中定义就是运行库和 MFC 方式不匹配。这两类错误不要去看代码回去改属性表。4.2 控制点布设与输入数据准备程序能跑不等于结果可信输入数据的质量决定了解算上限。控制点布设有几条硬经验点数不少于 8 个不要全部落在同一个平面上在图像上要铺开四个角各放一个、中间再补几个物方坐标的测量精度要明显高于像点量测精度否则控制点误差会直接吞掉解算精度。用一份仿真数据做首轮验证最稳妥先假定一组相机参数x0512、y0384、f1500 像素、摄站在物方上方若干米用共线方程正算出每个控制点的像点坐标再加一点随机噪声生成输入文件。这样解算结果应该回到假定值附近闭环验证比拿真实数据瞎猜强得多。1 0.0000 0.0000 0.0000 1024.5 768.2 2 2.0000 0.0000 0.0500 1180.3 770.1 3 2.0000 1.5000 0.0200 1176.8 512.4 4 0.0000 1.5000 0.0000 1020.7 509.9 5 0.0000 0.0000 1.2000 1026.1 402.3 6 2.0000 1.5000 1.2500 1182.0 398.7 7 1.0000 0.7500 0.6000 1101.2 601.5 8 0.5000 1.2000 1.1000 1060.4 455.8这组点的物方范围 2 米 × 1.5 米 × 1.25 米深度方向有变化不是共面布置。像点坐标分布在图像四角和中上部覆盖范围够。实际项目里如果只能在一个平面上布点DLT 会退化成平面单应11 个参数解不出来得换 8 参数模型这一点在避坑章展开。4.3 运行流程与结果自检完整跑一遍的步骤用 VC 打开DLT.sln切到 x64 配置生成解决方案把data/control_points.txt拷到 exe 所在目录或者点击界面上的打开按钮选绝对路径确认列表控件里读到了 8 行控制点点号、坐标显示正常点「解算」看结果框输出对照残差列表检查有没有单点残差明显偏大。正常输出大概长这样内方位元素: x0 511.87 y0 384.12 f 1499.63 p a hrefhttps://download.csdn.net/download/wlfeel1/2932412 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p