def returnProcessedImage(que,folder,img_flist):
	X = []
	for fname in img_flist:
		cur_img = imread(folder+'/'+fname , as_grey=True)
		cur_img = 1 - cur_img

		######## randomly add samples

		# random add contrast
		r_for_eq = random()
		cur_img = equalize_adapthist(cur_img,ntiles_x=8,ntiles_y=8,clip_limit=(r_for_eq+0.5)/3)

		
		#random morphological operation
		r_for_mf_1 = random()
		if 0.05 < r_for_mf_1 < 0.25: # small vessel
			selem1 = disk(0.5+r_for_mf_1)
			cur_img = dilation(cur_img,selem1)
			cur_img = erosion(cur_img,selem1)
		elif 0.25 < r_for_mf_1 < 0.5: # large vessel
			selem2 = disk(2.5+r_for_mf_1*3)
			cur_img = dilation(cur_img,selem2)
			cur_img = erosion(cur_img,selem2)
		elif 0.5 < r_for_mf_1 < 0.75: # exudate
			selem1 = disk(9.21)
			selem2 = disk(7.21)
			dilated1 = dilation(cur_img, selem1)
			dilated2 = dilation(cur_img, selem2)
			cur_img = np.subtract(dilated1, dilated2)
		
		cur_img = img_as_float(cur_img)
		X.append([cur_img.tolist()])
	# X = np.array(X , dtype = theano.config.floatX)
	que.put(X)
	return X

def get_segmentation_features(im):
    dilwindow = [4, 4]
    imthr = np.where(im > np.mean(im), 0.0, 1.0)
    imdil = morphology.dilation(imthr, np.ones(dilwindow))
    labels = measure.label(imdil)
    labels = imthr * labels
    labels = labels.astype(int)
    regions = measure.regionprops(labels)
    numregions = len(regions)
    while len(regions) < 1:
        dilwindow[0] = dilwindow[0] - 1
        dilwindow[1] = dilwindow[1] - 1
        if dilwindow == [0, 0]:
            regions = None
            break
        imthr = np.where(im > np.mean(im), 0.0, 1.0)
        imdil = morphology.dilation(imthr, np.ones(dilwindow))
        labels = measure.label(imdil)
        labels = imthr * labels
        labels = labels.astype(int)
        regions = measure.regionprops(labels)
    regionmax = get_largest_region(regions, labels, imthr)

    if regionmax is None:
        return (np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan)
    eccentricity = regionmax.eccentricity
    convex_area = regionmax.convex_area
    convex_to_total_area = regionmax.convex_area / regionmax.area
    extent = regionmax.extent
    filled_area = regionmax.filled_area
    return (eccentricity, convex_area, convex_to_total_area, extent,
            filled_area, numregions)

def removeNoise(img, r=7):

    # Creating a circle inside an array
    c = r # Center
    d = 2 * r + 1 # Diameter
    y, x = np.ogrid[-c:d-c, -c:d-c] # Create a True/False grid in numpy
    mask = x * x + y * y <= r * r # Circular shape
    structuringElement = np.zeros((d, d))
    structuringElement[mask] = 1 # Fill ones at the places with True

    # Applying erosion at the binary image
    eroded = erosion(img, structuringElement)

    # Dilate the remaining pixels from the eroded image
    dilated = dilation(eroded, structuringElement)

    # We have now opened the image. Now we need to close it:

    # We could have done this in the same step as the last, but we want to show all steps
    dilated2 = dilation(dilated, structuringElement)

    # Then we close by eroding back to normal
    eroded2 = erosion(dilated2, structuringElement)

    return eroded, dilated, dilated2, eroded2

def main():
    _bpcl = buffer_mask(bpcl)
    selem = morphology.disk(2)
    _bpcl = morphology.dilation(_bpcl, selem)
    pcs = make_potential_cloud_shadow_mask(nir, water, _bpcl, 0.02, 0.0005)
    pcs_buff = buffer_mask(pcs)
    # pcs_buffer = buffer_mask(pcs_buff)
    pcs_buffer = morphology.dilation(pcs_buff, selem)
    return _bpcl, pcs_buffer

