def circle_markers(blobs, pic_shape):
    '''Return array with circles around foci found'''

    markers_rad = np.zeros(pic_shape, dtype = np.bool)

    x_max, y_max = pic_shape

    for blob in blobs:

        x, y, r = blob

        r = r*np.sqrt(2)

        rr, cc = circle_perimeter(x, y, np.round(r).astype(int))
        rr_new, cc_new = [], []

        for x_c,y_c in zip(rr,cc):

            if (x_c >= 0) and (x_c < x_max) and (y_c >= 0) and (y_c < y_max):
                rr_new.append(x_c)
                cc_new.append(y_c)

        markers_rad[rr_new, cc_new] = True

    selem = np.array([0,1,0,1,1,1,0,1,0], dtype=bool).reshape((3,3))
    markers_rad = binary_dilation(binary_dilation(markers_rad, selem), selem)

    return markers_rad

def nlbin(im, threshold=0.5, zoom=0.5, escale=1.0, border=0.1, perc=80,
          range=20, low=5, high=90):
    """
    Performs binarization using non-linear processing.

    Args:
        im (PIL.Image):
        threshold (float):
        zoom (float): Zoom for background page estimation
        escale (float): Scale for estimating a mask over the text region
        border (float): Ignore this much of the border
        perc (int): Percentage for filters
        range (int): Range for filters
        low (int): Percentile for black estimation
        high (int): Percentile for white estimation

    Returns:
        PIL.Image containing the binarized image
    """
    if im.mode == '1':
        return im
    raw = pil2array(im)
    # rescale image to between -1 or 0 and 1
    raw = raw/np.float(np.iinfo(raw.dtype).max)
    if raw.ndim == 3:
        raw = np.mean(raw, 2)
    # perform image normalization
    if np.amax(raw) == np.amin(raw):
        raise KrakenInputException('Image is empty')
    image = raw-np.amin(raw)
    image /= np.amax(image)

    m = interpolation.zoom(image, zoom)
    m = filters.percentile_filter(m, perc, size=(range, 2))
    m = filters.percentile_filter(m, perc, size=(2, range))
    m = interpolation.zoom(m, 1.0/zoom)
    w, h = np.minimum(np.array(image.shape), np.array(m.shape))
    flat = np.clip(image[:w, :h]-m[:w, :h]+1, 0, 1)

    # estimate low and high thresholds
    d0, d1 = flat.shape
    o0, o1 = int(border*d0), int(border*d1)
    est = flat[o0:d0-o0, o1:d1-o1]
    # by default, we use only regions that contain
    # significant variance; this makes the percentile
    # based low and high estimates more reliable
    v = est-filters.gaussian_filter(est, escale*20.0)
    v = filters.gaussian_filter(v**2, escale*20.0)**0.5
    v = (v > 0.3*np.amax(v))
    v = morphology.binary_dilation(v, structure=np.ones((escale*50, 1)))
    v = morphology.binary_dilation(v, structure=np.ones((1, escale*50)))
    est = est[v]
    lo = np.percentile(est.ravel(), low)
    hi = np.percentile(est.ravel(), high)

    flat -= lo
    flat /= (hi-lo)
    flat = np.clip(flat, 0, 1)
    bin = np.array(255*(flat > threshold), 'B')
    return array2pil(bin)

def onoff(dc, attribute, roi, spectrum_name="spectrum"):

    # Extract data
    data = dc[attribute]

    # Get limits from ROI
    l, r = np.clip([roi.xmin, roi.xmax], 0, data.shape[2])
    b, t = np.clip([roi.ymin, roi.ymax], 0, data.shape[1])

    # make the 2d on spectrum mask
    mslice = np.zeros(data.shape[1:], dtype=bool)
    mslice[b:t+1,l:r+1] = True

    # make an off spectrum mask
    in_edge_pix = 4
    coreslice = binary_dilation(mslice, iterations=in_edge_pix)
    annular_width_pix = 4
    annslice = binary_dilation(coreslice, iterations=annular_width_pix)
    annslice = annslice-coreslice

    onoffspect = np.zeros(data.shape[0])

    for i, slc in enumerate(data):
        onoffspect[i] = (np.mean(slc[mslice]) -np.mean(slc[annslice]))

    # There is currently a bug that causes cubes to sometimes - but not always
    # have named coordinates.
    try:
        velocity = dc["Velocity", :,0,0]
    except IncompatibleAttribute:
        velocity = dc["World 0", :,0,0]

    return velocity, onoffspect

