def CannyFilter(time,sep,rawForce,Meta,tauMultiple=25,**kwargs):
        """
        Uses a canny filter (from image processing) to detect change points

        Args:
            time,sep,force,Meta,tauMultiple: see ZScoreByDwell
        Returns:
            0/1 matrix if we think there is or isnt an array
        """
        force = FilterToTau(time,tauMultiple,rawForce)
        n = force.size
        minV = np.min(force)
        maxV = np.max(force)
        # normalize the force 
        force -= minV
        force /= (maxV-minV)
        gradV = np.gradient(force)
        stdev = np.std(gradV)
        mu = np.mean(gradV)
        # convert the force to an image, to make Canny happy
        im = np.zeros((n,3))
        im[:,1] = force
        # set the 'sigma' value to our filtering value
        sigma = tauMultiple
        edges1 = feature.canny(im,sigma=sigma,low_threshold=0.8,
                               high_threshold=(1-10/n),use_quantiles=True)
        edges2 = feature.canny(im * -1,sigma=sigma,low_threshold=0.8,
                               high_threshold=(1-10/n),use_quantiles=True)
        # get where the algorithm thinks a transtition is happenening
        idx = np.where((edges1 == True) | (edges2 == True))[0]
        idx = WalkEventIdx(rawForce,idx)
        # switch canny to be between -0.5 and 0.5
        return idx - 0.5

def edge(ifile, ofile):
    img = io.imread(ifile, flatten=True)
    edges1 = feature.canny(img)
    edges2 = feature.canny(img, sigma=VALUE_SIGMA)
    out = np.uint8(edges2 * 255)

    io.imsave(ofile, out)

    # display results
    fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(8, 3), sharex=True, sharey=True)

    ax1.imshow(img, cmap=plt.cm.jet)
    ax1.axis('off')
    ax1.set_title('noisy image', fontsize=20)

    ax2.imshow(edges1, cmap=plt.cm.gray)
    ax2.axis('off')
    ax2.set_title('Canny filter, $\sigma=1$', fontsize=20)

    ax3.imshow(edges2, cmap=plt.cm.gray)
    ax3.axis('off')
    ax3.set_title('Canny filter, $\sigma=' + str(VALUE_SIGMA) + '$', fontsize=20)

    fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9,
                        bottom=0.02, left=0.02, right=0.98)

    plt.show()

def my_canny(img, fn = None, sigma=6, with_mask=False, save=False, show=False):
    height = img.shape[0]
    width = img.shape[1]
    if with_mask:
        import numpy as np
        mymask = np.zeros((height, width),'uint8')
        y1, x1 = 200, 150
        y2, x2 = 500, 350
        mymask[y1: y2, x1: x2] = 1
        ret = canny(img, sigma=sigma, mask=mymask)
    else:
        ret = canny(img, sigma)
        
    if show:
        from src.utils.io import showimage_pil
        showimage_pil(255*ret.astype('uint8'))
        
    if save:
        from src.utils.io import saveimage_pil
        if with_mask:
            feature = '_sigma' + str(sigma) + '_mask'
        else:
            feature = '_sigma' + str(sigma)


        saveimage_pil(255*ret.astype('uint8'), fn+feature+'.jpg',show=False)
    return ret

def compare(file1, file2):
    image1 = io.imread(file1, as_grey = True)
    image2 = io.imread(file2, as_grey = True)
    image1 = feature.canny(image1)
    image2 = feature.canny(image2)

    return ssim(image1, image2)

def worker(input_file_path, queue):
    int_values = []
    path_sections = input_file_path.split("/")

    for string in path_sections:
        if re.search(r'\d+', string) is not None:
            int_values.append(int(re.search(r'\d+', string).group()))

    file_count = int_values[-1]

    image = cv2.imread(input_file_path, cv2.CV_LOAD_IMAGE_GRAYSCALE)
    edges = img_as_ubyte(canny(image, sigma=canny_sigma))
    img_bw = cv2.threshold(edges, 250, 255, cv2.THRESH_BINARY)[1]

    point = _find_bottom_edge(img_bw)

    try:
        distance = len(img_bw) - point[1]
    except TypeError:
        try:
            edges = img_as_ubyte(canny(image, sigma=canny_sigma_closeup))
            img_bw = cv2.threshold(edges, 250, 255, cv2.THRESH_BINARY)[1]

            distance = len(img_bw) - point[1]
        except TypeError:
            distance = 0

    output = str(file_count) + ":" + str(distance) + "\n"
    queue.put(output)

    return output

