def test_square_image():
    im = np.zeros((50, 50)).astype(float)
    im[:25, :25] = 1.

    # Moravec
    results = peak_local_max(corner_moravec(im),
                             min_distance=10, threshold_rel=0)
    # interest points along edge
    assert len(results) == 57

    # Harris
    results = peak_local_max(corner_harris(im, method='k'),
                             min_distance=10, threshold_rel=0)
    # interest at corner
    assert len(results) == 1

    results = peak_local_max(corner_harris(im, method='eps'),
                             min_distance=10, threshold_rel=0)
    # interest at corner
    assert len(results) == 1

    # Shi-Tomasi
    results = peak_local_max(corner_shi_tomasi(im),
                             min_distance=10, threshold_rel=0)
    # interest at corner
    assert len(results) == 1

    def average_hough_detections(self, hough_radii, hough_res, num_best=5):
        """
        Smooths `num_best` hough detections with Gaussian and
        computes weighted average across the `num_best` hough
        detections to get more precise center_x, center_y and
        radius of circle
        """

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

        for radius, h in zip(hough_radii, hough_res):
            # For each radius, extract two circles
            h_smooth = skifilt.gaussian_filter(h, sigma=4)
            num_peaks = 1
            peaks = skif.peak_local_max(h, min_distance=40, num_peaks=num_peaks)
            centers.extend(peaks)
            accums.extend(h[peaks[:, 0], peaks[:, 1]])
            radii.extend([radius] * num_peaks)

        h_sum = np.sum([skifilt.gaussian_filter(x, sigma=2)
                        for x in hough_res[np.argsort(accums)[::-1][:num_best]]], axis=0)

        peaks = skif.peak_local_max(h_sum, min_distance=40, num_peaks=num_peaks)

        center_x, center_y = peaks[0]

        max_sel = [np.max(x.ravel()) for x in hough_res[np.argsort(accums)[::-1][:num_best]]]
        radii_sel = [radii[i] for i in np.argsort(accums)[::-1][:num_best]]

        radius = sum([m * r for m, r in zip(max_sel, radii_sel)]) / float(sum(max_sel))

        return center_x, center_y, int(radius)

def find_peaks(img, npeaks=float("inf")):
    threshold = lambda x: 10.0 ** (-x)
    log_th = __MIN_THRESHOLD
    peaks = peak_local_max(img, threshold_rel=threshold(log_th))

    while len(peaks) > npeaks and log_th > __MAX_THRESHOLD:
        log_th -= 0.05
        peaks = peak_local_max(img, threshold_rel=threshold(log_th))

    return peaks

def get_local_maxima(im_array, min_distance, threshold_rel):
    """Return the local maxima."""
    img =  peak_local_max(im_array,
                          indices=False,
                          min_distance=min_distance,
                          threshold_rel=threshold_rel)
    coords =  peak_local_max(im_array,
                             indices=True,
                             min_distance=min_distance,
                             threshold_rel=threshold_rel)
    return img, coords

def find_n_local_minima(image, n=0):
    """Takes a 3D image and creates a local minimum mask, and n minima mask

    """
    # Invert image so we are searching for maximas
    im_inv = np.invert(image)
    
    # Find peaks in each of r,g,b
    minimaR = peak_local_max(im_inv[:,:,0])
    minimaG = peak_local_max(im_inv[:,:,1])
    minimaB = peak_local_max(im_inv[:,:,2])
    
    # each of the minima will have an array of shape (n, 2)
    # where n is the number of minima

    # get 2D dimensions of the image
    y, x, z = image.shape
    zR = np.zeros((y, x))
    zG = np.zeros((y, x))
    zB = np.zeros((y, x))
    nMinima = np.zeros((y, x))

    # Make zR, zG, zB an accumulator space for the minima points
    zR[minimaR[:,0], minimaR[:,1]] = 1
    zG[minimaG[:,0], minimaG[:,1]] = 1
    zB[minimaB[:,0], minimaB[:,1]] = 1
    # Find where peaks in rgb overlap
    minima = zR * zG * zB

    if n != 0:
        i, j = np.nonzero(minima)

        # making sure the number of minima is not exceeded
        if n > np.sum(minima):
            n = np.sum(minima)

        ix = np.random.choice(len(i), n, replace=False)
        
        nMinima[i[ix], j[ix]] = 1

        # Label the masked image
        labelled = label(nMinima)

        #Make an nMinima shape of the image
        markers = np.zeros(image.shape)
        for i in range(z):
            markers[:,:,i] = labelled
        return minima, markers
    else:
        return minima