def add_neighbor_border_properties( segment1, segment2, vol):

	if segment1.label < segment2.label:
	# only process them once
		return

	bbx = get_bounding_box_intersection_two_segments(segment1, segment2)

	vol_crop = vol[bbx.xmin:bbx.xmax, bbx.ymin:bbx.ymax, bbx.zmin:bbx.zmax]

	CollectiveBoundingBox = namedtuple("CollectiveBoundingBox",["xmin", "xmax", "ymin", "ymax", "zmin", "zmax"])

	bbx1 = segment1.bounding_box
	bbx2 = segment2.bounding_box
	
	mask1 = binary_dilation(segment1.border_mask[bbx.xmin-bbx1.xmin : bbx.xmax-bbx1.xmin, bbx.ymin-bbx1.ymin : bbx.ymax-bbx1.ymin, bbx.zmin-bbx1.zmin : bbx.zmax-bbx1.zmin] ,structure=np.ones((3,3,1)))
	mask2 = binary_dilation(segment2.border_mask[bbx.xmin-bbx2.xmin : bbx.xmax-bbx2.xmin, bbx.ymin-bbx2.ymin : bbx.ymax-bbx2.ymin, bbx.zmin-bbx2.zmin : bbx.zmax-bbx2.zmin] ,structure=np.ones((3,3,1)))

	mask = mask1*mask2


	s1 = mask1.sum()
	s2 = mask2.sum()
	s = mask.sum()
	segment1.feature_dict["percent_border_with_neighbor"].append((segment2.label, s/float(s1) ))
	segment2.feature_dict["percent_border_with_neighbor"].append((segment1.label, s/float(s2) ))
	segment1.feature_dict["size_border_with_neighbor"].append((segment2.label, s))
	segment2.feature_dict["size_border_with_neighbor"].append((segment1.label, s))
	segment1.feature_dict["mean_intensity_border_with_neighbor"].append((segment2.label, np.sum(mask*vol_crop)/float(s)))
	segment2.feature_dict["mean_intensity_border_with_neighbor"].append((segment1.label, segment1.feature_dict["mean_intensity_border_with_neighbor"][-1][1]))
	
	score1 = segment1.feature_dict["mean_intensity_border_with_neighbor"][-1][1] * ( 1 - segment1.feature_dict["percent_border_with_neighbor"][-1][1])
	segment1.feature_dict["weighted_merge_score"].append((segment2.label, score1))
	segment2.feature_dict["weighted_merge_score"].append((segment1.label, score1))

def incorporate_cells(binary_image):    
    # invert input binary image
    inv_image = np.invert(binary_image)
    
    # matrix for binary_dilation
    struct = generate_binary_structure(2, 1)

    # do bunary dilation until the colony number even out
    plate_bin_dil = binary_dilation(inv_image, structure=struct)
    plate_dil_labels = label(plate_bin_dil)
    labels_number = len(np.unique(plate_dil_labels))  # initial number of colonies
    new_labels_number = labels_number - 1  # starting value
    cycle_number = 0  # starting value for dilation cycles
    while True:
        cycle_number += 1
        if cycle_number >= 30:
            break  # defence against infinite cycling
        else:
            if new_labels_number >= labels_number:
                break   # further dilation is useless (in theory)
            else:
                labels_number = new_labels_number
                plate_bin_dil = binary_dilation(plate_bin_dil, structure=struct)
                plate_dil_labels = label(plate_bin_dil)
                new_labels_number = len(np.unique(plate_dil_labels))
                
    return plate_bin_dil

