def test_blob_dog():
    r2 = math.sqrt(2)
    img = np.ones((512, 512))
    img3 = np.ones((5, 5, 5))

    xs, ys = circle(400, 130, 5)
    img[xs, ys] = 255

    xs, ys = circle(100, 300, 25)
    img[xs, ys] = 255

    xs, ys = circle(200, 350, 45)
    img[xs, ys] = 255

    blobs = blob_dog(img, min_sigma=5, max_sigma=50)
    radius = lambda x: r2 * x[2]
    s = sorted(blobs, key=radius)
    thresh = 5

    b = s[0]
    assert abs(b[0] - 400) <= thresh
    assert abs(b[1] - 130) <= thresh
    assert abs(radius(b) - 5) <= thresh

    b = s[1]
    assert abs(b[0] - 100) <= thresh
    assert abs(b[1] - 300) <= thresh
    assert abs(radius(b) - 25) <= thresh

    b = s[2]
    assert abs(b[0] - 200) <= thresh
    assert abs(b[1] - 350) <= thresh
    assert abs(radius(b) - 45) <= thresh

    assert_raises(ValueError, blob_dog, img3)

 def get_snr(cube, angle_list, y, x, mode, svd_mode, fwhm, ncomp, fmerit):
     if mode=='full':
         frame = pca(cube, angle_list, ncomp=ncomp, full_output=False, 
                     verbose=False, mask_center_px=mask_center_px, 
                     svd_mode=svd_mode)
     elif mode=='annular':
         y_cent, x_cent = frame_center(cube[0])
         annulus_radius = dist(y_cent, x_cent, y, x)
         frame = pca_annulus(cube, angle_list, ncomp, annulus_width, 
                             annulus_radius)
     else:
         raise RuntimeError('Wrong mode.')            
     
     if fmerit=='max':
         yy, xx = draw.circle(y, x, fwhm/2.)
         snr_pixels = [phot.snr_ss(frame, y_, x_, fwhm, plot=False, 
                                   verbose=False) for y_, x_ in zip(yy, xx)]
         return np.max(snr_pixels)
     elif fmerit=='px':
         return phot.snr_ss(frame, y, x, fwhm, plot=False, verbose=False)
     elif fmerit=='mean':
         yy, xx = draw.circle(y, x, fwhm/2.)
         snr_pixels = [phot.snr_ss(frame, y_, x_, fwhm, plot=False, 
                                   verbose=False) for y_, x_ in zip(yy, xx)]                                      
         return np.mean(snr_pixels)

def create_mask(img, num_circles, lo_thickness, hi_thickness, patch_size):
    im = rgb2gray(img)
    m = np.ones_like(im)

    np.random.seed(31415926)
    for i in range(num_circles):
        im_tmp = np.ones_like(m)
        yy = np.random.randint(0, m.shape[0])
        xx = np.random.randint(0, m.shape[1])
        r = np.random.randint(20, m.shape[0] / 2)
        t = np.random.randint(lo_thickness, hi_thickness)
        rro, cco = circle(yy, xx, r, shape=m.shape)
        rri, cci = circle(yy, xx, r - t, shape=m.shape)
        im_tmp[rro, cco] = 0
        im_tmp[rri, cci] = 1
        m[im_tmp == 0] = 0

    # Fix mask border.
    d = patch_size + 1
    m[:d, :] = 1
    m[-d:, :] = 1
    m[:, :d] = 1
    m[:, -d:] = 1

    return m

def test_blob_dog():
    img = np.ones((512, 512))

    xs, ys = circle(400, 130, 5)
    img[xs, ys] = 255

    xs, ys = circle(100, 300, 25)
    img[xs, ys] = 255

    xs, ys = circle(200, 350, 45)
    img[xs, ys] = 255

    blobs = blob_dog(img, min_sigma=5, max_sigma=50)
    area = lambda x: x[2]
    radius = lambda x: math.sqrt(x / math.pi)
    s = sorted(blobs, key=area)
    thresh = 5

    b = s[0]
    assert abs(b[0] - 400) <= thresh
    assert abs(b[1] - 130) <= thresh
    assert abs(radius(b[2]) - 5) <= thresh

    b = s[1]
    assert abs(b[0] - 100) <= thresh
    assert abs(b[1] - 300) <= thresh
    assert abs(radius(b[2]) - 25) <= thresh

    b = s[2]
    assert abs(b[0] - 200) <= thresh
    assert abs(b[1] - 350) <= thresh
    assert abs(radius(b[2]) - 45) <= thresh

    def get_snr(frame, y, x, fwhm, fmerit):
        """
        """
        if fmerit == 'max':
            yy, xx = draw.circle(y, x, fwhm / 2.)
            res = [snr(frame, (x_, y_), fwhm, plot=False, verbose=False,
                       full_output=True)
                   for y_, x_ in zip(yy, xx)]
            snr_pixels = np.array(res)[:, -1]
            fluxes = np.array(res)[:, 2]
            argm = np.argmax(snr_pixels)
            # integrated fluxes for the max snr
            return np.max(snr_pixels), fluxes[argm]

        elif fmerit == 'px':
            res = snr(frame, (x, y), fwhm, plot=False, verbose=False,
                      full_output=True)
            snrpx = res[-1]
            fluxpx = np.array(res)[2]
            # integrated fluxes for the given px
            return snrpx, fluxpx

        elif fmerit == 'mean':
            yy, xx = draw.circle(y, x, fwhm / 2.)
            res = [snr(frame, (x_, y_), fwhm, plot=False, verbose=False,
                       full_output=True) for y_, x_
                   in zip(yy, xx)]
            snr_pixels = np.array(res)[:, -1]
            fluxes = np.array(res)[:, 2]
            # mean of the integrated fluxes (shifting the aperture)
            return np.mean(snr_pixels), np.mean(fluxes)

