Visual LISP, AutoLISP and General Customization
cancel
Showing results for 
Show  only  | Search instead for 
Did you mean: 

Steiner Tree Problem

11 REPLIES 11
Reply
Message 1 of 12
owenp
1079 Views, 11 Replies

Steiner Tree Problem

I've got a big task ahead of me. Advanced help is gratefully received.

I want to build a Steiner tree in AutoCAD, based on a set of essentially random points on a plane.

For the uninitiated, it's essentially a network that uses the shortest amount of path to connect nodes on a plane, here's an example. I'm a drainage engineer, and it'll save me a great deal of time and guesswork when creating my network.

Once we've cracked the initial Steiner tree problem, I'd like to try incorporating obstacles (eg. buildings), that the network has to avoid. I realise that this is a complex problem and lisp may not be up to the task.

To begin with, we'll need to build on finding the Fermat Point of three points, which looks like a tree with 120 degree angles.

This is my second forray into Lisp, so I'm pretty basic, but I've made a stab at the code for the Fermat Point. Let me know what you think. 

 

 

(defun C:fp ()
  (setq fpoint 0)
  (setq pta (getpoint "\nPick first point"))
  (princ pta)
  (setq ptb (getpoint "\nPick second point"))
  (princ ptb)
  (setq ptc (getpoint "\nPick third point"))
  (princ ptc)
  
(defun Line (p1 p2)
  (entmakex (list (cons 0 "LINE")
                  (cons 10 p1)
                  (cons 11 p2))))

;;Find the average point of the three vertices
(setq lstAveragePoint
       (mapcar '(lambda (X Y Z)(/ (+ X Y Z) 3.0)) pta ptb ptc))
  (print "Average")
  (prin1 lstAveragePoint)
  
;; Get properties of the triangle a is angle, d is distance
  (setq aatob (angle pta ptb))
  (setq abtoc (angle ptb ptc))
  (setq actoa (angle ptc pta))
  (setq datob (distance pta ptb))
  (setq dbtoc (distance ptb ptc))
  (setq dctoa (distance ptc pta))

;;Find outer points of the equilateral triangles
  (defun outerpoint (pt ang dist x)
	(setq outside (polar pt (+ ang (/ pi 3.0)) dist))
    	(setq inside (polar pt (- ang (/ pi 3.0)) dist))

    (if (> (distance inside lstAveragePoint) (distance outside lstAveragePoint))
        (line inside x)
      (line outside x)
     )
   )
  (outerpoint pta aatob datob ptc)
  (outerpoint ptb abtoc dbtoc pta)
  (outerpoint ptc actoa dctoa ptb)
  
;;-------------=={ Intersections in Object List }==-----------;;
;;                                                            ;;
;;  Returns a list of all points of intersection between      ;;
;;  all objects in a list of VLA-Objects.                     ;;
;;------------------------------------------------------------;;
;;  Author: Lee Mac, Copyright © 2012 - www.lee-mac.com       ;;
;;------------------------------------------------------------;;
;;  Arguments:                                                ;;
;;  lst - List of VLA-Objects                                 ;;
;;------------------------------------------------------------;;
;;  Returns:  List of intersection points, or nil             ;;
;;------------------------------------------------------------;;


(defun fpoint ( lst / a l )
    (while (setq a (car lst))
        (foreach b (setq lst (cdr lst))
            (setq l (cons (LM:Intersections a b acextendnone) l))
        )
    )
    (apply 'append (reverse l))
)

(defun LM:Intersections ( obj1 obj2 mode / l r )
    (setq l (vlax-invoke obj1 'intersectwith obj2 mode))
    (repeat (/ (length l) 3)
        (setq r (cons (list (car l) (cadr l) (caddr l)) r)
              l (cdddr l)
        )
    )
    (reverse r)
)

;;  (point fpoint)

)

 I couldn't work out how to find the intersection point of the three lines, or trim them, but I've got a copy of Lee Mac's code in there to help.

 

11 REPLIES 11
Message 2 of 12
Kent1Cooper
in reply to: owenp


@Anonymous wrote:

... we'll need to build on finding the Fermat Point of three points, which looks like a tree with 120 degree angles.

.... I couldn't work out how to find the intersection point of the three lines....

 


An interesting problem.  I haven't dug too deep, and am having computer problems at the moment, but I have some suggestions that should simplify it a lot for you:  You don't need to figure all three Lines and find the intersection of three of them -- any two will do.  And you don't need to actually draw Lines so you can find their intersections as VLA objects -- you can use an ordinary (inters) function on the endpoints of two virtual lines.

Kent Cooper, AIA
Message 3 of 12
Kent1Cooper
in reply to: Kent1Cooper


@Kent1Cooper wrote:

... You don't need to figure all three Lines and find the intersection of three of them -- any two will do.  And you don't need to actually draw Lines so you can find their intersections as VLA objects -- you can use an ordinary (inters) function on the endpoints of two virtual lines.


For instance [minimally tested]:

 

(defun C:FP (/ tript pta ptb ptc outab outac); = Fermat Point of three designated points
  (defun tript (p1 p2 p3 / ptm ptn); = TRIangle-defining PoinT
    (setq
      ptm (polar p1 (+ (angle p1 p2) (/ pi 3)) (distance p1 p2))
      ptn (polar p1 (- (angle p1 p2) (/ pi 3)) (distance p1 p2))
    ); setq
    (if (> (distance ptm p3) (distance ptn p3)) ptm ptn)
  ); defun -- tript
  (setq
    pta (getpoint "\nPick first point: ")
    ptb (getpoint "\nPick second point: ")
    ptc (getpoint "\nPick third point: ")
    outab (tript pta ptb ptc); OUTside vertex of equilateral triangle on edge AB
    outbc (tript ptb ptc pta); same for edge BC
    fp (inters outab ptc outbc pta)
  ); setq
); defun -- C:FP

 

It leaves the Fermat Point for those three points in the variable called fp.

 

EDIT:  That doesn't test for whether any internal angle of the triangle between the three points is equal to or greater than 120 degrees, in which case the Fermat point for those three points is at that vertex.  Is that something you would want checked?

 

I'm pondering the complexity of the larger concept [even without the element of something like buildings to work around]....  Would it be necessary to calculate Fermat points for every possible combination of three points, and somehow compare them?  Do you have access to the algorithm for the applet on one of your links?

Kent Cooper, AIA
Message 4 of 12
owenp
in reply to: owenp

Thanks for your help Kent,

 

I'd just discovered the "inters" function last night, and have streamlined the code (looking for the intersect of 2 lines, not 3) so it works well for almost all situations. I'm quite pleased with it so far.

 

Now the only occasion it doesn't work is when there is an angle greater than 120 degrees.

I'll add a boolean to check for that. But in most cases, it works, so now I'm wondering how to make it work for 4 or more points.

I imagine it's a case of setting one of the points from some clever maths, instead of picking it out, but I'm in over my head here, and welcome any advice.

 

(defun C:fp (/ x)
  (setq fpoint 0)
  (setq pta (getpoint "\nPick first point"))
  (setq ptb (getpoint "\nPick second point"))
  (setq ptc (getpoint "\nPick third point"))


  (defun Line (p1 p2)
    (entmakex (list (cons 0 "LINE")
		    (cons 10 p1)
		    (cons 11 p2)
	      )
    )
  )
  ;;Find the average point of the three vertices
  (setq	lstAveragePoint
	 (mapcar '(lambda (x y z) (/ (+ x y z) 3.0)) pta ptb ptc)
  )
  ;; Get properties of the triangle where a is angle, d is distance
  (setq aatob (angle pta ptb))
  (setq abtoc (angle ptb ptc))
  (setq actoa (angle ptc pta))
  (setq datob (distance pta ptb))
  (setq dbtoc (distance ptb ptc))
  (setq dctoa (distance ptc pta))

  ;;Find outer points of the equilateral triangles
  (defun outerpoint (pt ang dist x)
    (setq outside (polar pt (+ ang (/ pi 3.0)) dist))
    (setq inside (polar pt (- ang (/ pi 3.0)) dist))
    (if	(> (distance inside lstAveragePoint)
	   (distance outside lstAveragePoint)
	)
      (setq outside inside)
    )
  )
  ;;Run "outerpoint" for each side
  (outerpoint pta aatob datob ptc)
  (setq oppc outside)
  (outerpoint ptb abtoc dbtoc pta)
  (setq oppa outside)
  (outerpoint ptc actoa dctoa ptb)
  (setq oppb outside)

  ;;Perform intersection
  (setq fpoint (inters pta oppa ptb oppb))
  ;;Draw lines
  (line pta fpoint)
  (line ptb fpoint)
  (line ptc fpoint)
  ;;  (point (inters pta oppa ptc oppc))
  ;;  (point (inters ptb oppb ptc oppc))
)
;; End of lisp routine

 

Tags (1)
Message 5 of 12
owenp
in reply to: Kent1Cooper

Kent - Sorry I just fully read your reply. Your code is very compact and neat! Mine is rather verbose and cobbled together from things I've found that work.

 

I managed to track down the java code, but it's a whole other language and I'm well out of my depth. See attached, I've had to change the extension from .java to .txt, as I can't upload it otherwise.

Message 6 of 12
Kent1Cooper
in reply to: owenp


@Anonymous wrote:

... now I'm wondering how to make it work for 4 or more points.

....


I think four points would be comparatively easy, using a similar approach as with three, except that the virtual lines that the layout springs from go between the far corners of triangles built on opposite sides of the virtual quadrilateral, and not to any of the original points as is done with three.  See the attached.  Lines go from those original points at 60-degree angles to the far-triangle-point [green] lines; it results in 120-degree angles in the "trees."  There are two trees possible [at least in this layout], red and yellow in the image.  In this case, the red one has a slightly smaller total length, but I'm not sure whether you can assume anything such as that the one with the longer center segment ["trunk"] will always have the shorter total length, or maybe that triangles need to be built only on the pair of opposite sides with the smaller total length.  It may be affected by the relationship between the four points -- the actual shape, proportions, etc., of the virtual quadrilateral.

 

I can only imagine that the complexity of the problem would increase exponentially as more points are added.  But does it look like it would be worth developing code for a four-point version?

Kent Cooper, AIA
Message 7 of 12
owenp
in reply to: Kent1Cooper

Great diagram!

You're right, the problem does get exponentially harder with more points, but the hard part is in choosing which four points to connect at any time. Lets give four points a go and we can take it from there. 

 

There's an algorithm for creating the four point tree, called the Melzak-Hwang algorithm, which this article describes quite well. I can't find any code equivalent.

 

In order to find the shortest path, a brute force approach is probably necessary; creating each path, measuring the total length and picking the shortest.

Message 8 of 12
owenp
in reply to: owenp

I just wanted to give an update, because I found a java based application that will provide the coordinates of steiner points for a list of points.   http://www.flywings.org.uk/FindSteinerTree/index.htm 

It may be worth checking out for anyone literate in java. 

Message 9 of 12
Kent1Cooper
in reply to: owenp


@Anonymous wrote:

.... Lets give four points a go and we can take it from there. 

....


A start, anyway [note limitation in comment near the top]:

 

(defun C:FT4 ; = Fermat Tree of 4 designated points
  ;; works, so far, only for convex virtual quadrilaterals not too extremely proportioned
  (/ tript treept pt1 pt2 pt3 pt4 pta ptb ptc ptd outab outcd inab incd)
  (defun tript (p1 p2 p3 / ptm ptn); = TRIangle-defining PoinT
    (setq
      ptm (polar p1 (+ (angle p1 p2) (/ pi 3)) (distance p1 p2))
      ptn (polar p1 (- (angle p1 p2) (/ pi 3)) (distance p1 p2))
    ); setq
    (if (> (distance ptm p3) (distance ptn p3)) ptm ptn)
  ); defun -- tript
  (defun treept (p1 p2 p3 / ptm ptn); = TREE-intersection PoinT
    (setq
      ptm (inters p1 p2 p3 (polar p3 (+ (angle p1 p2) (/ pi 3)) 1) nil)
      ptn (inters p1 p2 p3 (polar p3 (- (angle p1 p2) (/ pi 3)) 1) nil)
    ); setq
    (if (< (distance ptm p2) (distance ptn p2)) ptm ptn)
  ); defun -- treept
  (setq
    pt1 (getpoint "\nPick first point: ")
    pt2 (getpoint "\nPick second point: ")
    pt3 (getpoint "\nPick third point continuing in same direction around: ")
    pt4 (getpoint "\nPick fourth point: ")
  ); setq
  (if ; pair of opposite sides with smaller total length
    (<
      (+ (distance pt1 pt2) (distance pt3 pt4))
      (+ (distance pt2 pt3) (distance pt4 pt1))
    ); <
    (setq pta pt1 ptb pt2 ptc pt3 ptd pt4); then
    (setq pta pt2 ptb pt3 ptc pt4 ptd pt1); else
  ); if
  (setq
    outab (tript pta ptb ptc); OUTside vertex of equilateral triangle on edge A-B
    outcd (tript ptc ptd pta); same for edge C-D
    inab (treept outab outcd pta); INside tree-intersection from vertices A & B
    incd (treept outcd outab ptc); INside tree-intersection from vertices C & D
  ); setq
  (entmake (list '(0 . "LINE") (cons 10 pta) (cons 11 inab)))
  (entmake (list '(0 . "LINE") (cons 10 ptb) (cons 11 inab)))
  (entmake (list '(0 . "LINE") (cons 10 ptc) (cons 11 incd)))
  (entmake (list '(0 . "LINE") (cons 10 ptd) (cons 11 incd)))
  (entmake (list '(0 . "LINE") (cons 10 inab) (cons 11 incd)))
  (princ)
); defun -- C:FT4

Kent Cooper, AIA
Message 10 of 12
Kent1Cooper
in reply to: Kent1Cooper


@Kent1Cooper wrote:
.... That doesn't test for whether any internal angle of the triangle between the three points is equal to or greater than 120 degrees, in which case the Fermat point for those three points is at that vertex.  Is that something you would want checked?

....


Here's a version that accounts for that, and as in the 4-point routine, draws the Fermat Tree rather than just calculating the point:

 

(defun C:FT3 (/ tript pta ptb ptc ptlist outab outac fp dlist closept other2)
  ;; = Fermat Tree of 3 designated points
  ;; takes into account possible angle >= 120 degrees
  (defun tript (p1 p2 p3 / ptm ptn); = TRIangle-defining PoinT
    (setq
      ptm (polar p1 (+ (angle p1 p2) (/ pi 3)) (distance p1 p2))
      ptn (polar p1 (- (angle p1 p2) (/ pi 3)) (distance p1 p2))
    ); setq
    (if (> (distance ptm p3) (distance ptn p3)) ptm ptn)
  ); defun -- tript
  (setq
    pta (getpoint "\nPick first point: ")
    ptb (getpoint "\nPick second point: ")
    ptc (getpoint "\nPick third point: ")
    ptlist (list pta ptb ptc)
    outab (tript pta ptb ptc); OUTside vertex of equilateral triangle on edge AB
    outbc (tript ptb ptc pta); same for edge BC
    fp (inters outab ptc outbc pta nil); = [initial] Fermat Point
  ); setq
  (foreach pt ptlist (setq dlist (cons (list (distance pt fp) pt) dlist)))
    ;; = list of distances from initial fp paired with points
  (setq
    closept (cadar (vl-sort dlist '(lambda (x y) (< (car x) (car y))))); closest vertex to initial fp
    other2 (vl-remove closept ptlist)
  ); setq
  (if (>= (distance fp (car other2)) (distance closept (car other2))); 120 degrees or more
    (setq fp closept ptlist other2)
      ;; Fermat Point AT that vertex; draw Lines only from other 2 to there
  ); if
  (foreach pt ptlist ; original 3 points or reduced to 2 by above
    (entmake (list '(0 . "LINE") (cons 10 pt) (cons 11 fp)))
  ); foreach

  (princ)
); defun -- C:FT3

Kent Cooper, AIA
Message 11 of 12
owenp
in reply to: owenp

Kent,

 

Thanks for this, I always get confused trying to use "foreach", and I still can't get my head around that if statement that checks for the greater than 120 degree angle, but I'm pleased it works. Well done.

Any chance it could be used to solve the 4 point problem, no matter where the four points are?

 

I'm so grateful for your help, I bit off more than I could chew with this project.

Message 12 of 12
Kent1Cooper
in reply to: owenp


@Anonymous wrote:

.... 

Thanks for this, I always get confused trying to use "foreach", and I still can't get my head around that if statement that checks for the greater than 120 degree angle, but I'm pleased it works. Well done.

Any chance it could be used to solve the 4 point problem, no matter where the four points are?

....


The (if) about the 120-degree-or-larger angle doesn't do it by actually checking the angle.  It just checks whether the initially-calculated triangles-on-the-outside-based Fermat point falls outside the triangle, by comparing the distance from the vertex closest to it to one of the other vertices with the distance from the calculated Fermat point to that other vertex -- if the latter is greater, it's outside, or if it's equal, it's right on that closest vertex.  [To see where that calculated Fermat point falls outside, comment out that entire (if) function with the test, so it neither moves the fp point nor removes a vertex from the list of points before drawing the tree.]

 

I have a feeling the unlimited four-point problem is going to be a far more significant challenge for AutoLISP, involving a lot of comparisons of distances.  Maybe it can be done by solving the three-point case with each possible combination of three points, and somehow comparing the results.  Or it could be by solving for four points and testing whether the resulting tree extends outside the quadrilateral, though that could be a challenge for certain configurations, e.g. something like a Y arrangement of points.  But I'll play with it.

Kent Cooper, AIA

Can't find what you're looking for? Ask the community or share your knowledge.

Post to forums  

Autodesk Design & Make Report

”Boost