def fix_holes(im,out,isin):

    if isin.all():
        raise ValueError('Cannot fix holes. All pixels are labeled as holes.')

    isout = isin == False

    # strel for eroding
    se = num.ones((3,3),bool)

    # store dilated version here
    isout1 = num.zeros(isout.shape,dtype=bool)

    # loop, dilating known regions
    while not isout.all():
        
        # get pixels just inside the border
        morph.binary_dilation(isout,se,output=isout1)
        border = isout1 & num.logical_not(isout)
        (yb,xb) = num.nonzero(border)
        yn = num.vstack([yb-1,
                         yb-1,
                         yb-1,
                         yb,
                         yb,
                         yb+1,
                         yb+1,
                         yb+1])
        xn = num.vstack([xb-1,
                         xb,
                         xb+1,
                         xb-1,
                         xb+1,
                         xb-1,
                         xb,
                         xb+1])
        badidx = num.logical_or(yn >= im.shape[0],
                                num.logical_or(yn < 0,
                                               num.logical_or(xn >= im.shape[1],
                                                              xn < 0)))
        yn = yn[badidx == False]
        xn = xn[badidx == False]
        #print "xn = " + str(xn)
        #print "yn = " + str(yn)
        #print "isout[yn,xn] = " + str(isout[yn,xn].astype(int))
        out[yb,xb] = num.average(out[yn,xn],axis=0,weights=isout[yn,xn].astype(float))

#        plt.subplot(121)
#        plt.imshow(isout)
#        plt.subplot(122)
#        plt.imshow(out)
#        plt.show()
        
        # swap isout, isout1x
        isout2 = isout1
        isout1 = isout
        isout = isout2

    return isin

def removeLines(img, n):
    imfft = np.fft.fft2(imggray)
    imffts = np.fft.fftshift(imfft)
    
    mags = np.abs(imffts)
    angles = np.angle(imffts)
    
    visual = np.log(mags)

    #visual2 = (visual - visual.min()) / (visual .max() - visual.min())*255
    
    #print np.mean(visual)
    visual3 = np.abs(visual.astype(np.int16) - np.mean(visual))
    
    ret = houghLines(visual3)
    ret = morphology.binary_dilation(ret )
    ret = morphology.binary_dilation(ret )
    ret = morphology.binary_dilation(ret )
    ret = morphology.binary_dilation(ret )
    ret = morphology.binary_dilation(ret )
    w,h=ret.shape
    ret[w/2-3:w/2+3, h/2-3:h/2+3]=False
    
    
    delta = np.mean(visual[ret]) - np.mean(visual)
    imsave("visual_re" + str(n) + ".jpg", visual)
    
    visual_blured = ndimage.gaussian_filter(visual, sigma=5)
    
    #visual[ret] = np.minimum(visual[ret], np.mean(visual)) # visual[vismask] - delta
    
    visual[ret] =visual_blured[ret]
    imsave("visual_ret" + str(n) + ".jpg", visual)
    
    newmagsshift = np.exp(visual)
    
    newffts = newmagsshift * np.exp(1j*angles)
    
    newfft = np.fft.ifftshift(newffts)
    
    imrev = np.fft.ifft2(newfft)
    
    newim2 =  np.abs(imrev).astype(np.uint8)
    
    #newim2[newim2<20] = 255
    newim2 = np.maximum(newim2, img)
    
    
    
    #newim2 = morphology.grey_closing(newim2, 3 )
    
    
    return newim2

def foci_markers(blobs, pic_shape):
    '''Return array with foci markers'''

    markers = np.zeros(pic_shape, dtype = np.bool)

    for blob in blobs:

        x, y, r = blob

        markers[x,y] = True

    selem = np.array([0,1,0,1,1,1,0,1,0], dtype=bool).reshape((3,3))
    markers = binary_dilation(binary_dilation(markers, selem), selem)

    return markers

