def analyse(im, scales=[4,10,14,40], verbose=True):
    try:
        mask = im > im.mean()
    except:
        im = im.matrix
        mask = im > im.mean()

    granulo = granulometry(mask, sizes=np.arange(2, 19, 4))
    print 'granulo:', granulo

    plt.figure(figsize=(6, 2.2))
    plt.subplot(121)
    plt.imshow(mask, cmap=plt.cm.gray, origin='lower')
    opened_less = ndimage.binary_opening(mask, structure=disk_structure(scales[0]))
    opened = ndimage.binary_opening(mask, structure=disk_structure(scales[1]))
    opened_more = ndimage.binary_opening(mask, structure=disk_structure(scales[2]))
    opened_even_more = ndimage.binary_opening(mask, structure=disk_structure(scales[3]))
    plt.contour(opened_less, [0.5], colors='g', linewidths=2)
    plt.contour(opened, [0.5], colors='b', linewidths=2)
    plt.contour(opened_more, [0.5], colors='r', linewidths=2)
    plt.contour(opened_even_more, [0.5], colors='k', linewidths=2)

    plt.axis('off')
    plt.subplot(122)
    plt.plot(np.arange(2, 19, 4), granulo, 'ok', ms=8)
    plt.subplots_adjust(wspace=0.02, hspace=0.15, top=0.95, bottom=0.15, left=0, right=0.95)
    if verbose:
        plt.show()
    
    return opened_less, opened, opened_more, opened_even_more

def open_isl(isls, crms):
    import pylab as pl
    import scipy.ndimage as nd
    import numpy as N

    thr = crms
    ft1 = N.array(((1,0,1), (0,1,0), (1,0,1)), int)
    ft2 = N.array(((0,1,0), (1,1,1), (0,1,0)), int)
    ft3 = N.ones((3,3), int)
    ft5 = N.ones((5,5), int)
    for isl in isls:
      ma = ~isl.mask_active
      open1 = nd.binary_opening(ma, ft1)
      open2 = nd.binary_opening(ma, ft2)
      open3 = nd.binary_opening(ma, ft3)
      open5 = nd.binary_opening(ma, ft5)

      pl.figure()
      pl.suptitle('Island '+str(isl.island_id))
      pl.subplot(2,2,1); pl.imshow(N.transpose(isl.image), origin='lower', interpolation='nearest'); pl.title('Image')
      pl.subplot(2,2,2); pl.imshow(N.transpose(ma), origin='lower', interpolation='nearest'); pl.title('mask')
      pl.subplot(2,2,3); pl.imshow(N.transpose(open3), origin='lower', interpolation='nearest'); pl.title('open 3x3')
      pl.subplot(2,2,4); pl.imshow(N.transpose(open5), origin='lower', interpolation='nearest'); pl.title('open 5x5')
      #pl.subplot(2,2,3); pl.imshow(N.transpose(open1), origin='lower', interpolation='nearest'); pl.title('open diag')
      #pl.subplot(2,2,4); pl.imshow(N.transpose(open2), origin='lower', interpolation='nearest'); pl.title('open str')
      pl.savefig('cyga_p_w12_bigisl_'+str(isl.island_id)+'_open.png')

