OSGB格式三維模型投影坐標(biāo)轉(zhuǎn)換方法

程曉光，王婉秋，陸泉源

(飛燕航空遙感技術(shù)有限公司，南京 210001)

目前，基于傾斜攝影測(cè)量技術(shù)生產(chǎn)的三維模型已成為主要的測(cè)繪產(chǎn)品之一，其生產(chǎn)步驟包括密集點(diǎn)云提取、網(wǎng)格面模型構(gòu)建、紋理映射、成果編輯、成果導(dǎo)出等[1]，較為復(fù)雜。ContextCapture[2]是目前傾斜攝影測(cè)量領(lǐng)域主要的三維建模軟件之一[1, 3]，而OpenSceneGraph Binary(OSGB)是ContextCapture導(dǎo)出三維模型時(shí)最常用的格式之一，由開源三維圖形工具包OpenSceneGraph[4](OSG)推出，屬于二進(jìn)制文件格式，帶有嵌入式鏈接紋理數(shù)據(jù)，允許以細(xì)節(jié)層次(Level Of Detail，LOD)技術(shù)進(jìn)行多級(jí)金字塔顯示[4]。

三維模型使用的空間參考系統(tǒng)以投影坐標(biāo)系居多。根據(jù)國家標(biāo)準(zhǔn)[5-8]，在基本比例尺地形圖和正射影像地圖中，大地基準(zhǔn)采用2000國家大地坐標(biāo)系，高程基準(zhǔn)主要采用1985國家高程基準(zhǔn)，除1∶1 000 000外的其它比例尺的投影方法均為3°或6°分帶的高斯-克呂格投影。此外，各大城市和地區(qū)還建立了地方獨(dú)立坐標(biāo)系，目的是將高程歸化與投影變形的影響控制在微小范圍內(nèi)，保證城市范圍內(nèi)長度變形最小，以滿足當(dāng)?shù)匾?guī)劃建設(shè)需要[9]。大多數(shù)地方獨(dú)立坐標(biāo)系是在1954年北京坐標(biāo)系或1980年西安坐標(biāo)系的基礎(chǔ)上利用高斯-克呂格投影進(jìn)行改化后建立的。常見改化包括將中央子午線移至測(cè)區(qū)中央，將投影面改為測(cè)區(qū)平均高程面、對(duì)原點(diǎn)縱橫坐標(biāo)進(jìn)行偏移等[9]。在三維模型交付時(shí)，除了國家標(biāo)準(zhǔn)投影坐標(biāo)系外，還經(jīng)常需要以地方獨(dú)立坐標(biāo)系交付。

使用三維建模軟件在不同投影坐標(biāo)系下各生產(chǎn)一遍三維模型會(huì)造成巨大的工作量，且模型間可能存在偏差。如果ContextCapture在輸出模型時(shí)選擇的投影坐標(biāo)系與空三使用的坐標(biāo)系基于不同橢球體，則輸出模型的坐標(biāo)易有較大誤差。有商業(yè)軟件可對(duì)三維模型進(jìn)行投影坐標(biāo)轉(zhuǎn)換，但轉(zhuǎn)換后的模型經(jīng)常存在紋理丟失、顯示閃爍等問題。

針對(duì)ContextCapture導(dǎo)出的OSGB格式的三維模型的投影坐標(biāo)轉(zhuǎn)換問題，文中分析OSGB格式的三維模型的目錄和文件結(jié)構(gòu)，提出一種簡(jiǎn)單高效的投影坐標(biāo)轉(zhuǎn)換方法，先使用OSG獲得待轉(zhuǎn)換坐標(biāo)，再進(jìn)行坐標(biāo)轉(zhuǎn)換，最后利用文件坐標(biāo)數(shù)據(jù)替換法保存三維模型，并給出示例代碼。該方法擴(kuò)展性強(qiáng)，轉(zhuǎn)換后的模型無紋理丟失、無變形、無閃爍。

1 OSGB格式模型的目錄和文件結(jié)構(gòu)

