Author Topic: intersection  (Read 11619 times)

0 Members and 1 Guest are viewing this topic.

Fuccaro

  • Guest
intersection
« on: April 21, 2004, 08:41:34 AM »
Hello Lispers!

I use only AutoLisp for the moment. I try to write a program to calculate the illumination on surfaces.
I search for the intersection point between a 3dface and a ray (line) defined by two points. This involves a lot of mathematics: first I have to calculate the intersection between the line and the plane passing trough the 3 corners of the 3Dface, than I must check if that point is inside of the triangle (face).
It can be done, but probable you, the VL wizards, just take out from the hat a magic formula, something like vla_get_intersection_between_that_and_this...

hendie

  • Guest
intersection
« Reply #1 on: April 21, 2004, 09:17:09 AM »
hi Fuccaro
do you mean something like
Code: [Select]
(setq P1 (vla-intersectwith Object1 Object2 acextendnone))

 ~ Object1 & Object2 need to be valid VLA-objects, so you use something like
Code: [Select]
(setq Object1 (vlax-ename->vla-object (entlast)))

but I'm not sure how that would apply in your case since you are talking about the intersection of a line and a planar region. According to the help file it works with all objects except Except Pviewport and PolygonMesh.
If need be, I suppose you could always create a temporary "region" to test on.
If the two objects do not intersect, no data is returned.

I'm sure some of our more esteemed colleagues will have an answer for you

hendie

  • Guest
intersection
« Reply #2 on: April 21, 2004, 10:03:04 AM »
okay, I could be wrong.
I can't get it to work. I have my VLA object line and I have my VLA object region and I know for a fact that they do intersect.
I get a variant array returned but the value is always Nil.

bloody VL-alchemy  :genie:

anyone tell me where I'm going wrong ?


*edited* ~ tried it with a 3DFace also but still get the same error

hendie

  • Guest
intersection
« Reply #3 on: April 21, 2004, 10:35:39 AM »
Code: [Select]

_$ (setq MyLine (car (entsel "\nPlease select LINE : ")))
(setq Vline (vlax-ename->vla-object MyLine))

<Entity name: 4002cf78>
#<VLA-OBJECT IAcadLine 028db734>

_$ (setq MyRegion (car (entsel "\nPlease select REGION : ")))
(setq Vreg (vlax-ename->vla-object MyRegion))

<Entity name: 4002cf60>
#<VLA-OBJECT IAcadRegion 028d9664>

_$ (setq P1 (vla-intersectwith Vline Vreg acextendnone))

#<variant 8197 ...>



and when I do an inspect on p1 I get:
Code: [Select]

<type> Array of double
<value> # <safearray>


and the value is NIL


I even tried extruding the region but got an ; error: Automation Error. Not implemented yet message



and it all seemed so simple  :crazy:

hendie

  • Guest
intersection
« Reply #4 on: April 21, 2004, 10:51:16 AM »
I did a google search and found a few other Vlisps for the intersectwith function. All of them returned a variant array of Nil.
perhaps the object types are incompatible for this method.

Mark

  • Custom Title
  • Seagull
  • Posts: 28753
intersection
« Reply #5 on: April 21, 2004, 10:58:33 AM »
What do you get when you do this:
Code: [Select]
(vlax-safearray->list  (vlax-variant-value P1))
TheSwamp.org  (serving the CAD community since 2003)

Mark

  • Custom Title
  • Seagull
  • Posts: 28753
intersection
« Reply #6 on: April 21, 2004, 11:02:33 AM »
how about a dwg we can work with. upload to
http://theswamp.org/lilly.pond/
TheSwamp.org  (serving the CAD community since 2003)

hendie

  • Guest
intersection
« Reply #7 on: April 21, 2004, 11:04:49 AM »
Quote from: Mark Thomas
What do you get when you do this:
Code: [Select]
(vlax-safearray->list  (vlax-variant-value P1))


Code: [Select]

$ (vlax-safearray->list  (vlax-variant-value P1))
; error: ActiveX Server returned an error: Invalid index



 :horror:

DEVITG

  • Bull Frog
  • Posts: 479
finding intersections even with out lisp
« Reply #8 on: April 21, 2004, 10:08:52 PM »
On ACAD there is a marvelous tool

geomcal .
Quote


