diff --git a/exp/stone_detect.py b/exp/stone_detect.py
--- a/exp/stone_detect.py
+++ b/exp/stone_detect.py
@@ -3,7 +3,7 @@ sys.path.append("../src")
 
 import os
 import math
-import heapq
+import random
 
 import cv2 as cv
 import numpy as np
@@ -14,6 +14,8 @@ from annotations import DataFile,compute
 from hough import show
 from analyzer.epoint import EPoint
 
+random.seed(361)
+
 
 class NeighboringPoint(EPoint):
 	def __init__(self,x,y):
@@ -47,7 +49,7 @@ def filterStones(contours,bwImg,stoneDim
 		center=(moments["m10"]/(moments["m00"] or 1), moments["m01"]/(moments["m00"] or 1))
 		area=cv.contourArea(c)
 		(x,y,w,h)=cv.boundingRect(c)
-		if w>stoneDims[0] or h>stoneDims[1]*1.3 or w<2 or h<2:
+		if w>stoneDims[0] or h>stoneDims[1]*1.5 or w<2 or h<2:
 			cv.drawMarker(contourImg,tuple(map(int,center)),(0,0,255),cv.MARKER_TILTED_CROSS,12)
 			keep=False
 		coverage1=area/(w*h or 1)
@@ -65,50 +67,27 @@ def filterStones(contours,bwImg,stoneDim
 	return res
 
 
-def computeDistances(points):
-	n=len(points)
-	res=np.zeros((n,n),dtype=np.float32)
-
-	for (i,a) in enumerate(points):
-		for (j,b) in enumerate(points):
-			if j<i: continue
-			elif j==i: res[i,j]=0
-			else:
-				res[i,j]=res[j,i]=a.dist(b)
-
-	return res
+def point2lineDistance(a,b,p):
+	# https://en.wikipedia.org/wiki/Point-line_distance#Line_defined_by_two_points
+	ab=b-a
+	num=abs(ab.y*p.x - ab.x*p.y + b.x*a.y - a.x*b.y)
+	denum=math.sqrt(ab.y**2+ab.x**2)
+	return num/denum # double_area / side_length == height
 
 
-def collectNeighbours(points):
-	# we could do this in O(n log n)
-	# doi:10.1007/BF02187718
-	# https://link.springer.com/content/pdf/10.1007/BF02187718.pdf
-	""":param points: [EPoint, ...]
-	:return: [(k_i,dist_i) for i in range(4)]"""
-	neighborhood=[NeighboringPoint(p.x,p.y) for p in points]
-	dists=computeDistances(points)
-	for (i,a) in enumerate(points):
-		neighbours=heapq.nsmallest(4,enumerate(dists[i]),key=lambda x: x[1] if x[1]>0 else 1000)
-		neighborhood[i].neighbours=[neighborhood[j] for (j,d) in neighbours]
-	return neighborhood
-
-
-def growLine(r,s):
-	""":param r: NeighboringPoint
-	:param s: NeighboringPoint
-	:return: NeighboringPoint or None"""
-	u=s-r
-	alpha=math.atan2(u.y,u.x)
-	t=None
-	beta=alpha+math.pi # -alpha
-	for t_ in s.neighbours:
-		v=t_-s
-		beta_=math.atan2(v.y,v.x) # beta_ in (-pi, pi]
-		if abs(abs(alpha-beta_)-math.pi)>abs(abs(alpha-beta)-math.pi):
-			beta=beta_
-			t=t_
-	if abs(alpha-beta)<math.pi/12: return t
-	else: return None
+def groupLines(points,minCount,tolerance):
+	random.shuffle(points)
+	sample=points[:57]
+	for (i,a) in enumerate(sample):
+		for (j,b) in enumerate(sample):
+			if j<=i: continue
+			incidentPoints=[a,b]
+			for c in points:
+				if c is a or c is b: continue
+				if point2lineDistance(a,b,c)<=tolerance:
+					incidentPoints.append(c)
+			if len(incidentPoints)>=minCount:
+				yield incidentPoints
 
 
 if __name__=="__main__":
@@ -128,11 +107,12 @@ if __name__=="__main__":
 
 	rect=img[y1:y2,x1:x2]
 	unit=np.array([1,1,1],dtype=np.uint8)
+	kernel=np.ones((3,3),np.uint8)
 	maskB=cv.inRange(rect,centers[0]-unit,centers[0]+unit)
-	maskB=cv.dilate(maskB,np.ones((3,3),np.uint8),iterations=1)
-	maskB=cv.erode(maskB,np.ones((3,3),np.uint8),iterations=3)
+	maskB=cv.morphologyEx(maskB,cv.MORPH_OPEN,kernel,iterations=1)
+	maskB=cv.erode(maskB,kernel,iterations=2)
 	maskW=cv.inRange(rect,centers[1]-unit,centers[1]+unit)
-	maskW=cv.erode(maskW,np.ones((3,3),np.uint8),iterations=2)
+	maskW=cv.erode(maskW,kernel,iterations=2)
 
 	show(img,filename)
 	show(maskB,filename)
@@ -145,24 +125,32 @@ if __name__=="__main__":
 
 	(contours,hierarchy)=cv.findContours(stones,cv.RETR_LIST,cv.CHAIN_APPROX_SIMPLE)
 	stoneLocs=filterStones(contours,stones,stoneDims)
-	neighborhood=collectNeighbours([point for (point,contour) in stoneLocs])
 
 	linesImg=cv.cvtColor(np.zeros((h,w),np.uint8),cv.COLOR_GRAY2BGR)
 	cv.drawContours(linesImg,[c for (point,c) in stoneLocs],-1,(255,255,255),-1)
-	for p in neighborhood:
+	for (p,c) in stoneLocs:
 		cv.drawMarker(linesImg,(int(p.x),int(p.y)),(255,0,0),cv.MARKER_CROSS)
 
-	for p in neighborhood:
-		for q in p.neighbours:
-			length=1
-			while True:
-				q_=growLine(p,q)
-				if q_:
-					q=q_
-					length+=1
-				else: break
-			if length>1 or True:
-				print(length,p,q)
-			cv.line(linesImg,(int(p.x),int(p.y)),(int(q.x),int(q.y)),(255,255,0),min((length+1)//2,3))
+	lineSet=set()
+	minCount=min(max(math.sqrt(len(stoneLocs))-2,3),7)
+	print("min count:",minCount)
+	for points in groupLines([point for (point,contour) in stoneLocs],minCount,2):
+		points=tuple(sorted(tuple(p) for p in points))
+		lineSet.add(points)
+		obsolete=set()
+		pointSet=set(points)
+		for line in lineSet:
+			lineS=set(line)
+			if points is line: continue
+			if pointSet<lineS:
+				lineSet.remove(points)
+				break
+			if lineS<pointSet:
+				obsolete.add(line)
+		lineSet-=obsolete
+	for line in sorted(lineSet,key=len,reverse=True)[:16]:
+		print(len(line),line)
+		(xa,ya)=line[0]
+		(xb,yb)=line[-1]
+		cv.line(linesImg,(int(xa),int(ya)),((int(xb),int(yb))),(255,255,0),1)
 	show(linesImg)
-	cv.imwrite("/tmp/img.png",linesImg)