def fraction_positive(bin_im, positive_im, erode=2, overlap_thresh=0.9,
                      bin_name='nuclei', positive_name='tf'):
    """Compute fraction of objects in bin_im overlapping positive_im.

    The purpose of this function is to compute the fraction of nuclei
    that express a particular transcription factor. By providing the
    thresholded DAPI channel as `bin_im` and the thresholded TF channel
    as `positive_im`, this fraction can be computed.

    Parameters
    ----------
    bin_im : 2D array of bool
        The image of objects being tested.
    positive_im : 2D array of bool
        The image of positive objects.
    erode : int, optional
        Radius of structuring element used to smooth input images.
    overlap_thresh : float, optional
        The minimum amount of overlap between an object in `bin_im` and
        the `positive_im` to consider that object "positive".
    bin_name : string, optional
        The name of the objects being tested.
    positive_name : string, optional
        The name of the property being measured.

    Returns
    -------
    f : 1D array of float, shape (1,)
        The feature vector.
    name : list of string, length 1
        The name of the feature.

    Examples
    --------
    >>> bin_im = np.array([[1, 1, 0],
    ...                    [0, 0, 0],
    ...                    [1, 1, 1]], dtype=bool)
    >>> pos_im = np.array([[1, 0, 0],
    ...                    [0, 1, 1],
    ...                    [0, 1, 1]], dtype=bool)
    >>> f = fraction_positive(bin_im, pos_im, erode=0, overlap_thresh=0.6)
    >>> f[0]
    array([0.5])
    >>> f[1][0]
    'frac-nuclei-pos-tf-erode-0-thresh-0.60'
    """
    selem = skmorph.disk(erode)
    if erode > 0:
        bin_im = nd.binary_opening(bin_im, selem)
        positive_im = nd.binary_opening(positive_im, selem)
    lab_im = nd.label(bin_im)[0].ravel()
    pos_im = positive_im.ravel().astype(int)
    counts = sparse.coo_matrix((np.ones(lab_im.size),
                                (lab_im, pos_im))).todense()
    means = counts[:, 1] / np.sum(counts, axis=1)
    f = np.array([np.mean(means[1:] > overlap_thresh)])
    name = ['frac-%s-pos-%s-erode-%i-thresh-%.2f' %
            (bin_name, positive_name, erode, overlap_thresh)]
    return f, name

  def findCom(self,data): 
    data = data.astype(int)
    if self.update_counter >= 5:
      self.update_counter = 0

      ##########################################################################
      ## Update the background image, adding a new image and removing the oldest.
      ##########################################################################
      self.background_list.insert(0,data)
      self.background_list.pop()
      background = np.zeros((480, 640, 3), dtype=int)
      for b in self.background_list:
        background += b
      self.background = background/len(self.background_list)
    
    ############################################################################
    ## Detect foreground by looking at difference from mean.
    ############################################################################
    foreground = np.sum(np.abs(np.subtract(self.background,data)),axis=2)
    falseImage = foreground
    ## clean foreground image
    falseImage[falseImage > 100] = 255
    falseImage[falseImage < 101] = 0
    falseImage = ndimage.binary_opening(falseImage)
    falseImage = ndimage.binary_closing(falseImage)
    com = ndimage.measurements.center_of_mass(falseImage)
    self.update_counter += 1
    return com

def get_difference_spots(pix):
    bpix = pix > 20
    bpix = ndimage.binary_opening(bpix)
    bpix = ndimage.binary_closing(bpix)
    labels, n = ndimage.measurements.label(bpix)
    clicks = ndimage.measurements.center_of_mass(pix, labels, range(1, n+1))
    return clicks

def detect_vortices(cloud, radius=70, showplots=False):
    """
    Detects whether there are vortex-like features within a given radius
    of the peak density in the TOF image of an expanded BEC
    """
    OD = cloud.get_OD()    
    peak_coord = cloud.results['peak coordinates']
    center_region = ROI(center=peak_coord,
                        size=(1.5 * radius, 1.5 * radius)).slices
    smooth_cloud = ndi.median_filter(OD[center_region], size=4)
    minOD = smooth_cloud.min()
    maxOD = smooth_cloud.max()
    cloud_median = ndi.median_filter(smooth_cloud, size=10)
    belowthresh = where(smooth_cloud < cloud_median * 0.75, 1, 0)
    opened = ndi.binary_opening(belowthresh, iterations=1)
    closed = ndi.binary_closing(opened, iterations=1)
    vort_found = ndi.label(closed)[1]
    cloud.results['vort_found'] = vort_found
    if showplots == True:
        fig = plt.figure(1999)
        fig.add_subplot(221, xticks=[], yticks=[])   
        plt.imshow(smooth_cloud, interpolation='nearest', vmin=minOD,
                                   vmax=maxOD)
        fig.add_subplot(222, xticks=[], yticks=[]) 
        plt.imshow(cloud_median, interpolation='nearest', vmin=minOD,
                                   vmax=maxOD)
        fig.add_subplot(223, xticks=[], yticks=[]) 
        plt.imshow(closed, interpolation='nearest',
                   cmap=plt.cm.get_cmap('binary'))
        fig.add_subplot(224, xticks=[], yticks=[]) 
        plt.imshow(belowthresh, interpolation='nearest',
                   cmap=plt.cm.get_cmap('binary'))
    return vort_found

	def cell_to_image(self):
	
		# Find x, y coordinate bounds
		x_res = max(self.pix_list, key=itemgetter(0))[0]
		y_res = max(self.pix_list, key=itemgetter(1))[1]

		# Creating labeled_img
		self.cell_img = NP.zeros([x_res+2, y_res+2], dtype=NP.int_)
		
		for (x_pix, y_pix) in self.pix_list:
			self.cell_img[x_pix-1, y_pix-1] = 1

		# Find the pixels that make up the perimeter
		eroded_image = NDI.binary_erosion(self.cell_img)

		eroded_image_open = NDI.binary_opening(eroded_image, structure=NP.ones((3,3)))
		eroded_image_open2 = NDI.binary_erosion(eroded_image_open)

		# self.perim_img = self.cell_img - eroded_image
		self.eroded_img = eroded_image_open - eroded_image_open2
		self.perim_img = self.cell_img - eroded_image

		# Create a list of the coordinates of the pixels (use the center of the pixels)
		perim_image_ind = NP.where(self.perim_img == 1)
		perim_image_coord = NP.array([perim_image_ind[0], perim_image_ind[1]])
		self.perim_coord = NP.transpose(perim_image_coord)

		return

def segmentAlg(I,frame):
    n = 10
    np.random.seed(1)

    Nx, Ny = I.size

    im = np.asarray(I)

    mask = (im > im.mean()).astype(np.float)
    mask += 1 * im


    # Thresholding algorithm
    #binary_img = white.white(I,Nx,Ny,20,1.15)#mask < 80
    binary_img = im > 100
    #scipy.misc.imsave('images/outfile'+str(frame)+'_m.jpg', binary_img)

    binary_img = ndimage.binary_opening(binary_img, structure=np.ones((2,2))).astype(np.int)
    #scipy.misc.imsave('images/outfile'+str(frame)+'_m2.jpg', binary_img)

    label_im, nb_labels = ndimage.label(binary_img)

    sizes = ndimage.sum(mask, label_im, range(nb_labels + 1))
    mask_size = sizes > 150000
    remove_pixel = mask_size[label_im]
    label_im[remove_pixel] = 0

    #scipy.misc.imsave('images/outfile'+str(frame)+'_m3.jpg', label_im)

    return binary_img, label_im, nb_labels, sizes

def binaryOpening(binarydata, structure=None, iterations=3):
    '''
    Opening [R52] is a mathematical morphology operation [R53] that consists
    in the succession of an erosion and a dilation of the input with the same
    structuring element. Opening therefore removes objects smaller than the
    structuring element.

    Together with closing (binary_closing), opening can be used for noise removal.
    '''
    result = np.empty_like(binarydata)
    
    if structure is None:
        ndimage.binary_opening(binarydata, iterations=iterations, output=result)
    else:
        ndimage.binary_opening(binarydata, structure, iterations=iterations, output=result)
    return result

def og_cb(og_msg, param_list):
    global occupancy_difference_threshold, connected_comonents_size_threshold
    rospy.loginfo('og_cb called')
    diff_og = param_list[0]
    curr_og = rog.og_msg_to_og3d(og_msg, to_binary = False)

    if diff_og == None:
        param_list[0] = curr_og
        return

    pog.subtract(diff_og, curr_og)
    param_list[0] = curr_og

    diff_og.to_binary(occupancy_difference_threshold)
    # filter the noise
    connect_structure = np.zeros((3,3,3), dtype=int)
    connect_structure[1,1,:] = 1
#    connect_structure[1,1,0] = 0
    diff_og.grid = ni.binary_opening(diff_og.grid, connect_structure,
                                     iterations = 1)
    #    diff_og.grid, n_labels = diff_og.connected_comonents(connected_comonents_size_threshold)

    print 'np.all(diff_og == 0)', np.all(diff_og.grid == 0)
    diff_og_msg = rog.og3d_to_og_msg(diff_og)
    diff_og_msg.header.frame_id = og_msg.header.frame_id
    diff_og_msg.header.stamp = og_msg.header.stamp
    param_list[1].publish(diff_og_msg)

def smoothImage(image,disk_size):
    '''
    Smooths the outline of a binary object by using morphological opening with
    a user-defined kernel size.
    '''
    smoothed = ndimage.binary_opening(image,structure=np.ones((disk_size,disk_size))).astype('uint8') * 255
    return smoothed

def find_feature_mask_simple(s_msk, sigma=1, ax=None, x_values=None):
    """
    sigma is estimated globally.
    """
    # find emission features from observed spec.

    filtered_spec = s_msk - ni.median_filter(s_msk, 15)
    #filtered_spec = ni.gaussian_filter1d(filtered_spec, 0.5)
    #smoothed_std = get_smoothed_std(filtered_spec,
    #                                rad=3, smooth_length=3)
    std = np.nanstd(filtered_spec)
    for i in [0, 1]:
        std = filtered_spec[np.abs(filtered_spec)<3*std].std()

    emission_feature_msk_ = filtered_spec > sigma*std
    #emission_feature_msk_ = ni.binary_closing(emission_feature_msk_)
    emission_feature_msk = ni.binary_opening(emission_feature_msk_,
                                             iterations=1)

    if ax is not None:
        if x_values is None:
            x_values = np.arange(len(s_msk))

        #ax.plot(x_values, s_msk)
        ax.plot(x_values, filtered_spec)
        #ax.plot(x_values, smoothed_std)
        ax.axhline(sigma*std)
        ax.plot(x_values[emission_feature_msk],
                emission_feature_msk[emission_feature_msk],
                "ys", mec="none")

    return emission_feature_msk

def granulometry(data, sizes=None):
    s = max(data.shape)
    if sizes == None:
        sizes = range(1, s/2, 2)
    granulo = [ndimage.binary_opening(data, \
            structure=disk_structure(n)).sum() for n in sizes]
    return granulo

def thresholding(img, thresh, size=9):
    """
    Segment using a thresholding algorithm
    
    Input:
     - img  ndarray : Image array (ndim=2)
     - thresh float : Threshold value for pixels selectino
     - size     int : Minimum size a group of pixels must have
    
    Output:
     - regions : Binary array for each segmented region
    
    ---
    """

    logging.debug("Threshold: %.2f", thresh)
    logging.debug("Objects min size: %d", size)

    # Take the binary image thresholded
    img_bin = img > thresh

    # And use (MO) binary opening (erosion + dilation) for cleaning spurious Trues
    strct = ndi.generate_binary_structure(2, 2)
    img_bin = ndi.binary_opening(img_bin, strct)

    # Label each group/region (value==True) of pixels
    regions, nlbl = ndi.label(img_bin)
    for i in xrange(1, nlbl + 1):
        inds = np.where(regions == i)
        if inds[0].size < size:
            regions[inds] = 0

    logging.debug("Threshold labels: %s", np.unique(regions))

    return regions.astype(np.bool)

