def test_circle_perimeter_andres():
    img = np.zeros((15, 15), "uint8")
    rr, cc = circle_perimeter(7, 7, 0, method="andres")
    img[rr, cc] = 1
    assert np.sum(img) == 1
    assert img[7][7] == 1

    img = np.zeros((17, 15), "uint8")
    rr, cc = circle_perimeter(7, 7, 7, method="andres")
    img[rr, cc] = 1
    img_ = np.array(
        [
            [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
            [0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
            [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0],
            [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
            [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
            [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
            [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
            [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
            [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
            [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
            [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
            [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
            [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0],
            [0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
            [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        ]
    )
    assert_array_equal(img, img_)

def test_circle_perimeter_andres():
    img = np.zeros((15, 15), 'uint8')
    rr, cc = circle_perimeter(7, 7, 0, method='andres')
    img[rr, cc] = 1
    assert(np.sum(img) == 1)
    assert(img[7][7] == 1)

    img = np.zeros((17, 15), 'uint8')
    rr, cc = circle_perimeter(7, 7, 7, method='andres')
    img[rr, cc] = 1
    img_ = np.array(
        [[0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
         [0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
         [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0],
         [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
         [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
         [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
         [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
         [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
         [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
         [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
         [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
         [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
         [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0],
         [0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
         [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
    )
    assert_array_equal(img, img_)

def test_circle_perimeter_bresenham_shape():
    img = np.zeros((15, 20), 'uint8')
    rr, cc = circle_perimeter(7, 10, 9, method='bresenham', shape=(15, 20))
    img[rr, cc] = 1
    shift = 5
    img_ = np.zeros((15 + 2 * shift, 20), 'uint8')
    rr, cc = circle_perimeter(7 + shift, 10, 9, method='bresenham', shape=None)
    img_[rr, cc] = 1
    assert_array_equal(img, img_[shift:-shift, :])

def show_matching(img, img2, matching):
    print "matching..."
    ip_match = build_match_dic(img, img2, matching)

    print "Constructing intermediate image..."
    padding = 5 #padding around the edges

    bar = np.ndarray((img.shape[0], 5))
    bar.fill(1.0)
    viewer = ImageViewer(img)
    viewer.show()
    img3 = np.column_stack((img, bar, img2))
    viewer = ImageViewer(img3)
    viewer.show()

    img3 = img_as_ubyte(img3)

    viewer = ImageViewer(img3)
    viewer.show()

    img3 = np.pad(img3, pad_width=padding, mode='constant', constant_values=(0))


    viewer = ImageViewer(img3)
    viewer.show()
    print "Drawing lines..."

    colimg = color.gray2rgb(img3)
    for k,v in random.sample(ip_match.items(), int(float(len(ip_match.keys()))*0.005)):
        #Choose a random colour:
        col = [random.randint(0,255),random.randint(0,255),random.randint(0,255)]

        #Calculate coordinates after padding:
        x1 = k[0]+padding
        y1 = k[1]+padding
        x2 = v[0]+padding
        y2 = v[1] + img.shape[1]+bar.shape[1]+padding

        #Draw the points in both images:
        rr, cc = circle_perimeter(x1, y1, 3)
        colimg[rr, cc] = col
        rr, cc = circle_perimeter(x2, y2, 3)
        colimg[rr, cc] = col

        #Draw a line between the points:
        rr, cc = line(x1,y1,x2,y2)
        colimg[rr, cc] = col

    #Show the result:
    viewer = ImageViewer(colimg)
    viewer.show()

def find_iris(image, pupil, **kwargs):
    buffer = 20
    # run canny
    image = filter.canny(image, sigma=1, low_threshold=10, high_threshold=50)
    cx, cy, radius = pupil

    segments = get_segments(400, step=0.01)
    # get ray directions
    directions = zip(map(cos, segments[0]), map(sin, segments[0]))
    shape = image.shape
    points = []
    for d in directions:
        start = (cx + (radius + buffer) * d[0], cy + (radius + buffer)*d[1])
        ray = Ray(image, start, d)
        point = ray.fire()
        if point != None:
            points.append(point)

    for p in points:
        x, y = circle_perimeter(int(p[0]), int(p[1]), 3)
        x = x[x < rgb.shape[0]]
        y = y[y < rgb.shape[1]]
        rgb[x,y] = (220, 40, 40)

    e = Ellipse().fit_with_center(None, points)
    return image, points, e

def test_circles2():
    data = np.memmap("E:\\guts_tracking\\data\\fish202_aligned_masked_8bit_150x200x440.raw", dtype='uint8', shape=(440,200,150)).copy()

    i = 157

    hough_radii = np.arange(10, 100, 10)
    edges = feature.canny(data[i], sigma=3.0, low_threshold=0.4, high_threshold=0.8)
    hough_res = hough_circle(edges, hough_radii)

    centers = []
    accums = []
    radii = []

    for radius, h in zip(hough_radii, hough_res):
        peaks = feature.peak_local_max(h)
        centers.extend(peaks)
        accums.extend(h[peaks[:, 0], peaks[:, 1]])
        radii.extend([radius] * len(peaks))

    image = ski.color.gray2rgb(data[i])

    for idx in np.argsort(accums)[::-1][:5]:
        center_x, center_y = centers[idx]
        radius = radii[idx]
        cx, cy = circle_perimeter(center_y, center_x, radius)

        if max(cx) < 150 and max(cy) < 200:
            image[cy, cx] = (220, 20, 20)

    plt.imshow(image, cmap='gray')

    plt.show()

    def animate(i):
        print 'Frame %d' % i
        plt.title('Frame %d' % i)

        image = data[i]

        hough_radii = np.arange(10, 100, 10)
        edges = feature.canny(data[i], sigma=3.0, low_threshold=0.4, high_threshold=0.8)
        hough_res = hough_circle(edges, hough_radii)

        centers = []
        accums = []
        radii = []

        for radius, h in zip(hough_radii, hough_res):
            peaks = feature.peak_local_max(h)
            centers.extend(peaks)
            accums.extend(h[peaks[:, 0], peaks[:, 1]])
            radii.extend([radius] * len(peaks))

        image = ski.color.gray2rgb(data[i])

        for idx in np.argsort(accums)[::-1][:5]:
            center_x, center_y = centers[idx]
            radius = radii[idx]
            cx, cy = circle_perimeter(center_y, center_x, radius)

            if max(cx) < 150 and max(cy) < 200:
                image[cy, cx] = (220, 20, 20)

        im.set_data(image)

        return im,

def gen_ring_mask(votes, img, y0, x0, r0):
    import skimage.draw as draw
    #rad = np.sum(np.sum(votes[:,y0-1:y0+2,x0-1:x0+2], axis=1), axis=1)
    rad = votes[:, y0, x0]
    peaks = np.where(np.r_[False, rad[1:] > rad[:-1]] 
                & np.r_[rad[:-1] > rad[1:], False]
                & (rad >= 0.2*rad.max()))
    peaks = peaks[0]
    coeffs = []
    xx = range(len(rad))
    flag = np.zeros(len(rad), dtype='uint8')

    # fit guassian peak
    for p in peaks:
        if p > 1.3*r0:
            break
        yy = np.copy(rad)
        yy[:p-3] = 0
        yy[p+4:] = 0
        cc = fit_gauss(xx, yy, rad[p], p, 2)
        flag[np.floor(cc[1]-cc[2]):np.ceil(cc[1]+cc[2])+1] = 1
        coeffs.append(cc)
    #print flag, coeffs

    # create image mask
    mask = np.zeros(img.shape, dtype='uint8')
    for _r, _f in enumerate(flag):
        if _f == 1:
             _y, _x = draw.circle_perimeter(y0, x0, _r, 'andres')
             mask[_y, _x] = 1
    return mask

        def onclick(self,event):
            if event.inaxes==ax:
                if event.button==1:
                    if np.sum(MASKs[self.ind])==0:
                        
                        cx, cy = circle_perimeter(int(event.ydata), int(event.xdata), 3)
                        MASKs[self.ind][cx,cy] = (220, 80, 20, 1)
                        original_image = np.copy(MASKs[self.ind][:,:,3])
                        chull = convex_hull_image(original_image)
                        MASKs[self.ind][chull,:] = (255, 0, 0, .4)


                        ax.cla()
                        data =  areaFile.attrs['ROI_patches'][:,:,self.ind]
                        ax.imshow(data,cmap='binary_r')
                        ax.imshow(MASKs[self.ind])

                elif event.button==3:
                    MASKs[self.ind] = np.zeros(MASKs[self.ind].shape)

                    
                    data = areaFile.attrs['ROI_patches'][:,:,self.ind]

                    ax.cla()
                    ax.imshow(data,cmap='binary_r')
                    ax.imshow(MASKs[self.ind])

def draw_keypoints(img, keypoints, draw_prob):
    """Draws for each keypoint a circle (roughly matching the sigma of the scale)
    with a line for the orientation.
    Args:
        img    The image to which to add the keypoints (gets copied)
        keypoints The keypoints to draw
        draw_prob Probability of drawing a keypoint (the lower the less keypoints are drawn)
    Returns:
        Image with keypoints"""
    height, width = img.shape
    img = np.copy(img)

    # convert from grayscale image to RGB so that keypoints can be drawn in red
    img = img[:, :, np.newaxis]
    img = np.repeat(img, 3, axis=2)

    for (y, x), orientation, scale_idx, scale_size, kp_type in keypoints:
        if draw_prob < 1.0 and random.random() <= draw_prob:
            # draw the circle
            radius = int(scale_size)
            rr, cc = draw.circle_perimeter(y, x, radius, shape=img.shape)
            img[rr, cc, 0] = 1.0
            img[rr, cc, 1:] = 0

            # draw orientation
            orientation = util.quantize(orientation, [-135, -90, -45, 0, 45, 90, 135, 180])
            x_start = x
            y_start = y
            if orientation == 0:
                x_end = x + radius
                y_end = y
            elif orientation == 45:
                x_end = x + radius*(1/math.sqrt(2))
                y_end = y + radius*(1/math.sqrt(2))
            elif orientation == 90:
                x_end = x
                y_end = y + radius
            elif orientation == 135:
                x_end = x - radius*(1/math.sqrt(2))
                y_end = y - radius*(1/math.sqrt(2))
            elif orientation == 180:
                x_end = x - radius
                y_end = y
            elif orientation == -135:
                x_end = x - radius*(1/math.sqrt(2))
                y_end = y - radius*(1/math.sqrt(2))
            elif orientation == -90:
                x_end = x
                y_end = y - radius
            elif orientation == -45:
                x_end = x + radius*(1/math.sqrt(2))
                y_end = y - radius*(1/math.sqrt(2))
            x_end = np.clip(x_end, 0, width-1)
            y_end = np.clip(y_end, 0, height-1)
            rr, cc = draw.line(int(y_start), int(x_start), int(y_end), int(x_end))
            img[rr, cc, 0] = 1.0
            img[rr, cc, 1:] = 0
    img = np.clip(img, 0, 1.0)

    return img

def quantify_yeast_growth(input_filename, annotation_filename, 
                            profile_filename):
    

    downscaled = load_and_downscale(input_filename)
    annotation = AnnotatedImage.from_grayscale(downscaled)

    circle = fit_central_circle(downscaled)

    circle_coords = circle_perimeter(*circle)
    annotation[circle_coords] = 0, 255, 0

    x, y, r = circle
    center = (x, y)

    xdim, ydim, = downscaled.shape
    line_length = ydim - y

    mean_profile_line = find_mean_profile_line(downscaled, annotation,
                                center, -math.pi/4, math.pi/4, line_length)


    record_line_profile(profile_filename, mean_profile_line)

    with open(annotation_filename, 'wb') as f:
        f.write(annotation.png())

	def remove_deep_points(self, CIRCLE_RADIUS = 5, MINIMUM_BOUNDARY = 70):
		for point in self.POINTS_initial:
			segments = []
			circle_x, circle_y = circle_perimeter(point[0],point[1],CIRCLE_RADIUS)
			circle_points = np.array(list(sorted(set(zip(circle_x,circle_y)))))

			sortedpoints = np.empty([len(circle_points), 2], dtype=int)

			test1 = circle_points[circle_points[:,0] == point[0] - CIRCLE_RADIUS]
			start = len(test1)
			end = len_cpoints = len(sortedpoints)
			sortedpoints[0:start] = test1

			for x in xrange(point[0] - CIRCLE_RADIUS + 1, point[0] + CIRCLE_RADIUS):
				test1 = circle_points[circle_points[:,0] == x]
				testlen = len(test1)
				if x <= point[0]:
					sortedpoints[start:start+testlen/2] = test1[testlen/2:]
					sortedpoints[end-testlen/2:end] = test1[:testlen/2]
				else:
					sortedpoints[start:start+testlen/2] = test1[testlen/2:][::-1]
					sortedpoints[end-testlen/2:end] = test1[:testlen/2][::-1]
				start += testlen/2
				end -= testlen/2

			test1 = circle_points[circle_points[:,0] == point[0] + CIRCLE_RADIUS]
			sortedpoints[start:start + len(test1)] = test1[::-1]

			for c_perimeter in sortedpoints:
				segments.append(True)
				line_x, line_y = line(point[0], point[1], c_perimeter[0], c_perimeter[1])
				for line_points in zip(line_x,line_y)[1:]:
					# if original_image[line_points[0]][line_points[1]] != 0:
					if self.FOOTPRINT_cleaned_opening[line_points[0]][line_points[1]] != 0:
						segments[-1] = False
						break

			min_boundpoints = (MINIMUM_BOUNDARY / 360.0) * len_cpoints

			seg_sizes = []
			new_segment = True
			for segment in segments:
				if segment:
					if new_segment:
						seg_sizes.append(1)
						new_segment = False
					else:
						seg_sizes[-1] += 1
				else:
					new_segment = True

			if segments[0] == True and segments[-1] == True and len(seg_sizes) > 1:
				seg_sizes[0] = seg_sizes[0] + seg_sizes[-1]
				seg_sizes.pop()
			
			if(len(seg_sizes) == 0 or max(seg_sizes) < min_boundpoints):
				# boundary_image[point[0]][point[1]] = 0							#TAMA BANG TANGGALIN?

				if (point[0],point[1]) in self.POINTS_ordered:  ## IDENTIFY KUNG BAKIT NAGKA-ERRROR, EXAMPLE pt000120_merged4.py obj 1301
					self.POINTS_ordered.remove((point[0],point[1]))

def form_mask_BR_v2(start_temp, max_temp, temp_rate):
	dtype = [('temp', float), ('mask', np.ndarray), ('dist', np.ndarray)]
	masks = []
	#IMPROVEMENT!!!
	#MINUS LAST COORDINATES TO NEXT TEMPERATURE
	#LESS POINTS TO CHECK
	while (start_temp <= max_temp):
		coordinates = []
		coor_dist = []
		distance = int(round(start_temp))
		array_size = distance * 2 + 1
		img = np.zeros((array_size, array_size), dtype=np.uint8)
		rr, cc = circle_perimeter(distance, distance, distance)
		img[rr, cc] = 1
		rr,cc = np.nonzero(ndimage.binary_fill_holes(img).astype(int))
		img[rr, cc] = 1

		for idx in xrange(0, len(rr)):
			dist_temp = np.linalg.norm(np.array([rr[idx], cc[idx]]) - np.array([distance, distance]))
			if dist_temp <= start_temp:
				coordinates.append([rr[idx], cc[idx]])
				coor_dist.append(dist_temp)

		# coordinates.remove([distance, distance])
		coordinates = coordinates - np.array([distance, distance])
		masks.append((start_temp, coordinates, np.array(coor_dist)))
		start_temp += temp_rate

	return np.array(masks, dtype=dtype)

def circle_markers(blobs, pic_shape):
    '''Return array with circles around foci found'''

    markers_rad = np.zeros(pic_shape, dtype = np.bool)

    x_max, y_max = pic_shape

    for blob in blobs:

        x, y, r = blob

        r = r*2

        rr, cc = circle_perimeter(x, y, np.round(r).astype(int))
        rr_new, cc_new = [], []

        for x_c,y_c in zip(rr,cc):

            if (x_c >= 0) and (x_c < x_max) and (y_c >= 0) and (y_c < y_max):
                rr_new.append(x_c)
                cc_new.append(y_c)

        markers_rad[rr_new, cc_new] = True

    selem = np.array([0,1,0,1,1,1,0,1,0], dtype=bool).reshape((3,3))

    for i in range(4):
        markers_rad = binary_dilation(markers_rad, selem)

    return markers_rad

def hugh_circle_detection(image):
	# Load picture and detect edges
	edges = canny(image, sigma=3, low_threshold=10, high_threshold=50)

	fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(5, 2))

	# Detect two radii
	hough_radii = np.arange(15, 30, 2)
	hough_res = hough_circle(edges, hough_radii)

	centers = []
	accums = []
	radii = []

	for radius, h in zip(hough_radii, hough_res):
		# For each radius, extract two circles
		num_peaks = 2
		peaks = peak_local_max(h, num_peaks=num_peaks)
		centers.extend(peaks)
		accums.extend(h[peaks[:, 0], peaks[:, 1]])
		radii.extend([radius] * num_peaks)

	# Draw the most prominent 5 circles
	image = color.gray2rgb(image)
	for idx in np.argsort(accums)[::-1][:5]:
		center_x, center_y = centers[idx]
		radius = radii[idx]
		cx, cy = circle_perimeter(center_y, center_x, radius)
		image[cy, cx] = (220, 20, 20)

	ax.imshow(image, cmap=plt.cm.gray)
	plt.show()

def compute_singularity_measure(gradients):
    """ Compute singularity measure of data
    """
    # precompute closed curve coordinates
    circle_rad = 5

    grad_dim = gradients[0].shape
    closed_curves = np.empty(grad_dim, dtype=(list, list))
    for j in range(grad_dim[1]):
        for i in range(grad_dim[0]):
            rr, cc = circle_perimeter(
                i, j, circle_rad,
                method='andres', shape=grad_dim
            )
            closed_curves[i, j] = (rr, cc)

    # find singularities
    singularities = []
    for grad in gradients:
        width, height = grad.shape

        singularity = np.empty((width, height))
        for j in range(height):
            for i in range(width):
                res = np.sum(grad[closed_curves[i, j]])
                singularity[i, j] = res
        singularity = np.array(singularity)

        singularities.append(singularity)

    return np.array(singularities)

def get_edges(np_image):
    # Get the coordinates of the circle edges
    X = []
    Y = []
    Z = []
    
    circles = []  # to get the outlines of the circles
    C = []  # to get the centres of the circles, in relation to the different areas
    R = []  # to get radii
    
    coords = np.column_stack(np.nonzero(np_image))
    X = np.array(coords[:, 0])
    Y = np.array(coords[:, 1])

    # Fit a circle and compare measured circle area with
    # area from the amount of pixels to remove trash
    XC, YC, RAD, RESID = leastsq_circle(X, Y)

    # Hough radius estimate
    hough_radii = np.arange(RAD - RAD / 2., RAD + RAD / 2.)
    
    # Apply Hough transform
    hough_res = hough_circle(np_image, hough_radii)
         
    centers = []
    accums = []
    radii = []
    
    img = np.zeros_like(np_image)
    
    # For each radius, extract one circle
    for radius, h in zip(hough_radii, hough_res):
        peaks = peak_local_max(h, num_peaks=2)
        centers.extend(peaks)
        accums.extend(h[peaks[:, 0], peaks[:, 1]])
        radii.extend([radius, radius])
         
    for idx in np.argsort(accums)[::-1][:1]:
        center_x, center_y = centers[idx]
        C.append((center_x, center_y))
        radius = radii[idx]
        R.append(radius)
        cx, cy = circle_perimeter(int(round(center_x, 0)),
                                  int(round(center_y, 0)),
                                  int(round(radius, 0)))
        circles.append((cx, cy))
        
    if C:
        for cent in C:
            xc, yc = cent
            np_image[int(xc), int(yc)] = 255
            np_image[int(xc)-5:int(xc)+5, int(yc)-5:int(yc)+5] = 255
    
    if circles:
        for per in circles:
            e1, e2 = per
            np_image[e1, e2] = 255
    
    return [C, R, circles, []], np_image

def add_auto_masks_area(areaFile,addMasks=False):
    from skimage.morphology import binary_dilation, binary_erosion, disk
    from skimage import exposure
    from skimage.transform import hough_circle
    from skimage.morphology import convex_hull_image
    from skimage.feature import canny
    from skimage.draw import circle_perimeter

    import h5py


    MASKs = np.zeros(areaFile.attrs['ROI_patches'].shape)
    if 'ROI_masks' in (areaFile.attrs.iterkeys()):
        print 'Masks have already been created'
        awns = raw_input('Would you like to redo them from scratch: (answer y/n): ')
    else:
        awns = 'y'

    if awns=='y':
        MASKs = np.zeros(areaFile.attrs['ROI_patches'].shape)
        for i in range(areaFile.attrs['ROI_patches'].shape[2]):
            patch = areaFile.attrs['ROI_patches'][:,:,i]
            
            tt0 = exposure.equalize_hist(patch)
            tt = 255*tt0/np.max(tt0)
            thresh = 1*tt>0.3*255
            thresh2 = 1*tt<0.1*255

            tt[thresh] = 255
            tt[thresh2] = 0



            edges = canny(tt, sigma=2, low_threshold=20, high_threshold=30)
            try_radii = np.arange(3,5)
            res = hough_circle(edges, try_radii)
            ridx, r, c = np.unravel_index(np.argmax(res), res.shape)
            
            r, c, try_radii[ridx]
            image = np.zeros([20,20,4])
            cx, cy = circle_perimeter(c, r, try_radii[ridx]+2)
            
            try:
                image[cy, cx] = (220, 80, 20, 1)
                original_image = np.copy(image[:,:,3])
                chull = convex_hull_image(original_image)
                image[chull,:] = (255, 0, 0, .4)
                
            except:
                pass
            
            
            MASKs[:,:,i] = (1*image[:,:,-1]>0)
        if addMasks:
            areaFile.attrs['ROI_masks'] = MASKs
    else:
        pass    

    return np.array(MASKs)

def get_perimeter(img):
    """
    :param img:
    :return:
    """
    cx = img.shape[0]/2
    cc, rr = circle_perimeter(cx - 1, cx - 1, 256)
    return rr, cc

def auto_find_center_rings(avg_img, sigma=1, no_rings=4, min_samples=3,
                           residual_threshold=1, max_trials=1000):
    """This will find the center of the speckle pattern and the radii of the
    most intense rings.

    Parameters
    ----------
    avg_img : 2D array
        shape of the image
    sigma : float, optional
        Standard deviation of the Gaussian filter.
    no_rings : int, optional
        number of rings
    min_sample : int, optional
        The minimum number of data points to fit a model to.
    residual_threshold : float, optional
        Maximum distance for a data point to be classified as an inlier.
    max_trials : int, optional
        Maximum number of iterations for random sample selection.

    Returns
    -------
    center : tuple
        center co-ordinates of the speckle pattern
    image : 2D array
        Indices of pixels that belong to the rings,
        directly index into an array
    radii : list
        values of the radii of the rings

    Note
    ----
    scikit-image ransac
    method(http://www.imagexd.org/tutorial/lessons/1_ransac.html) is used to
    automatically find the center and the most intense rings.
    """

    image = img_as_float(color.rgb2gray(avg_img))
    edges = feature.canny(image, sigma)
    coords = np.column_stack(np.nonzero(edges))
    edge_pts_xy = coords[:, ::-1]
    radii = []

    for i in range(no_rings):
        model_robust, inliers = ransac(edge_pts_xy, CircleModel, min_samples,
                                       residual_threshold,
                                       max_trials=max_trials)
        if i == 0:
            center = int(model_robust.params[0]), int(model_robust.params[1])
        radii.append(model_robust.params[2])

        rr, cc = draw.circle_perimeter(center[1], center[0],
                                       int(model_robust.params[2]),
                                       shape=image.shape)
        image[rr, cc] = i + 1
        edge_pts_xy = edge_pts_xy[~inliers]

    return center, image, radii