def connect_extrema(im_pos, target, markers, visualize=False):
	'''
	im_pos : XYZ positions of each point in image formation (n x m x 3)
	'''
	height, width,_ = im_pos.shape
	centroid = np.array(target)

	im_pos = np.ascontiguousarray(im_pos.astype(np.int16))
	cost_map = np.ascontiguousarray(np.zeros([height, width], dtype=np.uint16))

	extrema = dgn.geodesic_map_MPI(cost_map, im_pos, np.array(centroid, dtype=np.int16), 1, 1)
	cost_map = extrema[-1]

	trails = []
	for m in markers:
		trail = dgn.geodesic_trail(cost_map.copy()+(32000*(im_pos[:,:,2]==0)).astype(np.uint16), np.array(m, dtype=np.int16))
		trails += [trail.copy()]
	if visualize:
		cost_map = deepcopy(cost_map)
		circ = circle(markers[0][0],markers[0][1], 5)
		circ = np.array([np.minimum(circ[0], height-1), np.minimum(circ[1], width-1)])
		circ = np.array([np.maximum(circ[0], 0), np.maximum(circ[1], 0)])
		cost_map[circ[0], circ[1]] = 0
		for i,t in enumerate(trails[1:]):
			# embed()
			cost_map[t[:,0], t[:,1]] = 0
			circ = circle(markers[i+1][0],markers[i+1][1], 5)
			circ = np.array([np.minimum(circ[0], height-1), np.minimum(circ[1], width-1)])
			circ = np.array([np.maximum(circ[0], 0), np.maximum(circ[1], 0)])
			cost_map[circ[0], circ[1]] = 0
		return trails, cost_map
	else:
		return trails

 def test_get_largest_region(self):
     test_img = np.zeros((1000,1000), dtype=np.uint8)
     rr, cc = circle(100,100,50)
     test_img[rr,cc] = 1
     rr, cc = circle(500,500,100)
     test_img[rr, cc] = 1
     largest_region = gzapi.get_largest_region(test_img)
     self.assertTupleEqual(largest_region.centroid, (500,500))

def create_ringed_spider_mask(im_shape, ann_out, ann_in=0, sp_width=10,
                              sp_angle=0):
    """
    Mask out information is outside the annulus and inside the spiders (zeros).

    Parameters
    ----------
    im_shape : tuple of int
        Tuple of length two with 2d array shape (Y,X).
    ann_out : int
        Outer radius of the annulus.
    ann_in : int
        Inner radius of the annulus.
    sp_width : int
        Width of the spider arms (3 branches).
    sp_angle : int
        angle of the first spider arm (on the positive horizontal axis) in
        counter-clockwise sense.

    Returns
    -------
    mask : numpy ndarray
        2d array of zeros and ones.

    """
    mask = np.zeros(im_shape)

    s = im_shape[0]
    r = s/2
    theta = np.arctan2(sp_width/2, r)

    t0 = np.array([theta, np.pi-theta, np.pi+theta, np.pi*2 - theta])
    t1 = t0 + sp_angle/180 * np.pi
    t2 = t1 + np.pi/3
    t3 = t2 + np.pi/3

    x1 = r * np.cos(t1) + s/2
    y1 = r * np.sin(t1) + s/2
    x2 = r * np.cos(t2) + s/2
    y2 = r * np.sin(t2) + s/2
    x3 = r * np.cos(t3) + s/2
    y3 = r * np.sin(t3) + s/2

    rr1, cc1 = polygon(y1, x1)
    rr2, cc2 = polygon(y2, x2)
    rr3, cc3 = polygon(y3, x3)

    cy, cx = frame_center(mask)
    rr0, cc0 = circle(cy, cx, min(ann_out, cy))
    rr4, cc4 = circle(cy, cx, ann_in)

    mask[rr0, cc0] = 1
    mask[rr1, cc1] = 0
    mask[rr2, cc2] = 0
    mask[rr3, cc3] = 0
    mask[rr4, cc4] = 0
    return mask

