Régression parabolique non triviale

[Ben314]

(Re)Salut,
C'est effectivement bien plus simple comme ça.
Quelques remarques :
- Je pense qu'il faut non seulement pas mal de points, mais je pense que pour que ça marche bien il faut un nombre à peu prés équivalent de points sur les deux branches de la parabole et qu'ils soient à peu prés répartis de la même façon par rapport au centre de la parabole (si d'un coté ils sont relativement proche du centre et pas de l'autre ça va déconner) vu que l'axe de la parabole, on le cherche en cherchant la droite qui maximise la somme des carrés des distances à la droite. Je sais pas si ça déconnerais pas aussi si la parabole est très aplati (ou il risque de trouver l'axe des Y comme droite maximisante).
- Les coeff de la matrice de Léon, c'est avec les notations précédentes et tu est pas du tout obligé d'utiliser des matrices pour les calculer.
- Le coeff. a de la parabole (), c'est bien le même que dans mon post du 8 Juin 2023 10:57, sauf que là, de le calculer avec les et les cosinus et sinus, c'est plus vraiment utile vu que, comme il y a un seul calcul à faire et plus des milliers, ça gagne plus rien de développer la formule.

Et sinon, je comprend pas ton message concernant la "parabole tournée" : les formules pour une rotation de centre (0,0), c'est celle de mon message du 15 juin à 13h :
et où et .
Et, évidement, pour changer le centre, il faut faire une translation de (0,0) vers le nouveau centre, c'est à dire ajouter les coordonnées du nouveau centre à X et Y.
Sauf que pour savoir ce que tu fait avec ça, il faudrait savoir sous quelle fore tu veut l'équation de la nouvelle parabole : soit sous forme cartésienne comme dans les différents messages de Léon, soit sous forme paramétrique comme semble le suggérer ton précédent message.

[leon1789]

Quelques remarques :
- Je pense qu'il faut non seulement pas mal de points, mais je pense que pour que ça marche bien il faut un nombre à peu prés équivalent de points sur les deux branches de la parabole et qu'ils soient à peu prés répartis de la même façon par rapport au centre de la parabole (si d'un coté ils sont relativement proche du centre et pas de l'autre ça va déconner) vu que l'axe de la parabole, on le cherche en cherchant la droite qui maximise la somme des carrés des distances à la droite. Je sais pas si ça déconnerais pas aussi si la parabole est très aplati (ou il risque de trouver l'axe des Y comme droite maximisante).

Je suis entièrement d'accord avec toi. Cette méthode a des limites d'utilisation.
Il faut calculer l'axe de la parabole de manière plus sûre.

[leon1789]

sous forme paramétrique comme semble le suggérer ton précédent message.

Ah ce serait bien plus simple pour dessiner la parabole !

Et pour assurer le sommet fixé (S1, S2), atteint disons pour t=0, il suffit de poser : D=B et E=-A


Il reste deux paramètres à déterminer, comme d'hab.

Juste pour info : si le sommet est en (0,0), l'équation implicite correspondante est
A = 4a et B = 4b où (a,b) est le foyer de la parabole.

[sylvain231]

Ok merci. Non je n'ai pas de librairie de calcul matriciel sous la main, il existe eigen mais je n'ai jamais réussi à l'installer.
Question : pour déduire le theta as-ton besoin du sin et du cos ou seulement de l'un des deux et si on a besoin des deux comment on fait ? car je sais que des fois les angles c'est compliqué

[sylvain231]

je vois que vous êtes partis sur une autre méthode, j'attends de voir

[sylvain231]

oui je viens de vérifier pour des paraboles presque plates ça ne marche pas du tout

[sylvain231]

je vous donne mon ensemble de points (sommet : x=200, y=300) :
[491.147, 227.655]
[462.454, 226.516]
[461.672, 221.156]
[460.218, 226.467]
[462.521, 230.463]
[444.12, 221.053]
[441.868, 226.97]
[431.279, 222.414]
[418.907, 231.843]
[405.601, 237.204]
[397.975, 241.513]
[394.867, 229.847]
[374.528, 231.399]
[352.032, 234.525]
[348.5, 239.83]
[339.034, 243.823]
[311.215, 245.222]
[313.473, 258.737]
[293.544, 253.98]
[281.78, 272.452]
[283.286, 259.333]
[269.725, 264.518]
[241.851, 283.36]
[247.822, 272.862]
[225.791, 288.793]
[236.361, 276.347]
[223.803, 282.513]
[209.34, 300.063]
[211.769, 302.712]
[208.739, 295.447]
[201.1, 303.362]
[188.385, 310.076]
[165.963, 317.603]
[167.426, 324.559]
[165.749, 325.947]
[152.545, 331.089]
[147.793, 336.425]
[142.001, 349.469]
[112.913, 358.083]
[103.538, 373.49]
[105.365, 364.987]
[87.9376, 387.12]
[81.5266, 389.234]
[76.6091, 390.403]
[70.4323, 403.594]
[77.5446, 407.94]
[66.4278, 402.054]
[47.4416, 419.374]
[48.0166, 431.14]
[25.2166, 458.267]
[24.2721, 442.828]
[27.8463, 479.668]
[5.51483, 477.176]
[6.52717, 490.404]
[0.633008, 483.01]
[12.2893, 499.955]
[-1.9699, 503.112]

