Author Topic: [Challenge] intersection of lines and planes  (Read 4381 times)

0 Members and 1 Guest are viewing this topic.

ElpanovEvgeniy

  • Water Moccasin
  • Posts: 1540
  • Moscow (Russia)
[Challenge] intersection of lines and planes
« on: January 20, 2011, 05:04:23 AM »
A new task to find the intersection point of plane and line.  :-)

arguments:
p1 p2 - the start and end points of the line
v - vector normal  planes (as dxf 210)
vp - point on  planes
f - the intersection between the test points or test extends line

Examples:
(test P1 P2 V VP FL) ;=>> point (intersection of lines and planes)
(test '(5. 5. 5.) '(10. 10. 10.) '(0. 0. 1.) '(0. 0. 2.) nil) ;=>> (2.0 2.0 2.0)
(test '(5. 5. 5.) '(10. 10. 10.) '(0. 0. 1.) '(0. 0. 2.) t) ;=>> nil
(test '(5. 5. 5.) '(10. 10. 10.) '(0. 0. 1.) '(0. 0. 2.) nil) ;=>> (2.0 2.0 2.0)
(test '(5. 5. 5.) '(10. 10. 10.) '(0. 0. 1.) '(0. 0. 2.) t) ;=>> (2.0 2.0 2.0)


« Last Edit: January 20, 2011, 07:15:16 AM by ElpanovEvgeniy »
蝸牛そろそろ登れ富士の山 /Kobayashi Issa/

qjchen

  • Bull Frog
  • Posts: 285
  • Best wishes to all
Re: [Challenge] intersection of lines and planes
« Reply #1 on: January 20, 2011, 07:12:22 AM »
THIS is my code
I am still testing

Code: [Select]
;;;;by qjchen@gmail.com
(defun test(P1 P2 V VP F / A B C D res x0 x1 x2 y0 y1 y2 z0 z1 z2)
(setq A (car V) B (cadr V) C (caddr V) x0 (car VP) y0 (cadr VP) z0 (caddr VP))
(setq x1 (car P1) y1 (cadr P1) z1 (caddr P1) x2 (car P2) y2 (cadr P2) z2 (caddr P2))
(if (/= (setq D (+ (* A (- x1 x2)) (* B (- y1 y2)) (* C (- z1 z2)))) 0.0)
  (progn
    (setq res
      (list (/ (+ (* A x0 (- x1 x2)) (* B (+ (* x1 (- y0 y2)) (* x2 (- y1 y0)) )) (* C (+ (* x1 (- z0 z2)) (* x2 (- z1 z0)) )) ) D)
            (/ (+ (* B y0 (- y1 y2)) (* C (+ (* y1 (- z0 z2)) (* y2 (- z1 z0)) )) (* A (+ (* y1 (- x0 x2)) (* y2 (- x1 x0)) )) ) D)
            (/ (+ (* C z0 (- z1 z2)) (* A (+ (* z1 (- x0 x2)) (* z2 (- x1 x0)) )) (* B (+ (* z1 (- y0 y2)) (* z2 (- y1 y0)) )) ) D)
      )
    )
    (if F (if (> (max (distance res P1) (distance res P2)) (distance P1 P2)) (setq res nil)))
  )
  (setq res nil)
)
res
)

(princ "\n (test '(5. 5. 5.) '(10. 10. 10.) '(0. 0. 1.) '(0. 0. 2.) nil) = \n")
(princ (test '(5. 5. 5.) '(10. 10. 10.) '(0. 0. 1.) '(0. 0. 2.) nil));Erase_DV

(princ "\n (test '(5. 5. 5.) '(10. 10. 10.) '(0. 0. 1.) '(0. 0. 2.) t) = \n")
(princ (test '(5. 5. 5.) '(10. 10. 10.) '(0. 0. 1.) '(0. 0. 2.) t));Erase_DV

(princ)

« Last Edit: January 20, 2011, 07:20:51 AM by qjchen »
http://qjchen.mjtd.com
My blog http://chenqj.blogspot.com (Chinese, can be translate into English)

gile

  • Water Moccasin
  • Posts: 2224
  • Marseille, France
Re: [Challenge] intersection of lines and planes
« Reply #2 on: January 20, 2011, 07:45:56 AM »
Hi,

