Discussion:Problème du cercle minimum

Titre

J'aurais intitulé l'article Problème du cercle minimal (au lieu de minimum ; Larousse : minimum en tant qu'adjectif est contesté ; Le Petit Robert : minimal tend à remplacer l'adjectif minimum) ou, plus précisément mais plus long, Problème du cercle de rayon minimal ou Problème du disque d'aire minimale ou, tout simplement, Problème du plus petit cercle/disque. JChG (d) 4 janvier 2013 à 16:39 (CET

J'ai repris un titre que j'avais trouvé quelque part. En cherchant à nouveau, il semble que ce soit la traduction faite en 1885, cf. [1]. Après, on peut expliciter dans l'intro.

cdang | m'écrire 4 janvier 2013 à 17:49 (CET)

On pourrait aussi éliminer le mot Problème, ce qui donnerait Cercle minimal, Cercle de rayon minimal, Disque d'aire minimale, ... JChG (d) 7 janvier 2013 à 14:13 (CET)

Notons qu'il existe aussi un Problème du cercle — faudrait-il donc éliminer le mote « problème » et renommer l'article en Cercle ? ;-)

Par ailleurs, ce n'est pas le cercle minimum en soi qui est intéressant, c'est la manière dont on le trouve, donc la résolution du problème. On peut aussi noter que l'article anglais contient le terme « Problem », et qu'il s'agit d'un problème d'optimisation, dans la même lignée que le problème du voyageur de commerce, dont le titre contient lui aussi le mot « problème ».

Ceci dit, je n'ai aucun attachement émotionnel au présent titre, renommez l'article si ça vous chante, moi, j'ai fait mon taf.

cdang | m'écrire 8 janvier 2013 à 12:08 (CET)

Ce n'est bien sûr pas le mot problème qui fait difficulté. Le but recherché est d'avoir un titre court, parce que les titres courts ont plus d'impact, tout en gardant de la précision. Ici le mot problème me semble superflu, ce qui n'est pas le cas des autres titres que vous citez. JChG (d) 12 janvier 2013 à 09:10 (CET)

Rapport avec la Régression circulaire

Ce message a été importé de Discussion:Régression circulaire#problème de Sylvester par cdang | m'écrire 7 janvier 2013 à 09:40 (CET)

Le seul rapport que je trouve est le suivant:

Soient n points dans R^2. Soit C(a,r) la solution du problème de régression circulaire(= le cercle le plus proche des n points, a= centre et r=rayon) et soit e la distance la plus grande entre les points extérieurs à C et le cercle C.

Alors la solution du problème de Sylvester , le cercle C'(b, r') , est telle que r' <= r+e. Autrement dit C(a, r+e) est une solution approchée du problème de Sylvester.Cordialement. Claudeh5 (d) 6 janvier 2013 à 10:12 (CET)

Programme Scilab

Bonjour,

voici un programme Scilab qui trouve le cercle minimum en utilisant la méthode de Chrystal, la recherche de l'enveloppe convexe se faisant avec la marche de Jarvis. Ce n'est pas la méthode la plus optimale, mais elle est simple à coder. Si n est le nombre de points des données et m le nombre de points de l'enveloppe convexe (m ≤ n), alors la complexité en temps de la marche de Jarvis est en O(m⋅n) et la complexité de l'algorithme de Chrystal est en O(m²), soit une complexité globale en O(m⋅n).

cdang | m'écrire 23 janvier 2017 à 15:05 (CET)

Programme Scilab chrystal_jarvis.sce

//============================================================================
// nom : chrystal_jarvis.sce
// auteur : Christophe Dang Ngoc Chan
// date de création : 2017-01-23
// dates de modification :
//    2017-01-24 : clarification de la fonction marche_jarvis()
//    2017-01-30 : calcul de la taille dans la fonction marche_jarvis()
//    2017-02-07 : détermination du cercle passant par trois points dans une
//        fonction séparée
//
//----------------------------------------------------------------------------
// version de Scilab : 6.0.0
// module Atoms requis : aucun
//----------------------------------------------------------------------------
// Objectif : Déterminer l'enveloppe convexe d'un ensemble de points par la
//    méthode de la marche de Jarvis, puis déterminer le cercle minimum par la
//    méthode de Chrystal.
// Entrées : aucune (points générés aléatoirement)
// Sorties : fenêtre graphique avec le résultat et affichage de la durée de
//    traitement
//============================================================================


// ******************
// * Initialisation *
// ******************


// **************
// * Constantes *
// **************

pisurdeux = 0.5*%pi; // évite de le calculer plusieurs fois

// ********** Génération aléatoire de points **********

n = 20; // nombre de points

