def ApplyThresholdToImageRegion(image2, Tb, Bb, Lb, Rb,shouldThresholdImage):

    image = rgb2gray(image2)

    global foregroundPixelValue
    global backgroundPixelValue

    thresholdValue = threshold_otsu(image)

    NumberOfRows = image.shape[0]
    NumberOfColumns = image.shape[1]

    numberOfBlackPixels = 0
    numberOfWhitePixels = 0
    selem = disk(3)

    # simpe thresholding
    for y in range(NumberOfRows):
        for x in range(NumberOfColumns):

            isWithinBoundary = IsWithinBoundary(y,x,image2, Tb, Bb, Lb, Rb,shouldThresholdImage)

            if (isWithinBoundary):
                if image[y,x] > thresholdValue:
                    #black
                    image[y,x] = 0
                    numberOfBlackPixels += 1
                else:
                    #white
                    image[y,x] = 1
                    numberOfWhitePixels += 1

    # assume foreground has more pixels in face region
    if (numberOfWhitePixels > numberOfBlackPixels):
        foregroundPixelValue = 1
        backgroundPixelValue = 0
        #print("foreground color is white")
    else:
        foregroundPixelValue = 0
        backgroundPixelValue = 1
        #print("foreground color is black")

    image = opening(image,selem)
    if (foregroundPixelValue == 0):
        image = opening(image,selem)
    else:
        image = closing(image,selem)


    if drawFaceImages:
        io.imshow(image)
        io.show()

    return image

def upsample_smooth(image):
  filename_split = os.path.splitext(image)
  filename_zero, fileext = filename_split
  basename = os.path.basename(filename_zero)
  upsampleRes = int(upsampleRes)
  upsampled_image = outPath+basename+'_cubicSpline.png'
  os.system("gdalwarp -tr", upsample_res, upsample_res," -r cubicspline ", image, upsampled_image)
  im = np.array(Image.open(upsampled_image))
  with rasterio.open(image) as r:
    im = r.read()
    p = r.profile
  im = im.squeeze()
  selem = disk(1)
  print("image min and max: ", im.min(),im.max())
  dilated = skimage.morphology.dilation(im, selem)
  print("dilated image min and max: ", dilated.min(),dilated.max())
  eroded = skimage.morphology.erosion(im, selem)
  print("eroded image min and max: ", eroded.min(),eroded.max())
  opened = opening(im, selem)
  print("opened image min and max: ", opened.min(),opened.max())
  closed = closing(im, selem)
  print("closed image min and max: ", closed.min(),closed.max())
  dilated = Image.fromarray(dilated)
  dilated.save(outPath+basename+'.png')
  with rasterio.open(outPath+basename+fileext, 'w', **p) as dst:
      dst.write(dilated, 1)
  color_outPath = outPath+'color/'
  if not os.path.exists(color_outPath):
        os.mkdir(color_outPath) 
  colored_image = color_outPath+basename+'.png'
  os.system("gdaldem color-relief", dilated, colorfile, colored_image)
  return im, dilated, eroded, opened, closed

def Mask_ROI_op(im,disk_size,thresh=None,black_spots=None,with_morph=False):
    l=np.array([0,0])
    if not isinstance(im,l.__class__):
        numpy_array=np.array(im)
    else:
        numpy_array=im 
    if len(numpy_array.shape)==3:
        numpy_array=numpy_array[:,:,0:3].mean(axis=2)        
    selem = disk(disk_size)
    openin = opening(numpy_array, selem)
    if thresh is None:
        thresh = threshold_otsu(openin)
    binary = openin > thresh
    if binary.dtype=='bool':
        binary=binary+0
    if black_spots is not None:    
        binary2 = openin > black_spots
        binary2 = binary2 + 0
        binary = binary - binary2 
    else:
        binary -=1
    binary=binary * -1
    if with_morph:
        return(binary,openin)
    else:
        return(binary)                     