def extract_binary_masks_from_structural_channel(Y, min_area_size=30, min_hole_size=15, gSig=5, expand_method='closing', selem=np.ones((3, 3))):
    """Extract binary masks by using adaptive thresholding on a structural channel

    Inputs:
    ------
    Y:                  caiman movie object
                        movie of the structural channel (assumed motion corrected)

    min_area_size:      int
                        ignore components with smaller size

    min_hole_size:      int
                        fill in holes up to that size (donuts)

    gSig:               int
                        average radius of cell

    expand_method:      string
                        method to expand binary masks (morphological closing or dilation)

    selem:              np.array
                        morphological element with which to expand binary masks

    Output:
    -------
    A:                  sparse column format matrix
                        matrix of binary masks to be used for CNMF seeding

    mR:                 np.array
                        mean image used to detect cell boundaries
    """

    mR = Y.mean(axis=0)
    img = cv2.blur(mR, (gSig, gSig))
    img = (img - np.min(img)) / (np.max(img) - np.min(img)) * 255.
    img = img.astype(np.uint8)

    th = cv2.adaptiveThreshold(img, np.max(
        img), cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY, gSig, 0)
    th = remove_small_holes(th > 0, min_size=min_hole_size)
    th = remove_small_objects(th, min_size=min_area_size)
    areas = label(th)

    A = np.zeros((np.prod(th.shape), areas[1]), dtype=bool)

    for i in range(areas[1]):
        temp = (areas[0] == i + 1)
        if expand_method == 'dilation':
            temp = dilation(temp, selem=selem)
        elif expand_method == 'closing':
            temp = dilation(temp, selem=selem)

        A[:, i] = temp.flatten('F')

    return A, mR

def get_symbols(image):
  dil_eros = bin_search(dilatation_cross_numb, [image], (1, 16), 1.0, "dec")
  block_size = 50
  binary_adaptive_image = erosion(dilation(threshold_adaptive(
    array(image.convert("L")), block_size, offset=10),
      square(dil_eros)), square(dil_eros))

  all_labels = label(binary_adaptive_image, background = True)
  objects = find_objects(all_labels)

  av_width = av_height = 0
  symbols = []

  for obj in objects:
    symb = (binary_adaptive_image[obj], (obj[0].start, obj[1].start))
    symbols.append(symb)
    av_height += symb[0].shape[0]
    av_width += symb[0].shape[1]

  av_width /= float(len(objects))
  av_height /= float(len(objects))

  symbols = [symb for symb in symbols
    if symb[0].shape[0] >= av_height and symb[0].shape[1] >= av_width]

  return symbols

def get_distorted(image, params, orient = "horizont"):
  shifts = []
  np_image = array(image.convert("L"))
  for el in params:
    if el[0] == "sin":
      shifts.append(lambda x: np_image.shape[0] / el[1] * \
        np.sin(x * el[2] / np_image.shape[1]))
    if el[0] == "cos":
      shifts.append(lambda x: np_image.shape[0] / el[1] * \
        np.cos(x * el[2] / np_image.shape[1]))
    if el[0] == "triang":
      lambda x: np_image.shape[0] / el[1] * \
        (x / el[2] / np_image.shape[1] - math.floor(x / (el[2] / np_image.shape[1])))
    if el[0] == "erosion":
      np_image = erosion(np_image, square(el[1]))
    if el[0] == "dilation":
      np_image = dilation(np_image, square(el[1]))

  if orient == "horizont":
    for idx in xrange(np_image.shape[0]):
      for shift in shifts:
        np_image[idx,:] = np.roll(np_image[idx,:], int(shift(idx)))
  if orient == "vert":
    for idx in xrange(np_image.shape[1]):
      for shift in shifts:
        np_image[:, idx] = np.roll(np_image[:, idx], int(shift(idx)))

  return Image.fromarray(np_image)