在ContextCapture的生產(chǎn)項(xiàng)目定義中設(shè)置輸出格式為OSGB，劃分瓦片(Tile)，使用LOD技術(shù)，則輸出的三維模型存儲(chǔ)文件夾的結(jié)構(gòu)見圖1。Data文件夾下為各瓦片的子文件夾，每個(gè)子文件夾存儲(chǔ)該瓦片不同金字塔級(jí)別的OSGB文件，其中一個(gè)OSGB文件和子文件夾同名。但是該模式并不便于模型瀏覽。為此，一種常見的做法是在三維模型存儲(chǔ)文件夾下，額外放置S3C文件以對(duì)構(gòu)成場(chǎng)景的瓦片文件進(jìn)行索引，方便瀏覽顯示。

圖1 OSGB格式的三維模型的目錄結(jié)構(gòu)

為降低內(nèi)存和磁盤使用，OSGB文件常采用4字節(jié)長度(單精度)而非8字節(jié)長度(雙精度)的浮點(diǎn)數(shù)來存儲(chǔ)頂點(diǎn)坐標(biāo)。對(duì)中國境內(nèi)的高斯-克呂格投影來說，如果坐標(biāo)以m為單位，則整數(shù)部分常會(huì)有6位或7位，但單精度浮點(diǎn)數(shù)只能保證7位十進(jìn)制有效數(shù)字[10]。如果在OSGB文件中存儲(chǔ)頂點(diǎn)原始投影坐標(biāo)，則容易造成坐標(biāo)小數(shù)部分有效位數(shù)不夠。為此，在采用投影坐標(biāo)系的情況下，OSGB文件一般對(duì)頂點(diǎn)坐標(biāo)做偏移，以保證小數(shù)點(diǎn)后有足夠的有效數(shù)字。坐標(biāo)偏移量存儲(chǔ)在metadata.xml的SRSOrigin標(biāo)簽中。在metadata.xml中，投影坐標(biāo)系信息存儲(chǔ)在SRS標(biāo)簽，支持歐洲石油調(diào)查組織[11](European Petroleum Survey Group，EPSG)定義的坐標(biāo)系編碼或眾知文本(Well-Known Text，WKT)。

單個(gè)OSGB文件存儲(chǔ)金字塔某級(jí)在特定空間范圍內(nèi)的模型，其基于OSG體系的典型組織結(jié)構(gòu)見圖2。下文如無特殊說明，則OSG類的默認(rèn)命名空間是osg。位于金字塔最底層的模型的根節(jié)點(diǎn)為一個(gè)Geode類對(duì)象；利用四叉樹或八叉樹等[12]技術(shù)做模型分割的根節(jié)點(diǎn)為一個(gè)Group類對(duì)象；除最底層外，未做分割的模型的根節(jié)點(diǎn)為一個(gè)PagedLOD類對(duì)象。這其中涉及到的主要節(jié)點(diǎn)類之間的繼承關(guān)系見圖3。

圖2 典型OSGB文件的組織結(jié)構(gòu)

圖3 主要節(jié)點(diǎn)類之間的繼承關(guān)系(箭頭起點(diǎn)為父類，終點(diǎn)為子類)

2 投影坐標(biāo)轉(zhuǎn)換方法

對(duì)OSGB格式的三維模型進(jìn)行投影坐標(biāo)轉(zhuǎn)換，可以先對(duì)模型頂點(diǎn)做坐標(biāo)轉(zhuǎn)換，再重新映射紋理，最后保存模型到文件。但該方法非常復(fù)雜，實(shí)現(xiàn)難度大。最簡(jiǎn)單直觀的思路是直接將OSGB文件中的投影坐標(biāo)替換為轉(zhuǎn)換并偏移后的坐標(biāo)，并修改metadata.xml和S3C文件中的坐標(biāo)系及坐標(biāo)偏移量信息。

2.1 待轉(zhuǎn)換坐標(biāo)獲取

