OSGB格式三维模型投影坐标转换方法

程晓光，王婉秋，陆泉源

(飞燕航空遥感技术有限公司，南京 210001)

目前，基于倾斜摄影测量技术生产的三维模型已成为主要的测绘产品之一，其生产步骤包括密集点云提取、网格面模型构建、纹理映射、成果编辑、成果导出等[1]，较为复杂。ContextCapture[2]是目前倾斜摄影测量领域主要的三维建模软件之一[1, 3]，而OpenSceneGraph Binary(OSGB)是ContextCapture导出三维模型时最常用的格式之一，由开源三维图形工具包OpenSceneGraph[4](OSG)推出，属于二进制文件格式，带有嵌入式链接纹理数据，允许以细节层次(Level Of Detail，LOD)技术进行多级金字塔显示[4]。

三维模型使用的空间参考系统以投影坐标系居多。根据国家标准[5-8]，在基本比例尺地形图和正射影像地图中，大地基准采用2000国家大地坐标系，高程基准主要采用1985国家高程基准，除1∶1 000 000外的其它比例尺的投影方法均为3°或6°分带的高斯-克吕格投影。此外，各大城市和地区还建立了地方独立坐标系，目的是将高程归化与投影变形的影响控制在微小范围内，保证城市范围内长度变形最小，以满足当地规划建设需要[9]。大多数地方独立坐标系是在1954年北京坐标系或1980年西安坐标系的基础上利用高斯-克吕格投影进行改化后建立的。常见改化包括将中央子午线移至测区中央，将投影面改为测区平均高程面、对原点纵横坐标进行偏移等[9]。在三维模型交付时，除了国家标准投影坐标系外，还经常需要以地方独立坐标系交付。

使用三维建模软件在不同投影坐标系下各生产一遍三维模型会造成巨大的工作量，且模型间可能存在偏差。如果ContextCapture在输出模型时选择的投影坐标系与空三使用的坐标系基于不同椭球体，则输出模型的坐标易有较大误差。有商业软件可对三维模型进行投影坐标转换，但转换后的模型经常存在纹理丢失、显示闪烁等问题。

针对ContextCapture导出的OSGB格式的三维模型的投影坐标转换问题，文中分析OSGB格式的三维模型的目录和文件结构，提出一种简单高效的投影坐标转换方法，先使用OSG获得待转换坐标，再进行坐标转换，最后利用文件坐标数据替换法保存三维模型，并给出示例代码。该方法扩展性强，转换后的模型无纹理丢失、无变形、无闪烁。

1 OSGB格式模型的目录和文件结构

在ContextCapture的生产项目定义中设置输出格式为OSGB，划分瓦片(Tile)，使用LOD技术，则输出的三维模型存储文件夹的结构见图1。Data文件夹下为各瓦片的子文件夹，每个子文件夹存储该瓦片不同金字塔级别的OSGB文件，其中一个OSGB文件和子文件夹同名。但是该模式并不便于模型浏览。为此，一种常见的做法是在三维模型存储文件夹下，额外放置S3C文件以对构成场景的瓦片文件进行索引，方便浏览显示。

图1 OSGB格式的三维模型的目录结构

为降低内存和磁盘使用，OSGB文件常采用4字节长度(单精度)而非8字节长度(双精度)的浮点数来存储顶点坐标。对中国境内的高斯-克吕格投影来说，如果坐标以m为单位，则整数部分常会有6位或7位，但单精度浮点数只能保证7位十进制有效数字[10]。如果在OSGB文件中存储顶点原始投影坐标，则容易造成坐标小数部分有效位数不够。为此，在采用投影坐标系的情况下，OSGB文件一般对顶点坐标做偏移，以保证小数点后有足够的有效数字。坐标偏移量存储在metadata.xml的SRSOrigin标签中。在metadata.xml中，投影坐标系信息存储在SRS标签，支持欧洲石油调查组织[11](European Petroleum Survey Group，EPSG)定义的坐标系编码或众知文本(Well-Known Text，WKT)。

单个OSGB文件存储金字塔某级在特定空间范围内的模型，其基于OSG体系的典型组织结构见图2。下文如无特殊说明，则OSG类的默认命名空间是osg。位于金字塔最底层的模型的根节点为一个Geode类对象；利用四叉树或八叉树等[12]技术做模型分割的根节点为一个Group类对象；除最底层外，未做分割的模型的根节点为一个PagedLOD类对象。这其中涉及到的主要节点类之间的继承关系见图3。

图2 典型OSGB文件的组织结构

图3 主要节点类之间的继承关系(箭头起点为父类，终点为子类)

2 投影坐标转换方法

对OSGB格式的三维模型进行投影坐标转换，可以先对模型顶点做坐标转换，再重新映射纹理，最后保存模型到文件。但该方法非常复杂，实现难度大。最简单直观的思路是直接将OSGB文件中的投影坐标替换为转换并偏移后的坐标，并修改metadata.xml和S3C文件中的坐标系及坐标偏移量信息。