It seems to me there's an error.
The intersection of a line defined by (5. 5. 5.) and (10. 10. 10.) or (0. 0. 0.) and (5. 5. 5.) (means its direction is (1. 1. 1.) and a plane which normal is (1. 0. 0.) and origin (0. 0. 2.) is (0. 0. 0.) rather than (2. 0. 0.).

IMO, this should be the results:

(test '(5. 5. 5.) '(10. 10. 10.) '(1. 0. 0.) '(2. 0. 0.) nil) ;=>> (2.0 2.0 2.0)
(test '(5. 5. 5.) '(10. 10. 10.) '(1. 0. 0.) '(2. 0. 0.) t) ;=>> nil
(test '(0. 0. 0.) '(5. 5. 5.) '(1. 0. 0.) '(2. 0. 0.) nil) ;=>> (2.0 2.0 2.0)
(test '(0. 0. 0.) '(5. 5. 5.) '(1. 0. 0.) '(2. 0. 0.) t) ;=>> (2.0 2.0 2.0)


Maybe I didn't understand the exercice ?...

EDIT: I undertand now. the VP point coordinates are defined about the plane coordinate system, I thaught it was defined in the WCS as the 2 other points.
Just replace (append p1 p2 norm org) by (append p1 p2 norm (trans org norm 0)) in the following code.

Anyway, here's a code which seems to work as I expect.

EDIT: corrected a typo (see here)

Code: [Select]
(defun ilp (p1 p2 norm org onSeg / x1 y1 z1 x2 y2 z2 xn yn zn xo yo zo scl)
  (mapcar 'set
  '(x1 y1 z1 x2 y2 z2 xn yn zn xo yo zo)
  (append p1 p2 norm org)
  )
  (if
    (and
      (not
(zerop
  (setq
    scl (+ (* xn (- x2 x1)) (* yn (- y2 y1)) (* zn (- z2 z1)))
  )
)
      )
      (or
(<= -1.
    (setq
      scl (/ (+ (* xn (- x1 xo)) (* yn (- y1 yo)) (* zn (- z1 zo))) scl)
    )
    0.
)
(not onSeg)
      )
    )
     (list
       (+ (* scl (- x1 x2)) x1)
       (+ (* scl (- y1 y2)) y1)
       (+ (* scl (- z1 z2)) z1)
     )
  )
)
« Last Edit: January 20, 2011, 10:17:12 AM by gile »
Speaking English as a French Frog

gile

  • Water Moccasin
  • Posts: 2224
  • Marseille, France
Re: [Challenge] intersection of lines and planes
« Reply #3 on: January 20, 2011, 08:46:12 AM »
Maybe it's more readable with 'factorised code' (but seems to run slower)

The same as upper: works with with all points and vector defined in the same coordinates system.

Code: [Select]
(defun DotProduct (v1 v2) (apply '+ (mapcar '* v1 v2)))