def open_image(img, mask_length):
    # Morphological opening on greyscale/binary image
    
    img = img.astype(np.uint8)
    img = opening(img, rectangle(mask_length,1))
    
    return(img)

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 smooth(image):
  filename_split = os.path.splitext(image)
  filename_zero, fileext = filename_split
  basename = os.path.basename(filename_zero)
  im = np.array(Image.open(image))
  with rasterio.open(image) as r:
    im = r.read()
    p = r.profile
  im = im.squeeze()
  selem = disk(1)
  print("image min and max: ", im.min(),im.max())
  dilated = skimage.morphology.dilation(im, selem)
  print("dilated image min and max: ", dilated.min(),dilated.max())
  eroded = skimage.morphology.erosion(im, selem)
  print("eroded image min and max: ", eroded.min(),eroded.max())
  opened = opening(im, selem)
  print("opened image min and max: ", opened.min(),opened.max())
  closed = closing(im, selem)
  print("closed image min and max: ", closed.min(),closed.max())
  #im[im==1]=0
  #im[im==2]=1
  median = cv2.medianBlur(im,9)
  average = cv2.blur(im,(9,9))
  #gaussian = cv2.GaussianBlur(im,(9,9),0)
  gaussian = cv2.GaussianBlur(dilated,(9,9),0)
  #bilateral = cv2.bilateralFilter(im,9,75,75)
  bilateral = cv2.bilateralFilter(gaussian,9,75,75)
  with rasterio.open(outPath+basename+fileext, 'w', **p) as dst:
      dst.write(bilateral, 1)
  color_outPath = outPath+'color/'
  if not os.path.exists(color_outPath):
        os.mkdir(color_outPath) 
  colored_image = color_outPath+basename+'.png'
  os.system("gdaldem color-relief", bilateral, colorfile, colored_image)
  return im, dilated, eroded, opened, closed, median, average, gaussian, bilateral

def opening(gray_img, kernel=None):
    """Wrapper for scikit-image opening functions. Opening can remove small bright spots (i.e. salt).

    Inputs:
    gray_img = input image (grayscale or binary)
    kernel   = optional neighborhood, expressed as an array of 1s and 0s. If None, use cross-shaped structuring element.

    :param gray_img: ndarray
    :param kernel = ndarray
    :return filtered_img: ndarray
    """

    params.device += 1

    # Make sure the image is binary/grayscale
    if len(np.shape(gray_img)) != 2:
        fatal_error("Input image must be grayscale or binary")

    # If image is binary use the faster method
    if len(np.unique(gray_img)) == 2:
        bool_img = morphology.binary_opening(gray_img, kernel)
        filtered_img = np.copy(bool_img.astype(np.uint8) * 255)
    # Otherwise use method appropriate for grayscale images
    else:
        filtered_img = morphology.opening(gray_img, kernel)

    if params.debug == 'print':
        print_image(filtered_img, os.path.join(params.debug_outdir, str(params.device) + '_opening' + '.png'))
    elif params.debug == 'plot':
        plot_image(filtered_img, cmap='gray')

    return filtered_img

def get_nuclei(img, opening_radius=6, block_size=80, threshold_offset=0):
    s = Sample(DOWNSAMPLE)
    binary = threshold_adaptive(s.downsample(img), int(block_size / s.rate), offset=threshold_offset)
    filled = fill_holes(binary)
    opened = opening(filled, selem=disk(opening_radius / s.rate))
    nuclei = apply_watershed(opened)
    nuclei = s.upsample(nuclei)
    return img_as_uint(nuclei)

