def log(self):
     """Laplacian of Gaussian."""
     # skimage.feature.blob_log(image, min_sigma=1, max_sigma=50,
     #     num_sigma=10, threshold=0.2, overlap=0.5, log_scale=False)
     blobs = feature.blob_log(self.image, **self.blob_ka, **self.log_ka)
     blobs[:, 2] = blobs[:, 2] * sqrt(2)  # Compute radii in 3rd column.
     return blobs

def psf_finder(stack, dominant_z, px = 100, wl = 700, patch_size = 48):
    '''
    find psfs and their centers from a stack.
    patch_size: the size of each psf small stack.
    '''
    hps = int(patch_size //2)
    NY, NX = stack[0].shape
    dominant_slice = stack[dominant_z]
    min_sig = 0.61*wl/px #theoretical diffraction-limited psf size at focus
    max_sig = 1.5* min_sig # this is kinda random
    blobs = blob_log(dominant_slice, min_sigma = min_sig, max_sigma = max_sig, threshold = 40 )
    centers = blobs[:,:2] # blob_centers
    lower_rest = np.logical_and(centers[:,0] > hps, centers[:,1] > hps)
    upper_rest = np.logical_and(centers[:,0] < NY - hps, centers[:,1] < NX-hps)
    total_rest = np.logical_and(lower_rest, upper_rest)
    centers = centers[total_rest]

    ind_accept = psf_isolation(centers, 30)
    centers = centers[ind_accept]
    print(centers)


    psf_collection = []
    for cc in centers:
        psf_patch = stack[:, cc[0]-hps:cc[0]+hps, cc[1]-hps:cc[1]+hps]
        psf_collection.append(psf_patch)

    return psf_collection, centers

def label_colonies(gray, min_foci_radius = 50, max_foci_radius = 200, \
        overlap=0, log_thres = 0.04, scale = 4):
    '''Label colonies on the image'''

    binary = (1 - binarize(gray))*255.

    min_sigma = ((min_foci_radius/3.)*2)/scale
    max_sigma = ((max_foci_radius/3.)*2)/scale
#    num_sigma = np.floor(max_sigma - min_sigma).astype(int)/10 + 1
    num_sigma = 10

    if scale != 1:
        new_shape = np.floor(np.array(gray.shape)/np.float(scale)).astype(np.int)
#        print new_shape, min_sigma, max_sigma
        small_im  = imresize(binary, new_shape)
    else:
        small_im = binary

    blobs_log = blob_log(small_im, min_sigma=min_sigma, max_sigma=max_sigma,\
            num_sigma=num_sigma, threshold=log_thres, overlap = overlap/100.)

    markers_num = blobs_log.shape[0]

    blobs_log = np.floor(blobs_log*np.float(scale)).astype(np.int)

    markers_fin = circle_markers(blobs_log, gray.shape)

    circles = np.copy(gray)
    circles[markers_fin] = 255

    return [markers_num, circles, blobs_log]

    def circles(self, filename):
        image = cv2.imread(filename, 0)
        image_gray = rgb2gray(image)

        blobs_log = blob_log(image_gray, max_sigma=30, num_sigma=10, threshold=.1)

        # Compute radii in the 3rd column.
        blobs_log[:, 2] = blobs_log[:, 2] * sqrt(2)

        blobs_dog = blob_dog(image_gray, max_sigma=30, threshold=.1)
        blobs_dog[:, 2] = blobs_dog[:, 2] * sqrt(2)

        blobs_doh = blob_doh(image_gray, max_sigma=30, threshold=.01)

        blobs_list = [blobs_log, blobs_doh]
        colors = ['yellow', 'red']
        titles = ['Laplacian of Gaussian', 'Determinant of Hessian']
        sequence = zip(blobs_list, colors, titles)

        fig, axes = plt.subplots(1, 2, sharex=True, sharey=True, subplot_kw={'adjustable': 'box-forced'})
        axes = axes.ravel()

        for blobs, color, title in sequence:
            ax = axes[0]
            axes = axes[1:]
            ax.set_title(title)
            ax.imshow(image, interpolation='nearest')
            for blob in blobs:
                y, x, r = blob
                c = plt.Circle((x, y), r, color=color, linewidth=2, fill=False)
                ax.add_patch(c)

        plt.savefig('output.png')
        plt.show()

def blob_detection(data, min_sigma=1, max_sigma=50, num_sigma=10, threshold=0.2, overlap=0.5):
    """Finds blobs in the given image.
    See also http://scikit-image.org/docs/dev/api/skimage.feature.html#skimage.feature.blob_log

    Args:
        data (ndarray): An image data.
        min_sigma (float, optional): The minimum standard deviation.
            Keep this low to detect smaller blobs. Defaults to 1.
        max_sigma (float, optional): The maximum standard deviation.
            Keep this high to detect larger blobs. Defaults to 50.
        num_sigma (int, optional): The number of intermediate values between `min_sigma` and `max_sigma`.
            Defaults to 10.
        threshold (float, optional): The absolute lower bound for scale space maxima.
            Reduce this to detect blobs with less intensities.
        overlap (float, optional): A value between 0 and 1.

    Returns:
        ndarray: Blobs detected.
            Each row represents coordinates and the standard deviation, `(x, y, r)`.

    """
    try:
        from skimage.feature import blob_log
    except ImportError:
        raise ImportError("No module named 'skimage'. 'spot_detection' requires 'scikit-image'")

    ## Laplacian of Gaussian
    blobs = blob_log(
        data, min_sigma=min_sigma, max_sigma=max_sigma, num_sigma=num_sigma, threshold=threshold, overlap=overlap)
    blobs[: , 2] = blobs[: , 2] * numpy.sqrt(2)
    _log.info('{} blob(s) were detected'.format(len(blobs)))
    return blobs

    def view_U_matrix(self,distance2=4,row_normalized='Yes',show_data='Yes',contooor='Yes',blob = 'Yes',save='Yes',save_dir = ''):
    	import scipy
    	from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show
    	umat = self.U_matrix(distance=distance2,row_normalized=row_normalized) 
    	data = getattr(self, 'data_raw')
    	proj = self.project_data(data)
    	msz =  getattr(self, 'mapsize')
    	coord = self.ind_to_xy(proj)
	#     freq = plt.hist2d(coord[:,1], coord[:,0], bins=(msz[1],msz[0]),alpha=1.0,cmap=cm.jet)[0]
	#     plt.close()
   	 
	#     fig, ax = plt.figure()
    	fig, ax= plt.subplots(1, 1)
    	im = imshow(umat,cmap=cm.RdYlBu_r,alpha=1) # drawing the function
    	# adding the Contour lines with labels`
    	# imshow(freq[0].T,cmap=cm.jet_r,alpha=1)
    	if contooor=='Yes':
        	mn = np.min(umat.flatten())
        	mx = np.max(umat.flatten())
        	std = np.std(umat.flatten())
        	md = np.median(umat.flatten())
        	mx = md + 0*std
#         	mn = md
#         	umat[umat<=mn]=mn
        	cset = contour(umat,np.linspace(mn,mx,15),linewidths=0.7,cmap=cm.Blues)
    
    	if show_data=='Yes':
        	plt.scatter(coord[:,1], coord[:,0], s=2, alpha=1.,c='Gray',marker='o',cmap='jet',linewidths=3, edgecolor = 'Gray')
        	plt.axis('off')
    
    	ratio = float(msz[0])/(msz[0]+msz[1])
    	fig.set_size_inches((1-ratio)*15,ratio*15)
    	plt.tight_layout()
    	plt.subplots_adjust(hspace = .00,wspace=.000)
    	sel_points = list()
    	if blob=='Yes':
        	from skimage.feature import blob_dog, blob_log, blob_doh
        	from math import sqrt
        	from skimage.color import rgb2gray
        	image = 1/umat
        	image_gray = rgb2gray(image)

        	#'Laplacian of Gaussian'
        	blobs = blob_log(image, max_sigma=5, num_sigma=4, threshold=.152)
        	blobs[:, 2] = blobs[:, 2] * sqrt(2)
        	imshow(umat,cmap=cm.RdYlBu_r,alpha=1)
        	sel_points = list()
        	for blob in blobs:
        		row, col, r = blob
        		c = plt.Circle((col, row), r, color='red', linewidth=2, fill=False)
        		ax.add_patch(c)
        		dist = scipy.spatial.distance_matrix(coord[:,:2],np.array([row,col])[np.newaxis,:])
        		sel_point = dist <= r
        		plt.plot(coord[:,1][sel_point[:,0]], coord[:,0][sel_point[:,0]],'.r')
        		sel_points.append(sel_point[:,0])

            
        if save=='Yes':
        	fig.savefig(save_dir, transparent=False, dpi=400) 
        return sel_points,umat

def detect_cells(image):
    image_gray = rgb2gray(image)

    blobs_log = blob_log(image_gray, max_sigma=30, num_sigma=10, threshold=.1)
    # Compute radii in the 3rd column.
    blobs_log[:, 2] = blobs_log[:, 2] * sqrt(2)

    blobs_dog = blob_dog(image_gray, max_sigma=30, threshold=.1)
    blobs_dog[:, 2] = blobs_dog[:, 2] * sqrt(2)

    blobs_doh = blob_doh(image_gray, max_sigma=30, threshold=.01)

    blobs_list = [blobs_log, blobs_dog, blobs_doh]
    colors = ['yellow', 'lime', 'red']
    titles = ['Laplacian of Gaussian', 'Difference of Gaussian',
          'Determinant of Hessian']
    sequence = zip(blobs_list, colors, titles)

    for blobs, color, title in sequence:
        fig, ax = plt.subplots(1, 1)
        ax.set_title(title)
        ax.imshow(image, interpolation='nearest')
        for blob in blobs:
            y, x, r = blob
            c = plt.Circle((x, y), r, color=color, linewidth=2, fill=False)
            ax.add_patch(c)

    plt.show()

 def get_blobs(self, image, **kwargs):
     '''
     Use Laplacian of Gaussian to find blobs in passed grayscale image.
     '''
     blobs = blob_log(image, **kwargs)
     # Compute radii in the 3rd column.
     blobs[:, 2] = blobs[:, 2] * sqrt(2)
     return blobs

def get_number_of_blobs(image):
    from skimage.feature import blob_dog, blob_log, blob_doh
    from math import sqrt
    from skimage.color import rgb2gray
    image_gray = rgb2gray(image)
    blobs_log = blob_log(image_gray, max_sigma=30, num_sigma=10, threshold=.1)

    return blobs_log.size

def blob_image_multiscale2(image, type=0,scale=2):
    # function that return a list of blob_coordinates, 0 = dog, 1 = doh, 2 = log
    list = []
    image = norm.normalize(image)
    for z, slice in tqdm(enumerate(image)):
        # init list of different sigma/zoom blobs
        featureblobs = []
        # x = 0,1,2,3,4
        if scale == 2:
            # for x in xrange(0,6):
            #     if type == 0:
            #         featureblobs.append(feature.blob_dog(slice, 2**x, 2**x))
            #     if type == 1:
            #         featureblobs.append(feature.blob_doh(slice, 2**x, 2**x))
            #     if type == 2:
            #         featureblobs.append(feature.blob_log(slice, 2**x, 2**x))
            for x in xrange(0,5):
                if type == 0:
                    featureblobs.append(feature.blob_dog(slice, 2**x, 2**(x+1)))
                if type == 1:
                    featureblobs.append(feature.blob_doh(slice, 2**x, 2**(x+1)))
                if type == 2:
                    featureblobs.append(feature.blob_log(slice, 2**x, 2**(x+1),16,.1))
        else:
            for x in xrange(0,4):
                if type == 0:
                    featureblobs.append(feature.blob_dog(slice, 3**x, 3**x))
                if type == 1:
                    featureblobs.append(feature.blob_doh(slice, 3**x, 3**x))
                if type == 2:
                    featureblobs.append(feature.blob_log(slice, 3**x, 3**x))
        # init list of blob coords
        blob_coords = []
        #print featureblobs
        # start at biggest blob size
        for featureblob in reversed(featureblobs):
            # for every blob found of a blobsize

            for blob in enumerate(featureblob):
                # if that blob is not within range of another blob, add it
                blob = blob[1]
                if not within_range(blob, blob_coords):
                    blob_coords.append([z, blob[0], blob[1], blob[2]])
        list.append(blob_coords)
    return list

    def show(self, som, distance2=1, row_normalized=False, show_data=True, contooor=True, blob=False, labels = False):
        umat = self.build_u_matrix(som, distance=distance2, row_normalized=row_normalized)
        msz = som.codebook.mapsize
        proj = som.project_data(som.data_raw)
        coord = som.bmu_ind_to_xy(proj)
        
        fig, ax = plt.subplots(1, 1)
        im = imshow(umat, cmap=plt.cm.get_cmap('RdYlBu_r'), alpha=1)

        if contooor:
            mn = np.min(umat.flatten())
            mx = np.max(umat.flatten())
            std = np.std(umat.flatten())
            md = np.median(umat.flatten())
            mx = md + 0*std
            cset = contour(umat, np.linspace(mn, mx, 15), linewidths=0.7, cmap=plt.cm.get_cmap('Blues'))

        if show_data:
            plt.scatter(coord[:, 1], coord[:, 0], s=2, alpha=1., c='Gray', marker='o', cmap='jet', linewidths=3, edgecolor='Gray')
            plt.axis('off')
            
        if labels:
            if labels == True:
              labels = som.build_data_labels()
            for label, x, y in zip(labels, coord[:, 1], coord[:, 0]):
                plt.annotate(str(label), xy = (x, y), horizontalalignment = 'center', verticalalignment = 'center')

        ratio = float(msz[0])/(msz[0]+msz[1])
        fig.set_size_inches((1-ratio)*15, ratio*15)
        plt.tight_layout()
        plt.subplots_adjust(hspace=.00, wspace=.000)
        sel_points = list()

        if blob:
            from skimage.color import rgb2gray
            from skimage.feature import blob_log

            image = 1/umat
            image_gray = rgb2gray(image)

            #'Laplacian of Gaussian'
            blobs = blob_log(image, max_sigma=5, num_sigma=4, threshold=.152)
            blobs[:, 2] = blobs[:, 2] * sqrt(2)
            imshow(umat, cmap=plt.cm.get_cmap('RdYlBu_r'), alpha=1)
            sel_points = list()

            for blob in blobs:
                row, col, r = blob
                c = plt.Circle((col, row), r, color='red', linewidth=2, fill=False)
                ax.add_patch(c)
                dist = scipy.spatial.distance_matrix(coord[:, :2], np.array([row, col])[np.newaxis, :])
                sel_point = dist <= r
                plt.plot(coord[:, 1][sel_point[:, 0]], coord[:, 0][sel_point[:, 0]], '.r')
                sel_points.append(sel_point[:, 0])

        plt.show()
        return sel_points, umat

def _detect_spots_blob_log(image, minSpotSize):
	
	maxSpotSize = minSpotSize
	spotSizeSteps = 1
	
	# find blobs starting from min_sigma to max_sigma in num_sigma steps
	blobs = blob_log(image, min_sigma=minSpotSize, max_sigma=maxSpotSize, num_sigma=spotSizeSteps, threshold=.005) # blobs_log = (y, x, r)
	
	return blobs

 def get_blobs(self, image):
   blobs_log = blob_log(image, max_sigma=30, num_sigma=10, threshold=.1)
   blobs_log[:, 2] = blobs_log[:, 2] * sqrt(2)
   blobs_dog = blob_dog(image, max_sigma=30, threshold=.1)
   blobs_dog[:, 2] = blobs_dog[:, 2] * sqrt(2)
   blobs_doh = blob_doh(image, max_sigma=30, threshold=.01)
   
   all_blobs = np.vstack([blobs_log, blobs_doh, blobs_dog])
   all_blobs = filter(lambda b: b[2] > 4, all_blobs)
   all_blobs = list(filter(lambda b: b[2] < 60, all_blobs))
   return all_blobs

def skimage_blob(frame):
    # gray_frm = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)

    detected_dogs = feature.blob_log(frame)

    for blob in detected_dogs:
        sigma = blob[-1]
        blob_rad = (2 ** 1 / 2) * sigma
        cv2.circle(frame, (blob[1], blob[0]), blob_rad, (255, 0, 0))

    return frame

 def predict_test(self, test_folder, destination_folder):
   test_image_file = self.get_image_names(test_folder)
   for i, im in enumerate(test_image_file):
     print 'Processing Test Image:', i
     file_name = os.path.join(test_folder, im)
     image = imread(file_name)
     im_final = self.apply_gaussian_filter(image)
     image_gray = rgb2gray(im_final)
     blobs = blob_log(image_gray, max_sigma=30, num_sigma=10, threshold=.55)
     blob_list = self.process_blobs(image, blobs)
     self.create_mask(image, im, blob_list, destination_folder)	  

def temblob(image,ind):

	"""
	Laplacian of gaussian blob detection for TEM images
	"""

	org = image[4:-256,4:-4]


	with warnings.catch_warnings():
		warnings.simplefilter("ignore")
		warnings.warn("user", UserWarning)

		igray = img_as_ubyte(rgb2gray(org))
		iinv = np.invert(igray)


	igaus = img_as_float(iinv)
	igaus = gaussian_filter(igaus, 1)
	h = 0.5
	sd = igaus - h
	msk = igaus
	dilat = reconstruction(sd, msk, method='dilation')
	hdome = igaus - dilat

	if ind == 'AgNP':
		kwargs = {}
		kwargs['threshold'] = 0.01
		kwargs['overlap'] = 0.4
		kwargs['min_sigma'] = 25
		kwargs['max_sigma'] = 50
		kwargs['num_sigma'] = 25
		calib = 500/(969-26)
		
	elif ind == 'AuNP':
		kwargs = {}
		kwargs['threshold'] = 0.01
		kwargs['overlap'] = 0.4
		kwargs['min_sigma'] = 18
		kwargs['max_sigma'] = 23
		kwargs['num_sigma'] = 5
		calib = 200/(777-23)
		
	else:
		warnmsg='Unable to identify keyword: {:}'.format(ind)
		warnings.warn(warnmsg, UserWarning)

	blobs = blob_log(hdome, **kwargs)
	diam = 2*sqrt(2)*blobs[:,-1]
	npdiam = [ind]
	npdiam.extend(calib*diam)

	return(npdiam)

    def image_blobs(self, n_frame):
        """
        input: number of frame
        """
        im0 = self.stack[n_frame]
        mx_sig = self.blobset[0]
        mi_sig = self.blobset[1]
        nsig = self.blobset[2]
#         th = (np.max(im0)-np.min(im0))/10. # threshold
        th = (np.mean(im0)-np.min(im0))/12.
        print("threshold:", th)
        self.c_list[n_frame] = blob_log(im0, 
            max_sigma = mx_sig, min_sigma = mi_sig, num_sigma=nsig, threshold = th, overlap = OL_blob)
        self.bl_flag[n_frame] = self.c_list[n_frame].shape[0]

def blob_feature(image, method='log'):
    if method == 'log':
        blobs = blob_log(image, )
    else:
        blobs = blob_doh(image, )

    blob_image = np.zeros(image.shape)

    #Draw the blobs to an image
    for blob in blobs:
        y,x,sigma = blob
        color = sigma
        size = int(np.sqrt(2*sigma)) if method == 'log' else sigma
        cv2.circle(blob_image, (x, y), size, sigma/1,-1)

    #dataset.show_image(blob_image)
    return blob_image

def run(min_s, max_s, sigma):
	# Run Blob Detection by Laplacian of Gaussian (LoG)
	blobs_log = blob_log(image, max_s, threshold=.01)

	# Compute radii in the 3rd column
	blobs_log[:, 2] = blobs_log[:, 2] * sqrt(2)

	# Configure the graph
	fig, ax = pyplot.subplots(1, 1)
	ax.set_title("Laplacian of Gaussian, min_s=" + str(min_s) + ", max_s="
			 + str(max_s) + ", sigma=" + str(sigma))
	ax.imshow(image, vmin=0, vmax=255, cmap=pyplot.cm.gray)

	# Load graph
	for blob in blobs_log:
		y, x, r = blob
		c = pyplot.Circle((x, y), r, color='yellow', linewidth=2, fill=False)
		ax.add_patch(c)

def Blobing(image, min_sigma, max_sigma, threshold):
    #image = mpimg.imread('stepping.png')[0:500, 0:500] 
    image_gray = rgb2gray(image)
    image_fltrd = canny(image_gray,
                        sigma=1.0,
                        low_threshold=None,
                        high_threshold=None,
                        mask=None)
    
    blobs_log = blob_log(image_fltrd,
                         min_sigma=min_sigma, # set to10
                         max_sigma=max_sigma, # set to 40
                         num_sigma=10,
                         threshold=threshold)
    
    #-----------------------------------------------------
    #----- Calculating the average of all blobs
    #      This will issue the location of each step
    lngth = len(blobs_log[0:])
    x_avg = sum(blobs_log[:,0])/lngth
    y_avg = sum(blobs_log[:,1])/lngth
    #print 'Step coordinates: (x,y):''(',x_avg,',', y_avg,')'
    
    # Compute belowradii in the 3rd column, this is just the radial of each circle
    blobs_log[:, 2] = blobs_log[:, 2]*sqrt(2)
    
    #-----------------------------------------------------
    #----- Ploting
    fig, ax = plt.subplots(1, 1)
    ax.imshow(image, interpolation='gaussian')
    #time.sleep(0.5)
    plt.close(fig) # save time by closing figure
 
    for circles in blobs_log:
        y, x, r = circles
        c = plt.Circle((x, y), r, color='blue', linewidth=2, fill=False)
        ax.add_patch(c)
        #print x, y
    plt.show()
    #time.sleep(0.5)
    plt.close(fig) # close figure
    
    return x_avg, y_avg # The Blobing returns the coordinates to be used for move command isuers