diff --git a/src/bactos_functions.py b/src/bactos_functions.py index cef843848a27560a8873091fee6d3a1b199648a4..2200a21ad10cc74a36a3bb6554d7d8cb31bb2734 100755 --- a/src/bactos_functions.py +++ b/src/bactos_functions.py @@ -9,15 +9,16 @@ import re import json import glob import numpy as np +from pathlib import Path from tifffile import imread -from skimage.feature import match_template +from skimage.feature import match_template, match_descriptors, SIFT from skimage.filters import gaussian, threshold_yen from skimage.restoration import denoise_wavelet, estimate_sigma from skimage.registration import phase_cross_correlation from skimage.morphology import disk, white_tophat, opening -from skimage.measure import label, regionprops, find_contours +from skimage.measure import label, regionprops, find_contours, ransac from skimage.segmentation import clear_border, expand_labels -from skimage import transform +from skimage import transform, exposure from skimage.util import img_as_float from scipy import ndimage as ndi @@ -25,10 +26,14 @@ import dateutil import h5py import datetime import os +import io from joblib import Parallel, delayed from csbdeep.utils import normalize from stardist.models import StarDist2D, Config2D +from lib_bactos_hdf import load_image, load_mask, load_intensity_profile +import matplotlib.pyplot as plt +from rich.progress import track def load_MM2_images(image_root_path, channel='*', position='*', time='*', z='*'): @@ -285,11 +290,57 @@ def register_task(image, refimage, gaussian_size=0, up_factor=20, show_info=Fals return dxdyt +def register_task_sift(image, refimage, remove_cosmics=True, gaussian_size=0): + """ + Use SIFT descriptor and RANSAC to obtain image shift + """ + + if remove_cosmics: + el = disk(2) + refimage = opening(refimage, el) + image = opening(image, el) + + if gaussian_size > 0: + image = gaussian(image, gaussian_size) + refimage = gaussian(refimage, gaussian_size) + + dstimg = exposure.rescale_intensity(image) + srcimg = exposure.rescale_intensity(refimage) + + # SIFT detector + detect1 = SIFT() + detect2 = SIFT() + + detect1.detect_and_extract(srcimg) + detect2.detect_and_extract(dstimg) + + # Match detected points between the two images + try: + matched = match_descriptors(detect1.descriptors, detect2.descriptors, + max_ratio=0.6) + src = detect1.keypoints[matched[:, 0]] + dst = detect2.keypoints[matched[:, 1]] + + # Use ransac to find the transform matrix + matrix, _ = ransac((src, dst), transform.EuclideanTransform, + min_samples=10, + residual_threshold=1) + + # print(matrix.translation) + out = -matrix.translation + except: + out = (0, 0) + + return out + + def register_images(imgs_stack, gaussian_size = 0, - register_filter=False, ref_filter=1, + register_filter=False, register_stack=True, ref_filter=1, move_filter=0, show_info=False, - up_factor=20, drift_filter=None, remove_cosmics=True, - scale_correction=False): + up_factor=20, drift_filter=None, + remove_cosmics=True, + scale_correction=False, + drift_algo='cor'): """ Register images to correct drift in time and between filter @@ -297,8 +348,20 @@ def register_images(imgs_stack, gaussian_size = 0, remove_cosmics: True, use an opening filter of disk size(2) to remove cosmics on images before registration - + register_filter: boolean, + Do we need to register the other filter to the ref_filter + + register_stack: boolean, + Do the registration of the whole stack + + scale_correction: boolean, + Try to correct the scaling between filters (default False) + + drift_algo: str in ['cor', 'sift'], + Define the algorithm to compute the drift between images. + 'cor': use crosscorrelation (the default) + 'sift': use SIFT point finding + RANSAC to find transform matrix """ # Do we need to do scale correction @@ -311,47 +374,65 @@ def register_images(imgs_stack, gaussian_size = 0, filters = [drift_filter] # A first drift correction in time - for pos in range(imgs_stack.shape[1]): - dxdy = [] - for img_filter in filters: - if show_info: - print('stabilisation de la position pos %i' % pos) - print('stabilisation du filtre %i' % img_filter) - - imgref = imgs_stack[0, pos, 0, 5:-5, 5:-5, img_filter] - - if remove_cosmics: - el = disk(2) - imgref = opening(imgref, el) - - if gaussian_size > 0: - imgref = gaussian(imgref, gaussian_size) - - # Non parrallel version - # dxdyt = [register_task(imgs_stack[t,pos,0,:,:,img_filter], imgref, up_factor=up_factor) for t in range(1, len(imgs_stack))] - # Parralel version - if remove_cosmics: - dxdyt = Parallel(n_jobs=-1)(delayed(register_task)(opening(imgs_stack[t,pos,0, 5:-5, 5:-5, img_filter], el), - imgref, up_factor=up_factor) for t in range(1, len(imgs_stack))) - else: - dxdyt = Parallel(n_jobs=-1)(delayed(register_task)(imgs_stack[t,pos,0, 5:-5, 5:-5, img_filter], imgref, up_factor=up_factor) for t in range(1, len(imgs_stack))) - dxdy += dxdyt - - if drift_filter is None: - for i, t in enumerate(range(1, len(imgs_stack))): - imgs_stack[t, pos, 0,:,:,img_filter] = ndi.shift(imgs_stack[t,pos,0,:,:,img_filter], dxdyt[i]) - - if drift_filter is not None: - for img_filter in range(imgs_stack.shape[5]): + if register_stack: + print('Start stack registration !') + for pos in range(imgs_stack.shape[1]): + dxdy = [] + for img_filter in filters: if show_info: print('stabilisation de la position pos %i' % pos) - print('stabilisation du filtre %i avec le déplacement filtre %i' % (img_filter, drift_filter)) - for ti, t in enumerate(range(1, len(imgs_stack))): - imgs_stack[t, pos, 0,:,:,img_filter] = ndi.shift(imgs_stack[t,pos,0,:,:,img_filter], dxdy[ti]) - + print('stabilisation du filtre %i' % img_filter) + + imgref = imgs_stack[0, pos, 0, 5:-5, 5:-5, img_filter] + + if drift_algo == 'sift': + print('Use SIFT to align stack') + try: + dxdyt = Parallel(n_jobs=-1, backend='multiprocessing')(delayed(register_task_sift)(imgs_stack[t,pos,0, 5:-5, 5:-5, img_filter], + imgs_stack[t-1,pos,0, 5:-5, 5:-5, img_filter], remove_cosmics=remove_cosmics) for t in range(1, len(imgs_stack))) + dxdyt = np.cumsum(dxdyt, axis=0).tolist() + + except: + print('Error in SIFT, failback to crosscorrelation registration') + drift_algo = 'cor' + + if drift_algo == 'cor': + print('Use crosscorrelation algorithm to align stack') + if remove_cosmics: + el = disk(2) + imgref = opening(imgref, el) + + if gaussian_size > 0: + imgref = gaussian(imgref, gaussian_size) + + # Non parrallel version + # dxdyt = [register_task(imgs_stack[t,pos,0,:,:,img_filter], imgref, up_factor=up_factor) for t in range(1, len(imgs_stack))] + # Parralel version + if remove_cosmics: + dxdyt = Parallel(n_jobs=-1, backend='multiprocessing')(delayed(register_task)(opening(imgs_stack[t,pos,0, 5:-5, 5:-5, img_filter], el), + imgref, up_factor=up_factor) for t in range(1, len(imgs_stack))) + else: + dxdyt = Parallel(n_jobs=-1, backend='multiprocessing')(delayed(register_task)(imgs_stack[t,pos,0, 5:-5, 5:-5, img_filter], + imgref, up_factor=up_factor) for t in range(1, len(imgs_stack))) + + dxdy += dxdyt + + if drift_filter is None: + for i, t in enumerate(range(1, len(imgs_stack))): + imgs_stack[t, pos, 0,:,:,img_filter] = ndi.shift(imgs_stack[t,pos,0,:,:,img_filter], dxdyt[i]) + + if drift_filter is not None: + for img_filter in range(imgs_stack.shape[5]): + if show_info: + print('stabilisation de la position pos %i' % pos) + print('stabilisation du filtre %i avec le déplacement filtre %i' % (img_filter, drift_filter)) + for ti, t in enumerate(range(1, len(imgs_stack))): + imgs_stack[t, pos, 0,:,:,img_filter] = ndi.shift(imgs_stack[t,pos,0,:,:,img_filter], dxdy[ti]) + # A second drift correction between filters, only on the first time if register_filter: + print('Register filters between two filters !') filter_ref = ref_filter filter_move = move_filter for pos in range(imgs_stack.shape[1]): @@ -369,8 +450,8 @@ def register_images(imgs_stack, gaussian_size = 0, if show_info: print('Error: %0.2e' % (error)) - for t in range(0, len(imgs_stack)): - imgs_stack[t,pos,0,:,:,filter_move] = ndi.shift(imgs_stack[t,pos,0,:,:,filter_move], dxdyt) + for t in range(0, len(imgs_stack)): + imgs_stack[t,pos,0,:,:,filter_move] = ndi.shift(imgs_stack[t,pos,0,:,:,filter_move], dxdyt) return imgs_stack @@ -449,7 +530,11 @@ def bacteria_bg(bacteria_contour, image): bg_values = np.empty(len(bcontours)) # Estimation of camera bg to fix a lower bound - camera_bg = np.percentile(np.ma.masked_equal(image, 0), (0.2,99.8))[0] + maimage = np.ma.masked_equal(image, 0) + if len(maimage[~maimage.mask]) > 0: + camera_bg = np.percentile(maimage[~maimage.mask], (0.2,99.8))[0] + else: + camera_bg = 0 for i in range(len(bcontours)): @@ -623,7 +708,7 @@ def segment_bactos(image, image_band=1, denoise=True, -def segment_bactos_stardist(image, image_band=1, model=None, prob_thresh): +def segment_bactos_stardist(image, image_band=1, model=None, prob_thresh=None): """ Use StarDist trained on TELEMOS bactos to segment bacteria on fluorescence images. @@ -644,7 +729,7 @@ def segment_bactos_stardist(image, image_band=1, model=None, prob_thresh): """ - # Select image + # Select image img = image[:,:,image_band] # Normalise image @@ -788,9 +873,9 @@ def compute_bactos_properties(labels, images, position=0, pix_per_micron=1, bactos_bbox = [r.bbox for r in regprops] bactos_labels = [r.label for r in regprops] + el = disk(2) if fluo_bg == 'auto': process_fluo = True - el = disk(2) bacteria_contours = enlarged_bacteria_contour(labels) else: process_fluo = False @@ -798,13 +883,16 @@ def compute_bactos_properties(labels, images, position=0, pix_per_micron=1, bactos_mean_intensity = np.empty((len(images),len(np.unique(labels))-1, len(channels))) for t in range(len(images)): for chan in channels: - regs = regionprops(labels, intensity_image=images[t, position, 0, :, :, chan]) + + # opening filter to remove cosmics! + imgt = opening(images[t, position, 0, :, :, chan], el) + + + regs = regionprops(labels, intensity_image=imgt) for i, r in enumerate(regs): bactos_mean_intensity[t,r.label-1, channels.index(chan)] = r.mean_intensity if process_fluo: - # opening filter to remove cosmics! - imgt = opening(images[t, position, 0, :, :, chan], el) fluo_bg = bacteria_bg(bacteria_contours, imgt) bactos_mean_intensity[t, :, channels.index(chan)] = bactos_mean_intensity[t, :, channels.index(chan)] - fluo_bg @@ -817,7 +905,8 @@ def extract_bactos(imgs_stack, mask_band=1, denoise=True, pix_per_micro=1, use_stardist=False, stardist_modelname = 'stardistBactosTELEMOS', stardist_modeldir = '../stardist/models/', - background=700, channels=[0, 1], prob_thresh=None): + background=700, channels=[0, 1], prob_thresh=None, + time_to_extract_bacteria=0): """ Make mask on selected band for the time=0 and extract intensity @@ -856,6 +945,9 @@ def extract_bactos(imgs_stack, mask_band=1, denoise=True, - prob_thresh, None or a float value btwn 0 1: Set the probability threshold for stadist predict (use None to use the best one from the training) + + - time_to_extract_bacteria, int: + The image to use in the time serie to extract the bacteria, the default is the first one so 0. """ data = {} @@ -871,10 +963,10 @@ def extract_bactos(imgs_stack, mask_band=1, denoise=True, tmpd = {} if use_stardist: - labels = segment_bactos_stardist(imgs_stack[0,p,0], image_band=mask_band, + labels = segment_bactos_stardist(imgs_stack[time_to_extract_bacteria,p,0], image_band=mask_band, model=stardist_model, prob_thresh=prob_thresh) else: - labels = segment_bactos(imgs_stack[0,p,0], image_band=mask_band, denoise=denoise, + labels = segment_bactos(imgs_stack[time_to_extract_bacteria,p,0], image_band=mask_band, denoise=denoise, min_size=bact_size[0], max_size=bact_size[1], min_eccentricity=bact_eccent[0], max_eccentricity=bact_eccent[1], pix_per_micron=pix_per_micro) @@ -1028,7 +1120,10 @@ def save_to_h5(output_file, imgs_stack, bactos_data, img_metadata, gp_name, def run_processing(data_list, nom_sortie, MM_version=2, background='auto', scale_correction=False, drug_filter=0, ctrl_filter=1, - stardist_path = None, bact_size=(0.5, 15), prob_thresh=None): + stardist_path = None, bact_size=(0.5, 15), prob_thresh=None, + time_to_extract_bacteria=0, build_bacteria_figures=True, + bacteria_extract_mask='TRP', drift_algo='cor', register_filter=True, + register_stack=True, register_mask=True, drift_filter=0): """ Loop over the given folder name to extract data and export them inside one hdf file (containing all data from each folder) @@ -1058,6 +1153,34 @@ def run_processing(data_list, nom_sortie, MM_version=2, background='auto', - prob_thresh, None or a float value bitween (0, 1) The probability threshold used by stardist predict (if None, use the best one from the training dataset) + - time_to_extract_bacteria, int: + The image to use in the time serie to extract the bacteria, the default is the first one so 0. + + - build_bacteria_figures, bool: + Add a figure for each bacteria to the hdf5 file in an svg format. This figure contains a summary for each bacteria. + + - bacteria_extract_mask, str in ['TRP' or 'DRUG']: + Give the mask to extract the mask of bacteria. The default 'TRP' is on the control filter + (the one of the tryptophan). If 'DRUG' is given the position of drug_filter is used + + - drift_algo, str in ['cor', 'sift']: + Define the algorithm to compute the drift between images. + - 'cor': use crosscorrelation (the default) + - 'sift': use SIFT point finding + RANSAC to find transform matrix + + - register_filter, bool: + Register the drug_filter on the ctrl_filter (default True) + + - register_stack, bool: + Register the image stack using *drift_algo*, (default True). If false the + raw stack is used. + + - register_mask, bool: + Register the two filter based on bacteria mask (default True) + + - drift_filter, bool: + The filter on which to register the images in time + Example: -------- @@ -1069,7 +1192,9 @@ def run_processing(data_list, nom_sortie, MM_version=2, background='auto', run_processing(expes, out_file) """ - for c in data_list: + #for i in track(range(len(data_list)), description="Processing data..."): + for i in range(len(data_list)): + c = data_list[i] datap = c last_dir = os.path.normpath(c).split(os.path.sep)[-1] print('Process %s, use name %s' % (c, last_dir)) @@ -1109,32 +1234,49 @@ def run_processing(data_list, nom_sortie, MM_version=2, background='auto', print('stabilisation') imgs = register_images(imgs, gaussian_size=0, move_filter=drug_filter, - ref_filter=ctrl_filter, register_filter=True, - show_info=False, up_factor=10, drift_filter=0, - scale_correction=scale_correction) + ref_filter=ctrl_filter, register_filter=register_filter, + register_stack=register_stack, show_info=False, up_factor=10, drift_filter=drift_filter, + scale_correction=scale_correction, + drift_algo=drift_algo) # Extract bacteria data = {} + mask_band = ctrl_filter + filter_ref = ctrl_filter + filter_target = drug_filter + if bacteria_extract_mask == 'DRUG': + mask_band = drug_filter + filter_ref = drug_filter + filter_target = ctrl_filter + bdata = extract_bactos(imgs, bact_size=bact_size, pix_per_micro=pix_per_micro, use_stardist=True, stardist_modeldir=stardist_path, - mask_band=ctrl_filter, + mask_band=mask_band, channels = sorted([drug_filter, ctrl_filter]), - prob_thresh=prob_thresh) + prob_thresh=prob_thresh, + time_to_extract_bacteria=time_to_extract_bacteria, + background=background) # Estimate position between two filters using bacteria pattern matching - for pos in range(imgs.shape[1]): - print('Resgister filter1 and filter0 for pos %i' % pos) - imgs = register_filter_from_bactos(bdata['POS%i'%pos]['bmask'], imgs, position=pos, - filter_ref=ctrl_filter, filter_target=drug_filter) - - # Update intensity data on re-registred images - bdata = update_bactos_data(bdata, imgs, pix_per_micro, background=background, - channels = sorted([drug_filter, ctrl_filter])) + if register_mask: + if ctrl_filter != drug_filter: + for pos in range(imgs.shape[1]): + print('Resgister filter1 and filter0 for pos %i' % pos) + imgs = register_filter_from_bactos(bdata['POS%i'%pos]['bmask'], imgs, position=pos, + filter_ref=filter_ref, filter_target=filter_target) + + # Update intensity data on re-registred images + bdata = update_bactos_data(bdata, imgs, pix_per_micro, background=background, + channels = sorted([drug_filter, ctrl_filter])) if len(bdata['POS0']['labels']) + len(bdata['POS1']['labels']) > 0: # Sauvegarde save_to_h5(nom_sortie, imgs, bdata, metas, last_dir, pix_per_micro, drug_filter=drug_filter, ctrl_filter=ctrl_filter) + + if build_bacteria_figures: + for pos in range(imgs.shape[1]): + save_bacterias_svg_figures(nom_sortie, last_dir, pos) else: print('No bacteria found, no need to save the data') @@ -1239,3 +1381,251 @@ def correct_scaling(img_stack, ref_image, moving_image): img_stack = register_stack_for_scaling(img_stack, moving_image, scaling) return img_stack + + +def get_bacteria_bbox(mask, bacteria_id, region_props=None, area_expend=0.5): + """ + Return the coordinate of the bounding box of bacteria + + Paremeters: + ----------- + + - mask, 2D array: + The mask of bacteria + + - bacteria_id, int: + The bacteria that we want on the mask + + - region_props, None or scikit image regionprops object: + The region props object for the mask could be passed to the function. + If its None, the region_props is evaluated using scikit image + + return xmin, xmax, ymin, ymax + """ + + if region_props is None: + region_props = regionprops(mask) + + r = region_props[bacteria_id] + + # Get crops around a bacteria + bxmin, bymin, bxmax, bymax = r.bbox + bheight = bxmax - bxmin + bwidth = bymax - bymin + bfcy, bfcx = r.centroid + offset = area_expend + maxsize = max((bwidth, bheight)) + new_size = maxsize + maxsize*offset + + cropxmin = int(bfcx - new_size/2) + cropxmax = int(bfcx + new_size/2) + cropymin = int(bfcy - new_size/2) + cropymax = int(bfcy + new_size/2) + + if cropxmin < 0: + cropxmin = 0 + if cropymin < 0: + cropymin = 0 + + if cropxmax > mask.shape[0]: + cropxmax = mask.shape[0] + if cropymax > mask.shape[1]: + cropymax = mask.shape[1] + + return cropxmin, cropxmax, cropymin, cropymax + + +def bacteria_image_wall(imgs, pos, crop, channels=None, percentil=(1, 99)): + """ + Create a unique image with all channels and all time of the bacteria + + Parameters: + ----------- + + - imgs, np array: + array of images in the format [time, pos, z, y, x, channel] + + - pos, int: + the position to use on the imgs stack + + - crop, list of int of size 4: + the spatial limit of the stack defined as [xmin, xmax, ymin, ymax] + + -channels, None or list of channel position: + Position of channel to stack vertically in the image_wall. If None all channel will be used. + + - percentil, tuple of size 2: + the min percentil and max percentil to crop the intensities of imgs + for autoscaling. + + return array + """ + + cropxmin, cropxmax, cropymin, cropymax = crop + + if channels is None: + channels = list(range(imgs.shape[5])) + + imghstack = np.hstack(imgs[:, pos, 0, cropymin:cropymax, cropxmin:cropxmax, channels[0]]) + ps, pe = np.percentile(imghstack, percentil) + imgwall = exposure.rescale_intensity(imghstack, in_range=(ps, pe)) + + if len(channels) > 1: + for ch in channels[1:]: + imghstack = np.hstack(imgs[:,pos,0, cropymin:cropymax, cropxmin:cropxmax, ch]) + ps, pe = np.percentile(imghstack, percentil) + imgwall = np.vstack((imgwall, exposure.rescale_intensity(imghstack, in_range=(ps, pe)))) + + return imgwall + + +def bacteria_summary_graph(bacteria_ids, h5data, experiment, position, + bacteria_offset=0.5, export_as_svg=True): + """ + Build a graphic thats shows pictures time series of a given bacteria + and the evolution of it's fluorescence. + + Parameters: + ----------- + + - bacteria_ids, (int or float) or 'all' or a list of int or float: + The bacteria id of list of bacteria ids to use. If 'all' take the unique + value of the mask to get all bacterias. + + - h5data, str: + The path to the data in hdf format + + - experiment, str: + The name of the experiment. Must be contained in the h5data + + - position, int: + The position, if multi stage position has been used on the microscope, to use + + - bacteria_offset, float: + The size offset (in percent) to add around the bacteria bounding box + + - export_as_svg, bool: + Export the figure as svg format if True (default). If False return a list of + matplotlib figure objects. + + return a list of figures + """ + + + if isinstance(bacteria_ids, (int, float)): + bacteria_ids = [bacteria_ids] + + h5data = Path(h5data) + + # Load the data + pos = position + exp = experiment + imgs = load_image(h5data, exp) + masks = load_mask(h5data, exp, pos) + bactosct = bacteria_contour(masks) + Ib_chan0 = load_intensity_profile(h5data, exp, 0, pos) + Ib_chan1 = load_intensity_profile(h5data, exp, 1, pos) + + if isinstance(bacteria_ids, str) and bacteria_ids == 'all': + bacteria_ids = np.unique(masks)[1:]-1 + + # Compute the properties of bacteria masks + rprops = regionprops(masks) + + # This list will contain all figures (one for each bacteria) + figures = [] + + for bid in bacteria_ids: + r = rprops[bid] + + # Get crops around a bacteria + cropxmin, cropxmax, cropymin, cropymax = get_bacteria_bbox(masks, bid, rprops, bacteria_offset) + # Compute the new width and height for this bacteria + nwidth = cropxmax - cropxmin + nheight = cropymax - cropymin + + # Get the centroid of bacteria + bfcy, bfcx = r.centroid + + # Start to build the figure + fig = plt.figure(figsize=(13, 1.5)) + plt.title(f'Exp: {exp}, Pos: {pos}, bacteria-{bid+1} (in file {h5data.name})') + G = plt.GridSpec(1, 9) + + # The main image to locate the bacteria on the FOV + plt.subplot(G[0, 0]) + p2, p98 = np.percentile(imgs[0,pos,0, ::2, ::2, 0], (2, 98)) + plt.imshow(exposure.rescale_intensity(imgs[0,pos,0, ::2, ::2, 0], in_range=(p2, p98)), 'gray') + plt.plot(bfcx/2, bfcy/2, 'ro', fillstyle='none') + plt.axis('off') + plt.text(-8, imgs.shape[4]/4, f'Bacteria-{bid+1}', rotation=90, ha='right', va='center') + plt.text(imgs.shape[4]/4, 0, f'POS - {pos}', ha='center', va='bottom') + + # Create the timewall of all images of this bacteria through time and filters + plt.subplot(G[0, 1:4]) + imgwall = bacteria_image_wall(imgs, pos, (cropxmin, cropxmax, cropymin, cropymax)) + plt.imshow(imgwall, 'gray') + + # Add channel names + for ch in range(imgs.shape[5]): + plt.text(-5, nheight/2+ch*nheight, f'C-{ch}', ha='right', va='center') + + # Plot the contour of bacteria on first row of the timewall and for all time steps + for j in range(len(imgs)): + plt.plot((bactosct[bid][1] + j*nwidth)-cropxmin, bactosct[bid][0]-cropymin, 'C3--') + + plt.axis('off') + + # Plot the mean intensity of each filters inside the segmented bacteria + plt.subplot(G[0, 4:]) + plt.plot(Ib_chan1[:, bid], 'C3o-') + ax = plt.gca() + ax.spines["right"].set_visible(False) + ax.spines["top"].set_visible(False) + ax.set_ylabel('I-C-0', color='C3') + ax.tick_params(axis='y', labelcolor='C3') + + ax2 = ax.twinx() + ax2.plot(Ib_chan0[:, bid], 'C0o-') + ax2.spines["left"].set_visible(False) + ax2.spines["top"].set_visible(False) + ax2.tick_params(axis='y', labelcolor='C0') + ax2.set_ylabel('I-C-1', color='C0') + + + plt.tight_layout(pad=0.9) + + if export_as_svg: + fsvg = io.BytesIO() + plt.savefig(fsvg, format='svg') + figures += [fsvg.getvalue().strip().replace(b'\n', b'')] + plt.close(fig) + del fig + else: + figures += [fig] + + return figures + + +def save_bacterias_svg_figures(h5file, exp, pos): + """ + Save figures in svg format for each bacteria, the figure shows a closed view of the + bacteria it's location and the signal inside the mask with time. + """ + + # create all the figures + print('Create svg figures in Parralel') + nbactos = np.unique(load_mask(h5file, exp, pos))[1:]-1 + figs = Parallel(n_jobs=-1, backend='multiprocessing')(delayed(bacteria_summary_graph)(int(b), h5file, exp, pos) for b in nbactos) + figs = [f[0] for f in figs] + + # Not in parralel we can use the following function + # figs = bacteria_summary_graph('all', h5file, exp, pos) + + print(f'Save svg figures to {exp}/{pos}/bacteria_svg_figures in {h5file}') + with h5py.File(h5file, 'a') as f: + pos = 'POS%i' % pos + if 'bacteria_svg_figures' in f[exp][pos]: + del f[exp][pos]['bacteria_svg_figures'] + + f[exp][pos].create_dataset('bacteria_svg_figures', data=figs, compression="gzip") diff --git a/src/bactos_viewver.py b/src/bactos_viewver.py index 8c7b67f03894e900f5b0530e721bcc26224214a5..faa693db9f33ab07c63fefa66d5de02e1a88ab37 100755 --- a/src/bactos_viewver.py +++ b/src/bactos_viewver.py @@ -3,7 +3,7 @@ Use ipywidgets and ipympl to create a graphical user interface inside jupyter notebook to check bacteria drug uptake on DISCO UV fluorescence image stack. """ -import ipywidgets as widgets +from ipywidgets import widgets, GridspecLayout, Layout import numpy as np import matplotlib.pyplot as plt from matplotlib.collections import PatchCollection @@ -18,13 +18,14 @@ import os import shutil import logging from matplotlib.backend_bases import NavigationToolbar2, Event +from base64 import b64encode OLDHOME = NavigationToolbar2.home def new_home(self, *args, **kwargs): # register a new event OLDHOME(self, *args, **kwargs) - event = Event('home_event', self) - self.canvas.callbacks.process('home_event', event) + # event = Event('home_event', self) + # self.canvas.callbacks.process('home_event', event) # redefine home command NavigationToolbar2.home = new_home @@ -73,7 +74,7 @@ def get_all_pos_fluo(selected_file, selected_exp, limit_to_selected_bactos=True, I0 = I0[:, selected_bactos] I1 = I1[:, selected_bactos] - # LE TRP est sur le channel 1 et le filtre 2 sur le channel 0 + # LE TRP est sur le channel 1 et le filtre 2 sur le channel 0 IF1 += [I1] IF2 += [I0] @@ -459,6 +460,8 @@ class BactosViewver(): self.Wimgout = widgets.Output() self.WplotAll = widgets.Output() self.Woutbactos = widgets.Output(layout=widgets.Layout(width="330px")) + self.BacteriaPlots = widgets.Output(layout=Layout(min_height="600px")) + self.bacteria_grid = None self.BactosAreaOutput = widgets.Output() self.I1output = widgets.Output() self.I2output = widgets.Output() @@ -535,9 +538,9 @@ class BactosViewver(): # Create the tab pannels tabs = widgets.Tab() # layout=widgets.Layout(width="75%") - tabs.children = [wimgs, self.AllPosPlotOutput, self.AllExpPlotOutput, self.BactosAreaOutput, + tabs.children = [wimgs, self.AllPosPlotOutput, self.AllExpPlotOutput, self.BacteriaPlots, self.BactosAreaOutput, self.I1output, self.I2output, self.I1overI2output, self.DrugIoutput] - tabs_titles = ('Plot: data explorer', 'Plot: Idrug for all pos', 'Plot: all exp', 'Data: Bactos area', + tabs_titles = ('Plot: data explorer', 'Plot: Idrug for all pos', 'Plot: all exp', 'Plot: bacteria summary','Data: Bactos area', 'Data: F2 fluo', 'Data: F1 fluo', 'Data: F2/F1 control', 'Data: Drug fluo') for i, title in enumerate(tabs_titles): @@ -560,6 +563,7 @@ class BactosViewver(): footer=None) tabs.observe(self.on_tab_changed, names='selected_index') + self.tabs = tabs display(self.final_lay) @@ -646,10 +650,18 @@ class BactosViewver(): # Ajust min-max if adjust_minmax: - self.F2minmax.value = list(np.percentile(np.ma.masked_equal(self.selected_imgs[0, 0, 0, :, :, self.drug_filter], 0), - (0.2, 99.9))) - self.F1minmax.value = list(np.percentile(np.ma.masked_equal(self.selected_imgs[0, 0, 0, :, :, self.ctrl_filter], 0), + try: + self.F2minmax.value = list(np.percentile(np.ma.masked_equal(self.selected_imgs[0, 0, 0, :, :, self.drug_filter], 0), + (0.2, 99.9))) + except: + pass + + try: + self.F1minmax.value = list(np.percentile(np.ma.masked_equal(self.selected_imgs[0, 0, 0, :, :, self.ctrl_filter], 0), (0.2, 99.9))) + except: + pass + else: self.plimg1.set_clim(self.F2minmax.value) self.plimg2.set_clim(self.F1minmax.value) @@ -665,7 +677,8 @@ class BactosViewver(): # Connect the update of profiles based on bacteria currently displayed on the image self.figure.canvas.mpl_connect('button_release_event', self.update_bactos_inview) - self.figure.canvas.mpl_connect('home_event', self.update_bactos_inview) + # TODO: MPL find new way to triger home event (it's no more available) + # self.figure.canvas.mpl_connect('home_event', self.update_bactos_inview) with self.WplotAll: plt.tight_layout() @@ -676,7 +689,7 @@ class BactosViewver(): drawing each contour individually). """ - # get figure axes + # get figure axes ax1, ax2, ax3, ax4, ax5, ax6 = self.figure.get_axes() if force_redraw: @@ -1106,7 +1119,7 @@ class BactosViewver(): self.figure_all_exp = plt.figure('Plot all experiments', figsize=(11,5)) plt.ion() - # Calcul du nombre d'expérience avec un contrôle + # Calcul du nombre d'expérience avec un contrôle # TODO !!! ax1 = plt.subplot(111) @@ -1259,6 +1272,16 @@ class BactosViewver(): sel = len([i for i, val in enumerate(self.selected_bactos) if val]) self.accordion.set_title(1, 'Select bactos (%i/%i)' % (sel, tot) ) + # Update the bacteria plot summary border if needed + if self.bacteria_grid is not None: + cb = self.bactos_colors[id_bactos] * 255 + if new_state: + bstyle = f"2px solid rgb({cb[0]}, {cb[1]}, {cb[2]})" + else: + bstyle = "2px dashed gray" + + self.bacteria_grid[id_bactos, 0].layout.border = bstyle + #Update tables self.update_area_table() self.update_intensity_tables() @@ -1323,7 +1346,7 @@ class BactosViewver(): self.H = len(self.selected_imgs[0,0,0,0,:,0]) self.bactos_dist_slider.max = min(self.W, self.H)/2 - 10 - # Need to update the extent + # Need to update the extent self.EXTENT = (0, self.H, self.W, 0) self.plimg1.set_extent(self.EXTENT) self.plimg2.set_extent(self.EXTENT) @@ -1350,7 +1373,7 @@ class BactosViewver(): self.ListPos.options = self.positions self.ListPos.value = self.selected_pos - # Load the control experiment if it as been given + # Load the control experiment if it as been given ctrl_exp = load_control_exp(self.selected_file, selected) # print(ctrl_exp) if ctrl_exp is not None: @@ -1367,7 +1390,7 @@ class BactosViewver(): # update bactos contours self.extract_bactos_contour_coordinate() - # print(selected_file, selected, selected_pos) + # print(selected_file, selected, selected_pos) self.i1 = load_intensity_profile(self.selected_file, selected, 0, self.selected_pos) self.i2 = load_intensity_profile(self.selected_file, selected, 1, self.selected_pos) @@ -1380,32 +1403,34 @@ class BactosViewver(): self.update_bactos_area() self.update_bactos_area_minmax() - try: - self.display_stack(self.selected_time, self.selected_pos, 0, 0, True) - except Exception as e: - print('Error to display images') - print(e) + if self.tabs.selected_index == 0: + try: + self.display_stack(self.selected_time, self.selected_pos, 0, 0, True) + except Exception as e: + print('Error to display images') + print(e) - try: - self.display_bactos_contour_fast(force_redraw=True) - except Exception as e: - print('Error to display bactos contour') - print(e) - - if sum(self.selected_bactos): try: - self.plot_profiles() + self.display_bactos_contour_fast(force_redraw=True) except Exception as e: - print('Error in plot fluo profiles') + print('Error to display bactos contour') print(e) + + if sum(self.selected_bactos): + try: + self.plot_profiles() + except Exception as e: + print('Error in plot fluo profiles') + print(e) - with self.AllPosPlotOutput: - try: - self.update_all_pos_plot() - except Exception as e: - print('Error') - print(e) + if self.tabs.selected_index == 1: + with self.AllPosPlotOutput: + try: + self.update_all_pos_plot() + except Exception as e: + print('Error') + print(e) #Update tables @@ -1413,8 +1438,12 @@ class BactosViewver(): self.update_intensity_tables() self.update_control_ratio_table() - self.Wimgout.clear_output() + if self.tabs.selected_index == 3: + self.plot_bacteria_svg_figures() + # Message container need to bee cleaned + self.Wimgout.clear_output() + def on_pos_changed(self, evt): """ Action triggered when the user select a position of the current experiement. @@ -1443,10 +1472,15 @@ class BactosViewver(): self.update_bactos_area_minmax() # Display new data - self.display_stack(self.selected_time, self.selected_pos, 0, 0, True) - self.display_bactos_contour_fast(force_redraw=True) - self.plot_profiles() + if self.tabs.selected_index == 0: + self.display_stack(self.selected_time, self.selected_pos, 0, 0, True) + self.display_bactos_contour_fast(force_redraw=True) + self.plot_profiles() + # Update single bacteria figure + if self.tabs.selected_index == 3: + self.plot_bacteria_svg_figures() + #Update tables self.update_area_table() self.update_intensity_tables() @@ -1458,7 +1492,9 @@ class BactosViewver(): new_ctrl_exp = evt['new'] save_control_exp(self.selected_file, self.ListExp.value, new_ctrl_exp) - self.plot_profiles() + + if self.tabs.selected_index == 0: + self.plot_profiles() # update the table self.update_control_ratio_table() @@ -1594,6 +1630,37 @@ class BactosViewver(): for a in (ax4, ax5, ax6): ax.set_xlim(0, self.selected_imgs.shape[0]) + def plot_bacteria_svg_figures(self): + """ + Update the plot of each bacteria summary + """ + with self.BacteriaPlots: + self.BacteriaPlots.clear_output() + try: + figs = load_bacterias_svg_figures(self.selected_file, self.ListExp.value, self.selected_pos) + except Exception as e: + print('error', e) + figs = None + self.bacteria_grid = None + + if figs is not None: + # figs = bacteria_summary_graph('all', fpath, exp, pos) + self.bacteria_grid = GridspecLayout(len(figs), 1, layout=Layout(height="800px")) + for i, fig in enumerate(figs): + dd = b64encode(fig) + url = f'data:image/svg+xml;base64,{dd.decode("utf8")}' + img_url = f'<img src="{url}" style="width: 100%; height: 100%;" >' + self.bacteria_grid[i, 0] = widgets.HTML(value=img_url) + cb = self.bactos_colors[i] * 255 + if self.selected_bactos[i]: + bstyle = f"2px solid rgb({cb[0]}, {cb[1]}, {cb[2]})" + else: + bstyle = "2px dashed gray" + + self.bacteria_grid[i, 0].layout.border = bstyle + + display(self.bacteria_grid) + def on_time_changed(self, evt): """ update the time varible @@ -1655,6 +1722,9 @@ class BactosViewver(): self.AllExpPlotOutput.clear_output(wait=True) self.update_all_exp_plot() plt.show(self.figure_all_exp) + + if tabid == 3: + self.plot_bacteria_svg_figures() def on_export(self, btn): """ @@ -1672,7 +1742,7 @@ class BactosViewver(): # Output file name (should be the same as the imput source) fout = os.path.abspath(self.source_file) fout = os.path.splitext(fout)[0] + '.xlsx' - + prevsname = [] with p.ExcelWriter(fout, engine='xlsxwriter') as writer: for expe in self.ListExp.options: ctrl_exp = load_control_exp(self.selected_file, expe) diff --git a/src/lib_bactos_hdf.py b/src/lib_bactos_hdf.py index 22caca9044a430d487f84ccb1d9447f51e5a0f00..1b40c844c5929bf0f09825baa523858be6982111 100755 --- a/src/lib_bactos_hdf.py +++ b/src/lib_bactos_hdf.py @@ -47,7 +47,7 @@ def load_image(data_in, expe_name): with h5py.File(data_in, 'r') as f: output = f[expe_name]['images'][:] - + return output @@ -74,7 +74,7 @@ def load_mask(data_in, expe_name, pos): with h5py.File(data_in, 'r') as f: pos = 'POS%i' % pos output = f[expe_name][pos]['bacteria_mask'][:] - + return output def get_stack_shape(data_in, expe_name): @@ -94,7 +94,7 @@ def get_stack_shape(data_in, expe_name): """ with h5py.File(data_in, 'r') as f: output = f[expe_name]['images'].shape - + return output def load_intensity_profile(data_in, expe_name, channel, pos): @@ -123,7 +123,7 @@ def load_intensity_profile(data_in, expe_name, channel, pos): with h5py.File(data_in, 'r') as f: pos = 'POS%i' % pos output = f[expe_name][pos]['ch%i_bacteria_intensities' % channel][:] - + return output @@ -149,8 +149,8 @@ def load_area(data_in, expe_name, pos): with h5py.File(data_in, 'r') as f: pos = 'POS%i' % pos output = f[expe_name][pos]['bacteria_area'][:] - - return output + + return output def save_selected_bactos(data_in, expe_name, pos, selected_bactos): @@ -205,7 +205,7 @@ def save_drug_signal(data_in, expe_name, pos, drug_data): pos = 'POS%i' % pos if 'drug_signal' in f[expe_name][pos]: del f[expe_name][pos]['drug_signal'] - + f[expe_name][pos].create_dataset('drug_signal', data=drug_data) @@ -228,16 +228,16 @@ def load_drug_signal(data_in, expe_name, pos): """ output = None - + with h5py.File(data_in, 'r') as f: pos = 'POS%i' % pos - + if 'drug_signal' in f[expe_name][pos]: output = f[expe_name][pos]['drug_signal'][:] - + return output - - + + def load_selected_bactos(data_in, expe_name, pos): """ Get the list of user selected bacteria save in h5file for a given experiment @@ -264,7 +264,7 @@ def load_selected_bactos(data_in, expe_name, pos): output = f[expe_name][pos]['selected_bactos'][:] else: output = None - + return output @@ -339,18 +339,17 @@ def load_filter_positions(data_in, expe_name): print('Use pos 0 for drug and pos 1 for controle') return ctrl_filter, drug_filter - - + def load_bacterias_svg_figures(h5file, exp, pos): """ Load svg figures of bacterias """ - + with h5py.File(h5file, 'r') as f: pos = 'POS%i' % pos if 'bacteria_svg_figures' in f[exp][pos]: output = f[exp][pos]['bacteria_svg_figures'][:] else: output = None - + return output