影像瓦片切割

#include "gdal_priv.h"  
#include "ogrsf_frmts.h"  
#include <algorithm>

bool ImageTile::processByscale(std::string srcimagename, std::string outpath, int tilesize, std::string formattype, double scale, int threadSize, int levelNo)
{
	//初始化成员变量
	GDALAllRegister();
	m_format = formattype;
	m_srcimagename = srcimagename;
	m_outputstr = outpath;
	m_levelNo = levelNo;
	getOGRGeoTransform(); // 获取图层参数信息
	this->tilesize = tilesize;
	double level = scale2Level(scale);
	int xNum = 1.0*m_srcWidth / tilesize / level + 1;
	int yNum = 1.0*m_srcHeight / tilesize / level + 1;

	std::string levelNoPath;
	int x = 0;
	int count = GetSpliteCols();
	while (x < xNum) //列数
	{
		int xCount = count;
		if (x + xCount > xNum) {
			xCount = xNum - x;
		}
		
		if (m_levelNo <= 0) {
			if (!spliteNCols(xNum, yNum, x, xCount, level, m_outputstr, m_format, m_srcimagename, m_bandcount, threadSize)) {
				std::cout << "Failed spliteNCols" << std::endl;
				return false;
			}
		}
		else {			
			if (0 == x) {
				std::stringstream ss;
				ss << "L" << m_levelNo;
				levelNoPath = TDAPICommon::JoinPath(m_outputstr, ss.str());
				if (!CreatePath(levelNoPath, yNum, xNum)) {
					std::cout << "Failed CreatePath," << levelNoPath << std::endl;
					return false;
				}
			}
			if (!spliteElectronicMap(xNum, yNum, x, xCount, level, levelNoPath, m_format, m_srcimagename, m_bandcount, threadSize)) {
				std::cout << "Failed spliteElectronicMap" << std::endl;
				return false;
			}
		}

		x += xCount;
	}

	return true;
}

