Publicité
+ Répondre à la discussion
Affichage des résultats 1 à 15 sur 15
  1. #1
    Rédacteur/Modérateur
    Avatar de pseudocode
    Homme Profil pro Xavier Philippeau
    Architecte système
    Inscrit en
    décembre 2006
    Messages
    9 960
    Détails du profil
    Informations personnelles :
    Nom : Homme Xavier Philippeau
    Âge : 41
    Localisation : France, Hérault (Languedoc Roussillon)

    Informations professionnelles :
    Activité : Architecte système
    Secteur : Industrie

    Informations forums :
    Inscription : décembre 2006
    Messages : 9 960
    Points : 15 759
    Points
    15 759

    Par défaut [java] Fitting de courbe (moindres carrés)

    Le "fitting" de courbe (ou "ajustement" en francais) a pour objectif de trouver la meilleure courbe passant par un ensemble de points.

    Par "meilleure courbe" on entend celle qui minimise une erreur. Dans notre cas, cette erreur est l'ecart quadratique entre la courbe et les points.

    Lorsque la courbe recherchée est une forme linéaire, la methode des "moindres carrés" permet de calculer les coefficients de la courbe en résolvant un système linéraire.

    Le code:

    La classe Matrix, qui s'occupe de stocker et inverser une matrice:
    Code java :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
     
    /**
     * Structure et Operations sur une matrice
     * 
     * @author Xavier PHILIPPEAU
     *  
     */
    public class Matrix {
     
    	// les données de la matrice ordonnnée en [x=colonne][y=ligne]
    	public double value[][];
     
    	// les dimensions de la matrice
    	public final int rows, cols;
     
    	// constructeur de matrice quelconque (vide)
    	public Matrix(int cols, int rows) {
    		this.cols = cols;
    		this.rows = rows;
    		this.value = new double[this.cols][this.rows];
    	}
     
    	// constructeur de matrice carré (pré-remplies)
    	public Matrix(int size, double[][] squarearray) {
    		this.cols = size;
    		this.rows = size;
    		this.value = squarearray;
    	}
     
    	/**
    	 * inversion par pivot de Gauss
    	 *
    	 * @return l'inverse de la matrice (ou null si non inversible)
    	 */
    	public Matrix inverse() {
     
    		//Une matrice ne peut etre inversée que si elle est carrée.
    		if (this.rows != this.cols) {
    			System.err.println("Inversion : La matrice n'est pas carré !!");
    			return null;
    		}
     
    		// Création de la matrice de travail T = [ this | Identité ]
    		Matrix T = new Matrix(this.cols * 2, this.cols);
    		for (int y = 0; y < this.rows; y++) {
    			for (int x = 0; x < this.cols; x++) {
    				T.value[x][y] = this.value[x][y];
    				if (x == y)	T.value[this.cols + x][y] = 1;
    			}
    		}
     
    		// Pour chaque ligne de la matrice T
    		for (int x = 0; x < T.rows; x++) {
     
    			// Cherche la ligne avec le pivot max (en valeur absolue)
    			int bestline = x;
    			double pivot = T.value[x][bestline];
    			for (int y = x + 1; y < T.rows; y++) {
    				if (Math.abs(T.value[x][y]) > Math.abs(pivot)) {
    					bestline = y;
    					pivot = T.value[x][bestline];
    				}
    			}
     
    			if (pivot == 0) {
    				System.err.println("Inversion : Le pivot est nul,inversion impossible !!");
    				return null;
    			}
     
    			// Echange des lignes (si necessaire)
    			if (bestline != x) {
    				double tmp;
    				for (int t = 0; t < T.cols; t++) {
    					tmp = T.value[t][x]; T.value[t][x] = T.value[t][bestline]; T.value[t][bestline] = tmp;
    				}
    			}
     
    			// Normalisation de la ligne du pivot
    			for (int t = 0; t < T.cols; t++) T.value[t][x] /= pivot;
     
    			// elimination des autres lignes
    			for (int y = 0; y < T.rows; y++) {
    				if (y==x) continue;
    				double coef = T.value[x][y];
    				for (int t = 0; t < T.cols; t++) T.value[t][y] -= coef * T.value[t][x];
    			}
    		}
     
    		// recupere la partie droite de T qui contient l'inverse de la matrice
    		// (la partie gauche devrait contenir l'identité)
    		Matrix inverse = new Matrix(this.cols,this.rows);
    		for (int y = 0; y < this.rows; y++)
    			for (int x = 0; x < this.cols; x++)
    				inverse.value[x][y] = T.value[this.cols+x][y];
     
    		return inverse;
    	}
     
    	/**
    	 * Multiplication de la matrice par un vecteur
    	 *
    	 * @param X le vecteur
    	 * @return le vecteur résultat: this[][] * X[]
    	 */
    	public double[] multiply(double[] X) {
    		double[] Y = new double[this.rows];
    		for (int y = 0; y < this.rows; y++)
    			for (int x = 0; x < this.cols; x++)
    				Y[y] += this.value[x][y]*X[x];
    		return Y;
    	}	
     
    }

    La classe LinearEquation qui s'occupe de modeliser une equation lineaire (coefs + valeur)
    Code java :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
     
    /**
     * Structure d'une equation lineaire
     * 
     *  coef[0]*X0 +  coef[1]*X1 + ... + coef[n]*Xn = residu
     *  
     *  ou X=(X0,X1,...,Xn) est le vecteur "inconnu"
     *
     * @author Xavier PHILIPPEAU
     */
    public class LinearEquation {
     
    	/**
    	 * Les coefficients de l'equation
    	 */
    	public double[] coefficents;
     
     
    	/**
    	 * Le résidu (= le résultat) de l'equation 
    	 */
    	public double residual;
     
    	/**
    	 * Constructeur
    	 * 
    	 * @param unknownsCount Le nombre d'inconnues
    	 */
    	public LinearEquation(int unknownsCount) {
    		this.coefficents = new double[unknownsCount];
    	}
     
    	/**
    	 * Constructeur
    	 * 
    	 * @param coefficents Les coefficients de l'equation
    	 * @param residual Le résidu (= le résultat) de l'equation 
    	 */
    	public LinearEquation(double[] coefficents, double residual) {
    		this.coefficents = coefficents;
    		this.residual = residual;
    	}
     
    }

    La classe LeastSquare qui s'occupe de construire et résoudre le système d'equations lineaires par la méthode des "moindres carrés":
    Code java :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
     
    import java.util.List;
     
    /**
     * résolution du système d'equations par la methode des moindres carrés
     *
     * @author Xavier PHILIPPEAU
     */
    public class LeastSquare {
     
    	/**
    	 * Resolution d'un systeme lineaire
    	 * 
    	 * @param systemsize Nombre d'inconnues dans le systeme
    	 * @param system Le systeme d'equations lineaires
    	 * @return le vecteur solution du système
    	 */
    	public static double[] linearFit(int systemsize, List<LinearEquation> system) {
     
    		// taille du systeme a resoudre (= nombre d'inconnues)
    		int n = systemsize;
     
    		// On cherche a trouver X qui resoud le systeme:
    		//
    		// J.X = R
    		//
    		// on multiplie par Jt, la transposée de J:
    		//
    		// Jt.J.X = Jt.R
    		//
    		// on pose:
    		//                    n
    		// A[c][l] = Jt.J = Somme[ coefficient[c][k] * coefficient[l][k] ]
    		//                   k=0
    		//
    		//                 n
    		// Y[l] = Jt.R = Somme[ residual[k] * coefficient[l][k] ]
    		//                k=0
    		//
    		// et "n", le nombre d'equation dans le systeme.
    		//
    		// Le systeme a resoudre devient ainsi:
    		//
    		// A.X = Y
     
    		Matrix A = new Matrix(n,n);
    		double Y[] = new double[n];
    		for(LinearEquation equation:system) {
    			for (int l = 0; l < n; l++) {
    				for (int c = 0; c < n; c++) {
    					A.value[c][l] += equation.coefficents[c] * equation.coefficents[l];
    				}
    				Y[l] += equation.residual * equation.coefficents[l];
    			}
    		}
     
    		// Resolution par inversion:
    		//
    		// A . X = Y  ==>  X = A^-1 . Y
     
    		double X[] = (A.inverse()).multiply(Y);		
     
    		return X;
    	}
     
    }

    PRomu@ld : sources
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

  2. #2
    Rédacteur/Modérateur
    Avatar de pseudocode
    Homme Profil pro Xavier Philippeau
    Architecte système
    Inscrit en
    décembre 2006
    Messages
    9 960
    Détails du profil
    Informations personnelles :
    Nom : Homme Xavier Philippeau
    Âge : 41
    Localisation : France, Hérault (Languedoc Roussillon)

    Informations professionnelles :
    Activité : Architecte système
    Secteur : Industrie

    Informations forums :
    Inscription : décembre 2006
    Messages : 9 960
    Points : 15 759
    Points
    15 759

    Par défaut Exemples d'utilisation

    Trouver la "meilleure ligne" au sens des moindres carrés:
    Code java :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
     
    public static double[] bestLine(List<double[]> data) {
     
    	List<LinearEquation> system = new ArrayList<LinearEquation>();
     
    	for(double[] d : data) {
    		double x = d[0];
    		double y = d[1];
     
    		// a.1 + b.X = Y
    		//
    		// (X,Y) sont connus (données du probleme)
    		// (a,b) sont les 2 inconnues
     
    		LinearEquation equation = new LinearEquation(2);
    		equation.coefficents = new double[] { 1, x };
    		equation.residual = y;
     
    		system.add(equation);
    	}
     
    	// resolution du systeme lineaire
    	double[] C = LeastSquare.linearFit(2, system);
    	return C;
     
    	// La solution s'ecrit: Y = C[0] + C[1]*X
    }

    Trouver la "meilleure parabole" au sens des moindres carrés:
    Code java :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
     
    public static double[] BestParabol(List<double[]> data) {
     
    	List<LinearEquation> system = new ArrayList<LinearEquation>();
     
    	for(double[] d : data) {
    		double x = d[0];
    		double y = d[1];
     
    		// a.1 + b.X + c.X^2 = Y
    		//
    		// (X,Y) sont connus (données du probleme)
    		// (a,b,c) sont les 3 inconnues
     
    		LinearEquation equation = new LinearEquation(3);
    		equation.coefficents = new double[] { 1, x, x*x };
    		equation.residual = y;
     
    		system.add(equation);
    	}
     
    	// resolution du systeme lineaire
    	double[] C = LeastSquare.linearFit(3, system);
    	return C;
     
    	// La solution s'ecrit: Y = C[0] + C[1]*X + C[2]*X^2
    }

    Trouver le "meilleure polynome de degré n" au sens des moindres carrés:
    Code java :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
     
    public static double[] BestPolynom(int order, List<double[]> data) {
     
    	List<LinearEquation> system = new ArrayList<LinearEquation>();
     
    	for(double[] d : data) {
    		double x = d[0];
    		double y = d[1];
     
    		// a.1 + b.X + c.X^2 + ... + z.X^ordre = Y
    		//
    		// (X,Y) sont connus (données du probleme)
    		// (a,b,c,...z) sont les inconnues
     
    		LinearEquation equation = new LinearEquation(order+1);
     
    		equation.residual = y;
    		for(int n=0;n<=order;n++)
    			equation.coefficents[n] = Math.pow(x, n);
     
    		system.add(equation);
    	}
     
    	// resolution du systeme lineaire
    	double[] C = LeastSquare.linearFit(order+1, system);
    	return C;
     
    	// La solution s'ecrit: Y = C[0] + C[1]*X + C[2]*X^2 + ... + C[n]*X^n
    }

    un exemple/explication du cas "non lineaire" appliqué à un cercle:
    http://www.orbitals.com/self/least/least.htm
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

  3. #3
    Expert Confirmé Sénior

    Inscrit en
    janvier 2007
    Messages
    10 183
    Détails du profil
    Informations personnelles :
    Âge : 57

    Informations forums :
    Inscription : janvier 2007
    Messages : 10 183
    Points : 14 353
    Points
    14 353

    Par défaut

    je la traduirai en C (Java moi yen a pas connaître) et je posterais le texte du papier cette semaine (peut-être aujourdhui ?) .. (je viens de rentrer de voyage..)
    "Un homme sage ne croit que la moitié de ce qu’il lit. Plus sage encore, il sait laquelle".

    Consultant indépendant.
    Architecture systèmes complexes. Programmation grosses applications critiques. Ergonomie.
    C, Fortran, XWindow/Motif, Java

    Je ne réponds pas aux MP techniques

  4. #4
    Rédacteur/Modérateur
    Avatar de pseudocode
    Homme Profil pro Xavier Philippeau
    Architecte système
    Inscrit en
    décembre 2006
    Messages
    9 960
    Détails du profil
    Informations personnelles :
    Nom : Homme Xavier Philippeau
    Âge : 41
    Localisation : France, Hérault (Languedoc Roussillon)

    Informations professionnelles :
    Activité : Architecte système
    Secteur : Industrie

    Informations forums :
    Inscription : décembre 2006
    Messages : 9 960
    Points : 15 759
    Points
    15 759

    Par défaut

    Trouver la "meilleure ellipse" au sens des moindres carrés, en utilisant les formules de mathworld.wolfram.com:

    Code java :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
     
    public static double[] BestEllipsoid(List<double[]> data) {
     
    	List<LinearEquation> system = new ArrayList<LinearEquation>();
     
    	for(double[] d : data) {
    		double x = d[0];
    		double y = d[1];
     
    		// a.x^2 + 2b.xy + c.y^2 + 2d.x + 2f.y + g = 0
    		//
    		// (X,Y) sont connus (données du probleme)
    		// on pose g = -1
    		// (a,b,c,d,f) sont les 5 inconnues
     
    		LinearEquation equation = new LinearEquation(5);
    		equation.coefficents = new double[] { x*x, x*y, y*y, x, y };
    		equation.residual = 1.0;
     
    		system.add(equation);
    	}
     
    	// resolution du systeme lineaire
    	double[] C = LeastSquare.linearFit(5, system);
     
    	// evaluation des 5 inconnues
    	double a = C[0];
    	double b = C[1]/2;
    	double c = C[2];
    	double d = C[3]/2;
    	double f = C[4]/2;
    	double g = -1;
     
    	// calculs des caracteristiques de l'ellipse:
    	//
    	// d'après http://mathworld.wolfram.com/Ellipse.html
    	//
    	// formules (15), (16), ..., (23)
     
    	// verifie que l'ellipse est valide
    	double delta=a*(c*g-f*f)-b*(b*g-d*f)+d*(b*f-d*c);
    	double J=a*c-b*b;
    	double I=a+c;
    	if(delta==0) return null;
    	if(J<=0) return null;
    	if((delta/I)>=0) return null;
     
    	/* calcul du centre de l'ellipse */
    	double xcentre = (c*d-b*f)/(b*b-a*c);
    	double ycentre = (a*f-b*d)/(b*b-a*c);
     
    	/* angle de rotation (en radian) */
    	double angle = 0.5*Math.atan(2*b/(c-a));
     
    	/* longueurs des 2 demi-axes */
    	double num = 2*(a*f*f + c*d*d + g*b*b - 2*b*d*f - a*c*g);
    	double sq = Math.sqrt( 1 + 4*b*b/((a-c)*(a-c)) );
     
    	double denom1 = (b*b-a*c) * ( (c-a) * sq - (c+a) );
    	double xhalflength = Math.sqrt(num/denom1);
     
    	double denom2 = (b*b-a*c) * ( (a-c) * sq - (c+a) );
    	double yhalflength = Math.sqrt(num/denom2);
     
    	/* retourne toujours un angle positif  */
    	if (angle<0) {
    		angle+=Math.PI/2;
    		// echange les longueurs des demi-axes
    		double tmp=xhalflength;
    		xhalflength=yhalflength;
    		yhalflength=tmp;
    	}
     
    	return new double[] { xcentre, ycentre, angle, xhalflength, yhalflength};
    }

    Et un petit test pour s'assurer que ca marche...
    Code java :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
     
    // construction des points d'une ellipse
    List<double[]> data = new ArrayList<double[]>();
    for(int i=0;i<100;i++) {
    	double t = i/100.0;
    	double a = t*2*Math.PI;
     
    	// ellipse de demi-axe 17 x 10, centrée en (0,0)
    	double xe = 17*Math.cos(a);
    	double ye = 10*Math.sin(a);
     
    	// + une rotation de 45 degrés
    	double r = Math.toRadians(45);
    	double xr = xe*Math.cos(r)+ye*Math.sin(r);
    	double yr = -xe*Math.sin(r)+ye*Math.cos(r);
     
    	// + une translation de (15,8)
    	double x = xr+15;
    	double y = yr+8;
     
    	data.add( new double[] {x,y});
    }
     
    double[] E = BestEllipsoid(data);
     
    System.out.println("centre: "+E[0]+","+E[1]);
    System.out.println("angle: "+Math.toDegrees(E[2])+" (degrés)");
    System.out.println("demi-axe: "+E[3]+" x "+E[4]);
     
    // affiche:
    // centre: 14.999999999999954,8.000000000000052
    // angle: 44.99999999999916 (degrés)
    // demi-axe: 17.00000000000003 x 10.000000000000034
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

  5. #5
    Rédacteur/Modérateur

    Homme Profil pro Jean-Marc Blanc
    Comme retraité, des masses
    Inscrit en
    avril 2007
    Messages
    2 958
    Détails du profil
    Informations personnelles :
    Nom : Homme Jean-Marc Blanc
    Âge : 73
    Localisation : Suisse

    Informations professionnelles :
    Activité : Comme retraité, des masses
    Secteur : Industrie

    Informations forums :
    Inscription : avril 2007
    Messages : 2 958
    Points : 4 901
    Points
    4 901

    Par défaut [java] Fitting de courbe (moindres carrés)

    Deux remarques à l'intention de pseudocode:

    1) Si la tâche à accomplir consiste à résoudre un système linéaire, ce qui est très souvent le cas, il est inutilement coûteux en temps d'exécution d'inverser la matrice (bien que ce soit ce que font la plupart des gens); le plus efficace est de la factoriser par la méthode LU dans le cas général, ou par la méthode de Cholesky, deux fois plus rapide, si elle est symétrique définie positive.

    2) Si tu cherches le vecteur qui satisfait le moins mal possible (au sens des moindres carrés) un système linéaire surdéterminé, regarde la méthode des valeurs singulières (singular value decomposition); tu trouveras toute la documentation à ce sujet dans Numerical Recipes ou dans la description de la fonction svd de MatLab.

    A ta disposition si tu as encore des questions à ce sujet.

  6. #6
    Invité de passage
    Inscrit en
    décembre 2006
    Messages
    3
    Détails du profil
    Informations forums :
    Inscription : décembre 2006
    Messages : 3
    Points : 1
    Points
    1

    Par défaut

    bonjour,
    je souhaite demander si je peux utiliser cette méthode pour résoudre un système de la forme : X^t*A=o , avec A matrice non carrée(n*m)
    merci d'avance

  7. #7
    Rédacteur/Modérateur
    Avatar de pseudocode
    Homme Profil pro Xavier Philippeau
    Architecte système
    Inscrit en
    décembre 2006
    Messages
    9 960
    Détails du profil
    Informations personnelles :
    Nom : Homme Xavier Philippeau
    Âge : 41
    Localisation : France, Hérault (Languedoc Roussillon)

    Informations professionnelles :
    Activité : Architecte système
    Secteur : Industrie

    Informations forums :
    Inscription : décembre 2006
    Messages : 9 960
    Points : 15 759
    Points
    15 759

    Par défaut

    Citation Envoyé par souheyeb Voir le message
    bonjour,
    je veus demander si je peut utiliser cette methode pour resoudre un systeme de la forme : X^t*A=o , avec A matrice non carrée(n*m)
    merci d'avance
    Heu... Non. La méthode des moindres carrés a pour objectif de trouver les "meilleurs" coefficients d'une courbe. La résolution du système linéaire surdéterminé est l'une des étapes de cette méthode, mais nullement le but en soi. D'ailleurs, comme l'a fait remarquer FR119492, elle n'est meme pas obligatoire.

    Pour la résolution d'un système linéaire surdéterminé, je te conseille plutôt la factorisation QR ou la décomposition SVD.
    ALGORITHME (n.m.): Méthode complexe de résolution d'un problème simple.

  8. #8
    Invité de passage
    Inscrit en
    mai 2012
    Messages
    1
    Détails du profil
    Informations forums :
    Inscription : mai 2012
    Messages : 1
    Points : 1
    Points
    1

    Par défaut

    Pour passer en 3d et chercher l'equation d'un plan , il suffit de déclarer
    double z = d[2];

    et
    de chercher l'equation lineaire z=aX+bY+c ?

    UN petit coup de main pour un informaticien non mathematicien ? (pitié!)

  9. #9
    Membre du Club
    Inscrit en
    mars 2009
    Messages
    54
    Détails du profil
    Informations forums :
    Inscription : mars 2009
    Messages : 54
    Points : 44
    Points
    44

    Par défaut Une librairie ?

    Bonjour,

    Merci pour ce sujet, pseudocode.
    On peut aussi utiliser une librairie d'algèbre linéaire pour créer le système linéaire et le résoudre... (ça fait moins de code à écrire et il y a peut-être des améliorations pour la résolution du système quand il est mal conditionné, par rapport à tA.A x = tA.b )
    Jama, par exemple ?
    http://math.nist.gov/javanumerics/jama/

    Exemple :
    Code :
    1
    2
    3
    4
    5
    6
    7
          double[][] array = {{1.,2.,3},{4.,5.,6.},{7.,8.,10.}};
          Matrix A = new Matrix(array);
          Matrix b = Matrix.random(3,1);
          Matrix x = A.solve(b);
          Matrix Residual = A.times(x).minus(b);
          double rnorm = Residual.normInf();

  10. #10
    Membre du Club
    Inscrit en
    mars 2009
    Messages
    54
    Détails du profil
    Informations forums :
    Inscription : mars 2009
    Messages : 54
    Points : 44
    Points
    44

    Par défaut Version avec Jama

    Citation Envoyé par FR119492 Voir le message
    1) Si la tâche à accomplir consiste à résoudre un système linéaire, ce qui est très souvent le cas, il est inutilement coûteux en temps d'exécution d'inverser la matrice (bien que ce soit ce que font la plupart des gens); le plus efficace est de la factoriser par la méthode LU dans le cas général, ou par la méthode de Cholesky, deux fois plus rapide, si elle est symétrique définie positive.
    A priori, pour une résolution "au mieux" (matrice rectangulaire, nb lignes > nb colonnes), Jama utilise une décomposition QR...
    Code Java :
    1
    2
    3
    4
    5
     
       public Matrix solve (Matrix B) {
          return (m == n ? (new LUDecomposition(this)).solve(B) :
                           (new QRDecomposition(this)).solve(B));
       }

    J'ai écrit une version dérivée du code de Xavier, adaptée à l'utilisation de Jama, voilà ce que ça donne :

    Code Java :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
     
    import Jama.Matrix;
    import java.awt.geom.Point2D;
    import java.util.AbstractList;
     
    [...]
     
    	/**
    	 * Renvoie les paramètre de l'ellipse qui s'ajuste mieux au nuage de
    	 *  points 2D reçu (au sens des moindres carrés)
    	 * @param lstPts2D liste de points 2D (Point2D.Double) du contour de l'ellipse
    	 * @return un tableau qui contient les paramètres de l'ellipse
    	 */
    	public static double[] getParametres_ellipse_mc( AbstractList lstPts2D ) {
    		// construction du système linéaire
    		int nNb_pts = lstPts2D.size();
    		Matrix mA = new Matrix( nNb_pts, 5 );
    		Matrix vB = new Matrix( nNb_pts, 1 );
    		for( int i = 0; i < nNb_pts; ++i ) {
    			Point2D.Double pt = (Point2D.Double)lstPts2D.get( i );
    			mA.set( i, 0, pt.x * pt.x );
    			mA.set( i, 1, pt.x * pt.y );
    			mA.set( i, 2, pt.y * pt.y );
    			mA.set( i, 3, pt.x );
    			mA.set( i, 4, pt.y );
    			vB.set(i, 0, -1.0 );
    		}
    		// resolution du systeme lineaire A.x = b
    		Matrix vX = mA.solve( vB );
     
    		// evaluation des 5 inconnues
    		double a = vX.get( 0, 0 );
    		double b = vX.get( 1, 0 ) / 2.0;
    		double c = vX.get( 2, 0 );
    		double d = vX.get( 3, 0 ) / 2.0;
    		double f = vX.get( 4, 0 ) / 2.0;
    		double g = 1.0;
     
    		// calcule les caracteristiques de l'ellipse:
    		//
    		// d'après http://mathworld.wolfram.com/Ellipse.html
    		//
    		// formules (15), (16), ..., (23)
     
    		// verifie que l'ellipse est valide
    		double delta = a * ( c * g - f * f) - b * ( b * g - d * f ) + d * ( b * f - d * c );
    		double nJ = a * c - b * b;
    		double nI = a + c;
    		if ( delta == 0 ) return null;
    		if ( nJ <= 0 ) return null;
    		if ( ( delta / nI ) >= 0) return null;
     
    		// calcul du centre de l'ellipse
    		double nX_centre = ( c*d - b*f )/( b*b - a*c );
    		double nY_centre = ( a*f - b*d )/( b*b - a*c );
     
    		// angle de rotation (en radian)
    		double nAngle_rad = 0.5 * Math.atan( 2.0 * b / ( c-a ));
     
    		// longueurs des 2 demi-axes
    		double num = 2.0 * ( a * f*f + c * d*d + g* b*b - 2*b*d*f - a*c*g);
    		double sq = Math.sqrt( (a-c)*(a-c) + 4.0*b*b );
     
    		double denom1 = ( b*b - a*c ) * ( sq - ( a + c ) );
    		double nDemi_axe_x = Math.sqrt( num / denom1 );
     
    		double denom2 = ( b*b - a*c ) * ( -1.0 * sq - ( a + c ) );
    		double nDemi_axe_y = Math.sqrt( num / denom2 );
     
    		// retourne un angle entre [0 et 2\Pi[
    		if ( nAngle_rad < 0.0 ) {
    			nAngle_rad += 2.0 * Math.PI;
    		}
    		return new double[]{nX_centre, nY_centre, nAngle_rad, nDemi_axe_x, nDemi_axe_y };
    	}// getParametres_ellipse_mc

    Avec la même procédure de test que Xavier j'obtiens aussi les bons résultats :
    Code :
    1
    2
    3
    4
    centre: 14.999999999999991,7.999999999999996
    angle: 44.99999999999999 (degrés)
    demi-axe: 17.000000000000025 x 9.999999999999982
    Merci encore pour ce sujet, il m'a bien aidé...

  11. #11
    Membre du Club
    Inscrit en
    mars 2009
    Messages
    54
    Détails du profil
    Informations forums :
    Inscription : mars 2009
    Messages : 54
    Points : 44
    Points
    44

    Par défaut

    La version "EllipseFitting" de Michael Doube est intéressante (voir le code ici) :
    http://github.com/mdoube/BoneJ/blob/...itEllipse.java

    Voir notamment la fonction "varToDimensions"...

  12. #12
    Expert Confirmé Sénior

    Inscrit en
    janvier 2007
    Messages
    10 183
    Détails du profil
    Informations personnelles :
    Âge : 57

    Informations forums :
    Inscription : janvier 2007
    Messages : 10 183
    Points : 14 353
    Points
    14 353

    Par défaut

    tu peux regarder j'en ai mise une ici-même en Fortran, facilement convertible..
    "Un homme sage ne croit que la moitié de ce qu’il lit. Plus sage encore, il sait laquelle".

    Consultant indépendant.
    Architecture systèmes complexes. Programmation grosses applications critiques. Ergonomie.
    C, Fortran, XWindow/Motif, Java

    Je ne réponds pas aux MP techniques

  13. #13
    Membre Expert
    Homme Profil pro
    Chercheur
    Inscrit en
    mars 2010
    Messages
    1 204
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Chercheur

    Informations forums :
    Inscription : mars 2010
    Messages : 1 204
    Points : 1 645
    Points
    1 645

    Par défaut

    Bonjour,

    inverser une matrice pour ensuite faire
    // Resolution par inversion:
    //
    // A . X = Y ==> X = A^-1 . Y

    double X[] = (A.inverse()).multiply(Y);
    c'est vraiment très maladroit!

  14. #14
    Membre éclairé

    Inscrit en
    janvier 2008
    Messages
    686
    Détails du profil
    Informations forums :
    Inscription : janvier 2008
    Messages : 686
    Points : 395
    Points
    395

    Par défaut

    C'est exactement ce que je cherchais.

    (enfin je cherchais un plane fitting)

    Merci.

  15. #15
    Modérateur
    Avatar de ToTo13
    Homme Profil pro Guillaume
    Ingénieur de Recherche
    Inscrit en
    janvier 2006
    Messages
    5 216
    Détails du profil
    Informations personnelles :
    Nom : Homme Guillaume
    Âge : 35
    Localisation : Etats-Unis

    Informations professionnelles :
    Activité : Ingénieur de Recherche
    Secteur : Santé

    Informations forums :
    Inscription : janvier 2006
    Messages : 5 216
    Points : 8 735
    Points
    8 735

    Par défaut

    En revanche, je pense que cette méthode doit pouvoir être utilisée pour trouver la meilleure gaussienne 2D passant par un ensemble de points.
    Ou ai je tort ?

    [EDIT] : J'ai définitivement tort car l'équation de la Gaussienne n'est pas linéaire à cause de l'exponentielle.
    Donc comment modifier la méthode pour résoudre dans le cas Gaussien ?
    Consignes aux jeunes padawans : une image vaut 1000 mots !
    - Dans ton message respecter tu dois : les règles de rédaction et du forum, prévisualiser, relire et corriger TOUTES les FAUTES (frappes, sms, d'aurteaugrafe, mettre les ACCENTS et les BALISES) => ECRIRE clairement et en Français tu DOIS.
    - Le côté obscur je sens dans le MP => Tous tes MPs je détruirai et la réponse tu n'auras si en privé tu veux que je t'enseigne.(Lis donc ceci)
    - ton poste tu dois marquer quand la bonne réponse tu as obtenu.

Liens sociaux

Règles de messages

  • Vous ne pouvez pas créer de nouvelles discussions
  • Vous ne pouvez pas envoyer des réponses
  • Vous ne pouvez pas envoyer des pièces jointes
  • Vous ne pouvez pas modifier vos messages
  •