def InsideOutside(s):
    '''
    Create inside-outside function for slice and extract nodes, values
    just inside, just outside and on the boundary
    
    Arguments
    ----
    s : 2D numpy integer array
        Extracted slice of label volume
    '''
    
    nx, ny = s.shape

    # Create boundary layer mask from difference between dilation
    # and erosion of label. The mask represents the layers of
    # voxels immediately inside and outside the boundary.
    bound_mask = binary_dilation(s) - binary_erosion(s)
    
    # Inside-outside function from complement Euclidian distance transforms
    # Positive outside, negative inside
    io = EDT(1-s) - EDT(s)
    
    # Extract x, y coordinates and IO function values boundary layers
    xy = np.argwhere(bound_mask) # N x 2 coordinates of non-zero voxels
    
    # Random downsample by 3
    # Every third point(ish) on boundary should be sufficient for accurate RBF
    n = xy.shape[0]
    samp = np.random.choice(np.arange(n), int(n/3.0))
    xy = xy[samp,:]
    
    io_xy = io[xy[:,0], xy[:,1]] 
    
    return io_xy, xy

	def get_neighbors(self, label_map):

		cropped_map = label_map[self.bounding_box.xmin : self.bounding_box.xmax+1, self.bounding_box.ymin : self.bounding_box.ymax +1, self.bounding_box.zmin:self.bounding_box.zmax +1]
		neighbors = np.unique( binary_dilation(self.border_mask).astype("int")* cropped_map)
		for x in neighbors:
			if not x in set((0,1,self.label)):
				self.neighbor_labels.add(x)

 def plot_maxima(self, display = False):
     vis_maxima = binary_dilation(self.maxima,
                                  structure = np.ones(shape = (1,5,5)))
     masked_maxima = np.ma.masked_where(vis_maxima == 0, vis_maxima)
     self.plot_maxima_stack(masked_maxima, self.smooth_dist_map)
     if display:
         plt.show()

