无人机影像配准并发布(共线方程)

news2024/11/16 7:30:50

无人机影像 + DEM 计算四个角点坐标(刚性变换)

  1. 像空间坐标(x,y,-f)

  2. 像空间坐标畸变纠正 deltax,deltay
    在这里插入图片描述

  3. 已知(x,y),求解(X,Y, Z)或者(Lat,Lon)
    在这里插入图片描述
    这里的Z是DEM上获取的坐标和Zs为相机坐标的高程,如果均为已知的情况下,则可以求解(X,Y),这里的(X,Y,Z)为地固地心坐标,单位为米。平地的情况只需要获取行高即可求解(X,Y),接着使用proj库将地固地心坐标转化为经纬度坐标即可。

  4. 地理配准
    这里直接采用**gdal_translategdal_wrap**,gdal_translate转换过程如下,大概就是将jpg进行地理配准。请注意,GDAL的影像起点是左上角,但是我们的相机模型是左下角,所以需要变换Y轴。转换之后的tif只有专业的QGIS之类的软件才能读取。
    在这里插入图片描述

    下面是QGIS读取的效果。但是为了geoserver能够识别,还需要转换,这时候需要gdal_wrap,这个时候就很关键,我们需要设置其为透明。

    gdalwarp -r cubic -ovr AUTO -dstalpha D:\code\roadProj\public\demo\DSC00002\test.tif D:\code\roadProj\public\demo\DSC00002\test3_geotiff.tif
    

    在这里插入图片描述

    gdal_translate.exe -of GTiff -gcp 0 5304 102.1265090139 29.6453703982 -gcp 7952 5304 102.1164515460 29.6474820822 -gcp 7952 0 102.1131839750 29.6445496193 -gcp 0 0 102.1233217949 29.6424804051 -ovr AUTO -co GCPs_Creation=YES D:\code\roadProj\public\demo\DSC00002\DSC00002.jpg D:\code\roadProj\public\demo\DSC00002\test.tif
    
    1. geoserver发布,具体操作比较简单

代码逻辑

  1. 下面是求解影像四个角点经纬度的简单思路,主要还是共线方程,代码中的1000还是得根据距离地面的高度,即需要DEM的高程值才能求解得到较为精确的精度。

  2. 具体实现分为相机模型(固定的参数不部分),大疆无人机是WGS84椭球,EPSG:4978是地心地固的转换参数。

struct CameraModel {
    double f = 7538.508;  // 像素为单位

    double u0 = 3982.417; // 像素为单位
    double v0 = 2671.637;

    double pixelSize = 4.5e-6; // 米为单位

    double k1 = 2.470920e-9;   // 径向畸变系数
    double k2 = -2.767172e-16;
    double k3 = 2.479935e-23;
    double k4 = -6.583598e-31;

    double p1 = 1.388595e-8; // 偏心畸变系数
    double p2 = 1.781812e-7;

    double alpha = -4.697031e-4; // CCD非正方形比例系数
    double beta = -1.300023e-4;

    double width = 7952; // 影像的高度和宽度
    double height = 5304;
};```

```cpp
///
/// \brief The ComputeBoundingBox class
/// 计算影像的包围盒的经纬度坐标 + 高程,然后贴地
///
class ComputeBoundingBox {
public:
    ComputeBoundingBox() {
        transTool.init("EPSG:4326",
                       "EPSG:4978"); // 椭球坐标->地心坐标 XYZ
    }

    ///
    /// \brief resetR
    /// \param phi 俯仰角
    /// \param omega 横滚角
    /// \param kappa 旋转角
    ///
    void resetR(double phi, double omega, double kappa) {
        this->phi = degreeToRadian(phi);
        this->omega = degreeToRadian(omega);
        this->kappa = degreeToRadian(kappa);
    }

    ///
    /// \brief resetCamera
    /// \param lat 维度
    /// \param lon 经度
    /// \param height 高程
    ///
    void resetCamera(const QString &imgNumber, double lat, double lon, double height) {
        this->imgNumber = imgNumber;
        this->cameraLat = lat;
        this->cameraLon = lon;
        this->height = height;
    }

    ///
    /// \brief compute 计算相机位置 + 旋转矩阵
    ///
    void compute();

private:
    ///
    /// \brief computeLatlon 根据像点坐标计算经纬度坐标(共线方程的逻辑)
    /// @param vector2d uv x轴,y轴坐标
    ///
    Eigen::Vector2d computeLatlon(Eigen::Vector2d &uv);

    ///
    /// \brief degreeToRadian 度转弧度
    /// \param degree
    /// \return
    ///
    double degreeToRadian(double degree) { return degree * M_PI / 180.0; }

    //    double computeDeltaX();

    //    double computeDeltaY();