[sylvain231]

j'ai trouvé un code d'inversion de matrice dans le livre "Numerical Recipes Second Edition" un livre qui regroupe des algorithmes de calcul page 48, il est disponible ici :
https://github.com/zhufengGNSS/Numerica ... dition.pdf

[sylvain231]

ah je viens de voir qu'on peut calculer l'inverse d'une matrice avec OpenCV donc c'est bon je suis OK sur l'inversion de matrices :

Code: Tout sélectionner

#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>
#include <iostream>

int main()
{
    // Création d'une matrice carrée
    cv::Mat mat = cv::Mat::zeros(3, 3, CV_32FC1);  // Exemple : matrice nulle

    try
    {
        // Calcul de l'inverse de la matrice
        cv::Mat inverse;
        cv::invert(mat, inverse);

        // Affichage de l'inverse de la matrice
        std::cout << "Inverse de la matrice :\n" << inverse << std::endl;
    }
    catch (cv::Exception& e)
    {
        // Gestion de l'exception
        std::cout << "Erreur :  la matrice n'est pas inversible" << std::endl;
    }

    return 0;
}


[leon1789]

Question : pour déduire le theta as-ton besoin du sin et du cos ou seulement de l'un des deux

Si tu connais le signe de t , alors le cos suffit à la déterminer.
On peut imposer t compris entre 0 et Pi , quitte à prendre a < 0 .

[sylvain231]

ok merci j'attends donc votre autre méthode car il m'arrivera dans 50% des cas d'avoir des ouvertures quasi-plates

[sylvain231]

mais je pense que pour gérer à la fois les paraboles "serrées" et les paraboles "plates" il faut séparer en entrée les points "à gauche" du sommet des points "à droite" et ça c'est tout à fait possible car dans mon programme je les connaitrais

[leon1789]

mais je pense que pour gérer à la fois les paraboles "serrées" et les paraboles "plates" il faut séparer en entrée les points "à gauche" du sommet des points "à droite" et ça c'est tout à fait possible car dans mon programme je les connaitrais

Il n'y a pas lieu de scinder les points :
D'une part, à quoi cela pourrait servir concrètement ? on va obtenir deux paraboles qu'on aura plus ou moins de mal à recoller en une seule.
D'autre part, la méthode par rotation et affinité marche bien pour tes données, que la parabole soit plate ou serrée ou entre les deux, qu'il ait davantage de points d'un coté que de l'autre.

Pour la liste de 57 points au dessus, on trouve t = 0.524 et a' = 0.955 / 10^3

Et pour tracer la parabole simplement, il suffit de poser A = sin(t) / a' et B = cos(t) / a'
on a alors la parabole sous forme paramétrée avec r parcourant les nombres réels, c'est très pratique pour dessiner son graphe.

Finalement, c'est justifié mathématiquement, et c'est simple mis à part l'obtention de t qui demande une recherche numérique, tout le reste est direct.

X (r) = 524*r^2 + 906*r + 200
Y (r) = 906*r^2 - 524*r + 300
r parcourant l'intervalle [-0.3 ; 0.3 ]



Je crois bien qu'on a répondu à toutes tes questions, non ?


PS. Notre clown proclame que tout ça est du bricolage sans intérêt. Il va sans dire que lui ne produit aucun résultat, comme d'habitude. En réalité, il a des idées loufoques et ne comprend évidemment rien aux maths en général : ici, il ne sait ni obtenir une parabole à coup sûr, ni traiter l'information du sommet fixé (dixit). Trop fort.

[sylvain231]

je ne trouve pas comme toi et ma parabole ne colle pas dans les 2 cas

[sylvain231]

donc si je comprend bien tu as repris l'ancienne méthode ? mais moi elle ne marchait ni dans la première liste de points ni dans la deuxième, j'ai dû manquer qqch, peux-tu mettre les résultats intermédiaires que je vérifie mes calculs STP ?

[sylvain231]

voici mon code y-a-t-til des erreurs (vérifier notamment le code de rotation de points) :

Code: Tout sélectionner

std::vector<Point2d> rotatePoint2ds(const std::vector<Point2d>& points, double theta) {
    std::vector<Point2d> rotatedPoint2ds;
    for (const auto& point : points) {
        rotatedPoint2ds.push_back({ point.x * cos(theta) - point.y * sin(theta), point.x * sin(theta) + point.y * cos(theta) });
    }
    return rotatedPoint2ds;
}


std::vector<Point2d> translatePoint2ds(const std::vector<Point2d>& points, const Point2d& sommet) {
    std::vector<Point2d> translatedPoint2ds;
    for (const auto& point : points) {
        translatedPoint2ds.push_back({point.x - sommet.x, point.y - sommet.y});
    }
    return translatedPoint2ds;
}

