一、實驗任務與內容

利用VC++撰寫一個3*3的均值濾波或中值濾波程式

（1）采用引數傳遞方式，并用一個變數輸出運算結果

（2）要有說明部分，說明輸入引數的格式，說明輸出資料方式

（3）將濾波后的結果用檔案保存下來，命名為Noise_result.tif

二、實驗原理

2.1 均值濾波

在遙感影像上選擇一個3×3的矩陣，如下式所示：
[ a 1 a 2 a 3 b 1 b 2 b 3 c 1 c 2 c 3 ] \left[\begin {matrix} a_1&a_2&a_3\\ b_1&b_2&b_3\\ c_1&c_2&c_3 \end{matrix}\right] ???a1?b1?c1??a2?b2?c2??a3?b3?c3?????
令：
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的矩陣，如下式所示：
[ a 1 a 2 a 3 b 1 b 2 b 3 c 1 c 2 c 3 ] \left[\begin {matrix}a_1&a_2&a_3\\b_1&b_2&b_3\\c_1&c_2&c_3\end{matrix}\right] ???a1?b1?c1??a2?b2?c2??a3?b3?c3?????
將9個元素按從大到小或者從小到大排列：
[ x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 , x 8 , x 9 ] = s o r t ( a 1 , a 2 , a 3 , b 1 , b 2 , b 3 , c 1 , c 2 , c 3 ) [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) [x1?,x2?,x3?,x4?,x5?,x6?,x7?,x8?,x9?]=sort(a1?,a2?,a3?,b1?,b2?,b3?,c1?,c2?,c3?)
令：
b 2 = x 5 b_2=x_5 b2?=x5?
遍歷影像的所有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 - 1)) + j + k % 3 - 1];
		}
		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 - 1)) + j + k % 3 - 1];
			}
			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;
}