def find_feature_mask(s_msk, sigma=1, ax=None, x_values=None):
    # find emission features from observed spec.

    filtered_spec = s_msk - ni.median_filter(s_msk, 15)
    filtered_spec = ni.gaussian_filter1d(filtered_spec, 0.5)
    smoothed_std = get_smoothed_std(filtered_spec,
                                    rad=3, smooth_length=3)
    emission_feature_msk_ = filtered_spec > sigma*smoothed_std
    #emission_feature_msk_ = ni.binary_closing(emission_feature_msk_)
    emission_feature_msk = ni.binary_opening(emission_feature_msk_,
                                             iterations=1)

    if ax is not None:
        if x_values is None:
            x_values = np.arange(len(s_msk))

        #ax.plot(x_values, s_msk)
        ax.plot(x_values, filtered_spec)
        ax.plot(x_values, smoothed_std)

        ax.plot(x_values[emission_feature_msk],
                emission_feature_msk[emission_feature_msk],
                "ys", mec="none")

    return emission_feature_msk

def artifact_mask(imdata, airdata, distance):
    """Computes a mask of artifacts found in the air region"""
    import nibabel as nb

    if not np.issubdtype(airdata.dtype, np.integer):
        airdata[airdata < .95] = 0
        airdata[airdata > 0.] = 1

    bg_img = imdata * airdata
    # Find the background threshold (the most frequently occurring value
    # excluding 0)
    # CHANGED - to the 75 percentile
    bg_threshold = np.percentile(bg_img[airdata > 0], 75)

    # Apply this threshold to the background voxels to identify voxels
    # contributing artifacts.
    qi1_img = np.zeros_like(bg_img)
    qi1_img[bg_img > bg_threshold] = 1
    qi1_img[distance < .10] = 0

    # Create a structural element to be used in an opening operation.
    struc = nd.generate_binary_structure(3, 1)
    qi1_img = nd.binary_opening(qi1_img, struc).astype(np.uint8)
    qi1_img[airdata <= 0] = 0

    return qi1_img

def artifact_mask(imdata, airdata, distance, zscore=10.):
    """Computes a mask of artifacts found in the air region"""
    from statsmodels.robust.scale import mad

    if not np.issubdtype(airdata.dtype, np.integer):
        airdata[airdata < .95] = 0
        airdata[airdata > 0.] = 1

    bg_img = imdata * airdata
    if np.sum((bg_img > 0).astype(np.uint8)) < 100:
        return np.zeros_like(airdata)

    # Find the background threshold (the most frequently occurring value
    # excluding 0)
    bg_location = np.median(bg_img[bg_img > 0])
    bg_spread = mad(bg_img[bg_img > 0])
    bg_img[bg_img > 0] -= bg_location
    bg_img[bg_img > 0] /= bg_spread

    # Apply this threshold to the background voxels to identify voxels
    # contributing artifacts.
    qi1_img = np.zeros_like(bg_img)
    qi1_img[bg_img > zscore] = 1
    qi1_img[distance < .10] = 0

    # Create a structural element to be used in an opening operation.
    struc = nd.generate_binary_structure(3, 1)
    qi1_img = nd.binary_opening(qi1_img, struc).astype(np.uint8)
    qi1_img[airdata <= 0] = 0

    return qi1_img

def findNeuron(bgImage, threshold, xNeuron,  yNeuron):
    """find bright object in small roi image."""
    mask = np.where(bgImage > threshold[0], 1, 0)
    mask = ndimage.binary_opening(mask,structure = np.ones((2,2)))
    mask = ndimage.binary_closing(mask)
        # --- Individually label all connected regions and get their center of mass
    label_im, nb_labels = ndimage.label(mask)
    centroids = ndimage.measurements.center_of_mass(bgImage, label_im, xrange(1,nb_labels+1))
    # --- select brightest object by default (mean brightness)
    meanBrightness = ndimage.measurements.mean(bgImage, label_im, xrange(1,nb_labels+1))