def get_stomata(max_proj_image, min_obj_size=200, max_obj_size=1000):
    """Performs image segmentation from a max_proj_image.
     Disposes of objects in range min_obj_size to
    max_obj_size

    :param max_proj_image: the maximum projection image
    :type max_proj_image: numpy.ndarray, uint16
    :param min_obj_size: minimum size of object to keep
    :type min_obj_size: int
    :param max_obj_size: maximum size of object to keep
    :type max_obj_size: int
    :returns: list of [ [coordinates of kept objects - list of slice objects],
                        binary object image - numpy.ndarray,
                        labelled object image - numpy.ndarray
                     ]

    """

    # pore_margin = 10
    # max_obj_size = 1000
    # min_obj_size = 200
    # for prop, value in segment_options:
    #     if prop == 'pore_margin':
    #         pore_margin = value
    #     if prop == 'max_obj_size':
    #         max_obj_size = value
    #     if prop == 'min_obj_size':
    #         min_obj_size = value
    #
    # print(pore_margin)
    # print(max_obj_size)
    # print(min_obj_size)

    #rescale_min = 50
    #rescale_max= 100
    #rescaled = exposure.rescale_intensity(max_proj_image, in_range=(rescale_min,rescale_max))
    rescaled = max_proj_image
    seed = np.copy(rescaled)
    seed[1:-1, 1:-1] = rescaled.max()
    #mask = rescaled
    #if gamma != None:
    #    rescaled = exposure.adjust_gamma(max_proj_image, gamma)
    #filled = reconstruction(seed, mask, method='erosion')
    closed = dilation(rescaled)
    seed = np.copy(closed)
    seed[1:-1, 1:-1] = closed.max()
    mask = closed


    filled = reconstruction(seed, mask, method='erosion')
    label_objects, nb_labels = ndimage.label(filled)
    sizes = np.bincount(label_objects.ravel())
    mask_sizes = sizes
    mask_sizes = (sizes > min_obj_size) & (sizes < max_obj_size)
    #mask_sizes = (sizes > 200) & (sizes < 1000)
    mask_sizes[0] = 0
    big_objs = mask_sizes[label_objects]
    stomata, _ = ndimage.label(big_objs)
    obj_slices = ndimage.find_objects(stomata)
    return [obj_slices, big_objs, stomata]

def find_edges(img, sigma = 4):
    img = feature.canny(img, sigma)

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

    return img

def getMinorMajorRatio(image):
# First we threshold the image by only taking values greater than the mean to reduce noise in the image
# to use later as a mask

    image = image.copy()
    # Create the thresholded image to eliminate some of the background
    imagethr = np.where(image > np.mean(image),0.,1.0)

    #Dilate the image
    imdilated = morphology.dilation(imagethr, np.ones((4,4)))

    # Create the label list
    label_list = measure.label(imdilated)
    label_list = imagethr*label_list
    label_list = label_list.astype(int)
    # calculate common region properties for each region within the segmentation

    region_list = measure.regionprops(label_list)
    maxregion = getLargestRegion(region_list, label_list, imagethr)
    
    # guard against cases where the segmentation fails by providing zeros
    ratio = 0.0
    if ((not maxregion is None) and  (maxregion.major_axis_length != 0.0)):
        ratio = 0.0 if maxregion is None else  maxregion.minor_axis_length*1.0 / maxregion.major_axis_length
    return ratio

def segment_lowest_cluster(img_k):
    """Returns a binarized image with only the smallest 
    cluster of the kmeans."""

    mini_img_k = np.amin(img_k)
    darkest_cluster = dilation(img_k < mini_img_k + 1, disk(3))
    return darkest_cluster