Command line: cal (or 'cal for transparent use)

CAL is an online geometry calculator that evaluates point (vector), real, or integer expressions. The expressions can access existing geometry using the object snap functions such as CEN, END, and INS. You can insert AutoLISP® variables into the arithmetic expression and assign the value of the expression back to an AutoLISP variable. You can use these arithmetic and vector expressions in any AutoCAD® command that expects points, vectors, or numbers.

CAL topics:

Understanding Syntax of Expressions
Formatting Feet and Inches
Formatting Angles
Using Points and Vectors
Using AutoLISP Variables
Using AutoCAD System Variables
Converting Units of Measurement
Using Standard Numeric Functions
Calculating a Vector from Two Points
Calculating the Length of a Vector
Obtaining a Point by Cursor
Obtaining the Last-Specified Point
Using AutoCAD Snap Modes in Arithmetic Expressions
Converting Points between UCS and WCS
Calculating a Point on a Line
Rotating a Point About an Axis
Obtaining an Intersection Point
Calculating a Distance
Obtaining a Radius
Obtaining an Angle
Calculating a Normal Vector
Using Shortcut Functions

It has a lot of tools and one just fited to yuor need
Quote
The ill and ilp functions determine intersection points.

ill(p1,p2,p3,p4)

Determines the intersection point between two lines (p1,p2) and (p3,p4). AutoCAD considers all points three-dimensional.

ilp(p1,p2,p3,p4,p5)

Determines the intersection point between a line (p1,p2) and a plane passing through three points (p3,p4,p5).

 


It can be used on Lisp

Code: [Select]


 (Setq ach (c:cal "pld(pa1,pa2,dst)")) ; locate a point on line [pa1 pa2] at dist

(SETQ ANPab (C:CAL "ang(INAB,pA2,pb2)")) ; calc a angle between two line's END point ant it's apex or vertice, Does not matter wher it is located






Quote
dist(p1,p2)

Determines the distance between two points, p1 and p2. This is the same as the vector expression abs(p1-p2).

dpl(p,p1,p2)

Determines the shortest distance between point p and the line passing through points p1 and p2.

 

dpp(p,p1,p2,p3)

Determines the distance from a point p to a plane defined by three points (p1,p2,p3).

dist(p1,p2)

Determines the distance between two points p1 and p2. This is the same as the vector expression abs(p1-p2).

 

The following example returns half the distance between the centers of two selected objects:

dist(cen,cen)/2

The following example finds the distance between the point 3,2,4 and a plane you define by selecting three endpoints:

dpp([3,2,4],end, end, end)





Really I do not know why it is not used more often.
Location @ Córdoba Argentina Using ACAD 2019  at Window 10

Fuccaro

  • Guest
intersection
« Reply #9 on: April 22, 2004, 02:12:40 AM »
GREAT!
Thank you DEVITG.

The ILP function fit exactly to my needs. It will return the intersection between the plane and the ray (let’s call X that point). Having the intersection point it is a 2D problem to check if it is inside of the triangle ABC. With the ILL I will have the intersection point (A1) between AX and BC. If the distance AX is greater than AA1, the point is outside. Repeating 3 times this test (for each point A, B and C) I can decide if the point is inside.

They are so smart people in Argentina... Thank you again!
And thanks for the other answers too!

hendie

  • Guest
intersection
« Reply #10 on: April 22, 2004, 03:24:36 AM »
brilliant ! glad you got it sorted Fuccaro, and thanks to Devitg for the solution.

it's still buggin me why I couldn't get there with the vla-intersectwith method though

DEVITG

  • Bull Frog
  • Posts: 479
the invisible it is allways at our reach
« Reply #11 on: April 22, 2004, 07:46:01 AM »
I discover the CAL funtion just browsing , I can not remember when I find
it but it is a simple tool , it can be used as a direct acad command from the command line and it can be used transparently.

just 'cal

and remember the famous acronym

KISS : KEEP IT SIMPLE SIR



 :D  :)
Location @ Córdoba Argentina Using ACAD 2019  at Window 10

hendie

  • Guest
intersection
« Reply #12 on: April 22, 2004, 08:04:04 AM »
I knew about the 'CAL tool, but I was not aware of the ILL & ILP options within it.

I must remember to store that in my brain somewhere in case I ever need it again

CAB

  • Global Moderator
  • Seagull
  • Posts: 10401
intersection
« Reply #13 on: April 22, 2004, 08:21:48 AM »
Fuccaro
You may find some use for these functions that were in my collection from the forums.

