数字图像处理：遥感图像的中值滤波与均值滤波

一、实验任务与内容

利用VC++编写一个3*3的均值滤波或中值滤波程序

（1）采用参数传递方式，并用一个变量输出运算结果

（2）要有说明部分，说明输入参数的格式，说明输出数据方式

（3）将滤波后的结果用文件保存下来，命名为Noise_result.tif

二、实验原理

2.1 均值滤波

在遥感图像上选择一个3×3的矩阵，如下式所示：
$$\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]$$
令：
b 2 = m e a n { a 1 + a 2 + a 3 + b 1 + b 2 + b 3 + c 1 + c 2 + c 3 } b_2=mean\{a_1+a_2+a_3+b_1+b_2+b_3+c_1+c_2+c_3\} b2​=mean{a1​+a2​+a3​+b1​+b2​+b3​+c1​+c2​+c3​}
遍历影像所有的3*3子阵，即进行了均值滤波

2.2 中值滤波

在遥感图像上选择一个3×3的矩阵，如下式所示：
$$\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]$$
将9个元素按从大到小或者从小到大排列：
$$[x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9]=sort(a_1,a_2,a_3,b_1,b_2,b_3,c_1,c_2,c_3)$$
令：
$$b_2=x_5$$
遍历影像的所有3*3子阵，即进行了中值滤波

三、编程思路

首先创建Filter类，其含有五个私有字段，分别为：Dataset，用于存储读入的图像数据集；Xsize、Ysize、Bandnum、dataType，分别用于储存读入图像的列数、行数、波段数和数据类型。

private:
	GDALDataset* Dataset;
	int Xsize;
	int Ysize;
	int Bandnum;
	GDALDataType dataType;

然后创建公有的构造函数和析构函数，构造函数中传入指向图像路径的字符指针，便于通过GDAL库读取图像。

public:
	Filter(const char*);
	~Filter();

然后创建公有的方法MeanFilter和MedianFilter，并传入指向保存路径的字符指针，便于通过GDAL库保存图像

void MeanFilter(const char*) const;
void MedianFilter(const char*) const;

以中值滤波为例，首先创建两个数组，用于储存读入的影像数据和处理后的影像数据：

unsigned char* OldBuf = new unsigned char[Xsize * Ysize * dataType];
unsigned char* NewBuf = new unsigned char[Xsize * Ysize * dataType];

然后将用GDAL读入的数据存入OldBuf中：

GDALRasterBand* band1 = Dataset->GetRasterBand(1);
band1->RasterIO(GF_Read, 0, 0, Xsize, Ysize, OldBuf, Xsize, Ysize, dataType, 0, 0);

由于使用的是待处理像元的八邻域，所以最外圈的像元不做处理，直接复制到NewBuf中去：

for(int i = 0; i < this->Ysize; i++) {
	for(int j = 0; j < this->Xsize; j++) {
		if(i == 0 || j == 0 || i == Ysize - 1 || j == Xsize - 1) {
			//影像中第i行第j列的像素在一维数组中的索引是Xsize*(i-1)+j
			NewBuf[Xsize * (i - 1) + j] = OldBuf[Xsize * (i - 1) + j];
		}
	}
}

创建一个长度为9的一维数组用于储存待处理像元及其八邻域的像素值：

int temp[9] = { 0 };//3*3的

遍历除最外圈像元以外的所有像元，将待处理像元的DN值与其八邻域的DN值存入temp数组，再将数组排序，取数组最中间的值，即 *(temp+4)，将其赋值给待处理像元，即完成了中值滤波。若将DN值存入数组后，将数组元素的平均值赋值给待处理像元，则完成了均值滤波。

for(int i = 1; i < this->Ysize - 1; i++) {
	for(int j = 1; j < this->Xsize - 1; j++) {
		for(int k = 0; k < 9; k++) {
			*(temp + k) = OldBuf[Xsize * (i + k / 3) + j + k % 3];
		}
		sort(temp, temp + 9);
		NewBuf[Xsize * i + j] = (unsigned char)*(temp + 4);
	}
}

最后将滤波后的影像保存即可。

使用方法：首先要注册GDAL，然后创建Filter类的实例，再调用方法即可：

int main() {
	GDALAllRegister();

	const char* imagepath = ".\\image\\nosieImg.tif";
	const char* savepath = ".\\image\\Nosie_result1.tif";

	Filter *filter=new Filter(imagepath);
	filter->MeanFilter(savepath);
	//filter->MedianFilter(savepath);
	delete filter;
	return 0;
}

四、效果展示

滤波前的影像：

中值滤波后的影像：

均值滤波后的影像：

五、源代码

Filter.h

#pragma once
#include<iostream>
#include<gdal_priv.h>
#include<gdal.h>

using namespace std;

class Filter {
private:
	GDALDataset* Dataset;
	int Xsize;
	int Ysize;
	int Bandnum;
	GDALDataType dataType;
public:
	Filter(const char*);
	~Filter();
	int getXsize() const;
	int getBandnum() const;
	void MeanFilter(const char*) const;
	void MedianFilter(const char*) const;
};

Filter.cpp

#include "Filter.h"
#include<iostream>
#include<gdal_priv.h>
#include<gdal.h>
#include<assert.h>
#include<algorithm>

using namespace std;