    Eigen::Vector3d cameraGeo;

    Eigen::Matrix3d rMatrix;

    Transform transTool;

    CameraModel intrinsic; // 内参数矩阵

    QString imgNumber;     // 影像编号

    double cameraLat;      // 相机外参数矩阵
    double cameraLon;
    double height;

    double phi;
    double omega;
    double kappa;
};
Eigen::Vector2d ComputeBoundingBox::computeLatlon(Eigen::Vector2d &uv) {
    double u = uv.x();
    double v = uv.y();
    Eigen::Vector3d cameraSpace(u, v, -this->intrinsic.f); // 像空间坐标

    double r = qSqrt(pow(u - this->intrinsic.u0, 2) + pow(v - this->intrinsic.v0, 2));

    // (x-x0) * (k1*r^2 + k2*r^4 + k3*r^6 + k4*r8) + p1*(r^2 + 2*(x-x)^2) + 2p2*(x-x0)(y-y0) + alpha*(x-x0) +
    // beta*(y-y0)
    double deltaX = (u - this->intrinsic.u0) * (this->intrinsic.k1 * pow(r, 2) + this->intrinsic.k2 * pow(r, 4) +
                                                this->intrinsic.k3 * pow(r, 6) + this->intrinsic.k4 * pow(r, 8)) +
                    this->intrinsic.p1 * (pow(r, 2) + 2 * pow((u - this->intrinsic.u0), 2)) +
                    2 * this->intrinsic.p2 * (u - this->intrinsic.u0) * (v - this->intrinsic.v0) +
                    this->intrinsic.alpha * (u - this->intrinsic.u0) + this->intrinsic.beta * (v - this->intrinsic.v0);
    double deltaY = (v - this->intrinsic.v0) * (this->intrinsic.k1 * pow(r, 2) + this->intrinsic.k2 * pow(r, 4) +
                                                this->intrinsic.k3 * pow(r, 6) + this->intrinsic.k4 * pow(r, 8)) +
                    this->intrinsic.p2 * (pow(r, 2) + 2 * pow((v - this->intrinsic.v0), 2)) +
                    2 * this->intrinsic.p1 * (u - this->intrinsic.u0) * (v - this->intrinsic.v0);

    Eigen::Vector3d cameraOffset(deltaX - this->intrinsic.u0, deltaY - this->intrinsic.v0,
                                 0);                              // 像点坐标偏移
    Eigen::Vector3d cameraSpaceTrue = cameraSpace + cameraOffset; // 实际的像点坐标[最后一位该如何求解]

    Eigen::Vector3d worldCoordBa =
        this->rMatrix * cameraSpaceTrue * this->intrinsic.pixelSize; // (xBa, yBa, zBa) pixelSize这个参数多余

    worldCoordBa =
        Eigen::Vector3d(worldCoordBa.x() / worldCoordBa.z() * 1000, worldCoordBa.y() / worldCoordBa.z() * 1000, 0);

    //    qDebug() << worldCoordBa.x() << " " << worldCoordBa.y() << " " << worldCoordBa.z();

    Eigen::Vector3d worldCoord = worldCoordBa + this->cameraGeo; // 真正的坐标
    //    std::cout << worldCoord.x() << " " << worldCoord.y() << " " << worldCoord.z();

    //    PJ_COORD latlonh = proj_coord(cameraLon, cameraLat, height, 0);
    PJ_COORD geoxyz = proj_coord(worldCoord.x(), worldCoord.y(), worldCoord.z(), 0); // 地心坐标
    PJ_COORD latlon = this->transTool.backward(geoxyz);
    //    std::cout << "lat:" << latlon.lpz.lam << " ,lon:" << latlon.lpz.phi << ", " << latlon.lpz.z;
    return Eigen::Vector2d(latlon.lp.phi, latlon.lp.lam); // lat and lon
}
void ComputeBoundingBox::compute() {

    PJ_COORD latlonh = proj_coord(cameraLon, cameraLat, height, 0);
    PJ_COORD geoxyz = transTool.forward(latlonh);

    this->cameraGeo = Eigen::Vector3d(geoxyz.xyz.x, geoxyz.xyz.y, geoxyz.xyz.z); // 地心坐标

    // 计算绕X轴旋转的旋转矩阵 Rx(φ)
    Eigen::Matrix3d Rx;
    Rx << 1, 0, 0, 0, std::cos(phi), -std::sin(phi), 0, std::sin(phi), std::cos(phi);

    // 计算绕Y轴旋转的旋转矩阵 Ry(ω)
    Eigen::Matrix3d Ry;
    Ry << std::cos(omega), 0, std::sin(omega), 0, 1, 0, -std::sin(omega), 0, std::cos(omega);

    // 计算绕Z轴旋转的旋转矩阵 Rz(κ)
    Eigen::Matrix3d Rz;
    Rz << std::cos(kappa), -std::sin(kappa), 0, std::sin(kappa), std::cos(kappa), 0, 0, 0, 1;

    // 计算总的旋转矩阵 R_total = Rz(κ) * Ry(ω) * Rx(φ)
    this->rMatrix = Rz * Ry * Rx;

    // 像素点畸变纠正

    Eigen::Vector2d lb(0, 0);                                          // 左下角为起点
    Eigen::Vector2d rb(this->intrinsic.width, 0);                      // 右下角坐标
    Eigen::Vector2d rt(this->intrinsic.width, this->intrinsic.height); // 右上角坐标
    Eigen::Vector2d lt(0, this->intrinsic.height);                     // 左上角成果

    std::vector<Eigen::Vector2d> latlonVec = {lb, rb, rt, lt};
    qDebug() << "####################" << this->imgNumber << "########################";
    for (Eigen::Vector2d &pixelCoord : latlonVec) {
        pixelCoord = this->computeLatlon(pixelCoord);
        qDebug() << "lat: " << QString::number(pixelCoord.x(), 'f', 10)
                 << ",lon:" << QString::number(pixelCoord.y(), 'f', 10);
    }
}

