Régression parabolique non triviale

[leon1789]

on part sur la nouvelle idée de Ben ...

[sylvain231]

et si vous cherchiez la distance d'un point avec la parabole d'équation 0 = a * (cos(theta) * (x-sommet.x) + sin(theta) * (y-sommet.y))^2- sin(theta) * (x - sommet.x) - cos(theta) * (y - sommet.y) ?

[sylvain231]

il peut être utile aussi de chercher sur google si d'autres n'ont pas eu le même problème :
un problème semblable : https://fr.mathworks.com/matlabcentral/ ... ate-system
c'est tout ce que j'ai trouvé ! il n'y a rien sur Google ! le fitting de cercles ou d'ellipses est bien documenté mais pas celui de paraboles bizarre !

[leon1789]

et si vous cherchiez la distance d'un point avec la parabole d'équation 0 = a * (cos(theta) * (x-sommet.x) + sin(theta) * (y-sommet.y))^2- sin(theta) * (x - sommet.x) - cos(theta) * (y - sommet.y) ?

C'est pas simple du tout ! Si on savait faire, on l'aurait fait depuis longtemps.
Voire un majorant "pas trop mauvais" de la distance, c'est ce que j'ai voulu faire juste avant.

[sylvain231]

ok donc on laisse tomber cette façon de faire, peux-tu lire le lien que je t'ai envoyé c'est une solution apparemment fonctionnelle pour fitter une parabole à un ensemble de points mais le sommet n'est pas connu c'est la seule différence avec nous mais peut-être que ça peut vous inspirer

[Ben314]

Concernant la distance minimale d'un point à une parabole donnée, au niveau mathématique, c'est béte comme choux à écrire, sauf que ça conduit à la résolution d'une équation du troisième degré qui demande pas mal de temps de calcul (et vu que c'est pas linéaire du tout, pas trop moyen d'utiliser les moindres carrés . . .)
Donc si on doit le faire pour tout les points, toutes un panel de valeurs du coeff. a (de y=ax^2) et tout un panel d'angles, c'est plus que super lourd . . .

[sylvain231]

oui c'est du brute force mais je peux essayer peux-tu me donner la formule de cette distance STP ?

[sylvain231]

on a pris 1000 theta, je pense qu'on peut en prendre moins et avec une échelle logarithmique des a je pense qu'on diminue déjà le temps de calcul

[sylvain231]

une équation du troisième degré ça se résout facilement avec les équations de Cardan
https://fr.wikipedia.org/wiki/M%C3%A9thode_de_Cardan

[Ben314]

Bon, sinon, la méthode ci dessus, au niveau implémentation, ça donnerais ça :
A minimiser sous contrainte (en considérant theta connu) :

avec, (comme précédemment),



Sauf erreur, le minimum (sous contrainte) de est alors