def measure_of_chaos(img,nlevels): 
    ''' Function for calculating a measure of image noise using level sets
    # Inputs: image - numpy array of pixel intensities
             nlevels - number of levels to calculate over (note that we approximating a continuious distribution with a 'large enough' number of levels)
    # Outputs: measure_value
     
    This function calculates the number of connected regions above a threshold at an evenly spaced number of threshold
    between the min and max values of the image.

    There is some morphological operations to deal with orphaned pixels and heuristic parameters in the image processing.

    # Usage 
    img = misc.imread('/Users/palmer/Copy/ion_image.png').astype(float)
    print measure_of_chaos(img,20)
    '''
    # Image in/preparation
    sum_notnull = np.sum(img > 0)
    if sum_notnull == 0:
        return 0
    im=clean_image(img)

    return float(np.sum([
        ndimage.label(
            morphology.binary_erosion(
                morphology.binary_dilation(im > lev,structure=dilate_mask)
            , structure=erode_mask)
        )[1] for lev in np.linspace(np.amin(im),np.amax(im),nlevels)]))/(sum_notnull*nlevels)

    def _prep(self):
        self._nforeground = self._bimg.sum()        
        # Dilate bimg to make it less strict for the big gap criteria
        # It is needed since sometimes the tracing goes along the
        # boundary of the thin fibre in the binary img
        self._dilated_bimg = binary_dilation(self._bimg)

        if not self._silent:
            print('(2) --Boundary DT...')
        self._make_dt()
        if not self._silent:
            print('(3) --Fast Marching with %s quality...' % ('high' if self._quality else 'low'))
        self._fast_marching()
        if not self._silent:
            print('(4) --Compute Gradients...')
        self._make_grad()

        # Make copy of the timemap
        self._tt = self._t.copy()
        self._tt[self._bimg <= 0] = -2

        # Label all voxels of soma with -3
        self._tt[self._soma.mask > 0] = -3

        # For making a large tube to contain the last traced branch
        self._bb = np.zeros(shape=self._tt.shape)

	def _get_border_pixels_from_mask(self, mask_in):



		mask = np.zeros((mask_in.shape[0]+2 , mask_in.shape[1] + 2, mask_in.shape[2]+2))
		mask[1:-1,1:-1, 1:-1] = mask_in

		dilated = binary_dilation (mask)


		old_x, old_y, old_z = np.nonzero(dilated - mask)
			
		old_x = old_x - 1
		old_y = old_y - 1
		old_z = old_z - 1

		new_x = np.clip(old_x, 0, mask_in.shape[0]-1)
		new_y = np.clip(old_y, 0, mask_in.shape[1]-1)
		new_z = np.clip(old_z, 0, mask_in.shape[2]-1)


		adjusted_boundary = zip(new_x, new_y, new_z)


		return adjusted_boundary

 def __init__(self, geometric_model='affine', tps_grid_size=3, tps_reg_factor=0, h_matches=15, w_matches=15, use_conv_filter=False, dilation_filter=None, use_cuda=True, normalize_inlier_count=False, offset_factor=227/210):
     super(WeakInlierCount, self).__init__()
     self.normalize=normalize_inlier_count
     self.geometric_model = geometric_model
     self.geometricTnf = GeometricTnf(geometric_model=geometric_model,
                                      tps_grid_size=tps_grid_size,
                                      tps_reg_factor=tps_reg_factor,
                                      out_h=h_matches, out_w=w_matches,
                                      offset_factor = offset_factor,
                                      use_cuda=use_cuda)
     # define dilation filter
     if dilation_filter is None:
         dilation_filter = generate_binary_structure(2, 2)
     # define identity mask tensor (w,h are switched and will be permuted back later)
     mask_id = np.zeros((w_matches,h_matches,w_matches*h_matches))
     idx_list = list(range(0, mask_id.size, mask_id.shape[2]+1))
     mask_id.reshape((-1))[idx_list]=1
     mask_id = mask_id.swapaxes(0,1)
     # perform 2D dilation to each channel 
     if not use_conv_filter:
         if not (isinstance(dilation_filter,int) and dilation_filter==0):
             for i in range(mask_id.shape[2]):
                 mask_id[:,:,i] = binary_dilation(mask_id[:,:,i],structure=dilation_filter).astype(mask_id.dtype)
     else:
         for i in range(mask_id.shape[2]):
             flt=np.array([[1/16,1/8,1/16],
                              [1/8, 1/4, 1/8],
                              [1/16,1/8,1/16]])
             mask_id[:,:,i] = scipy.signal.convolve2d(mask_id[:,:,i], flt, mode='same', boundary='fill', fillvalue=0)
         
     # convert to PyTorch variable
     mask_id = Variable(torch.FloatTensor(mask_id).transpose(1,2).transpose(0,1).unsqueeze(0),requires_grad=False)
     self.mask_id = mask_id
     if use_cuda:
         self.mask_id = self.mask_id.cuda();