def get_segmented_lungs(im):

    binary = im < -320
    cleared = clear_border(binary) 
    cleared=morph(cleared,5)
    label_image = label(cleared)
  
    areas = [r.area for r in regionprops(label_image)]
    areas.sort()
    if len(areas) > 2:
        for region in regionprops(label_image):
            if region.area < areas[-2]:
                for coordinates in region.coords:
                       label_image[coordinates[0], coordinates[1]] = 0
    binary = label_image > 0  
    selem = disk(2)
    binary = binary_erosion(binary, selem)
 
    selem = disk(10)
    binary = binary_closing(binary, selem)
    edges = roberts(binary)
    binary = ndi.binary_fill_holes(edges)
 
    get_high_vals = binary == 0
    im[get_high_vals] = 0
  
    binary = morphology.dilation(binary,np.ones([5,5]))
    return binary

def cluster_process(labels, original, activations):
    rbase = np.zeros(labels.shape)
    rubase = np.zeros(labels.shape)
    rubase[range(0,20),:] = 1
    rubase[:,range(0,20)] = 1
    rubase[range(-20,-1),:] = 1
    rubase[:,range(-20,-1)] = 1
    for i in range(1, int(np.max(labels))+1):
        base = np.zeros(labels.shape)
        base[labels==i] = 1
        li = len(base.nonzero()[0])
        if li>0:
            hull = convex_hull_image(base)
            lh =len(hull.nonzero()[0])
            sel_org = base*original
            sel_act = base*activations
            cond = (li > 4000 and float(lh) / float(li) < 1.07 and perimeter(base)**2.0 / li < 30) or np.max(base * rubase) > 0.5
            # print li>4000 and float(lh)/float(li)<1.07, perimeter(base)**2.0/li<30, np.max(base*rubase)>0.5, np.min(original[base>0])
            hard_array =[li > 4000, float(lh) / float(li) < 1.07]
            optional_array = [perimeter(base)**2.0/li < 25,
                              np.percentile(sel_org[sel_org>0], 5) > 0.2,
                              np.percentile(sel_act, 90) - np.percentile(sel_act, 90)]
            print hard_array, optional_array
            if debug and li>1000:
                rs(base,'subspread cluster')
            if cond:
                rbase = rbase + base
    rbase[rubase.astype(np.bool)] = 1
    return dilation(rbase, selem)

def load_scenes(filename):
    zipped_scenes = []
    print 'Working on: ' + filename
    img = data.imread('scenes/' + filename, as_grey=True)
    tmp = img
    tmp = filter.canny(tmp, sigma=2.0)
    tmp = ndimage.binary_fill_holes(tmp)
    #tmp = morphology.dilation(tmp, morphology.disk(2))
    tmp = morphology.remove_small_objects(tmp, 2000)
    contours = measure.find_contours(tmp, 0.8)
    ymin, xmin = contours[0].min(axis=0)
    ymax, xmax = contours[0].max(axis=0)
    if xmax - xmin > ymax - ymin:
        xdest = 1000
        ydest = 670
    else:
        xdest = 670
        ydest = 1000
    src = np.array(((0, 0), (0, ydest), (xdest, ydest), (xdest, 0)))
    dst = np.array(((xmin, ymin), (xmin, ymax), (xmax, ymax), (xmax, ymin)))
    tform3 = tf.ProjectiveTransform()
    tform3.estimate(src, dst)
    warped = tf.warp(img, tform3, output_shape=(ydest, xdest))
    tmp = filter.canny(warped, sigma=2.0)
    tmp = morphology.dilation(tmp, morphology.disk(2))
    descriptor_extractor.detect_and_extract(tmp)
    obj_key = descriptor_extractor.keypoints
    scen_desc = descriptor_extractor.descriptors
    zipped_scenes.append([warped, scen_desc, obj_key, filename])
    return zipped_scenes

def buffer_pcl(pcl):

    from skimage.morphology import dilation, disk
    selem = disk(3)
    dilated = dilation(pcl, selem)

    return dilated

