Topic: Triangulation

[Didge]

Hi Guys, I've been watching this site for a while now but this is my first topic.

Basically I have a triangular 3Dface aligned in an arbitary 3D plane,  corners defined as:
(X1,Y1,Z1) (X2,Y2,Z2) & (X3,Y3,Z3)

I would like to know what formula I should be using to obtain a Z elevation for a given point (X,Y) picked somewhere within the triangle.

Any suggestions would be most welcome.

Didge.

[CAB]

Glad to see you came out of the closet.
Your first one is a doozy.

I don't know the answer.

To clarify you pick a point in the current UCS & then want to project up or down
to a point that intersects the plane of the 3d triangle and return that point
with the resulting z value.

You may have to wait till one of the math wiz's wakes up.

[Didge]

Many thanks CAB, yes you defined the problem better than I did

I've trawled through many geometric web sites and discovered little more than my own ignorance.

[Jürg Menzi]

Hi Didge

Grab my matrix function library from HERE
and check the sample.
A solution can be:
- Select the triangle object
- Extract the points from object
- Create an infinite plane (setq PlaneList (MePlane3p P1 P2 P3))
- Select the desired point
- Calculate the distance by (MePlanePerpDist PlaneList Pnt)

Cheers

[Didge]

Thank-you very much Jurg,

I've been to your web sites many times over the years but I've never noticed that one before.

I seem to recall an easy formula from my school days but I guess I should've been listening in those maths classes.  It'll take me the weekend to digest your matrix functions, but it sounds very promising.

Once again, many thanks to you both guys.

[Jürg Menzi]

Didge, you're welcome...

You may have a look to following possible problems:
- the triangle points should always turn in the same direction (ccw) - otherwise you get the wrong sign for distance... (see code below)
- you've to check for 'inside triangle' because the plane is infinite... (see MeInsideBPlane from the matrix library)

Code: [Select]

;
; == Function MeGetCcw
; Determines the direction of 3 points.
; Copyright:
;   ©2001 MENZI ENGINEERING GmbH, Switzerland
; Arguments [Type]:
;   Pt1 = First point [LIST]
;   Pt1 = Second point [LIST]
;   Pt3 = Third point [LIST]
; Return [Type]:
;   >  1 = ccw [INT]
;   > -1 = cw [INT]
;   >  0 = Colinear [INT]
; Notes:
;   - None
;
(defun MeGetCcw (Pt0 Pt1 Pt2 / X1_Dif X1_Sqr X2_Dif X2_Sqr Y1_Dif Y1_Sqr
                                Y2_Dif Y2_Sqr)
 (setq X1_Dif (- (car Pt1) (car Pt0))
       Y1_Dif (- (cadr Pt1) (cadr Pt0))
       X2_Dif (- (car Pt2) (car Pt0))
       Y2_Dif (- (cadr Pt2) (cadr Pt0))
       X1_Sqr (* X1_Dif X1_Dif)
       Y1_Sqr (* Y1_Dif Y1_Dif)
       X2_Sqr (* X2_Dif X2_Dif)
       Y2_Sqr (* Y2_Dif Y2_Dif)
 )
 (cond
  ((> (* X1_Dif Y2_Dif) (* Y1_Dif X2_Dif)) 1)
  ((< (* X1_Dif Y2_Dif) (* Y1_Dif X2_Dif)) -1)
  ((or (< (* X1_Dif X2_Dif) 0) (< (* Y1_Dif Y2_Dif) 0)) -1)
  ((< (+ X1_Sqr Y1_Sqr) (+ X2_Sqr Y2_Sqr)) 1)
  (0)
 )
)
Oops... have seen that one function is missing in the matrix library:

Code: [Select]