// Pour une surface carrée
xmin = 0; xmax = 10; // bornes
ymin = 0; ymax = 10;
largeurx = xmax - xmin;
largeury = ymax - ymin;
// Commenter les deux lignes suivantes pour avoir un disque
points = rand(n, 2, "uniform").*[ones(n, 1)*largeurx, ones(n, 1)*largeury]...
   + [ones(n, 1)*xmin, ones(n, 1)*ymin];

// Pour un disque
r = 10; // rayon
x0 = 5; y0 = 5; // centre du cercle
rayons = r*rand(n, 1, "uniform");
angles = 2*%pi*rand(n, 1, "uniform");
// points = [rayons, rayons].*[cos(angles), sin(angles)]; // décommenter
// cette ligne pour avoir un disque


// *************
// * Fonctions *
// *************

// ********** marche_jarvis **********

function [ind_env_sens] = marche_jarvis_sens_unique(pts, sens)
    // détermine la moitié de l'enveloppe convexe
    //     par l'algorithme de la marche de Jarvis
    // entrées : pts : matrice n × 2 de réels, ensemble de points [X, Y]
    //    n : entier, nombre de points (évite d'avoir à le recalculer)
    // sens : entier, +1 ou -1 selon le sens de parcours
    // sorties : ind_env_sens : vecteur d'entiers,
    //    indices des points de l'enveloppe convexe
    n = size(pts, "r"); // nombre de points
    if (sens == 1)
        then indice_courant = 1;
            indice_stop = n;
            ind_env_sens = [1];
        else indice_courant = n;
            indice_stop = 1
            ind_env_sens = [n];
    end
    taille_env = 1; // nombre de points dans l'enveloppe
    condition = %t;
    while condition
        taille_env = taille_env + 1;
        pts_foo = pts - ones(n, 1)*pts(indice_courant, :);
        pts_foo(indice_courant, :) = [];
        angle_polaire = atan(-pts_foo(:, 1)*sens, pts_foo(:, 2)*sens);
        [angle_extr, indice_extr] = min(angle_polaire);
        if indice_extr >= indice_courant
            then indice_extr = indice_extr + 1;
        end
        condition = (indice_extr <> indice_stop);
        ind_env_sens(taille_env) = indice_extr;
        indice_courant = indice_extr;
    end
endfunction

function [env_cvx] = marche_jarvis(pts)
    // détermine l'enveloppe convexe
    //     par l'algorithme de la marche de Jarvis
    // entrées : pts : matrice n × 2 de réels, ensemble de points [X, Y]
    // sorties : env_cvx : matrice ? × 2 de réels, points de l'enveloppe convexe
    
    n = size(pts, "r"); // nombre de points
    pts = gsort(pts, "lr", "i"); // tri
    
    // progression de xmin vers xmax
    ind_env_positif = marche_jarvis_sens_unique(pts, +1)
    
    // progression de xmax vers xmin
    ind_env_negatif = marche_jarvis_sens_unique(pts, -1)
    
    // Bilan
    ind_env = [ind_env_positif', ind_env_negatif']
    env_cvx = pts(ind_env, :);
endfunction


// ********** Géométrie : cercle passant par trois points **********

function [C, r] = cercle_par_trois_points(P1, P2, P3, l12, l23, l31, ...
    angle1, angle2, angle3)
    // Détermine le cercle passant par trois points non-alignés.
    // Entrées : P1, P2 et P3 : points, représentés par des vecteurs 1 × 2 de
    //         réels, coordonnées (x, y).
    //     l12, l23, l31 : longueurs des côtés du triangle, réels
    //     angle1, angle2, angle3 : angles aux sommets du triangle, réels
    // Sorties : C : centre du cercle, vecteurs 1 × 2 de réels, coordonnées
    //         (x, y) ;
    //     r : réel, rayon du cercle.
    
    // vecteurs des côtés
    r = 0.5*mean([l23/sin(angle1), l31/sin(angle2), l12/sin(angle3)]);
    // moyenne de 3 valeurs normalement égales pour minimiser les erreurs
    invD = 0.5*inv(P1(1)*(P2(2) - P3(2)) + ...
        P2(1)*(P3(2) - P1(2)) + ...
        P3(1)*(P1(2) - P2(2)));
    C = invD*[norm(P1) norm(P2) norm(P3)].^2* ...
        [P2(2) - P3(2), P3(1) - P2(1)
        P3(2) - P1(2), P1(1) - P3(1)
        P1(2) - P2(2), P2(1) - P1(1)];
endfunction


// ********** chrystal **********

