def blackout_outside(self, img, sigma=3):
        img_g = skic.rgb2gray(img)
        edges = skif.canny(img_g, sigma=sigma)

        hough_radii = np.arange(180, 210, 2)
        hough_res = skit.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 = 1
            peaks = skif.peak_local_max(h, min_distance=40, num_peaks=num_peaks)
            if peaks != []:
                centers.extend(peaks)
                accums.extend(h[peaks[:, 0], peaks[:, 1]])
                radii.extend([radius] * num_peaks)

#                print radius, np.max(h.ravel()), len(peaks)

        if accums == [] and sigma==3:
            return self.blackout_outside(img, sigma=3)

    #     Draw the most prominent 5 circles
        image = (img.copy() / 4.0).astype(np.uint8)
        cx, cy = skid.circle(*self.average_hough_detections(hough_radii, hough_res))
        image[cy, cx] = img[cy, cx]

        return image

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 find_inner_circle_parameters(plane_array, rmin=200, rmax=250):
    """Given a single planar image (e.g. a section from CT data), find the
    locations of the inner of two circles in the image."""

    xdim, ydim = plane_array.shape

    edges = find_edges_sobel(plane_array)

    hough_radii = np.arange(rmin, rmax, 3)
    hough_res = hough_circle(edges, hough_radii)

    # Find the two clearest circles
    c1, c2 = find_n_best_hough_circles(hough_radii, hough_res, 2)

    # Work out which is the inner circle
    r1 = c1[2]
    r2 = c2[2]
    if r1 > r2:
        inner_circle_radius = r2
        cx, cy, r = c2
    else:
        inner_circle_radius = r1
        cx, cy, r = c1

    return cx, cy, r

def test_hough_circle():
    # Prepare picture
    img = np.zeros((120, 100), dtype=int)
    radius = 20
    x_0, y_0 = (99, 50)
    y, x = circle_perimeter(y_0, x_0, radius)
    img[x, y] = 1

    out1 = tf.hough_circle(img, radius)
    out2 = tf.hough_circle(img, [radius])
    assert_equal(out1, out2)
    out = tf.hough_circle(img, np.array([radius], dtype=np.intp))
    assert_equal(out, out1)
    x, y = np.where(out[0] == out[0].max())
    assert_equal(x[0], x_0)
    assert_equal(y[0], y_0)

def detect_Hough(data):
    image = data.copy()
    edges = canny(image, sigma=10, low_threshold=60, high_threshold=90)

    # Detect circles between 80% and 100% of image semi-diagonal
    lx, ly = data.shape
    sizex, sizey = lx/2., ly/2.
    max_r = numpy.sqrt(sizex**2 + sizey**2) 
    hough_radii = numpy.linspace(0.5*max_r, 0.9 * max_r, 20)
    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)

    # Use the most prominent circle
    idx = numpy.argsort(accums)[::-1][:1]
    center_x, center_y = centers[idx]
    radius = radii[idx]
    return center_x, center_y, radius

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 process(self, im):
        (width, height, _) = im.image.shape

        img_adapted = im.prep(self.transform)

        if width > self.max_resized or height > self.max_resized:
            scale_height = self.max_resized / height
            scale_width = self.max_resized / width
            scale = min(scale_height, scale_width)
            img_adapted = resize(img_adapted, (int(width * scale), int(height * scale)))

        edges = canny(img_adapted, sigma=self.sigma)

        # Detect two radii
        # Calculate image diameter
        shape = im.image.shape
        diam = math.sqrt(shape[0] ** 2 + shape[1] ** 2)
        radii = np.arange(diam / 3, diam * 0.8, 2)
        hough_res = hough_circle(edges, radii)

        accums = []
        for radius, h in zip(radii, hough_res):
            # For each radius, extract two circles
            peaks = peak_local_max(h, num_peaks=1, min_distance=1)
            if len(peaks) > 0:
                accums.extend(h[peaks[:, 0], peaks[:, 1]])

        if len(accums) == 0:  # TODO: fix, should not happen
            return [0]

        idx = np.argmax(accums)
        return [accums[idx]]

