使用gdal直接读取某坐标上的位置像素值,

于 2025-03-20 14:22:17 发布
阅读量221 收藏
点赞数 3
文章标签： 人工智能 计算机视觉
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/chijingjing/article/details/146395879
版权

需求：gdal直接跟据某个范围读取像素值会有缩放和重采样，如果二次重采样可能会产生锯齿或者模糊，为了绕开gdal自动采样和金字塔，可读取某个点位得像素值，然后自己重采样

初始化：

		GDALDataset* m_Image=img;
		GDALRasterIOExtraArg m_Arg;
void init(GDALDataset*img,GDALRIOResampleAlg alg)
	{
		 m_Image=img;
	
		INIT_RASTERIO_EXTRA_ARG(m_Arg);
		m_Arg.eResampleAlg = alg;
	}

1：直接临近值读取：

	/// \brief 根据经纬度读取对应位置的像素值
	/// \param longitude 经度
	/// \param latitude 纬度
	/// \param bandCount 波段数量
	/// \param Values 返回的像素值数组
	virtual bool IO(double longitude, double latitude, int bandCount, std::vector<double>& Values)
	{
		if (!m_Image || bandCount < 1 || bandCount > m_Image->GetRasterCount())
			return false;

		// 获取地理变换参数
		double geoTransform[6];
		if (m_Image->GetGeoTransform(geoTransform) != CE_None)
			return false;

		// 计算经纬度对应的行列号
		double col = (longitude - geoTransform[0]) / geoTransform[1];
		double row = (latitude - geoTransform[3]) / geoTransform[5];

		// 检查行列号是否在有效范围内
		if (col < 0 || row < 0 || col >= m_Image->GetRasterXSize() || row >= m_Image->GetRasterYSize())
			return false;

		// 读取每个波段的像素值
		Values.clear();
		for (int band = 1; band <= bandCount; band++)
		{
			GDALRasterBand* pBand = m_Image->GetRasterBand(band);
			if (!pBand)
				return false;
			auto datatype = pBand->GetRasterDataType();
			double value;
			CPLErr err = pBand->RasterIO(GF_Read, static_cast<int>(std::round(col)), static_cast<int>(std::round(row)),
				1, 1, &value, 1, 1, GDT_Float64, 0, 0, &m_Arg);
			if (err != CE_None)
				return false;

			Values.push_back(value);
		}

		return true;
	}

2：boxfilter滤波读取

	/// \brief 根据经纬度读取对应位置的像素值（BoxFilter）
	/// \param longitude 经度
	/// \param latitude 纬度
	/// \param bandCount 波段数量
	/// \param Values 返回的像素值数组
	/// \param kernelSize 滤波核大小（默认 3x3）
	virtual bool IO_BoxFilter(double longitude, double latitude, int bandCount, std::vector<double>& Values, int kernelSize = 3)
	{
		if (!m_Image || bandCount < 1 || bandCount > m_Image->GetRasterCount())
			return false;

		// 获取地理变换参数
		double geoTransform[6];
		if (m_Image->GetGeoTransform(geoTransform) != CE_None)
			return false;

		// 计算经纬度对应的行列号
		double col = (longitude - geoTransform[0]) / geoTransform[1];
		double row = (latitude - geoTransform[3]) / geoTransform[5];

		// 检查行列号是否在有效范围内
		if (col < 0 || row < 0 || col >= m_Image->GetRasterXSize() || row >= m_Image->GetRasterYSize())
			return false;

		// 计算滤波核的半径
		int radius = kernelSize / 2;

		// 读取每个波段的像素值
		Values.clear();
		for (int band = 1; band <= bandCount; band++)
		{
			GDALRasterBand* pBand = m_Image->GetRasterBand(band);
			if (!pBand)
				return false;

			double sum = 0.0;
			int count = 0;

			// 遍历滤波核范围内的像素
			for (int i = -radius; i <= radius; i++)
			{
				for (int j = -radius; j <= radius; j++)
				{
					int x = static_cast<int>(col) + j;
					int y = static_cast<int>(row) + i;

					if (x >= 0 && x < m_Image->GetRasterXSize() && y >= 0 && y < m_Image->GetRasterYSize())
					{
						double value;
						CPLErr err = pBand->RasterIO(GF_Read, x, y, 1, 1, &value, 1, 1, GDT_Float64, 0, 0);
						if (err == CE_None)
						{
							sum += value;
							count++;
						}
					}
				}
			}

			if (count > 0)
				Values.push_back(sum / count);
			else
				Values.push_back(0.0); // 如果没有有效像素，填充默认值
		}

		return true;
	}