function [C, r] = chrystal(poly_cvx)
    // détermine le cercle minimum d'un polygone convexe
    //     par l'algorithme de Chrystal
    // entrées : poly_cvx : matrice n × 2 de réels, ensemble de points [X, Y]
    //    classés dans le sens trigo ou anti-trigo
    //    (on utilise ici le fait que les deux premiers points sont contigus).
    // sorties : C : matrice 2 × 2 de réels, coordonnées du centre du cercle ;
    //    r : réel, rayon du cercle
    n = size(poly_cvx, "r"); // nombre de points
    indices = 1:n;
    i1 = 1; // premier sommet
    i2 = 2; // deuxième sommet
    condition = %t;
    
    while condition
        point1 = poly_cvx(i1, :);
        point2 = poly_cvx(i2, :);
        P12 = ones(n-2, 1)*(point2 - point1); // vecteur du côté [P1P2]
        nP12 = sqrt(P12(:, 1).^2 + P12(:, 2).^2); // longueur du côté P1P2
        ind_loc = indices;
        ind_loc([i1, i2]) = [];
        point3 = poly_cvx(ind_loc, :); // 3e sommet
        P23 = point3 - ones(n-2, 1)*point2; // vecteurs des autres côtés
        nP23 = sqrt(P23(:, 1).^2 + P23(:, 2).^2);
        P31 = ones(n-2, 1)*point1 - point3;
        nP31 = sqrt(P31(:, 1).^2 + P31(:, 2).^2);
        invnnn = 1./(nP12.*nP23.*nP31); // optimise les calculs
        // angle à partir du produit scalaire
        angle1 = acos(-(nP23.*invnnn).*sum(P12.*P31, "c"));
        angle2 = acos(-(nP31.*invnnn).*sum(P12.*P23, "c"));
        angle3 = acos(-(nP12.*invnnn).*sum(P23.*P31, "c"));
        [minangle3, k3] = min(angle3); // angle le plus petit en P3
        minangle1 = angle1(k3);
        minangle2 = angle2(k3);
        minpoint3 = point3(k3, :);
        // opposé à [P1P2] : [P1P2] est un diamètre
        if (minangle3 >= pisurdeux) then // triangle obtus en P3
            // [P1P2] est un diamètre
            C = mean([point1 ; point2], "r");
            r = 0.5*nP12(1);
            condition = %f;
        elseif and([minangle1 minangle2 minangle3] < pisurdeux) then
            [C, r] = cercle_par_trois_points(point1, point2, minpoint3, ...
               nP12(k3), nP23(k3), nP31(k3), minangle1, minangle2, minangle3);
            condition = %f;
        elseif minangle1 > pisurdeux then // triangle obtus en P1
            i1 = i2;
            i2 = ind_loc(k3);
        else // triangle obtus en P2
            i2 = ind_loc(k3);
        end
    end
endfunction


// ********** cercle_mini **********

function [C, r] = cercle_mini(points)
    // détermine le cercle minimum d'un polygone convexe dans les cas dégénérés
    //     et appelle l'algorithme de Chrystal dans le cas général
    // entrées : points : matrice n × 2 de réels, ensemble de points [X, Y]
    // sorties : C : matrice 1 × 2 de réels, coordonnée sdu centre du cercle ;
    //    r : réel, rayon du cercle
    n = size(points, "r"); // nombre de points
    select n
    case 1 then // un seul point
        C = points;
        r = 0;
    case 2 then // deux points
        C = mean(points, "r");
        r = 0.5*norm(points(2, :) - points(1, :));
    else
        enveloppe_convexe = marche_jarvis(points);
        [C, r] = chrystal(enveloppe_convexe)
    end
endfunction


// **********************
// * Programe principal *
// **********************

// ********** Détermination du cercle minimum **********

tic();
[C, r] = cercle_mini(points);
duree = toc();

// Si l'on veut tracer l'enveloppe convexe,
// remplacer les lignes ci-dessus par :
// enveloppe = marche_jarvis(points); 
// [C, r] = chrystal(enveloppe);


// ********** Affichage du résultat **********

figure(0);
clf();
plot(points(:, 1), points(:, 2), "+"); // points des données

// Si l'on veut tracer l'enveloppe convexe, ajouter la ligne suivante
// plot([enveloppe(:, 1) ; enveloppe(1, 1)], [enveloppe(:, 2) ; enveloppe(1, 2)]);

plot(C(1), C(2), "or"); // centre du cercle

anglescercle = linspace(0, pisurdeux, 20);
cosinus = cos(anglescercle);
sinus = sin(anglescercle);
plot(C(1) + r*[cosinus, -sinus, -cosinus, sinus], ...
    C(2) + r*[sinus, cosinus, -sinus, -cosinus], "r"); // cercle

axe = gca();
axe.isoview="on"; // axes isométriques

disp("Durée du calcul : "+string(duree)+" s. pour "+string(n)+" points");