手把手教你实现条纹结构光三维重建(3)——相机投影仪标定

我们都知道,投影仪其实就是个反向相机,如果我们了解双目标定的原理,那么相机和投影仪的标定就不难,关键是我们怎么得到投影仪在图像特征点(比如棋盘格角点)上的像素位置。

投影仪也类似于一个cmos,图像有像素位置(u,v),那么通过我们上一讲的条纹解码,给图像添加水平方向和垂直方向的投影,就可以通过解码,得到图像对应的投影相位值,此相位值就是投影的像素坐标(xp,yp)。如下图所示,具体的原理可以参考论文《Calibration of fringe projection profilometry A comparative review》。

这里,为了方便大家理解,可以参考中文文献《基于数字光栅投影的结构光三维测量技术与系统研究》 ,我将重点截图如下:

第四点,绝对相位映射到DMD图像坐标,如下所示:

say easy,show me code。。。

代码如下:

#include <iostream>
#include <vector>
#include <opencv2/opencv.hpp>
//#define PROJECTOR_WIDTH_1080P    1920          //这里定义的是1080P投影仪的分辨率
//#define PROJECTOR_HEIGHT_1080P   1080

#define PROJECTOR_WIDTH_1080P    912          //这里定义的是1080P投影仪的分辨率
#define PROJECTOR_HEIGHT_1080P   1140

#define PROJECTOR_WIDTH_720P    1280          //这里定义的是720P投影仪的分辨率
#define PROJECTOR_HEIGHT_720P   720

#define PROJECTOR_WIDTH_480P    854          //这里定义的是480P投影仪的分辨率
#define PROJECTOR_HEIGHT_480P   480

#define PROJECTOR_WIDTH_4500    912          //这里定义的是TI的4500投影仪的分辨率
#define PROJECTOR_HEIGHT_4500   1140
#define PI 3.1415926

std::vector<float> h_freq_array, v_freq_array;

// Absolute phase from 4 frames
cv::Mat get_phase4(const cv::Mat I1, const cv::Mat I2, const cv::Mat I3, const cv::Mat I4) {
	cv::Mat_<float> I1_(I1);
	cv::Mat_<float> I2_(I2);
	cv::Mat_<float> I3_(I3);
	cv::Mat_<float> I4_(I4);

	//获取wrap相位
	int m_nHeight = I1.rows;
	int m_nWidth = I1.cols;
	cv::Mat phase = cv::Mat::zeros(m_nHeight, m_nWidth, CV_32FC1);
	//#pragma omp parallel for
	for (int i = 0; i < m_nHeight; i++)
	{
		for (int j = 0; j < m_nWidth; j++)
		{
			int a1 = I1_.at<float>(i, j);
			int a2 = I2_.at<float>(i, j);
			int a3 = I3_.at<float>(i, j);
			int a4 = I4_.at<float>(i, j);
			phase.at<float>(i, j) = (float)atan2((a2 - a4), (a1 - a3));
			if (phase.at<float>(i, j) < 0)
				phase.at<float>(i, j) += (2 * PI);
		}
	}
	return phase;
}
cv::Mat unwrap_with_cue(const cv::Mat up, const cv::Mat upCue, float nPhase)
{
	// Determine number of jumps
	cv::Mat P = ((upCue)*nPhase - up) / (2 * PI);
	cv::Mat tmp(P.rows, P.cols, CV_32FC1);
	for (int i = 0; i < up.rows; i++) {
		for (int j = 0; j < up.cols; j++) {
			tmp.at<float>(i, j) = round(P.at<float>(i, j));
		}
	}
	// Add to phase
	cv::Mat upUnwrapped = up + tmp * 2 * PI;
	// Scale to range [0; 2pi]
	upUnwrapped *= 1.0 / nPhase;
	return upUnwrapped;
}

