def detectOpticDisc(image):
    kernel = octagon(10, 10)
    thresh = threshold_otsu(image[:,:,1])
    binary = image > thresh
    print binary.dtype
    luminance = convertToHLS(image)[:,:,2]
    t = threshold_otsu(luminance)
    t = erosion(luminance, kernel)
    
    
    labels = segmentation.slic(image[:,:,1], n_segments = 3)
    out = color.label2rgb(labels, image[:,:,1], kind='avg')
    skio.imshow(out)
    
    x, y = computeCentroid(t)
    print x, y
    rows, cols, _ = image.shape
    p1 = closing(image[:,:,1],kernel)
    p2 = opening(p1, kernel)
    p3 = reconstruction(p2, p1, 'dilation')
    p3 = p3.astype(np.uint8)
    #g = dilation(p3, kernel)-erosion(p3, kernel)
    #g = rank.gradient(p3, disk(5))
    g = cv2.morphologyEx(p3, cv2.MORPH_GRADIENT, kernel)
    #markers = rank.gradient(p3, disk(5)) < 10
    markers = drawCircle(rows, cols, x, y, 85)
    #markers = ndimage.label(markers)[0]
    #skio.imshow(markers)
    g = g.astype(np.uint8)
    #g = cv2.cvtColor(g, cv2.COLOR_GRAY2RGB)
    w = watershed(g, markers)
    print np.max(w), np.min(w)
    w = w.astype(np.uint8)
    #skio.imshow(w)
    return w

def run(args):
    logging.basicConfig(level=logging.INFO)

    slide = openslide.OpenSlide(args.wsi_path)

    # note the shape of img_RGB is the transpose of slide.level_dimensions
    img_RGB = np.transpose(np.array(slide.read_region((0, 0),
                           args.level,
                           slide.level_dimensions[args.level]).convert('RGB')),
                           axes=[1, 0, 2])

    img_HSV = rgb2hsv(img_RGB)

    background_R = img_RGB[:, :, 0] > threshold_otsu(img_RGB[:, :, 0])
    background_G = img_RGB[:, :, 1] > threshold_otsu(img_RGB[:, :, 1])
    background_B = img_RGB[:, :, 2] > threshold_otsu(img_RGB[:, :, 2])
    tissue_RGB = np.logical_not(background_R & background_G & background_B)
    tissue_S = img_HSV[:, :, 1] > threshold_otsu(img_HSV[:, :, 1])
    min_R = img_RGB[:, :, 0] > args.RGB_min
    min_G = img_RGB[:, :, 1] > args.RGB_min
    min_B = img_RGB[:, :, 2] > args.RGB_min

    tissue_mask = tissue_S & tissue_RGB & min_R & min_G & min_B

    np.save(args.npy_path, tissue_mask)

 def evaluate(self,img):
   img = cv2.medianBlur(img,5)
   if self.auto_thresh:
     self.thresh = threshold_otsu(img)
   if self.white_spots:
     bw = (img > self.thresh).astype(np.uint8)
   else:
     bw = (img <= self.thresh).astype(np.uint8)
   #cv2.imshow(self.name,bw*255)
   #cv2.waitKey(5)
   if not .1* img.size < np.count_nonzero(bw) < .8*img.size:
     print("reevaluating threshold!!")
     print("Ratio:",np.count_nonzero(bw)/img.size)
     print("old:",self.thresh)
     self.thresh = threshold_otsu(img)
     print("new:",self.thresh)
   m = cv2.moments(bw)
   r = {}
   try:
     r['x'] = m['m10']/m['m00']
     r['y'] = m['m01']/m['m00']
   except ZeroDivisionError:
     return -1
   x,y,w,h = cv2.boundingRect(bw)
   #if (h,w) == img.shape:
   #  return -1
   r['bbox'] = y,x,y+h,x+w
   return r

 def duel(frame):
     threshold = threshold_otsu(frame)
     binary1 = frame > threshold
     frame = frame - binary1 * frame
     threshold = threshold_otsu(frame)
     binary2 = frame > threshold
     binary = binary1 & binary2
     return binary