(defun ilp (p1 p2 norm org onSeg / scl)
  (if
    (and
      (/= 0 (setq scl (DotProduct norm (mapcar '- p2 p1))))
      (or
(<= -1. (setq scl (/ (DotProduct norm (mapcar '- p1 org)) scl)) 0.)
(not onSeg)
      )
    )
     (mapcar (function (lambda (x1 x2) (+ (* scl (- x1 x2)) x1))) p1 p2)
  )
)
« Last Edit: January 20, 2011, 09:02:42 AM by gile »
Speaking English as a French Frog

SOFITO_SOFT

  • Guest
Re: [Challenge] intersection of lines and planes
« Reply #4 on: January 20, 2011, 08:46:32 AM »
Hello:

easier with 3 points on the plane?

(Setq point-inter  (c:cal  "ilp (p1,p2,p3,p4,p5)"  )) <<<  THE POINT WHERE THE STRAIGHT P1 P2 intercepted the plan P3 P4  P5  
                
using the geo-cal.

Regards.                                                

gile

  • Water Moccasin
  • Posts: 2224
  • Marseille, France
Re: [Challenge] intersection of lines and planes
« Reply #5 on: January 20, 2011, 09:06:45 AM »
The geom cal runs much more slower than vector calculus an do not allow to evaluate if the returned point belongs to the segment.

With the plane defined by 3 points:

EDIT: no need to normalize the vector

Code: [Select]
(defun DotProduct (v1 v2) (apply '+ (mapcar '* v1 v2)))

(defun CrossProduct (v1 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)))
  )
)

(defun ilp3pts (p1 p2 p3 p4 p5 onSeg / norm)
  (if
    (and
      (setq norm (CrossProduct (mapcar '- p4 p3) (mapcar '- p5 p3)))
      (/= 0 (setq scl (DotProduct norm (mapcar '- p2 p1))))
      (or
(<= -1. (setq scl (/ (DotProduct norm (mapcar '- p1 p3))) scl) 0.)
(not onSeg)
      )
    )
     (mapcar (function (lambda (x1 x2) (+ (* scl (- x1 x2)) x1))) p1 p2)
  )
)
« Last Edit: January 20, 2011, 10:25:19 AM by gile »
Speaking English as a French Frog

SOFITO_SOFT

  • Guest
Re: [Challenge] intersection of lines and planes...One Refinement:
« Reply #6 on: January 20, 2011, 09:24:23 AM »
Refinement:
well, you can "shoot" from p1 to p2 and know very quickly if you "hit" in the triangle p3 p4 p5 (any UCS, you just have to use 3D points)
Code: [Select]
; return 0 or NIL or T ( vertex , out , in )
; work in all ucs with points ( x y z )
; toler-igual = tolerance little computation error ( 1e-8 ? ) 

(DEFUN PUNTO_DENTR_TRIAN (PP PA PB PC / )  <<  point PP in triangle PA PB PC ?

  (SETQ RADIO-A (MODUL_VECTO (RESTA_2VECTO PA PP)))  <<< module vector substraction (pa - pp)
  (SETQ RADIO-B (MODUL_VECTO (RESTA_2VECTO PB PP)))
  (SETQ RADIO-C (MODUL_VECTO (RESTA_2VECTO PC PP)))

  (IF (AND
(NOT (ZEROP RADIO-A))    <<< the point is in vertex triangle
(NOT (ZEROP RADIO-B))
(NOT (ZEROP RADIO-C))
      )
    (PROGN

      (SETQ TMP1 (DIVID_VECTO_NUMER (RESTA_2VECTO PA PP) RADIO-A)) <<< divide the vector substraction (pa - pp)  / radio-a
      (SETQ TMP2 (DIVID_VECTO_NUMER (RESTA_2VECTO PB PP) RADIO-B))
      (SETQ TMP3 (DIVID_VECTO_NUMER (RESTA_2VECTO PC PP) RADIO-C))

      (SETQ ANG1 (ARCO_COS (PRODU-ESCAL_2VECTO TMP1 TMP2)))  <<< arc_cos of  cross prod  of 2 vectors
      (SETQ ANG2 (ARCO_COS (PRODU-ESCAL_2VECTO TMP2 TMP3)))
      (SETQ ANG3 (ARCO_COS (PRODU-ESCAL_2VECTO TMP3 TMP1)))

      (EQUAL *PI2* (+ ANG1 ANG2 ANG3) TOLER-IGUAL)  <<< if adition of angles is 2pi is IN else if OUT

    )
    0.0  <<<< is in vertex
  )
)

SOFITO_SOFT

  • Guest
Re: [Challenge] intersection of lines and planes
« Reply #7 on: January 20, 2011, 09:48:41 AM »
Code: [Select]
[quote author=gile link=topic=36756.msg417558#msg417558 date=1295532405]
 much more slower  ??  :-o :-o ???

to evaluate if the returned point belongs to the segment. <<<<   is  olny one line ( 4 instructios ) verifi it !

( = ( + ( distance p1 point-inter ) ( distance p2 point-inter ) )  ( distance p1 p2 )  )

)
[/quote]

ElpanovEvgeniy

  • Water Moccasin
  • Posts: 1540
  • Moscow (Russia)
Re: [Challenge] intersection of lines and planes
« Reply #8 on: January 20, 2011, 10:01:32 AM »
  :oops:
I apologize, I really made mistakes in the test examples!

new examples:
Code: [Select]
(setq p1 '(1. 2. 3.)p2 '(4. 5. 6.) vp '(7. 8. 9.) v '(0. 0. 1.))
(test p1 p2 v vp nil); (7.0 8.0 9.0)
(test p1 p2 v vp t); nil

(setq p1 '(1. 2. 3.)p2 '(4. 5. 6.) vp '(7. 8. 9.) v '(0.57735 0.57735 0.57735))
(test p1 p2 v vp nil); (7.0 8.0 9.0)
(test p1 p2 v vp t); nil

(setq p1 '(1. 2. 3.)p2 '(4. 5. 6.) vp '(1. 3. 5.) v '(0.57735 0.57735 0.57735))
(test p1 p2 v vp nil); (2.0 3.0 4.0)

(setq p1 '(1. 2. 3.)p2 '(4. 5. 6.) vp '(1. 5. 9.) v '(0.57735 0.57735 0.57735))
(test p1 p2 v vp nil); (4.0 5.0 6.0)

(setq p1 '(1. 2. 3.)p2 '(4. 5. 6.) vp '(1. 3. 5.) v '(0. 0. 1.))
(test p1 p2 v vp nil); (3.0 4.0 5.0)

(setq p1 '(1. 2. 3.)p2 '(4. 5. 6.) vp '(1. 3. 5.) v '(0. 1. 0.))
(test p1 p2 v vp nil); (2.0 3.0 4.0)

(setq p1 '(1. 2. 3.)p2 '(4. 5. 6.) vp '(1. 3. 5.) v '(1. 0. 0.))
(test p1 p2 v vp nil); (1.0 2.0 3.0)

(setq p1 '(1. 2. 3.)p2 '(4. 5. 6.) vp '(1. 3. 5.) v '(0.707107 0.707107 0.0))
(test p1 p2 v vp nil); (1.5 2.5 3.5)
蝸牛そろそろ登れ富士の山 /Kobayashi Issa/

ElpanovEvgeniy

  • Water Moccasin
  • Posts: 1540
  • Moscow (Russia)
Re: [Challenge] intersection of lines and planes
« Reply #9 on: January 20, 2011, 10:16:58 AM »
THIS is my code
I am still testing


Hi Chen!
Your code works fine!
but my code is shorter...   :-)
蝸牛そろそろ登れ富士の山 /Kobayashi Issa/

ElpanovEvgeniy

  • Water Moccasin
  • Posts: 1540
  • Moscow (Russia)
Re: [Challenge] intersection of lines and planes
« Reply #10 on: January 20, 2011, 10:35:58 AM »
Hi gile! :)

Your program works great!


ps. I was worried by the number of variables in the code:
Code: [Select]
(mapcar 'set
  '(x1 y1 z1 x2 y2 z2 xn yn zn xo yo zo)
  (append p1 p2 norm org)
  )
蝸牛そろそろ登れ富士の山 /Kobayashi Issa/

gile

  • Water Moccasin
  • Posts: 2224
  • Marseille, France
Re: [Challenge] intersection of lines and planes
« Reply #11 on: January 20, 2011, 10:40:26 AM »
much more slower  ??  :-o :-o ???

to evaluate if the returned point belongs to the segment. <<<<   is  olny one line ( 4 instructios ) verifi it !

( = ( + ( distance p1 point-inter ) ( distance p2 point-inter ) )  ( distance p1 p2 )  )

)

I was just saying the geomcal is slow.
Look at this little benchmark using the upper ilp3pts (which evaluates if the point belongs to the segment) and only the geomcal (without evaluating where is the point).
Using 5 points:
(setq p1 '(1. 2. 3.) p2 '(4. 5. 6.) p3 '(0. 0. 0.) p4 '(0. 1. 0.) p5 '(0. 0. 1.))
(c:cal "ilp(p1,p2,p3,p4,p5)") returns (0.0 1.0 2.0)
(ilp3pts p1 p2 p3 p4 p5 nil) returns (0.0 1.0 2.0) as well
(ilp3pts p1 p2 p3 p4 p5 T) returns nil

Code: [Select]
_$ (benchmark '((ilp3pts p1 p2 p3 p4 p5 nil) (c:cal "ilp(p1,p2,p3,p4,p5)")))
Benchmarking ..................Elapsed milliseconds / relative speed for 32768 iteration(s):

    (ILP3PTS P1 P2 P3 P4 P5 nil).......3105 / 4.21 <fastest>
    (C:CAL "ilp(p1,p2,p3,p4,p5)").....13057 / 1.00 <slowest>
Speaking English as a French Frog

qjchen

  • Bull Frog
  • Posts: 285
  • Best wishes to all
Re: [Challenge] intersection of lines and planes
« Reply #12 on: January 20, 2011, 10:41:43 AM »
THIS is my code
I am still testing


Hi Chen!
Your code works fine!
but my code is shorter...   :-)

Hi, Evgeniy, I write a more long code, but I think it is more clear as like Gile's code (Certainly Gile's code is much better than mine) , with Vector calculation

Base on this wiki
http://en.wikipedia.org/wiki/Line-plane_intersection

Code: [Select]
;;;;vec plus
(defun q:vec:+(v1 v2)
  (mapcar '+ v1 v2)
)
;;;;vec substract
(defun q:vec:-(v1 v2)
  (mapcar '- v1 v2)
)
;;;;vec plus constant
(defun q:vec:*c(v a)
  (mapcar '(lambda(x) (* x a)) v)
)
;;;;vec dot product
(defun q:vec:dot*(v1 v2)
  (apply '+ (mapcar '* v1 v2))
)
;;;;Normalize a vec
(defun q:vec:Norm(v / l)
  (if (not (zerop (setq l (distance '(0 0 0) v))))
  (mapcar '(lambda(x) (/ x l)) v))
)
;;;; by qjchen@gmail.com To calculate a line and a plane intersection
(defun q:geo:line-int-plane(P1 P2 V VP F / d l n p)
 (setq n (q:vec:Norm V) l (q:vec:Norm (q:vec:- P2 P1)))
 (cond
  ((zerop (q:vec:dot* l n)) nil)
  (T (setq d (/ (q:vec:dot* (q:vec:- VP P1) n) (q:vec:dot* l n)))
     (setq p (q:vec:+ P1 (q:vec:*c l d)))
     (if (and F (< (q:vec:dot* (q:vec:- p P1) (q:vec:- P2 p)) 0))  nil p))
 )
)
 
(q:geo:line-int-plane '(5. 5. 5.) '(10. 10. 10.) '(1. 1. 1.) '(0 0 2.0) NIL)
« Last Edit: January 20, 2011, 11:29:03 AM by qjchen »
http://qjchen.mjtd.com
My blog http://chenqj.blogspot.com (Chinese, can be translate into English)

gile

  • Water Moccasin
  • Posts: 2224
  • Marseille, France
Re: [Challenge] intersection of lines and planes
« Reply #13 on: January 20, 2011, 11:01:32 AM »
Hi gile! :)

Your program works great!


ps. I was worried by the number of variables in the code:
Code: [Select]
(mapcar 'set
  '(x1 y1 z1 x2 y2 z2 xn yn zn xo yo zo)
  (append p1 p2 norm org)
  )

Tanks.

This way is less concise and readable the solution I posted above (using DotProduct routine and mapcar statements) which uses the same algorithm but it seems to run a little faster.
Speaking English as a French Frog

chlh_jd

  • Guest
Re: [Challenge] intersection of lines and planes
« Reply #14 on: January 20, 2011, 02:11:36 PM »
Hi all !
If use following method ,why take out a wrong result ?

Code: [Select]
;;;Use Vetor
(defun ilp (p1 p2 v vp org / ms lob pob pt)
  (setq ms (vla-get-ModelSpace
     (vla-get-ActiveDocument (vlax-get-acad-object))
   )
  )
  (setq lob (apply 'vla-AddLine
   (cons ms
(mapcar 'vlax-3d-point (list p1 p2))
   )
    )
  )
  (setq
    pob (apply 'vla-add3dface
       (cons ms
     (mapcar 'vlax-3d-point (list '(0. 0. 0.) v vp vp))
       )
)
  )
  (setq pt (Vla-IntersectWith
     lob
     pob
     (if org
       acExtendOtherEntity
       acExtendBoth
     )
   )
  )
  (vla-delete lob)
  (vla-delete pob)
  (vlax-safearray->list (vlax-variant-value pt)) 
)
;;;Use 3pts
(defun ilp3pts (p1 p2 p3 p4 p5 org / ms lob pob pt)
  (setq ms (vla-get-ModelSpace
     (vla-get-ActiveDocument (vlax-get-acad-object))
   )
  )
  (setq lob (apply 'vla-AddLine
   (cons ms
(mapcar 'vlax-3d-point (list p1 p2))
   )
    )
  )
  (setq
    pob (apply 'vla-add3dface
       (cons ms
     (mapcar 'vlax-3d-point (list p3 p4 p5 p5))
       )
)
  )
  (setq pt (Vla-IntersectWith
     lob
     pob
     (if org
       acExtendOtherEntity
       acExtendBoth
     )
   )
  )
  (vla-delete lob)
  (vla-delete pob)
  (vlax-safearray->list (vlax-variant-value pt))
)