cv::Mat decode_pattern(const std::vector<cv::Mat>& encode_images, const std::vector<float>& phases, const int projector_lens)
{
	//前面四组图案最低频率的编码(频率为1),所以不需要进行相位展开
	std::vector<cv::Mat> frames_low_freq(encode_images.begin(), encode_images.begin() + 4);
	cv::Mat upCue = get_phase4(frames_low_freq[0], frames_low_freq[1], frames_low_freq[2], frames_low_freq[3]);
	for (int index = 1; index < phases.size(); index++)  //两两求解双频
	{
		std::vector<cv::Mat> frames_high_freq(encode_images.begin() + 4 * (index), encode_images.begin() + 4 * (index + 1));
		cv::Mat unPhase = get_phase4(frames_high_freq[0], frames_high_freq[1], frames_high_freq[2], frames_high_freq[3]);
		upCue = unwrap_with_cue(unPhase, upCue, phases[index]);
	}
	cv::Mat decode_phase_img = projector_lens * ((upCue) / (2 * PI));
	return decode_phase_img;
}

void generate_freqs(std::vector <float>& freq_array, int length, int min_T)
{
	freq_array[4] = (double)length / min_T;     //我们需要生成五个频率,第五个频率为[投影宽度/周期]
	double x = sqrtf(sqrtf(freq_array[4]));    //第二个频率定义为第五个频率的开四次根号
	freq_array[3] = x * x * x; //第四个频率  
	freq_array[2] = x * x;     //第三个频率
	freq_array[1] = x;         //第二个频率
	freq_array[0] = 1;         //第一个频率
}

void init()
{
	v_freq_array.resize(5);
	generate_freqs(v_freq_array, PROJECTOR_HEIGHT_4500, 10);        //图像垂直方向——横条纹频率
	h_freq_array.resize(5);
	generate_freqs(h_freq_array, PROJECTOR_WIDTH_4500, 10);         //图像水平方向——竖条纹频率
}

void save_calibrate_xml(std::string filename,const cv::Mat& Kc, const cv::Mat& kc, double cam_error,
	const cv::Mat& Kp, const cv::Mat& kp, double proj_error,
	const cv::Mat& Rpc, const cv::Mat& Tpc, double stereo_error)
{
	cv::FileStorage fs(filename, cv::FileStorage::WRITE);
	if (!fs.isOpened())
		return;

	fs << "Kc" << cv::Mat(Kc) << "kc" << cv::Mat(kc) 
		<< "Kp" << cv::Mat(Kp) << "kp" << cv::Mat(kp) 
		<< "Rpc" << cv::Mat(Rpc) << "Tpc" << cv::Mat(Tpc) << "cam_error" << cam_error << "proj_error"
		<< proj_error << "stereo_error" << stereo_error;
	fs.release();
}