效果图

在这里插入图片描述

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

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

水文章——推荐一个视频播放器和一个图片查看器

视频播放器——PotPlayer http://www.potplayercn.com/ 图片查看器——JPEGVIEW https://www.bilibili.com/video/BV1ZY411P7fX/?spm_id_from333.337.search-card.all.click&vd_sourceab35b4ab4f3968642ce6c3f773f85138

PHP数组转对象和对象转数组

PHP数组转对象和对象转数组 <?php function array_to_object($arr){$obj new stdClass();foreach ($arr as $key > $val) {if (is_array($val) || is_object($val)) {$obj->$key array_to_object($val);} else {$obj->$key $val;}}return $obj; } function o…

pdf怎么转换成word 文档?这几种方法收藏一下

pdf怎么转换成word 文档&#xff1f;PDF和Word是我们平时工作和学习中最常用的两种文档格式。PDF文档格式通常用于电子书籍、合同、申请表等需要保持原样式的文档。而Word文档则通常用于编辑和修改。但是&#xff0c;有时我们需要将PDF文档转换为可编辑的Word文档&#xff0c;以…

【Docker】在Docker大火的背后,究竟隐藏着未来科技发展的哪些大趋势

这里写目录标题 在docker大火的背后是什么新科技的发展呢&#xff1f;1.容器化技术的普及2.云原生应用的兴起3.边缘计算的发展4.容器编排和管理平台的演进5.混合云和多云架构的普及 docker三大特性轻量化可移植性可扩展性 docker被谁抢了风头呢1.风头被Kubernetes (K8S)抢了2.缺…

日撸代码300行:第54天(基于 M-distance 的推荐)

代码来自闵老师”日撸 Java 三百行&#xff08;51-60天&#xff09;“&#xff0c;链接&#xff1a;日撸 Java 三百行&#xff08;51-60天&#xff0c;kNN 与 NB&#xff09;_闵帆的博客-CSDN博客 算法是基于M-distance的推荐&#xff0c;通过用户评分矩阵对用户进行电影推荐。…

如果你在选型低代码平台,可以从这5个角度去分析抉择

研究低代码平台已有3年&#xff0c;也算是个低代码资深用户了&#xff0c;很多企业面临低代码选型上的困难&#xff0c;选平台容易&#xff0c;换平台难。下面基于个人理解给大家做一份千字的注意事项&#xff01;希望对大家在选型低代码方面有一定帮助。最终&#xff0c;正确且…

[AWD靶场搭建]

文章目录 [AWD靶场搭建]前言AWD平台搭建靶机搭建Cadinal添加靶机 连接Asteroid大屏默认ssh账号密码参考 [AWD靶场搭建] 前言 觉得好玩搭建了一下AWD靶场&#xff0c;使用了vidar-team编写的 Cardinal AWD平台搭建 这里我是在kali搭建的&#xff0c;所以我下载了这个压缩包&…

centos7搭建k8s环境并部署springboot项目

之前看了很多文章&#xff0c;都是部署后一直报错&#xff0c;百度解决后下次又忘了&#xff0c;这次决定把从头到尾的过程记录下来方便下次再看&#xff0c;部署参考文章尚硅谷Kubernetes&#xff08;k8s&#xff09;视频学习笔记_尚硅谷k8s笔记_溯光旅者的博客-CSDN博客 1、…

13年测试老鸟,接口性能测试总结整理,据说这是全网最全的...

