for (i,c) in enumerate(contours):

		keep=True

		moments=cv.moments(c)

		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.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)

		hull=cv.convexHull(c)

		coverage2=area/(cv.contourArea(hull) or 1)

		if keep:

			res.append((EPoint(*center),c))

			cv.drawMarker(contourImg,tuple(map(int,center)),(255,0,0),cv.MARKER_CROSS)

	log.debug("accepted: %s",len(res))

	log.debug("rejected: %s",len(contours)-len(res))

	show(contourImg,"accepted and rejected stones")

	return res

def groupLines(points,minCount,tolerance):

	random.shuffle(points)

	sample=points[:]

		gray=cv.cvtColor(rect,cv.COLOR_BGR2GRAY)

		edges=cv.Canny(gray,70,130)

		show(edges,"edges")

		quantized=quantized & (255-cv.cvtColor(edges,cv.COLOR_GRAY2BGR))

		show(quantized,"quantized, edges separated")

		# detect black and white stones

		stones=self._detectStones(quantized,colors)

		# detect lines from edges and stones

		edgeImg=prepareEdgeImg(rect)

		hough=HoughTransform(edgeImg)

		stonesImg=np.zeros((self._rectH,self._rectW),np.uint8)

		for (point,c) in stones:

			cv.circle(stonesImg,(int(point.x),int(point.y)),2,255,-1)

		# cv.drawContours(stonesImg,[c for (point,c) in stones],-1,255,-1)

		show(stonesImg,"detected stones")

		hough.extract()

		# # detect lines passing through the stones

		# lines=self._constructLines(stones)

		#

		# # detect vanishing points of the lines

		# imgCenter=EPoint(w//2-x1, h//2-y1)

		# (a,b,c,d)=(p-EPoint(x1,y1) for p in self._annotations[filename][0])

		# (p,q,r,s)=(Line(a,b),Line(b,c),Line(c,d),Line(d,a))

		# v1=p.intersect(r)

		# v2=q.intersect(s)

		# log.debug("true vanishing points: %s ~ %s, %s ~ %s",v1,v1.toPolar(imgCenter),v2,v2.toPolar(imgCenter))

	def __init__(self,img):

		self._angleBandwidth=30 # degrees

		(h,w)=img.shape[:2]

		self._diagLen=int(np.sqrt(h**2+w**2))+1

		self._center=(w//2,h//2)

		self._acc=np.zeros((180,self._diagLen),dtype=np.int32)

		self.update(img)

	def extract(self):

		shift=self._diagLen//2

		peaks=[(alpha,d-shift) for (alpha,d) in peaks]

		peaks=self._filterClose(peaks)

		peaks.sort(key=lambda rc: rc[0])

		log.debug("detected peaks: %s",peaks)

		(alpha,beta)=self._detectDominantAngles(peaks)

		img=self._createImg()

		doublePeaks=peaks+[(alpha+180,-d) for (alpha,d) in peaks]

		params=self._computeGridParams([(gamma,d+shift) for (gamma,d) in doublePeaks if alpha<=gamma<=alpha+self._angleBandwidth])

		self._drawGridCurve(img,params,0)

		params=self._computeGridParams([(gamma,d+shift) for (gamma,d) in doublePeaks if beta<=gamma<=beta+self._angleBandwidth])

		self._drawGridCurve(img,params,1)

		self.show(img)

	def update(self,img,weight=1):

		start=datetime.now().timestamp()

		for (r,row) in enumerate(img):

			for (c,pix) in enumerate(row):

				if pix==0: continue

	def _computeDist(self,x,y,alphaDeg):

		alphaRad=alphaDeg*math.pi/180

		(x0,y0)=self._center

		(dx,dy)=(x-x0,y-y0)

		d=dx*math.cos(alphaRad)+dy*math.sin(alphaRad)

		return round(d)

	def _filterClose(self,peaks): # a naive implementation

		"""Discard points with Euclidean distance on the original image lower than 10.

		From such pairs keep only the one with a higher value in the accumulator.

		This can delete a series of points. If a-b and b-c are close and a>b>c, only a is kept."""

		center=EPoint(*self._center)

		res=[]

		for (alphaDeg,d) in peaks:

			alphaRad=alphaDeg*math.pi/180

			point=EPoint.fromPolar((alphaRad,d),center)

			ctrl=True

			for (betaDeg,e) in peaks:

				betaRad=betaDeg*math.pi/180

				point_=EPoint.fromPolar((betaRad,e),center)

				if point.dist(point_)<minDist and self._acc[(alphaDeg,d)]<self._acc[(betaDeg,e)]:

					ctrl=False

			if ctrl: res.append((alphaDeg,d))

		bandwidth=self._angleBandwidth

		k1=0

		k2=1

		histogram=[]

		while k1<n:

			while (k2<n and angles[k1]+bandwidth>angles[k2]) or (k2>=n and angles[k1]+bandwidth>angles[k2%n]+180):

				k2+=1

			histogram.append((angles[k1],k2-k1))

			k1+=1

		log.debug("angles histogram: %s",histogram)

		dominantAngles=sorted(histogram,key=lambda xy: xy[1],reverse=True)

		alpha=dominantAngles[0]

		beta=dominantAngles[0]

		log.debug("dominant angles: %s, %s",alpha,beta)

		return (alpha[0],beta[0])

	def _computeGridParams(self,lines):

		log.debug("computing grid parameters for: %s",lines)

		angles=[alpha for (alpha,d) in lines]

		dists=[d for (alpha,d) in lines]

		curve=lambda x,a,b,c,d: a*x**3+b*x**2+c*x+d

		(params,cov)=scipy.optimize.curve_fit(curve,dists,angles)

		log.debug("result: %s",params)

		return params

		arr=np.expand_dims(np.uint8(255*self._acc//maxVal),axis=2)

		img=np.concatenate((arr,arr,arr),axis=2)

		(h,w)=img.shape[:2]

		for x in range(0,w,4): # y axis

			img[h//2,x]=[255,255,255]

		for y in range(0,h,4):

			img[y,w//2]=[255,255,255]

		return img

	def _drawGridCurve(self,img,params,colorKey=0):

		colors=[[0,255,255],[255,0,255],[255,255,0]]

		color=colors[colorKey]

		(a,b,c,d)=params

		(h,w)=img.shape[:2]

		curve=lambda x: a*x**3+b*x**2+c*x+d

		for x in range(0,w,3):

			y=int(curve(x))

			if y<0 or y>=2*h: continue

			if y<h:	img[y,x]=color

			else: img[y%h,-x]=color