def test_squared_dot():
    im = np.zeros((50, 50))
    im[4:8, 4:8] = 1
    im = img_as_float(im)

    # Moravec fails

    # Harris
    results = peak_local_max(corner_harris(im))
    assert (results == np.array([[6, 6]])).all()

    # Shi-Tomasi
    results = peak_local_max(corner_shi_tomasi(im))
    assert (results == np.array([[6, 6]])).all()

def stk_to_rois(stk, threshold, min_size, max_window=8, downscale_factor=2):
    thresholded_stk = stk > threshold
    thresholded_stk = remove_small_objects(thresholded_stk, min_size)
    distance = ndi.distance_transform_edt(thresholded_stk)
    cropped_stk = stk.copy()
    cropped_stk[np.logical_not(thresholded_stk)] = 0
    combined_stk = cropped_stk + distance/distance.max()
    local_max = peak_local_max(combined_stk, indices=False, 
                               footprint=np.ones((max_window, max_window)), 
                               labels=thresholded_stk)
    markers = ndi.label(local_max)[0]
    labels = watershed(-combined_stk, markers, mask=thresholded_stk)
    new_markers = markers.copy()
    for i in set(labels.flatten()):
        if i == 0: continue
        if np.sum(labels==i) < min_size:
            new_markers[markers==i] = 0
    labels = watershed(-combined_stk, new_markers, mask=thresholded_stk)
    labels_set = set(labels.flatten())
    rois = []
    for label in labels_set:
        if label == 0: continue
        if np.sum((labels==label).astype(int)) < min_size: continue
        nroi = np.zeros((stk.shape[0], stk.shape[1]))
        cx,cy = np.where(labels==label)
        cx,cy = int(cx.mean()), int(cy.mean())
        x,y = np.ogrid[0:nroi.shape[0], 0:nroi.shape[1]]
        r = 4
        mask =  (cx-x)**2 + (cy-y)**2 <= r*r
        nroi[mask] = 1
        #nroi[labels==label] = 1
        rois.append(zoom(nroi, downscale_factor, order=0))
    rois = np.array(rois)
    return rois, thresholded_stk, labels

def find_storms(z):
	thresh_z = 20 #min z value for finding peaks
	min_distance = 5 #minimum distance between peaks
	#find peaks, using local maxima
	peaks = feature.peak_local_max(z, threshold_abs=thresh_z,
				min_distance=min_distance, indices=False)

	#uniquely label each peak
	markers, num_markers = ndimage.label(peaks)

	#use watershed algorithm to split image into basins, starting at markers
	labels = watershed(-z, markers, mask=z>thresh_z/2)

	#compute region properties
	props = measure.regionprops(labels, z)

	props = filter_storms(props)
	
	storms = []
	for p in props:
		s = {}
		s['x'], s['y'] = p.centroid
		s['centroid'] = p.centroid[::-1]
		s['max_intensity'] = p.max_intensity
		s['majlen'] = p.major_axis_length
		s['minlen'] = p.minor_axis_length
		s['angle'] = 180 - np.rad2deg(p.orientation)
		s['area'] = p.area
		storms.append(s)

	return storms

