Régression parabolique non triviale

[leon1789]

Si on place le sommet de la parabole pas trop loin, en (191 , 291) par exemple, alors la dernière méthode de Ben fonctionne bien.

La méthode ici :
superieur/regression-parabolique-non-triviale-t278846-240.html#p1549234

On obtient cette équation de parabole :

[leon1789]

Graphe de cette parabole serrée


Et la méthode de Ben fonctionne aussi pour la parabole plate, comme on l'a vue hier soir.

[sylvain231]

ah bonne nouvelle je teste

[sylvain231]

non moi ça ne marche pas 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]

et moi ça ne marche pas non plus pour la parabole plate

[sylvain231]

j'ai dû faire une erreur qq part mais je vois pas où

[leon1789]

Tourne d'abord les points d'un angle theta

puis calcule
alpha = E40
beta = E21
gamma = E02

puis cherche t entre 0 et Pi pour minimiser alpha + gamma - sqrt((alpha- gamma)^2+4 *beta^2) :

Ensuite, pour cet angle t minimisant, pose
p := 1/2/beta *(-alpha + gamma + sqrt((alpha- gamma)^2+4 *beta^2)) ;

et la parabole est
0 = p*( cos(t) * (x-S[1]) - sin(t) * (y-S[2]) ) ^2 - ( sin(t) * (x-S[1]) + cos(t) * (y-S[2])) ) ;

[sylvain231]

c'est pas les formules de Ben ?!

[sylvain231]

non ça ne marche pas non plus voici mon code :

Code: Tout sélectionner

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) {
       
        points = rotatePoint2ds(points, theta);
        double alpha = E40;
        double beta = E21;
       
        double gamma = E02;
        double Mmin = 1 / 2 / beta * (-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)) * (gamma - alpha + sqrt(pow(alpha - gamma, 2) + 4 * pow(beta, 2)));
    this->theta = 180 * theta_radians / CV_PI;
}


[sylvain231]

non mon code n'était pas bon voici le code corrigé mais ça ne marche pas non plus :

Code: Tout sélectionner

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) {
       
        vector<Point2d> points2 = rotatePoint2ds(points, theta);
        double alpha = E(points2, 4, 0);
        double beta = E(points2, 2, 1);
       
        double gamma = E(points2, 0, 2);
        double Mmin = 1 / 2 / beta * (-alpha + gamma + sqrt(pow(alpha - gamma, 2) + 4 * pow(beta, 2)));
        if (Mmin < minimum) {
            minimum = Mmin;
            this->theta_radians = theta;
        }
    }

    vector<Point2d> points2 = rotatePoint2ds(points, this->theta_radians);
    double alpha = E(points2, 4, 0);
    double beta = E(points2, 2, 1);

    double gamma = E(points2, 0, 2);
    this->a = (1 / (2 * beta)) * (gamma - alpha + sqrt(pow(alpha - gamma, 2) + 4 * pow(beta, 2)));
    this->theta = 180 * theta_radians / CV_PI;
}


[sylvain231]

il y avait encore une erreur j'ai corrigé cette ligne : double Mmin = alpha + gamma - sqrt(pow(alpha - gamma,2) + 4 * pow(beta ,2)); mais ça ne marhce toujours pas

[leon1789]

minimiser alpha + gamma - sqrt((alpha- gamma)^2+4 *beta^2)

donc

Code: Tout sélectionner

        double Mmin = alpha + gamma -  sqrt(pow(alpha - gamma, 2) + 4 * pow(beta, 2)) ;

[sylvain231]

regarde le message précédent je l'ai corrigé mais ça ne marche toujours pas

[sylvain231]

je trouve :
theta = 50.04, theta_radians= 0.873363 a = -1.98229 sommet = [191, 291]

[sylvain231]

tu veux qu'on vérifie alpha, beta, gamma pour theta =0 ?

[sylvain231]

j'obtiens pour theta =0 :
alpha = 1.7046e+06 beta = -26655.9 gamma = 509.903

[leon1789]

je trouve :
theta = 50.04, theta_radians= 0.873363 a = -1.98229 sommet = [191, 291]

impeccable !

[sylvain231]

ben non c'est pas bon ça colle pas aux points ! ma parabole est 134 degrés, a =2

[sylvain231]

là elle est perpendiculaire et beaucoup trop serrée !

[sylvain231]

ou alors c'est mon code de dessin qui n'est pas bon je te le donne :

Code: Tout sélectionner

void Parabole::draw(Mat& img, cv::Vec3b color) {
    // Itérer sur chaque pixel
    for (double y = 0; y < img.rows; y += 0.1) {
        for (double x = 0; x < img.cols; x += 0.1) {
            // Calculer X' et Y'
            double X_prime = x - sommet.x;
            double Y_prime = y - sommet.y;

            // Calculer la valeur de l'équation de la parabole
            double val;
            val = a * pow(cos(this->theta_radians) * X_prime + sin(this->theta_radians) * Y_prime, 2)
                 - sin(this->theta_radians) * X_prime - cos(this->theta_radians) * Y_prime;
           
            //  }

              // Si val est proche de 0, colorer le pixel
            if (fabs(val) < 0.5) { // Vous pouvez ajuster la tolérance
                img.at<cv::Vec3b>(y, x) = color; // Par exemple, en blanc
            }
        }
    }
}