diff --git a/exp/histogram.py b/exp/histogram.py
--- a/exp/histogram.py
+++ b/exp/histogram.py
@@ -4,8 +4,13 @@ import sys
 import cv2 as cv
 import numpy as np
 import scipy.cluster
+import scipy.ndimage
+from matplotlib import pyplot as plt
+import PIL.Image
+from PIL.ImageDraw import ImageDraw
 
 from annotations import DataFile,computeBoundingBox
+from hough import show,houghLines
 
 
 def createHistogram(img):
@@ -64,6 +69,27 @@ def computeClosest(x,centers):
 	return res
 
 
+def score(arr1,arr2):
+	try:
+		return (arr1&arr2).sum() / ((arr1|arr2).sum() or 1)
+	except TypeError:
+		print(type(arr1),type(arr2))
+		print(arr1.shape,arr2.shape)
+		print(arr1.dtype,arr2.dtype)
+		raise TypeError()
+
+def maxOp55(arr):
+	m=arr.max()
+	return 1 if m>127 and arr[2,2]==m else 0
+
+def ellipse(a,b):
+	img=PIL.Image.new("1",(a,b))
+	d=ImageDraw(img)
+	d.ellipse((1,1,a-1,b-1),fill=1)
+	img.save("/tmp/ellipse.png")
+	return np.array(img,dtype=np.uint8)
+
+
 filepath=sys.argv[1]
 annotations=DataFile(sys.argv[2])
 filename=os.path.basename(filepath)
@@ -75,25 +101,54 @@ print(x3,x4,y3,y4)
 rect=img[y3:y4,x3:x4,:]
 centers=quantize(rect)
 
-# colors=[(0,0,0),(255,255,255),(255,200,0)]
-# colorMap=colors[:]
-# for (i,c) in enumerate(centers):
-# 	colorMap[i]=computeClosest(c,colors)
-
 for x in range(x1,x2):
 	for y in range(y1,y2):
 		pix=img[y,x]
 		img[y,x]=computeClosest(pix,centers)
 
 rect=img[y1:y2,x1:x2]
-maskB=cv.inRange(rect,np.array([0,0,0]),np.array([60,60,60]))
+maskB=cv.inRange(rect,np.array([0,0,0]),np.array([80,80,80]))
 maskB=cv.erode(maskB,np.ones((3,3),np.uint8))
+# maskB=cv.erode(maskB,np.ones((3,3),np.uint8))
+# maskB=cv.erode(maskB,np.ones((3,3),np.uint8))
 maskW=cv.inRange(rect,np.array([160,160,160]),np.array([256,256,256]))
 maskW=cv.erode(maskW,np.ones((3,3),np.uint8))
+# maskW=cv.erode(maskW,np.ones((3,3),np.uint8))
+# maskW=cv.erode(maskW,np.ones((3,3),np.uint8))
 
-cv.imshow(filename,img)
-cv.waitKey()
-cv.imshow(filename,maskB)
-cv.waitKey()
-cv.imshow(filename,maskW)
-cv.waitKey()
+show(img,filename)
+show(maskB,filename)
+show(maskW,filename)
+stones=cv.bitwise_or(maskB,maskW)
+# houghLines(stones)
+
+(bh,bw)=stones.shape
+sw=bw//19
+sh=bh//19
+print(stones.shape,(sw,sh))
+ell=ellipse(sw,sh)*255
+# print(ell)
+hitMap=np.zeros_like(stones,dtype=np.uint8)
+for i in range(sw,bw):
+	for j in range(sh,bh):
+		region=stones[j-sh:j, i-sw:i]
+		hitMap[j,i]=255*score(region,ell)
+show(hitMap)
+
+gridMap=np.zeros_like(hitMap,dtype=np.uint8)
+for i in range(5,bw):
+	for j in range(5,bh):
+		region=hitMap[j-5:j, i-5:i]
+		gridMap[j,i]=255*maxOp55(region)
+show(gridMap)
+
+houghLines(gridMap)
+# ministones=cv.resize(stones,None,fx=0.25,fy=0.25,interpolation=cv.INTER_AREA)
+# dft = cv.dft(np.float32(ministones),flags = cv.DFT_COMPLEX_OUTPUT)
+# dft_shift = np.fft.fftshift(dft)
+# magnitude_spectrum = 20*np.log(cv.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))
+# plt.subplot(121),plt.imshow(stones, cmap = 'gray')
+# plt.title('Input Image'), plt.xticks([]), plt.yticks([])
+# plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
+# plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
+# plt.show()