GDAL中GDALDataset::RasterIO分块读取的实现

GDALDataset类中的RasterIO函数能够对图像任意指定区域、任意波段的数据按指定数据类型、指定排列方式读入内存和写入文件中，因此可以实现对大影像的分块读、写运算操作。针对特大的影像图像，有时为了减少内存消耗，对图像进行分块读取很有必要。在以下的测试代码中，给出了3种方式，每种方式的最终结果都是完全相同的，从内存占用情况来看：第一种大于第二种，第二种大于第三种。第三种消耗内存最小。

测试代码如下：

int test_gdal_GDALDataset()
{const char* image_name = "E:/GitCode/GDAL_Test/test_images/3.jpg";GDALAllRegister();GDALDataset* poDataset = (GDALDataset*)GDALOpen(image_name, GA_ReadOnly);if (poDataset == nullptr) {std::cout << "input image error" << std::endl;return -1;}int width = poDataset->GetRasterXSize();int height = poDataset->GetRasterYSize();int band_count = poDataset->GetRasterCount();size_t length = width * height * band_count;GDALDataType gdal_data_type = poDataset->GetRasterBand(1)->GetRasterDataType();int depth = GDALGetDataTypeSize((GDALDataType)gdal_data_type);//fprintf(stderr, "depth: %d\n", depth);assert(depth == 8 || depth == 16);int size_byte = 1;if (depth == 16) size_byte = 2;void* data1 = nullptr;void* data2 = nullptr;void* data3 = nullptr;if (depth == 8) {data1 = new unsigned char[length];memset(data1, 0, length);data2 = new unsigned char[length];memset(data2, 0, length);data3 = new unsigned char[length];memset(data3, 0, length);} else {data1 = new unsigned short[length];memset(data1, 0, length * 2);data2 = new unsigned short[length];memset(data2, 0, length * 2);data3 = new unsigned short[length];memset(data3, 0, length * 2);}GDALClose((GDALDatasetH)poDataset);{ // mode1GDALDataset* poDataset = (GDALDataset*)GDALOpen(image_name, GA_ReadOnly);int band_count = poDataset->GetRasterCount();fprintf(stderr, "mode=1: band_count = %d\n", band_count);int* pBandMap = new int[band_count];for (int i = 0; i < band_count; i++) {pBandMap[i] = i + 1;}GDALDataType gdal_data_type = poDataset->GetRasterBand(1)->GetRasterDataType();int depth = GDALGetDataTypeSize((GDALDataType)gdal_data_type);void* poSrcData = nullptr;if (depth == 8)poSrcData = new unsigned char[width * height * band_count];elsepoSrcData = new unsigned short[width * height * band_count];poDataset->RasterIO(GF_Read, 0, 0, width, height,poSrcData, width, height, gdal_data_type, band_count, pBandMap, 0, 0, 0);if (depth == 8) {unsigned char* p1 = (unsigned char*)poSrcData;for (int y = 0; y < height; y++) {unsigned char* p2 = (unsigned char*)data1 + width * band_count * y;for (int x = 0; x < width; x++) {for (int band = 0; band < band_count; band++) {p2[x * band_count + band] = p1[band * width * height + y * width + x];;}}}} else {unsigned short* p1 = (unsigned short*)poSrcData;for (int y = 0; y < height; y++) {unsigned short* p2 = (unsigned short*)data1 + width * band_count * y;for (int x = 0; x < width; x++) {for (int band = 0; band < band_count; band++) {p2[x * band_count + band] = p1[band * width * height + y * width + x];;}}}}GDALClose((GDALDatasetH)poDataset);delete[] pBandMap;delete[] poSrcData;}{ // mode2GDALDataset* poDataset = (GDALDataset*)GDALOpen(image_name, GA_ReadOnly);int crop_width = 200;int crop_height = 200;int loops_y = height / crop_height;int loops_x = width / crop_width;for (int y = 0; y < loops_y; y++) {for (int x = 0; x < loops_x; x++) {int band_count = poDataset->GetRasterCount();fprintf(stderr, "mode=2: band_count = %d\n", band_count);int* pBandMap = new int[band_count];for (int i = 0; i < band_count; i++) {pBandMap[i] = i + 1;}GDALDataType gdal_data_type = poDataset->GetRasterBand(1)->GetRasterDataType();int depth = GDALGetDataTypeSize((GDALDataType)gdal_data_type);void* poSrcData = nullptr;if (depth == 8)poSrcData = new unsigned char[crop_width * crop_height * band_count];elsepoSrcData = new unsigned short[crop_width * crop_height * band_count];int xOff = crop_width * x;int yOff = crop_height * y;poDataset->RasterIO(GF_Read, xOff, yOff, crop_width, crop_height,poSrcData, crop_width, crop_height, gdal_data_type, band_count, pBandMap, 0, 0, 0);if (depth == 8) {unsigned char* p1 = (unsigned char*)poSrcData;unsigned char* p2 = (unsigned char*)data2 + width * band_count * y * crop_height;for (int m = 0; m < crop_height; m++) {unsigned char* p3 = p2 + width * band_count * m + x * crop_width * band_count;for (int n = 0; n < crop_width; n++) {for (int band = 0; band < band_count; band++) {p3[n * band_count + band] = p1[band * crop_width * crop_height + m * crop_width + n];}}}}else {unsigned short* p1 = (unsigned short*)poSrcData;unsigned short* p2 = (unsigned short*)data2 + width * band_count * y * crop_height;for (int m = 0; m < crop_height; m++) {unsigned short* p3 = p2 + width * band_count * m + x * crop_width * band_count;for (int n = 0; n < crop_width; n++) {for (int band = 0; band < band_count; band++) {p3[n * band_count + band] = p1[band * crop_width * crop_height + m * crop_width + n];}}}}delete[] pBandMap;delete[] poSrcData;}}GDALClose((GDALDatasetH)poDataset);}{ // mode3int crop_width = 200;int crop_height = 200;int loops_y = height / crop_height;int loops_x = width / crop_width;for (int y = 0; y < loops_y; y++) {for (int x = 0; x < loops_x; x++) {GDALDataset* poDataset = (GDALDataset*)GDALOpen(image_name, GA_ReadOnly);int band_count = poDataset->GetRasterCount();fprintf(stderr, "mode=3: band_count = %d\n", band_count);int* pBandMap = new int[band_count];for (int i = 0; i < band_count; i++) {pBandMap[i] = i + 1;}GDALDataType gdal_data_type = poDataset->GetRasterBand(1)->GetRasterDataType();int depth = GDALGetDataTypeSize((GDALDataType)gdal_data_type);void* poSrcData = nullptr;if (depth == 8)poSrcData = new unsigned char[crop_width * crop_height * band_count];elsepoSrcData = new unsigned short[crop_width * crop_height * band_count];int xOff = crop_width * x;int yOff = crop_height * y;poDataset->RasterIO(GF_Read, xOff, yOff, crop_width, crop_height,poSrcData, crop_width, crop_height, gdal_data_type, band_count, pBandMap, 0, 0, 0);if (depth == 8) {unsigned char* p1 = (unsigned char*)poSrcData;unsigned char* p2 = (unsigned char*)data3 + width * band_count * y * crop_height;for (int m = 0; m < crop_height; m++) {unsigned char* p3 = p2 + width * band_count * m + x * crop_width * band_count;for (int n = 0; n < crop_width; n++) {for (int band = 0; band < band_count; band++) {p3[n * band_count + band] = p1[band * crop_width * crop_height + m * crop_width + n];}}}} else {unsigned short* p1 = (unsigned short*)poSrcData;unsigned short* p2 = (unsigned short*)data3 + width * band_count * y * crop_height;for (int m = 0; m < crop_height; m++) {unsigned short* p3 = p2 + width * band_count * m + x * crop_width * band_count;for (int n = 0; n < crop_width; n++) {for (int band = 0; band < band_count; band++) {p3[n * band_count + band] = p1[band * crop_width * crop_height + m * crop_width + n];}}}}GDALClose((GDALDatasetH)poDataset);delete[] pBandMap;delete[] poSrcData;}}}for (int i = 0; i < length * size_byte; i++) {unsigned char* p1 = (unsigned char*)data1;unsigned char* p2 = (unsigned char*)data2;unsigned char* p3 = (unsigned char*)data3;if (p1[i] != p2[i] || p1[i] != p3[i]) {fprintf(stderr, "error: data1 != data2 or data1 != data3\n");return -1;}}delete[] data1;delete[] data2;delete[] data3;return 0;
}

GitHub： https://github.com/fengbingchun/GDAL_Test