Skip to content
Snippets Groups Projects
Commit a2394430 authored by PIERRE Fabien's avatar PIERRE Fabien
Browse files

update readme

parent 1820c023
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 86 additions and 85 deletions
PB2_mex.cpp
+ 86
− 85
View file @ a2394430
...@@ -6,8 +6,8 @@ ...@@ -6,8 +6,8 @@
/************************* /*************************
* *
clc * clc
mex PB2_mex.cpp -output PB2 * mex PB2_mex.cpp -output PB2
* *
*/ */
...@@ -33,77 +33,77 @@ void projection_6(double * entree, double *sortie){ ...@@ -33,77 +33,77 @@ void projection_6(double * entree, double *sortie){
// Input:3x2 matrix A in double variables a11.. . a32 // Input:3x2 matrix A in double variables a11.. . a32
// Output : projection \Pi(A) in double variables a11.. . a32 // Output : projection \Pi(A) in double variables a11.. . a32
// Compute D= AˆT A // Compute D= AˆT A
double a11= entree[0];
double a12= entree[1];
double a21= entree[2];
double a22= entree[3];
double a31= entree[4];
double a32= entree[5];
double d11= a11 * a11 + a21 * a21 + a31 * a31 ;
double d12= a12 * a11 + a22 * a21 + a32 * a31 ;
double d22= a12 * a12 + a22 * a22 + a32 * a32 ;
// Compute Eigenvalues ( lmax , lmin )
//=squares of singular values ( smax , smin )
double trace = d11 + d22 ;
double det = d11 * d22 - d12 * d12 ;
double d = sqrt( max( 0.0 , 0.25 * trace * trace - det ) ) ;
double lmax = max( 0.0 , 0.5 * trace + d ) ;
double lmin = max( 0.0 , 0.5 * trace - d ) ;
double smax = sqrt ( lmax ) ;
double smin = sqrt ( lmin ) ;
// Check whether matrix lies outside K and need o be projed
if ( smax + smin > 1.0 ) {
// Compute orthonormalsy st emofEigenvectors , such t h a t
// ( v11 , v21 ) belongs to lmax and ( v12 , v22 ) belongs to lmin
double v11 , v12 , v21 , v22 ;
if ( d12 == 0.0 ) {
if ( d11 >= d22 ) { v11 = 1.0 ; v21 = 0.0 ; v12 = 0.0 ; v22 = 1.0 ; }
else
{ v11 = 0.0 ; v21 = 1.0 ; v12 = 1.0 ; v22 = 0.0 ; }
}
else {
v11 = lmax - d22 ; v21 = d12 ;
double l1 = hypot( v11 , v21 ) ;
v11 /= l1 ; v21 /= l1 ;
v12 = lmin - d22 ; v22 = d12 ;
double l2 = hypot( v12 , v22 ) ;
v12 /= l2 ; v22 /= l2 ;
}
// Project ( smax , smin ) toline ( 0 , 1 ) + tau * (1 , -1) , 0<=tau <=1.
double tau = 0.5 * ( smax - smin + 1.0 ) ;
double s1 = min ( 1.0 , tau ) ;
double s2 = 1.0 - s1 ;
// Compute \ Sigmaˆ+ * \ Sigma p
s1 /= smax ;
s2 = ( smin > 0.0 ) ? s2 / smin : 0.0 ;
// Compute double a11= entree[0];
double t11= s1 * v11 * v11 + s2 * v12 * v12 ; double a12= entree[1];
double t12= s1 * v11 * v21 + s2 * v12 * v22 ; double a21= entree[2];
double t21= s1 * v21 * v11 + s2 * v22 * v12 ; double a22= entree[3];
double t22= s1 * v21 * v21 + s2 * v22 * v22 ; double a31= entree[4];
double a32= entree[5];
// Result double d11= a11 * a11 + a21 * a21 + a31 * a31 ;
a11 = a11* t11 + a12 *t21 ; double d12= a12 * a11 + a22 * a21 + a32 * a31 ;
a12 = a11* t12 + a12 *t22 ; double d22= a12 * a12 + a22 * a22 + a32 * a32 ;
a21 = a21* t11 + a22 *t21 ;
a22 = a21* t12 + a22 *t22 ;
a31 = a31* t11 + a32 *t21 ; // Compute Eigenvalues ( lmax , lmin )
a32 = a31* t12 + a32 *t22 ; //=squares of singular values ( smax , smin )
sortie[0]= a11; double trace = d11 + d22 ;
sortie[1]= a12; double det = d11 * d22 - d12 * d12 ;
sortie[2]= a21; double d = sqrt( max( 0.0 , 0.25 * trace * trace - det ) ) ;
sortie[3]= a22; double lmax = max( 0.0 , 0.5 * trace + d ) ;
sortie[4]= a31; double lmin = max( 0.0 , 0.5 * trace - d ) ;
sortie[5]= a32; double smax = sqrt ( lmax ) ;
double smin = sqrt ( lmin ) ;
} // Check whether matrix lies outside K and need o be projed
if ( smax + smin > 1.0 ) {
// Compute orthonormalsy st emofEigenvectors , such t h a t
// ( v11 , v21 ) belongs to lmax and ( v12 , v22 ) belongs to lmin
double v11 , v12 , v21 , v22 ;
if ( d12 == 0.0 ) {
if ( d11 >= d22 ) { v11 = 1.0 ; v21 = 0.0 ; v12 = 0.0 ; v22 = 1.0 ; }
else
{ v11 = 0.0 ; v21 = 1.0 ; v12 = 1.0 ; v22 = 0.0 ; }
}
else {
v11 = lmax - d22 ; v21 = d12 ;
double l1 = hypot( v11 , v21 ) ;
v11 /= l1 ; v21 /= l1 ;
v12 = lmin - d22 ; v22 = d12 ;
double l2 = hypot( v12 , v22 ) ;
v12 /= l2 ; v22 /= l2 ;
}
// Project ( smax , smin ) toline ( 0 , 1 ) + tau * (1 , -1) , 0<=tau <=1.
double tau = 0.5 * ( smax - smin + 1.0 ) ;
double s1 = min ( 1.0 , tau ) ;
double s2 = 1.0 - s1 ;
// Compute \ Sigmaˆ+ * \ Sigma p
s1 /= smax ;
s2 = ( smin > 0.0 ) ? s2 / smin : 0.0 ;
// Compute
double t11= s1 * v11 * v11 + s2 * v12 * v12 ;
double t12= s1 * v11 * v21 + s2 * v12 * v22 ;
double t21= s1 * v21 * v11 + s2 * v22 * v12 ;
double t22= s1 * v21 * v21 + s2 * v22 * v22 ;
// Result
a11 = a11* t11 + a12 *t21 ;
a12 = a11* t12 + a12 *t22 ;
a21 = a21* t11 + a22 *t21 ;
a22 = a21* t12 + a22 *t22 ;
a31 = a31* t11 + a32 *t21 ;
a32 = a31* t12 + a32 *t22 ;
sortie[0]= a11;
sortie[1]= a12;
sortie[2]= a21;
sortie[3]= a22;
sortie[4]= a31;
sortie[5]= a32;
}
} }
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{ {
...@@ -129,8 +129,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) ...@@ -129,8 +129,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
height = dims[0]; height = dims[0];
width = dims[1]; width = dims[1];
int nb_chanel_dir= dims[2]; int nb_chanel_dir= dims[2];
Wr = (const double*)(mxGetPr(prhs[0])); Wr = (const double*)(mxGetPr(prhs[0]));
const int n_dim_result= 3; const int n_dim_result= 3;
...@@ -146,17 +145,20 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) ...@@ -146,17 +145,20 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
double *y= (double*)malloc(nb_chanel_dir* sizeof(double)); double *y= (double*)malloc(nb_chanel_dir* sizeof(double));
/* Copie des vecteurs pour traitement*/ /* Copie des vecteurs pour traitement*/
for(int lin=0; lin< height; lin++){ int N = height * width;
for(int col=0; col< width; col++){ for(int idx= 0; idx < N; idx++){
for(int w=0; w< nb_chanel_dir; w++){ int c=0;
x[w]= Wr[lin + height * col + w * width * height]; for(int w=0; w< nb_chanel_dir; w++){
} x[w]= Wr[idx + c];
c+=N;
projection_6(x,y); }
for(int w=0; w< nb_chanel_dir; w++){ projection_6(x,y);
Wl[lin + height * col + w * width * height]= y[w];
} c=0;
for(int w=0; w< nb_chanel_dir; w++){
Wl[idx + c]= y[w];
c+=N;
} }
} }
...@@ -165,4 +167,3 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) ...@@ -165,4 +167,3 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
return; return;
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment