def calculate_masked_stats():
    plate_no = "59798"
    parsed = get_plate_files(plate_no)
    for w in ['w2']:
        files = filter(lambda f: f.wave == w[1], parsed)
        # accum = np.zeros((2160, 2160), dtype=np.uint32)
        # files = filter(lambda x: 's1' not in x and 's7' not in x, all_files)
        nof = len(files)
        for i, frame in enumerate(files[0:5], 1):
            LogHelper.logText(frame.fullpath)
            img = imread(frame.fullpath)
            t = filters.threshold_yen(img)
            b1 = img > t
            b2 = binary_erosion(b1, square(2))
            b3 = binary_dilation(b2, square(10))
            b4 = binary_closing(b3, square(3))
            imm = np.ma.masked_where(b4, img)
            mn, mx = np.percentile(imm, (1, 99))
            LogHelper.logText(
                '%3d of %d, %4d-%4d-%4d-%5d, %.0f-%.0f'
                % (i, nof, imm.min(), mn, mx, imm.max(), imm.mean(), imm.std())
            )
            im2 = imm.filled(int(imm.mean()))
            out_name = "{0}\\{5}-{1}{2}-{3}-{4}.tif".format(ROOT_DIR, frame.row, frame.column, frame.site, LogHelper.init_ts, frame.experiment)
            imsave(out_name, im2)

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 morphoNoiseRemoval(img):
    "Removes noise by succession of 5 opening/closing morphological operators"
    for i in range(0,5):
        img = opening2(img, square(3))
        img = closing2(img, square(3))
        
    return img

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 removeChessboard(img):

    # Get the major lines in the image
    edges, dilatedEdges, (h, theta, d) = findLines(img)

    # Create image with ones to fill inn lines
    lines = np.ones(img.shape[:2])

    # Add lines to image as zeroes
    for _, angle, dist in zip(*hough_line_peaks(h, theta, d)):
        y0 = (dist - 0 * np.cos(angle)) / np.sin(angle)
        y1 = (dist - img.shape[1] * np.cos(angle)) / np.sin(angle)
        x, y = line(int(y1), 0, int(y0), img.shape[1] - 1)
        x = np.clip(x, 0, img.shape[0] - 1)
        y = np.clip(y, 0, img.shape[1] - 1)
        lines[x, y] = 0

    # Remove border edges from image with all edges
    w = 4
    edges = np.pad(edges[w:img.shape[0] - w, w:img.shape[1] - w], w, mode='constant')

    # Erode the lines bigger, such that they cover the original lines
    lines = erosion(lines, square(13))

    # Remove major lines and close shape paths
    removedChessboard = closing(edges * lines, square(8))

    return removedChessboard

def morph(img, tparams):
    ops = [mor.grey.erosion, mor.grey.dilation]
    t = ops[np.random.randint(2)] 
    if t == 0:    
        selem = mor.square(np.random.randint(1, tparams['selem_size'][0]))
    else:
        selem = mor.square(np.random.randint(1, tparams['selem_size'][1]))
    return t(img, selem)    

 def seg_sect(self, img):
     img_canny = canny(img, sigma=self.sigma,
                       low_threshold=self.low_threshold)
     
     img_dilate = binary_dilation(img_canny, square(3))
     img_erode = binary_erosion(img_dilate, square(3))
     img_fill = binary_fill_holes(img_erode)
     
     return img_fill

def process_cell(img):

    # la binariza en caso de que sea escala de grises
    if not img.dtype == 'bool':
        img = img > 0  # Binarizar

    # Calcular máscaras para limpiar lineas largas verticales
    h_k = 0.8
    sum0 = np.sum(img, 0)  # Aplastar la matriz a una fila con las sumas de los valores de cada columna.
    thr0 = sum0 < h_k * img.shape[0]
    thr0 = thr0.reshape(len(thr0), 1) # Convertirlo a vector de una dimensión

    # Calcular máscaras para limpiar lineas largas horizontales
    w_k = 0.5
    sum1 = np.sum(img, 1)
    thr1 = sum1 < w_k * img.shape[1]
    thr1 = thr1.reshape(len(thr1), 1)

    mask = thr0.transpose() * thr1 # Generar máscara final para la celda
    mask_lines = mask.copy()

    elem = morphology.square(5)
    mask = morphology.binary_erosion(mask, elem) # Eliminar ruido

    img1 = np.bitwise_and(mask, img) # Imagen filtrada

    # segmentación del bloque de números
    kerw = 5  # Kernel width
    thr_k = 0.8

    # Calcular mascara para marcar inicio y fin de región con dígitos horizontalmente
    sum0 = np.sum(img1, 0)
    sum0 = signal.medfilt(sum0, kerw)
    thr0 = sum0 > thr_k * np.median(sum0)
    thr0 = np.bitwise_and(thr0.cumsum() > 0, np.flipud(np.flipud(thr0).cumsum() > 0))
    thr0 = thr0.reshape(len(thr0), 1)

    # Calcular mascara para marcar inicio y fin de región con dígitos verticalmente
    sum1 = np.sum(img1, 1)
    sum1 = signal.medfilt(sum1, kerw)
    thr1 = sum1 > thr_k * np.median(sum1)
    thr1 = np.bitwise_and(thr1.cumsum() > 0, np.flipud(np.flipud(thr1).cumsum() > 0))
    thr1 = thr1.reshape(len(thr1), 1)

    # Mascara final para inicio y fin de caracteres (bounding box of digit region)
    mask = thr0.transpose() * thr1
    mask = morphology.binary_dilation(mask, morphology.square(2))


    img = np.bitwise_and(mask_lines.astype(img.dtype), img)  # Aplicar máscara para quitar lineas
    img = morphology.binary_dilation(img, morphology.disk(1)) # Dilatación para unir números quebrados por la máscara anterior
    img = morphology.binary_erosion(img, morphology.disk(1)) # Volver a la fomorma 'original' con los bordes unidos

    return np.bitwise_and(mask, img)

def _getPoseMask(peaks, height, width, radius=4, var=4, mode='Solid'):
    ## MSCOCO Pose part_str = [nose, neck, Rsho, Relb, Rwri, Lsho, Lelb, Lwri, Rhip, Rkne, Rank, Lhip, Lkne, Lank, Leye, Reye, Lear, Rear, pt19]
    # find connection in the specified sequence, center 29 is in the position 15
    # limbSeq = [[2,3], [2,6], [3,4], [4,5], [6,7], [7,8], [2,9], [9,10], \
    #            [10,11], [2,12], [12,13], [13,14], [2,1], [1,15], [15,17], \
    #            [1,16], [16,18], [3,17], [6,18]]
    # limbSeq = [[2,3], [2,6], [3,4], [4,5], [6,7], [7,8], [2,9], [9,10], \
    #            [10,11], [2,12], [12,13], [13,14], [2,1], [1,15], [15,17], \
    #            [1,16], [16,18]] # , [9,12]
    # limbSeq = [[3,4], [4,5], [6,7], [7,8], [9,10], \
    #            [10,11], [12,13], [13,14], [2,1], [1,15], [15,17], \
    #            [1,16], [16,18]] # 
    limbSeq = [[2,3], [2,6], [3,4], [4,5], [6,7], [7,8], [2,9], [9,10], \
                         [10,11], [2,12], [12,13], [13,14], [2,1], [1,15], [15,17], \
                         [1,16], [16,18], [2,17], [2,18], [9,12], [12,6], [9,3], [17,18]] #
    indices = []
    values = []
    for limb in limbSeq:
        p0 = peaks[limb[0] -1]
        p1 = peaks[limb[1] -1]
        if 0!=len(p0) and 0!=len(p1):
            r0 = p0[0][1]
            c0 = p0[0][0]
            r1 = p1[0][1]
            c1 = p1[0][0]
            ind, val = _getSparseKeypoint(r0, c0, 0, height, width, radius, var, mode)
            indices.extend(ind)
            values.extend(val)
            ind, val = _getSparseKeypoint(r1, c1, 0, height, width, radius, var, mode)
            indices.extend(ind)
            values.extend(val)
        
            distance = np.sqrt((r0-r1)**2 + (c0-c1)**2)
            sampleN = int(distance/radius)
            # sampleN = 0
            if sampleN>1:
                for i in xrange(1,sampleN):
                    r = r0 + (r1-r0)*i/sampleN
                    c = c0 + (c1-c0)*i/sampleN
                    ind, val = _getSparseKeypoint(r, c, 0, height, width, radius, var, mode)
                    indices.extend(ind)
                    values.extend(val)

    shape = [height, width, 1]
    ## Fill body
    dense = np.squeeze(_sparse2dense(indices, values, shape))
    ## TODO
    # im = Image.fromarray((dense*255).astype(np.uint8))
    # im.save('xxxxx.png')
    # pdb.set_trace()
    dense = dilation(dense, square(5))
    dense = erosion(dense, square(5))
    return dense

    def plot_preprocessed_image(self):
        """
        plots pre-processed image. The plotted image is the same as obtained at the end
        of the get_text_candidates method.
        """
        image = restoration.denoise_tv_chambolle(self.image, weight=0.1)
        thresh = threshold_otsu(image)
        bw = closing(image > thresh, square(2))
        cleared = bw.copy()

        label_image = measure.label(cleared)
        borders = np.logical_xor(bw, cleared)

        label_image[borders] = -1
        image_label_overlay = label2rgb(label_image, image=image)

        fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(12, 12))
        ax.imshow(image_label_overlay)

        for region in regionprops(label_image):
            if region.area < 10:
                continue

            minr, minc, maxr, maxc = region.bbox
            rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr,
                                      fill=False, edgecolor='red', linewidth=2)
            ax.add_patch(rect)

        plt.show()

def median_filter(image, selem=None):
    if selem is None:
        # default mask is 5x5 square
        selem = square(5)
    depth = image.shape[2]
    return np.dstack(median(channel[...,0], selem)
                     for channel in np.dsplit(image, depth)) / 255.

 def get_rough_detection(self, img, bigsize=40.0, smallsize=4.0, thresh = 0):
     diff = self.difference_of_gaussian(-img, bigsize, smallsize)
     diff[diff>thresh] = 1
     
     se = morphology.square(4)
     ero = morphology.erosion(diff, se)
     
     labimage = label(ero)
     #rec = morphology.reconstruction(ero, img, method='dilation').astype(np.dtype('uint8'))
     
     # connectivity=1 corresponds to 4-connectivity.
     morphology.remove_small_objects(labimage, min_size=600, connectivity=1, in_place=True)
     #res = np.zeros(img.shape)
     ero[labimage==0] = 0
     ero = 1 - ero
     labimage = label(ero)
     morphology.remove_small_objects(labimage, min_size=400, connectivity=1, in_place=True)
     ero[labimage==0] = 0
     res = 1 - ero
     res[res>0] = 255
     
     #temp = 255 - temp
     #temp = morphology.remove_small_objects(temp, min_size=400, connectivity=1, in_place=True)
     #res = 255 - temp
     
     return res

def squareMask(maskImg, square_width): #both odd and even square_with are allowed
    boxsize = maskImg.get_xsize()
    maskArray = EMNumPy.em2numpy(maskImg)
    
    if (boxsize <= square_width):
        print "ERROR: the width of the square cannot be larger than the boxsize of particles."
        sys.exit()
        
    #from skimage.morphology import square
    #Generates a flat, square-shaped structuring element.
    #Every pixel along the perimeter has a chessboard distance no greater than radius (radius=floor(width/2)) pixels.
    squareArray = square(square_width, dtype=np.uint8)
    m, n = squareArray.shape
    assert m==n
    
    if (m%2 == 0):
        pad_before = (boxsize - m)/2
        pad_after = (boxsize - m)/2
    else:
        pad_before = (boxsize - m)/2
        pad_after = (boxsize - m)/2+1
    #pad_width = (boxsize - square_width)/2
    #print "m, n, pad_before, pad_after", m, n, pad_before, pad_after
    #squareArrayPad = np.pad(squareArray, pad_width, mode='constant')
    squareArrayPad = np.pad(squareArray, (pad_before, pad_after), mode='constant')
    
    squareImg = EMNumPy.numpy2em(squareArrayPad)
    return squareImg

def process_image(image):
    tic = time.clock()
    # rescale intensity
    p2, p98 = np.percentile(image, (1, 99.9))
    image = rescale_intensity(1.0*image, in_range=(p2, p98))

    # do simple filter based on color value
    thresh = 0.5*threshold_func(image)
    filtered_image = np.zeros_like(image,dtype=np.uint8) # set up all-zero image
    filtered_image[image > thresh] = 1 # filtered values set to 1

    # perform watershed transform to split clusters
    distance = ndi.distance_transform_edt(filtered_image)
    local_maxi = peak_local_max(distance, indices=False, footprint=morphology.square(7),
                            labels=filtered_image, exclude_border=False)
    markers = ndi.label(local_maxi)[0]

    # segment and label particles
    labels = morphology.watershed(-distance, markers, mask=filtered_image)
    backup_labels = labels.copy()

    # remove boundaries and restore any small particles deleted in this process
    labels[find_boundaries(labels)] = 0
    for i in np.unique(backup_labels)[1:]:
        if np.count_nonzero(labels[backup_labels == i]) == 0:
            labels[backup_labels == i] = i
    toc = time.clock()
    procTime = toc - tic
    return image, labels, procTime

    def run3(self):
        """ Cette fonction test des alternatives à SIFT et ORB. Ne fonctionne pas."""
        for x in xrange(len(self.stack)-1):
            print('Traitement image ' + str(x+1))
            im1,im2 = 255.*gaussian_filter(self.stack[x,...], sqrt(self.initial_sigma**2 - 0.25)), 255.*gaussian_filter(self.stack[x+1,...], sqrt(self.initial_sigma**2 - 0.25))
            im1,im2 = enhance_contrast(normaliser(im1), square(3)), enhance_contrast(normaliser(im2), square(3))
            im1, im2 = normaliser(im1), normaliser(im2)
            
            b = cv2.BRISK()
            #b.create("Feature2D.BRISK")
            
            k1,d1 = b.detectAndCompute(im1,None)
            k2,d2 = b.detectAndCompute(im2,None)
            
            bf = cv2.BFMatcher(cv2.NORM_HAMMING)
            matches = bf.match(d1,d2)
            
            g1,g2 = [],[]
            for i in matches:
                g1.append(k1[i.queryIdx].pt)
                g2.append(k2[i.trainIdx].pt)

            model, inliers = ransac((np.array(g1), np.array(g2)), AffineTransform, min_samples=3, residual_threshold=self.min_epsilon, max_trials=self.max_trials, stop_residuals_sum=self.min_inlier_ratio)
            
            self.stack[x+1,...] = warp(self.stack[x+1,...], AffineTransform(rotation=model.rotation, translation=model.translation), output_shape=self.stack[x+1].shape)

        self.stack = self.stack.astype(np.uint8)

def segmentation(image):
	"""Executes image segmentation based on various features of the video stream"""
	gray = cv2.cvtColor(image, cv2.cv.CV_RGB2GRAY)
	blurred = cv2.GaussianBlur(gray, (5, 5), 0)

	ret, bw = cv2.threshold(blurred, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
	
	# Close binary image
	result = cv2.dilate(bw, square(10), iterations = 1)
	result = cv2.erode(result, square(10), iterations = 1)
	
	# Open binary image
	result = cv2.erode(result, square(10), iterations = 1)
	result = cv2.dilate(result, square(10), iterations = 1)
	
	return label(result) * 50

    def updateRasterInfo(self, **kwargs):
        kwargs['output_info']['statistics'] = ()
        kwargs['output_info']['histogram'] = ()

        self.window = square(int(kwargs.get('size', 3)))
        m = kwargs.get('measure', 'Mean').lower()
        if m == 'minimum':
            self.func = rank.minimum
        elif m == 'maximum':
            self.func = rank.maximum
        elif m == 'mean':
            self.func = rank.mean
        elif m == 'bilateral mean':
            self.func = rank.mean_bilateral
        elif m == 'median':
            self.func = rank.median
        elif m == 'sum':
            self.func = rank.sum
        elif m == 'entropy':
            self.func = rank.entropy
        elif m == 'threshold':
            self.func = rank.threshold
        elif m == 'autolevel':
            self.func = rank.autolevel
        return kwargs

def edge_confidence(epi, window=9, threshold=0.02):
    """
    Calculates the edge confidence according to eq. 2.
    This is a simple measurement which designates an edge if pixel intensities
    of gray value changes. An edge is determined by the sum of difference in
    pixel intensity of a central pixel to all pixels in a 1D window of size
    window. It is assumed that there is an edge if the sum is greater than a
    given threshold.

    Parameters
    ----------
    epi : numpy.array [v,u]
        Set of all gray-value epis for scanline s_hat.
    window_size : int, optional
        The 1D window siz in pixels. As the window should be centered
        symmetrically around the pixel to process the value shpuld be odd. For
        even numbers the next higher odd number is chosen.
    threshold : float, optional
        The threshold giving the smallest difference in EPI luminescence which
        must be overcome to designate an edge.

    Returns
    -------
    Ce : numpy.array [v,u]
         Edge confidence values for each EPI pixel.
    Me : numpy.array [v,u] of boolean.
        True means an edge was discovered for at that pixel.
    """

    v_dim = epi.shape[0]
    u_dim = epi.shape[1]
    # Check dimensions of input data
    assert epi.shape == (v_dim, u_dim,), 'Input EPI has wrong shape in function \'edge_confidence\'.'

    # Make window size odd
    if window % 2 == 0:
        warnings.warn(
            'window should be an odd number in function \'edge_confidence\'. Window size {g} was given but {u} is used instead.'.format(
                g=window, u=window + 1))
        window += 1

    # We avoid the border problem by padding the epi.
    padded_epi = np.pad(epi, ((0, 0), (int(window // 2), int(window // 2))), 'edge')
    assert padded_epi.shape == (v_dim, u_dim + window - 1,), 'Padded epi has wrong shape in function \'edge_confidence\'.'

    # Calculate confidence values
    Ce = np.zeros(epi.shape, dtype=np.float32)  # initiate array
    for k in range(window):
        Ce += (epi[...] - padded_epi[:, k:epi.shape[1] + k]) ** 2
    Me = Ce > threshold  # create confidence Mask
    Me = binary_opening(Me, selem=square(2), out=Me)  # work with square to avoid aliasing

    # Let's see if our results have reasonable meaning'
    assert np.all(
        Ce >= 0), 'Negative edge confidence found in function \'edge_confidence\'.'
    assert Ce.shape == (v_dim,
                        u_dim,), 'Ce output has incorrect shape in fucntion \'edge_confidence\'.'
    assert Me.shape == (v_dim,
                        u_dim,), 'Me output has incorrect shape in fucntion \'edge_confidence\'.'
    return Ce, Me

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 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