Code: [Select]
;|Here are a set of LISP functions that do the job for old-style polylines
(I did it a few years ago and it may not be perfect LISP). You will have
to modify the function PLINEVERTICES to be able to use new-style plines.
INPOLY does not calculate bulges, but gives correct results even for
self-intersecting plines. I used the fact that the sum of angles from
the checkpoint to all vertices is equal 0 (360) degrees if the
checkpoint is inside and any othre value if it is outside. The FUZZY
value is to avoid rounding errors.
Code is free to use by anybody

Tom
|;

(setq imr_fuzzy 0.0001)

(defun inpoly (pt en /)
  (if (equal 0.0 (anglesum pt (plinevertices en)) imr_fuzzy)
    nil
    t
  )
)


(defun c:inpoly (/ pt en)
  (setvar "cmdecho" 0)
  (if (= "POLYLINE" (cdr (assoc 0 (entget (setq en (car (entsel)))))))
    (while (setq pt (getpoint "\npoint to check: "))
      (if (inpoly pt en)
        (princ "INSIDE  polyline!")
        (princ "OUTSIDE polyline!")
      )
    )
  )
  (prin1)
)

(defun angd (scheitel p1 p2) (r2d (angr scheitel p1 p2)))

(defun r2d (winkel) (* (/ (float winkel) pi) 180.0))

(defun angr (scheitel p1 p2 / alf bet)
  (setq bet (angle scheitel p1)
        alf (angle scheitel p2)
        alf (- alf bet)
  )
  (if (< alf 0)
    (setq alf (+ (* 2 pi) alf))
  )
  alf
)

(defun plistisplanar (plist / isplanar)
  (setq isplanar t
        pl (cddr plist)
  )
  (while (and isplanar (< 3 (length plist)))
    (while (< 1 (length pl))
      (if (inters (car plist) (cadr plist) (car pl) (cadr pl))
        (setq isplanar nil
              pl nil
              plist nil
        )
      )
      (setq pl (cdr pl))
    )
    (if isplanar
      (setq plist (cdr plist)
            pl    (cddr plist)
      )
      (setq plist nil)
    )
  )
  isplanar
)

(defun anglesum (pt plist / as p1 p2 an1 an2 an)
  (setq as    0.0
        stp   (car plist)
        count 1
  )
  (while (< 1 (length plist))
    (setq p1    (car plist)
          p2    (cadr plist)
          plist (cdr plist)
          an    (angd pt p1 p2)
          an    (if (< 180.0 an)
                  (- an 360.0)
                  an
                )
          as    (+ as an)
    )
    (setq count (1+ count))
  )
  (setq
    an (angd pt p2 stp)
    an (if (< 180.0 an)
         (- an 360.0)
         an
       )
  )
  (+ as an)
)

(defun plinevertices (en / el vl)
  (while (= "VERTEX"
            (cdr (assoc 0
                        (setq el (entget (setq en (entnext
                                                    en
                                                  )
                                         )
                                 )
                        )
                 )
            )
         )
    (setq vl (cons (cdr (assoc 10 el)) vl))
  )
  (if (= (car vl) (car (reverse vl)))
    (reverse (cdr (reverse vl)))
    vl
  )
)




Code: [Select]
;|Re: Point inside the perimeter of a closed polyline ?
Subject: Re: Point inside the perimeter of a closed polyline ?
From: "Michael Doerr" <mdoerr@cebacus.de>
Date: Tue, 14 Dec 1999 06:47:25 +0000
Newsgroups: autodesk.autocad.customization
Hi Dominique,

the following code calculates the sum of the angles from a given point to
the points of the polyline and decides wether the point is inside or outside
in that way. 'PointInQuestion' means the point which has to be tested and
'Point_List' is the list of Points which represent the polyline. If the
function returns 'T' the given point is inside your polyline.

Greetings, Michael|;

(defun punktinpolylinie (pointinquestion
                         point_list
                         /
                        )
  (if (equal 0.0 (pipwinkelsumme pointinquestion point_list) 0.0001)
    nil
    t
  )
)

(defun pipwinkelsumme (pointinquestion           point_list
                       /            count        p1
                       p2           scheitel     winkeleins
                       winkelzwei
                      )
  (setq winkeleins 0.0
        scheitel   (car point_list)
        count      1
  )
  (while (< 1 (length point_list))
    (setq p1         (car point_list)
          p2         (cadr point_list)
          point_list (cdr point_list)
          winkelzwei (pipwinkelhilfe pointinquestion p1 p2)
          winkelzwei (if (< 180.0 winkelzwei)
                       (- winkelzwei 360.0)
                       winkelzwei
                     )
          winkeleins (+ winkeleins winkelzwei)
    )
    (setq count (1+ count))
  )
  (setq winkelzwei (pipwinkelhilfe pointinquestion p2 scheitel)
        winkelzwei (if (< 180.0 winkelzwei)
                     (- winkelzwei 360.0)
                     winkelzwei
                   )
  )
  (+ winkeleins winkelzwei)
)