目录&#xff1a;导读 前言一、Python编程入门到精通二、接口自动化项目实战三、Web自动化项目实战四、App自动化项目实战五、一线大厂简历六、测试开发DevOps体系七、常用自动化测试工具八、JMeter性能测试九、总结&#xff08;尾部小惊喜&#xff09; 前言 性能测试按照不同…

支持多种通信方式和协议方便接入第三方服务器或云平台

2路RS485串口是一种常用的通信接口&#xff0c;可以支持Modbus Slave协议&#xff0c;并可接入SCADA、HMI、DSC、PLC等上位机。它还支持Modbus RTU Master协议&#xff0c;可用于扩展多达48个Modbus Slave设备&#xff0c;如Modbus RTU远程数据采集模块、电表、水表、柴油发电机…

GAN论文精读

标题:Generative Adversarial Nets 摘要: 简写:作者提出了一个framework通过一个对抗的过程&#xff0c;在这里面会同时训练两个模型。 第一个模型为生成模型G&#xff0c;是用来抓住整个数据的分布 第二个模型为辨别模型D&#xff0c;是用来估计一个样本是否从G中产生。 …

BD Biosciences通过使用Liquid UI优化SAP QM,节省了80%的处理时间,提高了 95% 的数据准确性

背景 BD 生物科学公司成立于 1897 年&#xff0c;致力于改善患者的治疗效果&#xff0c;并在一个多世纪的时间里始终坚持这一理念&#xff0c;现已涉足诊断、生物科学以及各种医疗设备和仪器系统。 挑战 手动验证数据 原因&#xff1a;使用非自动程序演示和验证数据&#xff0c…

FRR+VPP

安装 三者的结合&#xff0c;实际上编译安装好就行了&#xff0c;不需要做任何代码上的修改&#xff0c;只需要安装和配置&#xff0c;然后你就有了一台路由器。 FRRouting使用frr-8.5.2版本&#xff0c;VPP使用23.06版本&#xff0c;DPDK和lcpng是VPP的插件&#xff0c;安装…

【CAS6.6源码解析】源码构建时-默认service配置不生效解决方案

CAS6的源码提供了默认的HTTPSandIMAPS-10000001.json配置用于授权所有的https和imaps服务&#xff0c;但是当添加JsonServiceRegistry模块启动后&#xff0c;会发现service是没有被注册的&#xff0c;是由于json路径引起的错误&#xff0c;可以把路径修改为绝对路径以解决此问题…

在idea中添加try/catch的快捷键

在idea中添加try/catch的快捷键 在idea中添加try/catch的快捷键 ctrlaltt 选中想被try/catch包围的语句&#xff0c;同时按下ctrlaltt&#xff0c; 出现下图 选择try/catch即可。

QTDAY3

闹钟 头文件 #ifndef WIDGET_H #define WIDGET_H#include <QWidget> #include <QTimerEvent> //定时器事件处理函数 #include <QTime> //时间类 #include <QString> #include <QPushButton> #include <QTextToSpeech> #include …

spring中存储和获取bean对象

文章目录 1. 存储bean对象1.1 使用配置文件存储bean对象1.2 使用五大类注解存储bean对象1.2.1 类注解1.2.2 五大类注解的作用1.2.3 方法注解 2.获取bean对象2.1 属性注入2.2 构造器注入2.3 getter注入2.4 注入对象的时候有spring中有多个bean对象怎么办2.4.1 将类中属性的名字&…

ava版知识付费平台免费搭建 Spring Cloud+Spring Boot+Mybatis+uniapp+前后端分离实现知识付费平台

提供私有化部署&#xff0c;免费售后&#xff0c;专业技术指导&#xff0c;支持PC、APP、H5、小程序多终端同步&#xff0c;支持二次开发定制&#xff0c;源码交付。 Java版知识付费-轻松拥有知识付费平台 多种直播形式&#xff0c;全面满足直播场景需求 公开课、小班课、独…

数据结构基础:2.顺序表。

顺序表的介绍和实现 一.线性表1.基本概念&#xff1a; 二.顺序表&#xff1a;1.基本概念&#xff1a;分类&#xff1a;1.静态顺序表&#xff1a;分类&#xff1a;2.动态顺序表&#xff1a;2.动态顺序表的功能接口的实现&#xff1a;0.顺序表打印&#xff1a;1.初始化和删除&…

功率放大器在电光调制中的应用有哪些

电光调制是一种利用光电效应将电信号转化为光信号的技术。在实现电光调制的过程中&#xff0c;功率放大器作为一个重要的组件&#xff0c;具有对输入电信号进行放大和控制的功能。本文将介绍功率放大器的基本原理、特点以及在电光调制中的应用。 基本原理 功率放大器是一种能够…