diff --git a/exp/stone_detect.py b/exp/stone_detect.py
--- a/exp/stone_detect.py
+++ b/exp/stone_detect.py
@@ -12,7 +12,8 @@ import scipy.ndimage
 
 from annotations import DataFile,computeBoundingBox
 from hough import show
-from analyzer.epoint import EPoint
+from analyzer.epoint import EPoint,homogenize
+from analyzer.grid import transformPoint
 
 random.seed(361)
 
@@ -26,6 +27,31 @@ class Line():
 	def getSortedPoints(self):
 		return tuple(sorted(self.points))
 
+	def computeAngle(self,line):
+		ab=self.a-self.b
+		cd=line.a-line.b
+		alpha=math.atan(ab.y/ab.x)
+		gamma=math.atan(cd.y/cd.x)
+		fi=max(alpha,gamma)-min(alpha,gamma)
+		return min(fi,math.pi-fi)
+
+	def intersect(self,line):
+		p=self.toProjective()
+		q=line.toProjective()
+		return EPoint.fromProjective(np.cross(p,q))
+
+	def toProjective(self):
+		return homogenize(np.cross(self.a.toProjective(),self.b.toProjective()))
+
+	def transform(self,matrix):
+		a=EPoint.fromProjective(transformPoint(self.a.toProjective(),matrix))
+		b=EPoint.fromProjective(transformPoint(self.b.toProjective(),matrix))
+		if a is None or b is None: return None
+		return Line(a,b)
+
+	def __str__(self): return "({0},{1})".format(self.a,self.b)
+	def __repr__(self): return "Line({0},{1})".format(repr(self.a),repr(self.b))
+
 
 def kmeans(img):
 	arr=np.reshape(img,(-1,3)).astype(np.float)