bool calibrate_Impl(const std::vector<cv::Mat>& chess_imgs, const std::vector<cv::Mat>& vps, const std::vector<cv::Mat>& ups)
{
	unsigned int img_height = chess_imgs[0].rows;
	unsigned int img_width = chess_imgs[0].cols;
	cv::Size pattern_size(8, 11);
	float checkerSize = 10;
	std::vector<cv::Point3f> Qi;
	for (int h = 0; h < pattern_size.height; h++)
		for (int w = 0; w < pattern_size.width; w++)
			Qi.push_back(cv::Point3f(checkerSize * w, checkerSize * h, 0.0));
	
	std::vector<std::vector<cv::Point2f> > qc, qp;
	std::vector<std::vector<cv::Point3f> > object_points;
	int nFrameSeq = chess_imgs.size();
	for (unsigned int i = 0; i < nFrameSeq; i++)
	{
		std::vector<cv::Point2f> qci;
		bool success =
			cv::findChessboardCorners(chess_imgs[i], pattern_size, qci, cv::CALIB_CB_ADAPTIVE_THRESH);
		if (!success)
		{
			std::cout << "图像 " << i << " 不能找到角点\n";
			continue;
		}
		else {
			// Refine corner locations
			cv::cornerSubPix(chess_imgs[i], qci, cv::Size(5, 5), cv::Size(1, 1),
				cv::TermCriteria(CV_TERMCRIT_EPS + CV_TERMCRIT_ITER, 30, 0.1));
		}
		// Draw colored chessboard
		cv::Mat chess_imgs_color;
		cv::cvtColor(chess_imgs[i], chess_imgs_color, cv::COLOR_GRAY2RGB);
		cv::drawChessboardCorners(chess_imgs_color, pattern_size, qci, success);
#if 1
		cv::resize(chess_imgs_color, chess_imgs_color, cv::Size(chess_imgs_color.cols * 0.5, chess_imgs_color.rows * 0.5));
		cv::imshow("chess_imgs_color", chess_imgs_color);
		cv::waitKey(0);
#endif
		// Vectors of accepted points for current view
		std::vector<cv::Point2f> qpi_a;
		std::vector<cv::Point2f> qci_a;
		std::vector<cv::Point3f> Qi_a;

		// Loop through checkerboard corners
		for (unsigned int j = 0; j < qci.size(); j++) {
			const cv::Point2f& qcij = qci[j];

			// Collect neighbor points
			const unsigned int WINDOW_SIZE = 10;
			std::vector<cv::Point2f> N_qcij, N_qpij;

			// avoid going out of bounds
			unsigned int starth = std::max(int(qcij.y + 0.5) - WINDOW_SIZE, 0u);
			unsigned int stoph = std::min(int(qcij.y + 0.5) + WINDOW_SIZE, img_height - 1);
			unsigned int startw = std::max(int(qcij.x + 0.5) - WINDOW_SIZE, 0u);
			unsigned int stopw = std::min(int(qcij.x + 0.5) + WINDOW_SIZE, img_width - 1);

			for (unsigned int h = starth; h <= stoph; h++) 
			{
				for (unsigned int w = startw; w <= stopw; w++) 
				{
						N_qcij.push_back(cv::Point2f(w, h));

						float upijwh = ups[i].at<float>(h, w);
						float vpijwh = vps[i].at<float>(h, w);
						N_qpij.push_back(cv::Point2f(upijwh, vpijwh));
				}
			}
			// std::cout << i << " findHomography " << N_qcij.size() << " " << N_qpij.size() <<
			// std::endl;
			// if enough valid points to build homography
			if (N_qpij.size() >= 50) {
				//                    std::cout << i << " findHomography" << std::endl;
				// translate qcij into qpij using local homography
				cv::Mat H = cv::findHomography(N_qcij, N_qpij, cv::RANSAC);
				/*std::cout << "H:\n" << H << endl;*/
				if (!H.empty()) {
					cv::Point3d Q =
						cv::Point3d(cv::Mat(H * cv::Mat(cv::Point3d(qcij.x, qcij.y, 1.0))));
					cv::Point2f qpij = cv::Point2f(Q.x / Q.z, Q.y / Q.z);

					qpi_a.push_back(qpij);
					qci_a.push_back(qci[j]);
					Qi_a.push_back(Qi[j]);
				}
			}
		}

		if (!Qi_a.empty()) {
			// Store projector corner coordinates
			qp.push_back(qpi_a);

			// Store camera corner coordinates
			qc.push_back(qci_a);

			// Store world corner coordinates
			object_points.push_back(Qi_a);
		}
	}
	if (object_points.size() <= 4) {
		std::cerr << "没有足够的标定数据!" << std::endl;
		return false;
	}
	// 标定相机
	cv::Mat Kc, kc;
	std::vector<cv::Mat> cam_rvecs, cam_tvecs;
	cv::Size frameSize(img_width, img_height);

	int cal_flags = 0
		+ cv::CALIB_FIX_K4
		+ cv::CALIB_FIX_K5;

	double cam_error = cv::calibrateCamera(
		object_points, qc, frameSize, Kc, kc, cam_rvecs, cam_tvecs, cal_flags,
		cv::TermCriteria(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, 50, DBL_EPSILON));

	// 标定投影仪
	cv::Mat Kp, kp;
	std::vector<cv::Mat> proj_rvecs, proj_tvecs;
	cv::Size screenSize(PROJECTOR_WIDTH_4500, PROJECTOR_HEIGHT_4500);
	double proj_error = cv::calibrateCamera(
		object_points, qp, screenSize, Kp, kp, proj_rvecs, proj_tvecs, cal_flags,
		cv::TermCriteria(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, 50, DBL_EPSILON));

	// 双目标定
	cv::Mat Rp, Tp, E, F;
	double stereo_error = cv::stereoCalibrate(
		object_points, qc, qp, Kc, kc, Kp, kp, frameSize, Rp, Tp, E, F,
		cv::CALIB_FIX_INTRINSIC + cal_flags,
		cv::TermCriteria(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, 100, DBL_EPSILON));

	//显示标定结果
	std::cout << "Kc = \n" << Kc;
	std::cout << "\nkc = \n" << kc;
	std::cout << "\ncam_error = " << cam_error << std::endl;

	std::cout << "\nKp = \n" << Kp;
	std::cout << "\nkp = \n" << kp;
	std::cout << "\nproj_error = " << proj_error << std::endl;

	std::cout << "\nRp = \n" << Rp;
	std::cout << "\nTp = \n" << Tp;
	std::cout << "\nstereo_error = " << stereo_error;

	std::string filename = "calibration_0.xml";
	save_calibrate_xml(filename, Kc,kc, cam_error,Kp,kp, proj_error,Rp,Tp, stereo_error);
}