3：双线性插值读取：

	/// \brief 根据经纬度读取对应位置的像素值（双线性插值）
	/// \param longitude 经度
	/// \param latitude 纬度
	/// \param bandCount 波段数量
	/// \param Values 返回的像素值数组
	virtual bool IO_BilinearInterpolation(double longitude, double latitude, int bandCount, std::vector<double>& Values)
	{
		if (!m_Image || bandCount < 1 || bandCount > m_Image->GetRasterCount())
			return false;

		// 获取地理变换参数
		double geoTransform[6];
		if (m_Image->GetGeoTransform(geoTransform) != CE_None)
			return false;

		// 计算经纬度对应的行列号
		double col = (longitude - geoTransform[0]) / geoTransform[1];
		double row = (latitude - geoTransform[3]) / geoTransform[5];

		// 检查行列号是否在有效范围内
		if (col < 0 || row < 0 || col >= m_Image->GetRasterXSize() || row >= m_Image->GetRasterYSize())
			return false;

		// 获取周围 4 个像素的行列号
		int col0 = static_cast<int>(col);
		int row0 = static_cast<int>(row);
		int col1 = col0 + 1;
		int row1 = row0 + 1;

		// 边界检查
		col1 = std::min(col1, m_Image->GetRasterXSize() - 1);
		row1 = std::min(row1, m_Image->GetRasterYSize() - 1);

		// 计算插值权重
		double dx = col - col0;
		double dy = row - row0;

		// 读取每个波段的像素值
		Values.clear();
		for (int band = 1; band <= bandCount; band++)
		{
			GDALRasterBand* pBand = m_Image->GetRasterBand(band);
			if (!pBand)
				return false;

			// 读取 4 个像素的值
			double v00, v01, v10, v11;
			CPLErr err00 = pBand->RasterIO(GF_Read, col0, row0, 1, 1, &v00, 1, 1, GDT_Float64, 0, 0, &m_Arg);
			CPLErr err01 = pBand->RasterIO(GF_Read, col0, row1, 1, 1, &v01, 1, 1, GDT_Float64, 0, 0, &m_Arg);
			CPLErr err10 = pBand->RasterIO(GF_Read, col1, row0, 1, 1, &v10, 1, 1, GDT_Float64, 0, 0, &m_Arg);
			CPLErr err11 = pBand->RasterIO(GF_Read, col1, row1, 1, 1, &v11, 1, 1, GDT_Float64, 0, 0, &m_Arg);

			if (err00 != CE_None || err01 != CE_None || err10 != CE_None || err11 != CE_None)
				return false;

			// 双线性插值计算
			double value = (1 - dx) * (1 - dy) * v00 +
				(1 - dx) * dy * v01 +
				dx * (1 - dy) * v10 +
				dx * dy * v11;

			Values.push_back(value);
		}

		return true;
	}