2.1 待转换坐标获取

在Geometry类对象的属性中，一般只有顶点数组含有需要转换的坐标数据。PagedLOD类对象根据视点到物体中心的距离决定显示哪种精细度的模型，而确定物体中心的方法由属性“中心模式”给出，既可以是包围球中心，也可以是用户自定义中心，还可以是二者的结合。对于后两种模式，需要对用户自定义中心的坐标进行转换，否则可能会造成模型显示时被误裁剪。文中定义三维模型需要转换的坐标包括Geometry类对象的顶点坐标和PagedLOD类对象的用户自定义中心坐标。

(1)

2.2 投影坐标转换

将投影坐标转换分为相同椭球体的投影坐标转换和不同椭球体间的投影坐标转换。前者的转换框架可见图4。后者采用七参数法[9]进行转换，其中未假设3个旋转角很小，转换框架可见图5。转换前后高程坐标不变。

图4 相同椭球体的投影坐标转换框架

图5 不同椭球体间的投影坐标转换框架

2.3 转换模型保存

(2)

在坐标转换和偏移完成后，需要将转换后模型保存到OSGB文件。理论上，可以调用osgDB::Registry类实例的writeNode函数将模型保存为OSGB文件。但这样做存在格式兼容性风险。如果ContextCapture不能支持所有版本的OSGB格式，则该方法生成的OSGB文件可能不能被ContextCapture准确读取。

一种简单的方法是直接对原始OSGB文件进行坐标数据替换。在OSGB文件的数据流中，一个Geometry类对象的顶点坐标连续放置，可被视为一个浮点型数组，而PagedLOD类对象的用户自定义中心坐标可被视为一个具有3个元素的浮点型数组。设转换前的数组为Asrc，存储了N个三维坐标，元素数目为3N，对应的转换和偏移后的数组为Adst，则以字节为单位的Asrc和Adst的大小nLength均为12N(单精度浮点型)或24N(双精度浮点型)。通过在文件数据流中比较长度为nLength的连续内存的值是否和Asrc相同来确定数据替换位置，在找到替换位置后，用Adst的值覆盖原始数据。

3 实现方法

利用C++和开源软件开发OSGB模型的投影坐标转换软件，其中Qt[13]用作编程框架，OpenSceneGraph 3.7[4]用于三维模型读取，PROJ 6.0[14]用于地图投影，GDAL 3.0[15]用于坐标系设置，Eigen3[16]用于矩阵运算。

3.1 坐标转换流程

图 6为设计的坐标转换流程。

图6 模型投影坐标转换流程

3.2 转换前模型解析

3.2.1 源坐标系解析

如果模型源坐标系以EPSG编码的形式提供，可调用OGRSpatialReference类的importFromEPSG函数导入编码对应的坐标系；如果模型源坐标系以WKT文本的形式提供，可调用OGRSpatialReference类的importFromWkt函数导入WKT对应的坐标系。

3.2.2 OSGB文件解析

要读取OSGB文件中的模型，可先调用osgDB::readNodeFile函数获得Node类型的模型对象。之后，要获取待转换坐标，需要对OSG定义的类NodeVisitor进行派生并重写相应的apply虚函数。示例代码如下：

class CTVisitor: public NodeVisitor {

public:

virtual void apply(Geode &node);

virtual void apply(Group &node);

}

对Geometry类对象的解析可放在void CTVisitor::apply(Geode &node)中，示例如下：

for(int n = 0; n < node.getNumDrawables(); n++){

Drawable *pDb = node.getDrawable(n);

if(pDb->className()== string("Geometry")){

ref_ptr pGm =(Geometry*)pDb;

//pVt为顶点坐标数组指针

Vec3Array *pVt =dynamic_cast

(pGm->getVertexArray().get());

}

}

traverse(node);

如果模型采用双精度浮点数存储顶点坐标，则需要将上述代码中的Vec3Array替换为Vec3dArray。对PagedLOD类对象的解析可放在void CTVisitor::apply(Group &node)中，示例如下：

if(node.className()== string("PagedLOD")){

PagedLOD lod =(PagedLOD)node;

if((lod.getCenterMode()== LOD::USER_ DEFINED_CENTER)||(lod.getCenterMode()== LOD::UNION_OF_BOUNDING_SPHERE_AND _USER_DEFINED)){

//cent为用户自定义中心坐标

Vec3 cent = lod.getCenter();

}

}

traverse(node);

3.3 投影坐标转换

3.3.1 大地坐标与投影坐标间的转换

大地坐标与投影坐标间的转换采用pj_transform函数，该函数定义如下：

int pj_transform(projPJ src, projPJ dst, long pt_count, int pt_offset, double *x, double *y, double *z)

其中，src为源坐标系对象，dst为目标坐标系对象，pt_count为待转换点数，pt_offset为首个待转换点在数组中的序号，x，y，z分别为存储三维坐标的数组的地址。在PROJ中，大地和投影坐标系对象需用PROJ字符串定义。大地坐标系的PROJ字符串为：