def ring_mask(i,dim=_DIM):
    # ring masks are a series of 50 concentric rings, each dim/100
    # pixels thick, around the center of the dim x dim array
    s = dim/100.
    rmin = (i * s) + _eps
    rmax = (i+1) * s
    c = (dim//2)+1
    mask = np.zeros((dim,dim),dtype=np.bool)
    mask[circle(c,c,rmax)]=True
    mask[circle(c,c,rmin)]=False
    return mask

def add_circles(image, circles, cmap=plt.cm.cool):
    output = gray2rgb(image)
    for r, x, y, accum in circles:
        shape = list(reversed(image.shape))
        xs_outer, ys_outer = circle(y, x, r, shape=shape)
        xs_inner, ys_inner = circle(y, x, r - 3, shape=shape)
        c_img = zeros(image.shape, dtype=bool)
        c_img[ys_outer, xs_outer] = 1
        c_img[ys_inner, xs_inner] = 0
        color_value = cmap(accum, bytes=True)[:-1]
        output[c_img] = color_value
    return output

def test_blob_dog():
    r2 = math.sqrt(2)
    img = np.ones((512, 512))

    xs, ys = circle(400, 130, 5)
    img[xs, ys] = 255

    xs, ys = circle(100, 300, 25)
    img[xs, ys] = 255

    xs, ys = circle(200, 350, 45)
    img[xs, ys] = 255

    blobs = blob_dog(img, min_sigma=5, max_sigma=50)
    radius = lambda x: r2 * x[2]
    s = sorted(blobs, key=radius)
    thresh = 5

    b = s[0]
    assert abs(b[0] - 400) <= thresh
    assert abs(b[1] - 130) <= thresh
    assert abs(radius(b) - 5) <= thresh

    b = s[1]
    assert abs(b[0] - 100) <= thresh
    assert abs(b[1] - 300) <= thresh
    assert abs(radius(b) - 25) <= thresh

    b = s[2]
    assert abs(b[0] - 200) <= thresh
    assert abs(b[1] - 350) <= thresh
    assert abs(radius(b) - 45) <= thresh

    # Testing no peaks
    img_empty = np.zeros((100,100))
    assert blob_dog(img_empty).size == 0

    # Testing 3D
    r = 10
    pad = 10
    im3 = ellipsoid(r, r, r)
    im3 = util.pad(im3, pad, mode='constant')

    blobs = blob_dog(im3, min_sigma=3, max_sigma=10,
                          sigma_ratio=1.2, threshold=0.1)
    b = blobs[0]

    assert b[0] == r + pad + 1
    assert b[1] == r + pad + 1
    assert b[2] == r + pad + 1
    assert abs(math.sqrt(3) * b[3] - r) < 1

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

        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 add_LV(self, scale, drift=10):
        radius = self.lv_radius = self.base_lv_radius * scale
        thick  = self.lv_thick  = self.base_lv_thick * (4+scale)/5
        
        # get center
        shift = np.array([np.random.randn()*drift, radius])
        
        # print shift
        center = np.array([self.v/2, 10]) + shift
        self.lv_center = center
        self.lv_radius_vec = np.array([0,radius])

        # get points
        self.lv_big = d.circle(center[0], center[1], radius)
        self.lv_small = d.circle(center[0], center[1], radius-thick)

	def getMask(self, world, rad):
		'''Get a mask for the robot's world'''
		# Get the center of the robot pose
		x, y = self.rect.center
		# Draw a circle around that location
		xx, yy = draw.circle(x, y, rad, world.shape)
		# Set the index bounds of the world		
		xMin = 0
		xMax = world.shape[0]-1
		yMin = 0
		yMax = world.shape[1]-1
		# Find the points of the circle that exced the index bounds
		xxMin = np.asarray(np.where(xx<xMin))[0]
		xxMax = np.asarray(np.where(xx>xMax))[0]
		yyMin = np.asarray(np.where(yy<yMin))[0]
		yyMax = np.asarray(np.where(yy>yMax))[0]
		# Clip the shape of the circle to the bounds of the world
		xyd = np.concatenate((xxMin,xxMax,yyMin,yyMax))
		xx = np.delete(xx,xyd)
		yy = np.delete(yy,xyd)
		# Make an empty mask same size as the world
		mask = np.zeros(world.shape)
		# Apply the circle to mask of the world
		mask[xx, yy] = 1
		return mask

def pick_bounding_box(im_rgb):
    # global bounding_box
    bounding_box = []

    im = im_rgb.copy()
    im_display = im_rgb.copy()

    # Display instructions
    txt = "Click on 4 points."
    cv2.putText(im_display, txt, (10, 25), cv2.FONT_HERSHEY_DUPLEX, 1, (255, 255, 255))

    # Have user specify region
    cv2.namedWindow("pick_region")
    cv2.setMouseCallback("pick_region", mouse_event, bounding_box)
    while len(bounding_box) < 4:
        # Display circles on corners
        for pt in bounding_box:
            circ = circle(pt[0], pt[1], 5)
            im_display[circ[0], circ[1]] = 255
        # Display lines between points
        # if len(bounding_box) > 1:
            # for i in range(len(bounding_box)-1):
                # pt1 = boundin_box[i]
                # pt2 = boundin_box[i+1]

        cv2.imshow("pick_region", im_display)
        ret = cv2.waitKey(30)
    cv2.destroyWindow("pick_region")

    return np.array(bounding_box)

    def count_with_edge(dot):
        r, c = dot
        mask = np.zeros((height, width), dtype=np.uint8)
        rr, cc = draw.circle(r, c, radius)

        # min_r = np.where(rr >= 0)[0][0]
        # min_c = np.where(cc >= 0)[0][0]
        # min_i = max(min_r, min_c)

        # print(r, c)
        # print(height, width)
        # print(rr)
        # print(cc)
        # max_r = np.where(rr < height)[0][-1]
        # max_c = np.where(cc < width)[0][-1]
        # max_i = min(max_r, max_c)

        # # print(width, height)
        # # print(rr.shape)
        # # print(cc.shape)
        # # print(rr[min_i: max_i])
        # # print(cc[min_i: max_i])
        # rr = rr[min_i: max_i]
        # cc = cc[min_i: max_i]
        # print(rr)
        # print(cc)

        for _r, _c in zip(rr, cc):
            if 0 <= _r < height and 0 <= _c < width:
                mask[_r, _c] = 1

        # mask[rr, cc] = 1
        binary_for_cal = np.where(mask == 1, binary, False)
        return np.count_nonzero(binary_for_cal)

def add_cm(image, centers, color):
    return_image = image.copy()
    for center in centers:
        rr, cc = circle(center[0], center[1], 10)
        return_image[rr, cc, 0] = 1
        return_image[rr, cc, color] = 1
    return return_image

    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 test_sets_all(labels_all, x_center, y_center, radius, feature_filename, testindex, rot=False):# n_rot=1, x_center=0, y_center=0):
#returns (X_test, y_test)
#labels_all: labels [img, x, y] c: class labels c=0 no label
# x_center, y_center, radius: parameters for interesting region trainingspxls are used from
#feature_filename: filename of hdf5 file with dataset 'data' containing feature_array [x,y,img/imgRot,f] 
#testindex: between 0 and k indicate fold used for testing all others are used for training
#rot: True if featuresarray with imgRot is used
    
    #open feature file
    featuresF = h5py.File(feature_filename,'r')

    #global parameters
    n_features = featuresF['data'].shape[3]
    xDim = featuresF['data'].shape[0]
    yDim = featuresF['data'].shape[1]
    
    #array for features of one frame
    features = np.zeros((xDim, yDim, n_features), dtype=np.float32)

    x, y = circle(y_center, x_center, radius)

    X_test = np.zeros((len(x), n_features))
    y_test = np.zeros((len(x),))
    if rot==False:
        features[:,:,:] = featuresF['data'][:,:,testindex,:]
        X_test[:, :] = features[x,y,:]
    else:
        features[:,:,:] = featuresF['data'][:,:,testindex*n_rot,:]
        X_test[:, :] = features[x,y,:]
    y_test[:] = labels_all[testindex, x, y] 
    featuresF.close()
    return (X_test, y_test)

def display_edges(image, g, threshold):
    """Draw edges of a RAG on its image
 
    Returns a modified image with the edges drawn.Edges are drawn in green
    and nodes are drawn in yellow.
 
    Parameters
    ----------
    image : ndarray
        The image to be drawn on.
    g : RAG
        The Region Adjacency Graph.
    threshold : float
        Only edges in `g` below `threshold` are drawn.
 
    Returns:
    out: ndarray
        Image with the edges drawn.
    """
    image = image.copy()
    for edge in g.edges_iter():
        n1, n2 = edge
 
        r1, c1 = map(int, rag.node[n1]['centroid'])
        r2, c2 = map(int, rag.node[n2]['centroid'])
 
        line  = draw.line(r1, c1, r2, c2)
        circle = draw.circle(r1,c1,2)
 
        if g[n1][n2]['weight'] < threshold :
            image[line] = 0,1,0
        image[circle] = 1,1,0

    return image