def segment_image(image, plantid, threshold, locations, opening_size=5):
    """
        Segments an image based on simple thresholding
        
        Inputs:
        - image : the image (as a numpy array)
        - plantid : which plant on the image (int)
        - threshold : the threshold to use
        - locations : the coordinates of the plants in the
        
        Output:
        - mask of the image
        """
    image_part = get_piece(image, locations[plantid])
    mask = median_filter(image_part.mean(2), (4,4)) > threshold
    opening(mask, selem=square(opening_size), out=mask)  # open image, remove
    # clutter
    return mask  # gaussian_filter(image_part[mask, :], sigma=convolve_sigma)

 def func(frame):
     _dtype = frame.dtype
     kernel = mor.disk(3)
     frameWP = frame - mor.white_tophat(frame, kernel) * (mor.white_tophat(frame, kernel) > 1000).astype(float)
     kernel = mor.rectangle(25, 1)
     closed = mor.closing(frameWP, kernel)
     opened = mor.opening(closed, kernel)
     result = ((frameWP.astype(float) / opened.astype(float)) * 3000.0)
     return result.astype(_dtype)

def ProcessImage(im, targetDim = 250, doDenoiseOpening = True):

	#Resize to specified pixels max edge size
	scaling = 1.
	if im.shape[0] > im.shape[1]:
		if im.shape[0] != targetDim:
			scaling = float(targetDim) / im.shape[0]
			im = misc.imresize(im, (targetDim, int(round(im.shape[1] * scaling))))
	else:
		if im.shape[1] != targetDim:
			scaling = float(targetDim) / im.shape[1]
			im = misc.imresize(im, (int(round(im.shape[0] * scaling)), targetDim))
	#print "scaling", scaling

	greyim = 0.2126 * im[:,:,0] + 0.7152 * im[:,:,1] + 0.0722 * im[:,:,2]

	#Highlight number plate
	imnorm = np.array(greyim, dtype=np.uint8)
	se = np.ones((3, 30), dtype=np.uint8)
	opim = morph.opening(imnorm, se)
	diff = greyim - opim + 128.

	misc.imsave("diff.png", diff)

	#Binarize image
	vals = diff.copy()
	vals = vals.reshape((vals.size))

	meanVal = vals.mean()
	stdVal = vals.std()
	threshold = meanVal + stdVal

	#print "Threshold", threshold

	binIm = diff > threshold
	misc.imsave("threshold.png", binIm)
	#print vals.shape
	#plt.plot(vals)
	#plt.show()

	#Denoise
	diamond = morph.diamond(2)
	if doDenoiseOpening:
		currentIm = morph.binary_opening(binIm, diamond)
	else:
		currentIm = binIm
	denoiseIm2 = morph.binary_closing(currentIm, np.ones((3, 13)))

	#print "currentIm", currentIm.min(), currentIm.max(), currentIm.mean()
	#print "denoiseIm2", denoiseIm2.min(), denoiseIm2.max(), currentIm.mean()
	#misc.imsave("denoised1.png", currentIm * 255)
	#misc.imsave("denoised2.png", denoiseIm2 * 255)

	#Number candidate regions
	#print "Numbering regions"
	numberedRegions, maxRegionNum = morph.label(denoiseIm2, 4, 0, return_num = True)
	return numberedRegions, scaling

