C# GDAL 数字图像处理Part11 Brovey图像融合

最新推荐文章于 2024-07-10 10:54:49 发布

A_RSswjtuer

最新推荐文章于 2024-07-10 10:54:49 发布

阅读量1.1k

点赞数 2

分类专栏：数字图像处理文章标签： c# 开发语言图像处理

本文链接：https://blog.csdn.net/im_nobody33/article/details/128103646

版权

数字图像处理专栏收录该内容

13 篇文章 14 订阅

订阅专栏

这个...时隔很多天，终于又来更新了。

今天是图像融合中的Brovey方法。先将理论，后讲实现。

一、Brovey图像融合

Brovey图像融合方法+Python代码_煎糕丶的博客-CSDN博客_brovey融合

这是Brovey图像融合的公式，我们单独拎出一个来看看：

$R = \frac{r\times Pan}{r+g+b}$

讲一下：图像融合后新的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。

今天就到这~