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

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

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

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<double>& mins, 
    std::vector<double>& 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=" << S_MIN
            << " min=" << min << " max=" << max << std::endl;
        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, int szInPixel, 
    QString outDir,
    QList<QPoint> validBlocks,           /*only valid blocks will be saved*/
    std::function<void(int&)> cb)
{
    //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);
    unsigned int 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 " 
        << outDir.toStdString() 
        <<" data type :" << dataType
        <<" dataType size:" << dataTypeSizeInBytes
        << std::endl;
    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;
        GDALDatasetH out = 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