def watershed_counter(path_image):
    # load the image and convert it to a floating point data type
    rgb_image = img_as_float(io.imread(path_image))
    image = rgb2grey(rgb_image)

    bin_image = image > threshold_otsu(image)
    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 2.5))
    ax1.imshow(image, cmap=plt.cm.gray)
    ax1.set_title('Original Image')
    ax1.axis('off')

    ax2.hist(image)
    ax2.set_title('Otsu Thresholded Histogram')
    ax2.axvline(threshold_otsu(image), color='r')

    ax3.imshow(bin_image, cmap=plt.cm.gray)
    ax3.set_title('Thresholded Image')
    ax3.axis('off')

    # Now we want to separate the two objects in image
    # Generate the markers as local maxima of the distance to the background

    distance = ndi.distance_transform_edt(bin_image)
    local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((3, 3)), labels=bin_image)
    markers = ndi.label(local_maxi)[0]
    labels = watershed(distance, markers, mask=bin_image)

    regions = regionprops(labels)
    regions = [r for r in regions if r.area > 50]
    num = len(regions)

    fig, axes = plt.subplots(ncols=4, figsize=(8, 2.7))
    ax0, ax1, ax2, ax3 = axes

    ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
    ax0.set_title('Overlapping objects')
    ax1.imshow(-distance, cmap=plt.cm.jet, interpolation='nearest')
    ax1.set_title('Distances')
    ax2.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest')
    ax2.set_title(str(num) + ' Total Objects')
    ax3.imshow(rgb_image, cmap=plt.cm.gray, interpolation='nearest')
    ax3.contour(labels, [0.5], linewidths=1.2, colors='y')
    ax3.axis('off')

    for ax in axes:
        ax.axis('off')

    fig.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
                        right=1)

    print num

    plt.show()

 def initial_segmentation(self, emphasized):
     binary_img = np.zeros(emphasized.shape)
     for z in range(0, emphasized.shape[0]):
         th = threshold_otsu(emphasized[z, :, :])
         binary_img[z, :, :] = emphasized[z, :, :] > th
         tmp = np.multiply(emphasized[z, :, :], binary_img[z, :, :])
         tmp[tmp < 0] = 0
         th = threshold_otsu(tmp)
         binary_img[z, :, :] = tmp > th
         binary_img[z, :, :] = remove_small_objects(binary_img[z, :, :].astype(bool), 5)
     binary_img = binary_img.astype('i4')
     return binary_img

    def CleanImgage(self, thresh = 30):

        # Image processing

        footprint = np.array([[-1,-1,-1],[-1,8,-1], [-1,-1,-1]])
        self.clean_image = cv2.medianBlur(self.image, 5)
        self.clean_image = cv2.filter2D(self.clean_image,-1,footprint)
        self.clean_image = cv2.medianBlur(self.clean_image, 5)
        self.markers = np.zeros_like(self.image)
        self.markers[self.clean_image < threshold_otsu(self.image)] = 1
        self.markers[self.clean_image >= ((threshold_otsu(self.image)*thresh)/100)] = 2
        self.markers[self.clean_image >= ((threshold_otsu(self.image)*50)/100)] = 3

def load_data(size):
    """
    Loads data from /chars and /notChars to train binary Classifier
    """
    char_files = glob.glob("/users/derekjanni/pyocr/chars/*")
    non_files = glob.glob("/users/derekjanni/pyocr/notChars/*")
    # pre-process character data
    X_char = []
    for i in char_files:
        img = Image.open(i)
        img = img.convert("L")
        img = img.resize((size, size))
        image = np.asarray(img).astype("int32")
        image.setflags(write=True)
        thresh = threshold_otsu(image)
        binary = image > thresh
        X_char.append(binary)

    Y_char = [1 for i in range(len(X_char))]
    # pre-process non-char data
    X_non = []
    for i in char_files:
        img = Image.open(i)
        img = img.convert("L")
        img = img.resize((size, size))
        image = np.asarray(img).astype("int32")
        image.setflags(write=True)
        thresh = threshold_otsu(image)
        binary = image > thresh
        X_non.append(binary)

    # add 1000 circles of varying radius and position
    for i in range(1000):
        radius = np.random.randint(10, 15) + np.random.randint(5)
        offset = np.random.randint(-5, 5)
        img = np.ones((size, size), dtype=np.int32)
        cv2.circle(img, ((offset + size) / 2, (offset + size) / 2), radius, (0, 0, 0), -1)
        X_non.append(img)

    Y_non = [0 for i in range(len(X_non))]
    X = np.asarray(X_char + X_non).reshape(-1, 1, 50, 50).astype(np.float32)
    Y = np.asarray(Y_char + Y_non).astype(np.int32)
    # some trickery to emulate a train-test split
    indices = list(zip(X, Y))  # name to indicate what's going on here, not the actual data struct
    random.shuffle(indices)
    X, Y = zip(*indices)
    X = np.asarray(X)
    Y = np.asarray(Y)
    return X, Y