;
; == Function MeMaxInLst
; Searches for the biggest number in a list of numbers.
; Arguments [Type]:
;   Lst  = List to search [LIST]
; Return [Type]:
;   > Biggest Number [REAL] or [INT]
; Notes:
;   - None
;
(defun MeMaxInLst (Lst)
 (apply 'max Lst)
)

[Serge J. Gianolla]

Let me have a wild guess. You wanna do draping on a TIN mesh.

[CAB]

OK, How about a half baked geometry solution?
I couldn't make any headway with the formulas either.
Here is my test code.
User picks the points of the 3d triangle. p1 p2 p3
Then picks a point within the triangle. p4
Using the inters function I get the 2d intersect with p1-p4 and p2-p3 -> p1>p2p3
Then using a z value I get the 3d intersect point with p1>p2p3 and p2-p3
now we have points p1 and p1>p2p3 in 3d which pass over p4
Using the inters again to get p4,0 to p4,10 and p1 to p1>p2p3 we get the 3d point
in the triangle plane. If you get the points from objects in the drawing
you will have to use the trans on them

Code: [Select]

(defun c:test ()
  ;;  (setq p1 nil p2 nil p3 nil p4 nil)
  (or p1 (setq p1 (getpoint "\nFirst"))) ; get the points only once
  (or p2 (setq p2 (getpoint "\nSecond")))
  (or p3 (setq p3 (getpoint "\nThird")))

  (or p4 (setq p4 (getpoint "\nForth")))
  (setq p4 (list (car p4) (cadr p4))) ; make 2d

  ;;  2d intersect point of p1-p4 with p2-p3
  (setq p1>p2p3 (inters p1 p4 p2 p3 nil))
  ;;  3d intersect point of p1-p4* with p2-p3
  (setq p1>p2p3 (inters
               (list (car p1>p2p3) (cadr p1>p2p3) 0.0)
               (list (car p1>p2p3) (cadr p1>p2p3) 10.0)
               p2 p3 nil)
  )
  ;;  get the pick point interset with the plane
  (setq px (inters p1 p1>p2p3
                   (list (car p4) (cadr p4) 0.0)
                   (list (car p4) (cadr p4) 10.0)
                   nil)
  )
)

[CAB]

Here is a test version that you must select the 3d lines that form the triangle.
This is a test an no error checking.

Code: [Select]

(defun c:test3 (/ e1 e2 p1 p2 p3 p4 p5 p1>p2p3 px)
  ;;  (setq p1 nil p2 nil p3 nil p4 nil)
  (setq e1 (car (entsel "\nSelect side of triangle")))
  (setq e2 (car (entsel "\nSelect Second side of triangle")))
  ;(setq p3 (entsel "\nThird"))

  (or p4 (setq p4 (getpoint "\nSelect point in triangle")))
  (setq p4 (trans p4 1 0) ;  World
        p4 (list (car p4) (cadr p4)) ; make 2d
        p1 (cdr (assoc 10 (entget e1)))
        p2 (cdr (assoc 11 (entget e1)))
        p3 (cdr (assoc 10 (entget e2)))
        p5 (cdr (assoc 11 (entget e2))))
  ;;  get the 3 triangle points
  (cond
    ((< (distance p1 p3) 0.0001) (setq p3 p5))
    ((< (distance p2 p3) 0.0001) (setq p3 p5))
    )

  ;;  2d intersect point of p1-p4 with p2-p3
  (setq p1>p2p3 (inters p1 p4 p2 p3 nil))
  ;;  3d intersect point of p1-p4* with p2-p3
  (setq p1>p2p3 (inters
               (list (car p1>p2p3) (cadr p1>p2p3) 0.0)
               (list (car p1>p2p3) (cadr p1>p2p3) 10.0)
               p2 p3 nil)
  )
  ;;  get the pick point interset woth the plane
  (setq px (inters p1 p1>p2p3
                   (list (car p4) (cadr p4) 0.0)
                   (list (car p4) (cadr p4) 10.0)
                   nil)
  )
  (setq px (trans px 0 1))
  (command ".point" "non" px)
  px
)

[CAB]

BTW, this is the only formula I could come up with but can't seem to
figure how to apply it.

Take the equation of the plane and the equation of the line
(parametrically), and solve them as a system of equations.

For instance...

plane: 2x +4y +5z = 0
line:
x=t
y=2t
z=3t+5

then you plug in x y and z into the plane's equation, and solve for t...

2*t + 4*2*t + 5*3*t + 5*5 = 0
25*t = -25
t = -1

then, once you know t, just plug it in to the line's equation to find the
point of intersection:

x = -1
y = 2*-1 = -2
z = 3*-1 + 5 = 5 - 3 = 2

so the intersection is at <-1,-2,2>. Unless I messed something up, which is
entirely possible...

[Serge J. Gianolla]

Hi CAB,

User picks the points of the 3d triangle. p1 p2 p3
Then picks a point within the triangle. p4
Using the inters function I get the 2d intersect with p1-p4 and p2-p3 -> p1>p2p3 snip...

Typically the 3DFace will be up in the air [with Z values], and the arbitrary pt P4 picked inside the triangle will have Z=0. Therefore there is no way P1 P4 ever intersecting P2 P3. So we've got go back to drawing board

[CAB]

Dis you test it?

[Serge J. Gianolla]

No I did not.

Let say that Z for P1 is 46.6, P2 is 39.4 and for P3 43.9 and picked P4 Z to Zero [which may not be true according to the ELEV command - that is another parameter to consider] The line P1 to P4 will not intersect P2 to P3, at least not inside the said 3DFace.

Granted Didge did not give us much to go by! Is P4 user selected or like you have the fence line of a property [just one line for explanation sake] in planimetry sitting on 0 and we have a TIN mesh made of 3DFaces [could be 100's]. Now to drape the line to the mesh, you have to evaluate the endpoints to each 3DFace to see if it belongs to that triangle? Does it make sense?

[CAB]

I didn't say p1 & p4 intersect.
Here is an explanation.

[Serge J. Gianolla]

I didn't say p1 & p4 intersect.

Neither have I, Alan.

Isn't it working in this case because P1 and P4 - looking at your image - are conveniently located on 0 for Z value? What if P1 is anywhere but z=0? Can't test, away from AutoCAD!

[CAB]

OK, I moved my triangle, same result using the 'test.lsp'.
p1 is no longer at 0,0,0

p1 (34.3598 13.5639 10.0)
p2 (44.3598 23.5639 20.0)
p3 (46.3598 16.5639 14.0)
p4 (39.8598 16.8139)

p1>p2p3 [2d] (45.3598 20.0639)
p1>p2p3 [3d] (45.3598 20.0639 17.0)

px  (39.8598 16.8139 13.5)

I could be missing something as I'm not a 3d guy.
I tested the test3.lsp in a rotated arouns z  UCS & it workd fine but did not test any other rotation.

[Serge J. Gianolla]

You keep working with the assumption that P4 is user-selected. If it is true then, isn't it easier to add the Z's for P1, P2 and P3 then average them:
(Z1+Z2+Z3)/3 then the obtained value is replacing Z4?

[CarlB]

I didn't say p1 & p4 intersect.
Here is an explanation.

CAB, that's a nice illustration.  I think your theory is sound.  It doesn't matter that Pt1 is at 0,0,0 that is just to make it easier to follow - the math uses elevation differences anyway.  As long as you followed a sign convention it seems you should get the right answer.  And if it turns out that the z of the pick point is outside the range of the P1>Px line, then you could also tell whether the point falls inside the triangle.  Much easier to follow than the matrix math

I don't know about your line-plane intersection example, it may be correct but I think several steps/calculations are needed to go from a 3d face and a point to get those equations.

[CAB]

Serge
I don't see why the p4, picked or sent by another routine has any bearing but I'm not that good with 3d stuff.
The average z would only apply to the center of the triangle. Well maybe not called the center but where all
three bisector lines intersect. Otherwise the average z is not correct.

[Serge J. Gianolla]

OOps

The average z would only apply to the center of the triangle.

 

[Serge J. Gianolla]

Yes, it does matter if the pt is external to 3DFace. 1st case scenario, P1 P4 could be parallel to P2 P3 never intersecting. Even though if could be a fluke, it may happen. Another case is P4 anywhere outside the 3DFace when intersecting the lines the new point projected ends up to another 3DFace if we speak of a Triangular Irregular Network mesh, thereabout giving the wrong elevation.

[CAB]

After further testing, the routine only works if you rotate about the z axis.
Rotation about the x or y axis don't work.
Perhaps that is what Serge was trying to tell me all along.
If the triangle was 5 feet in the air and you were 10 feet above that looking straight down then
any point on the ground would be in the triangle. But when you rotate the x and or y axis and you are
looking down but at an angle to the xy and the point you are looking at on the ground may not be under the
triangle relative to the z axis. Not sure if may description was accurate.
This is why I shy away from the 3d stuff.

[CAB]

You can see where he triangle projects into the xy pland along the z axis.
So "dropping" a point through what appears to be  the triangle will nit the xy plane outside the triangle.

[CAB]

Another formula if anyone wants to take a crack at it.
http://astronomy.swin.edu.au/~pbourke/geometry/planeline/

[CarlB]

CAB,

I'm not sure why you say it "didn't work".  Your routine is designed to project a point onto a plane, projecting along the z axis.  So it won't appear to look right when your view isn't perpendicular to the z axis.  Your solution, I think, solved what the original poster was trying to do.  Iintersecting an arbitrary line with a plane is more complex and will require the matrix math solutions, then "point in polygon" algorithm to find whether the point-plane intersection is inside the polygon.

[CAB]

Thanks Carl.
It does work as long as you rotate only the z-axis.
I updated the picture. Not exactly accurate but represents the problem with other rotations.
Picking the magenta X will result in a z-axis projection outside the triangle. Yellow vector.

[Kerry]

I don't believe the original posters stated intent is decisive enough to provide a solution.

Is the requirement to "use" the < perhaps twisted > view on the triangle as a "viewport" onto the current UCS. ?

If so,
because the problem is expressed using AutoCAD parameters, I see no problem using a Geometric solution, rather than a formulae.

OP feedback is required.

[Serge J. Gianolla]

Still no news from original poster, so I shall assume that there is only one triangle and indeed the point is inside - no need to check for P4 validity.
In the old GeomCal there is the ilp function. Can be accessed in Topics AutoLISP Reference, in Cal as a header: Obtaining an intersection point

intersection of a line p1-p2 and a plane defined by p3-p4-p5

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

In our case the picked pt to find Z value would be p1 and P2 would be at any user-defined distance along Z axis, assuming or setting UCS as WCS.

[Kerry]

Hi Serge,
Yes, I looked at the CAL ILP function too. It will work, given my understanding of the OP's post.
 .... but I've been wrong about interpretation before, so I held off posting a solution.

kwb
 

[Didge]

Sorry for the delay in responding Guys, I've spent a long weekend moving home and sadly no internet connection at the new place as yet.

I wasn't expecting such a detailed response, thank-you all very much.

As for CAB's input, what can I say except that I'm still trying to take it all in, not the easiest of tasks on a return to work - thank-you very much, but please let me sleep on some of those examples. (Beginning to think I've bitten-off more than I can chew

Serge is almost correct, it's not quite draping but interpolating levels from a TIN. I've found a great triangulation lisp routine on the web, "Triangulate" written by DANIELE PIAZZA and based upon some C code by the same name, credited to PAUL BOURKE.

I've written a small lisp app that creates 3 dimensional pipe & chamber networks using ACAD's solid modelling. My employer's GIS record system provides the Pipe Invert Levels and diameters, I just need to interpolate the Chamber Cover Levels to complete a totally automatic system.

With regard to the P4 issue, it was my intention to create a (Get_Z) function that takes a numeric X,Y argument, then scans through a selection set of 3D Faces before returning an interpolated Z elevation for that X,Y point.

Hope that makes sense.

Didge.

[Didge]

I think I've found the mathmatical formula I was looking for, namely the "Plane Equation":

http://www.ems-i.com/smshelp/Data_Module/Interpolation/Linear_Interpolationsms.htm

Many thanks to all those that assisted.

[CAB]

Didge,
Hope you share the (Get_Z) function when you get it working.

[Didge]

Here's a snippet I've come up with to check the formula (its a command, not a function).
All seems ok with the interpolation, but I still need to add a check to ensure that the specified point is within the triangular boudary of the 3Dface.



;*********************************************************************
; GET-Z  -  Lisp command to interpolate Z elevations from a 3D face                           *
; =====     using the plane equation.                                                                   *
;                                                                                                                       *
;*********************************************************************
(defun c:GET-Z ( / ENT P1 P1 P3 X Y X1 Y1 Z1 X2 Y2 Z2 X3 Y3 Z3 A B C D Z)

         (setq ENT (car (entsel "\nSelect 3DFace ")))
         (setq P (getpoint "\nSelect location > "))

         (setq P1 (cdr (assoc 10 (entget ENT)))
                 P2 (cdr (assoc 11 (entget ENT)))
                 P3 (cdr (assoc 12 (entget ENT)))
         )
         (setq  X (car P)
                  Y (cadr P)
                 X1 (car P1)
                 Y1 (cadr P1)
                 Z1 (caddr P1)
                 X2 (car P2)
                 Y2 (cadr P2)
                 Z2 (caddr P2)
                 X3 (car P3)
                 Y3 (cadr P3)
                 Z3 (caddr P3)
         )
 
; Add check here to determine if the point (P) is within the triangular
; 3dface, if so process the following code.
 
         (setq A (+ (* Y1 (- Z2 Z3)) (* Y2 (- Z3 Z1)) (* Y3 (- Z1 Z2)) ))
         (setq B (+ (* Z1 (- X2 X3)) (* Z2 (- X3 X1)) (* Z3 (- X1 X2)) ))
         (setq C (+ (* X1 (- Y2 Y3)) (* X2 (- Y3 Y1)) (* X3 (- Y1 Y2)) ))
         (setq D (-  (- (* A X1)) (* B Y1) (* C Z1) ))
         (setq Z (-  (- (* (/ A C) X))  (* (/ B C) Y)  (/ D C)  ))
         (prompt (strcat "\nElevation= " (rtos Z 2 3) "              "))

  (princ)
)

[CAB]

Didge,
Nice work.

Here is a quickie that uses Doug Boards function to test for which side.
I wrapped it & modified the return value.

Code: [Select]

;;  return t if p is inside of 3 points
(defun inside (p p1 p2 p3 / sideof)

  ;;D. C. Broad, Jr.
  ;;(sideof <ray-origin> <another-point-on-ray> <point-to-be-tested>)
  (defun sideof (p1 p2 p / r)
    (setq r
           (cond
             ((equal p1 p 1e-10) 0)
             (t (sin (- (angle p1 p) (angle p1 p2))))
           )
    )
    (if (equal r 0 1e-10)
      0
      (if (minusp r) -1 1) ; CAB modified
    )
  )
  ;;return values
  ;;negative = point is to the right side of the ray
  ;;0 = point is on the ray
  ;;otherwise point is on the left side of the ray.
  ;;P1 should not equal P2 for meaningful results.

  (defun 2d (pt) (list (car pt) (cadr pt)))

  (setq p  (2d p)
        p1 (2d p1)
        p2 (2d p2)
        p3 (2d p3)
  )
  (= (sideof p1 p2 p) (sideof p2 p3 p) (sideof p3 p1 p))
)

[Didge]

That works like a dream, thanks CAB, and thanks to Doug Boards for sharing that function.

Just for the challenge, I was hoping to stick with pure mathematics using Ceva's theorem;

http://en.wikipedia.org/wiki/Ceva%27s_theorem

Unfortunately I failed to get that formula working last night (probably too much whiskey on my part).

It should return "1" if and only if the point is within the triangle, my attempt returned "1" regardless of the points location.

[Jürg Menzi]

With less Single Malt  :

Code: [Select]

;
; -- MeIsInTriangle
; Checks a point for inside a triangle.
; Arguments [Typ]:
;   Pnt = Point to test [LIST]
;   Pta = First triangle point [LIST]
;   Ptb = Second triangle point [LIST]
;   Ptc = Third triangle point [LIST]
; Return [Typ]:
;   > True if inside [BOOLEAN]
;   > False if outside [BOOLEAN]
; Notes:
;   - Triangle points must be ccw
;   - Points on edge are inside
;
(defun MeIsInTriangle (Pnt Pta Ptb Ptc / CroAbp CroBcp CroCap)
 (setq CroAbp (-
               (* (- (car Ptc) (car Ptb)) (- (cadr Pnt) (cadr Ptb)))
               (* (- (cadr Ptc) (cadr Ptb)) (- (car Pnt) (car Ptb)))
              )
       CroCap (-
               (* (- (car Ptb) (car Pta)) (- (cadr Pnt) (cadr Pta)))
               (* (- (cadr Ptb) (cadr Pta)) (- (car Pnt) (car Pta)))
              )
       CroBcp (-
               (* (- (car Pta) (car Ptc)) (- (cadr Pnt) (cadr Ptc)))
               (* (- (cadr Pta) (cadr Ptc)) (- (car Pnt) (car Ptc)))
              )
 )
 (and (>= CroAbp 0) (>= CroCap 0) (>= CroBcp 0))
)

[Didge]

Well there we have it, conclusive proof that single malt is detrimental to lisp coding - that function works very well Jurg, thank-you.

I generally prefer using pure maths wherever possible as the coding usually runs cleaner and faster.
In this case however, I've tried running both of the "inside-checking" functions against a large data set and I can't find any measurable difference between them in terms of speed. They both work very well.

Thank-you both very much

[Jürg Menzi]

Glad to help...
You're right, alcohol and programming are not compatible... :doa:
In my free time I prefer scotch single malt (my favorite: Talisker from Isle of Skye)...

[Didge]

Talisker is indeed a mighty fine brew Jurg and one of my favourites too, in this case however the following was responsible. A nice peaty little number from Islay

http://www.missionliquor.com/Store/Qstore/Qstore.cgi?CMD=011&PROD=1080182500&PNAME=Laphroaig+10+Years+Original+Cask+Strength+750ml

[bullah]

please share with me,  i use this

;*********************************************************************
; GET-Z  -  Lisp command to interpolate Z elevations from a 3D face                           *
; =====     using the plane equation.                                                                   *
;                                                                                                                       *
;*********************************************************************
(defun c:GET-Z ( / ENT P1 P1 P3 X Y X1 Y1 Z1 X2 Y2 Z2 X3 Y3 Z3 A B C D Z)

         (setq ENT (car (entsel "\nSelect 3DFace ")))
         (setq P (getpoint "\nSelect location > "))

         (setq P1 (cdr (assoc 10 (entget ENT)))
                 P2 (cdr (assoc 11 (entget ENT)))
                 P3 (cdr (assoc 12 (entget ENT)))
         )
         (setq  X (car P)
                  Y (cadr P)
                 X1 (car P1)
                 Y1 (cadr P1)
                 Z1 (caddr P1)
                 X2 (car P2)
                 Y2 (cadr P2)
                 Z2 (caddr P2)
                 X3 (car P3)
                 Y3 (cadr P3)
                 Z3 (caddr P3)
         )
 
; Add check here to determine if the point (P) is within the triangular
; 3dface, if so process the following code.
 
         (setq A (+ (* Y1 (- Z2 Z3)) (* Y2 (- Z3 Z1)) (* Y3 (- Z1 Z2)) ))
         (setq B (+ (* Z1 (- X2 X3)) (* Z2 (- X3 X1)) (* Z3 (- X1 X2)) ))
         (setq C (+ (* X1 (- Y2 Y3)) (* X2 (- Y3 Y1)) (* X3 (- Y1 Y2)) ))
         (setq D (-  (- (* A X1)) (* B Y1) (* C Z1) ))
         (setq Z (-  (- (* (/ A C) X))  (* (/ B C) Y)  (/ D C)  ))
         (prompt (strcat "\nElevation= " (rtos Z 2 3) "              "))

  (princ)
)

but some time ok, some time not ok.. sorry my english is bad. Im use autocad 2002

[ymg]

Try this function:

Code - Auto/Visual Lisp: [Select]

1. ;;****************************************************************************;
2. ;; (getz p t1 t2 t3)                                                          ;
3. ;; Given point p and triangle defined by points t1, t2, t3                    ;
4. ;; Returns: (x y z) where z is on face of triangle.                           ;
5. ;;                                                                            ;
6. ;; By ymg  August 2013                                                        ;
7. ;;****************************************************************************;
8.  
9. (defun getz (p t1 t2 t3 /  v1 v2)
10.      
11.    ; Calculating a normal vector with gile's functions in line.               ;
12.    ; Note that I do not calculate the unit vector for the normal n.           ;
13.    (setq v1 (mapcar '- t2 t1)
14.          v2 (mapcar '- t3 t1)
15.           n (list (- (* (cadr v1) (caddr v2)) (* (caddr v1) (cadr v2)))
16.                   (- (* (caddr v1) (car v2)) (* (car v1) (caddr v2)))  
17.                   (- (* (car v1) (cadr v2)) (* (cadr v1) (car v2)))    
18.             )
19.    )
20.    (list (car p)(cadr p)(/ (apply '+ (mapcar '* n (mapcar '- t1 p)))(caddr n)))
21.  )
22.  

And if you need to check if your point is in a  given  triangle:

Code - Auto/Visual Lisp: [Select]

1. (defun IsInTriangle_p (p p1 p2 p3 / x x1 x2 x3 y y1 y2 y3)
2.    (setq x  (car p)  y  (cadr p)
3.          x1 (car p1) y1 (cadr p1)
4.          x2 (car p2) y2 (cadr p2)
5.          x3 (car p3) y3 (cadr p3)
6.    )
7.    (=  (minusp (- (* (- x1 x) (- y2 y)) (* (- y1 y) (- x2 x))))
8.        (minusp (- (* (- x2 x) (- y3 y)) (* (- y2 y) (- x3 x))))
9.        (minusp (- (* (- x3 x) (- y1 y)) (* (- y3 y) (- x1 x))))
10.    )  
11. )
12.  

[bullah]

im use DTM.VLX to make 3D FACE the command is "DTM" and i use "DTMZ" command displays estimated elevations (Z-coordinates) of X-Y points (plan) picked in the area of the previously generated DTM model. but i want point xyz and text in plan. so im used

;*********************************************************************
; GET-Z  -  Lisp command to interpolate Z elevations from a 3D face                           *
; =====     using the plane equation.                                                                   *
;                                                                                                                       *
;*********************************************************************
(defun c:GETZ3 ( / ENT P1 P1 P3 X Y X1 Y1 Z1 X2 Y2 Z2 X3 Y3 Z3 A B C D Z)

         (setq ENT (car (entsel "\nSelect 3DFace ")))
         (setq P (getpoint "\nSelect location > "))

         (setq P1 (cdr (assoc 10 (entget ENT)))
                 P2 (cdr (assoc 11 (entget ENT)))
                 P3 (cdr (assoc 12 (entget ENT)))
         )
         (setq  X (car P)
                  Y (cadr P)
                 X1 (car P1)
                 Y1 (cadr P1)
                 Z1 (caddr P1)
                 X2 (car P2)
                 Y2 (cadr P2)
                 Z2 (caddr P2)
                 X3 (car P3)
                 Y3 (cadr P3)
                 Z3 (caddr P3)
         )
 
; Add check here to determine if the point (P) is within the triangular
; 3dface, if so process the following code.
 
         (setq A (+ (* Y1 (- Z2 Z3)) (* Y2 (- Z3 Z1)) (* Y3 (- Z1 Z2)) ))
         (setq B (+ (* Z1 (- X2 X3)) (* Z2 (- X3 X1)) (* Z3 (- X1 X2)) ))
         (setq C (+ (* X1 (- Y2 Y3)) (* X2 (- Y3 Y1)) (* X3 (- Y1 Y2)) ))
         (setq D (-  (- (* A X1)) (* B Y1) (* C Z1) ))
         (setq Z (-  (- (* (/ A C) X))  (* (/ B C) Y)  (/ D C)  ))
         (command "point" P)
         (command "change" "l" "" "p" "e" (rtos Z 2 3) "" )
         (command "text" "j" "bl" P "0.20" "90" (rtos z 2 3) )
         (command "change" "l" "" "p" "e" (rtos Z 2 3) "" )

  (princ)
)

;*********************************************************************

I have a problem to choose select 3dFace and point out 3dface can still have value. should have no value.  thanks 

[ymg]

Bullah,

You will always have problem trying to select a 3dFace in a mesh,
because every face is always double.

Ideally you would need to have a structure to keep all your points and 3dFaces and check
the location of the points within it.

Short of this, another technique is to select a point, with it create a fenced selection set in the +x direction,
do the same  in the -x direction.  Now the face that is common to both selection set is the face you are in.

Idea was given to me by robierzogg at Hispacad.

Following is an implementation:

Code - Auto/Visual Lisp: [Select]

1. ;; Given the ename of a 3DFACE                                                ;
2. ;; Returns List of 3 unique points forming that face                          ;
3. (defun 3df->pl (en)
4.   (defun distinct (l) (if l (cons (car l) (distinct (vl-remove (car l) l)))))
5.   (setq ent (entget en)
6.           l (distinct (mapcar '(lambda (a) (cdr (assoc a ent))) '(10 11 12 13)))
7.   )
8. )
9.  
10. ;; Given 2 List                                                               ;
11. ;; Returns List of Elements Common to Both List.                              ;
12. (defun common (l1 l2)
13.   (if l1
14.      (if (member (car l1) l2)
15.         (cons (car l1) (common (cdr l1) l2))
16.         (common (cdr l1) l2)
17.      )
18.   )
19. )
20.  
21. ;; Given a Selection Set of Entities                                          ;
22. ;; Returns List of ename                                                      ;
23. (defun ss->enl ( ss / i l )
24.     (if ss
25.        (repeat (setq i (sslength ss))
26.             (setq l (cons (ssname ss (setq i (1- i))) l))
27.        )
28.     )
29. )
30.  
31.  
32. (defun crossproduct (u v)
33.   (list
34.     (- (* (cadr u) (caddr v)) (* (cadr v) (caddr u)))
35.     (- (* (car  v) (caddr u)) (* (car  u) (caddr v)))
36.     (- (* (car  u) (cadr  v)) (* (car  v) (cadr  u)))
37.   )
38. )
39.  
40. (defun getz (p p1 p2 p3 / n)  
41.    (setq p (list (car p) (cadr p) 0.)
42.          n (crossproduct (mapcar '- p2 p1) (mapcar '- p3 p1))        
43.    )
44.    (list (car p)(cadr p)(/ (apply '+ (mapcar '* n (mapcar '- p1 p)))(caddr n)))
45. )
46.  
47. ;; Return the maximum extent along the X of Selection Set of 3DFACE           ;
48. (defun maxtr (ss /  i ent tmp bb rtn)
49.   (setq  rtn 0)      
50.   (repeat (setq i (sslength ss))
51.      (setq ent (entget (ssname ss (setq i (1- i))))
52.            tmp (mapcar '(lambda (a) (cdr (assoc a ent))) '(10 11 12 13))
53.             bb (list (apply 'mapcar (cons 'min tmp)) (apply 'mapcar (cons 'max tmp)))
54.            rtn (max rtn (- (caadr bb) (caar bb)))
55.      )
56.   )
57. )    
58.  
59. ;|                                                          
60.                   Función zpto                              
61.                                                           |;
62.  
63. (defun c:zpto (/ p enl enl1 enl2 tr)
64.   (setvar "osmode" 0)
65.   (setq maxt (maxtr (ssget "_X"  '((0 . "3DFACE")))))  
66.   (while (setq p (getpoint "\nPick a Point: "))
67.    
68.     ;buscamos la 3dface repetida
69.     (setq enl1 (ss->enl (ssget "_F" (list p (list (+ (car p) maxt) (cadr p))) '((0 . "3DFACE"))))
70.           enl2 (ss->enl (ssget "_F" (list p (list (- (car p) maxt) (cadr p))) '((0 . "3DFACE"))))
71.     )
72.     (cond
73.        ((setq enl (common enl1 enl2)) (setq tr (3df->pl (car enl)))
74.                                       (princ "\nPunto intersección: ")
75.                                       (princ (getz p (car tr) (cadr tr) (caddr tr)))
76.        )
77.        (t (princ"\nPoint is Outside of  TIN"))
78.     )
79.   )
80. )
81.  

Note that in order to keep the Fenced Selection Set as small as possible,
we first scan all the faces and limit the search distance to the biggest
x distance of the faces.

To use issue command ZPTO, then simply pick points.  No need to select
the 3DFace.

ymg

[ymg]

Bullah,

And here is a GETZ optimized for speed:

Code - Auto/Visual Lisp: [Select]

1. (defun getz (p t1 t2 t3 / n1 n2 n3 x x1 x21 x31 y y1 y21 y31 z1 z21 z31)
2.    (setq x  (car  p)  y  (cadr p)
3.          x1 (car t1) y1 (cadr t1) z1 (caddr t1)        
4.         x21 (- (car t2) x1)  y21 (- (cadr t2) y1) z21 (- (caddr t2) z1)
5.         x31 (- (car t3) x1)  y31 (- (cadr t3) y1) z31 (- (caddr t3) z1)
6.          n1 (- (* y21 z31) (* z21 y31))
7.          n2 (- (* z21 x31) (* x21 z31))  
8.          n3 (- (* x21 y31) (* y21 x31))                      
9.    )  
10.    (list x y (/ (+ (* (- x1 x) n1) (* (- y1 y) n2) (* z1 n3)) n3))
11.  )
12.  

This is around 3 times faster that the one with all the mapcar.

And similarly for IsInTriangle_p :

Code - Auto/Visual Lisp: [Select]

1. (defun isintriangle_p (p p1 p2 p3 / x x1x x2x x3x y y1y y2y y3y)
2.    (setq  x (car  p)          y (cadr  p)        
3.         x1x (- (car p1) x)  y1y (- (cadr p1) y)
4.         x2x (- (car p2) x)  y2y (- (cadr p2) y)
5.         x3x (- (car p3) x)  y3y (- (cadr p3) y)  
6.    )
7.    (=  (minusp (- (* x1x y2y) (* y1y x2x)))
8.         (minusp (- (* x2x y3y) (* y2y x3x)))
9.         (minusp (- (* x3x y1y) (* y3y x1x)))
10.    )  
11. )
12.  

Should give you approx. 15% speed improvement.

ymg

[bullah]

Thank YMG.. i try put point its ok, when i try create text for the z value its not ok please help me.. Thank you very much

;; Given the ename of a 3DFACE                                                ;
;; Returns List of 3 unique points forming that face                          ;
(defun 3df->pl (en)
  (defun distinct (l) (if l (cons (car l) (distinct (vl-remove (car l) l)))))
  (setq ent (entget en)
          l (distinct (mapcar '(lambda (a) (cdr (assoc a ent))) '(10 11 12 13)))
  )
)
 
;; Given 2 List                                                               ;
;; Returns List of Elements Common to Both List.                              ;
(defun common (l1 l2)
  (if l1
     (if (member (car l1) l2)
        (cons (car l1) (common (cdr l1) l2))
        (common (cdr l1) l2)
     )
  )
)
 
;; Given a Selection Set of Entities                                          ;
;; Returns List of ename                                                      ;
(defun ss->enl ( ss / i l )
    (if ss
       (repeat (setq i (sslength ss))
            (setq l (cons (ssname ss (setq i (1- i))) l))
       )
    )
)
 
 
(defun crossproduct (u v)
  (list
    (- (* (cadr u) (caddr v)) (* (cadr v) (caddr u)))
    (- (* (car  v) (caddr u)) (* (car  u) (caddr v)))
    (- (* (car  u) (cadr  v)) (* (car  v) (cadr  u)))
  )
)
 
(defun getz (p p1 p2 p3 / n)   
   (setq p (list (car p) (cadr p) 0.)
         n (crossproduct (mapcar '- p2 p1) (mapcar '- p3 p1))         
   )
   (list (car p)(cadr p)(/ (apply '+ (mapcar '* n (mapcar '- p1 p)))(caddr n)))
)
 
;; Return the maximum extent along the X of Selection Set of 3DFACE           ;
(defun maxtr (ss /  i ent tmp bb rtn)
  (setq  rtn 0)       
  (repeat (setq i (sslength ss))
     (setq ent (entget (ssname ss (setq i (1- i))))
           tmp (mapcar '(lambda (a) (cdr (assoc a ent))) '(10 11 12 13))
            bb (list (apply 'mapcar (cons 'min tmp)) (apply 'mapcar (cons 'max tmp)))
           rtn (max rtn (- (caadr bb) (caar bb)))
     )
  )
)   
 
;|                                                         
                  Función zpto                             
                                                          |;
 
(defun c:zpto1 (/ p enl enl1 enl2 tr bulah)
  (setvar "osmode" 0)
  (setq maxt (maxtr (ssget "_X"  '((0 . "3DFACE"))))) 
  (while (setq p (getpoint "\nPick a Point: "))
 
    ;buscamos la 3dface repetida
    (setq enl1 (ss->enl (ssget "_F" (list p (list (+ (car p) maxt) (cadr p))) '((0 . "3DFACE"))))
     enl2 (ss->enl (ssget "_F" (list p (list (- (car p) maxt) (cadr p))) '((0 . "3DFACE"))))
    )
    (cond
       ((setq enl (common enl1 enl2)) (setq tr (3df->pl (car enl)))                           
                                      (setq putpoint (getz p (car tr) (cadr tr) (caddr tr)))
                                      (command "point" putpoint )

       )

       (t (princ"\nPoint is Outside of  TIN"))
    )
  )
)

[ymg]

Bullah,

Your modified code was working.  However a point does not show much unless you
change DDPTYPE.

Here I have modified so that the points are entmake'd and kept the list in the command line:

Code - Auto/Visual Lisp: [Select]

1. ;; Given the ename of a 3DFACE                                                ;
2. ;; Returns List of 3 unique points forming that face                          ;
3. (defun 3df->pl (en)
4.   (defun distinct (l) (if l (cons (car l) (distinct (vl-remove (car l) l)))))
5.   (setq ent (entget en)
6.           l (distinct (mapcar '(lambda (a) (cdr (assoc a ent))) '(10 11 12 13)))
7.   )
8. )
9.  
10. ;; Given 2 List                                                               ;
11. ;; Returns List of Elements Common to Both List.                              ;
12. (defun common (l1 l2)
13.   (if l1
14.      (if (member (car l1) l2)
15.         (cons (car l1) (common (cdr l1) l2))
16.         (common (cdr l1) l2)
17.      )
18.   )
19. )
20.  
21. ;; Given a Selection Set of Entities                                          ;
22. ;; Returns List of ename                                                      ;
23. (defun ss->enl ( ss / i l )
24.     (if ss
25.        (repeat (setq i (sslength ss))
26.             (setq l (cons (ssname ss (setq i (1- i))) l))
27.        )
28.     )
29. )
30.  
31.  
32. (defun crossproduct (u v)
33.   (list
34.     (- (* (cadr u) (caddr v)) (* (cadr v) (caddr u)))
35.     (- (* (car  v) (caddr u)) (* (car  u) (caddr v)))
36.     (- (* (car  u) (cadr  v)) (* (car  v) (cadr  u)))
37.   )
38. )
39.  
40. (defun getz (p p1 p2 p3 / n)  
41.    (setq p (list (car p) (cadr p) 0.)
42.          n (crossproduct (mapcar '- p2 p1) (mapcar '- p3 p1))        
43.    )
44.    (list (car p)(cadr p)(/ (apply '+ (mapcar '* n (mapcar '- p1 p)))(caddr n)))
45. )
46.  
47. ;; Return the maximum extent along the X of Selection Set of 3DFACE           ;
48. (defun maxtr (ss /  i ent tmp bb rtn)
49.   (setq  rtn 0)      
50.   (repeat (setq i (sslength ss))
51.      (setq ent (entget (ssname ss (setq i (1- i))))
52.            tmp (mapcar '(lambda (a) (cdr (assoc a ent))) '(10 11 12 13))
53.             bb (list (apply 'mapcar (cons 'min tmp)) (apply 'mapcar (cons 'max tmp)))
54.            rtn (max rtn (- (caadr bb) (caar bb)))
55.      )
56.   )
57. )    
58.  
59. ;|                                                          
60.                   Función zpto                              
61.                                                           |;
62.  
63. (defun c:zpto (/ p enl enl1 enl2 tr)
64.   (setvar "osmode" 0)
65.   (setq maxt (maxtr (ssget "_X"  '((0 . "3DFACE")))))  
66.   (while (setq p (getpoint "\nPick a Point: "))
67.  
68.     ;Searching For The Repeated 3DFACE
69.     (setq enl1 (ss->enl (ssget "_F" (list p (list (+ (car p) maxt) (cadr p))) '((0 . "3DFACE"))))
70.           enl2 (ss->enl (ssget "_F" (list p (list (- (car p) maxt) (cadr p))) '((0 . "3DFACE"))))
71.     )
72.     (cond
73.        ((setq enl (common enl1 enl2)) (setq tr (3df->pl (car enl))
74.                                              p (getz p (car tr) (cadr tr) (caddr tr))
75.                                       )      
76.                                       (entmakex (list '(0 . "POINT") (cons 10 p)))
77.                                       (princ "\nPoint on 3dFace: ")
78.                                       (princ p)
79.        )
80.        (t (princ"\nPoint is Outside of  TIN"))
81.     )
82.   )
83. )
84.  

[bullah]

"terima kasih tuan"  thank you sir 

[ymg]

Terima kasih kembali, Bullah