python遥感图像开发小软件_遥感影像深度学习标注软件的开发要点

最近一个月开发了一款遥感影像深度学习标注软件。经过一个月的艰苦编码，基本已经稳定， 讲开发过程做一个简单记录，以备后用。

一、遥感影像的标注与图片影像有何不同

1、遥感影像文件尺寸大，单幅影像动辄达到几百M甚至上T的数据量。 用labelme之类的软件无法打开。

2、遥感影像很多是16Bit  32BiT 的数据，不经过数据类型转换和拉伸处理无法正常显示。

3、数据成果的分幅问题， 普通图像可以直接导出标注结果， 但遥感图像的标注结果想要被深度学习算法利用，必须要经过分幅处理。

4、多波段数据需要选择波段。

针对这些遥感影像本身的特点， 我们基于labelme， 做了针对性的改进。

首先，大幅的影像，我们需要给影像做金字塔文件， 在GDAL提供了很方便的金字塔函数:

CPLErr CPL_DLL CPL_STDCALL

GDALBuildOverviews( GDALDatasetH,const char *, int, int *,int, int *, GDALProgressFunc, void * );

经过这样的处理，我们就拥有了多层次的图像数据，单这仅仅是第一步，为了平滑显示遥感影像， 需要对各层次影像进行切块，也就是LOD技术(Levels of Detail的简称，意为多细节层次)。

实现LOD技术，不但在让标注软件在浏览影像时候更平滑顺畅，而且方便用户在不同尺度上对物体进行标注。

其次，要解决16BIT 和 32 bit 影像的显示问题。 我们都知道，计算机显示图像时，只能用RGB三原色模型， 每个波段只能用8BIT显示，  那么在显示图像的时候必然会涉及到scalar type的转换，

虽然GDAL 在 rasterIO 的时候，也可自动实现颜色scalar type转换， 但是这往往还不够，为了使得显示的颜色更加饱满，通常需要去掉 直方图中前后各0.2% 的像素值。

/**/

void getHistogramMinMax(GDALDatasetH dataset, std::vector&mins,

std::vector& maxs, GDALProgressFunc callback = 0L) {

GDALDataType datatype= GDALGetRasterDataType(GDALGetRasterBand(dataset, 1));

ScalarType scalarType;if (datatype == GDT_UInt16) scalarType =UINT16;if (datatype == GDT_Int16) scalarType =SINT16;if (datatype == GDT_UInt32) scalarType =UINT32;if (datatype == GDT_Int32) scalarType =SINT32;int bands =GDALGetRasterCount(dataset);for (int i = 0; i < bands; i++)

{//scalar min

const float64 S_MIN =defaultMin(scalarType);//scalar max

const float64 S_MAX =defaultMax(scalarType);int bins = S_MAX - S_MIN + 1;int *panHistogram = new int[bins];double *accuHistogram = new double[bins];

GDALGetRasterHistogram(GDALGetRasterBand(dataset,1),

S_MIN- 0.5, S_MAX + 0.5,

bins, panHistogram,false, false, 0, 0);/*acc histo*/

for (int j = 0; j < bins; j++) {

accuHistogram[j]=panHistogram[j];

}for (int j = 1; j < bins; j++) {

accuHistogram[j]+= accuHistogram[j - 1];

}

size_t samples=GDALGetRasterXSize(dataset);

size_t lines=GDALGetRasterYSize(dataset);for (int j = 0; j < bins; j++) {

accuHistogram[j]/= accuHistogram[bins - 1];

}int min = 0, max = 0;for (int u = 0; u < bins; u++) {if (accuHistogram[u] > 0.02) {

min=u;break;

}

}for (int u = 0; u < bins; u++) {if (accuHistogram[u] > 0.98) {

max=u;break;

}

}

std::cout<< "= getHistogramMinMax:" << "S_MIN=" <

mins.push_back(S_MIN+min);

maxs.push_back(S_MIN+max);delete[] panHistogram;delete[] accuHistogram;

}return;

}

对于大幅的遥感影像而言， 如何将标注成果导出为CNN深度模型的输入，也是一个必要要解决的问题，在GDAL中， 有一个很好用的python脚本工具, 但是这个分块工具的用途是给WMS服务用的， 涉及许多的坐标系变换，用起来相当复杂，

我们采用自己分幅的方式， 将成果输出为固定大小的TILE.    这个时候要注意， 分出来的每个小块还是需要保留该区域的地理坐标信息， 即保存它的 GeoTrans[6]  这个属性。  在这个过程中有一个地方要注意，  GeoTrans[5]的值一般是一个负数。

bool Utils::gdal2Tile(QString image, intszInPixel,

QString outDir,

QList validBlocks, /*only valid blocks will be saved*/std::functioncb)