def analyze_ra_rotation(rotate_fn):
    """ Get RA axis center of rotation XY coordinates

    Args:
        rotate_fn (str): FITS file of RA rotation

    Returns:
        tuple(int): A tuple of integers corresponding to the XY pixel position
        of the center of rotation around RA axis
    """
    d0 = fits.getdata(rotate_fn)  # - 2048

    # Get center
    position = (d0.shape[1] // 2, d0.shape[0] // 2)
    size = (1500, 1500)
    d1 = Cutout2D(d0, position, size)

    d1.data = d1.data / d1.data.max()

    # Get edges for rotation
    rotate_edges = canny(d1.data, sigma=1.0)

    rotate_hough_radii = np.arange(100, 500, 50)
    rotate_hough_res = hough_circle(rotate_edges, rotate_hough_radii)
    rotate_accums, rotate_cx, rotate_cy, rotate_radii = \
        hough_circle_peaks(rotate_hough_res, rotate_hough_radii, total_num_peaks=1)

    return d1.to_original_position((rotate_cx[-1], rotate_cy[-1]))

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 _detect_spots_hough_circle(image, radius):
	edges = canny(image)
	imshow(edges)
	show()
	hough_radii = np.arange(radius/2, radius*2, 10)
	hough_circles = hough_circle(edges, hough_radii)
	
	print(hough_circles)

def count_dishes(out, sink):
    edges = canny(
        sink,
        sigma=2,
        #low_threshold=10,
        high_threshold=0.3
    )
    
    hough_radii = np.arange(25, 70, 1)
    hough_res = hough_circle(edges, hough_radii)
    centers = []
    accums = []
    radii = []

    for radius, h in zip(hough_radii, hough_res):
        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)

    sink = color.gray2rgb(sink)
    hits = {}

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

        if is_sink_hole(center_x, center_y, radius):
            continue

        for d in hits.keys():
            dx, dy, dr = d

            dt = distance(center_x, dx, center_y, dy)
            
            if dt <= 40 and abs(dr - radius) < 50:
                hits[d] += 1
                break
        else:
            hits[(center_x, center_y, radius)] = 1


    dishes = [k for k,v in hits.iteritems()]

    for dish in dishes:
        center_x, center_y, radius = dish

        cx, cy = circle_perimeter(center_y, center_x, radius)

        try:
            sink[cy, cx] = (220, 250, 20)
        except IndexError:
            continue


    draw_res(out, sink, edges)
    return len(dishes)

def __tutorial_hough_circle_detection_skiimage(img_path, min_dim=40, max_dim=60):
    """
    This algorithm is crap, the one using opencv is much much better
    :param img_path:
    :param min_dim:
    :param max_dim:
    :return:
    """

    # load picture and detect edges
    img = cv2.imread(img_path)
    img = cv2.cvtColor(img, cv2.COLOR_BGRA2GRAY)
    img_edges = canny(img, sigma=3, low_threshold=10, high_threshold=50)

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

    # detect all radii within the range
    min_radius = int(min_dim / 2)
    max_radius = int(max_dim / 2)
    hough_radii = np.arange(start=min_radius, stop=max_radius, step=1)
    hough_res = hough_circle(img_edges, hough_radii)

    for radius, h in zip(hough_radii, hough_res):
        # for each radius, extract 2 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)

    img_width = img.shape[1]
    img_height = img.shape[0]

    # get the sorted accumulated values
    # accums_sorted = np.asarray(accums)
    # accums_sorted = accums_sorted[idx_sorted]

    # don't consider circles with accum value less than the threshold
    accum_threshold = 0.3
    idx_sorted = np.argsort(accums)[::-1]

    # draw the most prominent n circles, i.e those
    # with the highest n peaks
    img_color = color.gray2rgb(img)
    for idx in idx_sorted:
        if accums[idx] < accum_threshold:
            continue
        center_y, center_x = centers[idx]
        radius = radii[idx]
        cx, cy = circle_perimeter(center_x, center_y, radius)
        cx[cx > img_width - 1] = img_width - 1
        cy[cy > img_height - 1] = img_height - 1
        img_color[cy, cx] = (220, 20, 20)

    # save the result
    skimage.io.imsave("D://_Dataset//GTSDB//Test_Regions//_img2_.png", img_color)

