Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Close

Visual LISP, AutoLISP and General Customization

- Autodesk Community
- >
- AutoCAD Customization
- >
- Visual LISP, AutoLISP and General Customization
- >
- Steiner Tree Problem

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic to the Top
- Bookmark
- Subscribe
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

619 Views, 11 Replies

11-12-2012 08:59 AM

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.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

11-12-2012 09:42 AM in reply to:
owenp

owenp 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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

11-12-2012 11:43 AM in reply to:
Kent1Cooper

Kent1Cooper wrote:... You don't need to figure

all threeLines and find the intersection ofthreeof them -- anytwowill 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 twovirtuallines.

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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

11-13-2012 05:48 AM 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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

11-13-2012 06:18 AM 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.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

11-13-2012 12:32 PM in reply to:
owenp

owenp 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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

11-14-2012 04:46 AM 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.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

11-14-2012 08:00 AM 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.h

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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

11-14-2012 08:32 AM in reply to:
owenp

owenp 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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

11-21-2012 12:51 PM in reply to:
Kent1Cooper

Kent1Cooper wrote:.... Thatdoesn'ttest for whether any internal angle of the triangle between the three points isequal to or greater than 120 degrees, in which case the Fermat point for those three points isatthat 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

Log into access your profile, ask and answer questions, share ideas and more. Haven't signed up yet? Register

Announcements

Start with some of our most frequented solutions to get help installing your software.

Upgrading to a 2015 product? Make sure to check these out 1st!

- Privacy | Legal Notices & Trademarks | Report Noncompliance | Site map | © Copyright 2014 Autodesk Inc. All rights reserved

Except where otherwise noted, this work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License. Please see the Autodesk Creative Commons FAQ for more information.