TheSwamp
Code Red => AutoLISP (Vanilla / Visual) => Topic started by: ElpanovEvgeniy 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)
-
THIS is my code
I am still testing
;;;;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)
-
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 (http://www.theswamp.org/index.php?topic=36756.msg417547#msg417547))
(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)
)
)
)
-
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.
(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)
)
)
-
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.
-
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
(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)
)
)
-
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)
; 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
)
)
-
[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]
-
:oops:
I apologize, I really made mistakes in the test examples!
new examples:
(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)
-
THIS is my code
I am still testing
Hi Chen!
Your code works fine!
but my code is shorter... :-)
-
Hi gile! :)
Your program works great!
ps. I was worried by the number of variables in the code:
(mapcar 'set
'(x1 y1 z1 x2 y2 z2 xn yn zn xo yo zo)
(append p1 p2 norm org)
)
-
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
_$ (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>
-
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
;;;;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)
-
Hi gile! :)
Your program works great!
ps. I was worried by the number of variables in the code:
(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.
-
Hi all !
If use following method ,why take out a wrong result ?
;;;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))
)
-
Don't have time to play :-(
But have posted similar code here (http://www.cadtutor.net/forum/showthread.php?56064-Intersection-Line-amp-Rectangle)
-
(defun ilp ...
Very nice Gile :-)
-
my variant:
(defun inters_w-line-2p (P1 P2 V VP fl / P)
(if (setq p ((lambda (a b c)
(inters a b (list (car a) (cadr a) c) (list (car b) (cadr b) c) fl)
)
(trans p1 0 v)
(trans p2 0 v)
(caddr (trans vp 0 v))
)
)
(trans p v 0)
)
)
-
my variant:
(defun inters_w-line-2p (P1 P2 V VP fl / P)
(if (setq p ((lambda (a b c)
(inters a b (list (car a) (cadr a) c) (list (car b) (cadr b) c) fl)
)
(trans p1 0 v)
(trans p2 0 v)
(caddr (trans vp 0 v))
)
)
(trans p v 0)
)
)
very funny idea !
I learn a lot from this code for trans
-
Thank you!
This code is faster than the mathematical algorithms... :-)
-
About line-plane intersection from stig madsen
http://www.theswamp.org/index.php?topic=2505.0 (http://www.theswamp.org/index.php?topic=2505.0)
mardi
-
Very nice trick Evgeniy 8-)
-
Very nice trick Evgeniy 8-)
x2 :-)
-
Excellent solution Evgeniy !!
-
I'm glad you liked this theme!