def punchhole_removal(im):
    import numpy as np
    from PIL import Image
    from skimage import io
    from skimage.color import rgba2rgb, rgb2gray
    from skimage.transform import hough_circle, hough_circle_peaks
    from skimage.feature import canny
    from skimage.draw import circle
    from skimage.util import img_as_ubyte

    ''' check for punch holes and remove  '''
    
    max_peaks =  24 #maximum number of peaks to be found. changed from 99 to 24 for reducing the unnecessary punch holes being filled.

    img = np.array(im)# Load picture .
    img_rgb = rgba2rgb(img)# convert to RGB
    img_gray = rgb2gray(img_rgb)# convert to gray
    image = img_as_ubyte(img_gray)
    width, height = image.shape
    x1 =  punchhole_margin
    x2 =  (int)(width - punchhole_margin)
    y1 =  (int)(height - punchhole_margin)
    y2 =  punchhole_margin

    edges = canny(image, 3, 10, 40) # perform canny to detect the edges
    hough_radii = np.arange(31, 34, 1) #get the radius range with step as 1.
    hough_res = hough_circle(edges, hough_radii) # detect the circles centres coordinates

    # Select the most prominent circles based on the max_peaks
    accums, cx, cy, radii = hough_circle_peaks(hough_res, hough_radii,total_num_peaks=max_peaks)

    for center_y, center_x, radius in zip(cy, cx, radii):

        #if the circles centres fall in the border regions, 
        #get the dominant color near the hole and fill the hole with a linear gradient of the dominant color
        if(((0 < center_y < width) and (0 < center_x < y2)) or \
           ((0 < center_y < width) and (y1 < center_x < height)) or\
           ((0 < center_y < x1) and (0 < center_x < height)) or \
           ((x2 < center_y < width) and (0 < center_x < height))):

            index=0
            rr, cc= circle(center_y, center_x, radius+1, img.shape)
            dominantpix = dominantcolor(center_x, center_y, radius, img)           
            dark_grad = [dominantpix[0], dominantpix[1],dominantpix[2]]
            light_grad = [dominantpix[0]+1, dominantpix[1]+1, dominantpix[2]+1]
            #white_grad = [255,255,255]
            RGBA_list = lineargradient(dark_grad,light_grad,len(list(rr)))   
          
            for i , j in zip(list(rr), list(cc)):
                pixlist = RGBA_list[index]
                pixtuple = tuple(pixlist)
                img[i,j]= (pixtuple[0], pixtuple[1], pixtuple[2], 255)
                index += 1
           
    finalimage=Image.fromarray(img)

    return finalimage

def find_boundary(img, min_radius, max_radius):
    edges = canny(img, sigma=0.)
    radii = np.arange(min_radius, max_radius)
    h = hough_circle(edges, radii)
    h_max = np.max(h, axis=(1, 2))
    k = np.min(np.argsort(h_max)[-2:])
    radius = radii[k]
    h = h[k, ...]
    row, col = np.unravel_index(np.argmax(h), h.shape)
    return row, col, radius

    def find_circle(self, image, frame, dim, **kw):
        dx, dy = None, None

        pframe = self._preprocess(frame, blur=0)
        edges = canny(pframe, sigma=3)
        hough_radii = arange(dim * 0.9, dim * 1.1, 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)

        # for idx in argsort(accums)[::-1][:1]:
        try:
            idx = argsort(accums)[::-1][0]
        except IndexError:
            return dx, dy

        center_y, center_x = centers[idx]
        radius = radii[idx]

        draw_circle_perimeter(frame, center_x, center_y, radius, (220, 20, 20))
        # cx, cy = circle_perimeter(int(center_x), int(center_y), int(radius))

        # draw perimeter
        # try:
        #     frame[cy, cx] = (220, 20, 20)
        # except IndexError:
        #     pass

        # draw center
        # cx, cy = circle(int(center_x), int(center_y), int(2))
        # frame[cy, cx] = (220, 20, 20)
        draw_circle(frame, center_x, center_y, 2, (220, 20, 20))

        h, w = frame.shape[:2]

        ox, oy = w / 2, h / 2
        dx = center_x - ox
        dy = center_y - oy

        cx, cy = circle(int(ox), int(oy), int(2))
        frame[cy, cx] = (20, 220, 20)

        image.set_frame(frame)
        return float(dx), -float(dy)

