一、实验任务与内容
利用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;
}