void read_calib_images(std::string folder_path, int img_count,
	std::vector<std::vector<cv::Mat>>& multi_calib_images_H, std::vector<std::vector<cv::Mat>>& multi_calib_images_V)
{
	for (int group = 0; group < img_count; group++)   //总共的图像组数
	{
		std::vector<cv::Mat> calib_images_H, calib_images_V;
		std::string folder_path_index = folder_path + std::to_string(group) + "//*.bmp";
		std::vector<cv::String> image_files;
		cv::glob(folder_path_index, image_files);
		int img_index = 0;
		for (const auto& file : image_files) {
			cv::Mat image = cv::imread(file,0);
			if (image.empty()) {
				std::cerr << "Failed to read image: " << file << std::endl;
				continue;
			}
			if (img_index < 20)
				calib_images_V.push_back(image);
			else
				calib_images_H.push_back(image);
			img_index++;
		}
		multi_calib_images_H.push_back(calib_images_H);
		multi_calib_images_V.push_back(calib_images_V);
	}
}

int main()
{
	init();

	//1. 读取标定图像
	std::vector<std::vector<cv::Mat>> multi_calib_images_H, multi_calib_images_V;
	std::string folder_path = ".//calibImage//";
	int img_count = 7;
	read_calib_images(folder_path, img_count, multi_calib_images_H, multi_calib_images_V);
	if (multi_calib_images_V.size() <= 4)
	{
		std::cout << "至少需要4组图像!\n";
		return 0;
	}
	std::vector<cv::Mat> chess_imgs, vps, ups;
	int encode_img_count = multi_calib_images_V[0].size();    //获取每一组编码图像的数量,5个频率,4步相移,其实就是20张图像
	for (int i = 0; i < img_count; i++)
	{
		cv::Mat chess_img = multi_calib_images_H[i][encode_img_count];   //获取最后一张没有条纹的棋盘格图像
		cv::Mat decode_img_V = decode_pattern(multi_calib_images_V[i], v_freq_array, PROJECTOR_HEIGHT_4500);
		cv::Mat decode_img_H = decode_pattern(multi_calib_images_H[i], h_freq_array, PROJECTOR_WIDTH_4500*2);
	
		chess_imgs.push_back(chess_img.clone());
		vps.push_back(decode_img_V.clone());
		ups.push_back(decode_img_H.clone());
	}

	calibrate_Impl(chess_imgs, vps, ups);

	return 0;
}

 部分运行效果(检测到的棋盘格角点)如下所示:

我们将最后一组图像作为测试,进行三维重建,效果如下:

如果大家有什么疑问,欢迎留言,我将一一解答。

相关的测试数据集可以下面链接中下载

https://download.csdn.net/download/laiyinping/89464283

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mfbz.cn/a/731103.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

IO读取properties文件实现JDBC连接池实战