在Geometry類對(duì)象的屬性中，一般只有頂點(diǎn)數(shù)組含有需要轉(zhuǎn)換的坐標(biāo)數(shù)據(jù)。PagedLOD類對(duì)象根據(jù)視點(diǎn)到物體中心的距離決定顯示哪種精細(xì)度的模型，而確定物體中心的方法由屬性“中心模式”給出，既可以是包圍球中心，也可以是用戶自定義中心，還可以是二者的結(jié)合。對(duì)于后兩種模式，需要對(duì)用戶自定義中心的坐標(biāo)進(jìn)行轉(zhuǎn)換，否則可能會(huì)造成模型顯示時(shí)被誤裁剪。文中定義三維模型需要轉(zhuǎn)換的坐標(biāo)包括Geometry類對(duì)象的頂點(diǎn)坐標(biāo)和PagedLOD類對(duì)象的用戶自定義中心坐標(biāo)。

(1)

2.2 投影坐標(biāo)轉(zhuǎn)換

將投影坐標(biāo)轉(zhuǎn)換分為相同橢球體的投影坐標(biāo)轉(zhuǎn)換和不同橢球體間的投影坐標(biāo)轉(zhuǎn)換。前者的轉(zhuǎn)換框架可見圖4。后者采用七參數(shù)法[9]進(jìn)行轉(zhuǎn)換，其中未假設(shè)3個(gè)旋轉(zhuǎn)角很小，轉(zhuǎn)換框架可見圖5。轉(zhuǎn)換前后高程坐標(biāo)不變。

圖4 相同橢球體的投影坐標(biāo)轉(zhuǎn)換框架

圖5 不同橢球體間的投影坐標(biāo)轉(zhuǎn)換框架

2.3 轉(zhuǎn)換模型保存

(2)

在坐標(biāo)轉(zhuǎn)換和偏移完成后，需要將轉(zhuǎn)換后模型保存到OSGB文件。理論上，可以調(diào)用osgDB::Registry類實(shí)例的writeNode函數(shù)將模型保存為OSGB文件。但這樣做存在格式兼容性風(fēng)險(xiǎn)。如果ContextCapture不能支持所有版本的OSGB格式，則該方法生成的OSGB文件可能不能被ContextCapture準(zhǔn)確讀取。

一種簡(jiǎn)單的方法是直接對(duì)原始OSGB文件進(jìn)行坐標(biāo)數(shù)據(jù)替換。在OSGB文件的數(shù)據(jù)流中，一個(gè)Geometry類對(duì)象的頂點(diǎn)坐標(biāo)連續(xù)放置，可被視為一個(gè)浮點(diǎn)型數(shù)組，而PagedLOD類對(duì)象的用戶自定義中心坐標(biāo)可被視為一個(gè)具有3個(gè)元素的浮點(diǎn)型數(shù)組。設(shè)轉(zhuǎn)換前的數(shù)組為Asrc，存儲(chǔ)了N個(gè)三維坐標(biāo)，元素?cái)?shù)目為3N，對(duì)應(yīng)的轉(zhuǎn)換和偏移后的數(shù)組為Adst，則以字節(jié)為單位的Asrc和Adst的大小nLength均為12N(單精度浮點(diǎn)型)或24N(雙精度浮點(diǎn)型)。通過在文件數(shù)據(jù)流中比較長度為nLength的連續(xù)內(nèi)存的值是否和Asrc相同來確定數(shù)據(jù)替換位置，在找到替換位置后，用Adst的值覆蓋原始數(shù)據(jù)。

3 實(shí)現(xiàn)方法

利用C++和開源軟件開發(fā)OSGB模型的投影坐標(biāo)轉(zhuǎn)換軟件，其中Qt[13]用作編程框架，OpenSceneGraph 3.7[4]用于三維模型讀取，PROJ 6.0[14]用于地圖投影，GDAL 3.0[15]用于坐標(biāo)系設(shè)置，Eigen3[16]用于矩陣運(yùn)算。

3.1 坐標(biāo)轉(zhuǎn)換流程

圖 6為設(shè)計(jì)的坐標(biāo)轉(zhuǎn)換流程。

圖6 模型投影坐標(biāo)轉(zhuǎn)換流程

3.2 轉(zhuǎn)換前模型解析

3.2.1 源坐標(biāo)系解析

如果模型源坐標(biāo)系以EPSG編碼的形式提供，可調(diào)用OGRSpatialReference類的importFromEPSG函數(shù)導(dǎo)入編碼對(duì)應(yīng)的坐標(biāo)系；如果模型源坐標(biāo)系以WKT文本的形式提供，可調(diào)用OGRSpatialReference類的importFromWkt函數(shù)導(dǎo)入WKT對(duì)應(yīng)的坐標(biāo)系。