def segmentationize(imageSe):
    """
    Divides coherent forms of an image in smaller groups of type integer.
    """
    
    # create an matrix of distances to the next sourrounding area
    distance = ndimage.distance_transform_edt(imageSe, sampling=3)
    erosed = ndimage.binary_erosion(imageSe, iterations=8).astype(imageSe.dtype)
    distanceE = ndimage.distance_transform_edt(erosed, sampling=3)
    distance += (2 * distanceE)
    labels, num = label(imageSe, background=0, return_num='True')
    sizes_image = ndimage.sum(imageSe, labels, range(num))
    sizes_image = np.sort(sizes_image, axis=None)
    pos = int(0.4 * num)
    areal = int(sizes_image[pos] ** 0.5)
    if areal <= 10:
        areal = 10
    elif (areal % 2) != 0:
        areal += 1
    footer = circarea(areal) # draw circle area
    
    # find the positions of the maxima from the distances
    local_maxi = peak_local_max(distance, indices=False, footprint=footer, labels=imageSe)
    markers = label(local_maxi)
    
    # watershed algorithm starts at the maxima and returns labels of particles
    simplefilter("ignore", FutureWarning)   # avoid warning in watershed method
    labels_ws = watershed(-distance, markers, mask=imageSe)
    simplefilter("default", FutureWarning)
    
    return labels, labels_ws, local_maxi

def label_nuclei(binary, min_size):
    '''Label, watershed and remove small objects'''

    distance = medial_axis(binary, return_distance=True)[1]

    distance_blured = gaussian_filter(distance, 5)

    local_maxi = peak_local_max(distance_blured, indices=False, labels=binary, min_distance = 30)

    markers = measure_label(local_maxi)

#    markers[~binary] = -1

#    labels_rw = segmentation.random_walker(binary, markers)

#    labels_rw[labels_rw == -1] = 0

#    labels_rw = segmentation.relabel_sequential(labels_rw)

    labels_ws = watershed(-distance, markers, mask=binary)

    labels_large = remove_small_objects(labels_ws,min_size)

    labels_clean_border = clear_border(labels_large)

    labels_from_one = relabel_sequential(labels_clean_border)

#    plt.imshow(ndimage.morphology.binary_dilation(markers))
#    plt.show()

    return labels_from_one[0]

def corner_peaks(image, min_distance=10, threshold_abs=0, threshold_rel=0.1,
                 exclude_border=True, indices=True, num_peaks=np.inf,
                 footprint=None, labels=None):
    """Find corners in corner measure response image.

    This differs from `skimage.feature.peak_local_max` in that it suppresses
    multiple connected peaks with the same accumulator value.

    Parameters
    ----------
    See `skimage.feature.peak_local_max`.

    Returns
    -------
    See `skimage.feature.peak_local_max`.

    Examples
    --------
    >>> from skimage.feature import peak_local_max, corner_peaks
    >>> response = np.zeros((5, 5))
    >>> response[2:4, 2:4] = 1
    >>> response
    array([[ 0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  1.,  1.,  0.],
           [ 0.,  0.,  1.,  1.,  0.],
           [ 0.,  0.,  0.,  0.,  0.]])
    >>> peak_local_max(response, exclude_border=False)
    array([[2, 2],
           [2, 3],
           [3, 2],
           [3, 3]])
    >>> corner_peaks(response, exclude_border=False)
    array([[2, 2]])
    >>> corner_peaks(response, exclude_border=False, min_distance=0)
    array([[2, 2],
           [2, 3],
           [3, 2],
           [3, 3]])

    """

    peaks = peak_local_max(image, min_distance=min_distance,
                           threshold_abs=threshold_abs,
                           threshold_rel=threshold_rel,
                           exclude_border=exclude_border,
                           indices=False, num_peaks=num_peaks,
                           footprint=footprint, labels=labels)
    if min_distance > 0:
        coords = np.transpose(peaks.nonzero())
        for r, c in coords:
            if peaks[r, c]:
                peaks[r - min_distance:r + min_distance + 1,
                      c - min_distance:c + min_distance + 1] = False
                peaks[r, c] = True

    if indices is True:
        return np.transpose(peaks.nonzero())
    else:
        return peaks