def pestFeatureExtraction(filename):
	selem = disk(8)
	image = data.imread(filename,as_grey=True)
	thresh = threshold_otsu(image)
	elevation_map = sobel(image)
	markers = np.zeros_like(image)

	if ((image<thresh).sum() > (image>thresh).sum()):
		markers[image < thresh] = 1
		markers[image > thresh] = 2
	else:
		markers[image < thresh] = 2
		markers[image > thresh] = 1

	segmentation = morphology.watershed(elevation_map, markers)
	segmentation = dilation(segmentation-1, selem)
	segmentation = ndimage.binary_fill_holes(segmentation)

	segmentation = np.logical_not(segmentation)
	image[segmentation]=0;

	hist = np.histogram(image.ravel(),256,[0,1])

	hist = list(hist[0])
	hist[:] = [float(x) / (sum(hist) - hist[0]) for x in hist]
	hist.pop(0)

	features = np.empty( (1, len(hist)), 'float' )
	
	a = np.array(list(hist))
	f = a.astype('float')
	features[0,:]=f[:]

	return features

def blobs(image, remove_mb = None, val = 160, size = 100):
    """ Convolve a kernel on the image and a gaussian filter to highligh blobs. Find blobs using the
    Difference of Gaussian. Remove from the list of blobs the blobs that are at the membrane.
    return 3 different list
    """

    thresh = threshold_otsu(image)

    #Find all the blobs in the image using Difference of Gaussian
    blobs_in_image = feature.blob_dog(image, min_sigma=0.01,
                        max_sigma=3, threshold=thresh)
    blob_list = []
    for blob in blobs_in_image:
        y, x, r = blob
        blob_list.append((y, x))



    if remove_mb == None:
        blob_in_image_after_binary = set(blob_list)

    else:
        #Create a mask to remove blobs that are at the membrane and surrounded
        #by bright big object
        binary = image >= val*thresh/100
        binary = dilation(binary, square(3))
        binary = remove_small_objects(binary, min_size=size)
        # Create a list of coordinate with the binary image
        coor_binary = np.nonzero(binary)
        list_blob_masked = zip(*coor_binary)
        #Substract the list of coordinate from the binary image to the list of blobs
        blob_in_image_after_binary = (set(blob_list) - set (list_blob_masked))

    return blob_in_image_after_binary

 def morpho_rec2(self, img, size=10):
     # internal gradient of the cells: 
     se = morphology.diamond(size)
     dil = morphology.dilation(img, se)
     rec = morphology.reconstruction(dil, img, method='erosion').astype(np.dtype('uint8'))
             
     return rec

def make_thicker(img, filter_size):
    # Make greyscale image thicker
    # Parameter
    
    img = dilation(img, disk(filter_size))

    return img

def mask_to_objects_2d(mask, background=0, offset=None):
    """Convert 2D (binary or label) mask to polygons. Generates borders fitting in the objects.
    Parameters
    ----------
    mask: ndarray
        2D mask array. Expected shape: (height, width).
    background: int
        Value used for encoding background pixels.
    offset: tuple (optional, default: None)
        (x, y) coordinate offset to apply to all the extracted polygons.
    Returns
    -------
    extracted: list of AnnotationSlice
        Each object slice represent an object from the image. Fields time and depth of AnnotationSlice are set to None.
    """
    if mask.ndim != 2:
        raise ValueError("Cannot handle image with ndim different from 2 ({} dim. given).".format(mask.ndim))
    if offset is None:
        offset = (0, 0)
    # opencv only supports contour extraction for binary masks: clean mask and binarize
    mask_cpy = np.zeros(mask.shape, dtype=np.uint8)
    mask_cpy[mask != background] = 255
    # create artificial separation between adjacent touching each other + clean
    contours = dilation(mask, square(3)) - mask
    mask_cpy[np.logical_and(contours > 0, mask > 0)] = background
    mask_cpy = clean_mask(mask_cpy, background=background)
    # extract polygons and labels
    polygons = _locate(mask_cpy, offset=offset)
    objects = list()
    for polygon in polygons:
        # loop for handling multipart geometries
        for curr in flatten_geoms(polygon.geoms) if hasattr(polygon, "geoms") else [polygon]:
            x, y = get_polygon_inner_point(curr)
            objects.append((polygon, mask[y - offset[1], x - offset[0]]))
    return objects