Filter::Filter(const char* path) {
	this->Dataset = (GDALDataset*)GDALOpen(path, GA_ReadOnly);
	this->Xsize = Dataset->GetRasterXSize();
	this->Ysize = Dataset->GetRasterYSize();
	this->Bandnum = Dataset->GetRasterCount();
	this->dataType = Dataset->GetRasterBand(1)->GetRasterDataType();
}

Filter::~Filter() {}

int Filter::getXsize() const {
	return this->Xsize;
}

int Filter::getBandnum() const {
	return this->Bandnum;
}

void Filter::MeanFilter(const char* savepath) const {
	unsigned char* OldBuf = new unsigned char[Xsize * Ysize * dataType];
	unsigned char* NewBuf = new unsigned char[Xsize * Ysize * dataType];

	GDALRasterBand* band1 = Dataset->GetRasterBand(1);
	band1->RasterIO(GF_Read, 0, 0, Xsize, Ysize, OldBuf, Xsize, Ysize, dataType, 0, 0);

	//外围的一圈不修改，所以直接赋值过去
	for(int i = 0; i < this->Ysize; i++) {
		for(int j = 0; j < this->Xsize; j++) {
			if(i == 0 || j == 0 || i == Ysize - 1 || j == Xsize - 1) {
				//影像中第i行第j列的像素在一维数组中的索引是Xsize*(i-1)+j
				NewBuf[Xsize * (i - 1) + j] = OldBuf[Xsize * (i - 1) + j];
			}
		}
	}

	/*均值滤波实现*/
	int temp[9] = { 0 };//3*3的
	for(int i = 1; i < this->Ysize - 1; i++) {
		for(int j = 1; j < this->Xsize - 1; j++) {
			for(int k = 0; k < 9; k++) {
				*(temp + k) = OldBuf[Xsize * (i + k / 3) + j + k % 3];
			}
			int sum = 0;
			for(int t = 0; t < 9; t++) {
				sum += *(temp + t);
			}

			NewBuf[Xsize * i + j] = sum / 9;
		}
	}

	/*保存滤波之后的影像*/
	GDALDriver* Driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	int BandMap[1] = { 1 };
	char** papszOption = nullptr;
	papszOption = CSLSetNameValue(papszOption, "INTERLEAVE", "BAND");
	GDALDataset* saveDataset = Driver->Create(savepath, this->Xsize, this->Ysize, this->Bandnum, this->dataType, papszOption);

	if(!saveDataset) {
		assert(!saveDataset);
	}

	saveDataset->RasterIO(GF_Write, 0, 0, this->Xsize, this->Ysize, NewBuf, this->Xsize, this->Ysize, this->dataType, this->Bandnum, BandMap, 0, 0, 0);

	GDALClose(Dataset);
	GDALClose(saveDataset);
}

void Filter::MedianFilter(const char* savepath) const {
	unsigned char* OldBuf = new unsigned char[Xsize * Ysize * dataType];
	unsigned char* NewBuf = new unsigned char[Xsize * Ysize * dataType];

	GDALRasterBand* band1 = Dataset->GetRasterBand(1);
	band1->RasterIO(GF_Read, 0, 0, Xsize, Ysize, OldBuf, Xsize, Ysize, dataType, 0, 0);

	//外围的一圈不修改，所以直接赋值过去
	for(int i = 0; i < this->Ysize; i++) {
		for(int j = 0; j < this->Xsize; j++) {
			if(i == 0 || j == 0 || i == Ysize - 1 || j == Xsize - 1) {
				//影像中第i行第j列的像素在一维数组中的索引是Xsize*(i-1)+j
				NewBuf[Xsize * (i - 1) + j] = OldBuf[Xsize * (i - 1) + j];
			}
		}
	}

	/*中值滤波实现*/
	int temp[9] = { 0 };//3*3的
	for(int i = 1; i < this->Ysize - 1; i++) {
		for(int j = 1; j < this->Xsize - 1; j++) {
			for(int k = 0; k < 9; k++) {
				*(temp + k) = OldBuf[Xsize * (i + (k / 3 - 1)) + j + k % 3 - 1];
			}
			sort(temp, temp + 9);
			NewBuf[Xsize * i + j] = (unsigned char)*(temp + 4);
		}
	}

	/*保存滤波之后的影像*/
	GDALDriver* Driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	int BandMap[1] = { 1 };
	char** papszOption = nullptr;
	papszOption = CSLSetNameValue(papszOption, "INTERLEAVE", "BAND");
	GDALDataset* saveDataset = Driver->Create(savepath, this->Xsize, this->Ysize, this->Bandnum, this->dataType, papszOption);

	if(!saveDataset) {
		assert(!saveDataset);
	}

	saveDataset->RasterIO(GF_Write, 0, 0, this->Xsize, this->Ysize, NewBuf, this->Xsize, this->Ysize, this->dataType, this->Bandnum, BandMap, 0, 0, 0);

	GDALClose(Dataset);
	GDALClose(saveDataset);
}

main.cpp

#include<gdal_priv.h>
#include<gdal.h>
#include<iostream>
#include"Filter.h"

using namespace std;

int main() {
	GDALAllRegister();

	const char* imagepath = ".\\image\\nosieImg.tif";
	const char* savepath = ".\\image\\Nosie_result1.tif";

	Filter *filter=new Filter(imagepath);
	//filter->MeanFilter(savepath);
	filter->MedianFilter(savepath);
	delete filter;
	return 0;
}