def detect_sources(snmap, threshold):
    hot = (snmap > threshold)
    hot = binary_dilation(hot, iterations=2)
    hot = binary_fill_holes(hot)
    blobs,nblobs = label(hot)
    print(nblobs, 'blobs')
    #print('blobs min', blobs.min(), 'max', blobs.max())
    slices = find_objects(blobs)
    px,py = [],[]
    for i,slc in enumerate(slices):
        blob_loc = blobs[slc]
        sn_loc = snmap[slc]
        imax = np.argmax((blob_loc == (i+1)) * sn_loc)
        y,x = np.unravel_index(imax, blob_loc.shape)
        y0,x0 = slc[0].start, slc[1].start
        px.append(x0+x)
        py.append(y0+y)
        #if i == 0:
        #    plt.subplot(2,2,1)
        #    plt.imshow(blob_loc, interpolation='nearest', origin='lower')
        #    plt.colorbar()
        #    plt.subplot(2,2,2)
        #    plt.imshow((blob_loc==(i+1))*sn_loc, interpolation='nearest', origin='lower')
        #    plt.subplot(2,2,3)
        #    plt.plot(x, y, 'ro')
    return np.array(px),np.array(py)

    def fix_hole(self,im,out,polygon,isin=None):

        s = num.ones((3,3),bool)        

        # get points on the inside of the polygon
        if isin is None:
            isin = point_inside_polygon(self.X,self.Y,polygon)
        (y_isin,x_isin) = num.nonzero(isin)

        # get points on the outside boundary of the polygon
        isboundary = num.logical_and(morph.binary_dilation(isin,s),
                                     ~isin)

        (y_isboundary,x_isboundary) = num.nonzero(isboundary)

        # for each point inside the polygon, get an average from
        # all the boundary points
        for i in range(len(y_isin)):
            x = x_isin[i]
            y = y_isin[i]
            d = num.sqrt((y_isboundary-y)**2+(x_isboundary-x)**2)
            w = num.exp(-d)
            out[y,x] = num.sum(im[isboundary] * w) / num.sum(w)

        return isin

    def erodeDilate(self, img, method=None):
        """
        Use morphological operators to erode or dilate the ridge structure.
        Dilate uses a recursive call to first dilate then erode.  Dilation
        alone produces ridge structures that are too thick to look
        authentic. Recursive call introduces random spurious minutiae when
        some valley structures are bridged.
        """
        img = np.array(img)
        if not method:
            method = random.choice(('erode', 'dilate', 'none'))
        inkIndex = np.where(img < 250)
        imgBin = np.zeros(np.shape(img))
        imgBin[inkIndex] = 1
        
        strel = morphology.generate_binary_structure(2,2)
        if method == 'erode':
            imgBin = morphology.binary_erosion(imgBin, strel)
        elif method == 'dilate':
            imgBin = morphology.binary_dilation(imgBin, strel)
        else:
            return img

        inkIndex = np.where(imgBin == 1)
        returnImg = 255*np.ones(np.shape(img))
        returnImg[inkIndex] = 0
        
        # Recursive call to erode after dilation to give more authentic
        # appearance.  Erode after dilate introduces spurious minutiae
        # but does not make the ridge structure too thick
        if method == 'dilate':
            self.erodeDilate(returnImg, method='erode')
        
        return returnImg

def getPara(predict, true, threshold, resolution, windowsize):
    (TP, FP, TN, FN, class_lable) = perf_measure(true, predict, threshold)
    if((TP + FN) == 0):
        TPR = 0
    else:
        TPR = np.float(TP) / (TP + FN)

    class_lable = class_lable.astype(bool).reshape(250,  130)
    true = true.astype(bool).reshape((250,  130))

    num = 2
    x = np.arange( -num , num+1, 1)
    xx, yy  = np.meshgrid( x, x )
    struc = (xx * xx + yy * yy)<= num * num
    class_lable = binary_dilation(class_lable, struc)
    class_lable = binary_erosion(class_lable, struc)

    # predict2 = remove_small_objects(class_lable, windowsize * resolution, in_place=False)
    predict2 = remove_small_objects(class_lable, windowsize, in_place=False)
    labeled_array1, num_features1 = label(predict2)
    labeled_array2, num_features2 = label(true)
    FP_num = num_features1 - num_features2
    if FP_num < 0:
        FP_num = 0
    return TPR, FP_num

 def plot_boundary_map(X):
     bounds = binary_dilation(X) - X
     H,W = X.shape
     rgba = np.zeros((H,W,4), np.uint8)
     rgba[:,:,1] = bounds*255
     rgba[:,:,3] = bounds*255
     plt.imshow(rgba, interpolation='nearest', origin='lower')