遥感图像mnf算法

     

 

 

 

 1. 使用右和上方做差比较方法，求取噪声协方差。右边界和上边界需要使用镜面反弹

2. Noise协方差矩阵，0.5的要求来自于Green，1988文章

3. eig(CovI, CovN); 这里的CovN注意方向

 4. 求特征值以后特征值大小排列正好反了

5. 去均值

6. 与特征向量相乘

求协方差：

template <typename T> bool BILMnfCovAndMean(HANDLE hReadFile, ulong Width, ulong Height, int Bands, 

    double* noiseMean, double* noiseCov, double* mean, double* cov)//传入的矩阵值的初始化必须为0

{

    /*

        MNF算法, cov(x,y) = 1/(n-1)*( 求和xi*yi - n*x*y) 协方差通过分解得到，速度比原始的公式要块,并且方便编程

        使用右和上方做差比较方法，求取噪声协方差。右边界和上边界需要使用镜面反弹

        噪声协方差 = 0.5*cov(噪声)= 0.5 * (1/(n-1)) * (求和Ni*Nj - n*N*N)    Noise协方差矩阵，

        date: 20170708

        author: 南京small山炮

    */

    if (hReadFile == NULL || noiseMean == NULL || noiseCov == NULL || mean == NULL || cov == NULL)

    {

        printf("Fuck, NULL Value Has Been Found!\n");

        return false;

    }

    if (Width <= 1 || Height <= 1 || Bands <= 1)

    {

        printf("高或宽或波段小于等于1，玩个蛋蛋!\n");

        return false;

    }

    ulong readBytes = 0;

    int k = 0;  T tTmp = 0;

    ulong BWF = Bands * Width * sizeof(double);

    ulong BWT = Bands * Width  * sizeof(T);

    ulong proSize = BWT;

    if (BWT < BWF) proSize = BWF;

    double* pNMean = noiseMean; double* pMean = mean;

    double* pNFBuf = NULL; T* pTPreBuf = NULL; T* pTImg = NULL;

    double* pCov = cov; double* pNCov = noiseCov;

    if (proSize <= DistributeMemorySize())

    {

        if (proSize <= SIZE16M)     k = SIZE16M  / proSize;

        else if(proSize <= SIZE32M) k = SIZE32M  / proSize;

        else if(proSize <= SIZE64M) k = SIZE64M  / proSize;

        else                        k = SIZE128M / proSize;

        double dNVal = 0;  ulong KBWT = k * BWT;

        pTPreBuf = (T*)malloc(BWT); pTImg = (T*)malloc(KBWT); pNFBuf = (double*)malloc(Width * sizeof(double));

        //第一次，求镜像的数据

        ReadFile(hReadFile,(LPVOID)pTPreBuf,BWT,&readBytes,NULL); //存放上一次的数据

        ReadFile(hReadFile,(LPVOID)pTImg,   BWT,&readBytes,NULL);

        T* pTPre  = pTPreBuf; T* pTCur = pTImg; 

        T* pJBuf = NULL; T* pJCur = NULL; T* pJPre = NULL;

        for (int i = 0; i < Bands; i++)

        {

            pJBuf = pTPre; pJCur = pTCur; pJPre = pTPre; pCov += i;   pNCov += i; 

            for (ulong x = 0; x < Width - 1; x++)

            {

                tTmp = *pJCur++;

                dNVal = tTmp - (*pJCur + *pJPre++)/2.0;

                pNFBuf[x] = dNVal;  

                pNMean[i] += dNVal; 

                *pNCov += dNVal * dNVal; 

                *pCov +