{//open file

char *strImage =image.toLocal8Bit().data();

GDALDatasetH dataset=GDALOpen(strImage, GA_ReadOnly);if (!dataset) {return false;

}double geoTrans[6];

GDALGetGeoTransform(dataset, geoTrans);/*get wkt*/std::string wkt = (char*)GDALGetProjectionRef(dataset);//get image info

unsigned int height =GDALGetRasterYSize(dataset);

unsignedint width =GDALGetRasterXSize(dataset);int bsNum =GDALGetRasterCount(dataset);

GDALDataType dataType= GDALGetRasterDataType(GDALGetRasterBand(dataset, 1));int dataTypeSizeInBytes = GDALGetDataTypeSize(dataType) / 8;char* buf = new char[szInPixel * szInPixel * bsNum *dataTypeSizeInBytes];

std::cout<< "= gdal2Tile: outDir"

<

ossimString driverName;char **option = 0;if (image.endsWith("img", Qt::CaseInsensitive)) {

driverName= "HFA";

option= CSLSetNameValue(option, "BLOCKSIZE", "256");

}else if (image.endsWith("ecw", Qt::CaseInsensitive)){

driverName= "ECW";

}else if (image.endsWith("gta", Qt::CaseInsensitive)){

driverName= "GTA";

}else if (image.endsWith("pix", Qt::CaseInsensitive)){

driverName= "PCIDSK";

}else if (image.endsWith("hdr", Qt::CaseInsensitive)){

driverName= "ENVI";

option= CSLSetNameValue(option, "INTERLEAVE", "bip");

option= CSLSetNameValue(option, "SUFFIX", "REPLACE");

}else{

driverName= "GTIFF";

option= CSLSetNameValue(option, "INTERLEAVE", "PIXEL");

option= CSLSetNameValue(option, "COMPRESS", "NONE");

option= CSLSetNameValue(option, "TILED", "YES");

option= CSLSetNameValue(option, "BLOCKXSIZE", "256");

option= CSLSetNameValue(option, "BLOCKYSIZE", "256");

option= CSLSetNameValue(option, "BIGTIFF", "IF_NEEDED");

}

GDALDriverH driver=GDALGetDriverByName(driverName.c_str());int tile_x_num = std::ceil(1.0 * width /szInPixel);int tile_y_num = std::ceil(1.0 * height /szInPixel);for (int u = 0; u < validBlocks.size(); u++) {int i =validBlocks[u].y();int j =validBlocks[u].x();int xsize = (j != tile_x_num - 1) ? szInPixel : width - j*szInPixel;int ysize = (i != tile_y_num - 1) ? szInPixel : height - i*szInPixel;/*read all bands*/CPLErr err= GDALDatasetRasterIO(dataset, GF_Read, j*szInPixel,

i*szInPixel, xsize, ysize, buf, xsize, ysize, dataType,

bsNum,0,

bsNum* dataTypeSizeInBytes, bsNum * dataTypeSizeInBytes *szInPixel, dataTypeSizeInBytes);/*create new dataset, first get the output filename*/image= image.replace('\\', '/');int extPos = image.lastIndexOf('.');int lastSlashPos = image.lastIndexOf('/');

QString baseName= image.mid(lastSlashPos, extPos -lastSlashPos);

QString ext= image.mid(extPos, -1);if (outDir.endsWith('/') || outDir.endsWith('\\')) {

outDir= outDir.replace('\\', '/');

}else{

outDir+= '/';

outDir= outDir.replace('\\', '/');

}

QString outImage= outDir +baseName+ QString("_")+QString::number(i)+ QString("_") +QString::number(j)+ext;

GDALDatasetHout =GDALCreate(driver,

outImage.toLocal8Bit().data(),

xsize,

ysize,

bsNum,

dataType,

option);/*write to new dataset*/GDALDatasetRasterIO(out, GF_Write, 0,0, xsize, ysize, buf, xsize, ysize, dataType,

bsNum,0,

bsNum* dataTypeSizeInBytes, bsNum * dataTypeSizeInBytes *szInPixel, dataTypeSizeInBytes);if (!wkt.empty()) {

GDALSetProjection(out, wkt.c_str());double newGeoTrans[6];

memcpy(newGeoTrans, geoTrans,sizeof(double) * 6);

newGeoTrans[0] = geoTrans[0] + geoTrans[1] * j *szInPixel;

newGeoTrans[3] = geoTrans[3] + geoTrans[5] * i *szInPixel;

GDALSetGeoTransform(out, newGeoTrans);

}else{

GDALSetProjection(out, wkt.c_str());double newGeoTrans[6];

newGeoTrans[0] = 0;

newGeoTrans[1] = 1;

newGeoTrans[2] = 0;

newGeoTrans[3] =ysize;

newGeoTrans[4] = 0;

newGeoTrans[5] = -1;

}

GDALClose(out);int percent = (u)*100.0 /validBlocks.size();

cb(percent);

}

GDALClose(dataset);

CSLDestroy(option);delete[] buf;int percent = 100;

cb(percent);return true;

}

5、 倾斜地物的标注

许多遥感影像中的地物，如球场， 房屋  等，  是长方形的几何体， 但是由于拍摄的问题，我们得到的是倾斜的长方体， 这种情况，如果用矩形框进行标注，会框选大部分的无效区域， 如果用

多边形标注，又会显得麻烦， 这时，如果用倾斜多边形标注会非常简洁。

这款遥感影像标注软件是基于labelme的改进，目前在开发中，不过已经很好地解决了上述问题.

6.  标记软件与数字地球的结合

数字地球提供了一个全局的环境，能帮助标注者辅助判读图像，CESIUM是一款基于webgl技术的，BS架构的开源数字地球软件， 相比worldwind java 和OSGEARTH ， 这个社区更为

活跃，有丰富的案例和代码可供参考。在RSLABEL标注软件，提供了插件机制， 将CESIUM作为一个模块嵌入到工作环境中，用户可将正在编辑的影像嵌入到数字地球上，在三维地形

的配合下进行判读。 网页嵌入到应用程序中，需要用到QWebView ,  该控件可以让桌面应用程序装载web页面，   通过evaluejavascript 函数让网页执行应用程序的指令。

https://github.com/enigma19971/RSLabel