基于Eigen的遙感多光譜影像主成分分析

主成分分析簡介

在多元統(tǒng)計(jì)分析中凶赁,主成分分析(英語:Principal components analysis魁淳,PCA)是一種分析盆驹、簡化數(shù)據(jù)集的技術(shù)。主成分分析經(jīng)常用于減少數(shù)據(jù)集的維數(shù)聊记,同時(shí)保持?jǐn)?shù)據(jù)集中的對(duì)方差貢獻(xiàn)最大的特征撒妈。這是通過保留低階主成分,忽略高階主成分做到的排监。這樣低階成分往往能夠保留住數(shù)據(jù)的最重要方面狰右。

計(jì)算步驟

輸入:n維樣本集D=(x(1),x(2),...,x(m)),要降維到的維數(shù)n.

輸出:降維后的樣本集D′    1) 對(duì)所有的樣本進(jìn)行中心化:        2) 計(jì)算樣本的協(xié)方差矩陣X*XT

3) 對(duì)矩陣X*XT進(jìn)行特征值分解

4)取出最大的n個(gè)特征值(從大到小排序)對(duì)應(yīng)的特征向量(w1,w2,...,wn), 將所有的特征向量標(biāo)準(zhǔn)化后舆床,組成特征向量矩陣W棋蚌。

5)對(duì)樣本集中的每一個(gè)樣本x(i),轉(zhuǎn)化為新的樣本z(i)=WT*x(i)    6) 得到輸出樣本集D′=(z(1),z(2),...,z(m)) 即前m成分

博客園有一篇博客詳細(xì)介紹了PCA 原理

代碼實(shí)現(xiàn)

主要用到了gdal和Eigen庫

gdal用于讀寫遙感多光譜影像

Eigen則便于各種矩陣運(yùn)算

#include"gdal_priv.h"
#include"cpl_conv.h" // for CPLMalloc()
#include<iostream>
#include<fstream>
#include<string>
#include"Eigen/Dense"
using namespace std;
using namespace Eigen;

//eigen實(shí)現(xiàn)主成分分析
void featurenormalize(MatrixXd &X)
{
    //計(jì)算每一維度均值
    MatrixXd meanval = X.colwise().mean();
    RowVectorXd meanvecRow = meanval;
    //樣本均值化為0
    X.rowwise() -= meanvecRow;
}
void computeCov(MatrixXd &X, MatrixXd &C)
{
    //計(jì)算協(xié)方差矩陣C = XTX / n-1;
    C = X.adjoint() * X;
    C = C.array() / (X.rows() - 1);
}
void computeEig(MatrixXd &C, MatrixXd &vec, MatrixXd &val)
{
    //計(jì)算特征值和特征向量嫁佳,使用selfadjont按照對(duì)陣矩陣的算法去計(jì)算,可以讓產(chǎn)生的vec和val按照有序排列(默認(rèn)從大到泄饶骸)
    SelfAdjointEigenSolver<MatrixXd> eig(C);

    vec = eig.eigenvectors();
    val = eig.eigenvalues();
}
int computeDim(MatrixXd &val)
{
    //輸出信息量達(dá)到95%的前n主成分
    /*int dim;
    double sum = 0;
    for (int i = val.rows() - 1; i >= 0; --i)
    {
    sum += val(i, 0);
    dim = i;

    if (sum / val.sum() >= 0.95)
    break;
    }
    return val.rows() - dim;*/
    return 7;//這里設(shè)置輸出7個(gè)主成分
}

