Author Topic: intersection  (Read 11618 times)

0 Members and 1 Guest are viewing this topic.

CAB

  • Global Moderator
  • Seagull
  • Posts: 10401
intersection
« Reply #15 on: April 22, 2004, 11:49:37 AM »
You are welcome.

Not that it helped you much, sounds like you are way ahead of me. :)

CAB
I've reached the age where the happy hour is a nap. (°¿°)
Windows 10 core i7 4790k 4Ghz 32GB GTX 970
Please support this web site.

SMadsen

  • Guest
intersection
« Reply #16 on: April 22, 2004, 03:25:34 PM »
If anyone's interested, here's some non-matrix math behind plane calculations.

Some of the utility functions are good to have  :)

To see what calculations are available, run C:TEST. I've tried to comment it.

Code: [Select]

;;; VECTORS
;;; Asks for origin and 2 points on a plane. Just a utility
;;; function to get started. Returns origin and unit vectors
;;; for the 2 points.
(defun vectors (/ p0 p1 p2 u v)
  (and (setq p0 (getpoint "\nOrigin: "))
       (setq p1 (getpoint p0 "\nVector 1: ")
             p2 (getpoint p0 "\nVector 2: ")
       )
       (setq u (mapcar '- p1 p0)
             v (mapcar '- p2 p0)
       )
  )
  (list p0 u v)
)

;;; DOT
;;; Find dot product of two vectors, u and v
(defun dot (u v)
  (+ (* (car u) (car v))
     (* (cadr u) (cadr v))
     (* (caddr u) (caddr v))
  )
)

;;; CROSS
;;; Find cross product of two vectors, u and v
(defun cross (u v / x y z)
  (setq x (- (* (cadr u) (caddr v))
             (* (caddr u) (cadr v))
          )
        y (- (* (caddr u) (car v))
             (* (car u) (caddr v))
          )
        z (- (* (car u) (cadr v))
             (* (cadr u) (car v))
          )
  )
  (list x y z)
)

;;; ONPLANE
;;; Determines if point pt is on the plane specified by
;;; the planes normal vector n
;;; If ONPLANE returns 0.0 then pt is on plane
(defun onPlane (n pt d)
  (+ (* (car n) (car pt))
     (* (cadr n) (cadr pt))
     (* (caddr n) (caddr pt))
     (- d)
  )
)

;;; DISTTOPLANE
;;; Calculates Z distance from a point to a plane specified
;;; by the planes normal vector n
(defun distToPlane (n pt d)
  (/ (+ (* (car n) (car pt))
        (* (cadr n) (cadr pt))
        (* (caddr n) (caddr pt))
        d
     )
     (sqrt (+ (expt (car n) 2.0)
              (expt (cadr n) 2.0)
              (expt (caddr n) 2.0)
           )
     )
  )
)

;;; GETINTERS
;;; Calculates intersection point based on unit vector pn, point pt
;;; and intersection parameter si
(defun getInters (pt pn si)
  (mapcar '+ (mapcar (function (lambda (n) (* n si))) pn) pt)
)


(defun C:TEST (/      dnum   dp1    dp2    n      p0     p01    p1
               p1-plane      p2     p2-plane      pint   sd     sden
               si     sn     snum   snumalt       tden   tnum   tnumalt
               u      v      vects  w
              )
  (setq vects (vectors)                 ; get origin and two unit vectors of a plane
        p0    (car vects)               ; assign origin to p0
        u     (cadr vects)              ; assign 1st unit vector to u
        v     (caddr vects)             ; assign 2nd unit vector to v
        n     (cross u v)               ; find normal to plane
        dnum  (dot n p0)                ; calculate d in ax+by+cz+d
  )

  ;; find denominators for s and t in the parametric plane equation
  ;; V(s,t)=Vo+su+tv. Check for non-zero values before proceding.
  (cond ((and (not (equal (setq sden (dot u (cross n v))) 0.0 1E-6))
              (not (equal (setq tden (dot v (cross n u))) 0.0 1E-6))
         )
         ;; get points for line in space
         (and (setq p1 (getpoint "\nSpecify first point on line: "))
              (setq p2 (getpoint p1 "\nSecond point on line: "))
              (setq w (mapcar '- p2 p1)) ; find unit vector for p1-p2
         )

         ;; check if p1 and p2 is on the plane p0-u-v
         (setq p1-plane (onPlane n p1 dnum)
               p2-plane (onPlane n p2 dnum)
         )

         ;; calculate s and t in the parametric equation
         (setq snum (/ (dot w (cross n v)) sden)
               tnum (/ (dot w (cross n u)) tden)
         )

         ;; just an alternative way of calc'ing s and t:
         (setq snumalt (/ (- (* (dot u v) (dot w v)) (* (dot v v) (dot w u)))
                          (- (expt (dot u v) 2.0) (* (dot u u) (dot v v)))
                       )
               tnumalt (/ (- (* (dot u v) (dot w u)) (* (dot u u) (dot w v)))
                          (- (expt (dot u v) 2.0) (* (dot u u) (dot v v)))
                       )
         )

         ;; distance from p1 to plane p0-u-v:
         (setq dp1 (distToPlane n p1 dnum))
         ;; distance from p2 to plane p0-u-v:
         (setq dp2 (distToPlane n p2 dnum))

         ;; just some output:
         (mapcar 'princ
                 (list "\nDistance from p1 to plane: "
                       dp1
                       "\nDistance from p2 to plane: "
                       dp2
                 )
         )

         ;; is p1-p2 parallel to plane?
         (and (equal (dot n w) 0.0 1e-6)
              (princ "\nLine is parallel with plane")
         )

         ;; get unit vector p0 to p1
         (setq p01 (mapcar '- p1 p0))
         ;; check denominator for intersection parameter si
         (cond
           ((not (zerop (setq sn (- (dot n p01))
                              sd (dot n w)
                        )
                 )
            )
            ;; calculate intersection parameter si
            (setq si (/ sn sd))
            ;; get intersection point pi based on p1, unit vector w and si
            (setq pint (getInters p1 w si))
            ;; investigate parameter si:
            (cond ((or (equal si 0.0 1E-6)
                       (equal si 1.0 1E-6)
                   )
                   (princ "\nLine touches on plane")
                  )
                  ((< 0.0 si 1.0) (princ "\nIntersection is physical"))
                  (T (princ "\nIntersection is non-physical"))
            )
           )
           ((princ "\nNo intersection"))
         )
        )
        (T (princ "\nVectors do not describe a plane"))
  )
  (terpri)
  (if pint pint (princ))
)

Fuccaro

  • Guest
intersection
« Reply #17 on: April 23, 2004, 03:21:06 AM »
I like the matrix methods. If you write a routine to solve the 3x3 and the 4x4 matrix, you can use it in almost all your matrix-geometry programs...
SMadsen
I can't completely understand how GETINTERS works, but still searching...
Probably I will follow you and I will create a file with all the geometric routines I will ewer need.
Have a nice day!

SMadsen

  • Guest
intersection
« Reply #18 on: April 23, 2004, 05:56:03 AM »
CAB,
Just had some time to see what you posted. It's really cool. Thanks for putting it up here.

Fuccaro,
I'm not good at matrix math so if I can avoid it, well ...
I have been writing some transformation and stuff and it's fun to mess with once it gets going but I have to reflush the brain cells before it gets going :roll:

The GETINTERS is just a shortcut for getting the point that si describes along the vector p1p2: p1+s*w (with the variable names used in C:TEST).
I'm sure you know this but si is the ratio along the vector at which the vector intersects the plane. That's why 0<=si<=1 describes an intersection physically on the vector where si is the distance ratio from p1 to p2. E.g. if si=0.5 then it intersects at the midpoint of p1p2. si>1 means that the intersection lies at an extension of p1p2. All that GETINTERS does is to find that distance and move it back into p1:

w = p2 minus p1 = unit vector, i.e. p1p2 transposed into origin (it's not really a unit vector but the actual vector! I just call them unit vectors)

Distance to intersection = si * w. Moving it back to p1 is merely an addition of coordinates p1+(si*w), which is what MAPCAR does in GETINTERS. Here's an illustration of the process:


SMadsen

  • Guest
intersection
« Reply #19 on: April 23, 2004, 01:03:16 PM »
Just for fun, I tried to implement the ray-point-in-poly method to determine if an intersection point between a line and a 3dface is within the 3dface.

There are some problems with points directly on the edges (as always with this method), but it should work.
It only shoots out one ray in a distance of 100 times one edge but that should be enough for a rectangular object.

Code: [Select]
(defun C:TEST (/ a    dnum vint ent  entl tnum i    n    p0   p01  p1
                 w    uint plst sden pint tden v    u    si   sn   sd
                 p2   snum)
  (and (setq ent (car (entsel "\nSelect 3Dface: ")))
       (setq entl (entget ent))
       (= (cdr (assoc 0 entl)) "3DFACE")
       (foreach n '(10 11 12 13 10)
         (setq plst (cons (cdr (assoc n entl)) plst))
       )
       (setq ent (car (entsel "\nSelect line: ")))
       (setq entl (entget ent))
       (= (cdr (assoc 0 entl)) "LINE")
       (setq p1   (cdr (assoc 10 entl))
             p2   (cdr (assoc 11 entl))
             w    (mapcar '- p2 p1)
             p0   (cadr plst)
             u    (mapcar '- (car plst) p0)
             v    (mapcar '- (caddr plst) p0)
             n    (cross u v)
             dnum (dot n p0)
       )
  )

  ;; find denominators for s and t in the parametric plane equation
  ;; V(s,t)=Vo+su+tv. Check for non-zero values before proceding.
  (cond
    ((and u v n
          (not (equal (setq sden (dot u (cross n v))) 0.0 1E-6))
          (not (equal (setq tden (dot v (cross n u))) 0.0 1E-6))
     )

     ;; calculate s and t in the parametric equation
     (setq snum (/ (dot w (cross n v)) sden)
           tnum (/ (dot w (cross n u)) tden)
     )

     ;; is p1-p2 parallel to plane?
     (and (equal (dot n w) 0.0 1e-6)
          (princ "\nLine is parallel with plane")
     )

     ;; get unit vector p0 to p1
     (setq p01 (mapcar '- p1 p0))
     ;; check denominator for intersection parameter si
     (cond
       ((not (zerop (setq sn (- (dot n p01))
                          sd (dot n w)
                    )
             )
        )
        ;; calculate intersection parameter si
        (setq si (/ sn sd))
        ;; get intersection point pi based on p1, unit vector w and si
        (setq pint (getInters p1 w si))
        ;; investigate parameter si:
        (cond ((or (equal si 0.0 1E-6)
                   (equal si 1.0 1E-6)
               )
               (princ "\nLine touches on plane")
              )
              ((< 0.0 si 1.0) (princ "\nIntersection is physical"))
              (T (princ "\nIntersection is non-physical"))
        )
        (and pint
             (setq uint (getInters pint v 100.0)
                   vint (getInters pint u 100.0)
             )
             (setq a 0
                   i 0
             )
             (repeat (1- (length plst))
               (if (inters pint uint (nth i plst) (nth (1+ i) plst))
                 (setq a (1+ a))
               )
               (setq i (1+ i))
             )
             (if (or (zerop a)(zerop (rem a 2.0)))
               (princ "\nIntersection point is outside 3DFACE")
               (princ "\nIntersection point is inside 3DFACE"))
        )
       )
       ((princ "\nNo intersection"))
     )
    )
    (T (princ "\nVectors do not describe a plane"))
  )
  (terpri)
  (if pint pint (princ))
)