def segment(image, thresh):
    #preprocess image
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    
    #perform euclidean distance transform
    distances = ndimage.distance_transform_edt(thresh)
    localMax = peak_local_max(distances, indices = False, min_distance = 3, labels = thresh)

    #perform connected component analysis on local peaks
    markers = ndimage.label(localMax, structure = np.ones((3, 3)))[0]
    labels = watershed(-distances, markers, mask = thresh)

    #loop over labels returned from watershed to mark them
    for label in np.unique(labels):
        if label == 0:
            continue

        mask = np.zeros(gray.shape, dtype="uint8")
        mask[labels == label] = 255

        #find contours in mask and choose biggest contour by area
        contours = cv2.findContours(mask.copy(), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)[-2]
        contour = max(contours, key = cv2.contourArea)

        #draw circle around max size contour
        ((x, y), r) = cv2.minEnclosingCircle(contour)
        cv2.circle(image, (int(x), int(y)), int(r), (0, 255, 0), 2)
    
    #show final image
    cv2.imshow("Output", image)
    return len(np.unique(labels) - 1)

    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 test_subpix():
    img = np.zeros((50, 50))
    img[:25,:25] = 255
    img[25:,25:] = 255
    corner = peak_local_max(corner_harris(img), num_peaks=1)
    subpix = corner_subpix(img, corner)
    assert_array_equal(subpix[0], (24.5, 24.5))

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 black_background(image, kernel):

    shifted = cv2.pyrMeanShiftFiltering(image, 10, 39)
    gray = cv2.cvtColor(shifted, cv2.COLOR_BGR2GRAY)
    thresh = cv2.threshold(gray, 0, 255,
                           cv2.THRESH_BINARY | cv2.THRESH_OTSU)[1]
    D = ndimage.distance_transform_edt(thresh)
    localMax = peak_local_max(D, indices=False, min_distance=10,
                              labels=thresh)

    # perform a connected component analysis on the local peaks,
    # using 8-connectivity, then appy the Watershed algorithm
    markers = ndimage.label(localMax, structure=np.ones((3, 3)))[0]
    labels = watershed(-D, markers, mask=thresh)
    # create a mask
    mask2 = np.zeros(gray.shape, dtype="uint8")
    #  loop over the unique labels returned by the Watershed  algorithm for
    for label in np.unique(labels):
        # if the label is zero, we are examining the 'background' so simply ignore it
        if label == 0:
            continue
        # otherwise, allocate memory for the label region and draw
        # it on the mask
        mask2[labels == label] = 255
    return mask2

def LargestWatershedRegion(shapes,dims,skipBias): 
    L=len(shapes)-skipBias     
    shapes=shapes.reshape((-1,) + dims[1:]) 
    D=len(dims)
    num_peaks=4
#    structure=np.ones(tuple(3*np.ones((np.ndim(shapes)-1,1))))
    for ll in range(L): 
        temp=shapes[ll]
        local_maxi = peak_local_max(gaussian_filter(temp,[1]*(D-1)), exclude_border=False, indices=False, num_peaks=num_peaks)
        markers,junk = label(local_maxi)
        nonzero_mask=temp>0
        if np.sum(nonzero_mask)>(3**3)*num_peaks:
            labels = watershed(-temp, markers, mask=nonzero_mask) #watershed regions
            ind = 1
            temp2 = np.copy(temp)
            temp2[labels!=1]=0
            total_intensity = sum(temp2.reshape(-1,))
            for kk in range(2,labels.max()+1):
                temp2 = np.copy(temp)
                temp2[labels!=kk]=0
                total_intensity2 = sum(temp2.reshape(-1,))
                if total_intensity2>total_intensity:
                    ind = kk
                    total_intensity=total_intensity2
            temp[labels!=ind]=0
            shapes[ll]=temp
    shapes=shapes.reshape((len(shapes),-1))
    return shapes

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_subpix_dot():
    img = np.zeros((50, 50))
    img[25, 25] = 255
    corner = peak_local_max(corner_harris(img),
                            min_distance=10, threshold_rel=0, num_peaks=1)
    subpix = corner_subpix(img, corner)
    assert_array_equal(subpix[0], (25, 25))

    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