这个...时隔很多天,终于又来更新了。
今天是图像融合中的Brovey方法。先将理论,后讲实现。
一、Brovey图像融合
这是Brovey图像融合的公式,我们单独拎出一个来看看:
讲一下:图像融合后新的R值 = 图像融合前旧的R值 x 全色波段对应位置的值 ÷ 旧的RGB和
对于另外的新的G、B值也是一样的。
好,Brovey的理论就只有这么多。
二、代码实现
提前说一句,在以下讲解中,我们的已有变量是:
Pan_dt: 全色波段的dataset
RGB_dt: 多光谱波段的dataset
Bands_idx: 从多光谱的多个波段中获取的3个波段的波段索引,长度为3的int数组
第一步,先从两个dataset中获取相同大小的数据。在图像融合中,我们输入的是全色波段和多光谱波段,意要把两幅图像融合为更清晰的多光谱图像。所以读取dataset的时候我们获取的数据长度就应该是 全色dataset的宽 x 全色dataset的高。代码如下:
int Xsize = this.Pan_dt.RasterXSize, Ysize = this.Pan_dt.RasterYSize;
int[] Pan = new int[Xsize * Ysize];
int[] R = new int[Xsize * Ysize];
int[] G = new int[Xsize * Ysize];
int[] B = new int[Xsize * Ysize];
this.Pan_dt.GetRasterBand(1).ReadRaster(0, 0, Xsize, Ysize, Pan, Xsize, Ysize, 0, 0);
this.RGB_dt.GetRasterBand(this.Bands_idx[0]).ReadRaster(0, 0, this.RGB_dt.RasterXSize, this.RGB_dt.RasterYSize, R, Xsize, Ysize, 0, 0);
this.RGB_dt.GetRasterBand(this.Bands_idx[1]).ReadRaster(0, 0, this.RGB_dt.RasterXSize, this.RGB_dt.RasterYSize, G, Xsize, Ysize, 0, 0);
this.RGB_dt.GetRasterBand(this.Bands_idx[2]).ReadRaster(0, 0, this.RGB_dt.RasterXSize, this.RGB_dt.RasterYSize, B, Xsize, Ysize, 0, 0);
int[][] PRGB_data = new int[4][] { Pan, R, G, B };
我们可以看到这里获得了一个名为PRGB_data的二维数组,第一层的含义是 [Pan波段、R波段、G波段、B波段] ,其中每一个波段(第二层)是一个长度为 Xsize*Ysize 的一维数组。
提取出来后就可以进行第二步,Brovey的计算,其实也很简单,代码如下:
for (int i = 1; i <= 3; i++)
{
for (int y = 0; y < this.Pan_dt.RasterYSize; y++)
for (int x = 0; x < this.Pan_dt.RasterXSize; x++)
{
PRGB_data[i][x + y * Pan_dt.RasterXSize] = (int)System.Math.Round(PRGB_data[0][x + y * Pan_dt.RasterXSize] *
((double)PRGB_data[i][x + y * Pan_dt.RasterXSize] /
(double)(PRGB_data[1][x + y * Pan_dt.RasterXSize] + PRGB_data[2][x + y * Pan_dt.RasterXSize] + PRGB_data[3][x + y * Pan_dt.RasterXSize])));
}
}
可以看到就是做一个个循环然后套公式就可以了。这一步做完了之后PRGB_data的RGB波段(第一层索引为1,2,3)就已经是更新完毕了,后面的第三步就是把这三个数组的数据组合成一个TIFF输出出去。
在C# GDAL 数字图像处理Part9 辐射定标与输出Dataset_A_RSswjtuer的博客-CSDN博客中我们已经提到过了输出dataset。
我们知道,从 dataset 中 ReadRaster 时,是读出了一个指定长度的数组,那我们现在对这个数组修改完了,想要把它塞回去,怎么办呢?在这里,就不用像Part9一样一个个像素地塞回去了,代码如下:
for (int i = 1; i <= OutputDt.RasterCount; i++)
{
OutputDt.GetRasterBand(i).WriteRaster(0, 0, OutputDt.RasterXSize, OutputDt.RasterYSize, PRGB_data[i], OutputDt.RasterXSize, OutputDt.RasterYSize, 0, 0);//整段写入
}
OutputDt.FlushCache();
可以看到直接用WriteRaster直接整个写进去就OK。
今天就到这~