void writePcaImg(const char* path, int width, int height, double *pBuff, double *adfGeo, const char *prj, int bandNum, int imageSize, int pcaInd)
{
    GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); //圖像驅(qū)動(dòng)
    char** ppszOptions = NULL;
    int depth = 8;//圖像位深
    int dim = 1;//每個(gè)圖像波段數(shù)蒿往,這里將每個(gè)主成分存儲(chǔ)到一個(gè)單波段圖像
    GDALDataset* dst = pDriver->Create(path, width, height, dim, GDT_Float64, ppszOptions);//創(chuàng)建圖像
    if (dst == nullptr)
        printf("Can't Write Image!");
    dst->SetGeoTransform(adfGeo);//設(shè)置坐標(biāo)
    dst->SetProjection(prj);//設(shè)置投影
    dst->RasterIO(GF_Write, 0, 0, width, height, &pBuff[(bandNum - pcaInd)*imageSize], width, height,
        GDT_Float64, dim, nullptr, dim*depth, width*dim*depth, depth);//寫入圖像
    GDALClose(dst);
}
int main(int argc, char *argv[])
{   //讀取影像
    char* pszFilename = "D:/gdalData/pca/before.img";
    char *outPath = "D:/pca_temp/pca";
    GDALDataset  *poDataset;
    GDALAllRegister();
    poDataset = (GDALDataset *)GDALOpen(pszFilename, GA_ReadOnly);
    if (poDataset == NULL)
    {
        printf_s("read failed!\n");
    }
    else
    {
        printf_s("read successful!\n");
    }
    double adfGeoTransform[6];
    if (poDataset->GetGeoTransform(adfGeoTransform) == CE_Failure)//讀取坐標(biāo)信息
    {
        printf("獲取參數(shù)失敗");
    }
    const char *prj = poDataset->GetProjectionRef();//讀取投影信息


    int iWidth = poDataset->GetRasterXSize();//圖像寬度
    int iHeight = poDataset->GetRasterYSize();//圖像高度
    int iBandCount = poDataset->GetRasterCount();//波段數(shù)
    int iImageSize = iWidth * iHeight;//圖像像元數(shù)
    
    double *pBuff1 = new double[iImageSize*iBandCount];//開辟空間存儲(chǔ)原始圖像
    
    poDataset->RasterIO(GF_Read, 0, 0, iWidth, iHeight, pBuff1,
        iWidth, iHeight, GDT_Float64, iBandCount, 0, 0, 0, 0);//讀取原始圖像
    
    MatrixXd staMat = Map<MatrixXd>(pBuff1, iImageSize, iBandCount);//將圖像讀入eigen矩陣


    MatrixXd X(iImageSize, iBandCount), C(iBandCount, iBandCount);//按波段存儲(chǔ)至X矩陣,構(gòu)建協(xié)方差矩陣C
    MatrixXd vec, val;//構(gòu)建特征向量湿弦、特征值矩陣vec瓤漏、val
    X = MatrixXd(staMat);
    
    //零均值化
    featurenormalize(X);
    
    //計(jì)算協(xié)方差
    computeCov(X, C);
    
    //計(jì)算特征值和特征向量
    computeEig(C, vec, val);
    
    //計(jì)算損失率,確定降低維數(shù)
    int dim = computeDim(val);
    //計(jì)算結(jié)果
    MatrixXd res = X * vec.rightCols(dim);

    //將主成分分量存儲(chǔ)至pBuff2
    double *pBuff2 = new double[iImageSize*iBandCount];


    for (int i = 0; i < dim; ++i)
    {
        for (int j = 0; j < iImageSize; ++j)
        {
            pBuff2[i*iImageSize + j] = res(j, i);
    
        }
    }


    //各個(gè)主成分寫入圖像(包含坐標(biāo)及投影信息)
    for (int i = 0; i < iBandCount; i++)
    {
        char x[]=" ";
        strcpy(x, outPath);
        char dstPath[10] = {};
        sprintf(dstPath, "%d.tif", i + 1);
        strcat(x, dstPath);
        writePcaImg(x, iWidth, iHeight, pBuff2, adfGeoTransform, prj, 7, iImageSize, i + 1);
        cout << "pca " << i + 1 << " complete" << endl;
            
    }

    cout << "pca complete!" << endl;
    cin.get();
    return 0;


}
最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
  • 序言:七十年代末颊埃,一起剝皮案震驚了整個(gè)濱河市蔬充,隨后出現(xiàn)的幾起案子,更是在濱河造成了極大的恐慌班利,老刑警劉巖饥漫,帶你破解...
    沈念sama閱讀 212,383評(píng)論 6 493
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場離奇詭異罗标,居然都是意外死亡庸队,警方通過查閱死者的電腦和手機(jī),發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 90,522評(píng)論 3 385
  • 文/潘曉璐 我一進(jìn)店門闯割,熙熙樓的掌柜王于貴愁眉苦臉地迎上來彻消,“玉大人,你說我怎么就攤上這事纽谒≈づ颍” “怎么了如输?”我有些...
    開封第一講書人閱讀 157,852評(píng)論 0 348
  • 文/不壞的土叔 我叫張陵鼓黔,是天一觀的道長。 經(jīng)常有香客問我不见,道長澳化,這世上最難降的妖魔是什么? 我笑而不...
    開封第一講書人閱讀 56,621評(píng)論 1 284
  • 正文 為了忘掉前任稳吮,我火速辦了婚禮缎谷,結(jié)果婚禮上,老公的妹妹穿的比我還像新娘灶似。我一直安慰自己列林,他們只是感情好,可當(dāng)我...
    茶點(diǎn)故事閱讀 65,741評(píng)論 6 386
  • 文/花漫 我一把揭開白布酪惭。 她就那樣靜靜地躺著希痴,像睡著了一般。 火紅的嫁衣襯著肌膚如雪春感。 梳的紋絲不亂的頭發(fā)上砌创,一...
    開封第一講書人閱讀 49,929評(píng)論 1 290
  • 那天虏缸,我揣著相機(jī)與錄音,去河邊找鬼嫩实。 笑死刽辙,一個(gè)胖子當(dāng)著我的面吹牛,可吹牛的內(nèi)容都是我干的甲献。 我是一名探鬼主播宰缤,決...
    沈念sama閱讀 39,076評(píng)論 3 410
  • 文/蒼蘭香墨 我猛地睜開眼,長吁一口氣:“原來是場噩夢(mèng)啊……” “哼晃洒!你這毒婦竟也來了撵溃?” 一聲冷哼從身側(cè)響起,我...
    開封第一講書人閱讀 37,803評(píng)論 0 268
  • 序言:老撾萬榮一對(duì)情侶失蹤锥累,失蹤者是張志新(化名)和其女友劉穎缘挑,沒想到半個(gè)月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體桶略,經(jīng)...
    沈念sama閱讀 44,265評(píng)論 1 303
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡语淘,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 36,582評(píng)論 2 327
  • 正文 我和宋清朗相戀三年,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了际歼。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片惶翻。...
    茶點(diǎn)故事閱讀 38,716評(píng)論 1 341
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡,死狀恐怖鹅心,靈堂內(nèi)的尸體忽然破棺而出吕粗,到底是詐尸還是另有隱情,我是刑警寧澤旭愧,帶...
    沈念sama閱讀 34,395評(píng)論 4 333
  • 正文 年R本政府宣布颅筋,位于F島的核電站,受9級(jí)特大地震影響输枯,放射性物質(zhì)發(fā)生泄漏议泵。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 40,039評(píng)論 3 316
  • 文/蒙蒙 一桃熄、第九天 我趴在偏房一處隱蔽的房頂上張望先口。 院中可真熱鬧,春花似錦瞳收、人聲如沸碉京。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,798評(píng)論 0 21
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽谐宙。三九已至,卻和暖如春血崭,著一層夾襖步出監(jiān)牢的瞬間卧惜,已是汗流浹背厘灼。 一陣腳步聲響...
    開封第一講書人閱讀 32,027評(píng)論 1 266
  • 我被黑心中介騙來泰國打工, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留咽瓷,地道東北人设凹。 一個(gè)月前我還...
    沈念sama閱讀 46,488評(píng)論 2 361
  • 正文 我出身青樓,卻偏偏與公主長得像茅姜,于是被迫代替她去往敵國和親闪朱。 傳聞我的和親對(duì)象是個(gè)殘疾皇子,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 43,612評(píng)論 2 350

推薦閱讀更多精彩內(nèi)容