bool ImageTile::spliteElectronicMap(int xNum, int yNum, int x, int xCount, double level, std::string filepath, std::string formattype, std::string srcimagename, int bandCount, int threadSize)
{
	std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();

	std::string dstTileSuffixName = TDRasterCommon::GetSuffixByDriver(formattype);
	if (dstTileSuffixName.empty()) {
		std::cout << "No file suffix corresponding to " << formattype << " driver was found" << std::endl;
		return false;
	}

	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
	auto armDriver = GetGDALDriverManager()->GetDriverByName(formattype.c_str());
	auto poDriver = GetGDALDriverManager()->GetDriverByName("MEM");

	if (nullptr == armDriver) {
		std::cout << "Gdal driver " << formattype << "does not exist" << std::endl;
		return false;
	}

	auto srcDataset = (GDALDataset*)GDALOpen(srcimagename.c_str(), GA_ReadOnly); //获取原始数据句柄
	if (nullptr == srcDataset) {
		std::cout << " Failed to open " << srcimagename << std::endl;
		return false;
	}
	auto pProjection = srcDataset->GetProjectionRef();
	if (nullptr == pProjection) {
		std::cout << "GetProjectionRef returns nullptr" << std::endl;
		return false;
	}

	//GDALRasterIOExtraArg exterArg;
	//INIT_RASTERIO_EXTRA_ARG(exterArg);
	//exterArg.eResampleAlg = GDALRIOResampleAlg::GRIORA_Bilinear;

	//读取影响的位置
	double startX = x * tilesize *level;
	double startY = 0;
	double sw = (m_srcWidth - (xNum - 1)*tilesize * level) / (tilesize * level);
	double w = m_srcWidth - (xNum - 1)*tilesize * level;
	double sh = (m_srcHeight - (yNum - 1)*tilesize * level) / (tilesize * level);
	double h = m_srcHeight - (yNum - 1)*tilesize * level;

	//
	size_t storeWidth = 0;
	//读取影像的宽度
	size_t readWidth = 0;

	if (x + xCount == xNum) {
		readWidth = tilesize * (xCount - 1)*level + w;
		storeWidth = tilesize * (xCount - 1) + tilesize * sw;
	}
	else {
		readWidth = tilesize * level* xCount;
		storeWidth = tilesize * xCount;
	}

	//读取影像的高度
	size_t readHeight = m_srcHeight;
	size_t storeHeight = tilesize * (yNum - 1) + tilesize * sh;

	//开辟空间
	size_t imgBufNum = storeWidth * storeHeight * bandCount;
	GByte *imgBuf = new GByte[imgBufNum];
	memset(imgBuf, 0, imgBufNum);

	//读取N列瓦片影像
	srcDataset->RasterIO(GF_Read, startX, startY, readWidth, readHeight, imgBuf, storeWidth, storeHeight,
		GDT_Byte, bandCount, nullptr, 0, 0, 0);
	std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
	std::cout << "read raster time : " << std::chrono::duration_cast<std::chrono::seconds>(end - start).count() << " s" << std::endl;
	ThreadPool::ThreadPool threadPool(1, threadSize);
	for (size_t y = 0; y < yNum; ++y) {
		size_t tileSizeX = tilesize;
		size_t tileSizeY = tilesize;
		if ((yNum - 1) == y) {
			tileSizeY = tilesize * sh;
		}

		for (size_t i = 0; i < xCount; ++i) {
			size_t curX = x + i;
			if ((xNum - 1) == curX) {
				tileSizeX = tilesize * sw;
			}

			threadPool.CommitOneTask([this, curX, i, y, level, tileSizeX, tileSizeY, storeWidth, storeHeight, armDriver, poDriver, filepath, dstTileSuffixName, bandCount, srcDataset, imgBuf, pProjection]()
			{
				std::stringstream ss;
				ss << filepath << "/R" << curX << "/C" << y;
				std::string outputFilePath = ss.str();
				TDAPICommon::CreateDir(outputFilePath);

				ss.str("");
				ss << outputFilePath << "/" << curX << y << m_levelNo;
				std::string outputBaseFileName = ss.str();
				std::string  memName = outputBaseFileName + ".tif"; // 配置输出参数
				std::string  target_name = outputBaseFileName + dstTileSuffixName; // 配置输出参数				
				GDALDataset* poDstDS;
				poDstDS = poDriver->Create(memName.c_str(), tilesize, tilesize, bandCount, GDT_Byte, nullptr);
				if (nullptr == poDstDS) {
					return;
				}
				//设置放射变换参数		
				if (nullptr != pProjection) {
					poDstDS->SetProjection(pProjection);
				}
				
				double startX = curX * tilesize*level;
				double startY = y * tilesize*level;
				double xx, yy;
				TDRasterCommon::ImageRowCol2Projection(m_adfGeoTransform, startX, startY, xx, yy);
				double updateGeoTransfrom[6] = { 0.0 };
				poDstDS->SetGeoTransform(UpdateGeoTransfrom(m_adfGeoTransform, xx, yy, level, updateGeoTransfrom)); //设置图像参数

				GByte *pStart = new GByte[tileSizeX*tileSizeY*bandCount];

				memset(pStart, 0, tileSizeX*tileSizeY*bandCount);

				//计算原始影像一行瓦片的长度
				size_t rowTileWidth = storeWidth * tilesize;
				//计算N行瓦片的长度
				int nRowTileWidth = rowTileWidth * y;

				for (size_t iBand = 0;iBand < bandCount;++iBand) {
					//计算数据瓦片的波段起始位置 
					size_t dstStart = tileSizeX * tileSizeY * iBand;
					//计算读取的原始影像的影像的起始位置
					size_t srcStart = storeWidth * storeHeight * iBand;

					for (size_t j = 0; j < tileSizeY; ++j) {
						memcpy(pStart + dstStart + j * tileSizeX, imgBuf + srcStart + nRowTileWidth + j * storeWidth + i * tilesize, tileSizeX);
					}
				}

				poDstDS->RasterIO(GF_Write, 0, 0, tileSizeX, tileSizeY, pStart, tileSizeX, tileSizeY,
					GDT_Byte, bandCount, nullptr, 0, 0, 0);

				create_copy(poDstDS, target_name, armDriver);
				GDALClose(poDstDS);
				poDstDS = NULL;

				delete[] pStart;
			});
		}
	}

	threadPool.StopThreadPool();
	delete[] imgBuf;
	imgBuf = nullptr;
	GDALClose(srcDataset);
	srcDataset = NULL;
	end = std::chrono::steady_clock::now();
	std::cout << "handle one time is : " << std::chrono::duration_cast<std::chrono::seconds>(end - start).count() << " s" << std::endl;
	return true;
}