def findSigma(image):
	sig = 3.5
	total = len(image)*len(image[0])
	flag = True
	cnt = 0
	while(flag):
		cnt = cnt+1
		edges = feature.canny(image, sigma=sig)
		edSum = np.sum(edges)
		tmp = total/edSum
		print sig

		###if there are too many pixels, increase sig
		if tmp<200:
			sig = sig + .13
		###too few pixels, decr sig
		if tmp>401:
			sig = sig - .13
		elif tmp>200 and tmp <401:
			return edges

		##sometimes any sigma we put in will be incorct so we let feature decide after some trying
		elif cnt>10 and tmp == 0:
			edges = feature.canny(image)
			return edges

def edge(ifile, ofile):
    img = io.imread(ifile, flatten=True)
    edges1 = feature.canny(img)
    edges2 = feature.canny(img, sigma=VAL_SIGMA)
    out = np.uint8(edges2 * 255)

    io.imsave(ofile, out)
    return ofile

def manipulate():
    edges_sigma = 1.75

    while len(framesToConvert) > 0:
        currentFrame = framesToConvert.pop(-1)

        if not os.path.isdir(currentFrame.basePath + "/raw/"):
            os.mkdir(currentFrame.basePath + "/raw/")
            os.mkdir(currentFrame.basePath + "/grayscale/")
            os.mkdir(currentFrame.basePath + "/cropped/")
            os.mkdir(currentFrame.basePath + "/resize150/")
            os.mkdir(currentFrame.basePath + "/edges_" + str(edges_sigma) + "/")

        cv2.imwrite(currentFrame.basePath + "/raw/FRAME_" + str(currentFrame.count) + ".jpg", currentFrame.image)

        grayscaleImage = cv2.cvtColor(currentFrame.image, cv2.COLOR_BGR2GRAY)
        cv2.imwrite(currentFrame.basePath + "/grayscale/FRAME_" + str(currentFrame.count) + ".jpg", grayscaleImage)

        croppedImage = grayscaleImage[0:RAW_IMAGE_HEIGHT, ((RAW_IMAGE_WIDTH - RAW_IMAGE_HEIGHT) / 2):RAW_IMAGE_WIDTH - ((RAW_IMAGE_WIDTH - RAW_IMAGE_HEIGHT) / 2)]
        cv2.imwrite(currentFrame.basePath + "/cropped/FRAME_" + str(currentFrame.count) + ".jpg", croppedImage)

        resizedImage = cv2.resize(croppedImage, (FINAL_IMAGE_SIZE, FINAL_IMAGE_SIZE), interpolation=cv2.INTER_AREA)
        cv2.imwrite(currentFrame.basePath + "/resize150/FRAME_" + str(currentFrame.count) + ".jpg", resizedImage)

        edges = img_as_ubyte(canny(resizedImage, sigma=1.75))
        cv2.imwrite(currentFrame.basePath + "/edges_" + str(edges_sigma) + "/FRAME_" + str(currentFrame.count) + ".jpg", edges)

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 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 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 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 find_edges(img, sigma = 4):
    img = feature.canny(img, sigma)

    selem = disk(10)
    img = dilation(img, selem)

    return img

    def test_tif_series_2otsu_Y23(self):
        frame_list = glob("./test_resources/Y23_0.4X_raw_tifs/*.tif")
        for idx in range(400,450):
            frame = frame_list[idx]
            self.image = np.array(Image.open(frame))
            self.thresholds = otsu(self.image, nclasses=2)
            if idx==425:
                # Plot intensity histogram and thresholds
                plt.hist(self.image.flatten(), bins=256, log=True)
                for t in self.thresholds:
                    plt.axvline(t, color='r')
                plt.title('Y23 using class=2')
                plt.draw()
                plt.savefig('Y23_2-class.png', dpi=300)
                plt.clf()
                Image.open('Y23_2-class.png').show()

                # Plot actual data and thresholded
                edge = feature.canny(self.image>self.thresholds[0])
                highlight = np.copy(self.image)
                highlight[edge] = 2*self.image.max() - self.image.min()
                fig = plt.figure()
                ax11 = fig.add_subplot(221)
                ax12 = fig.add_subplot(222)
                ax21 = fig.add_subplot(223)
                ax22 = fig.add_subplot(224)
                ax11.imshow(self.image)
                ax12.imshow(self.image>self.thresholds[0])
                ax21.imshow(highlight)
                ax22.imshow(edge)
                plt.show()

    def _calc_crispness(self, grey_array):
        """Calculate three measures of the crispness of an channel.

        PARAMETERS
        ----------
        grey_array : 2D numpy array
            Raw data for the grey channel.

        PRODUCES
        --------
        crispnesses : list
            Three measures of the crispness in the grey channel of types:

            - ``sobel``, ``canny``, and ``laplace``
        """
        grey_array = grey_array/255
        sobel_var = filters.sobel(grey_array).var()
        canny_array = feature.canny(grey_array, sigma=1).var()
        canny_ratio = np.sum(canny_array == True)/float(
                                                    len(canny_array.flatten()))
        laplace_var = filters.laplace(grey_array, ksize=3).var()
        self.feature_data.extend([sobel_var, canny_ratio, laplace_var])
        if self.columns_out:
            self.column_names.extend(['crisp_sobel', 'crisp_canny',
                                      'crisp_laplace'])

