Hi guys.
To avoid checking each point against any other, I thought to use a grid of ‘quads’, that is a rectangular area with its list of points.
And then check distances inside each quad and also from the points into a quad to the points of adjacent quads.
… Don’t know if the logic makes sense, and also if the script follows this logic without errors ( I guess not … it has been written pretty quickly … )
The ‘split’ option is the number of quads used on the larger side of the bounding box of the points.
The script prints an alert message if the quads are too small and no distance gets calculated.
It also tries to use Rhino libraries for input/output only and not for data structure and calculations.
Anyway:
import Rhino
import math
def main():
# input
split = 10
while True:
gob = Rhino.Input.Custom.GetObject()
gob.AcceptNothing( True )
gob.SetCommandPrompt( 'Points ?' )
gob.GeometryFilter = Rhino.DocObjects.ObjectType.Point
spind, spopt = gob.AddOptionInteger( 'Split',
Rhino.Input.Custom.OptionInteger( split ) )
gob.GetMultiple( 0, 0 )
if gob.CommandResult() == Rhino.Commands.Result.Cancel:
return
res = gob.Result()
if res == Rhino.Input.GetResult.Option:
opind = gob.OptionIndex()
if opind == spind:
split = spopt.CurrentValue
elif res == Rhino.Input.GetResult.Object:
obrefs = gob.Objects()
break
else:
return
# turn points into lists [X,Y]
plist = []
for obref in obrefs:
pt = obref.Point().Location
plist.append( [ pt.X, pt.Y ] )
print '%d points' % len( plist )
# get bounding box: points min & max
mi = plist[ 0 ][ : ]
ma = plist[ 0 ][ : ]
for pt in plist[ 1 : ]:
if pt[ 0 ] < mi[ 0 ]:
mi[ 0 ] = pt[ 0 ]
if pt[ 0 ] > ma[ 0 ]:
ma[ 0 ] = pt[ 0 ]
if pt[ 1 ] < mi[ 1 ]:
mi[ 1 ] = pt[ 1 ]
if pt[ 1 ] > ma[ 1 ]:
ma[ 1 ] = pt[ 1 ]
# bbox size
lenx = ma[ 0 ] - mi[ 0 ]
leny = ma[ 1 ] - mi[ 1 ]
# size of 'quads' containing point lists
lenma = max( lenx, leny )
stp = lenma / split
xcnt = int( round( lenx / stp ) )
ycnt = int( round( leny / stp ) )
stx = lenx / xcnt
sty = leny / ycnt
# Rhino objects drawn for user info
ids = []
pla = Rhino.Geometry.Plane.WorldXY
p0 = Rhino.Geometry.Point3d( mi[ 0 ], mi[ 1 ], 0 )
p1 = Rhino.Geometry.Point3d( ma[ 0 ], ma[ 1 ], 0 )
rec = Rhino.Geometry.Rectangle3d( pla, p0, p1 )
ids.append( Rhino.RhinoDoc.ActiveDoc.Objects.AddCurve(
rec.ToNurbsCurve() ) )
Rhino.RhinoDoc.ActiveDoc.Views.Redraw()
interx = Rhino.Geometry.Interval( 0, stx )
intery = Rhino.Geometry.Interval( 0, sty )
# pla = Rhino.Geometry.Plane.WorldXY
# pla.OriginX = mi[ 0 ]
# for xind in range( xcnt ):
# pla.OriginY = mi[ 1 ]
# for yind in range( ycnt ):
# rec = Rhino.Geometry.Rectangle3d( pla, interx, intery )
# ids.append( Rhino.RhinoDoc.ActiveDoc.Objects.AddCurve(
# rec.ToNurbsCurve() ) )
# pla.OriginY += sty
# pla.OriginX += stx
# Rhino.RhinoDoc.ActiveDoc.Views.Redraw()
# initialize quads grid
qad = [ [ [] for ind in range( ycnt ) ] for jnd in range( xcnt ) ]
# storing point into 'quads'
print '%d X %d = %d quads will be used' % ( xcnt * ycnt, xcnt, ycnt )
for pt in plist:
qx = int( ( pt[ 0 ] - mi[ 0 ] ) / stx - 0.001 )
qy = int( ( pt[ 1 ] - mi[ 1 ] ) / sty - 0.001 )
qad[ qx ][ qy ].append( pt )
for ix in range( xcnt ):
for iy in range( ycnt ):
print '%d points in quad %d,%d' % ( len( qad[ ix ][ iy ] ), ix, iy )
# tex = '%d,%d' % ( xind, yind )
# pnt = Rhino.Geometry.Point3d( pt[ 0 ], pt[ 1 ], 0 )
# ids.append( Rhino.RhinoDoc.ActiveDoc.Objects.AddTextDot( tex, pnt ) )
# Rhino.RhinoDoc.ActiveDoc.Views.Redraw()
# initializing values
mindis = lenx + leny
pointa = []
pointb = []
gid = None
pla = Rhino.Geometry.Plane.WorldXY
# main loop on 'quads'
for ix in range( xcnt ):
for iy in range( ycnt ):
pla.OriginX = mi[ 0 ] + stx * ix
pla.OriginY = mi[ 1 ] + sty * iy
# drawing
if gid:
Rhino.RhinoDoc.ActiveDoc.Objects.Delete( gid, True )
rec = Rhino.Geometry.Rectangle3d( pla, interx, intery )
gid = Rhino.RhinoDoc.ActiveDoc.Objects.AddCurve(
rec.ToNurbsCurve() )
Rhino.RhinoDoc.ActiveDoc.Views.Redraw()
# any point in this quad ?
if not qad[ ix ][ iy ]:
print 'skipped quad %d,%d' % ( ix, iy )
continue
# points in this quad
pts = qad[ ix ][ iy ]
# calc in locad quad
print 'checking local points in quad %d,%d' % ( ix, iy )
while len( pts ) > 1:
p0 = pts.pop()
for pt in pts:
dx = abs( p0[ 0 ] - pt[ 0 ] )
if dx > mindis:
continue
dy = abs( p0[ 1 ] - pt[ 1 ] )
if dy > mindis:
continue
dis = math.sqrt( dx * dx + dy * dy )
if dis < mindis:
mindis = dis
pointa = p0
pointb = pt
if pointa:
print 'min. distance so far: %.4f' % mindis
# indices of near quads to check
x0 = ix
x1 = ix + 1 if ix < ( xcnt -1 ) else ix
y0 = iy
y1 = iy + 1 if iy < ( ycnt -1 ) else iy
# calc against near quads
pts = qad[ ix ][ iy ]
for qx in range( x0, x1 + 1 ):
for qy in range( y0, y1 + 1 ):
if ( ix == qx ) and ( iy == qy ):
continue
print 'checking quad %d,%d against quad %d,%d' % ( ix, iy, qx, qy )
qpts = qad[ qx ][ qy ]
if not qpts:
continue
for pt in pts:
for qp in qpts:
dx = abs( qp[ 0 ] - pt[ 0 ] )
if dx > mindis:
continue
dy = abs( qp[ 1 ] - pt[ 1 ] )
if dy > mindis:
continue
dis = math.sqrt( dx * dx + dy * dy )
if dis < mindis:
mindis = dis
pointa = qp
pointb = pt
if pointa:
print 'min. distance so far: %.4f' % mindis
if gid:
Rhino.RhinoDoc.ActiveDoc.Objects.Delete( gid, True )
if ids:
for gid in ids:
Rhino.RhinoDoc.ActiveDoc.Objects.Delete( gid, True )
if not pointa:
print 'No distance calculated. Please try with a lower split factor'
return
Rhino.RhinoDoc.ActiveDoc.Objects.AddLine(
Rhino.Geometry.Point3d( pointa[ 0 ], pointa[ 1 ], 0 ),
Rhino.Geometry.Point3d( pointb[ 0 ], pointb[ 1 ], 0 ) )
tex = '%.4f' % mindis
pnt = Rhino.Geometry.Point3d(
( pointa[ 0 ] + pointb[ 0 ] ) / 2.0,
( pointa[ 1 ] + pointb[ 1 ] ) / 2.0,
0 )
Rhino.RhinoDoc.ActiveDoc.Objects.AddTextDot( tex, pnt )
Rhino.RhinoDoc.ActiveDoc.Views.Redraw()
main()
Cheers
Will see if I have more initiative later on today or tomorrow, and if my brain can actually comprehend something of all that… --Mitch