bool ImageTile::spliteNCols(int xNum, int yNum, int x, int xCount, double level, std::string filepath, std::string formattype, std::string srcimagename, int bandCount, int threadSize)
{
	std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();

	std::string dstTileSuffixName = TDRasterCommon::GetSuffixByDriver(formattype);
	if (dstTileSuffixName.empty()) {
		std::cout << "No file suffix corresponding to "<< formattype << " driver was found" << std::endl;
		return false;
	}

	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
	auto armDriver = GetGDALDriverManager()->GetDriverByName(formattype.c_str());

	if (nullptr == armDriver) {
		std::cout << "Gdal driver " << formattype << "does not exist" << std::endl;
		return false;
	}	
	
	auto srcDataset = (GDALDataset*)GDALOpen(srcimagename.c_str(), GA_ReadOnly); //获取原始数据句柄
	if (nullptr == srcDataset){
		std::cout << " Failed to open " << srcimagename << std::endl;
		return false;
	}
	auto pProjection = srcDataset->GetProjectionRef();
	if (nullptr == pProjection) {
		std::cout << "GetProjectionRef returns nullptr" << std::endl;
		return false;
	}

	//GDALRasterIOExtraArg exterArg;
	//INIT_RASTERIO_EXTRA_ARG(exterArg);
	//exterArg.eResampleAlg = GDALRIOResampleAlg::GRIORA_Bilinear;

	//读取影响的位置
	double startX = x * tilesize *level;
	double startY = 0;
	double sw = (m_srcWidth - (xNum - 1)*tilesize * level) / (tilesize * level);
	double w = m_srcWidth - (xNum - 1)*tilesize * level;
	double sh = (m_srcHeight - (yNum - 1)*tilesize * level) / (tilesize * level);
	double h = m_srcHeight - (yNum - 1)*tilesize * level;

	//
	size_t storeWidth = 0;
	//读取影像的宽度
	size_t readWidth = 0;

	if (x + xCount == xNum) {
		readWidth = tilesize*(xCount - 1)*level + w;
		storeWidth = tilesize * (xCount - 1) + tilesize * sw;
	}
	else {
		readWidth = tilesize*level* xCount;
		storeWidth = tilesize * xCount;
	}

	//读取影像的高度
	size_t readHeight = m_srcHeight;
	size_t storeHeight = tilesize * (yNum - 1) + tilesize * sh;

	//开辟空间
	size_t imgBufNum = storeWidth * storeHeight * bandCount;
	GByte *imgBuf = new GByte[imgBufNum];
	memset(imgBuf,0, imgBufNum);

	//读取N列瓦片影像
	srcDataset->RasterIO(GF_Read, startX, startY, readWidth, readHeight, imgBuf, storeWidth, storeHeight,
		GDT_Byte, bandCount, nullptr, 0, 0, 0);
	std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
	std::cout << "read raster time : " << std::chrono::duration_cast<std::chrono::seconds>(end - start).count() << " s" << std::endl;

	
	TDThreadPool::ThreadPool threadPool(1, threadSize);
	for (size_t y = 0; y < yNum; ++y) {
		size_t tileSizeX = tilesize;
		size_t tileSizeY = tilesize;
		if ((yNum - 1) == y) {
			tileSizeY = tilesize * sh;
		}

		for (size_t i = 0; i < xCount; ++i) {
			size_t curX = x + i;
			if ((xNum - 1) == curX) {
				tileSizeX = tilesize * sw;
			}
			std::string  target_name = filepath + "/R" + DecIntToHexStr(y) + "C" + DecIntToHexStr(curX) + dstTileSuffixName; // 配置输出参数

			threadPool.CommitOneTask([this, curX,i,y, level,tileSizeX, tileSizeY, storeWidth, storeHeight, armDriver,target_name, bandCount, srcDataset, imgBuf, pProjection]()
			{
				GDALDataset* poDstDS;
				poDstDS = armDriver->Create(target_name.c_str(), tilesize, tilesize, bandCount, GDT_Byte, nullptr);
				if (nullptr == poDstDS) {
					return;
				}
				//设置放射变换参数			
				if (nullptr != pProjection) {
					poDstDS->SetProjection(pProjection);
				}

				double startX = curX * tilesize*level;
				double startY = y * tilesize*level;
				double xx, yy;
				TDRasterCommon::ImageRowCol2Projection(m_adfGeoTransform, startX, startY, xx, yy);
				double updateGeoTransfrom[6] = {0.0};
				poDstDS->SetGeoTransform(UpdateGeoTransfrom(m_adfGeoTransform, xx, yy, level, updateGeoTransfrom)); //设置图像参数

				GByte *pStart = new GByte[tileSizeX*tileSizeY*bandCount];
				
				memset(pStart, 0, tileSizeX*tileSizeY*bandCount);

				//计算原始影像一行瓦片的长度
				size_t rowTileWidth = storeWidth * tilesize;
				//计算N行瓦片的长度
				size_t nRowTileWidth = rowTileWidth * y;

				for (size_t iBand = 0;iBand < bandCount;++iBand) {
					//计算数据瓦片的波段起始位置 
					size_t dstStart = tileSizeX * tileSizeY * iBand;
					//计算读取的原始影像的影像的起始位置
					size_t srcStart = storeWidth * storeHeight * iBand;

					for (size_t j = 0; j < tileSizeY; ++j) {
						memcpy(pStart + dstStart + j * tileSizeX, imgBuf + srcStart + nRowTileWidth + j* storeWidth + i * tilesize, tileSizeX);
					}
				}

				poDstDS->RasterIO(GF_Write, 0, 0, tileSizeX, tileSizeY, pStart, tileSizeX, tileSizeY,
					GDT_Byte, bandCount, nullptr, 0, 0, 0);

				GDALClose(poDstDS);
				poDstDS = NULL;

				delete[] pStart;
			});
		}
	}

	threadPool.StopThreadPool();
	delete[] imgBuf;
	imgBuf = nullptr;
	GDALClose(srcDataset);
	srcDataset = NULL;
	end = std::chrono::steady_clock::now();
	std::cout << "handle one time is : " << std::chrono::duration_cast<std::chrono::seconds>(end - start).count() << " s" << std::endl;
	return true;
}