@@ -94,11 +120,61 @@ def groupLines(points,minCount,tolerance
 				yield ab
 
 
+def computeRectiMatrix(p,q,r,s):
+	# p || q, r || s
+	vanish1=homogenize(np.cross(p.toProjective(),q.toProjective()))
+	vanish2=homogenize(np.cross(r.toProjective(),s.toProjective()))
+	horizon=homogenize(np.cross(vanish1,vanish2))
+
+	return np.matrix([horizon,[0,1,0],[0,0,1]])
+
+
+def scoreMatrix(matrix,p,r,lines):
+	p_=p.transform(matrix)
+	r_=r.transform(matrix)
+	if p_ is None or r_ is None:
+		return math.inf
+	score=0
+	for ab in lines:
+		if ab is p or ab is r: continue
+		ab_=ab.transform(matrix)
+		if ab_ is None:
+			score+=math.pi/2
+			continue
+		alpha=min(ab_.computeAngle(p_), ab_.computeAngle(r_))
+		if alpha<=math.pi/30:
+			score+=alpha
+		else: score+=math.pi/2
+	return score
+
+
+def groupParallels(lines,tolerance,w):
+	ctrl=False
+	for (i,p) in enumerate(lines):
+		for (j,q) in enumerate(lines[i+1:]):
+			a=p.intersect(q)
+			if a.y>0 and a.x<w: continue
+			for r in lines[i+1+j+1:]:
+				b=r.intersect(p) # !! ideal points
+				c=r.intersect(q)
+				ab=b-a
+				ac=c-a
+				if abs(ab.x*ac.y-ab.y*ac.x)/2<=tolerance: # area
+					yield (p,q,r)
+					ctrl=True
+					break
+			if ctrl: break
+		if ctrl:
+			ctrl=False
+			continue
+
+
 if __name__=="__main__":
 	filepath=sys.argv[1]
 	annotations=DataFile(sys.argv[2])
 	filename=os.path.basename(filepath)
-	(x1,y1,x2,y2)=computeBoundingBox(annotations[filename][0])
+	corners=annotations[filename][0]
+	(x1,y1,x2,y2)=computeBoundingBox(corners)
 	(w,h)=(x2-x1,y2-y1)
 	img=cv.imread(filepath)
 	(x3,x4,y3,y4)=(x1+w//4,x1+3*w//4,y1+h//4,y1+3*h//4)
@@ -134,9 +210,10 @@ if __name__=="__main__":
 	cv.drawContours(linesImg,[c for (point,c) in stoneLocs],-1,(255,255,255),-1)
 	for (p,c) in stoneLocs:
 		cv.drawMarker(linesImg,(int(p.x),int(p.y)),(255,0,0),cv.MARKER_CROSS)
+	matsImg=np.copy(linesImg)
 
 	lineDict=dict()
-	minCount=min(max(math.sqrt(len(stoneLocs))-2,3),7)
+	minCount=min(max(math.sqrt(len(stoneLocs))-4,3),7)
 	print("min count:",minCount)
 	for line in groupLines([point for (point,contour) in stoneLocs],minCount,2):
 		key=line.getSortedPoints()
@@ -147,15 +224,109 @@ if __name__=="__main__":
 		for ab in lineDict.values():
 			if ab is line: continue
 			if line.points<ab.points: # == impossible
-				del lineDict[line.getSortedPoints()]
+				del lineDict[key]
 				break
 			if ab.points<line.points:
 				obsolete.add(ab.getSortedPoints())
 		for key in obsolete: del lineDict[key]
 
-	for line in sorted(lineDict, 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)
+	print("valid lines:",len(lineDict))
+	lines=sorted(lineDict.values(), key=lambda ab: len(ab.points), reverse=True)
+	res=[]
+	for line in lines:
+		v=line.b-line.a
+		alpha=math.atan(v.y/v.x)
+		res.append((round(math.pi/2-alpha if alpha>0 else math.pi/2+alpha,3),repr(line)))
+		(xa,ya)=line.a
+		(xb,yb)=line.b
+		cv.line(linesImg,(int(xa),int(ya)),(int(xb),int(yb)),(255,255,0),1)
+	res.sort()
+	# for row in res: print(row)
+	# linePack=[
+	# 	Line(EPoint(174.457,166.127),EPoint(174.96,27.253)),
+	# 	Line(EPoint(191.333,38.075),EPoint(191.485,71.227)),
+	# 	Line(EPoint(210.117,167.092),EPoint(205.0,50.0)),
+	# 	Line(EPoint(127.7,25.6),EPoint(120.172,179.405)),
+	# 	Line(EPoint(127.809,58.481),EPoint(124.324,127.427)),
+	# 	Line(EPoint(85.964,191.64),EPoint(97.68,14.477)),
+	# 	Line(EPoint(56.447,124.662),EPoint(54.889,137.918))
+	# ]
+	# linePack=[
+	# 	Line(EPoint(56.447,124.662),EPoint(139.695,128.104)),
+	# 	Line(EPoint(288.267,74.433),EPoint(191.485,71.227)),
+	# 	Line(EPoint(252.926,29.71),EPoint(174.96,27.253)),
+	# 	Line(EPoint(274.412,120.07),EPoint(41.065,112.759)),
+	# 	Line(EPoint(289.674,108.019),EPoint(26.17,100.538)),
+	# 	Line(EPoint(240.702,107.818),EPoint(41.065,112.759)),
+	# 	Line(EPoint(174.457,166.127),EPoint(88.192,164.5))
+	# ]
+	# for (i,p) in enumerate(linePack):
+	# 	for q in linePack[i+1:]:
+	# 		print(p.intersect(q))
+	# show(linesImg)
+
+	# parallels=[]
+	# for (i,ab) in enumerate(lines):
+	# 	for cd in lines[i+1:]:
+	# 		if ab.computeAngle(cd)>math.pi/4:
+	# 			continue
+	# 		parallels.append((ab,cd))
+	# print("parallel lines candidate pairs:",len(parallels))
+
+	# parallels=list(groupParallels(lines,50,w))
+	# print("parallel triples:",len(parallels))
+	# for (p,q,r) in parallels:
+	# 	print(p,q,r)
+	# 	print(p.intersect(q))
+	# 	print(p.intersect(r))
+	# 	print(q.intersect(r))
+	# 	print()
+	# 	matsImg_=np.copy(matsImg)
+	# 	for pi in (p,q,r):
+	# 		cv.line(matsImg_,(int(pi.a.x),int(pi.a.y)),(int(pi.b.x),int(pi.b.y)),(0,255,255),1)
+	# 	show(matsImg_)
+
+	# p=Line(corners[0],corners[1])
+	# q=Line(corners[2],corners[3])
+	# r=Line(corners[1],corners[2])
+	# s=Line(corners[3],corners[0])
+	#
+	# matrix=computeRectiMatrix(p,q,r,s)
+	# print(matrix)
+	# for (point,contour) in stoneLocs:
+	# 	point_=EPoint.fromProjective(transformPoint(point.toProjective(),matrix))
+	# 	# print(point,"->",point_)
+	# 	cv.line(matsImg,(int(point.x),int(point.y)),(int(point_.x),int(point_.y)),(0,255,255),1)
+	# points1=np.float32([[p.x-x1,p.y-y1] for p in corners])
+	# points2=np.float32([[0,0],[0,400],[400,400],[400,0]])
+	# print(points1)
+	# print(points2)
+	# mat=cv.getPerspectiveTransform(points1,points2)
+	# print(mat)
+	# show(rect)
+	# warped=cv.warpPerspective(rect,matrix,(400,400))
+	# show(warped)
+
+	# mats=[]
+	# for (i,(p,q)) in enumerate(parallels):
+	# 	for (r,s) in parallels[i+1:]:
+	# 		if p is r or p is s or q is r or q is s:
+	# 			continue
+	# 		matrix=computeRectiMatrix(p,q,r,s)
+	# 		score=scoreMatrix(matrix,p,r,lines)
+	# 		mats.append((score,p,q,r,s,matrix))
+	# mats.sort(key=lambda x: x[0])
+	# for (score,p,q,r,s,matrix) in mats[:4]:
+	# 	print(score,p,q,r,s,matrix)
+	# 	matsImg_=np.copy(matsImg)
+	# 	for ab in (p,q,r,s):
+	# 		((xa,ya),(xb,yb))=(ab.a,ab.b)
+	# 		cv.line(matsImg_,(int(xa),int(ya)),(int(xb),int(yb)),(255,255,0),1)
+	# 	show(matsImg_)
+	# for (score,p,q,r,s,matrix) in mats[-4:]:
+	# 	print(score,p,q,r,s,matrix)
+	# 	matsImg_=np.copy(matsImg)
+	# 	for ab in (p,q,r,s):
+	# 		((xa,ya),(xb,yb))=(ab.a,ab.b)
+	# 		cv.line(matsImg_,(int(xa),int(ya)),(int(xb),int(yb)),(0,0,255),1)
+	# 	show(matsImg_)
diff --git a/src/analyzer/epoint.py b/src/analyzer/epoint.py
--- a/src/analyzer/epoint.py
+++ b/src/analyzer/epoint.py
@@ -1,6 +1,12 @@
 import math
 
 
+def homogenize(v):
+	for k in v:
+		if k!=0: return v/k
+	return v
+
+
 ## Euclidean 2D plane point: (x,y).
 class EPoint:
 	def __init__(self,x,y):
@@ -21,6 +27,12 @@ class EPoint:
 	def toProjective(self):
 		return (1,self._x,self._y)
 
+	def toPolar(self,center):
+		v=self-center
+		alpha=math.atan2(v.y,v.x)
+		k=self.dist(center)
+		return (alpha,k)
+
 	def dist(self,a):
 		return math.sqrt((self._x-a._x)**2+(self._y-a._y)**2)