def find_circles(image, min_r=30, max_r=300, cutoff=0.5, step=4, blur_sigma=3):
    rs = arange(min_r, max_r, step)
    threshold_image = threshold(image)
    hough_space = hough_circle(threshold(image), rs)
    sigma_3d = (blur_sigma,
                blur_sigma,
                float(blur_sigma) / step)
    blurred_hough_space = gaussian_filter(hough_space, sigma=sigma_3d)
    local_maxima = peak_local_max(blurred_hough_space, exclude_border=False)
    circles = column_stack((local_maxima, hough_space[tuple(local_maxima.T)]))
    circles[:,0] = rs[list(circles[:,0])]
    return circles[circles[:,3] > cutoff]

def test_hough_circle():
    # Prepare picture
    img = np.zeros((120, 100), dtype=int)
    radius = 20
    x_0, y_0 = (99, 50)
    x, y = circle_perimeter(y_0, x_0, radius)
    img[y, x] = 1

    out = tf.hough_circle(img, np.array([radius]))

    x, y = np.where(out[0] == out[0].max())
    # Offset for x_0, y_0
    assert_equal(x[0], x_0 + radius)
    assert_equal(y[0], y_0 + radius)

def extract_hough_circle(img_rgb, img_gray, out_filepath):
    # Canny
    img = img_as_ubyte(img_gray)
    edges = canny(img, sigma=3, low_threshold=10, high_threshold=50)

    # fig, ax = plt.subplots(nrows=1, ncols=1)
    # ax.imshow(edges, cmap=plt.cm.gray)
    # ax.axis('off')    
    # ax.set_title('Canny Edges for Hough Circle', fontsize=18)
    # plt.tight_layout()
    # plt.savefig('canny_edges_for_hough_circle.png')

    # Detect
    min_radii = 15; max_radii = 30; step_radii = 1
    plausible_radii = np.arange(min_radii, max_radii, step_radii)

    hough_circles = hough_circle(edges, plausible_radii)

    centers = []; accums = []; radii = []
    for radius, h in zip(plausible_radii, hough_circles):
        n_extracted_circle = 1 # ...for each radius
        peaks = peak_local_max(h, num_peaks=n_extracted_circle)
        centers.extend(peaks)
        accums.extend(h[peaks[:, 0], peaks[:, 1]])
        radii.extend([radius] * n_extracted_circle)

    # Draw the most prominent circles
    n_top_circle = 15
    fig, ax = plt.subplots(ncols=1, nrows=1)
    for idx in np.argsort(accums)[::-1][:n_top_circle]:
        center_x, center_y = centers[idx]
        center_color = (0, 225, 0)
        img_rgb[center_x, center_y] = center_color

        radius = radii[idx]
        perim_color = (255, 0, 0)
        perim_x_list, perim_y_list = circle_perimeter(center_y, center_x, radius)
        # if all(i < img_rgb.shape[1] for i in perim_x_list) and all(i < img_rgb.shape[0] for i in perim_y_list):
        #     img_rgb[perim_y_list, perim_x_list] = perim_color
        for perim_x, perim_y in zip(perim_x_list, perim_y_list):
            if perim_x < img_rgb.shape[1] and perim_y < img_rgb.shape[0]:
                img_rgb[perim_y, perim_x] = perim_color

    ax.imshow(img_rgb, cmap=plt.cm.gray)
    ax.axis('off')
    ax.set_title('Hough Circle', fontsize=18)
    plt.tight_layout()
    plt.savefig(out_filepath)
    plt.close(fig)