程式碼
試試看更長的code
#include <iostream>
#include <fstream>
#include <math.h>
#include <Dense>
#include <SVD>
using Eigen::MatrixXf;
using Eigen::JacobiSVD;
using namespace std;
#define POINT 4
int main()
{
fstream input;
int length;
int i,j,p,pixel;
int avg[POINT];
int sum;
float cov[POINT][POINT]={};
float x,y=0;
float KLT_eff,DCT_eff;
float r;
unsigned char X[POINT][512*512/POINT];
input.open("baboon.bmp", ios::in | ios::binary);
//get the file size.
input.seekg(0, ios::end);
length = input.tellg();
input.seekg(0, ios::beg);
char *data = new char[length];
input.read(data, length); //read the input file and write into the array.
//build the four vector X1, X2, X3, X4
for (i=0, pixel=1078; i<512*512/POINT; i++){
for (p=0; p<POINT;p++){
X[p][i]=data[pixel++];
}
}
//find each mean of X1, X2, X3, X4
for (p=0; p<POINT; p++){
for (i=0, sum=0; i<512*512/POINT; i++){
sum=X[p][i]+sum;
}
avg[p]=sum/(512*512/POINT);
}
//Covariance Matrix
for (i=0; i<POINT; i++){
for (j=0; j<POINT; j++){
for (p=0; p<512*512/POINT; p++){
cov[i][j]=(X[i][p]-avg[i])*(X[j][p]-avg[j])+cov[i][j];
}
cov[i][j]=cov[i][j]/((512*512/POINT)-1);
}
}
//print the matrix
cout<<"Covariance Matrix:"<<endl;
for (i=0;i<POINT;i++){
for (j=0;j<POINT;j++){
cout <<cov[i][j]<<" ";
}
cout <<endl;
}
//KLT basis vectors
MatrixXf m(POINT,POINT);
for (i=0; i<POINT; i++){
for (j=0; j<POINT; j++){
m(i,j)=cov[i][j];
}
}
//std::cout << "Here is the matrix m:" << std::endl << m << std::endl;
JacobiSVD<MatrixXf> svd(m, Eigen::ComputeThinU | Eigen::ComputeThinV);
std::cout <<endl<< "Its singular values are:" << std::endl << svd.singularValues() << std::endl;
//Tdct //inv_Tdct
MatrixXf Tdct(POINT,POINT);
MatrixXf inv_Tdct(POINT,POINT);
for (i=0; i<POINT; i++){
for (j=0; j<POINT; j++){
if (i==0)
Tdct(i,j)=1/(sqrt((float)POINT));
else
Tdct(i,j)=sqrt(2/(float)POINT)*cos((M_PI*(2*j+1)*i)/(2*POINT));
}
}
//cout<<endl<<"Tdct matrix:"<<endl<<Tdct<<endl;
//build the matrix to calculate inverse matrix
MatrixXf tmp(POINT,2*POINT);
for (i=0; i<POINT; i++){
for (j=0; j<POINT; j++){
tmp(i,j)=Tdct(i,j);
}
}
for (i=0; i<POINT; i++){
for (j=POINT; j<2*POINT; j++){
if ((j-POINT)==i) tmp(i,j)=1.0;
else tmp(i,j)=0;
}
}
//calculating the inverse matrix
for (i=0; i<POINT; i++){ //which row is standard to be one and elimate the other rows to be zero
float beone=tmp(i,i);
//let standard row (i,i) to be one, other col elements need to be resized.
for (j=0; j<2*POINT; j++){ //col
tmp(i,j)=tmp(i,j)/beone;
}
//elimating others row
for (int m=0; m<POINT; m++){ //row
r=tmp(m,i); //ref
for (int l=0; l<2*POINT; l++){ //col
if (m!=i){
tmp(m,l)=tmp(m,l)-tmp(i,l)*r;
}
}
}
}
//build inv_Tdct matrix
for (i=0;i<POINT;i++){
for (j=0;j<POINT;j++){
inv_Tdct(i,j)=tmp(i,j+POINT);
}
}
MatrixXf tmp2(POINT,POINT);
MatrixXf D(POINT,POINT);
//DCT matrix: D=Tdct*Cov*inv_Tdct
for (i=0; i<POINT; i++){
for (j=0; j<POINT; j++){
float buf=0;
for (int k=0; k<POINT; k++){
buf=buf+Tdct(i,k)*m(k,j); //tmp2 matrix = Tdct matrix* Covariance matrix; tmp2=Tdct*Czm
}
tmp2(i,j)=buf;
}
}
for (i=0; i<POINT; i++){
for (j=0; j<POINT; j++){
float buf=0;
for (int k=0; k<POINT; k++){
buf=buf+tmp2(i,k)*inv_Tdct(k,j); //D matrix = tmp2 matrix* inverse Tdct matrix; D=tmp2*inv_Tdct
}
D(i,j)=buf;
}
}
cout<<endl<<"D matrix:"<<endl<<D<<endl<<endl;
//De-correlation Efficiency for KLT
//Y
for (i=0,y=0;i<POINT;i++)
y=0;
//X
for (i=0,x=0;i<POINT;i++){
for (j=0;j<POINT;j++){
if (i!=j)
x=x+m(i,j);
}
}
KLT_eff=1-y/x;
cout <<"Decorrelation Efficiency: "<<endl<<"KLT: "<<KLT_eff*100<<"%"<<endl;
//De-correlation Efficiency for DCT
//Y
for (i=0,y=0;i<POINT;i++){
for (j=0;j<POINT;j++){
if (i!=j)
y=y+abs(D(i,j));
}
}
//X
for (i=0,x=0;i<POINT;i++){
for (j=0;j<POINT;j++){
if (i!=j)
x=x+abs(m(i,j));
}
}
DCT_eff=1-y/x;
cout<<"DCT: "<<DCT_eff*100<<"%"<<endl;
return 0;
}
參考:http://www.ewdna.com/2012/02/css-block.html
HTML encoder:http://formatmysourcecode.blogspot.tw/
會需要用到HTML encoder主要的原因就是一些符號在HTML裡是有特殊意義的(像是 > < &)
都需要轉換才能夠使用 否則code內的< > 就會被當成HTML的特殊意義使用
網頁設計(HTML學習):http://www.phd.com.tw/knowledge/html/typesetting/
將code轉成HTML:http://formatmysourcecode.blogspot.tw/
8/31 2015