void Parabole::fit2(std::vector<cv::Point2d> points, cv::Point2d sommet) {
    this->sommet = sommet;
    int n = points.size();
    cout << "sommet = " << sommet << endl;
    points = translatePoint2ds(points, sommet);
    double bestTheta = 0;
    double minimum = std::numeric_limits<double>::infinity();
    for (double theta = 0; theta <= CV_PI; theta += CV_PI / 1000) {

        vector<Point2d> rotated = rotatePoint2ds(points, theta);
        double E1 = 0;
        for (int i = 0; i < n; i++) {
            E1 += pow(rotated[i].y, 2);
        }
        double E2 = 0;
        for (int i = 0; i < n; i++) {
            E2 += pow(rotated[i].x, 2) * rotated[i].y;
        }
        double E3 = 0;
        for (int i = 0; i < n; i++) {
            E3 += pow(rotated[i].x, 4);
        }
        double delta = E1 - (E2 / E3);
        if (delta < minimum) {
            minimum = delta;
            bestTheta = theta;
        }
    }
    this->theta_radians = bestTheta;
    this->theta = bestTheta * 180 / CV_PI;
    vector<Point2d> rotated = rotatePoint2ds(points, bestTheta);
    double E2 = 0;
    for (int i = 0; i < n; i++) {
        E2 += pow(rotated[i].x, 2) * rotated[i].y;
    }
    double E3 = 0;
    for (int i = 0; i < n; i++) {
        E3 += pow(rotated[i].x, 4);
    }
    this->a = E2 / E3;
    this->aa = sin(theta_radians) / (4 * a);
    this->b = cos(theta_radians) / (4 * a);
    cout << "aa" << aa << "b" << b << endl;
    cout << "a=" << this->a << endl;

}

[leon1789]

donc si je comprend bien tu as repris l'ancienne méthode ? mais moi elle ne marchait ni dans la première liste de points ni dans la deuxième, j'ai dû manquer qqch, peux-tu mettre les résultats intermédiaires que je vérifie mes calculs STP ?

ok, on va y aller pas à pas. Il n'y a pas de raison que tu n'y arrives pas. Je reviens...
Et on va modifier un truc au passage.

[sylvain231]

merci, j'ai déjà vu un bug j'avais mis double delta = E1 - E2/ E3 au lieu de double delta = E1 - (pow(E2,2) / E3); mais ça ne change rien, ça ne marche toujours pas

[sylvain231]

fais le sur le premier exemple alors où la différence avec ce qu'on devrait trouver est beaucoup plus importante

[leon1789]

fais le sur le premier exemple alors où la différence avec ce qu'on devrait trouver est beaucoup plus importante

on va reprendre l'exemple avec 3 points (10;1) , (18;10) , (2;10) et le sommet en (0;0)

En modifiant les calculs.

On va chercher t pour minimiser E(Xt^4) - E(Xt^2 * Yt)^2 / E(Yt^2)
et poser c = E(Xt^2 * Y) / E(Y^2)
puis A = c * sin(t) et B = c * cos(t)

Une équation cartésienne de la parabole sera alors
(B * X - A *Y )^2 / (A^2 + B^2) - (A*X + B*Y) = 0
et une équation paramétrique (pour le dessin)
X (r) = A*r^2 + B*r
Y (r) = B*r^2 - A*r

Etape 1 : trouver t et c

mets à jour ton code avec ceci (tu mets l'indentation que tu veux, pas souci évidement)

Code: Tout sélectionner

(...)
        vector<Point2d> rotated = rotatePoint2ds(points, theta);
        double E1 = 0;
        for (int i = 0; i < n; i++) {
            E1 += pow(rotated[i].x, 4);
        }
        double E2 = 0;
        for (int i = 0; i < n; i++) {
            E2 += pow(rotated[i].x, 2) * rotated[i].y;
        }

        double E3 = 0;
        for (int i = 0; i < n; i++) {
            E3 += pow(rotated[i].y, 2);
        }

        double delta = E1 - (pow(E2,2) / E3);
        if (delta < minimum) {
            minimum = delta;
            bestTheta = theta;
        }
}

    this->theta_radians = bestTheta;
    this->theta = bestTheta * 180 / CV_PI;
    vector<Point2d> rotated = rotatePoint2ds(points, bestTheta);
    double E2 = 0;
    for (int i = 0; i < n; i++) {
        E2 += pow(rotated[i].x, 2) * rotated[i].y;
    }
    double E3 = 0;
    for (int i = 0; i < n; i++) {
        E3 += pow(rotated[i].y, 2);
    }
    this->c = E2 / E3;
    this->a = sin(theta_radians) *c;
    this->b = cos(theta_radians)  *c;
    cout << "a" << a << "b" << b << endl;


Et maintenant donne moi les t, c, a, b obtenus
avec les 3 points (10;1) , (18;10) , (2;10) et le sommet toujours en (0;0)