def read_input_otsu():
	print("Otsu")
	alphabet=["a","b","c","d","e","f","g","h","i","j","k","l","m","n","o","p","q","r","s","t","u","v","w","x","y","z"]
	#alphabet=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,17,19,20,21,2]
	n=70
	train_data=[]
	train_solutions=[]
	for i in range(n):
		for letter in alphabet:
			image = imread("chars74k-lite/chars74k-lite/"+letter+"/"+letter+"_"+str(i)+".jpg")
			img = np.asarray(image)
			img.setflags(write=True)
			thresh = threshold_otsu(img)
			binary = image > thresh
			binary = np.reshape(binary,(1,400))
			train_data.append(np.asarray(binary[0]).astype(float))
			#
			t=[0]*26#init all vectors to zero
			t[ord(letter)-97] = 1
			t = np.asarray(t)
			train_solutions.append(t)
	#print(train_data[0])
	#train_data = threshold_otsu(train_data)
	#print(train_data[0])
	#
	m=88
	test_data=[]
	test_solutions=[]
	for j in range(n,m):
		for letter in alphabet:
			image = imread("chars74k-lite/chars74k-lite/"+letter+"/"+letter+"_"+str(j)+".jpg")
			img = np.asarray(image)
			img.setflags(write=True)
			thresh = threshold_otsu(img)
			binary = image > thresh
			binary = np.reshape(binary,(1,400))
			test_data.append(np.asarray(binary[0]).astype(float))
			#
			t=[0]*26#init all vectors to zero
			t[ord(letter)-97] = 1
			t = np.asarray(t)
			test_solutions.append(t)
	#
	'''print(len(train_data))
	print(len(train_solutions))
	print(len(test_data))
	print(len(test_solutions))'''
	return train_data, train_solutions, test_data, test_solutions

    def locate(i):
        """
        Median subtract each hologram, convolve with Mexican hat kernel,
        then smooth the absolute value of the convolution, and use
        Otsu's thresholding to segment the image into specimens.
        Record the time, x and y centroids, and some intensity features
        within each segment.
        """
        median_sub_holo = hologram_cube[i, ...] - median_holo
        conv_holo = convolve_fft(median_sub_holo,
                                 MexicanHat2DKernel(convolution_kernel_radius),
                                 fftn=fft2, ifftn=ifft2)
        smooth_abs_conv = gaussian_filter(np.abs(conv_holo),
                                          gaussian_kernel_radius)

        thresh = threshold_otsu(smooth_abs_conv - np.median(smooth_abs_conv))
        # thresh = threshold_yen(smooth_abs_conv - np.median(smooth_abs_conv))

        masked = np.ones_like(smooth_abs_conv)
        masked[smooth_abs_conv <= thresh] = 0

        label_image = label(masked)
        regions = regionprops(label_image, smooth_abs_conv)
        for region in regions:
            centroid = region.weighted_centroid
            pos = (i, centroid[0], centroid[1],
                   region.max_intensity, region.mean_intensity)
            positions.append(pos)

def cont(imagen, depth=2**16, gaussian=3, screenpercent=0.7,t=0):

    imagen = gaussian_filter(imagen, gaussian)

    if t==0:
        otsu = threshold_otsu(imagen, depth)
    elif t==1:
        otsu = filters.threshold_isodata(imagen, depth)
    else:
        otsu = filters.threshold_li(imagen)
    imagen = binarizar(imagen, otsu)
    imagen = gaussian_filter(imagen, gaussian)

    contours = measure.find_contours(imagen, 1)
    centro = np.asanyarray([1280*0.5, 960*0.5])
    while len(contours) > 1:
        if sum(np.abs(centro - contours[1].mean(axis=0)) < [1280*screenpercent*0.5, 960*screenpercent*0.5]) != 2:
            del contours[1]
        elif sum(np.abs(centro - contours[0].mean(axis=0)) < [1280*screenpercent*0.5, 960*screenpercent*0.5]) != 2:
            del contours[0]
        else:
            if contours[1].size < contours[0].size:
                del contours[1]
            else:
                del contours[0]
    return imagen, contours[0]

    def read_img_bw(self, image_file, array_type="one_zero"):
        '''
        Read an image from a file and return its matrix
        inputs:
            image_file: full path to the image file including extension
            array_type: String with two values:
                "one_zero": matrix is uint8 with only 1s and 0s
                "bool": matrix is bool with True and False
        '''
        #read image and convert to matrix
        image = Image.open(image_file).convert('L')