#        # --- Calculate the distance of each new cms to the old neuron position
#        # --- and select the new neuron position to be the object closest to
#        # --- the old location
#    dist = []
#    for coords in centroids:
#        dist.append((coords[0]-yNeuron)**2 + (coords[1]-xNeuron)**2)
#    if len(dist)==0:
#        yNewNeuron,xNewNeuron = yNeuron,  xNeuron
#    else:
#        loc = np.argmin(dist)
#        yNewNeuron,xNewNeuron = centroids[loc]
    if nb_labels >1:
        loc = np.argmax(meanBrightness)
        yNewNeuron,xNewNeuron = centroids[loc]
    else:
        yNewNeuron,xNewNeuron = yNeuron,  xNeuron
        loc = -1  
    neuronObject = np.where(label_im == loc+1,0,1)
    neuronArea = np.sum(neuronObject)
        # --- Get average of the neuron fluoresence --- 
    tmp_neuron = np.ma.masked_array(bgImage, neuronObject)
    newNeuronAverage = np.ma.average(tmp_neuron[tmp_neuron>threshold[1]])
        
    return yNewNeuron,xNewNeuron, newNeuronAverage,neuronArea,  neuronObject

def detect_current(cloud, showplots=False):
    """
    Detects whether there is a vortex-like signature of persistent
    current in the center of a TOF image of an expanded ring BEC
    """
    OD = cloud.get_OD()    
    peak_coord = cloud.results['peak coordinates']
    center_region = ROI(center=peak_coord, size=(40, 40)).slices
    cloud_center = ndi.median_filter(OD[center_region], size=2)
    minOD = cloud_center.min()
    maxOD = cloud_center.max()
    cloud_median = ndi.median_filter(cloud_center, size=10)
    belowthresh = where(cloud_center < cloud_median * 0.75, 1, 0)
    opened = ndi.binary_opening(belowthresh, iterations=1)
    closed = ndi.binary_closing(opened, iterations=3)
    current_found = ndi.label(closed)[1]
    cloud.results['current_found'] = current_found
    if showplots == True:
        fig = plt.figure(1999)
        fig.add_subplot(221, xticks=[], yticks=[])   
        plt.imshow(cloud_center, interpolation='nearest', vmin=minOD,
                                   vmax=maxOD)
        fig.add_subplot(222, xticks=[], yticks=[]) 
        plt.imshow(cloud_median, interpolation='nearest', vmin=minOD,
                                   vmax=maxOD)
        fig.add_subplot(223, xticks=[], yticks=[]) 
        plt.imshow(closed, interpolation='nearest',
                   cmap=plt.cm.get_cmap('binary'))
        fig.add_subplot(224, xticks=[], yticks=[]) 
        plt.imshow(belowthresh, interpolation='nearest',
                   cmap=plt.cm.get_cmap('binary'))
    return current_found, asum(closed)     

def post_process_prediction(prediction):
    # find connected components and remove small ones
    # create structure element (ball)
    str_el = create_ball_3d(3)
    img_2 = ndimage.binary_opening(prediction >= 2, str_el)
    connected_components, n_components = ndimage.label(img_2) # 0 and 1 are background/brain
    discard_components = []
    component_sizes = []
    all_components = np.arange(n_components)+1
    for component in all_components:
        size_of_component = np.sum(connected_components == component)
        if size_of_component < 3000:
            discard_components.append(component)
        component_sizes.append(size_of_component)
    if len(discard_components) == n_components:
        discard_components = discard_components[discard_components!=np.argmax(component_sizes)]
    keep_components = [i for i in all_components if i not in discard_components]
    new_mask = np.zeros(prediction.shape, dtype=bool)
    for keep_me in keep_components:
        mask = ndimage.binary_dilation(connected_components == keep_me, create_ball_3d(5))
        new_mask = (new_mask | mask)
    prediction_cleaned = np.zeros(prediction.shape, dtype=np.int32)
    prediction_cleaned[new_mask] = prediction[new_mask]
    prediction_cleaned[prediction_cleaned == 1] = 0
    prediction_cleaned[prediction_cleaned > 0] -= 1
    return prediction_cleaned