用数字描绘科学发展
为决策服务创先争优——国家统计局江西调查总队--《理论导报》2010年10期
用数字描绘科学发展
为决策服务创先争优——国家统计局江西调查总队
【摘要】：正国家统计局江西调查总队是国家统计局的派出机构,为正厅级机构,内设17个处室,在全省11个设区市设置正县处级调查队、38个县(市、区)设置正科局级调查队。主要职责是:组织实施主要农产品产量、农作物播种面积、城乡居民收入、农户固定资产投资、工业品价格、固定资产投资价格、居民消费价格、农产品价格、房地产价格、规模以下工业、企业景气、采购经理、畜牧业监测、退耕还林监测、农村贫困监测、农民
【关键词】：
【分类号】：C829.2【正文快照】：
国家统计局江西调查总队是国家统计局的派出机构，为正厅级机构，内设17个处室，在全省n个设区市设置正县处级调查队、38个县(市、区)设置正科局级调查队。主要职责是:组织实施主要农产品产量、农作物播种面积、城乡居民收入、农户固定资产投资、工业品价格、固定资产投资价
欢迎：、、)
支持CAJ、PDF文件格式，仅支持PDF格式
【相似文献】
中国期刊全文数据库
田兆静;;[J];中国统计;2011年07期
孔德全;杨建军;;[J];内蒙古统计;2011年03期
吴波;;[J];统计与咨询;2011年03期
杨汉东;;[J];统计与决策;2011年15期
本刊编辑部;;[J];统计与管理;2011年03期
胡继红;何雯;;[J];数据;2011年09期
杨斌;谢慧;;[J];中国统计;2011年04期
高蓝星;;[J];调研世界;2011年07期
;[J];共产党人;2011年01期
潘晓明;;[J];中国性科学;2011年08期
中国博士学位论文全文数据库
王克林;[D];暨南大学;2011年
&快捷付款方式
&订购知网充值卡
400-819-9993
《中国学术期刊（光盘版）》电子杂志社有限公司
同方知网数字出版技术股份有限公司
地址：北京清华大学 84-48信箱 大众知识服务
出版物经营许可证 新出发京批字第直0595号
订购热线：400-819-82499
服务热线：010--
在线咨询：
传真：010-
京公网安备75号查看: 3350|回复: 6
DLG、DRG、DEM、DOM概述
马上注册地信网，享受更多功能，学习更多知识，成就人生精彩！
才可以下载或查看，没有帐号？
本帖最后由 丈量大地 于
00:57 编辑
DLG、GRG、DEM、DOM都是数字化的产品，其概述如下：
数字线划地图（Digital Line Graphic，缩写DLG）
数字线划地图（Digital Line Graphic，缩写DLG）
1.jpg (273.63 KB, 下载次数: 33)
00:53 上传
数字线划地图是地形图上现有核心要素信息的矢量格式数据集，可保存各要素间的空间关系和相关的属性信息，全面地描述地表目标。内容包括行政界线、地名、水系及水利设施工程、交通网和地图数学基础（高斯坐标系和地理坐标系）。数字线划地图是满足地理信息分析要求的数据结构，可视为带有智能的数据集，不但含有几何数据，还有社会人文信息。
数字栅格地图（Digital Raster Graphic，缩写DRG）
2.jpg (123.11 KB, 下载次数: 30)
00:54 上传
数字栅格地图是现有纸质地形图经计算机处理后得到的栅格数据文件。每一幅地形图在扫描数字化后，经几何纠正，并进行内容更新和数据压缩处理，彩色地形图还应经色彩校正，使每幅图像的色彩基本一致。数字栅格地图在内容上、几何精度和色彩上与国家基本比例尺地形图保持一致。
数字高程模型（Digital Elevation Model，缩写DEM）
3.jpg (231.44 KB, 下载次数: 32)
00:55 上传
　　数字高程模型是在某一投影平面（如高斯投影平面）上规则格网点的平面坐标（X，Y）及高程（Z）的数据集。DEM的格网间隔应与其高程精度相适配，并形成有规则的格网系列。数字高程模型是描述地表起伏形态特征的空间数据模型，根据不同的高程精度，可分为不同类型。为完整反映地表形态，还可增加离散高程点数据。
　数字正射影像（Digital Orthophoto Map，缩写DOM）
4.jpg (107.05 KB, 下载次数: 33)
00:56 上传
看看，，，，，了解下
好资料！！谢谢提供
4D产品概述
随着测绘技术和计算机技术的结合与不断发展，地图不在局限与以往的模式，现代数字地图主要由DOM(数字正射影像图)、DEM(数字高程模型)、DRG(数字栅格地图)、DLG(数字线划地图)以及复合模式组成。
DOM (数字正射影像图)
　　利用航空相片、遥感影像，经象元纠正，按图幅范围裁切生成的影像数据。 它的信息丰富直观，具有良好的可判读性和可量测性，从中可直接提取自然地理和社会经济信息。
DEM (数字高程模型)
　　数字高程模型是以高程表达地面起伏形态的数字集合。可制作透视图、断面图，进行工程土石方计算、表面覆盖面积统计，用于与高程有关的地貌形态分析、通视条件分析、洪水淹没区分析。
DRG (数字栅格地图)
　　数字栅格地图是纸制地形图的栅格形式的数字化产品。 可作为背景与其他空间信息相关，用于数据采集、评价与更新，与DOM、DEM集成派生出新的可视信息。
DLG (数字线划地图)
　　现有地形图上基础地理要素分层存储的矢量数据集。数字线划图既包括空间信息也包括属性信息，可用于建设规划、资源管理、投资环境分析等各个方面以及作为人口、资源、环境、交通、治安等各专业信息系统的空间定位基础。
4D数字产品采集方法
　　数字高程模型(DEM)
　　数字高程模型的数据采集通常包括以下几种方法：
　　利用自动记录的测距经纬仪（常用电子速测经纬仪或全丫经纬仪）的野外实测。这种速测经纬仪一般都有微处理器，可以自动记录和显示有关数据，还能进行多种测站上的计算工作。其记录的数据可以通过串行通讯，输入计算机中进行处理。
现有地图数字化
　　利用数字化仪对已有地图上的信息(如等高线)进行数字化的方法,目前常用字数字化仪有手扶跟踪数字化仪和扫描数字化。
空间传感器
　　利用全球定位系统GPS,结合雷达和激光测高仪等进行数据采集。早在2000年，美国“奋进”号航天飞机在结束了9天的绕地飞行后，采用星载成像雷达和合成孔径雷达等高新技术，采集了地球上人类所能正常活动地区（约占地表总面积的80%）的地面高程信息，经处理可制成数字高程模型和三维地形图。此次计划所取得的测绘成果，覆盖面大、精度高、有统一的基准，不但在民用方面应用广泛，而且在导弹发射、战场管理、后勤规划等军事活动中具有重要价值，因此引起了各国军界和传媒的广泛关注。
数字摄影测量方法
　　这是DEM数据采集最常用最有效的方法之一。利用附有的自动记录装置（接口）的立体测图仪或立体坐标仪、解析测图仪及数字摄影系统，进行人工、半自动或全自动的量测来获取数据。
LIDAR +CCD相机
　　LIDAR也叫机载激光雷达，是一种安装在飞机上的机载激光探测和测距系统，是由GPS（全球卫星定位系统）、INS（惯性导航系统）和激光测距三大技术的集成应用系统。如加拿大OPTECH公司生产的ALTM3100系统和德国 IGI公司生产的LiteMapper5600系统，ALTM3100和LiteMapper5600机载激光扫描遥感系统同时还集成了CCD相机，它与激光探测与测距系统协同作业，同步纪录探测点位的影像信息，因此它可直接获取一个地区高精度的数字高程模型（DEM）、数字地表模型（DTM）、数字正射影像图（DOM）, 由于这种方法可以直接获取高精度的正射影像数据，免去了影像处理的环节，它的成果可以广泛应用于城市测绘、规划、林业、交通、电力、灾害等部门。
　　数字栅格地图（DRG）
　　数字栅格地图是通过一张纸质或其他质地的模拟地形图，由扫描仪扫描生成一维阵列影像，同时对每一系统的灰度（或分色）进行量化，再经二值处理、图形定向、几何校正即形成一幅数字栅格地图，需要经过以下几个步骤：
　　（1）图形扫描：采用扫描分辨率不低于500dpi的单色或彩色扫描仪扫描。
　　（2）图幅定向：将栅格图幅由扫描仪人材变换为高斯投影平面直角坐标。
　　（3）几何校正：消除图底及扫描产生的几何畸变。可以采用相关软件对栅格图像的畸变进行纠正，纠正时要按公里格网进行，通过仿射变换及双线性变换，实现图幅纠正。
　　（4）色彩纠正：用PHOTOSHOP等软件进行栅格图编辑轰动单色图按要素人工设色，对彩色图作色彩校正，为使色彩统一，应按规定的RGB比例选择所用的几种色调。
　　数字正射影图（DOM）
　　数字正射影像图是利用数字高程模型（DEM）对经扫描处理的数字化航空像片或高空采集的卫星影像数据，逐像元进行投影差改正、镶嵌，按国家基本比例尺地形图图幅范围剪裁生成的数字正射影数据集。
　　对于航空像片,利用全数字摄影系统,恢复航摄时的摄影姿态,建立立体模型,在系统中对DEM进行检测、编辑和生成，最后制作出精度较高的DOM。
　　对于卫星影像数据，可利用已有DEM数据，通过单片数字微分纠正生成DOM数据。
　　数字线划地图（DLG）
　　数字线划地图是基础地理信息系统的核心数字产品，它采集方法很多，主要包括以下几个方面：
　　（1）平板仪测量：夹板仪测量采集的是非数字产品，它最终生成的成果是纸质或薄膜地图。它要生成DLG还需要经过内业数字化、编辑处理。目前，平板仪测量已经不是GIS野外数据采集的主要手段，它正逐渐被全野外数字测量所取代。
　　（2）全野外数字测量：利用电子手簿、便携机或掌上电脑与全站仪相连，测量结果直接以数字形式存储，不需要经过内业数字化处理。
　　（3）GPS测量：采用实时动态GPS测量系统，用两台或更多台GPS接收机来协同工作，将一台接收机作为基站，放在已知点上，其他接收机对空间目标测量，采集的数据存放便携电脑或掌上电脑中。
　　（4）地图数字化：地图数字化有两种作业方式：手扶跟踪数字化和扫描矢量化。手扶跟踪数字化是使用数字化仪进行地图数字化，扫描矢量化是通过专用软件对扫描处理后的数字栅格地图进行屏幕跟踪矢量化。
　　（5）摄影测量：摄影测量经历了模拟摄影测量和解析摄影测量阶段，随着计算机技术及其应用的发展以及数字图像处理、模式识别、计算机视觉等学科的发展，现已进入数字摄影测量阶段。长期以来，摄影测量在基本比例图生产中占据着胡里非常重要的位置，特别是发展到今天的数字摄影量阶段，摄影测量以其高效快速、生成数据产品齐全而发挥着着其他测量手段无法比拟的作用。
4D产品的表示方法
数字高程模型（DEM）
　　数学方法
　　用数学方法来表达，可以采用整体拟合方法即根据区域所有的高程点数据，用傅立叶级数和高次多项式拟合地面高程曲面。也是用局部拟合方法，将地表复杂表面正方形规则区域或面积大致相关的不规则区域，根据有限个点拟合形成高程曲面。
　　图形方法
　　线模式
　　等高线是表示地形最常见的形式。其他的地形特征线也是表达地面高程的重要信息源，如山脊线、谷底线、海线及坡度变换线等。
　　点模式
　　用离散采样数据点建立DEM是常用的方法之一。数据采样可以按规则格网采样，可以是密度一致的或不一致的；可以是不规则采样，如不规则三角网、邻近网模型等；也可以有选择性采样，采集山峰、洼坑、隘口、边界等重要特征点。
数字栅格地图（DRG）
　　数字栅格地图的表示方法一般用数学矩阵表示，以简单的文件格式存储在计算机中。栅格地图每一位置点应矩阵中的一个元素，栅格地图的最终产品以LZW压缩存储。
数字正射影像图（DOM）
　　DOM是利用数字化的航空像片影像或数字遥感影像，经几何改正和镶嵌，按一定图幅范围生成的数字正射影像集。DOM包含三个数据层：注记层、格网修饰层、正射影像数据层，其中前两层为矢量数据、最后一层为栅格数据。正射影像地图使用的是一中矢量——栅格混合数据模型表示，这使得它既具有矢量格式数据的高精可量测特征，又具有景仰数据直观、信息丰富等特征。
数字线划地图（DLG）
　　数字线划地图采用矢量模型表示，矢量模型有两种：面片（Spaghetti）数据模型和拓扑数据模型。
　　面片模型中点用空间坐标（x,y,z）表示，线用组空间坐标表示，而面有由封闭的线（即多边形）来表示。矢量数据的属性信息在计算机中利用附加不同的编码来区分在表现号来描绘；面可以用不同图案、不同颜色来表示。
　　在拓扑模型中，多边形的被分割成一系列的弧和节点，节点、弧及多边形之间空间关系在属性表中定义。各种空间要素的拓扑关系包括点—点、点—线、点—面、线—线、线—面、面—面等多种形式。
4D产品及GIS应用
　　DEM的应用，在于通过计算机采用一定的算法，能够很方便地将DEM数据转换为等高线、透视图、断面图、坡度图以及专题图等各种衅解产品，或者按照用户的需求计算出体积、空间距离、表面覆盖面积等工程数据和统计数据以及进行通视分析、域特征地貌与地形自动分割等。
　　DRG可作为背景用于数据参照或修测其他地理相关信息，用于数字线划地图（DLG）的数据采集、评价和更新，还可与数字高程模型（DEM）、数字正射影像图（DOM）等数据信息集成使用，派生出新的可视信息，从而提取、更新地图数据，绘制纸质地图。
　　DOM具有精度高、信息丰富、直观真实等优点，可用作为背景控制信息，评价其他数据的精度、现实性和完整性；可从中撮自然资源和社会发展信息，为防止灾害和公共设施建设规划等提供可靠依据；还可从中提取和派生新的信息，实现地图的修测更新。在城市测绘领域，DOM被广泛应用于城市规划设计、交通规划设计、城市绿化覆盖率调查、城市建成区发展调查、风景名胜区规划、城市发展与生态环境调查与可持续发展研究等诸多方面，取得了显著的社会与经济效益。
　　DLG满足地理信息系统进行各种空间分析要求，视为带有智能的数据。可随机地进行数据选取和显示，与其他几种产品叠加，便于分析、决策。
　　4D产品构成了地理信息系统的基础数据框架，是其他信息空间载体，用户可依据自身的要求，选择适合自己的基础数据产品，研制各种专题地理信息系统。如在电力管理信息系统中引入适当的GIS系统，可以为电力管理提供行之有效的辅助决策方法；地理信息系统应用于地名数据库管理，提高人们对城市的监控能力；地理信息系统应用于房地产管理，将空间数据与大量的非空间数据（属性数据）结合起来，为维护主地产市场正常、高效运行发挥重要作用；地理信息系统应用于规划管理，提高办事效率，而且利用矢量数据（DLG）与栅格数据（DOM）相结合，使GIS的信息表达更加丰富，形象生动，而且为系统交通、管线、通信、银行、土地等部门。
　　随着国家将基础测绘列入国民经济和社会发展计划，全国许多省市政府把基础地理数据产品建设作为省级基础测绘的重点，目前，全国1：5万基础地理数据的更新和建库已经完成，1：1万基础地理数据的更新和建库正在建设中。
　　做为城市测绘部门，在不断完善大比例尺基础地理数据产品的更新和建库，建立一个良性数据更新，维护体系的同时，应建立基于中小比例尺数字产品的基础地理产品库，并在此基础上做深层次的开发应用，最终纳入到基础地理信息系统的管理中。目前数字城市作为城市建设的一个热点，已得到各级政府的广泛重视，有些地区已进入前期的实施阶段，基础地理信息数据库作为数字城市的基础框架，在数字城市的建设中发挥着重要作用。
太好了&&太感动了
值得一看，好东西，学习学习
论坛推荐话题 /1
地信论坛一直是行业同类平台中的领先者，为了更好的领跑和更适应移动互联，我们决定重新打造地信论坛的微信公众号，将长期以来广大网友为论坛所奉献的丰富知识与经验精选给大家，帮助大家开心的学习，快乐的工作。
热线: 400-
Powered byGIS中DEM的表示方法
１引言数字高程模型——ＤＥＭ是表示空间连续起伏的数字模型。从数学上讲，高程模型是一个二维的连续函数，因此，ＤＥＭ是高程模型的离散化表示。一般，地理信息系统中的数据结构只具有二维的意义，高程是地理空间的第三维坐标，因此ＤＥＭ的建立是一个必要的补充。它是建立不同层次的资源与环境地理信息系统、军事地理信息系统不可缺少的组成部分。ＤＥＭ的应用十分广泛，如铁路、公路勘测设计中的选线，坡度计算、土方计算和隧道设计；水利工程中的水库设计，运河和渠道的布线；对不同地形作静态分析与比较；生成坡度图、坡向图和晕渲图等；用于地貌分析，计算侵蚀和径流等；精确制导武器的地形匹配，无人驾驶飞机的导航；军事目的的三维立体地形显示；战场环境的分析计算，如通视概率的计算，坡度、坡向的计算；此外，若用其它特征代替高程后，还可进行诸如人口、降水量等的分析。还有一个与ＤＥＭ意义非常相近的词：数字地形模型。数字地形模型包括地面起伏和属性两个含义，它不仅有地表的高程数据，...&
(本文共4页)
权威出处：
对实际地面进行描述与分析总是首先以单点的形式获取地形数据，并从点到线、从线到片、再由面片镶嵌而表达整个地形表面。每一个数据点包括３个空间坐标向量（Ｘ，Ｙ，Ｚ）和相应的属性信息。数据点在（Ｘ，Ｙ）坐标平面内的分布不外乎两种基本形式：规则的正方形栅格形式和随机形式。规则分布的数据一般来自自动或半自动摄影测量与遥感方式或者随机分布的数据内插。而实际测绘生产中的绝大部分原始数据都是随机分布的，如地面测量数据、解析测图仪作业和既有地图数字化产生的数字等高线数据等。由于每一个栅格数据点都有明确的邻域（如四邻域），因而栅格数据可以自动连成四边形网络。对于随机分布的数据则存在一个选择最邻近的点组成多边形格网的问题。由于基于三角形的描述是刚体变换不变的，适合于各种数据分布密度，有利于更新和直接利用各种地形特征信息进行数据分析，因此随机三角形格网（ＴＩＮ）被广泛用于随机分布数据的ＤＴＭ的建立。特别是Ｄｅｌａｕｎａｙ三角形格网，由于其唯一性和良好的三...&
(本文共4页)
权威出处：
1引言不规则三角网(Triangulated IrregularNetwork,TIN)在地理信息系统的发展与应用中有着相当重要的作用和地位。地学领域中,野外采集的离散数据在很多情况下不能有效地满足构建DTM(数字地面模型)的需要,特别是在地形复杂的地区,需要进行加密内插,才能使构建的DTM真实地反映地表。使用不规则三角网来构建DTM能有效地解决这一问题,不规则三角网通过从不规则分布的数据点生成的连续三角面来逼近地形表面,其不仅能以不同层次的分辨率来描述地形表面,而且在某一特定分辨率下能用更少的空间和时间更精确地表示更加复杂的表面,特别当地形包含有大量特征线如断裂线、构造线等时,TIN模型能更好的顾及这些特征,从而能更精确合理地表达地表形态[1]。Delaunay三角网在地形拟合方面表现相当出色,它是相互邻接且相互不重叠的三角形的集合,因其具有空外接圆性质以及最小角最大的性质,可最大限度的保证网中三角形满足近似等边(角)性,避免...&
(本文共3页)
权威出处：
1传统不规则三角网等值线的追踪算法及不足常见的不规则三角网等值线的追踪算法为:(1)设定等值特征值Z,判断三角形是否有等值点通过。(2)先搜索边界三角形的边界边,查找是否有等值线经过,然后搜索内部三角形,判断是否有等值线通过。这种算法的局限性有:a、建立不规则三角网时必需设置边界三角形以及边界边的标识,增加存储空间;b、需对边界三角形和内部三角形分别追踪等值线,增加程序结构;c、需对追踪处理过的三角形进行标识,当有多条等值线时处理较为麻烦,增加判断条件;d、在追踪等值线时必需识别等值线穿过三角形的入口和出口,以及共用一边的两个三角形的联系,需要大量的指针,增加程序复杂程度。2新的等值线追踪算法2·1基本思路利用C++语言,根据成熟的不规则三角网建立算法建立三角网。对给定的等值特征点Z求三角形边上的等值点的平面坐标,并将三角形的ID赋予等值点的属性值。将等值点存储在向量容器中,如果向量容器中有与该等值点平面坐标相等的点,则将这两点...&
(本文共2页)
权威出处：
三角网是由一系列连续三角形构成的网状的平面控制图形,是三角测量中布设连续三角形的主要扩展形式,同时向各方向扩展而构成网状。目前,建立不规则三角网TIN(Triangulated IrregularNetwork,TIN)的算法主要有递归生长法、凸闭包收缩法、逐点插入法。优点为点位分布均匀、各点之间相互牵制、图形强度较高,缺点是扩展较缓慢。当数据量较大时,比如机载激光雷达数据,少则几万,多则上千万个点,以上各种算法都已无法满足建网的需求,效率低下且效果不佳[1]。本文将利用虚拟三角网对机载激光雷达数据进行虚拟分割并重新组织选取地面种子点,由地面种子点建立初始三角网,再结合逐点插入法建立不规则三角网,为构建高精度数字地面模型提供较好的数据源。1虚拟三角网Cho等[2]在利用机载激光雷达数据进行建筑物探测和重建时,将虚拟网格的概念引入数据组织,克服了运算效率低和信息损失的问题。虚拟网格就是将原始的点云数据进行虚拟“分割”,按照点云的空...&
(本文共5页)
权威出处：
0引言数字地面模型(digital terrain models,简称DTM)是以数字的形式存储地球表面上所有信息的总和,是描述地面诸特征空间分布的数值的集合。主要用于描述地面起伏的状况,可用于提取各种地形参数,如坡度、坡向、粗糙度等,并进行通视分析、流域结构生成等应用分析[1-2]。DTM的表示方法主要有3种:等高线、规则格网和不规则三角网(TIN)。其中,TIN具有高精度、高效率和易于处理地形线等特点,它是按一定的规则将离散点连接成覆盖整个区域且互不重叠、结构最佳的三角形,建立离散点之间的空间关系。其目的是克服利用离散点计算地球表面上任意点高程的困难,近年来得到较快的发展[3-4]。TIN的构建方法已经成为G IS的一个研究热点,一些改进算法也陆续的被提出。1 D elaunay三角网及传统算法Delaunay三角网具有良好的形态,在表达地表形态方面表现较为出色[5],是目前公认的最优的三角网,它是一系列相连的但不重叠的三角...&
(本文共3页)
权威出处：
扩展阅读：
CNKI手机学问
有学问，才够权威！
出版：《中国学术期刊（光盘版）》电子杂志社有限公司
地址：北京清华大学 84-48信箱 大众知识服务
互联网出版许可证 新出网证(京)字008号
京ICP证040431号
服务咨询：400-810--9993
订购咨询：400-819-9993
传真：010-4558人阅读
GIS底层开发（38）
OpenCL（7）
基于GPU的高分一号影像正射校正的设计与实现一 RPC正射校正的原理影像正射校正的方法有很多，主要包含两大类：一类是严格的几何纠正模型，另一类是近似几何纠正模型。当遥感影像的成像模型和有关参数已知时，可以根据严格的成像模型来校正图像，这种方法属于严格几何纠正，最具代表的是共线方程法。当传感器成像模型未知或者无法获取相关的辅助参数时，可以用假定的数学模型模拟成像模型，对影像实现校正，这种方法属于近似几何纠正，主要有：几何多项式纠正、有理函数法、局部区域校正等模型。本文将主要对RPC正射校正模型进行展开讨论。RPC模型将像点坐标d（line，sample）表示为以地面点大地坐标D（Latitude，Longitude，Height）为自变量的比值。为了减小计算过程中舍入误差，增强参数求解的稳定性，需要把地面坐标和影像坐标正则化到（-1，1）之间。对于一个影像，定义如下比值多项式：&&&&&&&&&&&&&&&&&& && &&& &&& &&& &&& &&& &&& &(1-1)式中， &&其中，b1和d1通常为1，（P,L，H）为正则化的地面坐标，（X，Y）为正则化的影像坐标，a1,a2,a3,...a20,b1,b2,b3,...b20,c1,c2,c3,...c20,d1,d2,d3,..d20为RPC参数。正则化规则如下：&这里，LAT_OFF,LAT_SCALE,LONG_OFF,LONG_SCALE,HEIGHT_OFF,HEIGHT_SCALE为物方坐标的正则化参数。SAMP_OFF,SAMP_SCALE,LINE_OFF,LINE_SCALE为像方坐标的正则化参数。它们与RPC模型中4个多项式的80个系数共同保存于卫星厂家提供给用户的RPC文件中。在RPC模型中，可用一阶多项式来表示由光学投影引起的畸变误差模型，由二阶多项式趋近由地球曲率、投影折射、镜头倾斜等因素引起的畸变，可用三阶多项式来模拟高阶部分的其他未知畸变。&& RPC模型的优化方法有两种：分为直接优化法和间接优化法。直接优化法是利用足够的地面控制点信息重新计算RPC系数。因此优化后的RPC系数可以直接代入用于影像的定位和计算，而不需要对已经建立起来的模型进行改变。间接优化法是在像方或物方空间进行一些补偿换算，以减小一些系统误差。在间接优化法中，仿射变换法是最常用的方法。仿射变换公式如下：式中，(x,y )是在影像上量测得到的GCP的影像坐标。(sample,line)是利用RPC模型计算得到的GCP的影像坐标。(e0,e1,e2,f0,f1,f2)是两组坐标之间的仿射变换参数。利用少量的地面控制点即可解算出仿射变换的参数。这样就能利用仿射变换对由RPC模型计算得到的影像坐标解求出更为精确的行列坐标，从而达到优化RPC模型的目的。&本博客暂时先不考虑控制点和DEM，待以后完善。根据RPC系数，RPC正解变换就是从影像的行列号坐标到地面坐标点的变换，在RPC模型中，正解变换需要用到仿射变换六个参数以及RPC系数。RPC反变换是地面到影像行列号坐标的变换，只需要RPC系数就可以。二、RPC正射校正的CPU上实现RPC正射校正需要解决几个基本的问题，一个问题是第一部分提到的RPC正解和反解变换，另外还需要根据RPC系数填充多项式的系数。1、计算各个多项式的项，即计算公式中(1-1)中的乘积项，也就是L,P,H相乘的那部分，其代码如下：void RPCComputeCoeffTerms( double dfLong, double dfLat, double dfHeight, double *padfTerms )
padfTerms[0] = 1.0;
padfTerms[1] = dfL
padfTerms[2] = dfL
padfTerms[3] = dfH
padfTerms[4] = dfLong * dfL
padfTerms[5] = dfLong * dfH
padfTerms[6] = dfLat * dfH
padfTerms[7] = dfLong * dfL
padfTerms[8] = dfLat * dfL
padfTerms[9] = dfHeight * dfH
padfTerms[10] = dfLong * dfLat * dfH
padfTerms[11] = dfLong * dfLong * dfL
padfTerms[12] = dfLong * dfLat * dfL
padfTerms[13] = dfLong * dfHeight * dfH
padfTerms[14] = dfLong * dfLong * dfL
padfTerms[15] = dfLat * dfLat * dfL
padfTerms[16] = dfLat * dfHeight * dfH
padfTerms[17] = dfLong * dfLong * dfH
padfTerms[18] = dfLat * dfLat * dfH
padfTerms[19] = dfHeight * dfHeight * dfH
}2、计算各个RPC项和L,P,H等的乘积和，代码如下：double RPCEvaluateSum( double *padfTerms, double *padfCoefs )
double dfSum = 0.0;
for( i = 0; i & 20; i++ )
dfSum += padfTerms[i] * padfCoefs[i];
return dfS
}3、正解变换，即影像到地面坐标的变换，这个函数需要用到迭代，具体代码如下：void RPCInverseTransform( stRPCInfo *pRPC,double *dbGeoTran,
double dfPixel, double dfLine, double dfHeight,
double *pdfLong, double *pdfLat )
double dfResultX, dfResultY;
//初始值使用放射变换系数求解
dfResultX = dbGeoTran[0] + dbGeoTran[1] * dfPixel + dbGeoTran[2] * dfL
dfResultY = dbGeoTran[3] + dbGeoTran[4] * dfPixel + dbGeoTran[5] * dfL
//开始用正变换的函数进行迭代
double dfPixelDeltaX=0.0, dfPixelDeltaY=0.0;
for( nIter = 0; nIter & 20; nIter++ )
double dfBackPixel, dfBackL
//反算像点坐标，计算阈值
RPCTransform( pRPC, dfResultX, dfResultY, dfHeight,
&dfBackPixel, &dfBackLine );
dfPixelDeltaX = dfBackPixel - dfP
dfPixelDeltaY = dfBackLine - dfL
dfResultX = dfResultX
- dfPixelDeltaX * dbGeoTran[1] - dfPixelDeltaY * dbGeoTran[2];
dfResultY = dfResultY
- dfPixelDeltaX * dbGeoTran[4] - dfPixelDeltaY * dbGeoTran[5];
if( fabs(dfPixelDeltaX) & 0.001
&& fabs(dfPixelDeltaY) & 0.001 )
//printf(&迭代%d次\n&,nIter);
*pdfLong = dfResultX;
*pdfLat = dfResultY;
}4、反解变换，即地面坐标点到影像行列号坐标的变换，可由RPC系数直接求取，具体代码如下：void RPCTransform( stRPCInfo *pRPC,
double dfLong, double dfLat, double dfHeight,
double *pdfPixel, double *pdfLine )
double dfResultX, dfResultY;
double adfTerms[20];
RPCComputeCoeffTerms(
(dfLong - pRPC-&dfLONG_OFF) / pRPC-&dfLONG_SCALE,
(dfLat - pRPC-&dfLAT_OFF) / pRPC-&dfLAT_SCALE,
(dfHeight - pRPC-&dfHEIGHT_OFF) / pRPC-&dfHEIGHT_SCALE,
adfTerms );
dfResultX = RPCEvaluateSum( adfTerms, pRPC-&adfSAMP_NUM_COEFF )
/ RPCEvaluateSum( adfTerms, pRPC-&adfSAMP_DEN_COEFF );
dfResultY = RPCEvaluateSum( adfTerms, pRPC-&adfLINE_NUM_COEFF )
/ RPCEvaluateSum( adfTerms, pRPC-&adfLINE_DEN_COEFF );
*pdfPixel = dfResultX * pRPC-&dfSAMP_SCALE + pRPC-&dfSAMP_OFF;
*pdfLine = dfResultY * pRPC-&dfLINE_SCALE + pRPC-&dfLINE_OFF;
}至此，一些基本的准备做好了，那么下面就开始真正的校正了，RPC校正的大概的流程如下：（1）根据输入未校正的影像得到输出影像的地理范围、地面分辨率大小以及仿射变换的留个参数。（2）在输出的校正影像上，求得每一个像素对应的地理坐标，然后变换到原始影像上的行列号，根据一定的插值算法求得输出新影像上的对应的像素值。（3）将对应的像素值写入到输出文件中。基本流程已经确定后，那么下面就开始讲如何对大数据进行分块，是按照输入影像分块还是输出影像分块，为了方便编程，一般采用输出影像进行分块，具体过程如下：首先将输出影像按照行进行分块，然后计算每一个输出分块所对应的原始影像的行列号范围，再进行重采样操作，最后得到输出影像的像素值。在本博客中，每多少行为一块分块只是本人自己的实验，实际产品和生产中需要动态检测内存等资源信息。还有一个需要注意的是，假设影像是多波段的影像，数据怎么排列，为了效率，我创建的输出影像是按照BIP规格存储，就一个像素的所有波段的数据存储在一起，每一个分块的数据也按照BIP方式读取，写入也按照BIP格式写入。综上所述，CPU上的正射校正可以比容易写出来，影像的IO还是一如既往的用GDAL，具体代码如下：bool ImageWarpRPC(const char * pszSrcFile, const char * pszDstFile, const char * pszDstWKT,
double dResX, double dResY,
const char * pszFormat)
GDALAllRegister();
CPLSetConfigOption(&GDAL_FILENAME_IS_UTF8&, &NO&);
// 打开原始图像
GDALDatasetH hSrcDS = GDALOpen(pszSrcFile, GA_ReadOnly);
if (NULL == hSrcDS)
GDALDataType eDataType = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1));
int nBandCount = GDALGetRasterCount(hSrcDS);
double dbGeonTran[6];
GDALGetGeoTransform(hSrcDS,dbGeonTran);
const char* pszItem = GDALGetMetadataItem(hSrcDS,&INTERLEAVE&,&IMAGE_STRUCTURE&);
//获得原始图像的行数和列数
int nXsize = GDALGetRasterXSize(hSrcDS);
int nYsize = GDALGetRasterYSize(hSrcDS);
// 提取RPC系数
char** pszMetadata = GDALGetMetadata(hSrcDS,&RPC&);
GDALRPCInfo sRPCI
if (pszMetadata != NULL)
GDALExtractRPCInfo( pszMetadata, &sRPCInfo );
stRPCInfo stI
memcpy(&stInfo,&sRPCInfo,sizeof(stInfo));
// 计算输出图像四至范围、大小、仿射变换六参数等信息
double adfGeoTransform[6] = {0};
double adfExtent[4] = {0};
//minx,maxx,miny,maxy
nPixels = 0, nLines = 0;
double dbX[4];
double dbY[4];
//计算四个角点的地理坐标
RPCInverseTransform(&stInfo,dbGeonTran,0,0,100,dbX,dbY);
RPCInverseTransform(&stInfo,dbGeonTran,nXsize,0,100,dbX+1,dbY+1);
RPCInverseTransform(&stInfo,dbGeonTran,nXsize,nYsize,100,dbX+2,dbY+2);
RPCInverseTransform(&stInfo,dbGeonTran,0,nYsize,100,dbX+3,dbY+3);
adfExtent[0] = min(min(dbX[0],dbX[1]) , min(dbX[2],dbX[3]));
adfExtent[1] = max(max(dbX[0],dbX[1]) , max(dbX[2],dbX[3]));
adfExtent[2] = min(min(dbY[0],dbY[1]) , min(dbY[2],dbY[3]));
adfExtent[3] = max(max(dbY[0],dbY[1]) , max(dbY[2],dbY[3]));
int nMinCellSize = min(dbGeonTran[1],fabs(dbGeonTran[5]));
nPixels = ceil( fabs(adfExtent[1] - adfExtent[0]) / dbGeonTran[1] );
nLines = ceil ( fabs(adfExtent[3] - adfExtent[2]) / fabs(dbGeonTran[5]) );
adfGeoTransform[0] = adfExtent[0];
adfGeoTransform[3] = adfExtent[3];
adfGeoTransform[1] = dbGeonTran[1];
adfGeoTransform[5] = dbGeonTran[5];
// 创建输出图像
GDALDriverH hDriver = GDALGetDriverByName( pszFormat );
if (NULL == hDriver)
char **papszOptions = NULL;
papszOptions = CSLSetNameValue( papszOptions, &INTERLEAVE&, &PIXEL& );
GDALDatasetH hDstDS = GDALCreate( hDriver, pszDstFile, nPixels, nLines, nBandCount, eDataType, papszOptions );
if (NULL == hDstDS)
GDALSetProjection( hDstDS,
pszDstWKT);
GDALSetGeoTransform( hDstDS, adfGeoTransform );
//然后是图像重采样
double dfPixel = 0;
double dfLine = 0;
int nFlag = 0;
float dfValue = 0;
CPLErr err = CE_F
int nOffset = 0;
//内存偏移量
int nBandList[4] = {1,2,3,4};
//分块的做法
int nSubHeight = 2000;
int nBlockNum = (nLines+nSubHeight-1)/nSubH
//块的数量
/*double *pLineXs = new double[nPixels*nSubHeight];
double *pLineYs = new double[nPixels*nSubHeight];*/
unsigned short *pValues = new unsigned short[nPixels*nBandCount*nSubHeight];
unsigned short *pSrcValues = NULL;
for (int nIndex = 0; nIndex & nBlockN nIndex ++)
//最后一块特殊处理
int nBottomIndex = (nIndex + 1) *nSubH
if (nBlockNum - 1 == nIndex)
nBottomIndex = nL
//分块的左上角地理坐标
double dbX1 = adfGeoTransform[0] + 0* adfGeoTransform[1]
+ nIndex*nSubHeight
* adfGeoTransform[2];
double dbY1 = adfGeoTransform[3] + 0 * adfGeoTransform[4]
+ nIndex*nSubHeight
* adfGeoTransform[5];
//分块的右下角地理坐标
double dbX2 = adfGeoTransform[0] + nPixels
* adfGeoTransform[1]
+ nBottomIndex
* adfGeoTransform[2];
double dbY2 = adfGeoTransform[3] + nPixels
* adfGeoTransform[4]
+ nBottomIndex
* adfGeoTransform[5];
//分块的右上角地理坐标
double dbX3 = adfGeoTransform[0] + nPixels * adfGeoTransform[1]
+ nIndex*nSubHeight
* adfGeoTransform[2];
double dbY3 = adfGeoTransform[3] + nPixels * adfGeoTransform[4]
+ nIndex*nSubHeight
* adfGeoTransform[5];
//分块的左下角地理坐标
double dbX4 = adfGeoTransform[0] + 0 * adfGeoTransform[1]
+ nBottomIndex
* adfGeoTransform[2];
double dbY4 = adfGeoTransform[3] + 0 * adfGeoTransform[4]
+ nBottomIndex
* adfGeoTransform[5];
//由输出的图像地理坐标系变换到原始的像素坐标系
RPCTransform(&stInfo,dbX1,dbY1,100,&dfPixel,&dfLine);
int nCol1 = (int)(dfPixel + 0.5);
int nRow1 = (int)(dfLine + 0.5);
RPCTransform(&stInfo,dbX2,dbY2,100,&dfPixel,&dfLine);
int nCol2 = (int)(dfPixel + 0.5);
int nRow2 = (int)(dfLine + 0.5);
RPCTransform(&stInfo,dbX3,dbY3,100,&dfPixel,&dfLine);
int nCol3 = (int)(dfPixel + 0.5);
int nRow3 = (int)(dfLine + 0.5);
RPCTransform(&stInfo,dbX4,dbY4,100,&dfPixel,&dfLine);
int nCol4 = (int)(dfPixel + 0.5);
int nRow4 = (int)(dfLine + 0.5);
int nMinRow = min( min( nRow1,nRow2 ), min( nRow3 , nRow4));
if (nMinRow & 0)
nMinRow = 0;
int nMaxRow = max( max( nRow1,nRow2 ), max( nRow3 , nRow4));
if (nMaxRow &= nYsize)
nMaxRow = nYsize-1;
int nHeight = nMaxRow-nMinRow+1;
//读取原始数据
pSrcValues = new unsigned short[nXsize*nHeight*nBandCount];
int nDataSize = GDALGetDataTypeSize(eDataType)/8;
((GDALDataset*)hSrcDS)-&RasterIO(GF_Read,0,nMinRow,nXsize,nHeight,
pSrcValues,nXsize,nHeight,eDataType,nBandCount,nBandList,nBandCount*nDataSize,
nXsize*nBandCount*nDataSize,1*nDataSize);
//这一块适合做并行
if (HasCUDA())
DWORD t1 = GetTickCount();
/*RPCWarpCL(stInfo,adfGeoTransform,pSrcValues,nXsize,nHeight,nMinRow,
nYsize,pValues,nPixels,nSubHeight,nIndex*nSubHeight,nBandCount,nIndex,nBlockNum);*/
RPCWarpCuda(stInfo,adfGeoTransform,pSrcValues,nXsize,nHeight,nMinRow,
nYsize,pValues,nPixels,nSubHeight,nIndex*nSubHeight,nBandCount,nIndex,nBlockNum);
DWORD t2 = GetTickCount();
double tt = (t2-t1);
printf(&%f毫秒\n&,tt);
for (int nRow = nIndex*nSubH nRow & nBottomI nRow ++)
for (int nCol = 0; nCol & nP nCol ++)
nOffset = ( (nRow - nIndex*nSubHeight)*nPixels + nCol ) * nBandC
double dbX = adfGeoTransform[0] + nCol * adfGeoTransform[1]
* adfGeoTransform[2];
double dbY = adfGeoTransform[3] + nCol * adfGeoTransform[4]
* adfGeoTransform[5];
RPCTransform(&stInfo,dbX,dbY,100,&dfPixel,&dfLine);
int nColIndex = (int)(dfPixel + 0.5);
int nRowIndex = (int)(dfLine + 0.5);
int nOffsetSrc = ( (nRowIndex - nMinRow)*nXsize + nColIndex ) * nBandC
//超出范围的用0填充
if (nColIndex & 0 || nColIndex &= nXsize || nRowIndex & 0 || nRowIndex &= nYsize)
for (int nBandIndex = 0; nBandIndex & nBandC nBandIndex ++)
pValues[nOffset++] = 0;
for (int nBandIndex = 0; nBandIndex & nBandC nBandIndex ++)
pValues[nOffset++] = pSrcValues[nOffsetSrc++];
delete []pSrcV
//写入数据
if (nIndex & nBlockNum - 1)
( (GDALDataset*)hDstDS )-&RasterIO(GF_Write,0,nIndex*nSubHeight,nPixels,nSubHeight,
pValues,nPixels,nSubHeight,(GDALDataType)eDataType,nBandCount,nBandList,nBandCount*nDataSize,
nPixels*nBandCount*nDataSize,1*nDataSize);
else if (nIndex == nBlockNum - 1)
int nRealIndex = nIndex*nSubH
nSubHeight = nLines - nIndex*nSubH
( (GDALDataset*)hDstDS )-&RasterIO(GF_Write,0,nRealIndex,nPixels,nSubHeight,
pValues,nPixels,nSubHeight,(GDALDataType)eDataType,nBandCount,nBandList,nBandCount*nDataSize,
nPixels*nBandCount*nDataSize,1*nDataSize);
/*delete []pLineXs;
delete []pLineYs;*/
delete []pV
GDALClose(hSrcDS );
GDALClose(hDstDS );
}上面的代码还是有很多瑕疵的，比如计算输出影像的范围，比较保险的做法是根据原始影像的四条边来计算，还有影像的数据类型还不支持泛型，这些都将在以后不断完善。细心的读者可能已经发现有RPCWarpCuda和RPCWrapCL的函数，不错，这正是CUDA实现和opencl上的实现，其实大家都知道影像几何校正的时间大部分都花费在影像的重采样上，所以本博客对于正射校正的GPU实现也主要针对坐标的正反算和重采样上。下面开始GPU之旅吧！&三、RPC正射校正的GPU实现自从NVIDIA推出CUDA，GPGPU的发展速度非常迅速，以前要编写GPU上处理非图形的编程，需要将问题域转化为图形问题，而GPGPU的出现，大大简化了GPU编程，随着技术的进一步发展，之后又出现opencl，这是跨平台的一个标准，最先由苹果公司联合几家公司推出的标准，opencl的具体实现交给各个硬件厂商。下面就着重介绍CUDA上的实现。CUDA上实现的函数声明如下：extern &C& void RPCWarpCuda(stRPCInfo& stInfo,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& double*pfGeoTransform,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& unsignedshort* poDataIn,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& intnWidthIn,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& intnHeightIn,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& intnMinRowIn,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& intnSrcHeight,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& unsignedshort* poDataOut,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& intnWidthOut,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& intnHeightOut,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& intnMinRowOut,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& intnBandCount,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& intnBlockIndex,&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& intnBlockNum);&CUDA的运行模型是由多个线程运行同一个计算过程，在一个计算设备上，划分为不同的线程块，每一个线程块由不同的线程组成，这样组成了一个两层的线程组织体系，对于图像处理来说，最适合划分为二维网格，每一个线程计算一个输出像素，在RPC正射校正模型中，波段数、RPC的参数以及归一化的参数以及输出影像的仿射变换系数在每一个线程都是一样的，所以只需拷贝这些参数到GPU上一次，这些参数在CUDA设备上需要声明为__constant__的，如下：__constant__ double dbGeoTrans[6];&__constant__&&& double&&&&& dfLINE_OFF;__constant__&&& double&&&&& dfSAMP_OFF;__constant__&&& double&&&&& dfLAT_OFF;__constant__&&& double&&&&& dfLONG_OFF;__constant__&&& double&&&&& dfHEIGHT_OFF;&&&& //缩放比例__constant__&&& double&&&&& dfLINE_SCALE;__constant__&&& double&&&&& dfSAMP_SCALE;__constant__&&& double&&&&& dfLAT_SCALE;__constant__&&& double&&&&& dfLONG_SCALE;__constant__&&& double&&&&& dfHEIGHT_SCALE;&&&& //系数__constant__&&& double&&&&& adfLINE_NUM_COEFF[20];__constant__&&& double&&&&& adfLINE_DEN_COEFF[20];__constant__&&& double&&&&& adfSAMP_NUM_COEFF[20];__constant__&&& double&&&&& adfSAMP_DEN_COEFF[20];&//波段数__constant__ int nBandCount;&在RPCWarpCuda中，只在第一个块的数据传递过来的时候将这些数据拷贝的CUDA设备的常量内存上。代码如下，if(0 == nBlockIndex)
cudaMalloc(&poDataOut_d,nDataSizeOut);
//给参数常量分配空间
cudaMemcpyToSymbol(adfLINE_DEN_COEFF,stInfo.adfLINE_DEN_COEFF,sizeof(double)*20);
cudaMemcpyToSymbol(adfLINE_NUM_COEFF,stInfo.adfLINE_NUM_COEFF,sizeof(double)*20);
cudaMemcpyToSymbol(adfSAMP_DEN_COEFF,stInfo.adfSAMP_DEN_COEFF,sizeof(double)*20);
cudaMemcpyToSymbol(adfSAMP_NUM_COEFF,stInfo.adfSAMP_NUM_COEFF,sizeof(double)*20);
cudaMemcpyToSymbol(dfHEIGHT_OFF,&stInfo.dfHEIGHT_OFF,sizeof(double));
cudaMemcpyToSymbol(dfHEIGHT_SCALE,&stInfo.dfHEIGHT_SCALE,sizeof(double));
cudaMemcpyToSymbol(dfLAT_OFF,&stInfo.dfLAT_OFF,sizeof(double));
cudaMemcpyToSymbol(dfLAT_SCALE,&stInfo.dfLAT_SCALE,sizeof(double));
cudaMemcpyToSymbol(dfLINE_OFF,&stInfo.dfLINE_OFF,sizeof(double));
cudaMemcpyToSymbol(dfLINE_SCALE,&stInfo.dfLINE_SCALE,sizeof(double));
cudaMemcpyToSymbol(dfLONG_OFF,&stInfo.dfLONG_OFF,sizeof(double));
cudaMemcpyToSymbol(dfLONG_SCALE,&stInfo.dfLONG_SCALE,sizeof(double));
cudaMemcpyToSymbol(dfSAMP_OFF,&stInfo.dfSAMP_OFF,sizeof(double));
cudaMemcpyToSymbol(dfSAMP_SCALE,&stInfo.dfSAMP_SCALE,sizeof(double));
//传递放射变换系数
cudaMemcpyToSymbol(dbGeoTrans,pfGeoTransform,sizeof(double)*6);
cudaMemcpyToSymbol(nBandCount,&nBands,sizeof(int));
}由于输出影像是均匀分块，所以只需要第一次给输出的数据申请设备存在，这样可以减少设备内存分配的时间，在最后一次的时候将输出数据设备内存释放掉。在RPCWarpCuda中，只在第一个块的数据传递过来的时候将这些数据拷贝的CUDA设备的常量内存上。还有调用内很函数的时候将线程块设为32*32的，这刚好是1024个线程，现在很多设备都至少支持每个线程块都支持1024个线程，但是真正产品中需要检测设备的参数。一切准备工作就绪后，那么下面就只需编写内核函数了，其实内核函数要做的工作主要是影像坐标和地面点的坐标之间的换算以及影像重采样，其核函数如下：&__global__ void RPCWarpKernel(
unsigned short* poDataIn,
int nWidthIn,
int nHeightIn,
int nMinRowIn,
int nSrcHeight,
unsigned short* poDataOut,
int nWidthOut,
int nHeightOut,
int nMinRowOut)
int idy = blockIdx.y*blockDim.y + threadIdx.y; //行
int idx = blockIdx.x*blockDim.x + threadIdx.x; //列
if (idx & nWidthOut && idy & nHeightOut)
//求出输出影像的实际行列号
int nRow = idy + nMinRowO
int nCol =
double dbX = dbGeoTrans[0] + nCol * dbGeoTrans[1]
* dbGeoTrans[2];
double dbY = dbGeoTrans[3] + nCol * dbGeoTrans[4]
* dbGeoTrans[5];
double dfPixel = 0.0;
double dfLine = 0.0;
RPCTransformCUDA(dbX,dbY,100,&dfPixel,&dfLine);
//求出原始影像所在坐标，然后最邻近采样
int nColIndex = (int)(dfPixel + 0.5);
int nRowIndex = (int)(dfLine + 0.5);
int nOffsetSrc = ( (nRowIndex - nMinRowIn)*nWidthIn + nColIndex ) * nBandC
int nOffset = ( (nRow - nMinRowOut)*nWidthOut + nCol ) * nBandC
//超出范围的用0填充
if (nColIndex & 0 || nColIndex &= nWidthIn || nRowIndex & 0 || nRowIndex &= nSrcHeight)
for (int nBandIndex = 0; nBandIndex & nBandC nBandIndex ++)
poDataOut[nOffset++] = 0;
for (int nBandIndex = 0; nBandIndex & nBandC nBandIndex ++)
poDataOut[nOffset++] = poDataIn[nOffsetSrc++];
}至此，CUDA上的实现就介绍到这里，关于OPENCL上的实现也和CUDA实现差不多，具体就不阐述了，我会上传代码，有兴趣的读者可以下载自己研究并改进。四、性能测试&&& 我测试环境是硬件环境是GT750m，CUDA是5.5版本，opencl是CUDA上的实现。测试数据时高分一号宽幅的数据，行数是13400，列数是12000,4个波段，数据类型是unsigned short，原始影像文件的大小是1.2GB。经过大量的测试，测试的结果如下：
包括IO时间
不包括IO时间
19.26 从上面的结论可以看出，如果不统计影像IO的时间，CUDA和OPENCL都可以获得很大的加速比，加速比能够达到20左右，并且CUDA和OPENCL的加速效果相当接近。影像IO的时间大约是80-100秒左右。如果计算影像IO的时间，加速效果大概是2.68倍左右，所以下一步该研究如何优化影像的IO操作。校正前后的图像比较如下：校正前校正后。&&五、结论与展望通过对RPC正射校正的的GPU实现，可以看出，GPU在遥感影像处理中可以发挥它巨大的优势，有GPU的用武之地。但是本文没对GPU优化进行进一步研究。下一步研究加入DEM以及加入控制点提高精度；对于影像IO部分，也需要深入的研究。代码下载地址是：对于相关的环境，可能需要读者自己配置了。也欢迎大家提出意见，共同进步。
参考知识库
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：577181次
积分：6523
积分：6523
排名：第2590名
原创：102篇
转载：20篇
评论：331条
毕业于FuZhou University，目前在北京从事GIS与遥感的研发工作，研究兴趣为GIS，并专注于C/C++。如果有什么事可以通过邮箱联系到我，如果需要转载原创的文章，也请添加链接，谢谢！
文章：34篇
阅读：203765
(1)(1)(1)(1)(2)(8)(1)(1)(2)(3)(2)(2)(3)(2)(3)(1)(2)(1)(5)(1)(4)(5)(1)(2)(1)(1)(4)(5)(5)(6)(5)(10)(5)(7)(9)(1)(2)(6)