四、C++和Opencv实现的基于PCA主成分分析的高一水体信息提取

#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/highgui.hpp"
#include <gdal.h>
#include <gdal_priv.h>  
#include <iostream>
#include <string>
#include <vector>
#include "opencv/cv.h"  
#include "opencv/highgui.h" 
#include<opencv2/opencv.hpp>  
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
#include <fstream>
#include <sstream>
using namespace std;

using namespace cv;

 

cv::Mat GDAL2Mat(const char* fileName);//转换成opencv能识别的类型
void RemoveSmallRegion(Mat &Src, Mat &Dst, unsigned int AreaLimit, int CheckMode, int NeihborMode);
void Operation(cv::InputArray input1, cv::InputArray input2, cv::OutputArray output);
int main()
{
int w_mean_var = 5;
int w_half = (w_mean_var - 1) / 2;
Mat img;
vector<Mat> band;
const char* filename = "yangben3.tif";
img = GDAL2Mat(filename);
split(img, band);
int channels = img.channels();

if (channels!=4){
cout << "请注意查看影像的波段个数" << endl;
return -1;
}
split(img, band);
Mat b1 = band[0];
Mat b2 = band[1];
Mat b3 = band[2];
Mat b4 = band[3];
int m = b1.rows;

int n = b1.cols;

 

Mat ndvi;
Operation(b4,b3,ndvi);
Mat ndwi;
Operation(b2, b4, ndwi);
double amin = 0, amax = 0;
cv::Point minPt, maxPt;
minMaxLoc(b4, &amin, &amax, &minPt, &maxPt);
Mat ni;
ni= 255.0f*(b4 - amin) / (amax - amin);
Mat nir;
ni.copyTo(nir);

int k = 3;
Mat r = Mat::zeros(k,k,CV_32F);
Mat v = Mat::zeros(1, k, CV_32F);
Mat a = Mat::zeros(1, k, CV_32F);

for (int i = 0; i < m; i++){
for (int j = 0; j < n; j++){
a.at<float>(0, 0) = ndvi.at<float>(i, j);
a.at<float>(0, 1) = ndwi.at<float>(i, j);
a.at<float>(0, 2) = nir.at<float>(i, j);
   Mat at = a.t();
r = r + at*a;
v = v + a;
}
}
r = r / (m*n);
v = v / (n*m);
Mat vt = v.t();
Mat C = r - vt*v;

Mat eValuesMat;
Mat eVectorsMat;
eigen(C, eValuesMat, eVectorsMat);
Mat y;
eValuesMat.copyTo(y);

y.reshape(1, 1).copyTo(y);
cv::sort(y, y, CV_SORT_EVERY_ROW + CV_SORT_ASCENDING);
sortIdx(y, yi, SORT_EVERY_ROW + SORT_ASCENDING);
cout << yi.type() << endl;

Mat newy;
eValuesMat.copyTo(newy);
newy.reshape(1, 1).copyTo(newy);
cv::sort(newy, newy, CV_SORT_EVERY_ROW + CV_SORT_DESCENDING);

Mat mb(1, 1, CV_32FC1, Scalar(0));
reduce(y, mb, 1, CV_REDUCE_SUM);
float ysum = mb.at<float>(0,0);
Mat rate = y / ysum;
Mat newrate = newy / ysum;
float sumrate = 0;
Mat newi=Mat::zeros(1,3,CV_32F);
int len = 0;
for (int k = 2; k >= 0; k--){
sumrate = sumrate + rate.at<float>(0,k);
newi.at<float>(0, 2 - k) = yi.at<float>(0,k);
len++;
if (sumrate > 0.85)break; 
}
Mat NT = Mat::zeros(3,1,CV_32F);
for (int p = 0; p < len; p++){
for (int q = 0; q < 3; q++){
int l =(int)(newi.at<float>(0,p));
NT.at<float>(q, p) = eVectorsMat.at<float>(q,l);
}
}
Mat aa = Mat::zeros(1, k, CV_32F);
Mat TT = NT.t();
Mat RN = Mat::zeros(m,n,CV_64F);
    for (int i = 0; i < m;i++){
for (int j = 0; j < n; j++){
RN.at<double>(i, j) = ndvi.at<float>(i, j)*TT.at<float>(0, 0) + ndwi.at<float>(i, j)*TT.at<float>(0, 1) + nir.at<float>(i, j)*TT.at<float>(0, 2);
}
    }
Mat gray = Mat::zeros(m,n,CV_8U);
for (int i = 0; i < m;i++){
for (int j = 0; j < n;j++){
if (RN.at<double>(i, j) < (-163)){
gray.at<uchar>(i, j) = 255;
}
else{
gray.at<uchar>(i, j) = 0;
}
}

}

 

Mat Gray;
gray.copyTo(Gray);
RemoveSmallRegion(Gray, Gray, 600, 1, 1);
RemoveSmallRegion(Gray, Gray, 100, 0, 0);
namedWindow("image", WINDOW_NORMAL);
imshow("image", Gray);
waitKey(0);
img.release();
return 0;
}




