| Message |
Here is a randomized algorithm from de Berg et al., Computational Geometry with Applications, Springer-Verlag. Section 4.7. It consists of FOUR scripts, so watch for the script dividers and name the scripts appropriately.
The algorithm is simple but ingenious. The first script extracts the nodes from a shape and randomly permutes them. It calls the second script, "MiniDisk," to do the work.
MiniDisk grows the minimum enclosing disk point by point. The first two distinct points define a unique disk--they are the ends of its diameter. Each additional point that lies outside the current disk makes it larger. The task of enlarging the disk to include the previous points, while necessarily containing the new point on its boundary, is the job of "MiniDisk1".
MiniDisk1 begins with the circle defined by the first point in the list together with the special point which must be on the boundary. Then, in the same way MiniDisk grew the enclosing disk sequentially, MiniDisk1 also grows it, but this time with the constraint that the special point must always be on the boundary. At each step, a new point is considered. If it lies outside the current disk, then the disk must be enlarged. The enlarged disk must include the special point _and_ the new point on its boundary. The task of finding this disk is relegated to "MiniDisk2."
MiniDisk2 actually does some geometry: for each point in a list, it must compute the smallest circle containing it and two special points. That, fortunately, is straightforward.
The algorithm has expected running time O(N * lg(N)) where N is the number of vertices in the original shape.
All circles are represented as their centers and radii. This avoids the errors introduced by approximating circles as regular N-gons (ArcView represents circles as anything from a regular 16-gon to a regular 360-gon, depending on how it was constructed). |
| |
'
' Smallest enclosing disk
'
' Returns the smallest disk enclosing a shape.
'
' SELF: Shape
'
' Returns: {ptCenter, xRadius}
'
' WAH @QD
' AV 3.2a, 17 March 2002
'
' Comments: de Berg et al., Computational Geometry with
' Applications, Springer-Verlag. Section 4.7
'=============================================================='
theShape = SELF
if (theShape.Is(Polyline)) then
lPts = {}
for each lPtShape in theShape.AsList
lPts = lPts + lPtShape
end
elseif (theShape.Is(Multipoint)
or theShape.Is(GeoCurve)) then
lPts = theShape.AsList
elseif (theShape.Is(Rect)) then
ptO = theShape.ReturnOrigin
ptS = theShape.ReturnSize
ptCenter = ptS / (2@2) + ptO
xRadius = ptCenter.Distance(ptO)
return {ptCenter, xRadius}
elseif (theShape.Is(Line)) then
pt0 = theShape.ReturnStart
pt1 = theShape.ReturnEnd
ptCenter = (pt0 + pt1)/(2@2)
return {ptCenter, ptCenter.Distance(pt0)}
elseif (theShape.Is(Circle)) then
return {theShape.ReturnCenter, theShape.GetRadius}
elseif (theShape.Is(Oval)) then
ptS = theShape.ReturnSize
ptCenter = theShape.ReturnCenter
return {ptCenter, ptS.GetX max ptS.GetY}
elseif (theShape.Is(Ellipse)) then
ptCenter = theShape.ReturnCenter
ptAxes = theShape.ReturnAxes
return {ptCenter, (ptAxes.GetX max ptAxes.GetY)/2}
elseif (theShape.Is(Point)) then
return {theShape.Clone, 0}
else
return {Point.MakeNull, Number.MakeNull}
end
'
' Randomly permute lPts.
'
for each i in lPts.Count-1..1 by -1
k = Number.MakeRandom(0, i+1) mod (i+1)
ptI = lPts.Get(i)
lPts.Set(i, lPts.Get(k))
lPts.Set(k, ptI)
end
'
' Apply the algorithm.
'
return av.Run("MiniDisk", lPts)
' end of script
********
'
' MiniDisk
'
' SELF: List of points, already randomly permuted.
'
' Returns: {ptCenter, xRadius} for the
' smallest disc enclosing the points.
'
' Called by "Smallest enclosing disk".
'
' SELF must have at least two points.
'=================================================='
lPts = SELF
if (lPts.Count > 1) then
ptCenter = (lPts.Get(0) + lPts.Get(1)) / (2@2)
xRadius = ptCenter.Distance(lPts.Get(0))
lPtVertices = {lPts.Get(0), lPts.Get(1)}
if (lPts.Count > 2) then
for each i in 2..(lPts.Count-1)
ptP = lPts.Get(i)
if (ptP.Distance(ptCenter) > xRadius) then
lDisk = av.Run("MiniDisk1", {lPts, i-1, ptP})
ptCenter = lDisk.Get(0)
xRadius = lDisk.Get(1)
lPtVertices = lDisk.Get(2)
end
end
end
else
ptCenter = lPts.Get(0)
xRadius = 0
lPtVertices = lPts.Clone
end
return {ptCenter, xRadius, lPtVertices}
' end of script
********
'
' MiniDisk1
'
' SELF: {
' List of points,
' Index of last point,
' Point
' }
'
' Returns: {ptCenter, xRadius} for the
' smallest disc enclosing the points in the
' list with the point on its boundary.
'
' Called by "Minidisk".
'=================================================='
lPts = SELF.Get(0)
N = SELF.Get(1)
ptQ = SELF.Get(2)
ptCenter = (lPts.Get(0) + ptQ) / (2@2)
xRadius = ptCenter.Distance(ptQ)
lPtVertices = {lPts.Get(0), ptQ}
for each i in 1..N
ptP = lPts.Get(i)
if (ptP.Distance(ptCenter) > xRadius) then
lDisk = av.Run("MiniDisk2", {lPts, i-1, ptQ, ptP})
ptCenter = lDisk.Get(0)
xRadius = lDisk.Get(1)
lPtVertices = lDisk.Get(2)
end
end
return {ptCenter, xRadius, lPtVertices}
' end of script
********
'
' MiniDisk2
'
' SELF: {
' List of points,
' Index of last point,
' PointQ,
' PointR
' }
'
' Returns: {ptCenter, xRadius} for the
' smallest disc enclosing the points in the
' list with PointR and PointQ on its boundary.
'
' Called by "Minidisk1".
'
' Modified 8/11/02: added if (xDelta <> 0.0).
'=================================================='
if (SELF = NIL) then
lPts = {1@1}
N = 0
ptQ =0@0
ptR = 1@0
else
lPts = SELF.Get(0)
N = SELF.Get(1)
ptQ = SELF.Get(2)
ptR = SELF.Get(3)
end
ptCenter = (ptR + ptQ) / (2@2)
xRadius = ptCenter.Distance(ptQ)
ptO = 0@0
ptB = ptR - ptQ
c2 = (ptR.Distance(ptO)^2 - (ptQ.Distance(ptO)^2))/2
for each i in 0..N
ptP = lPts.Get(i)
if (ptP.Distance(ptCenter) > xRadius) then
'
' Construct the disk through ptP, ptQ, and ptR
' by intersecting perpendicular bisectors of
' two sides of the triangle.
'
' (When the points are not distinct, do a
' special case.)
'
if (0@0 = ptB) then
ptCenter = (ptP + ptQ) / (2@2)
xradius = ptQ.Distance(ptCenter)
else
ptA = ptQ - ptP
xDelta = ptA.GetX * ptB.getY - (ptA.GetY * ptB.GetX)
if (xDelta <> 0.0) then
c1 = (ptQ.Distance(ptO)^2 - (ptP.Distance(ptO)^2))/2
x = (ptB.GetY * c1 - (ptA.GetY * c2)) / xDelta
y = (ptA.GetX * c2 - (ptB.GetX * c1)) / xDelta
ptCenter = x@y
xRadius = ptCenter.Distance(ptP)
end
end
end
end
'''MsgBox.Info(ptCenter.AsString + nL + xradius.AsString,"")
return {ptCenter, xRadius, {ptP, ptQ, ptR}}
' end of script |