#         image = image.resize((100,100), Image.LANCZOS)
        image_array = image.getdata()
               
        image_array = numpy.array(image_array).astype(numpy.uint8).reshape((image.size[0],image.size[1]))
#         print image_array.shape
        
        #Threshold the image according to array_type parameter
        if array_type == "bool":
            thresh = threshold_otsu(image_array)
            image_array = image > thresh
            image_array = ~image_array
        else:
            image_array[image_array < 128] = 1
            image_array[image_array >= 128] = 0
            image_array[image_array == 1] = 255
#         plt.imshow(image_array)
#         plt.show()
        return image_array

def mouse_centroid(im, previous_centroid):
    """
    Find mouse's centroid in a single image
    
    Parameters:
        im = image of analyze (numpy array)
        
        previous_centroid = coordinates of the mouse in the previous frame
        
    Returns:
        Coordinates of the mouse's centroid
    """
    original = copy(im)
    im = im < filters.threshold_otsu(im) * 0.2
    distance = ndimage.distance_transform_edt(im)
    centers = (distance > 0.8 * distance.max())
    if len(centers) == 0:
        return previous_centroid
    labels = label(centers)
    centroids = [r.weighted_centroid for r in regionprops(labels, original)]
    if len(centroids) == 1:
        return list(centroids[0])
    elif len(centroids) > 1:
        d = lambda a, b: ((a[0] - b[0])**2 + (a[1] - b[1])**2)**0.5
        dists = [d(c, previous_centroid) for c in centroids]
        d_idx = np.array(dists == np.min(dists))
        return list(np.array(centroids)[d_idx][0])

def max_mask_iter(fns, offset=0, close_radius=0, erode_radius=0):
    """Find masks for a set of images having brightness artifacts.

    Parameters
    ----------
    fns : list of string
        The images being examined.
    offset : int, optional
        Offset the threshold automatically found.
    close_radius : int, optional
        Perform a morphological closing of the mask of this radius.
    erode_radius : int, optional
        Perform a morphological erosion of the mask, after any closing,
        of this radius.

    Returns
    -------
    maxes : iterator of bool array
        The max mask image corresponding to each input image.
    """
    ms = maxes(fns)
    t = imfilter.threshold_otsu(ms)
    ims = it.imap(io.imread, fns)
    masks = ((im < t + offset) for im in ims)
    if close_radius > 0:
        masks = (morphop(mask, 'close', close_radius) for mask in masks)
    if erode_radius > 0:
        masks = (morphop(mask, 'erode', erode_radius) for mask in masks)
    return masks

def im_proc(im):
    """Apply series of morphological procedures on image."""
    th = threshold_otsu(im)
    im_bin = im > th
    return(ndi.binary_fill_holes(
                morphology.closing(
                im_bin,np.ones((3,3)))))

def nln(src_image):
    orig_image = io.imread(src_image, 0)
    if len(orig_image.shape) == 3:
        image = (orig_image.sum(axis=2) / 3).astype('ubyte')
    else:
        image = orig_image
    h, w = image.shape
    thresh = filters.threshold_otsu(image)
    binary = (image < thresh).astype('ubyte') # 二值化，黑色为1
    sum_of_black_pixels = np.sum(binary)
    if sum_of_black_pixels == 0:
        raise Exception("white or black image.")
    binary_line = binary.sum(axis=0) # 在X轴上投影
    hx = binary_line / (sum_of_black_pixels * 1.0)
    binary_line_v = binary.sum(axis=1) # 在Y轴上投影
    hy = binary_line_v / (sum_of_black_pixels * 1.0)

    H = 32
    W = 32
    binary_new = np.zeros( (H, W) )
    for y in range(h):
        for x in range(w):
            if binary[y,x] == 1:
                x2 = int((W-1) * np.sum(hx[0:x+1]))
                y2 = int((H-1) * np.sum(hy[0:y+1]))
                x2_end = int((W - 1) * np.sum(hx[0:x + 2]))
                y2_end = int((H - 1) * np.sum(hy[0:y + 2]))
                if y == h - 1:
                    y2 = H - 1
                    y2_end = H
                if x == w - 1:
                    x2 = W - 1
                    x2_end = W
                binary_new[y2:y2_end, x2:x2_end] = binary[y,x]
    return binary_new.astype('ubyte')