参考文章 Java中的池化思想 面试官&#xff1a;为什么数据库连接很消耗资源&#xff0c;资源都消耗在哪里&#xff1f; 池化思想是什么&#xff1f;连接池是什么&#xff1f; 在Java中&#xff0c;池化思想是一种通过创建和管理可重复使用的对象池来提高性能和资源利用率的编…

【图解IO与Netty系列】Netty编解码器、TCP粘包拆包问题处理、Netty心跳检测机制

Netty编解码器、TCP粘包拆包问题处理、Netty心跳检测机制 Netty编解码器编码器解码器编解码器Netty提供的现成编解码器 TCP粘包拆包问题处理Netty心跳检测机制 Netty编解码器 网络传输是以字节流的形式传输的&#xff0c;而我们的应用程序一般不会直接对字节流进行处理&#x…

建筑驱鸟设备 | 建筑专用超声波驱鸟器

从半夜的鸣叫到频繁的鸟粪污染&#xff0c;鸟类活动有时会成为城市居民不得不面对的小小困扰。通过合理的驱鸟方法&#xff0c;我们可以有效地减少鸟类对建筑物的侵扰&#xff0c;保护建筑物的完好和安全&#xff0c;同时维护城市居民的生活质量。 建筑专用超声波驱鸟器&#x…

理解 JTBD 框架和EJ 理念:深挖以用户为中心的设计

在与用户的交流中&#xff0c;我们发现对用户需求的精准洞察普遍困扰着产品经理、设计、企划人员&#xff0c;因为当今消费者行为已经由单品消费转向场景消费&#xff0c;千人千面的个性化需求出现&#xff0c;消费者数据维度极大丰富&#xff0c;这对把握用户体验造成了很大挑…

Android开发系列(五)Jetpack Compose之Icon Image

Icon是用于在界面上显示矢量图标的组件。它提供了很多内置的矢量图标&#xff0c;也支持自定义图标。要使用Icon组件&#xff0c;可以通过指定图标资源的名称或引用来创建一个Icon对象。例如&#xff0c;使用Icons.Default.Home来创建一个默认风格的首页图标。可以通过设置图标…

TrueNAS系统在ARM平台上的移植

随着家庭及中小型企业对存储和共享需求的日益增长&#xff0c;高效、可靠的文件存储系统成为支撑各类应用的关键。 在众多存储系统中&#xff0c;TrueNAS以其卓越的数据完整性与可靠性、简洁高效的应用程序部署和管理、灵活的虚拟化应用添加能力&#xff0c;以及出色的可用性&a…

【第24章】Vue实战篇之用户信息展示

文章目录 前言一、准备1. 获取用户信息2. 存储用户信息3. 加载用户信息 二、用户信息1.昵称2.头像 三、展示总结 前言 这里我们来展示用户昵称和头像。 一、准备 1. 获取用户信息 export const userInfoService ()>{return request.get(/user/info) }2. 存储用户信息 i…

【面试题】风险评估和应急响应的工作流程

风险评估和应急响应是网络安全管理中两个重要的环节。下面分别介绍它们的工作流程&#xff1a; 一、风险评估工作流程&#xff1a; 1.确定评估范围&#xff1a;明确需要评估的信息系统或资产的范围。 2.资产识别&#xff1a;识别并列出所有需要评估的资产&#xff0c;包括硬件…

美妆短剧撬动33亿市值后,爆款短剧有了新风向

6月1日起微短剧分级备案正式施行&#xff0c;所有短剧未经备案不得播出&#xff0c;该备案也是短剧行业的首个行业规范&#xff0c;让近两年来肆意增长的短剧迎来新一轮洗牌&#xff0c;在保障短剧质量的同时&#xff0c;也促进了行业的发展。 ▲ 图片来源&#xff1a;网络 面对…

Freertos-----任务之间的消息传递(使用消息队列信号量方法)

这次来分享任务之间的数据传递的方法&#xff0c;方法有很多种&#xff0c;我展示2种&#xff0c;让大家对freertos有更深刻的印象 目录 消息队列 信号量 消息队列 首先直接打开普中的例程&#xff0c;然后在里面加上ADC的驱动代码&#xff0c;先初始化外设先&#xff0c;我…