def find_border_centroids(hfp, keys, areas, largeobjkey, disttransfkey, resultkey):

    for k, bounds in keys.iteritems():

        # bounds = (shp[0],) * 2
        for lbl, lblim in hfp['faces', largeobjkey].label_image_iterator(key=k, background=0, area=areas[k]):

            hfp.logging('---\nLabel {} found in image {}', lbl, k)

            # Avoid very small artifacts
            lblim = morphology.opening(lblim)

            # Connected component analysis to detect when a label touches the border multiple times
            conncomp = vigra.analysis.labelImageWithBackground(lblim.astype(np.uint32), neighborhood=8, background_value=0)

            for l in np.unique(conncomp):
                # Ignore background
                if l == 0: continue

                # Get the current label object
                curobj = conncomp == l

                # Get disttancetransf of the object
                curdist = np.array(hfp['faces', disttransfkey, k])
                curdist[curobj == False] = 0

                # Detect the global maximum of this object
                amax = np.amax(curdist)
                curdist[curdist < amax] = 0
                curdist[curdist > 0] = lbl
                # Only one pixel is allowed to be selected
                bds = lib.find_bounding_rect(curdist)
                centroid = (int((bds[1][0] + bds[1][1]-1) / 2), int((bds[0][0] + bds[0][1]-1) / 2))

                # Now translate the calculated centroid to the position within the orignial 3D volume
                centroidm = (centroid[0] - bounds, centroid[1] - bounds)
                # hfp.logging('centroidxy = {}', centroidm)
                # Set the pixel
                try:
                    if centroidm[0] < 0 or centroidm[1] < 0:
                        raise IndexError
                    else:
                        if k == 'xyf':
                            hfp[resultkey][centroidm[0], centroidm[1], 0] = lbl
                        elif k == 'xyb':
                            hfp[resultkey][centroidm[0], centroidm[1], -1] = lbl
                        elif k == 'xzf':
                            hfp[resultkey][centroidm[0], 0, centroidm[1]] = lbl
                        elif k == 'xzb':
                            hfp[resultkey][centroidm[0], -1, centroidm[1]] = lbl
                        elif k == 'yzf':
                            hfp[resultkey][0, centroidm[0], centroidm[1]] = lbl
                        elif k == 'yzb':
                            hfp[resultkey][-1, centroidm[0], centroidm[1]] = lbl
                except IndexError:
                    pass

def GrayMM(img,thresh,size):
    structElement = morphology.rectangle(size,size)
  
    img[img == -9999] = 0
    img = img/5000
    img[img < 0] = 0
    img[img > 1] = 1
    
    outdata = morphology.opening(img, structElement)
    #threshold after bin
    imgOut = outdata > 255*thresh/5000
    
    return imgOut

def ROI_binary_mask(sample, size=5, ticket=(80, 80, 80)):
    #  very slow function at resolution 4
    PreprocRes = np.copy(sample)

    for i in range(3):  # RGB
        # this one is painfully slow..
        PreprocRes[:, :, i] = Preprocessing(sample[:, :, i])
    res = combining(PreprocRes)
    ticket = FindTicket(sample, ticket)
    res = res - ticket
    res[res > 0] = 1
    res[res < 0] = 0
    res = opening(res, disk(size))
    return res

def buffer_mask(stuff):
    stuff_shape = stuff.shape
    stuff_buff = np.array(np.copy(stuff))
    for x in range(1,stuff_shape[0]):
        for y in range(1,stuff_shape[1]):
            if stuff[x,y] == False:
                if np.count_nonzero(stuff[x-1:x+2,y-1:y+2]) > 4:
                    # print x,y
                    stuff_buff[x,y] = True
    selem = morphology.disk(1)
    stuff_buff = morphology.opening(stuff_buff, selem)
    selem = morphology.disk(1)
    stuff_buff = morphology.closing(stuff_buff, selem)
    return stuff_buff

def label_clouds(pcp, opening_selem, closing_selem):
    # img = buffer_pcl()
    # pcp = calc_pcp
    img = pcp
    _img = np.zeros((img.shape[0]+30, img.shape[1]+30))
    _img[15:-15, 15:-15] = img
    
    o_selem   = morphology.disk(opening_selem)
    _img = morphology.opening(_img, o_selem)
    c_selem   = morphology.disk(closing_selem)
    bin_im = morphology.closing(_img, c_selem)
    
    labels, nbr_objs = measurements.label(bin_im)
    print(nbr_objs)
    return labels[15:-15, 15:-15], nbr_objs

	def clean_footprint(self, OPENING_THRESH = 5):
		self.FOOTPRINT_cleaned_opening = opening(self.FOOTPRINT_added_boundary, square(OPENING_THRESH))

		labeled_image, num_of_labels = ndimage.measurements.label(self.FOOTPRINT_cleaned_opening, [[1,1,1],[1,1,1],[1,1,1]])
		
		if num_of_labels > 1:
			b_slices = ndimage.find_objects(labeled_image)
			area = []
			for idx,s in enumerate(b_slices):
				area.append(len(self.FOOTPRINT_cleaned_opening[self.FOOTPRINT_cleaned_opening[s] != 0]))
			self.FOOTPRINT_cleaned_opening[labeled_image != (area.index( max(area))) + 1] = 0


		# area = np.bincount(obj.flatten())[indexedObject[0]]
		return len(self.FOOTPRINT_cleaned_opening[self.FOOTPRINT_cleaned_opening!=0])