def isolate_single_seed(stack_path, output_path):
    """Load stack then isolate seed."""

    sigma = 10
    iterations = 1

    raw_stack = Image3D.from_path(stack_path)

    print raw_stack.shape

    smoothed = gaussian_filter(raw_stack, sigma).view(Image3D)

    edges = sobel_magnitude_nd(smoothed).view(Image3D)

    #edges.save('edges')

    labels = np.zeros(raw_stack.shape)
    cx, cy, cz = map(lambda x: x/2, raw_stack.shape)
    labels[cx, cy, cz] = 1

    threshold = threshold_otsu(smoothed)
    thresholded = smoothed > threshold

    #thresholded.view(Image3D).save('thresh')

    segmentation = watershed(edges, markers=labels, mask=thresholded)

    #segmentation.view(Image3D).save('seg')

    dilated = binary_dilation(segmentation, iterations=iterations)

    isolate = np.multiply(raw_stack, dilated)

    isolate.view(Image3D).save(output_path)

def image_prep(img):
	# apply filters for better contrast and noise removal
	img_gamma = adjust_gamma(img, 0.7)
	img_median = median(img_gamma, disk(1))

	# apply threshold
	val =  threshold_otsu(img_median)
	img_otsu = img_median > val

	# label image regions
	label_image = label(img_otsu)

	candidates = []
	for region in regionprops(label_image):
		minr, minc, maxr, maxc = region.bbox
		if (maxr - minr > maxc - minc):
			candidates.append(region)

	# find numbers
	areas = []
	for candidate in candidates:
		areas.append(candidate.area)
	areas.sort()
	areas.reverse()

	n = 1
	v = []
	for candidate in candidates:
		if (candidate.area == areas[0] or candidate.area == areas[1] or candidate.area == areas[2]):
			v.append(candidate.image)
			imsave('num%d.png' % n, candidate.image)
			n += 1
	return v

def label_stack(z_stack, parameters):
    '''
    Apply an automated threshold and a watershed algorithm to
    label cells in each planes of the stack

    Parameters:
    -----------
    z_stack: 3d array of shape (nz, nx, ny)

    Returns:
    --------
    labeled_stack: 3d array of the same shape as z_stack,
    with each detected region labeled
    '''
    segment_method = parameters['segment_method']
    correction = parameters['correction']
    labeled_stack = np.zeros(z_stack.shape, dtype=np.uint8)
    max_int_proj = z_stack.max(axis=0)
    thresh = None

    if segment_method == 'otsu':
        thresh = threshold_otsu(max_int_proj) * correction
    elif segment_method == 'naive':
        thresh = max_int_proj.max() * correction
    else:
        err_string = ('Segmentation method {}'
                      'is not implemented'.format(segment_method))
        raise NotImplementedError(err_string)
    while thresh > z_stack.max():
        log.warning('''Reducing threshold''')
        thresh *= 0.9
    for n, frame in enumerate(z_stack):
        labeled_stack[n] = label_from_thresh(frame, thresh, parameters)
    return labeled_stack

def testSkimage2():
    img = Image.open('../img/1.png')
    # img = Image.open('../img/2.png')
    img = np.array(img)
    imggray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    # (thresh, imgbw) = cv2.threshold(imggray, 128, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
    camera = imggray
    # camera = data.camera()
    val = filters.threshold_otsu(camera)

    hist, bins_center = exposure.histogram(camera)

    plt.figure(figsize=(9, 4))
    plt.subplot(131)
    plt.imshow(camera, cmap='gray', interpolation='nearest')
    plt.axis('off')
    plt.subplot(132)
    plt.imshow(camera < val, cmap='gray', interpolation='nearest')
    plt.axis('off')
    plt.subplot(133)
    plt.plot(bins_center, hist, lw=2)
    plt.axvline(val, color='k', ls='--')

    plt.tight_layout()
    plt.show()
    return