Vecteur, normale, angle 3D ...

[(gile)]

Pour certains calculs ou comparaisons en 3D, il peut être intéressant d'utiliser des vecteurs normés (vecteur de une unité).

 

Mis à jour le 07/10/06 (nouvelles routines)

 

Pour trouver le vecteur normé d'un segment défini par ses extrémités, on "déplace" le départ en (0.0 0.0 0.0) et l'extrémité à une unité du départ :

 

;;; VEC1 Retourne le vecteur normé (1 unité) de p1 à p2 (nil si p1 = p2)

(defun vec1 (p1 p2)
 (if (not (equal p1 p2 1e-009))
   (mapcar '(lambda (x1 x2)
       (/ (- x2 x1) (distance p1 p2))
     )
    p1
    p2
   )
 )
)

;; ou

(defun vec1 (p1 p2)
 (vunit (mapcar '- p2 p1))
)

 

Produit scalaire de deux vecteurs.

Le résultat est un nombre réel.

Si le produit scalaire de deux vecteur est égal à 0, les vecteurs sont perpendiculaires.

 

;;; VXV Retourne le produit scalaire (réel) de deux vecteurs

(defun vxv (v1 v2)
 (apply '+ (mapcar '* v1 v2))
)

 

Le produit scalaire d'un vecteur par lui même retourne le carré de sa norme (sa longueur).

 

;;; VLEN Retourne la longueur (norme) d'un vecteur

(defun vlen (v)
 (sqrt (vxv v v))
)

 

Pour trouver le vecteur unitaire de même sens et direction qu'un vecteur quelconque.

 

;;; VUNIT Retourne le vecteur unitaire d'un vecteur

(defun vunit (v / l)
 (if (/= 0 (setq l (vlen v)))
   (mapcar '(lambda (x) (/ x l)) v)
 )
)

 

Produit vectoriel de deux vecteurs.

Le résultat est un vecteur perpendiculaire aux deux vecteurs.

 

;;; V^V Retourne le produit vectoriel (vecteur) de deux vecteurs

(defun v^v (v1 v2)
 (if (inters '(0 0 0) v1 '(0 0 0) v2)
   (list (- (* (cadr v1) (caddr v2)) (* (caddr v1) (cadr v2)))
  (- (* (caddr v1) (car v2)) (* (car v1) (caddr v2)))
  (- (* (car v1) (cadr v2)) (* (cadr v1) (car v2)))
   )
   ;; '(0.0 0.0 0.0)
 )
)

 

Il peut être intéressant aussi de connaître le vecteur normal (ou la normale) d'un un plan défini par 3 points, pour déterminer le code dxf 210 ou la propriété Normal en Vlisp (les coordonnées des points 3d doivent alors être définies dans le SCG).

 

Norm_3dpts utilise la fonction LINEARP définie plus bas.

 

;;; NORM_3PTS retourne le vecteur normal du plan défini par 3 points 3D

(defun norm_3pts (p0 p1 p2 / norm)
 (foreach v '(p1 p2)
   (set v (mapcar '- (eval v) p0))
 )
 (mapcar '(lambda (x) (/ x (distance '(0 0 0) norm)))
  (setq norm (v^v p1 p2))
 )
)

;; ou

(defun norm_3pts (p0 p1 p2 / norm)
 (vec1	'(0 0 0)
(v^v (vec1 p0 p1) (vec1 p0 p2))
 )
) 

 

La direction du vecteur normal est fonction de l'ordre dans lequel sont spécifiés les points : p0 comme origine et p1 et p2 dans le sens trigonométrique.

Exemple :

(norm_3pts '(0 0 0) '(1 0 0) '(0 1 0)) retourne (0.0 0.0 1.0)

(norm_3pts '(0 0 0) '(0 1 0) '(1 0 0)) retourne (0.0 0.0 -1.0)

 

http://img347.imageshack.us/img347/3500/normal3as.png

 

Les deux routines ci-dessus retournent les mêmes résultats que les expressions vec1(dep,ext) et nor(p0,p1,p2) employés avec la calculatrice géométrique.

 

Exemples d'utilisation des vecteurs normés :

 

Ces routines utilisant les vecteurs fonctionnent quelque soit le système de coordonnées courant et quelque soit celui dans lequel les coordonnées des points sont définies (le même pour tous les points, quand même !).

 

Pour trouver la valeur d'un angle, définis par trois points 3D :

 

;;; ANGLE_3PTS Retourne l'angle (radians) défini par son sommet et deux points
;;; L'angle retourné est toujours positif et inférieur à pi radians.

(defun angle_3pts (som p1 p2 / v1 v2)
 (if (and
(setq v1 (vec1 som p1))
(setq v2 (vec1 som p2))
     )
   (cond
     ((equal v1 v2 1e-009) 0.0)
     ((equal v1 (mapcar '- v2) 1e-009) pi)
     (T (* 2 (asin (/ (distance v1 v2) 2))))
   )
   (Prompt "\nErreur: points confondus.")
 )
) 

;;; ou, par le "théorème d'Al Kashi"

(defun [b]angle_3pts[/b] (som p1 p2 / d1 d2 d3)
 (setq	d1 (distance som p1)
d2 (distance som p2)
d3 (distance p1 p2)
 )
 (if (and (< 0 d1) (< 0 d2))
   (acos (/ (+ (* d1 d1) (* d2 d2) (- (* d3 d3)))
     (* 2 d1 d2)
  )	
   )
 )
)

 

Pour trouver la valeur d'un angle, définis par quatre points 3D (les extrémités de 2 segments) :

 

;;; ANGLE_4PTS Retourne l'angle entre 2 segments définis par leurs extrémités respectives 
(defun angle_4pts (p1 p2 p3 p4 / som v1 v2)
 (cond
   ((or (equal p1 p2 1e-009)
 (equal p3 p4 1e-009)
    )
    (prompt "\nErreur: points confondus.")
   )
   ((linearp (remove_doubles (list p1 p2 p3 p4)))
    (cond
      ((equal (vec1 p1 p2) (vec1 p3 p4) 1e-009) 0.0)
      ((equal (vec1 p1 p2) (vec1 p4 p3) 1e-009) pi)
    )
   )
   ((cond
      ((setq som (inters p1 p2 p3 p4 nil))
(if (null (setq v1 (vec1 som p1)))
  (setq v1 (vec1 som p2))
)
(if (null (setq v2 (vec1 som p3)))
  (setq v2 (vec1 som p4))
)
(* 2 (asin (/ (distance v1 v2) 2)))
      )
      (T
(if (coplanp (list p1 p2 p3 p4))
  (prompt "\nErreur: segments parallèles.")
  (prompt "\nErreur: points non coplanaires.")
)
      )
    )
   )
 )
)

 

Les trois dernières routines utilisent les fonctions ASIN et ACOS pour retourner un arc sinus ou cosinus, curieusement asin et acos ne sont pas des fonctions AutoLISP prédéfinies :

 

;;; ASIN et ACOS Retournent l'arc sinus ou l'arc cosinus du nombre, en radians

(defun ASIN (num)
 (if (<= -1 num 1)
   (atan num (sqrt (- 1 (expt num 2))))
   (princ
     "\nErreur: L'argument pour ASIN doit être compris entre -1 et 1"
   )
 )
)

(defun ACOS (num)
 (if (<= -1 num 1)
   (atan (sqrt (- 1 (expt num 2))) num)
   (princ
     "\nErreur: L'argument pour ACOS doit être compris entre -1 et 1"
   )
 )
) 

 

L'utilisation de vecteurs permet aussi de réaliser des tests sur des listes de points 3D :

 

Evaluer si les points d'une liste sont colinéaires :

 

;;; LINEARP Retourne T si tous les points de la liste sont alignés

(defun linearp (lst)
 (cond
   ((= 2 (length lst)) T)
   ((or (equal	(vec1 (car lst) (cadr lst))
	(vec1 (car lst) (caddr lst))
	1e-009
 )
 (equal	(vec1 (car lst) (cadr lst))
	(vec1 (caddr lst) (car lst))
	1e-009
 )
    )
    (linearp (cdr lst))
   )
 )
)

 

Evaluer si tous les points d'une liste sont coplanaires :

 

;;; COPLANP Retourne T si tous les points de la liste sont coplanaires
(defun coplanp (lst)
 (cond
   ((= 3 (length lst)) T)
   ((or (inters (car lst) (cadr lst) (caddr lst) (cadddr lst) nil)
 (equal	(vec1 (car lst) (cadr lst))
	(vec1 (caddr lst) (cadddr lst))
	1e-009
 )
 (equal	(vec1 (car lst) (cadr lst))
	(vec1 (cadddr lst) (caddr lst))
	1e-009
 )
    )
    (coplanp (cdr lst))
   )
 )
)

 

PS : J'ai ajouté "angle_4pts" et des contrôles dans "angle_3pts" et "norm_3pts". [Edité le 5/3/2006 par (gile)][Edité le 17/3/2006 par (gile)]

[Edité le 20/3/2006 par (gile)]

[Edité le 7/10/2006 par (gile)]

[(gile)]

J'ai modifié ANGLE_4PTS pour corriger une inexactitude en cas de points alignés.

 

Pour les fonctions LINEARP et COPLANP il peut être prudent de supprimer les doublons dans les listes de points :

 

(linearp (remove_doubles (list p1 p2 p3 p4)))

 

avec la fonction REMOVE_DOUBLES :

 

;;; REMOVE_DOUBLES - Suprime tous les doublons d'une liste

(defun REMOVE_DOUBLES (lst)
 (cond
   ((atom lst) lst)
   (T
    (cons (car lst) (REMOVE_DOUBLES (vl-remove (car lst) lst)))
   )
 )
)

[(gile)]

Un petit exemple d'utilisation des routines ci-dessus (à placer, par exemple, dans la barre d'outils ou le menu Renseigements)

 

ANGLE_3D affiche sur la ligne de commande l'angle formé par 4 points 3D (respectivement les extrémité de deux segments) indépendamment du SCU courant.

 

Si le premier et le troisième points sont confondus (sommet de l'angle) il est possible de faire ENTER ou clic droit au lieu de le re-spécifier.

 

Si les segments sont parallèles ou non coplanaires il est retourné un message d'erreur.

 

Les routines ANGLE_4PTS, VEC1, LINEARP, COPLANP, ASIN, REMOVE_DOUBLES doivent être chargées.

 

(defun c:angle_3d (/ pt1 pt2 pt3 pt4 ang)
 (initget 1)
 (setq pt1 (getpoint "\nPremier point: "))
 (initget 1)
 (setq pt2 (getpoint pt1 "\nSecond point: "))
 (grdraw pt1 pt2 1)
 (if (not
(setq pt3 (getpoint "\nTroisième point ou :"))
     )
   (setq pt3 pt1)
 )
 (initget 1)
 (setq pt4 (getpoint pt3 "\nQuatrième point: "))
 (grdraw pt3 pt4 1)
 (setq ang (angle_4pts pt1 pt2 pt3 pt4))
 (cond
   ((numberp ang)
    (if (= (getvar "ANGDIR") 1)
      (setq ang (- ang))
    )
    (princ
      (strcat "\nAngle = "
       (angtos (+ (getvar "ANGBASE") ang))
       " ou "
       (angtos (+ (getvar "ANGBASE") (- (* 2 pi) ang)))
      )
    )
   )
 )
 (redraw)
 (princ)
) 

[(gile)]

Suite à ce sujet, deux autres routines qui peuvent être utiles en 3D.

 

Pour trouver l'élévation d'un point par rapport à un plan :

 

;;; [b]ELEV[/b] retourne l'élévation point pt par rapport au plan défini par p1 p2 p3

(defun [b]elev[/b] (pt p1 p2 p3)
 (* (cos ([b]angle_3pts[/b] p1 (mapcar '+ p1 ([b]norm_3pts[/b] p1 p2 p3)) pt))
    (distance p1 pt)
 )
) 

 

Pour obtenir les coordonnées de la projection orthogonale d'un point sur un plan :

 

;;; [b]PROJ_PT[/b] Retourne les coordonnées de la projection orthogonale du point pt
;;; sur le plan défini par les points p1 p2 p3.

(defun [b]proj_pt[/b] (pt p1 p2 p3)
 (mapcar '-
  pt
  (mapcar '(lambda (x) (* x ([b]elev[/b] pt p1 p2 p3)))
	  ([b]norm_3pts[/b] p1 p2 p3)
  )
 )
) 

 

Avec ces deux routines, on peut aisément obtenir les mêmes résultats qu'avec les fonctions dpp et ilp de la calculatrice géométrique

 

Distance entre le point pt et le plan p1 p2 p3 :

(abs (elev pt p1 p2 p3))

 

Intersection de la ligne p1 p2 et du plan p3 p4 p5 :

 

(inters p1

p2

(proj_pt p1 p3 p4 p5)

(proj_pt p2 p3 p4 p5)

nil

)