3.2.2 OSGB文件解析

要讀取OSGB文件中的模型，可先調(diào)用osgDB::readNodeFile函數(shù)獲得Node類型的模型對(duì)象。之后，要獲取待轉(zhuǎn)換坐標(biāo)，需要對(duì)OSG定義的類NodeVisitor進(jìn)行派生并重寫相應(yīng)的apply虛函數(shù)。示例代碼如下：

class CTVisitor: public NodeVisitor {

public:

virtual void apply(Geode &node);

virtual void apply(Group &node);

}

對(duì)Geometry類對(duì)象的解析可放在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為頂點(diǎn)坐標(biāo)數(shù)組指針

Vec3Array *pVt =dynamic_cast

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

}

}

traverse(node);

如果模型采用雙精度浮點(diǎn)數(shù)存儲(chǔ)頂點(diǎn)坐標(biāo)，則需要將上述代碼中的Vec3Array替換為Vec3dArray。對(duì)PagedLOD類對(duì)象的解析可放在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為用戶自定義中心坐標(biāo)

Vec3 cent = lod.getCenter();

}

}

traverse(node);

3.3 投影坐標(biāo)轉(zhuǎn)換

3.3.1 大地坐標(biāo)與投影坐標(biāo)間的轉(zhuǎn)換

大地坐標(biāo)與投影坐標(biāo)間的轉(zhuǎn)換采用pj_transform函數(shù)，該函數(shù)定義如下：

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

其中，src為源坐標(biāo)系對(duì)象，dst為目標(biāo)坐標(biāo)系對(duì)象，pt_count為待轉(zhuǎn)換點(diǎn)數(shù)，pt_offset為首個(gè)待轉(zhuǎn)換點(diǎn)在數(shù)組中的序號(hào)，x，y，z分別為存儲(chǔ)三維坐標(biāo)的數(shù)組的地址。在PROJ中，大地和投影坐標(biāo)系對(duì)象需用PROJ字符串定義。大地坐標(biāo)系的PROJ字符串為：

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

高斯-克呂格投影坐標(biāo)系的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為大地坐標(biāo)系使用的橢球體的字符串，1954年北京坐標(biāo)系是“krass”，1980年西安坐標(biāo)系是“IAU76”，2000國家大地坐標(biāo)系是“GRS80”；dCentLon為中央子午線以角度制表示的經(jīng)度，dFalseEast為以m為單位的東偏，dFalseNorth為以m為單位的北偏。之后，調(diào)用pj_init_plus函數(shù)，輸入?yún)?shù)為坐標(biāo)系的PROJ字符串，即可獲得projPJ類型的坐標(biāo)系對(duì)象。

3.3.2 大地坐標(biāo)與空間直角坐標(biāo)間的轉(zhuǎn)換

從大地坐標(biāo)轉(zhuǎn)換到空間直角坐標(biāo)采用pj_geodetic_to_geocentric函數(shù)，該函數(shù)定義如下：

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

從空間直角坐標(biāo)轉(zhuǎn)換到大地坐標(biāo)采用pj_geocentric_to_geodetic函數(shù)，該函數(shù)定義如下：

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為待轉(zhuǎn)換點(diǎn)數(shù)，pt_offset為首個(gè)待轉(zhuǎn)換點(diǎn)在數(shù)組中的序號(hào)，x，y，z分別為存儲(chǔ)三維坐標(biāo)的數(shù)組的地址。

3.4 轉(zhuǎn)換后模型生成

在生成轉(zhuǎn)換后模型時(shí)，需要確保轉(zhuǎn)換前后目錄結(jié)構(gòu)不變，OSGB文件名和相對(duì)路徑不變。

3.4.1 metadata.xml生成