(defun pipwinkelhilfe (pointinquestion p1 p2 / alpha beta)
  (setq beta  (angle pointinquestion p1)
        alpha (angle pointinquestion p2)
        alpha (- alpha beta)
  )
  (if (< alpha 0)
    (setq alpha (+ (* 2 pi) alpha))
  )
  (* (/ (float alpha) pi) 180.0)
)




Code: [Select]
;; Updated (03-07-03)
;; Given:
;;-------------------------------------------------------
;; Function to find or create an invisible ray VLA-Object
;; on Layer "0" an to create a global symbol $cv_ray.
;; It will leave the Ray object in the drawing to save time.
(defun @cv_ray ( / e ss)
  (if (setq ss (ssget "X" '((0 . "RAY")(8 . "0")(60 . 1))))
    (setq $cv_ray (vlax-ename->vla-object (ssname ss 0)))
    (setq $cv_ray
      (entmakex
        (list
         '(0 . "RAY")'(100 . "AcDbEntity")'(100 . "AcDbRay")
         '(8 . "0")'(60 . 1)'(10 0.0 0.0 0.0)'(11 1.0 0.0 0.0)
        )
      )
      $cv_ray (vlax-ename->vla-object $cv_ray)
    )
  )
)
;; And:
;;-------------------------------------------------
;; Function originated by Ken Alexander (03-05-03)
;; that is 10X faster than my @cv_parse_list,
;; to group data into triplets.
;; Thanks, Ken!
;;
(defun @cv_triple_up (old / new)
  (while
    (setq new (cons (list (car old)(cadr old)(caddr old)) new)
          old (cdddr old)
    )
  )
  (reverse new)
)
;; And:
;;-------------------------------------------------------------------
;; Function to determine if a point <PIQ> is inside a closed polyline
;; based on the number of intersections found between a ray whose
;; basepoint is the PIQ, and that PIQ is not one of the intersection
;; points.
;; Arguments:
;;   PIQ = 3D Point in WCS
;;   Outer = Outer Polyline VLA-Object
;; Returns:
;;   either T (inside) or nil (on or outside)
(defun @cv_inside (PIQ Outer / Points)
  (vl-load-com)
  (and
    (> (vl-list-length PIQ) 1)
    (vl-every 'numberp PIQ)
    (= (type Outer) 'VLA-Object)
    (vlax-property-available-p Outer 'Closed)
    (= (vla-get-closed Outer) :vlax-true)
    (or $cv_ray (@cv_ray))
    (or (vlax-put $cv_ray "Basepoint" PIQ) T)
    (setq Points (vlax-invoke Outer "IntersectWith" $cv_ray acExtendNone))
    (setq Points (@cv_triple_up Points))
    (= (rem (length Points) 2) 1)
    (not (equal PIQ (vlax-curve-getclosestpointto Object PIQ) 1e-11))
  )
)

;; Then:
;; ... define your point and Outer object ...
(setq Inside (@cv_inside PIQ Outer))



Code: [Select]
Re: Point inside the perimeter of a closed polyline ?
Subject: Re: Point inside the perimeter of a closed polyline ?
From: kboutora@francemel.com (Kamal Boutora)
Date: Thu, 16 Dec 1999 15:40:45 +0000
Newsgroups: autodesk.autocad.customization
In article <838ik5$kje17@adesknews2.autodesk.com>,
not.robertb@mwengineers.com says...
> Could you post code, or a link to some?
>

  From the faq of comp.graphics.algorithms (a very interesting faq that
  every programmer should have IMHO).

    you can look at the latest HTML version at either of these two sites:

http://www.exaflop.org/docs/cgafaq
http://www.cis.ohio-state.edu/hypertext/faq/usenet/graphics/algorithms-
faq/faq.html

    The code below is from Wm. Randolph Franklin <wrf@ecse.rpi.edu>

    int pnpoly(int npol, float *xp, float *yp, float x, float y)
    {
      int i, j, c = 0;
      for (i = 0, j = npol-1; i < npol; j = i++) {
        if ((((yp[i]<=y) && (y<yp[j])) ||
             ((yp[j]<=y) && (y<yp[i]))) &&
            (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))

          c = !c;
      }
      return c;
    }




Subject: Re: About Polyline
Author: Horst Kraemer <horst.kraemer@snafu.de>
Date Posted: Oct 2 1999 9:15:21:000AM


On Sat, 02 Oct 1999 09:34:42 GMT, philipma@ms5.hinet.net (Philip Ma)
wrote:

> Hi,
>
> Maybe the questions are stupid, but I need the answers. Thank you in
> advance!
>
> There are some questions about 2-D closed polyline as following:
>   1. How to know a polyline is crisscross (intersection itself) ?
>   2. If the polyline is not crisscross, how to know the polyline is
> clockwise or counterclockwise?
>   3. If the polyline is not crisscross, how to know a point is inside,
> outside or on the polyline?
>   4. If the polyline is not crisscross, how to know the area of the
> polyline?

In the order of ascending computational cost:

If

          P0=(x0,y0),P1=(x1,y1),....,PN=(xN,yN)=P0

is your closed polygon consisting of a sequence of N distinct points
P1,..,PN then the _oriented_ area is 1/2 of the sum of

        x_i*y_{i+1}-y_i*x_{i+1}

for i from 0 to N-1, i.e. for all adjacent pairs of vertices of the
closed polygon. This area may be positive or negative, according to
the orientation of the sequence relative to the coordinate system. In
any case its magnitude is the surface area.

If the coordinate system is oriented in a way that turning the
positive x-axis counterclockwise by 90deg will turn it into the
positive y-axis, then a positive area means that the polygon is
oriented counterclockwise. In fact "positive" doesn't mean
"counterclockwise". It just means "the same orientation as the
coordinate system". In fact in 2D geometry "clockwise" has no meaning.
If you look at a plane from the other side "clockwise" will turn into
"counterclockwise". You have to define additionally which side of the
plane is the "upside" and this definition can only be done by "looking
at it" from 3D. (Sorry for the side-step...)

There are many methods to test if a point of the plane is inside a
simple polygon. None of them is technical trivial.

The method I prefer is this one. It has been published by the german
mathematician Ahlemeyer, although he may not be the first who invented
it.  Let's assume for simplicity that the point to be tested is the
origin.

Divide the plane into 3 regions.

S0 consists only of the origin (0,0)

S+ is the "upper sector", consisting of all points with y>0 + all
points (x,0) with x>0

S- is the "lower sector", consisting of all points with y<0 + all
points (x,0) with x<0.

Every point of the plane belongs to exactly one sector.

Initialize a counter with 0.

Now make a round trip as before through all pairs of vertex points. If
any vertex is in S0 then (0,0) is a a vertex of the polygon and you
make your decision depending on the detail if an edge/vertex is
"inside" or "outside" the polygon.

Otherwise
  if two adjacent vertices (x_i,y_i),(x_{i+1},y_{i+1}) are in
  the same sector, do nothing.
  else if they are in two distinct sectors then calculate

      d =  x_i*y_{i+1) - y_i*x_{i+1)

  if d=0 then the origin is situated on the edge and you
  make your decision
  else if d>0 then add 1 to the counter
  else (d<0) subtract 1 from the counter

If no final decision could be taken in the middle then the point is
outside if the counter is ZERO, otherwise the point is outside. (Then
the counter is either 2 or -2).


Testing if a polygon is simple (not crossing itself) is probably the
most difficult task in terms of computation time. The crude and
straightforward approach is to test for each pair of edges if they
intersect or not. This will only be a problem for _very_ many
vertices.

Some time ago we hat a long thread about this problem in
sci.math.research, but nobody could come up with a _working_ algorithm
which was not O(n^2), i.e. which did need significantly less than
n^2/2 distinct tests. There were a lot of "convincing" proposals but
none of them was waterproof.


Regards
Horst
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.

Fuccaro

  • Guest
intersection
« Reply #14 on: April 22, 2004, 10:01:53 AM »
CAB
Receiving Devitg’s message I thought that I am on the right way and in short time I will finish my lisp. Now I have a lot of information to analyze... Just joking!
Before to ask help in this forum, I thought at two different solutions to find if a point X is inside of a triangle ABC. One was the method you sent me now: by summing the angles AXB+BXC+CXA. If the sum is =2*PI the point is inside. I anticipated some complications because in fact we are working in 3D and I prefer not to change the UCS.
Look at me: I write you long explanations, when in fact I wish to say you a single thing: thank you CAB!

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))
)