4：Mitchell-Netravali 插值

	/// \brief 根据经纬度读取对应位置的像素值（Mitchell-Netravali 插值）
	/// \param longitude 经度
	/// \param latitude 纬度
	/// \param bandCount 波段数量
	/// \param Values 返回的像素值数组
	virtual bool IO_MitchellNetravali(double longitude, double latitude, int bandCount, std::vector<double>& Values)
	{
		if (!m_Image || bandCount < 1 || bandCount > m_Image->GetRasterCount())
			return false;

		// 获取地理变换参数
		double geoTransform[6];
		if (m_Image->GetGeoTransform(geoTransform) != CE_None)
			return false;

		// 计算经纬度对应的行列号
		double col = (longitude - geoTransform[0]) / geoTransform[1];
		double row = (latitude - geoTransform[3]) / geoTransform[5];

		// 检查行列号是否在有效范围内
		if (col < 0 || row < 0 || col >= m_Image->GetRasterXSize() || row >= m_Image->GetRasterYSize())
			return false;

		// Mitchell-Netravali 插值核函数
		auto mitchell_netravali_kernel = [](double t) -> double {
			t = std::abs(t);
			if (t < 1.0) {
				return (7.0 * t * t * t - 12.0 * t * t + 5.333333) / 6.0;
			}
			else if (t < 2.0) {
				return (-2.333333 * t * t * t + 12.0 * t * t - 20.0 * t + 10.666667) / 6.0;
			}
			return 0.0;
		};

		// 读取每个波段的像素值
		Values.clear();
		for (int band = 1; band <= bandCount; band++)
		{
			GDALRasterBand* pBand = m_Image->GetRasterBand(band);
			if (!pBand)
				return false;

			double sum = 0.0;
			double weight_sum = 0.0;

			// 遍历 4x4 邻域
			for (int i = -1; i <= 2; i++)
			{
				for (int j = -1; j <= 2; j++)
				{
					int x = static_cast<int>(col) + j;
					int y = static_cast<int>(row) + i;

					if (x >= 0 && x < m_Image->GetRasterXSize() && y >= 0 && y < m_Image->GetRasterYSize())
					{
						double value;
						CPLErr err = pBand->RasterIO(GF_Read, x, y, 1, 1, &value, 1, 1, GDT_Float64, 0, 0,&m_Arg);
						if (err == CE_None)
						{
							double dx = col - x;
							double dy = row - y;
							double weight = mitchell_netravali_kernel(dx) * mitchell_netravali_kernel(dy);
							sum += value * weight;
							weight_sum += weight;
						}
					}
				}
			}

			if (weight_sum > 0)
				Values.push_back(sum / weight_sum);
			else
				Values.push_back(0.0); // 如果没有有效像素，填充默认值
		}

		return true;
	}

5：Catmull-Rom 插值

	/// \brief 根据经纬度读取对应位置的像素值（Catmull-Rom 插值）
	/// \param longitude 经度
	/// \param latitude 纬度
	/// \param bandCount 波段数量
	/// \param Values 返回的像素值数组
	virtual bool IO_CatmullRom(double longitude, double latitude, int bandCount, std::vector<double>& Values)
	{
		if (!m_Image || bandCount < 1 || bandCount > m_Image->GetRasterCount())
			return false;

		// 获取地理变换参数
		double geoTransform[6];
		if (m_Image->GetGeoTransform(geoTransform) != CE_None)
			return false;

		// 计算经纬度对应的行列号
		double col = (longitude - geoTransform[0]) / geoTransform[1];
		double row = (latitude - geoTransform[3]) / geoTransform[5];

		// 检查行列号是否在有效范围内
		if (col < 0 || row < 0 || col >= m_Image->GetRasterXSize() || row >= m_Image->GetRasterYSize())
			return false;

		// Catmull-Rom 插值核函数
		auto catmull_rom_kernel = [](double t) -> double {
			t = std::abs(t);
			if (t < 1.0) {
				return 1.0 - 2.5 * t * t + 1.5 * t * t * t;
			}
			else if (t < 2.0) {
				return 2.0 - 4.0 * t + 2.5 * t * t - 0.5 * t * t * t;
			}
			return 0.0;
		};

		// 读取每个波段的像素值
		Values.clear();
		for (int band = 1; band <= bandCount; band++)
		{
			GDALRasterBand* pBand = m_Image->GetRasterBand(band);
			if (!pBand)
				return false;

			double sum = 0.0;
			double weight_sum = 0.0;

			// 遍历 4x4 邻域
			for (int i = -1; i <= 2; i++)
			{
				for (int j = -1; j <= 2; j++)
				{
					int x = static_cast<int>(col) + j;
					int y = static_cast<int>(row) + i;

					if (x >= 0 && x < m_Image->GetRasterXSize() && y >= 0 && y < m_Image->GetRasterYSize())
					{
						double value;
						CPLErr err = pBand->RasterIO(GF_Read, x, y, 1, 1, &value, 1, 1, GDT_Float64, 0, 0, &m_Arg);
						if (err == CE_None)
						{
							double dx = col - x;
							double dy = row - y;
							double weight = catmull_rom_kernel(dx) * catmull_rom_kernel(dy);
							sum += value * weight;
							weight_sum += weight;
						}
					}
				}
			}

			if (weight_sum > 0)
				Values.push_back(sum / weight_sum);
			else
				Values.push_back(0.0); // 如果没有有效像素，填充默认值
		}

		return true;
	}

cejutue
关注
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值