qui est une fonction de theta qu'il faut minimiser (avec une boucle par exemple).
Et il est obtenu pour un rapport (la parabole ayant alors pour équation

[sylvain231]

et le a (l'ouverture de la parabole) il est égal à quoi ?

[Ben314]

et le a (l'ouverture de la parabole) il est égal à quoi ?

C'est le du message.

Et si ça t'intéresse, le carré de distance d'un point (variable) situé sur la parabole à un point (fixe) c'est :

Et pour trouver le minimum de cette fonction, il "suffit" de chercher les zéros de la dérivée, c'est à dire les tels que puis de regarder parmi les solutions réelles (sauf cas très particulier, il y a 1 ou 3 solutions) laquelle donne la plus petite valeur pour .
Et la racine carrée de cette valeur minimale, c'est la "vraie" distance du point à la parabole.

[leon1789]

une équation du troisième degré ça se résout facilement avec les équations de Cardan
https://fr.wikipedia.org/wiki/M%C3%A9thode_de_Cardan

Ces formules passent par les nombres complexes, ensuite on aura 3 candidats possibles, il faudra retenir le bon (cf le message de Ben juste avant), puis chercher à minimiser une somme de façon brute (sans formule pour éliminer des inconnues), avec un ou plusieurs inconnues dans un domaine réel non borné, etc. C'est difficilement envisageable.

[sylvain231]

à p si je comprend bien
non ça ne marche pas du tout du moins sur la parabole serrée voici mon code :

Code: Tout sélectionner

double E(vector <Point2d> points, int a, int b) {
    double res = 0;
    int n = points.size();
    for (int i = 0; i < n; i++) {
        res += pow(points[i].x, a) * pow(points[i].y, b);
    }
    return res / n;
}

void Parabole::fit(std::vector<cv::Point2d> points, cv::Point2d sommet) {
    this->sommet = sommet;
    int n = points.size();
    cout << "sommet = " << sommet << endl;
    points = translatePoint2ds(points, sommet);
    double E40 = E(points, 4, 0);
    double E31 = E(points, 3, 1);
    double E22 = E(points, 2, 2);
    double E13 = E(points, 1, 3);
    double E04 = E(points, 0, 4);
    double E21 = E(points, 2, 1);
    double E30 = E(points, 3, 0);
    double E12 = E(points, 1, 2);
    double E03 = E(points, 0, 3);
    double E02 = E(points, 0, 2);
    double E11 = E(points, 1, 1);
    double E20 = E(points, 2, 0);
    double minimum = std::numeric_limits<double>::infinity();
    for (double theta = 0; theta <= CV_PI; theta += CV_PI / 1000) {
        double c = cos(theta);
        double s = sin(theta);
        double alpha = E40 * pow(c, 4) - 4 * E31 * pow(c, 3) * s + 6 * E22 * pow(c, 2) * pow(s, 2)
            - 4 * E13 * c * pow(s, 3) + E04 * pow(s, 4);
        double beta = E21 * pow(c, 3) + (E30 - 2 * E12) * pow(c, 2) * s +
            (E03 - 2*E21) * c * pow(s, 2) + E12 * pow(s, 3);
        double gamma = E02 * pow(c, 2) + 2 * E11 * c * s + E20 * pow(s, 2);
        double Mmin = 0.5 * (alpha + gamma - sqrt(pow(alpha - gamma, 2) + 4 * pow(beta, 2)));
        if (Mmin < minimum) {
            minimum = Mmin;
            this->theta_radians = theta;
        }
    }
    double c = cos(this->theta_radians);
    double s = sin(this->theta_radians);
    double alpha = E40 * pow(c, 4) - 4 * E31 * pow(c, 3) * s + 6 * E22 * pow(c, 2) * pow(s, 2)
        - 4 * E13 * c * pow(s, 3) + E04 * pow(s, 4);
    double beta = E21 * pow(c, 3) + (E30 - 2 * E12) * pow(c, 2) * s +
        (E03 - 2*E21) * c * pow(s, 2) + E12 * pow(s, 3);
    double gamma = E02 * pow(c, 2) + 2 * E11 * c * s + E20 * pow(s, 2);
    this->a = (1 / (2 * beta)) * (alpha - gamma - sqrt(pow(alpha - gamma, 2) + 4 * pow(beta, 2)));
    this->theta = 180 * theta_radians / CV_PI;
}


[sylvain231]

une équation du troisième degré ça se résout facilement avec les équations de Cardan
https://fr.wikipedia.org/wiki/M%C3%A9thode_de_Cardan

Ces formules passent par les nombres complexes, ensuite on aura 3 candidats possibles, il faudra retenir le bon (cf le message de Ben juste avant), puis chercher à minimiser une somme de façon brute (sans formule pour éliminer des inconnues), avec un ou plusieurs inconnues dans un domaine réel non borné, etc. C'est difficilement envisageable.

OK le problème semble très difficile alors car cela fait trois jours qu'on est dessus et toujours pas de solution !

[sylvain231]

j'essaierai demain de faire le brute force mais là je fatigue

[leon1789]

Bon, sinon, la méthode ci dessus, au niveau implémentation, ça donnerais ça :
A minimiser sous contrainte (en considérant theta connu) :

avec, (comme précédemment),



Sauf erreur, le minimum (sous contrainte) de est alors

qui est une fonction de theta qu'il faut minimiser (avec une boucle par exemple).
Et il est obtenu pour un rapport (la parabole ayant alors pour équation

Est-ce que ce n'est pas l'opposé pour p = a/b ?

Et je viens de tester : ça fonctionne pour la parabole plate. OK.

Pour la serrée, il y a deux lieux pour le minimum ... pas facile de faire le bon choix.

Le minimum de gauche, c'est le bon t pour réaliser un bon ajustement.
Mais malheureusement, le minimum global est à droite, et on retombe dans le même souci qu'avant.

[Ben314]

Est-ce que ce n'est pas l'opposé pour p = a/b ?

Tout à fait, si !!!!
(sur mon papier, j'avais pris avec un pour le coeff. en ab)

Et j'ai modifié le message contenant l'erreur au cas où quelqu'un tombe dessus . . .

[Ben314]

Et si la courbe que tu as effacé représente le cas général (des fois où ça déconne), je me demande si le "bon" minimum, il n'est pas forcément encadré "d'assez prés" par des valeurs très grandes.
Peut-être qu'en repérant les min. locaux avec le "plus gros creux", ça marcherais.

Mais bon, le moins qu'on puisse dire, c'est que c'est empirique . . .

[leon1789]

Bon, pour résumer la journée :

minimiser fonctionne pour la parabole plate, pas la serrée.

minimiser fonctionne pour la parabole serrée, pas la plate.