在轉(zhuǎn)換后模型的metadata.xml中，對(duì)于EPSG已定義的投影坐標(biāo)系，可將SRS標(biāo)簽值設(shè)置為坐標(biāo)系的EPSG編碼；對(duì)于EPSG未定義的投影坐標(biāo)系，可先調(diào)用OGRSpatialReference類的SetWellKnownGeogCS函數(shù)設(shè)置大地坐標(biāo)系，再調(diào)用SetTM函數(shù)設(shè)置中央子午線經(jīng)度、中央緯線緯度、中央子午線縮放因子、東偏、北偏等參數(shù)，調(diào)用SetProjCS函數(shù)設(shè)置坐標(biāo)系名稱，調(diào)用exportToWkt函數(shù)獲得坐標(biāo)系的WKT，最后將SRS標(biāo)簽值設(shè)置為WKT。SRSOrigin標(biāo)簽值為Δxdst，Δydst,Δzdst。

3.4.2 OSGB文件生成

設(shè)Asrc的數(shù)組指針為pSrc，Adst的數(shù)組指針為pDst，以字節(jié)為單位的OSGB文件大小為nOsgbSize。先將轉(zhuǎn)換前的OSGB文件全部讀入unsigned char類型、大小為nOsgbSize的數(shù)組中。設(shè)數(shù)組指針為pszOsgb，則坐標(biāo)數(shù)據(jù)替換的示例代碼如下：

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內(nèi)的數(shù)據(jù)保存到轉(zhuǎn)換后目錄下的OSGB文件。

3.4.3 S3C文件坐標(biāo)系設(shè)置

雖然S3C文件也存儲(chǔ)坐標(biāo)系信息，但由于新版S3C文件的格式已做加密，難以用編程的方法修改其坐標(biāo)系信息。因此，先將轉(zhuǎn)換前的S3C文件復(fù)制到轉(zhuǎn)換后目錄下，再利用CC_S3CComposer.exe將S3C文件中的坐標(biāo)系信息替換為metadata.xml中的坐標(biāo)系信息。

4 軟件實(shí)例

圖7是軟件界面，支持中國常用大地坐標(biāo)系基礎(chǔ)上的高斯-克呂格投影；支持自定義和自動(dòng)計(jì)算輸出坐標(biāo)偏移量；支持手動(dòng)和文件輸入橢球體轉(zhuǎn)換七參數(shù)。

圖7 軟件界面

在Intel Core i5-2300的CPU、24GB內(nèi)存的臺(tái)式機(jī)上，軟件將具有6個(gè)瓦片、1179個(gè)OSGB文件的三維模型從2000國家大地坐標(biāo)系基礎(chǔ)上120°中央子午線的投影坐標(biāo)系轉(zhuǎn)換到1954年北京大地坐標(biāo)系基礎(chǔ)上120.75°中央子午線的投影坐標(biāo)系，耗時(shí)20.214 s，平均單個(gè)OSGB文件耗時(shí)17 ms。在利用CC_S3CComposer.exe更新S3C文件的坐標(biāo)系信息后，可以用ContextCapture Viewer[1]打開轉(zhuǎn)換后的模型，幾何無變形，紋理無缺失，顯示無閃爍。為了對(duì)比，由人工在目標(biāo)坐標(biāo)系下重新生產(chǎn)三維模型。通過在軟件轉(zhuǎn)換后模型和人工重新生產(chǎn)的模型上選取同名特征點(diǎn)并查詢其坐標(biāo)，發(fā)現(xiàn)坐標(biāo)一致。

5 結(jié)束語

文中提出的投影坐標(biāo)轉(zhuǎn)換方法簡(jiǎn)單高效，對(duì)OSGB文件的格式要求不高，不受是否采用LOD技術(shù)、有無紋理貼圖、紋理壓縮比例、切塊方式等因素的影響。由于OSG基本支持所有主流的三維模型格式，所以文中方法容易擴(kuò)展到其它二進(jìn)制的三維模型格式。此外，坐標(biāo)數(shù)據(jù)替換法也適用于其它需要對(duì)模型坐標(biāo)進(jìn)行修改的情況。

文中方法未修改Geometry類對(duì)象的法向量，但理論上，修改頂點(diǎn)坐標(biāo)之后需要重新計(jì)算法向量。雖然在第4節(jié)例子中不更新法向量在視覺上沒有明顯影響，但更新模型法向量仍是一個(gè)改進(jìn)方向。