前端模糊搜索关键字高亮

效果 代码 <template><view class"flexStart new-box"><view class"company"><!-- 输入框样式 --><view class"spaceBetween companyName" click.stop"isCompany true"><input type"text&quo…

xargs 传参

xargs的默认命令是 echo&#xff0c;空格是默认定界符。这意味着通过管道传递给 xargs的输入将会包含换行和空白&#xff0c;不过通过 xargs 的处理&#xff0c;换行和空白将被空格取代。xargs是构建单行命令的重要组件之一。 xargs -n1 // 一次输出一个参数到一行&#xf…

qmt量化交易策略小白学习笔记第47期【qmt编程之期货仓单】

qmt编程之获取期货数据 qmt更加详细的教程方法&#xff0c;会持续慢慢梳理。 也可找寻博主的历史文章&#xff0c;搜索关键词查看解决方案 &#xff01; 感谢关注&#xff0c;咨询免费开通量化回测与获取实盘权限&#xff0c;欢迎和博主联系&#xff01; 期货仓单 提示 1…

使用Python selenium爬虫领英数据,并进行AI岗位数据挖掘

随着OpenAI大火&#xff0c;从事AI开发的人趋之若鹜&#xff0c;这次使用Python selenium抓取了领英上几万条岗位薪资数据&#xff0c;并使用Pandas、matplotlib、seaborn等库进行可视化探索分析。 但领英设置了一些反爬措施&#xff0c;对IP进行限制封禁&#xff0c;因此会用到…

mysql 某个时间字段取值时间标识的字符串的值

SELECT STR_TO_DATE(substr(out_trade_no, 1,14), %Y-%m-%d %H:%i:%s) FROM o_order WHERE id 364457; UPDATE o_order SET created_time DATE_FORMAT(STR_TO_DATE(substr(out_trade_no, 1,14), %Y%m%d %H%i%s), %Y-%m-%d %H:%i:%s) WHERE id 364457; 举例&#xff1a; 1…

五种实用方法!手把手教你系统盘瘦身

随着电脑的使用时间变长&#xff0c;电脑硬盘会逐渐被各种类型的数据占满&#xff0c;其中系统盘的可用空间也在慢慢变小。这是因为系统在运行过程中会产生大量临时文件和缓存文件&#xff0c;同时&#xff0c;系统的每一次更新升级也都会生成相关的文件夹存放在系统盘中&#…

高阶图神经网络 (HOGNN) 的概念、分类和比较

图神经网络&#xff08;GNNs&#xff09;是一类强大的深度学习&#xff08;DL&#xff09;模型&#xff0c;用于对相互连接的图数据集进行分类和回归。它们已被用于研究人类互动、分析蛋白质结构、设计化合物、发现药物、识别入侵机器、模拟单词之间的关系、寻找有效的交通路线…

AI/ML 数据湖参考架构架构师指南

这篇文章的缩写版本于 2024 年 3 月 19 日出现在 The New Stack 上。 在企业人工智能中&#xff0c;主要有两种类型的模型&#xff1a;判别模型和生成模型。判别模型用于对数据进行分类或预测&#xff0c;而生成模型用于创建新数据。尽管生成式人工智能最近占据了新闻的主导地…

excel基本操作

excel 若要取消在数据表中进行的所有筛选 步骤操作&#xff1a; 单击“数据”选项卡。在“排序和筛选”组中&#xff0c;找到“清除”按钮。点击“清除”按钮。 图例&#xff1a; 将文本文件的数据导入到Excel工作表中进行数据处理 步骤&#xff1a; 在Excel中&#xff0c…

DevEco鸿蒙开发请求网络交互设置

首先&#xff0c;在鸿蒙项目下config.json中找到module项&#xff0c;在里面填写"reqPermissions": [{"name": "ohos.permission.INTERNET"}] 在页面对应js文件内&#xff0c;填写import fetch from system.fetch;。 GET和POST区别 GET将表单数…