def main():
    for file_path in glob.glob("/home/lucas/Downloads/Lucas/GSK 10uM/*.JPG"):

        img = data.imread(file_path, as_grey=True)

        img = transform.resize(img, [600, 600])
        img_color = transform.resize(data.imread(file_path), [600, 600])

        img[img >img.mean()-0.1] = 0

        # io.imshow(img)
        # io.show()
        #
        edges = canny(img)
        bordas_fechadas = closing(img > 0.1, square(15)) # fechando gaps
        fill_cells = ndi.binary_fill_holes(bordas_fechadas)
        # io.imshow(fill_cells)
        # io.show()
        img_label = label(fill_cells, background=0)
        n= 0
        for  x in regionprops(img_label):
            if x.area < 2000 and x.area > 300:
                n +=1
                print x.area
                minr, minc, maxr, maxc = x.bbox
                try:
                    out_path_name = file_path.split("/")[-1].rstrip(".JPG")
                    io.imsave("out/cell_{}_pic_{}_area_{}.png".format(n, out_path_name, str(round(x.area))),img_color[minr-3: maxr+3, minc-3: maxc+3])
                    #io.show()
                except:
                    pass

def getBoundary(img, debug=False, **kwds):
    from skimage import feature
    edge = feature.canny(img, **kwds)
    start_row = None
    middle_col = (edge.shape[1]-1)//2
    start_cols = np.ones(edge.shape[0], dtype=int)*middle_col
    stop_cols = np.ones(edge.shape[0], dtype=int)*middle_col
    for i, row in enumerate(edge):
        isedge = row > 0
        if isedge.any():
            w = np.where(isedge)[0]
            start_cols[i], stop_cols[i] = w[0], w[-1]
            # set the row that starts to have object to be measured
            if start_row is None:
                start_row = i
            stop_row = i+1
        continue
    if debug:
        print(start_row, stop_row)
        # for start, stop in zip(start_cols, stop_cols):
        #    print(start, stop)
        #    continue
        from matplotlib import pyplot as plt
        plt.figure()
        plt.plot(start_cols, '.')
        plt.plot(stop_cols, '.')
        plt.savefig("edge.png")
        plt.close()
    return start_row, stop_row, start_cols, stop_cols

    def edges(cls):
        from scipy import ndimage, misc
        import numpy as np
        from skimage import feature
        col = Image.open("f990.jpg")
        gray = col.convert('L')

        # Let numpy do the heavy lifting for converting pixels to pure black or white
        bw = np.asarray(gray).copy()

        # Pixel range is 0...255, 256/2 = 128
        bw[bw < 245]  = 0    # Black
        bw[bw >= 245] = 255 # White
        bw[bw == 0] = 254
        bw[bw == 255] = 0
        im = bw
        im = ndimage.gaussian_filter(im, 1)
        edges2 = feature.canny(im, sigma=2)
        labels, numobjects =ndimage.label(im)
        slices = ndimage.find_objects(labels)
        print('\n'.join(map(str, slices)))
        misc.imsave('f990_sob.jpg', im)
        return

        #im = misc.imread('f990.jpg')
        #im = ndimage.gaussian_filter(im, 8)
        sx = ndimage.sobel(im, axis=0, mode='constant')
        sy = ndimage.sobel(im, axis=1, mode='constant')
        sob = np.hypot(sx, sy)
        misc.imsave('f990_sob.jpg', edges2)