#pragma once
#include "DibImage.h"
#include "Exception.h"
//矩阵类-声明
class Matrix
{
private:
//矩阵宽度
long width;
//矩阵高度
long height;
protected:
//矩阵的首地址
double * p;
public:
static Matrix One(long n);
void info();//显示信息
long getHeight();
long getWidth();
bool isVector();//如果是false,那就是一个数
double Deter();//求行列式determinant
bool isPtv();//是否正定
Matrix Eigen();//特征值
Matrix Adjoint();//伴随矩阵
Matrix T();//转置
Matrix Inv();//逆矩阵
void Unity();//归一化
void Store(const char * filename);//把矩阵存储到文件里
Matrix BiCubic(long n);//构造如下类似矩阵
/*{{2,0.5,0,0,0,0},
{0.5,2,0.5,0,0,0},
{0,0.5,2,0.5,0,0},
{0,0,0.5,2,0.5,0},
{0,0,0,0.5,2,0.5,},
{0,0,0,0,0.5,2,}
};*/
Matrix subMatrix(long offset) ;//throw (Exception &);
Matrix subMatrix(long x,long y,long WidthOffset,long HeightOffset);//取整数位置子区
//Matrix subMatrix(double x,double y,long WidthOffset,long HeightOffset);//取小数位置子区
void Test();//测试函数
/***********重载部分-Overloaded Part*****************/
Matrix operator+(Matrix &);
Matrix operator-(Matrix &);
Matrix operator*(Matrix &);
friend Matrix operator*(double alpha,Matrix &);//实数与矩阵相乘
Matrix operator*(double alpha);//矩阵与实数相乘
Matrix operator/(Matrix &);//实际是实数相除或矩阵和实数相除
Matrix operator/(double sub);
Matrix operator+=(Matrix &);
Matrix operator-=(Matrix &);
Matrix operator*=(Matrix &);//矩阵与实数相乘
Matrix operator*=(double alpha);//矩阵与实数相乘
Matrix &operator=(Matrix &);//赋值
double * operator[](long heightPos);//用于实现用[][]操作矩阵元素
friend Matrix sqrt(Matrix m);//开方
friend double abs(Matrix &);//取绝对值(泛数)
friend double sum(Matrix &);//求和
friend Matrix multiply(Matrix &m1,Matrix & m2);//按元素相乘
friend double operator+(double dbl,Matrix &);
friend double operator+(Matrix &,double dbl);
friend double operator-(double dbl,Matrix &);
friend double operator-(Matrix &,double dbl);
friend ostream & operator<<(ostream &,Matrix &);//输出
Matrix(void);//默认构造函数
Matrix(long n);//单位矩阵
Matrix(DibImage & dibImage);//类型转换(把图像文件转换成矩阵,用于图像计算处理)
Matrix(double * arrAddress,long arrWidth);//构造一维矩阵(一行,arrWidth列)
Matrix(double * arrAddress,long arrHeight,long arrWidth);//构造二维矩阵
Matrix(long Height,long Width);//构造空矩阵
Matrix(const Matrix &);//复制构造函数
virtual ~Matrix(void);//默认析构函数
/***********重载部分-Overloaded Part*****************/
};
#include "StdAfx.h"
#include "./matrix.h"
//矩阵类-实现
#define SIGN(a,b) ((b) > 0 ? fabs(a) : -fabs(a))
Matrix::Matrix(void):width(1),height(1)
{
p=new double[1];
//width=1;
//height=1;
//cout<<"Default Constructor Invoked."<<endl;
}
Matrix::Matrix(long n)
{
height=width=n;
p=new double[n*n];
for(long i=0;i<n;i++){
for(long j=0;j<n;j++){
if(i==j) *(p+n*i+j)=1;
else *(p+n*i+j)=0;
}
}
}
Matrix::Matrix(DibImage & dibImage)
{
width=dibImage.getWidth();height=dibImage.getHeight();
p=new double[width*height];
if(p==NULL) assert(0);//内存分配失败
long lLineBytes = WIDTHBYTES(width * 8);
LPSTR Src = ::FindDIBBits((LPSTR)::GlobalLock((HGLOBAL) dibImage.GetHDIB()));
for(long i=0;i<height;i++){
for(long j=0;j<width;j++){
*(p+width*i+j)=double(*((unsigned char *)Src + lLineBytes * (height-1-i) + j));
// *(p+width*i+j)=double(*((unsigned char *)Src + lLineBytes * (i) + j));
}
}
::GlobalUnlock((HGLOBAL) dibImage.GetHDIB());
}
Matrix::Matrix(double * arrAddress,long arrWidth)
{
long arrHeight=1;
p=new double[arrWidth*arrHeight];
for(long i=0;i<arrHeight;i++){
for(long j=0;j<arrWidth;j++){
*(p+arrWidth*i+j)=*(arrAddress+arrWidth*i+j);
}
}
width=arrWidth;
height=arrHeight;
}
Matrix::Matrix(double * arrAddress,long arrHeight,long arrWidth)
{
p=new double[arrWidth*arrHeight];
for(long i=0;i<arrHeight;i++){
for(long j=0;j<arrWidth;j++){
*(p+arrWidth*i+j)=*(arrAddress+arrWidth*i+j);
}
}
width=arrWidth;
height=arrHeight;
}
Matrix::Matrix(long Height,long Width)
{
p=new double[Height*Width];
width=Width;
height=Height;
}
Matrix::Matrix(const Matrix & m)//copy constructor
{
height=m.height;
width=m.width;
p=new double[height*width];
for(long i=0;i<height;i++){
for(long j=0;j<width;j++){
*(p+width*i+j)=*(m.p+width*i+j);
}
}
}
long Matrix::getWidth(){
return width;
}
long Matrix::getHeight(){
return height;
}
Matrix Matrix::One(long n)
{
Matrix m(1,n);
for(long i=0;i<n;i++)
m[0][i]=i;
return m;
}
bool Matrix::isVector()
{
return !(width==1 && height==1);
}
Matrix Matrix::subMatrix(long offset) //throw (Exception &)
{
if(!(height==width && offset<=width && offset>=0))
throw Exception("Error at Matrix Matrix::subMatrix(long offset)./n NOT height==width && offset<=width && offset>=0");
double * t=new double[offset*offset];
for(long i=0;i<offset;i++){
for(long j=0;j<offset;j++){
*(t+offset*i+j)=*(p+width*i+j);
}
}
Matrix m(t,offset,offset);
delete [] t;
return m;
}
Matrix Matrix::subMatrix(long x,long y,long WidthOffset,long HeightOffset)
{
//cout<<"x="<<x<<" y="<<y<<" w="<<WidthOffset<<" h="<<HeightOffset<<" in subMatrix 1/n";
bool log=true;
if(x<0 || y<0){
cout<<"x<0 or y<0 in Matrix::subMatrix()!/n";
log=false;
}else if(WidthOffset<0 || HeightOffset<0){
cout<<"WidthOffset<0 || HeightOffset<0 in Matrix::subMatrix()!/n";
log=false;
}else if(WidthOffset*2+1>width || HeightOffset*2+1>height){
cout<<">Subarea is bigger than Matrix in Matrix::subMatrix()!/n";
log=false;
}else if((x-WidthOffset)<0 || (x+WidthOffset)>width){
cout<<"Width error in Matrix::subMatrix()!/n";
log=false;
}else if((y-HeightOffset)<0 || (y+HeightOffset)>height){
cout<<"Height error in Matrix::subMatrix()!/n";
log=false;
}
assert(log);
Matrix m(HeightOffset*2+1,WidthOffset*2+1);
for(long i=0;i<HeightOffset*2+1;i++){
for(long j=0;j<WidthOffset*2+1;j++){
m[i][j]=(*this)[i+y-HeightOffset][j+x-WidthOffset];
}
}
return m;
}
double Matrix::Deter()//矩阵的行列式
{
assert(width==height);
double result=1;
double k;
Matrix m=*this;
for(long i=0;i<height-1;i++){//Gauss消去法,变换成上三角矩阵
for(long j=i+1;j<height;j++){
k=m[j][i]/m[i][i];
m[j][i]=0;
for(long n=i+1;n<width;n++){
m[j][n]=m[j][n]-k*m[i][n];
}
}
}
for(long i=0;i<height;i++){
for(long j=0;j<width;j++){
if(i==j) result*=m[i][j];
}
}
return result;
}
bool Matrix::isPtv()
{
assert(width==height);//是方阵才可以计算
bool result=true;
Matrix m;
for(long i=1;i<=height;i++){
m=this->subMatrix(i);
if(m.Deter()<=0){
result=false;
break;
}
}
return result;
}
Matrix Matrix::T()
{
double * t=new double[width*height];
for(long i=0;i<height;i++){
for(long j=0;j<width;j++){
*(t+height*j+i)=*(p+width*i+j);
}
}
Matrix m(t,width,height);
delete [] t;
return m;
}
Matrix Matrix::Eigen()
{
Matrix matrix=*this;
assert(matrix.getHeight()==matrix.getWidth());
double **a,*wr,*wi;
long nn,m,l,k,j,its,i,mmin;
double z,y,x,w,v,u,t,s,r,q,p,anorm;
long n=matrix.getWidth();
a=new double*[n+1];
wr=new double[n+1];
wi=new double[n+1];
for(i=0;i<n+1;i++)
a[i]=new double[n+1];
for(i=0;i<n;i++){
for(j=0;j<n;j++){
a[i+1][j+1]=matrix[i][j];
}
}
anorm=fabs(a[1][1]);
for (i=2;i<=n;i++)
for (j=(i-1);j<=n;j++)
anorm += fabs(a[i][j]);
nn=n;
t=0.0;
while (nn >= 1) {
its=0;
do {
for (l=nn;l>=2;l--) {
s=fabs(a[l-1][l-1])+fabs(a[l][l]);
if (s == 0.0) s=anorm;
if (fabs(a[l][l-1]) + s == s) break;
}
x=a[nn][nn];
if (l == nn) {
wr[nn]=x+t;
wi[nn--]=0.0;
} else {
y=a[nn-1][nn-1];
w=a[nn][nn-1]*a[nn-1][nn];
if (l == (nn-1)) {
p=0.5*(y-x);
q=p*p+w;
z=sqrt(fabs(q));
x += t;
if (q >= 0.0) {
z=p+SIGN(z,p);
wr[nn-1]=wr[nn]=x+z;
if (z) wr[nn]=x-w/z;
wi[nn-1]=wi[nn]=0.0;
} else {
wr[nn-1]=wr[nn]=x+p;
wi[nn-1]= -(wi[nn]=z);
}
nn -= 2;
} else {
if (its == 30) printf("Too many iterations in HQR");
if (its == 10 || its == 20) {
t += x;
for (i=1;i<=nn;i++) a[i][i] -= x;
s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]);
y=x=0.75*s;
w = -0.4375*s*s;
}
++its;
for (m=(nn-2);m>=l;m--) {
z=a[m][m];
r=x-z;
s=y-z;
p=(r*s-w)/a[m+1][m]+a[m][m+1];
q=a[m+1][m+1]-z-r-s;
r=a[m+2][m+1];
s=fabs(p)+fabs(q)+fabs(r);
p /= s;
q /= s;
r /= s;
if (m == l) break;
u=fabs(a[m][m-1])*(fabs(q)+fabs(r));
v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));
if (u+v == v) break;
}
for (i=m+2;i<=nn;i++) {
a[i][i-2]=0.0;
if (i != (m+2)) a[i][i-3]=0.0;
}
for (k=m;k<=nn-1;k++) {
if (k != m) {
p=a[k][k-1];
q=a[k+1][k-1];
r=0.0;
if (k != (nn-1)) r=a[k+2][k-1];
if (x=fabs(p)+fabs(q)+fabs(r)) {
p /= x;
q /= x;
r /= x;
}
}
if (s=SIGN(sqrt(p*p+q*q+r*r),p)) {
if (k == m) {
if (l != m)
a[k][k-1] = -a[k][k-1];
} else
a[k][k-1] = -s*x;
p += s;
x=p/s;
y=q/s;
z=r/s;
q /= p;
r /= p;
for (j=k;j<=nn;j++) {
p=a[k][j]+q*a[k+1][j];
if (k != (nn-1)) {
p += r*a[k+2][j];
a[k+2][j] -= p*z;
}
a[k+1][j] -= p*y;
a[k][j] -= p*x;
}
mmin = nn<k+3 ? nn : k+3;
for (i=l;i<=mmin;i++) {
p=x*a[i][k]+y*a[i][k+1];
if (k != (nn-1)) {
p += z*a[i][k+2];
a[i][k+2] -= p*r;
}
a[i][k+1] -= p*q;
a[i][k] -= p;
}
}
}
}
}
} while (l < nn-1);
}
Matrix eigenvalue(1,n);
for(i=0;i<n;i++){
eigenvalue[0][i]=wr[i+1];
}
return eigenvalue;
}
Matrix Matrix::Adjoint()
{
return this->Inv()*this->Deter();
}
void Matrix::Unity()
{
double total;
total=sqrt(sum(multiply(*this,*this)));
for(long i=0;i<this->height;i++){
for(long j=0;j<this->width;j++){
(*this)[i][j]/=total;
}
}
return;
}
void Matrix::Store(const char * filename){
ofstream fout(filename);
int w=5;
for(long i=0;i<height;i++){
for(long j=0;j<width;j++){
fout.width(w);
fout<<i;
fout.width(w);
fout<<j;
fout.width(15);
fout<<(*this)[i][j]<<endl;/*
cout.width(width);
cout<<i;
cout.width(width);
cout<<j;
cout.width(width);
cout<<(*this)[i][j]<<endl;*/
}
}
fout.close();
}
Matrix Matrix::Inv()
{
assert(width==height && width>0);
Matrix m=*this;
double * pDbSrc=m.p;
long nLen=m.width;
int *is,*js,i,j,k;
double d,p;
is = new int[nLen];
js = new int[nLen];
for(k=0;k<nLen;k++)
{
d=0.0;
for(i=k;i<nLen;i++)
for(j=k;j<nLen;j++)
{
p=fabs(pDbSrc[i*nLen + j]);
if(p>d)
{
d = p;
is[k] = i;
js[k] = j;
}
}
if(d+1.0==1.0)
{
delete is;
delete js;
assert(0);
//return FALSE;
}
if(is[k] != k)
for(j=0;j<nLen;j++)
{
p = pDbSrc[k*nLen + j];
pDbSrc[k*nLen + j] = pDbSrc[(is[k]*nLen) + j];
pDbSrc[(is[k])*nLen + j] = p;
}
if(js[k] != k)
for(i=0; i<nLen; i++)
{
p = pDbSrc[i*nLen + k];
pDbSrc[i*nLen + k] = pDbSrc[i*nLen + (js[k])];
pDbSrc[i*nLen + (js[k])] = p;
}
pDbSrc[k*nLen + k]=1.0/pDbSrc[k*nLen + k];
for(j=0; j<nLen; j++)
if(j != k)
{
pDbSrc[k*nLen + j]*=pDbSrc[k*nLen + k];
}
for(i=0; i<nLen; i++)
if(i != k)
for(j=0; j<nLen; j++)
if(j!=k)
{
pDbSrc[i*nLen + j] -= pDbSrc[i*nLen + k]*pDbSrc[k*nLen + j];
}
for(i=0; i<nLen; i++)
if(i != k)
{
pDbSrc[i*nLen + k] *= -pDbSrc[k*nLen + k];
}
}
for(k=nLen-1; k>=0; k--)
{
if(js[k] != k)
for(j=0; j<nLen; j++)
{
p = pDbSrc[k*nLen + j];
pDbSrc[k*nLen + j] = pDbSrc[(js[k])*nLen + j];
pDbSrc[(js[k])*nLen + j] = p;
}
if(is[k] != k)
for(i=0; i<nLen; i++)
{
p = pDbSrc[i*nLen + k];
pDbSrc[i*nLen + k] = pDbSrc[i*nLen +(is[k])];
pDbSrc[i*nLen + (is[k])] = p;
}
}
delete is;
delete js;
return m;
}
Matrix Matrix::BiCubic(long n)
{
assert(n>1);
delete [] p;
width=n;
height=n;
p=new double[n*n];
//首行
*(p)=0.5; *(p+1)=2;
for(long k=2;k<n;k++) *(p+k)=0;
//中间行
for(long i=1;i<n-1;i++){//行
for(long j=0;j<n;j++){//列
if(j==(i-1)){
*(p+n*i+j)=0.5;
j++;
*(p+n*i+j)=2;
j++;
*(p+n*i+j)=0.5;
j++;
}
*(p+n*i+j)=0;
}
}
//末行
for(long l=0;l<n-2;l++) *(p+n*(n-1)+l)=0;
*(p+n*n-1)=2;
*(p+n*n-2)=0.5;
return *this;
}
Matrix Matrix::operator +(Matrix &m1)
{
assert(m1.height==height && m1.width==width);
long tmpHeight=m1.height;
long tmpWidth=m1.width;
double * t=new double[tmpWidth*tmpHeight];
for(long i=0;i<tmpHeight;i++){
for(long j=0;j<tmpWidth;j++){
*(t+tmpWidth*i+j)=*(m1.p+tmpWidth*i+j)+*(p+tmpWidth*i+j);
}
}
Matrix m(t,tmpHeight,tmpWidth);
delete [] t;
return m;
}
Matrix Matrix::operator -(Matrix &m1)
{
assert(m1.height==height && m1.width==width);
long tmpHeight=m1.height;
long tmpWidth=m1.width;
double * t=new double[tmpWidth*tmpHeight];
for(long i=0;i<tmpHeight;i++){
for(long j=0;j<tmpWidth;j++){
*(t+tmpWidth*i+j)=*(p+tmpWidth*i+j)-*(m1.p+tmpWidth*i+j);
}
}
Matrix m(t,tmpHeight,tmpWidth);
delete [] t;
return m;
}
Matrix Matrix::operator *(Matrix &m1)
{
if(!this->isVector() && m1.isVector()){//左为数,右为矩阵
Matrix m;
m=p[0]*m1;
return m;
}else if(this->isVector() && !m1.isVector()){//左为矩阵,右为数
Matrix m;
m=*this*m1[0][0];
return m;
}else if(!this->isVector() && m1.isVector()){//左右都为数
double * t=new double[1];
t[0]=p[0]*m1[0][0];
Matrix m(t,1,1);
delete [] t;
return m;
}else if(this->isVector() && m1.isVector() && width==m1.height){//左为矩阵,右为矩阵
double sum;
double * t=new double[height*m1.width];
for(long i=0;i<height;i++){
for(long j=0;j<m1.width;j++){
sum=0;
for(long k=0;k<width;k++){
sum+=(*(p+width*i+k))*(m1[k][j]);
}
*(t+m1.width*i+j)=sum;
}
}
Matrix m(t,height,m1.width);
delete [] t;
return m;
}else{
assert(0);//未知运算
return *this;
}
}
Matrix operator*(double alpha,Matrix & m1)
{
Matrix m=m1;
for(long i=0;i<m.height;i++){
for(long j=0;j<m.width;j++){
m[i][j]=alpha*m1[i][j];
}
}
return m;
}
Matrix Matrix::operator*(double alpha)
{
return alpha*(*this);
}
Matrix Matrix::operator+=(Matrix & m)
{
return *this+m;
}
Matrix Matrix::operator-=(Matrix & m)
{
return *this-m;
}
Matrix Matrix::operator *=(double alpha)
{
return *this*alpha;
}
Matrix Matrix::operator *=(Matrix & m1)
{
return *this*m1;
}
Matrix sqrt(Matrix m)
{
assert(m[0][0]>=0 && m.getHeight()==1 && m.getWidth()==1);
m[0][0]=sqrt(m[0][0]);
return m;
}
double abs(Matrix & m)
{
double sum=0;
for(long i=0;i<m.height;i++){
for(long j=0;j<m.width;j++){
sum+=m[i][j]*m[i][j];
}
}
return sqrt(sum);
}
double sum(Matrix & m)
{
double total=0;
for(long i=0;i<m.height;i++){
for(long j=0;j<m.width;j++){
total+=m[i][j];
}
}
return total;
}
Matrix multiply(Matrix & m1,Matrix & m2){
assert(m1.getHeight()==m2.getHeight() && m1.getWidth()==m2.getWidth());
Matrix m(m1.getHeight(),m1.getWidth());
for(long i=0;i<m1.getHeight();i++){
for(long j=0;j<m1.getWidth();j++){
m[i][j]=m1[i][j]*m2[i][j];
}
}
//cout<<m;
return m;
}
Matrix Matrix::operator /(Matrix &m1)
{
assert(m1.width==1 && m1.height==1);
assert(m1[0][0]!=0);
return *this/m1[0][0];
}
Matrix Matrix::operator /(double sub)
{
assert(sub!=0);
Matrix m=*this;
for(long i=0;i<height;i++){
for(long j=0;j<width;j++){
m[i][j]=*(p+width*i+j)/sub;
}
}
return m;
}
Matrix & Matrix::operator =(Matrix & m)
{
if(&m==this) return *this;
height=m.height;
width=m.width;
delete [] p;
p=new double[height*width];
for(long i=0;i<height;i++){
for(long j=0;j<width;j++){
*(p+width*i+j)=*(m.p+width*i+j);
}
}
return *this;
}
double * Matrix::operator [](long heightPos)
{
assert(heightPos>=0 && heightPos<=height);//报错
return p+width*heightPos;
}
double operator+(double dbl,Matrix & m)
{
assert(m.getHeight()==1 && m.getWidth()==1);
return dbl+m[0][0];
}
double operator+(Matrix & m,double dbl)
{
return dbl+m;
}
double operator-(double dbl,Matrix & m)
{
assert(m.getHeight()==1 && m.getWidth()==1);
return dbl-m[0][0];
}
double operator-(Matrix & m,double dbl)
{
return -(dbl-m);
}
ostream & operator<<(ostream & os,Matrix & m)
{
os<<"Matrix:height="<<m.height<<endl;
os<<"Matrix:width="<<m.width<<endl;
for(long i=0;i<m.height;i++){
for(long j=0;j<m.width;j++){
os.width(15);
os<<m[i][j];
}
os<<endl;
}
return os;
}
Matrix::~Matrix(void)
{
delete [] p;
//cout<<"Default Destructor Invoked."<<endl;
}
void Matrix::Test(){
double arr[4][4]={{1,1,1,1},{1,2,3,4},{1,3,6,10},{1,4,10,20}};
Matrix m1,m(&arr[0][0],4,4),m2(&arr[0][0],4,4);
m.isPtv();
/*
cout<<"arranger="<<m.Arg()<<endl;
cout<<m;*/
}
void Matrix::info(){
cout<<"height="<<height<<endl;
cout<<"width="<<width<<endl;
}