QString("+proj=longlat +ellps=%1 +no_defs +type= crs").arg(strEllps)

高斯-克吕格投影坐标系的PROJ字符串为：

QString("+proj=tmerc +lat_0=0 +lon_0=%1 +k=1 +x_0=%2 +y_0=%3 +ellps=%4 +units=m +no_defs +type=crs").arg(dCentLon).arg(dFalseEast).arg(dFalseNorth).arg(strEllps)

其中，strEllps为大地坐标系使用的椭球体的字符串，1954年北京坐标系是“krass”，1980年西安坐标系是“IAU76”，2000国家大地坐标系是“GRS80”；dCentLon为中央子午线以角度制表示的经度，dFalseEast为以m为单位的东偏，dFalseNorth为以m为单位的北偏。之后，调用pj_init_plus函数，输入参数为坐标系的PROJ字符串，即可获得projPJ类型的坐标系对象。

3.3.2 大地坐标与空间直角坐标间的转换

从大地坐标转换到空间直角坐标采用pj_geodetic_to_geocentric函数，该函数定义如下：

int pj_geodetic_to_geocentric(double a, double es, long pt_count, int pt_offset, double *x, double *y, double *z)

从空间直角坐标转换到大地坐标采用pj_geocentric_to_geodetic函数，该函数定义如下：

int pj_geocentric_to_geodetic(double a, double es, long pt_count, int pt_offset, double *x, double *y, double *z)

其中，a为椭球体半长轴，es为第一偏心率平方，pt_count为待转换点数，pt_offset为首个待转换点在数组中的序号，x，y，z分别为存储三维坐标的数组的地址。

3.4 转换后模型生成

在生成转换后模型时，需要确保转换前后目录结构不变，OSGB文件名和相对路径不变。

3.4.1 metadata.xml生成

在转换后模型的metadata.xml中，对于EPSG已定义的投影坐标系，可将SRS标签值设置为坐标系的EPSG编码；对于EPSG未定义的投影坐标系，可先调用OGRSpatialReference类的SetWellKnownGeogCS函数设置大地坐标系，再调用SetTM函数设置中央子午线经度、中央纬线纬度、中央子午线缩放因子、东偏、北偏等参数，调用SetProjCS函数设置坐标系名称，调用exportToWkt函数获得坐标系的WKT，最后将SRS标签值设置为WKT。SRSOrigin标签值为Δxdst，Δydst,Δzdst。

3.4.2 OSGB文件生成

设Asrc的数组指针为pSrc，Adst的数组指针为pDst，以字节为单位的OSGB文件大小为nOsgbSize。先将转换前的OSGB文件全部读入unsigned char类型、大小为nOsgbSize的数组中。设数组指针为pszOsgb，则坐标数据替换的示例代码如下：

for(int n = 0; n < nOsgbSize-nLength; n++){

int b = memcmp(pszOsgb + n, pSrc, nLength);

if(b == 0){

memcpy(pszOsgb + n, pDst, nLength);

break;

}

}

最后将pszOsgb内的数据保存到转换后目录下的OSGB文件。

3.4.3 S3C文件坐标系设置

虽然S3C文件也存储坐标系信息，但由于新版S3C文件的格式已做加密，难以用编程的方法修改其坐标系信息。因此，先将转换前的S3C文件复制到转换后目录下，再利用CC_S3CComposer.exe将S3C文件中的坐标系信息替换为metadata.xml中的坐标系信息。

4 软件实例

图7是软件界面，支持中国常用大地坐标系基础上的高斯-克吕格投影；支持自定义和自动计算输出坐标偏移量；支持手动和文件输入椭球体转换七参数。

图7 软件界面

在Intel Core i5-2300的CPU、24GB内存的台式机上，软件将具有6个瓦片、1179个OSGB文件的三维模型从2000国家大地坐标系基础上120°中央子午线的投影坐标系转换到1954年北京大地坐标系基础上120.75°中央子午线的投影坐标系，耗时20.214 s，平均单个OSGB文件耗时17 ms。在利用CC_S3CComposer.exe更新S3C文件的坐标系信息后，可以用ContextCapture Viewer[1]打开转换后的模型，几何无变形，纹理无缺失，显示无闪烁。为了对比，由人工在目标坐标系下重新生产三维模型。通过在软件转换后模型和人工重新生产的模型上选取同名特征点并查询其坐标，发现坐标一致。

5 结束语

文中提出的投影坐标转换方法简单高效，对OSGB文件的格式要求不高，不受是否采用LOD技术、有无纹理贴图、纹理压缩比例、切块方式等因素的影响。由于OSG基本支持所有主流的三维模型格式，所以文中方法容易扩展到其它二进制的三维模型格式。此外，坐标数据替换法也适用于其它需要对模型坐标进行修改的情况。

文中方法未修改Geometry类对象的法向量，但理论上，修改顶点坐标之后需要重新计算法向量。虽然在第4节例子中不更新法向量在视觉上没有明显影响，但更新模型法向量仍是一个改进方向。