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);

	{ // mode1
		GDALDataset* 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];
		else
			poSrcData = 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;
	}

	{ // mode2
		GDALDataset* 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];
				else
					poSrcData = 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);
	}

	{ // mode3
		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++) {
				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];
				else
					poSrcData = 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