遥感图像mnf算法

最新推荐文章于 2023-01-08 15:47:45 发布

small山炮

最新推荐文章于 2023-01-08 15:47:45 发布

阅读量574

点赞数

文章标签：算法 c++

本文链接：https://blog.csdn.net/jikechong/article/details/126040131

版权

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 += *pJBuf * *pJBuf;

pMean[i] += *pJBuf++;

}

pNFBuf[Width - 1] = dNVal;

*pNCov += dNVal * dNVal;

pNMean[i] += dNVal;

*pCov += *pJBuf * *pJBuf;

pMean[i] += *pJBuf++;

pCov++; pNCov++; pJPre++; pJCur++;

for (int j = i + 1; j < Bands; j++)

{

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

{

tTmp = *pJCur++;

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

*pNCov += dNVal * pNFBuf[x];

*pCov += pTPre[x] * *pJBuf++;

}

*pCov += pTPre[Width - 1] * *pJBuf++;

pCov++;

*pNCov += dNVal * pNFBuf[Width - 1];

pNCov++; pJPre++; pJCur++;

}

pTCur += Width; pTPre += Width;

}

//SetFilePointer(hReadFile, -BWT, NULL, FILE_CURRENT);

SetFilePointer(hReadFile, BWT, NULL, FILE_BEGIN);

ulong kH = (Height - 1) / k;

for (ulong y = 0; y < kH; y++)

{

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

pTCur = pTImg;

for (ulong ik = 0; ik < k; ik++)

{

pCov = cov; pNCov = noiseCov; pTPre = pTPreBuf;

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

{

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

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

{

tTmp = *pJCur++;

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

*pJPre++ = tTmp;

pNFBuf[x] = dNVal;

pNMean[i] += dNVal;

*pNCov += dNVal * dNVal;

*pCov += *pJBuf * *pJBuf;

pMean[i] += *pJBuf++;

}

pNFBuf[Width - 1] = dNVal;

*pNCov += dNVal * dNVal;

pNMean[i] += dNVal;

pNCov++; *pJPre++ = *pJCur++;

*pCov += *pJBuf * *pJBuf;

pMean[i] += *pJBuf++;

pCov++;

for (int j = i + 1; j < Bands; j++)

{

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

{

tTmp = *pJCur++;

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

*pNCov += dNVal * pNFBuf[x];

*pCov += pTCur[x] * *pJBuf++;

}

*pCov += pTCur[Width - 1] * *pJBuf++;

pCov++;

*pNCov += dNVal * pNFBuf[Width - 1];

pNCov++; pJPre++; pJCur++;

}

pTCur += Width; pTPre += Width;

}

k = (Height - 1) - k * kH;

if (k != 0)

{

KBWT = k * BWT;

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

pTCur = pTImg;

for (ulong ik = 0; ik < k; ik++)

{

pCov = cov; pNCov = noiseCov; pTPre = pTPreBuf;

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

{

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

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

{

tTmp = *pJCur++;

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

*pJPre++ = tTmp;

pNFBuf[x] = dNVal;

pNMean[i] += dNVal;

*pNCov += dNVal * dNVal;

*pCov += *pJBuf * *pJBuf;

pMean[i] += *pJBuf++;

}

pNFBuf[Width - 1] = dNVal;

*pNCov += dNVal * dNVal;

pNMean[i] += dNVal;

pNCov++; *pJPre++ = *pJCur++;

*pCov += *pJBuf * *pJBuf;

pMean[i] += *pJBuf++;

pCov++;

for (int j = i + 1; j < Bands; j++)

{

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

{

tTmp = *pJCur++;

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

*pNCov += dNVal * pNFBuf[x];

*pCov += pTCur[x] * *pJBuf++;

}

*pCov += pTCur[Width - 1] * *pJBuf++;

pCov++;

*pNCov += dNVal * pNFBuf[Width - 1];

pNCov++; pJPre++; pJCur++;

}

pTCur += Width; pTPre += Width;

}

else //分块或者内存映射去做

{

}

if (pNFBuf != NULL){free(pNFBuf);pNFBuf = NULL;}

if (pTPreBuf != NULL){free(pTPreBuf);pTPreBuf = NULL;}

if (pTImg != NULL){free(pTImg);pTImg = NULL;}

ulong HW = Width * Height;

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

{

pMean[i] /= HW;

pNMean[i] /= HW;

}

ulong HW2 = HW - 1;

pCov = cov; double* p2Cov = NULL; double temp = 0;

pNCov = noiseCov; double* p2NCov = NULL; double tmpN = 0;

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

{

temp = HW * pMean[i];

tmpN = HW * pNMean[i];

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

{

*pCov = (*pCov - temp * pMean[j]) / HW2 ; pCov++;

*pNCov = ((*pNCov - tmpN * pNMean[j]) / HW2) * 0.5; pNCov++; //0.5的要求来自于Green，1988文章

}

if (i >= Bands - 1) break;

p2Cov = cov + i + 1; p2NCov = noiseCov + i + 1;

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

{

*pCov++ = *p2Cov; p2Cov += Bands;

*pNCov++ = *p2NCov; p2NCov += Bands;

}

return true;

}

求广义特征值和特征向量：

template <typename _T1, typename _T2>

bool RsEig(_T1* A, _T2* B, int n, double* V, double* D, const char* flag = "chol")

{

if (A == NULL || B == NULL || n <= 0)

{

return false;

}

unsigned int NND= n * n * sizeof(double);

double* L = new double[NND];

double* Tmp = new double[NND];

double* dt = new double[NND];

if (L == NULL || Tmp == NULL || dt == NULL)

{

goto ERR;

}

if(!RsChol<_T2>(B, n, n, L)) goto ERR;

if(!RsMatInv<double>(L, n)) goto ERR; //R的逆是L的逆转置

_T1* mA = A;

double* pL = L; double* pL2 = NULL;

double* pDT = dt;

double* pTmp = Tmp;

unsigned int JN = 0;

double a = 0;

//因为L是下三角, R就是上三角 . L*A*R 也是对称矩阵

for (int j = 0; j < n; j++)

{

JN = j * n; pL2 = L;

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

{

a = 0; mA = A + i; pL = L + JN;

for (int k = 0; k <= j; k++)

{

a += *pL++ * *mA; mA += n;

}

pDT[i] = a; a = 0;

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

{

a += pDT[k] * pL2[k];

}

*pTmp++ = a; pL2 += n;

}

memset(dt, 0 , NND);

//太慢了

jacobi(Tmp, dt, n);

double* pDT2 = NULL;

for (int j = 0; j < n; j++)

{

JN = j * n; D[j] = Tmp[JN + j]; pL = L + JN + j; pDT2 = dt + JN;

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

{

pL2 = pL; pDT = pDT2 + i;

a = 0;

for (int k = j; k < n; k++)

{

a += *pL2 * *pDT;

pDT += n; pL2 += n;

}

V[i * n + j] = -1.0 * a;

}

delete dt;

delete L;

delete Tmp;

return true;

ERR:

if (L != NULL)

{

delete L;

L = NULL;

}

if (Tmp != NULL)

{

delete Tmp;

Tmp = NULL;

}

if (dt == NULL)

{

delete dt;

dt = NULL;

}

return false;

}

small山炮

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
遥感图像mnf算法

mnf算法c++
复制链接

扫一扫