Skip to content
Snippets Groups Projects
Commit fba928a1 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

New SVD decomposition for mean plane computation. This should solve numerous... 

New SVD decomposition for mean plane computation. This should solve numerous instability problems and bad 2D interpolations.
parent 112b0475
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Adapt/dsvdcmp.cpp
+ 24
− 49
View file @ fba928a1
// $Id: dsvdcmp.cpp,v 1.4 2001-01-08 08:05:39 geuzaine Exp $ // $Id: dsvdcmp.cpp,v 1.5 2001-11-01 09:36:59 geuzaine Exp $
#include <math.h>
#include "Gmsh.h"
#include "nrutil.h" #include "nrutil.h"
double dpythag(double a, double b) double dpythag(double a, double b)
...@@ -136,7 +136,9 @@ void dsvdcmp(double **a, int m, int n, double w[], double **v) ...@@ -136,7 +136,9 @@ void dsvdcmp(double **a, int m, int n, double w[], double **v)
} }
break; break;
} }
if (its == 30) nrerror("no convergence in 30 dsvdcmp iterations"); if (its == 1000){
Msg(GERROR, "SVD decomposition: no convergence in 1000 iterations");
}
x=w[l]; x=w[l];
nm=k-1; nm=k-1;
y=w[nm]; y=w[nm];
...@@ -190,51 +192,24 @@ void dsvdcmp(double **a, int m, int n, double w[], double **v) ...@@ -190,51 +192,24 @@ void dsvdcmp(double **a, int m, int n, double w[], double **v)
free_dvector(rv1,1,n); free_dvector(rv1,1,n);
} }
void dsvbksb(double **u, double w[], double **v, int m, int n, double b[], double x[])
{
int jj,j,i;
double s,*tmp;
/* cf. Numerical Recipes in C, p. 62 */ tmp=dvector(1,n);
#define PREC 1.e-16
void invert_singular_matrix(double **M, int n, double **I){
double **V, **T, *W;
int i, j, k;
V = dmatrix(1,n,1,n);
T = dmatrix(1,n,1,n);
W = dvector(1,n);
dsvdcmp(M, n, n, W, V);
for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){
I[i][j] = 0.0 ;
T[i][j] = 0.0 ;
}
}
for(i=1 ; i<=n ; i++){
for (j=1;j<=n;j++) { for (j=1;j<=n;j++) {
if(fabs(W[i]) > PREC){ s=0.0;
T[i][j] += M[j][i] / W[i] ; if (w[j]) {
} for (i=1;i<=m;i++) s += u[i][j]*b[i];
/* s /= w[j];
else{
T[i][j] += 0.0 ;
}
*/
} }
tmp[j]=s;
} }
for(i=1 ; i<=n ; i++){
for (j=1;j<=n;j++) { for (j=1;j<=n;j++) {
for(k=1 ; k<=n ; k++){ s=0.0;
I[i][j] += V[i][k] * T[k][j] ; for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj];
} x[j]=s;
} }
free_dvector(tmp,1,n);
} }
free_dmatrix(V,1,n,1,n);
free_dmatrix(T,1,n,1,n);
free_dvector(W,1,n);
}
#undef PREC
Mesh/Utils.cpp
+ 160
− 213
View file @ fba928a1
// $Id: Utils.cpp,v 1.2 2001-08-12 20:45:58 geuzaine Exp $ // $Id: Utils.cpp,v 1.3 2001-11-01 09:36:59 geuzaine Exp $
#include "Gmsh.h" #include "Gmsh.h"
#include "Numeric.h" #include "Numeric.h"
...@@ -11,6 +11,9 @@ ...@@ -11,6 +11,9 @@
extern Context_T CTX; extern Context_T CTX;
void dsvdcmp(double **a, int m, int n, double w[], double **v);
void dsvbksb(double **u, double w[], double **v, int m, int n, double b[], double x[]);
void direction (Vertex * v1, Vertex * v2, double d[3]){ void direction (Vertex * v1, Vertex * v2, double d[3]){
d[0] = v2->Pos.X - v1->Pos.X; d[0] = v2->Pos.X - v1->Pos.X;
d[1] = v2->Pos.Y - v1->Pos.Y; d[1] = v2->Pos.Y - v1->Pos.Y;
...@@ -28,240 +31,147 @@ void Projette (Vertex * v, double mat[3][3]){ ...@@ -28,240 +31,147 @@ void Projette (Vertex * v, double mat[3][3]){
v->Pos.Z = Z; v->Pos.Z = Z;
} }
/*
Le concept d'un plan moyen calcule au sens des moidres carres n'est
pas le bon. Imagine un quart de cercle extrude d'une faible
hauteur. Le plan moyen sera dans le plan du cercle! En attendant
mieux, il y a un test de coherence pour les surfaces non-planes. */
void MeanPlane(List_T *points, Surface *s){ void MeanPlane(List_T *points, Surface *s){
int i, j, ix, iy, iz, N; int i, j, min, ndata, na;
double det,sys[3][3],b[3],res[3],mod,t1[3],t2[3],ex[3],s2s[2][2],r2[2],X,Y,Z; double **U, **V, *W, res[4], ex[3], t1[3], t2[3];
Vertex *v; Vertex *v;
double xm=0., ym=0., zm=0.;
N = List_Nbr (points); ndata = List_Nbr(points);
na = 3;
for (i = 0; i < 3; i++){
b[i] = 0.0;
for (j = 0; j < 3; j++){
sys[i][j] = 0.0;
}
}
/* ax + by + cz = 1 */
ix = iy = iz = 0;
// TOLERANCE ! WARNING WARNING U = dmatrix(1,ndata,1,na);
double eps = 1.e-6 * CTX.lc; V = dmatrix(1,ndata,1,ndata);
W = dvector(1,na);
for (i = 0; i < N; i++){ for (i=0; i<ndata; i++){
List_Read(points, i, &v); List_Read(points, i, &v);
xm += v->Pos.X;
if (!i){ ym += v->Pos.Y;
X = v->Pos.X; zm += v->Pos.Z;
Y = v->Pos.Y;
Z = v->Pos.Z;
}
else{
if(fabs(X-v->Pos.X) > eps) ix = 1;
if(fabs(Y-v->Pos.Y) > eps) iy = 1;
if(fabs(Z-v->Pos.Z) > eps) iz = 1;
}
sys[0][0] += v->Pos.X * v->Pos.X;
sys[1][1] += v->Pos.Y * v->Pos.Y;
sys[2][2] += v->Pos.Z * v->Pos.Z;
sys[0][1] += v->Pos.X * v->Pos.Y;
sys[0][2] += v->Pos.X * v->Pos.Z;
sys[1][2] += v->Pos.Y * v->Pos.Z;
sys[2][1] = sys[1][2];
sys[1][0] = sys[0][1];
sys[2][0] = sys[0][2];
b[0] += v->Pos.X;
b[1] += v->Pos.Y;
b[2] += v->Pos.Z;
}
s->d = 1.0;
/* x = X */
if (!ix){
s->d = X;
res[0] = 1.;
res[1] = res[2] = 0.0;
Msg(DEBUG, "Mean plane of type 'x = c'");
}
/* y = Y */
else if (!iy){
s->d = Y;
res[1] = 1.;
res[0] = res[2] = 0.0;
Msg(DEBUG, "Mean plane of type 'y = c'");
}
/* z = Z */
else if (!iz){
s->d = Z;
res[2] = 1.;
res[1] = res[0] = 0.0;
Msg(DEBUG, "Mean plane of type 'z = c'");
}
/* by + cz = -x */
else if (!sys3x3_with_tol (sys, b, res, &det)){
s->d = 0.0;
s2s[0][0] = sys[1][1];
s2s[0][1] = sys[1][2];
s2s[1][0] = sys[1][2];
s2s[1][1] = sys[2][2];
b[0] = -sys[0][1];
b[1] = -sys[0][2];
if (sys2x2 (s2s, b, r2)){
res[0] = 1.;
res[1] = r2[0];
res[2] = r2[1];
Msg(DEBUG, "Mean plane of type 'by + cz = -x'");
} }
xm/=(double)ndata;
ym/=(double)ndata;
zm/=(double)ndata;
/* ax + cz = -y */ for (i=0; i<ndata; i++){
List_Read(points, i, &v);
else{ U[i+1][1] = v->Pos.X-xm;
s->d = 0.0; U[i+1][2] = v->Pos.Y-ym;
s2s[0][0] = sys[0][0]; U[i+1][3] = v->Pos.Z-zm;
s2s[0][1] = sys[0][2]; }
s2s[1][0] = sys[0][2]; dsvdcmp(U,ndata,na,W,V);
s2s[1][1] = sys[2][2]; if(W[1]<W[2] && W[1]<W[3]) min=1;
b[0] = -sys[0][1]; else if(W[2]<W[1] && W[2]<W[3]) min=2;
b[1] = -sys[1][2]; else min=3;
if (sys2x2 (s2s, b, r2)){ res[0] = V[1][min];
res[0] = r2[0]; res[1] = V[2][min];
res[1] = 1.; res[2] = V[3][min];
res[2] = r2[1]; norme(res);
Msg(DEBUG, "Mean plane of type 'ax + cz = -y'");
} free_dmatrix(U,1,ndata,1,na);
free_dmatrix(V,1,ndata,1,ndata);
/* ax + by = -z */ free_dvector(W,1,na);
else{ // check coherence of results for non-plane surfaces
s->d = 1.0; if(s->Typ != MSH_SURF_PLAN){
s2s[0][0] = sys[0][0]; double res2[3], c[3], cosc, sinc, angplan;
s2s[0][1] = sys[0][1]; double eps=1.e-3;
s2s[1][0] = sys[0][1]; Vertex v1, v2, v3;
s2s[1][1] = sys[1][1]; v1 = InterpolateSurface(s,0.5,0.5,0,0);
b[0] = -sys[0][2]; v2 = InterpolateSurface(s,0.5+eps,0.5,0,0);
b[1] = -sys[1][2]; v3 = InterpolateSurface(s,0.5,0.5+eps,0,0);
if (sys2x2 (s2s, b, r2)){ t1[0] = v2.Pos.X-v1.Pos.X;
res[0] = r2[0]; t1[1] = v2.Pos.Y-v1.Pos.Y;
res[1] = r2[1]; t1[2] = v2.Pos.Z-v1.Pos.Z;
res[2] = 1.; t2[0] = v3.Pos.X-v1.Pos.X;
Msg(DEBUG, "Mean plane of type 'ax + by = -z'"); t2[1] = v3.Pos.Y-v1.Pos.Y;
} t2[2] = v3.Pos.Z-v1.Pos.Z;
else{ norme(t1);
Msg(GERROR, "Problem in mean plane computation"); norme(t2);
} prodve(t1, t2, res2);
} norme(res2);
prodve (res, res2, c);
prosca (res, res2, &cosc);
sinc = sqrt (c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
angplan = myatan2 (sinc, cosc);
angplan = angle_02pi(angplan)*180./Pi;
if((angplan>70 && angplan<110) ||
(angplan>260 && angplan<280)){
Msg(WARNING, "SVD failed (angle=%g): using rough algo...", angplan);
res[0]=res2[0];
res[1]=res2[1];
res[2]=res2[2];
goto end;
} }
} }
s->a = res[0];
s->b = res[1];
s->c = res[2];
mod = sqrt (res[0] * res[0] + res[1] * res[1] + res[2] * res[2]);
for (i = 0; i < 3; i++)
res[i] /= mod;
/* L'axe n'est pas l'axe des x */
ex[0] = ex[1] = ex[2] = 0.0; ex[0] = ex[1] = ex[2] = 0.0;
if(res[0] == 0.0) if(res[0]==0.)
ex[0] = 1.0; ex[0] = 1.0;
else if(res[1] == 0.0) else if(res[1]==0.)
ex[1] = 1.0; ex[1] = 1.0;
else else
ex[2] = 1.0; ex[2] = 1.0;
prodve(res, ex, t1); prodve(res, ex, t1);
norme(t1);
mod = sqrt (t1[0] * t1[0] + t1[1] * t1[1] + t1[2] * t1[2]);
for (i = 0; i < 3; i++)
t1[i] /= mod;
prodve (t1, res, t2); prodve (t1, res, t2);
norme(t2);
mod = sqrt (t2[0] * t2[0] + t2[1] * t2[1] + t2[2] * t2[2]); end:
for (i = 0; i < 3; i++) res[3] = (xm*res[0]+ym*res[1]+zm*res[2]);
t2[i] /= mod;
for (i = 0; i < 3; i++) s->a = res[0];
s->plan[0][i] = t1[i]; s->b = res[1];
for (i = 0; i < 3; i++) s->c = res[2];
s->plan[1][i] = t2[i]; s->d = res[3];
for (i = 0; i < 3; i++) for(i=0; i<3; i++) s->plan[0][i] = t1[i];
s->plan[2][i] = res[i]; for(i=0; i<3; i++) s->plan[1][i] = t2[i];
for(i=0; i<3; i++) s->plan[2][i] = res[i];
Msg(DEBUG1, "Plane : (%g x + %g y + %g z = %g)", s->a, s->b, s->c, s->d); Msg(DEBUG1, "Surface: %d", s->Num);
Msg(DEBUG2, "Normal : (%g , %g , %g )", res[0], res[1], res[2]); Msg(DEBUG2, "SVD : %g,%g,%g (min=%d)", W[1],W[2],W[3],min);
Msg(DEBUG2, "Plane : (%g x + %g y + %g z = %g)", s->a, s->b, s->c, s->d);
Msg(DEBUG2, "Normal : (%g , %g , %g )", s->a, s->b, s->c);
Msg(DEBUG2, "t1 : (%g , %g , %g )", t1[0], t1[1], t1[2]); Msg(DEBUG2, "t1 : (%g , %g , %g )", t1[0], t1[1], t1[2]);
Msg(DEBUG3, "t2 : (%g , %g , %g )", t2[0], t2[1], t2[2]); Msg(DEBUG3, "t2 : (%g , %g , %g )", t2[0], t2[1], t2[2]);
/* Matrice orthogonale */
if (!iz){
for (i = 0; i < 3; i++){
for (j = 0; j < 3; j++){
s->invplan[i][j] = (i == j) ? 1. : 0.;
s->plan[i][j] = (i == j) ? 1. : 0.;
}
}
}
else{
for(i=0;i<3;i++){ for(i=0;i<3;i++){
for(j=0;j<3;j++){ for(j=0;j<3;j++){
s->invplan[i][j] = s->plan[j][i]; s->invplan[i][j] = s->plan[j][i];
} }
} }
}
//check coherence for plane surfaces
// this is the end of the algo as it was used for surface drawing: if(s->Typ==MSH_SURF_PLAN){
Curve *c;
#if 0 for(i=0;i<List_Nbr(s->Generatrices);i++) {
/* L'axe n'est pas l'axe des x */ List_Read(s->Generatrices,i,&c);
if(res[0] > res[1]){ if(c->Num>0){
ex[0] = 0.; List_Read(s->Generatrices,i,&c);
ex[1] = 1.; for(j=0;j<List_Nbr(c->Control_Points);j++){
ex[2] = 0.; List_Read(c->Control_Points,j,&v);
double d = s->a*v->Pos.X+s->b*v->Pos.Y+s->c*v->Pos.Z-s->d;
if(fabs(d)>CTX.lc*1.e-3){
Msg(GERROR1, "Plane surface %d (%gx+%gy+%gz+%g=0) is not plane!",
s->Num, s->a, s->b, s->c, s->d);
Msg(GERROR3, "Control point %d = (%g,%g,%g), val=%g",
v->Num, v->Pos.X, v->Pos.Y, v->Pos.Z, d);
return;
} }
else{
ex[0] = 1.;
ex[1] = 0.;
ex[2] = 0.;
} }
prodve(res,ex,t1);
mod = sqrt (t1[0] * t1[0] + t1[1] * t1[1] + t1[2] * t1[2] ) ;
for(i=0;i<3;i++) t1[i]/=mod;
prodve(t1,res,t2);
mod = sqrt (t2[0] * t2[0] + t2[1] * t2[1] + t2[2] * t2[2] ) ;
for(i=0;i<3;i++) t2[i]/=mod;
for(i=0;i<3;i++)s->plan[0][i] = t1[i];
for(i=0;i<3;i++)s->plan[1][i] = t2[i];
for(i=0;i<3;i++)s->plan[2][i] = res[i];
/* Matrice orthogonale */
for(i=0;i<3;i++){
for(j=0;j<3;j++){
s->invplan[i][j] = s->plan[j][i];
} }
} }
#endif
} }
}
#define Precision 1.e-10 #define Precision 1.e-10
...@@ -300,7 +210,45 @@ void find_bestuv (Surface * s, double X, double Y, ...@@ -300,7 +210,45 @@ void find_bestuv (Surface * s, double X, double Y,
*V = minv; *V = minv;
} }
void invert_singular_matrix(double **M, int n, double **I); // Numerical Recipes in C, p. 62
void invert_singular_matrix(double **M, int n, double **I){
double **V, **T, *W;
int i, j, k;
V = dmatrix(1,n,1,n);
T = dmatrix(1,n,1,n);
W = dvector(1,n);
dsvdcmp(M, n, n, W, V);
for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){
I[i][j] = 0.0 ;
T[i][j] = 0.0 ;
}
}
#define PREC 1.e-16 //singular value precision
for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){
if(fabs(W[i]) > PREC){
T[i][j] += M[j][i] / W[i] ;
}
}
}
#undef PREC
for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){
for(k=1 ; k<=n ; k++){
I[i][j] += V[i][k] * T[k][j] ;
}
}
}
free_dmatrix(V,1,n,1,n);
free_dmatrix(T,1,n,1,n);
free_dvector(W,1,n);
}
void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V) { void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V) {
double Unew,Vnew,err; double Unew,Vnew,err;
...@@ -529,7 +477,6 @@ double angle_plan (Vertex * V, Vertex * P1, Vertex * P2, double n[3]){ ...@@ -529,7 +477,6 @@ double angle_plan (Vertex * V, Vertex * P1, Vertex * P2, double n[3]){
prosca (PA, PB, &cosc); prosca (PA, PB, &cosc);
prosca (c, n, &sinc); prosca (c, n, &sinc);
angplan = myatan2 (sinc, cosc); angplan = myatan2 (sinc, cosc);
return angplan; return angplan;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment