[Challenge]Quickhull

[(gile)]

Salut,

En tombant sur cet article à propos d'un algorithme pour calculer l'enveloppe convexe d'une collection de points, je me suis dit que ça pourrait faire l'objet d'un nouveau challenge et que ça pourrait aussi être une base pour résoudre un problème plus complex posé par @Luna dans ce sujet : faire une enveloppe concave (concave hull).

Donc, il s'agit toujours d'écrire un fonction (dans le langage de son choix) qui prend en argument une liste de points et qui renvoie la liste des points qui constituent l'enveloppe convexe (convex hull) des points passés en argument. Il existe au moins un autre algorithme bien connu pour extraire l'enveloppe convexe d'une collection de points, mais comme il a déjà été souvent implémenté, je trouve plus intéressant de se pencher sur celui-là et je suppute qu'il se prête mieux à la généralisation vers une enveloppe concave (mais c'est une autre histoire...).

On doit pouvoir tester la fonction 'quickhull' via une sélection d'entités points (nodal) en dessinant une polyligne.
Exemple en LISP :
 

(defun c:test (/ ss i pts)
  (if (setq ss (ssget '((0 . "point"))))
    (progn
      (repeat (setq i (sslength ss))
	(setq pts
	       (cons (cdr (assoc 10 (entget (ssname ss (setq i (1- i))))))
		       pts
	       )
	)
      )
      (command "_pline")
      (foreach p (quickhull pts)
	(command "_non" p)
      )
      (command "_close")
    )
  )
  (princ)
)

 

[Fraid]

Bonjour,

Si il y a un argument en plus qui est la "résolution" (longueur max des faces du polygone) je sais faire.

Mais sans, je ne vois pas comment déterminer si on doit rejoindre ou non un point...surtout si on a une volée en forme de croissant.

Notamment avec l'algorithme que tu nous présente au début du message. 

Avant que tu poste ce challenge, je commençais à reprendre mon lisp qui me donne un profil dans un nuage de points.

(l'outil Autocad est vraiment inutile...) en essayant d'intégrer l'algo que tu avais posté sur theswamp.

du coup, il me semble que c'est la solution non?

[Curlygoth]

en 2D je saurais le faire en testant si le point est dans la polyligne ou non si en parle de 3D je passe mon tour XD

(j'utilise la technique du nombre d'intersection d'une demi droite a partir d'un point si le nombre de point est impair on est dedans sinon on est dehors ^^)

Je donne la source pour les curieux :

http://maxence.delannoy.pagesperso-orange.fr/pt_poly.htm

En 3D avec un nuage de points il faudrait relié les points en eux pour creer des face 3D et ensuite meme systeme tester les point a l'interieur et "compter le nombre d'intersection"

Quote

Si il y a un argument en plus qui est la "résolution" (longueur max des faces du polygone) je sais faire.

je suis d'accord ! moi je dis : je pense que se serait réalisable pour moi ^^ en VBA mais je le ferais pas demain XD

[Fraid]

21 minutes ago, Curlygoth said:

3D

non non, j'extrait les points 3d que je projette à plat ensuite

[Fraid]

je suis désolé, mais avec ce lisp pour extraire des profils d'un nuage, j’atteins les limites de partage...

déjà parce que ce n'est qu'une expérimentation encore en cours

et que cela répond à un cahiers des charges bien précis qui ferais gagner trop d'argent à des gent comme Vinchi, ou bouigche pour pas les citer...

[Curlygoth]

je demande rien ^^ juste des pistes ! moi ça me va !

mais du coup tu prends des deltas Z pour checker et voir si ton point est dedans ou non ^^ enfin c'est ce que je comprends dans ton message ^^

moi je pensais relié les point pour faire une triangulation de tous les points pour faire un mappage et tester le point via la solution que j'ai mis plus haut en calculant le nombre dintersection entre les face 3D (en union pour pas faire autant de boucle que de faces XD) et compter le nombre de point d'intersection

 

[(gile)]

Le but du challenge est de fournir une fonction qui prends en argument une collection de points 2D (je conseille des commencer avec une ou quelques dizaines de points pas plus) et qui renvoie les points qui constituent l'enveloppe convexe de cette collection.
J'ai implémenté ce type de fonction en utilisant l'algorithme dit du parcours de Graham (en C# et F# ainsi qu'en LISP : gc:ConvHull dans MathGeom sur cette page).
Je proposais de le faire avec ici avec QuickHull, mais après tout, toute solution est bienvenue.

[Curlygoth]

a donc si je comprends bien on parle d'un nuage de point dans le plan et il faut définir quels sont les points qui délimite l'extérieur ?

pour ensuite les comparer a d'autres ?

MMmm... donc je prends les 3 premiers points dasn le nuages et je tests s'il y a des points plus excentré ?

C'est faisable en VBA... mais à ma connaissance je vais partir de zéro ^^ il n'y a pas de date limite ? xD

[(gile)]

@CurlygothSi tu lis attentivement les articles de Wikipédia  que j'ai donné en liens, tu auras une description du problème dans celui qui concerne l'enveloppe convexe et la description de deux algorithmes : parcours de Graham et quickhull.

[Fraid]

Bonjour,

En attendant, un p'tit lispounet pour dessiner des points dans une zone de manière pseudo-aléatoire

(defun add-random-pts (nb zon / acdc cpt nbc)
        (setq acdc (vla-get-activedocument (vlax-get-acad-object)) cpt -1 nbc (strlen (itoa (1- (caadr zon)))))
        (while (< (setq cpt (1+ cpt)) nb)
            (vla-AddPoint (vla-get-modelspace acdc) (vlax-3d-point (+(random nbc nil)(caar zon)) (+ (random nbc nil)(cadar zon)) 0))
        )
)

;Random
;Retourne un Entier pseudo-aléatoire 
;(random arg1 arg2)
;Arg1 : Entier nombre de chiffres maxi. ex 2 -> 1 à 99, 5 -> 1 à 99999
;Arg2 : T ou nil si T 2 -> 10 à 99, 5 -> 10000 à 99999
;Retour : Entier
(defun Random (n m / rd)
    (if (null *RDM*)(setq *RDM* (getvar 'date)))
    (setq rd (fix (* (/ (setq *RDM* 
                            (rem (+ 1013904223.0   (*  1664525.0 *RDM*)) 
                                 4294967296.0 
                            )
                        )
                        4294967296.0
                    )
                    (expt 10.0 n)
                 )
            )
    )
    (if (and m (< rd (expt 10 (1- n))))
        (setq rd (+ rd (expt 10 (1- n))));Maintient
    )
    rd
)
;pour 50 points dans une zone 0,0 100,100
(add-random-pts 50 '((0 0)(100 100)))             

 

[Curlygoth]

je pensais réfléchir par moi même mais effectivement s'il y a des exemples ça ira plus vite XD

[Fraid]

T'inquiète, c'est juste pour éviter de dessiner une volée

cela ne résous en rien le challenge

[Curlygoth]

Verifier si un point est dans une polyligne :

Function Point_dans_Polyligne(PolyLigne, PT_TEST) As Boolean
Dim Pt1(0 To 2) As Double
Dim Pt2(0 To 2) As Double
Dim Droite As AcadRay
Dim Points As Variant
Pt2(0) = PT_TEST(0) + 1
Pt2(1) = PT_TEST(1)
Pt2(2) = PT_TEST(2)
Set Droite = ThisDrawing.ModelSpace.AddRay(PT_TEST, Pt2)
Droite.Update
PointsINT = PolyLigne.IntersectWith(Droite, acExtendNone)
Droite.Delete
Dim NbrePoints As Integer
NbrePoints = UBound(PointsINT) / 3
Point_dans_Polyligne = False
If NbrePoints Mod 2 = 0 Then
Point_dans_Polyligne = False
Else
Point_dans_Polyligne = True
End If
End Function

[(gile)]

Outre que je ne vois pas bien le rapport avec ce challenge, cet algorithme n'est pas fiable à 100%.
Tu n'envisages pas le cas où la droite passe par un sommet (point vert) ou deux sommets (point jaune), ni celui où elle est tangente à un arc (point rouge). Dans ces trois cas, les points sont à l'intérieur mais la droite a 3 intersections...

 

Ça m'a toujours intrigué cette tendance à la verbosité dans la culture VB :

Point_dans_Polyligne = False
If NbrePoints Mod 2 = 0 Then
Point_dans_Polyligne = False
Else
Point_dans_Polyligne = True
End If

Alors qu'on devrait pouvoir écrire tout simplement :

Point_dans_Polyligne = NbrePoints Mod 2 <> 0

 

[Curlygoth]

Wow ! c'est vrai tu as totalement raison ! (je suis en demi droite)

donc plus partie droite de ta ligne rouge et de ta ligne verte et partie gauche de ta ligne jaune !

mais oui je vois bien la problématique ! et bien écoute je te remercie de ton retour je vais trouver effectivement une alternative à cette problématique...

Vérifier si un point d'intersection est aussi un point de la polyligne ? et si c'est le cas je l'ignore .. attend dans le cas de la ligne jaune ça marche pas XD

Allez je vais reflechir à ça 😉

 

[VDH-Bruno]

Bonjour,

Rapidement un premier jet:

(defun Quickhull (lst / P Q position cut-pts getPtH HalfConHull)

  ;; Position d'un point C par rapport à une droite orienté (AB): 0 = sur la droite; 0 < au-dessus; 0 > au-dessous
  (defun position (A B C)
    (- (* (- (car B) (car A)) (- (cadr C) (cadr A)))
       (* (- (cadr B) (cadr A)) (- (car C) (car A)))
    )
  )

  ;; coupe une lite de points en 2 par rapport à un segment ->( (liste des points dessus)(liste des points dessous))
  (defun cut-pts (A B lst haut bas)
    (cond ((null lst) (list haut bas))
	  ((minusp (position A B (car lst)))
	   (cut-pts A B (cdr lst) haut (cons (car lst) bas))
	  )
	  ((zerop (position A B (car lst)))
	   (cut-pts A B (cdr lst) haut bas)
	  )
	  (T (cut-pts A B (cdr lst) (cons (car lst) haut) bas))
    )
  )

  ;; Recherche de point le plus éloigné en fonction de 2 pt et d'une liste de points
  (defun getPtH	(pt1 pt2 lpt / dist H)
    (setq dist 0.)
    (foreach x lpt
      (if (< dist (+ (distance pt1 x) (distance x pt2)))
	(setq H	   x
	      dist (+ (distance pt1 x) (distance x pt2))
	)
      )
    )
    H
  )

  ;; Demi enveloppe convexe
  (defun HalfConHull (P Q lst / H)
    (cond
      ((null (cdr lst)) lst)
      (T
       (setq H	 (getPtH P Q lst)
	     lst (cut-pts P H lst nil nil)
       )
       (append
	 (HalfConHull P H (car lst))
	 (cons H
	       (HalfConHull H Q (car (cut-pts H Q (cadr lst) nil nil)))
	 )
       )
      )
    )
  )

  ;; Traitement
  (if (< (length lst) 4)
    lst
    (setq P   (assoc (apply 'min (mapcar 'car lst)) lst)
	  lst (vl-remove P lst)
	  Q   (assoc (apply 'max (mapcar 'car lst)) lst)
	  lst (cut-pts P Q lst nil nil)
	  lst (append (cons P (HalfConHull P Q (car lst)))
		      (cons Q (HalfConHull Q P (cadr lst)))
	      )
    )
  )

)

Pas tout à fait la version que je voulais donner, mais celle-ci à l'air de faire le travail. Peut être si je peux (le temps) je regarderai pour d'éventuelle optimisation.

(Ps: J'ai pris le partie de supprimer les éventuels points aligné dans les segments droit😉)

A+ Bruno

[didier]

Bonjour

Ah zut de zut

Quelqu'un a répondu !
Pour une fois que je me mettais dans un challenge que j'avais compris !
(car les listes triées deux à deux dans l'ordre alphabétique hébreu médiéval avec quelques chiffres, très peu pour moi, ça ne me parle guère)
Allez j'attends encore un peu pour peaufiner ma version et je dépose.

Je pense que pour les prochains il faut donner une date de dépôt,
qu'on le veuille ou non on lit les réponses,
dans le cas présent j'ai sauté le message de @VDH-Bruno
mais il existe et ça devrait pas pour l'égalité des chances.

Amicalement

[VDH-Bruno]

Désolé Didier, mais comme je voyais pas mal de commentaires, j'ai cru à tord que le challenge était un peu bloqué... Alors je m'y suis attelé ce matin. 

Amicalement aussi

A+ Bruno

 

[didier]

Bonjour @VDH-Bruno

Il n'y a rien de personnel contre toi, c'est toi qui publies, c'est tout, du coup je te cite.
Je pense vraiment que pour ces challenges il faudrait définir une date de dépose à respecter par chacun.
Ce qui n'empêche pas les discussions, questions, clarifications, on est bien d'accord.

Avant de déposer ma solution (que j'ai du mal à la finaliser) j'attends un peu, car par manque de temps avec les autres occupations de la vie courante.

Amicalement

[didier]

Bonjour

En parlant de questions, j'en ai une :

Sur l'image suivant le résultat à obtenir est bien la solution ROUGE ?
La ligne CYAN dessinée à la main n'est pas à faire entrer dans le périmètre, c'est bien ça ?

Amicalement

[Fraid]

Bonjour,

Je me permet de poster

mais c'est tellement brut, que cela répond au challenge que part son résultat.

pas son algorithme "Viking"

(defun Quickhull (lpt / acdc mods oldosmode oldcmdecho lptmp ptrec zone sel ob ct)
    (vla-zoomextents (vlax-get-acad-object))
    (setq acdc (vla-get-activedocument (vlax-get-acad-object)) 
          mods (vla-get-modelspace acdc)
          oldosmode (getvar 'osmode)
          oldcmdecho (getvar 'cmdecho)        
          zone (coordaffichage)
          lptmp (cdr lpt)
    )
    (setvar 'osmode 0)
    (setvar 'cmdecho 0)
    (foreach pt1 lpt
        (foreach pt2 lptmp
            (vla-AddLine mods (vlax-3d-point (car pt1)(cadr pt1) 0) (vlax-3d-point (car pt2)(cadr pt2) 0))
        )
        (setq lptmp (cdr lptmp))
    )
    (setq ptrec (vlax-make-safearray vlax-vbDouble '(0 . 14)))
    (vlax-safearray-fill ptrec  (list   (1-(caar zone)) (1-(cadar zone)) 0
                                        (1-(caar zone)) (1+(cadadr zone)) 0 
                                        (1+(caadr zone)) (1+(cadadr zone)) 0
                                        (1+(caadr zone)) (1-(cadar zone)) 0
                                        (1-(caar zone)) (1-(cadar zone)) 0
                                )
    )
    (vla-AddPolyline mods ptrec)
    (vla-zoomextents (vlax-get-acad-object))
    (vla-regen acdc acActiveViewport)
    (vl-cmdf "_.-boundary" (car zone) "")
    (setq sel (ssget "_x" (list (cons 0 "*POLYLINE"))) ct -1)
    (while (setq ob (ssname sel (setq ct (1+ ct))))
        (if (/= ct 1)(vla-delete (vlax-ename->vla-object ob)))
    )
    (setq sel (ssget "_x" (list (cons 0 "LINE"))) ct -1)
    (while (setq ob (ssname sel (setq ct (1+ ct))))
        (vla-delete (vlax-ename->vla-object ob))
    )
    (setvar 'osmode oldosmode)
    (setvar 'cmdecho oldcmdecho)
    (princ)
)
(defun coordaffichage ( / viewctr viewsize screensize)  
  (mapcar 'set '(viewctr viewsize screensize) (mapcar 'getvar '(viewctr viewsize screensize)))
  ((lambda (dimxy) (list (mapcar '- viewctr dimxy) (mapcar '+ viewctr dimxy)))
    (list (* viewsize (/ (car screensize) (cadr screensize) 2.)) (/ viewsize 2) 0.0)
  )
)

et un test un peux différent, puisque je dessine la polyligne

(defun c:test (/ ss i pts)
  (if (setq ss (ssget '((0 . "point"))))
    (progn
      (repeat (setq i (sslength ss))
        (setq pts
               (cons (cdr (assoc 10 (entget (ssname ss (setq i (1- i))))))
                   pts
               )
        )
      )
      (quickhull pts)
    )
  )
  (princ)
)

 

[(gile)]

Super, ça démarre.

@didier C'est pas grave si des réponses commencent à arriver, ce n'est pas un concours il s'agit juste d'apprendre en comparant différentes façon de faire et il me semble qu'il reste de place pour bien d'autres réponses. J'attendrai la tienne pour poster ma réponse en LISP.
@VDH-Bruno Bravo, ça répond tout à fait à la demande. Je suis sûr que si tu as le temps tu pourras "optimiser" ta routine (ma réponse en LISP est proche de la tienne, mais plus concise).
@Fraid Ça ne répond pas tout à fait à la demande donc on ne peut pas comparer directement avec les autres réponses et on est assez loin de l'algorithme du quicksort. En tout cas, c'est pas vraiment 'quick'. J'ai eu un peu peur quand vous parliez de "nuage de points", avec ta solution, il y a intérêt à ce que ce soit un cirrus de beau temps.

Suivant la tâche à accomplir, certains langages sont plus ou moins bien adaptés. Dans ce cas là, parmi les langages que je connais F# est certainement le plus adapté. Je ne pense pas "divulgacher" le challenge en donnant une réponse en F# maintenant.

let quickhull (pts: Point2d list) =
    let area (p1: Point2d) (p2: Point2d) (p3: Point2d) =
        ((p2.X - p1.X) * (p3.Y - p1.Y)) - ((p2.Y - p1.Y) * (p3.X - p1.X))
    let rec rightPts p1 p2 min acc = function   
        | [] -> acc
        | x :: xs -> 
            match area p1 p2 x with
            | a when 0. <= a -> rightPts p1 p2 min acc xs
            | a -> 
                match acc with
                | [] -> rightPts p1 p2 a [x] xs
                | y :: ys -> 
                    if a < min
                    then rightPts p1 p2 a (x :: acc) xs
                    else rightPts p1 p2 min (y :: x :: ys) xs
    let rec loop p1 p2 pts =
        match rightPts p1 p2 0. [] pts with
        | [] -> [p1]
        | h :: t -> loop p1 h t @ loop h p2 t
    let minX = pts |> List.minBy (fun p -> p.X)
    let maxX = pts |> List.maxBy (fun p -> p.X)
    loop minX maxX pts @ loop maxX minX pts

 

[didier]

Bonjour @(gile)

Assurément le F# n'est pas vulgarisé... la lecture rend dubitatif ! (ce n'est pas un gros mot)

Je vais me remettre sur ma réponse, je ne travaille pas vite et je n’ai pas beaucoup de temps, tout cumulé je pense livrer (ce soir si les enfants me libèrent)
Demain matin on aura le fin mot.

Ce qui serait bien, c'est ce que je vais essayer de faire, ce serait de déposer des commentaires en entête de ligne pour expliquer ce qu'on fait
Histoire que ce soit un tantinet formateur pour les autres, non ?

Au niveau rapidité je pense être dans les clous, mais tu sais à quel point je suis impératif, je tente de me soigner.

Amicalement

[VDH-Bruno]

1 hour ago, (gile) said:

@VDH-Bruno Bravo, ça répond tout à fait à la demande. Je suis sûr que si tu as le temps tu pourras "optimiser" ta routine (ma réponse en LISP est proche de la tienne, mais plus concise).

@(gile)Je ne doute pas de la concision de ta version, que je détaillerai avec sourire en me disant "mais c'est bien sur", pour un second jet c'est compliqué, je suis entrain de pisser du trait en mode dessinateur d’exécution, pour la quarantaine de bonhommes qui attendent lundi dans l’atelier ferraillage/coffrage, les joies du télétravail (le travail qui s'invite à la maison, et non le travaille devant la télé comme certains seraient tenté de croire).

 

 

4 hours ago, didier said:

Il n'y a rien de personnel contre toi, c'est toi qui publies, c'est tout, du coup je te cite.

@didier Pas de soucis j'ai beaucoup d'humour, je ne me suis pas sentie visé d’ailleurs celle-ci:

8 hours ago, didier said:

... @VDH-Bruno
mais il existe et ça devrait pas pour l'égalité des chances.

lu comme ça, cela m'a beaucoup amusé😊

Après algorithme tel que proposé sur la page Wikipédia me semble pas incompatible avec un style impératif. Donc ça devrait aboutir bonne continuation à toi

Salutations

Bruno

[didier]

Coucou

Ça y est je livre et me délivre, c'est tellement rare que je participe à un challenge...

Exécution : le sablier flashe à partir de dix mille points

Rendu :

 

Code : avec commentaires (je demanderais volontiers aux autres de faire de même pour la formation des lecteurs)

Amicalement

(defun da:horaire (a b c / )
;| calcul de la superficie du triangle qui est négative si les sommets tournent en sens horaire
http://mathforum.org/library/drmath/view/54386.html  (0.5)(x1*y2 - y1*x2 -x0*y2 + y0*x2 + x0*y1 - y0*x1)
 |;
(minusp (* (- (+ (+ (- (- (* (car b) (cadr c))(* (cadr b) (car c)))(* (car a) (cadr c)))(* (cadr a) (car c)))(* (car a) (cadr b)))(* (cadr a) (car b)))0.5))
  );fin defun da:horaire

(defun da:makepoly2d (lstpt2d);dessin d'une poly 2D depuis une liste de points ((x1 y1)...(xn yn))
  (entmake
    (append
      (list
	'(0 . "LWPOLYLINE")
	'(100 . "AcDbEntity")
	'(100 . "AcDbPolyline")
	(cons 90 (length lstpt2d))
	'(70 . 1) ;1 pour fermée, 0 pour ouverte
	)
      (mapcar '(lambda (p) (cons 10 p)) lstpt2d)
      )
    )
  (princ)
  );fin de da:makepoly2d

;fonction principale qui attend en unique argument une liste de points ((x1 y1 z1)...(xn yn zn))
(defun daqh (lstpt / l3 l4 l5 l6 ldalsthull lstpt pt ptpi)

  ; tri par Y croissants
  (setq lst (vl-sort lstpt '(lambda(a b) (if (eq (cadr a) (cadr b))(< (car a) (car b))(< (cadr a) (cadr b))))))

  ; point avec le plus petit Y, point pivot
  (setq ptpi (car lst))

  ; liste des angles depuis pivot depuis l'axe des X
  (setq lda (mapcar '(lambda (X Y) (angle ptpi y)) (reverse (cdr (reverse lst))) (cdr lst)))
  
  ;création liste des points avec la valeur d'angle en premier (angle x y z)
  (setq l3 (mapcar 'cons lda (cdr lst)))

  ; tri par premier élément de la liste "angle"
  (setq l4 (vl-sort l3 '(lambda(a b) (if (eq (cdr a) (cdr b))(< (car a) (car b))(< (car a) (car b))))))

  ;liste des points triés en retirant l'angle en premier élément pour n'avoir que des points (x y z)
  (setq l5 (mapcar 'cdr l4))

  ;rajout du point pivot dans la liste des points
  (setq l6 (cons ptpi l5))

  (setq lsthull (list (cadr l6) (car l6)))
  ;liste sans les deux premiers points
  (setq l6 (cddr l6))

  ;test sens horaire 
  (foreach ptemp l6
    (while (and (cdr lsthull) (da:horaire (cadr lsthull) (car lsthull) ptemp))
      (setq lsthull (cdr lsthull))
      )
    (setq lsthull (cons ptemp lsthull))
    )

  ;(setq zaza (reverse lsthull))

  ;suppression des z dans la liste des points
  (setq lsthull (mapcar '(lambda (pt) (list (car pt) (cadr pt))) (reverse lsthull)))
  ;dessin de la poly 2D d'enveloppe
  (da:makepoly2d lsthull)
);fin defun daqh


;fonction pour test qui demande un choix de points filtré sur entité POINT
(defun c:TestdaQH ( / ss lstpt n)
  (prompt "\nChoix des points à cerner...\n")
  (setq ss (ssget (list '(0 . "POINT"))))
  (repeat (setq n (sslength ss))
    (setq lstpt (cons (cdr (assoc 10 (entget (ssname ss (setq n (1- n))))))lstpt))
  )
(daqh lstpt)
  );fin de c:TestdaQH

(prompt "\nTaper Testdaqh pour trouver l'enveloppe convexe des points\n")
(princ)


      

 

 

 

 

 

Amicalement