void Operation(cv::InputArray input1, cv::InputArray input2, cv::OutputArray output){
Mat b1 = input1.getMat();
Mat b2 = input2.getMat();
Mat up,down,bz;
up=b1-b2;
down=b1+b2;
bz=up/down;
double amin = 0, amax = 0;
cv::Point minPt, maxPt;
minMaxLoc(bz, &amin, &amax, &minPt, &maxPt);
Mat gray;
gray = 255.0f*(bz - amin) / (amax - amin);
gray.copyTo(output);
}
cv::Mat GDAL2Mat(const char* fileName)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
GDALDataset * poDataset;
poDataset = (GDALDataset*)GDALOpen(fileName, GA_ReadOnly); 
int Cols = poDataset->GetRasterXSize();
int Rows = poDataset->GetRasterYSize(); 
int BandSize = poDataset->GetRasterCount();
double *adfGeoTransform = new double[6];
poDataset->GetGeoTransform(adfGeoTransform);


std::vector <cv::Mat> imgMat;  
float* pafScan = (float*)CPLMalloc(sizeof(float)* Cols);
for (int i = 0; i< BandSize; i++)
{
GDALRasterBand *pBand = poDataset->GetRasterBand(i + 1);
pafScan = new float[Cols*Rows];
pBand->RasterIO(GF_Read, 0, 0, Cols, Rows, pafScan, Cols, Rows, GDT_Float32, 0, 0);
cv::Mat A = cv::Mat(Rows, Cols, CV_32FC1, pafScan);
imgMat.push_back(A.clone());
delete[]pBand;
A.release();
}
delete[] pafScan;
cv::Mat img;
img.create(Rows, Cols, CV_32FC(BandSize));
merge(imgMat, img);
//GDALClose((GDALDatasetH)poDataset);
imgMat.clear();
return img;
}
void RemoveSmallRegion(Mat &Src, Mat &Dst, unsigned int AreaLimit, int CheckMode, int NeihborMode)
{
int RemoveCount = 0;  
Mat PointLabel = Mat::zeros(Src.size(), CV_8UC1);
if (CheckMode == 1) 
{
cout << "去除小连通域.";
for (int i = 0; i < Src.rows; i++)
{
for (int j = 0; j < Src.cols; j++)
{
if (Src.at<uchar>(i, j) < 10)
{
PointLabel.at<uchar>(i, j) = 3;
}
}
}
}
else
{
cout << "去除孔洞";
for (int i = 0; i < Src.rows; i++)
{
for (int j = 0; j < Src.cols; j++)
{
if (Src.at<uchar>(i, j) > 10)
{
PointLabel.at<uchar>(i, j) = 3;//如果原图是白色区域，标记为合格，像素为3  
}
}
}
}
vector<Point2i>NeihborPos;
NeihborPos.push_back(Point2i(-1, 0));
NeihborPos.push_back(Point2i(1, 0));
NeihborPos.push_back(Point2i(0, -1));
NeihborPos.push_back(Point2i(0, 1));
if (NeihborMode == 1)
{
cout << "Neighbor mode: 8邻域." << endl;
NeihborPos.push_back(Point2i(-1, -1));
NeihborPos.push_back(Point2i(-1, 1));
NeihborPos.push_back(Point2i(1, -1));
NeihborPos.push_back(Point2i(1, 1));
}
else cout << "Neighbor mode: 4邻域." << endl;
int NeihborCount = 4 + 4 * NeihborMode;
int CurrX = 0, CurrY = 0;
for (int i = 0; i < Src.rows; i++)
{
for (int j = 0; j < Src.cols; j++)
{
if (PointLabel.at<uchar>(i, j) == 0)
{   
vector<Point2i>GrowBuffer;
GrowBuffer.push_back(Point2i(j, i));
PointLabel.at<uchar>(i, j) = 1;
int CheckResult = 0;
for (unsigned int z = 0; z < GrowBuffer.size(); z++)
{
for (int q = 0; q < NeihborCount; q++)
{
CurrX = GrowBuffer.at(z).x + NeihborPos.at(q).x;
CurrY = GrowBuffer.at(z).y + NeihborPos.at(q).y;
if (CurrX >= 0 && CurrX<Src.cols&&CurrY >= 0 && CurrY<Src.rows)  
{
if (PointLabel.at<uchar>(CurrY, CurrX) == 0)
{
GrowBuffer.push_back(Point2i(CurrX, CurrY));    
PointLabel.at<uchar>(CurrY, CurrX) = 1;          
}
}
}
}
if (GrowBuffer.size()>AreaLimit)   
CheckResult = 2;
else
{
CheckResult = 1;
RemoveCount++;
}
for (unsigned int z = 0; z < GrowBuffer.size(); z++)
{
CurrX = GrowBuffer.at(z).x;
CurrY = GrowBuffer.at(z).y;
PointLabel.at<uchar>(CurrY, CurrX) += CheckResult;
}
}
}
}
CheckMode = 255 * (1 - CheckMode);
for (int i = 0; i < Src.rows; ++i)
{
for (int j = 0; j < Src.cols; ++j)
{
if (PointLabel.at<uchar>(i, j) == 2)
{
Dst.at<uchar>(i, j) = CheckMode;
}
else if (PointLabel.at<uchar>(i, j) == 3)
{
Dst.at<uchar>(i, j) = Src.at<uchar>(i, j);

}
}
}
cout << RemoveCount << " objects removed." << endl;

}

结果如下图：可以看到后面用c++实现的水体信息提取没有在matlab中实现的效果好