eigensystem.cpp

基于小波的SAR斑点处理

//////////////////////////////////////////////////////////////////////
// 求实对称矩阵的特征值和特征向量									 //
// 使用HouseHolder约化和因子分解方法								 //
// 开发: 余家忠														//
// 时间: 2000年6月27日												//
/////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "math.h"
#include "EigenSystem.h"

// 辅助函数,根据参数b的值来确定返回值a的符号
double Sign(double a,double b)
{
	return ( b>=0.0 ? fabs(a) : -fabs(a) );
}

// 辅助函数,求参数a和b的平方和的平方根
double Pythag(double a,double b)
{
	return sqrt( a*a + b*b );
}

// 辅助函数,利用HouseHolder变换将实对称矩阵约化为三对角矩阵.
// symMatrix为实对称矩阵,大小为n*n.输出时,被产生变换的正交矩阵取代,
// d[n]返回三对角矩阵的对角元素,
// e[n]返回非对角元素,且有e[0]=0.
void Simplify(double * symMatrix,int n, double * d, double * e)
{
	int i,j,k,m;
	double scale,hh,h,g,f;

	ASSERT(symMatrix!=NULL && d!=NULL && e!=NULL);

	for( i=n-1; i>0; i-- )
	{
		m = i-1;
		h = scale = 0.0;
		if( m > 0 )
		{
			for( k=0; k<=m; k++)
				scale += fabs(symMatrix[i*n+k]);
			if( scale == 0.0 )	// 跳过变换
				e[i] = symMatrix[i*n+m];
			else
			{
				for( k=0; k<=m; k++ )
				{
					// 将标定的a用于变换
					symMatrix[i*n+k] /= scale;
					// 在h中形成delta
					h += symMatrix[i*n+k]*symMatrix[i*n+k];
				}
				f = symMatrix[i*n+m];
				g = ( f>=0.0 ? -sqrt(h) : sqrt(h) );
				e[i] = scale*g;
				h -= f*g;
				symMatrix[i*n+m] = f-g;	// 将u存入symMatrix的第i行
				f = 0.0;
				for(j=0; j<=m; j++)
				{
					// 将u/H存入symMatrix的第i列
					symMatrix[j*n+i] = symMatrix[i*n+j]/h;
					g = 0.0;	// 在g中A*u的一个元素
					for(k=0; k<=j; k++)
						g += symMatrix[j*n+k]*symMatrix[i*n+k];
					for(k=j+1; k<=m; k++)
						g += symMatrix[k*n+j]*symMatrix[i*n+k];
					e[j] = g/h;	// 在e的暂时不用的元素中形成p的元素

					f += e[j]*symMatrix[i*n+j];
				}

				hh = f/(h+h);	// 形成K
				for(j=0; j<=m; j++)	// 形成q,并存入e中p的位置上
				{
					f = symMatrix[i*n+j];
					e[j] = g = e[j]-hh*f;	// 注意:有e[m]=e[i-1]
					for(k=0; k<=j; k++)		// 约化symMatrix
						symMatrix[j*n+k] -= (f*e[k]+g*symMatrix[i*n+k]);
				}
			}
		}
		else
			e[i] = symMatrix[i*n+m];
		d[i] = h;
	}

	d[0] = 0.0;
	e[0] = 0.0;
	for(i=0; i<n; i++)		// 开始矩阵变换的积累
	{
		m = i-1;
		if(d[i])			// 当i=0时跳过这一块
		{
			for(j=0; j<=m; j++)
			{
				g = 0.0;
				
				// 利用symMatrix中存储的u和u/H来形成P*Q
				for(k=0; k<=m; k++)
					g += symMatrix[i*n+k]*symMatrix[k*n+j];
				for(k=0; k<=m; k++) 
					symMatrix[k*n+j] -= g*symMatrix[k*n+i];
			}
		}

		d[i] = symMatrix[i*n+i];

		// 为下一次迭代将symMatrix重新置成单位矩阵
		symMatrix[i*n+i] = 1.0;
		for(j=0; j<=m; j++)	symMatrix[j*n+i] = symMatrix[i*n+j] = 0.0;
	}
}

// 利用因子分解的方法求三对角矩阵的特征值和特征向量.
// 采用隐含位移的QL算法(正交、下三角)。输入时，d[1...n]
// 为三对角矩阵的对角元素，输出时，它返回特征值。e[1...n]
// 为输入三对角矩阵的次对角元素值，e[1]为任意值，输出时，
// e中的值已背破坏。如果需要求一个三对角矩阵的特征向量，则
// 矩阵tMatrix[1...n][1...n]中输入单位矩阵;如果需要求一个由tred2
// 已约化的矩阵之特征向量,则tMatrix由tred2输出的矩阵作为输入. 不管
// 在哪种情况下,tMatrix的第k列返回与d[k]相对应的归一化特征向量.
//
void TQLI(double * tMatrix,int n,double * d,double * e)
{
	int i,j,k,m,iter;
	double s,r,p,g,f,dd,c,b;

	ASSERT( tMatrix!=NULL && d!=NULL && e!=NULL );

	// 为方便,重编e中的元素
	for(i=1; i<n; i++)
		e[i-1] = e[i];
	e[n-1] = 0.0;

	for( j=0; j<n; j++)
	{
		iter = 0;
		do{
			// 寻找一个单一的小的次对角元素用来分裂矩阵
			for(m=j; m<n-1; m++)
			{
				dd = fabs(d[m])+fabs(d[m+1]);
				if((float)(fabs(e[m])+dd) == dd) break;
			}
			if( m != j )
			{
				if(iter++ == 30)
				{
					AfxMessageBox("Too many iterations!\nThe result is false!");
					return;
				}
				g = (d[j+1]-d[j])/(2.0*e[j]);	// 形成位移
				r = Pythag(g,1.0);
				g = d[m]-d[j]+e[j]/(g+Sign(r,g));	// 这是dm-ds
				s = c = 1.0;
				p = 0.0;
				for(i=m-1; i>=j; i--)	// 如同原来QL方法一样的平面转动,
				{						// 随后吉文斯旋转,用以恢复三对角形式.
					f = s*e[i];
					b = c*e[i];
					e[i+1] = ( r = Pythag(f,g));
					if(r == 0.0)		// 从下溢处恢复
					{
						d[i+1] -= p;
						e[m] = 0.0;
						break;
					}
					s = f/r;
					c = g/r;
					g = d[i+1]-p;
					r = (d[i]-g)*s+2.0*c*b;
					d[i+1] = g+(p = s*r);
					g = c*r-b;
					for(k=0; k<n; k++)	// 形成特征向量
					{
						f = tMatrix[k*n+i+1];
						tMatrix[k*n+i+1] = s*tMatrix[k*n+i]+c*f;
						tMatrix[k*n+i] = c*tMatrix[k*n+i]-s*f;
					}
				}
				if(r == 0.0 && i >= j) continue;
				d[j] -= p;
				e[j] = g;
				e[m] = 0.0;
			}
		}while(m != j);
	}
}