利用VC++编写计算植被指数的程序
（1）采用参数传递的方式，并使用一个变量输出结果
（2）要有说明部分，说明输入参数的格式，说明输出数据的方式
（3）将计算结果用文件保存下来，命名为vegetable.tif
众所周知，近红外（NIR）的波段为0.78到3微米，而红光（R）的波段为0.625到0.74微米。如图1，在植物的反射波谱特性曲线中，近红外和红光波段中间有较大的落差，这是植物在反射波谱特性曲线中特有的性质，利用这个性质可以在遥感图像中提取植被。
定义归一化植被指数（NDVI）为：
N
D
V
I
=
N
I
R
−
R
N
I
R
+
R
NDVI= {{NIR-R}\over{NIR+R}}
NDVI=NIR+RNIR−R
NDVI越大，表示近红外和红外波段之间的落差越大，则越有可能是植被。通过计算遥感影像的归一化植被指数，再均匀地映射到[0,255]，遥感影像上的像元越白，则越有可能是植被。
首先创建ImageProcessing类，有六个私有成员，分别为：Dataset，存储读入影像的数据集；xSize、ySize、imageSize、bandNum、dataType分别存储读入影像的列数、行数、像元数、波段数、DN值的数据类型。
private:
GDALDataset* Dataset;
int xSize;
int ySize;
int imageSize;
int bandNum;
GDALDataType dataType;
然后添加公有成员，首先是构造函数和析构函数。构造函数通过传入指向待处理图像路径的字符指针来构造实例。
ImageProcessing(const char*);
~ImageProcessing();
然后创建两个函数calcuNDVI和saveImage分别用来计算NDVI值和保存影像。calcuNDVI函数传入近红外的波段数和红外的波段数，返回一个一维数组，其中的元素是将所有像元的NDVI值由[-1,1]平均映射到[0,255]后的值。saveImage函数传入保存的路径和calcuNDVI返回的一维数组。
unsigned char* calcuNDVI(int, int) const;
void saveImage(const char*, unsigned char*);
在计算NDVI之前，先创建一个函数NDVI2DN，这个函数可以将[-1,1]的值均匀地映射到[0,255]
int NDVI2DN(float);
首先创建四个一维数组，分别用来存储近红外波段的DN值、红外波段的DN值、[-1,1]的NDVI、[0,255]的NDVI
unsigned char* NIRBuffer = new unsigned char[imageSize * sizeof(dataType)];
unsigned char* RBuffer = new unsigned char[imageSize * sizeof(dataType)];
float* NDVIFloatBuffer = new float[imageSize];
unsigned char* NDVIBuffer = new unsigned char[imageSize * sizeof(dataType)];
然后用RasterIO函数将红外波段和近红外波段读入数组
GDALRasterBand* NIRBand = Dataset->GetRasterBand(NIRBandnum);
NIRBand->RasterIO(GF_Read, 0, 0, xSize, ySize, NIRBuffer, xSize, ySize, dataType, 0, 0);
GDALRasterBand* RBand = Dataset->GetRasterBand(RBandnum);
RBand->RasterIO(GF_Read, 0, 0, xSize, ySize, RBuffer, xSize, ySize, dataType, 0, 0);
然后逐像元计算NDVI，然后映射到[0,255]，最后存入NDVIBuffer，删除堆内存中的数组，将NDVIBuffer返回，即完成了遥感影像植被指数的计算。
for(int i = 0; i < ySize; i++) {
for(int j = 0; j < xSize; j++) {
int index = xSize * i + j;
if(int(NIRBuffer[index] + RBuffer[index]) == 0) {
NDVIFloatBuffer[index] = 0;
}
else {
NDVIFloatBuffer[index] = float(NIRBuffer[index] - RBuffer[index]) / float(NIRBuffer[index] + RBuffer[index]);
}
NDVIBuffer[index] = NDVI2DN(NDVIFloatBuffer[index]);
}
}
delete[] NDVIFloatBuffer;
delete[] NIRBuffer;
delete[] RBuffer;
return NDVIBuffer;
最后将影像保存即可。
使用示例：
首先注册GDAL，然后定义待处理图片路径、保存路径。接着创建ImageProcessing类的实例，调用saveImage方法，在saveImage方法中传出参数，最后删除实例释放内存。
int main() {
GDALAllRegister();
const char* imagepath = "..\\image\\test2.tif";
const char* savepath = "..\\image\\vegetable.tif";
ImageProcessing* imageProcessing = new ImageProcessing(imagepath);
//波段4是近红外，波段3是红外
imageProcessing->saveImage(savepath, imageProcessing->calcuNDVI(4,3));
delete imageProcessing;
return 0;
}
归一化植被指数均匀映射到[0,255]之后输出的遥感影像如图2所示：
#pragma once
#include<iostream>
#include<gdal.h>
#include<gdal_priv.h>
using namespace std;
class ImageProcessing {
private:
GDALDataset* Dataset;
int xSize;
int ySize;
int imageSize;
int bandNum;
GDALDataType dataType;
public:
ImageProcessing(const char*);
~ImageProcessing();
int getBandNum() const;
int getImageSize() const;
void calcuNDVI(const char*, int, int) const;
};
#include<iostream>
#include"ImageProcessing.h"
#include<gdal.h>
#include<gdal_priv.h>
#include<assert.h>
using namespace std;
int NDVI2DN(float);
ImageProcessing::ImageProcessing(const char* path) {
this->Dataset = (GDALDataset*)GDALOpen(path, GA_ReadOnly);
this->xSize = Dataset->GetRasterXSize();
this->ySize = Dataset->GetRasterYSize();
this->imageSize = this->xSize * this->ySize;
this->bandNum = Dataset->GetRasterCount();
this->dataType = Dataset->GetRasterBand(1)->GetRasterDataType();
}
ImageProcessing::~ImageProcessing() {}
int ImageProcessing::getBandNum() const {
return this->bandNum;
}
int ImageProcessing::getImageSize() const {
return this->imageSize;
}
void ImageProcessing::calcuNDVI(const char* savepath, int NIRBandnum, int RBandnum) const {
unsigned char* NIRBuffer = new unsigned char[imageSize * sizeof(dataType)];
unsigned char* RBuffer = new unsigned char[imageSize * sizeof(dataType)];
float* NDVIFloatBuffer = new float[imageSize];
unsigned char* NDVIBuffer = new unsigned char[imageSize * sizeof(dataType)];
GDALRasterBand* NIRBand = Dataset->GetRasterBand(NIRBandnum);
NIRBand->RasterIO(GF_Read, 0, 0, xSize, ySize, NIRBuffer, xSize, ySize, dataType, 0, 0);
GDALRasterBand* RBand = Dataset->GetRasterBand(RBandnum);
RBand->RasterIO(GF_Read, 0, 0, xSize, ySize, RBuffer, xSize, ySize, dataType, 0, 0);
for(int i = 0; i < ySize; i++) {
for(int j = 0; j < xSize; j++) {
int index = xSize * i + j;
if(int(NIRBuffer[index] + RBuffer[index]) == 0) {
NDVIFloatBuffer[index] = 0;
}
else {
NDVIFloatBuffer[index] = float(NIRBuffer[index] - RBuffer[index]) / float(NIRBuffer[index] + RBuffer[index]);
}
NDVIBuffer[index] = NDVI2DN(NDVIFloatBuffer[index]);
}
}
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, 1, this->dataType, papszOption);
if(!saveDataset) {
assert(!saveDataset);
}
saveDataset->RasterIO(GF_Write, 0, 0, this->xSize, this->ySize, NDVIBuffer, this->xSize, this->ySize, this->dataType, 1, BandMap, 0, 0, 0);
GDALClose(Dataset);
GDALClose(saveDataset);
delete[] NDVIFloatBuffer;
delete[] NDVIBuffer;
delete[] NIRBuffer;
delete[] RBuffer;
}
int NDVI2DN(float NDVI) {
NDVI += 1;
int DN;
if(NDVI == 1) {
DN = 127;
}
else {
NDVI = NDVI / 2 * 255.0;
DN = int(NDVI);
}
return DN;
}
#include<iostream>
#include<gdal.h>
#include<gdal_priv.h>
#include"ImageProcessing.h"
using namespace std;
int main() {
GDALAllRegister();
const char* imagepath = "..\\image\\test2.tif";
const char* savepath = "..\\image\\vegetable.tif";
ImageProcessing* imageProcessing = new ImageProcessing(imagepath);
//cout << imageProcessing->getBandNum() << endl;
//cout << imageProcessing->getImageSize();
//波段4是近红外，波段3是红外
imageProcessing->calcuNDVI(savepath, 4, 3);
delete imageProcessing;
return 0;
}