def FindTicket(RGB_image, _3tuple=(80, 80, 80)):
    # Find the "black ticket on the images"
    temp_image_3 = np.copy(RGB_image)
    temp_image_3[:, :, :] = 0
    for i in range(3):
        temp_image_1 = np.zeros(shape=RGB_image.shape[0:2])
        temp_image_1[np.where(RGB_image[:, :, i] < _3tuple[i])] = 1
        temp_image_3[:, :, i] = temp_image_1

    temp_resultat = temp_image_3.sum(axis=2)

    temp_resultat[temp_resultat > 2] = 3
    temp_resultat[temp_resultat < 3] = 0
    temp_resultat[temp_resultat == 3] = 1
    
    #temp_resultat = Filling_holes_2(temp_resultat)
    temp_resultat = closing(temp_resultat, disk(20))
    temp_resultat = opening(temp_resultat, disk(20))
    temp_resultat = RemoveBorder(temp_resultat)
    return temp_resultat

    def _floodfill(self, img):
        back = self._scharr(img)
        # Binary thresholding.
        back = back > 0.05

        # Thin all edges to be 1-pixel wide.

        back = skm.skeletonize(back)
        # Edges are not detected on the borders, make artificial ones.
        back[0, :] = back[-1, :] = True
        back[:, 0] = back[:, -1] = True

        # Label adjacent pixels of the same color.
        labels = label(back, background=-1, connectivity=1)

        # Count as background all pixels labeled like one of the corners.
        corners = [(1, 1), (-2, 1), (1, -2), (-2, -2)]
        for l in (labels[i, j] for i, j in corners):
            back[labels == l] = True

        # Remove remaining inner edges.
        return skm.opening(back)

def change_bg(orig, background=None,
              threshold_yellow=[(80, 145), (65, 145), (0, 50)],
              threshold_black=[(0, 30), (0, 30), (0, 30)],
              threshold_face=[(42, 90), (20, 70), (0, 42)],
              show_pic=False,
              ):
    thresholds = [threshold_yellow, threshold_black, threshold_face]
    if background is None:
        background = orig * 0
    background = array(Image.fromarray(background).resize(orig.shape[:2][::-1]))
    masks = [orig.copy() for i in range(len(thresholds))]
    for mask_i, threshold in enumerate(thresholds):
        for i, (th_lo, th_hi) in enumerate(threshold):
            masks[mask_i][:, :, i] = (th_lo <= masks[mask_i][:, :, i]).astype(int) * (masks[mask_i][:, :, i] <= th_hi).astype(int)
        mask = masks[mask_i][:, :, 0] * masks[mask_i][:, :, 1] * masks[mask_i][:, :, 2]
        masks[mask_i][:, :, 0], masks[mask_i][:, :, 1], masks[mask_i][:, :, 2] = [mask for i in range(3)]
    mask = 1-reduce(lambda x, y: x * y, [1-w for w in masks])
    mask_orig = mask.copy()
    mask = morphology.dilation(mask, ones([10, 3, 1]))
    mask = morphology.opening(mask, ones([18, 12, 1]))
    mask = morphology.dilation(mask, ones([3, 6, 1]))

    result = (orig * mask + background * (1 - mask))
    if show_pic:
        subplot(221)
        imshow(orig)
        title('orig')
        subplot(222)
        imshow(mask*255)
        title('mask')
        subplot(223)
        imshow(mask_orig*255)
        title('mask_orig')
        subplot(224)
        imshow(result)
        title('result')
        show()
    return result, mask, mask_orig