diff --git a/setup.cfg b/setup.cfg index 4aff539ee4f8b45e96c59fc9c13139e170f7cbba..acaf86953eb94d0ce55fc7bef3f3ecc4c69cbcde 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = epicure -version = 0.2.5 +version = 0.2.6 description = Napari plugin to manually correct epithelia segmentation in movies long_description = file: README.md long_description_content_type = text/markdown @@ -35,12 +35,15 @@ install_requires = numpy magicgui qtpy + pyqtwebengine scikit-image + scipy opencv_python_headless roifile xlsxwriter laptrack - matplotlib + plotly + kaleido imagecodecs edt diff --git a/src/epicure/Utils.py b/src/epicure/Utils.py index 597725aff9ba0bbae5a7cf92b5ab79e6f25498f2..3c92529c67538fa9af611e5ed55b5d30f71747df 100644 --- a/src/epicure/Utils.py +++ b/src/epicure/Utils.py @@ -4,12 +4,18 @@ import time import math from skimage.measure import regionprops, find_contours, regionprops_table from skimage.segmentation import find_boundaries, expand_labels -from napari.utils.translations import trans -from napari.utils.notifications import show_info -from napari.utils import notifications as nt +from napari.utils.translations import trans # type: ignore +from napari.utils.notifications import show_info # type: ignore +from napari.utils import notifications as nt # type: ignore from skimage.morphology import skeletonize, binary_dilation, disk +from scipy.ndimage import center_of_mass, find_objects +from scipy.ndimage import label as ndlabel +from scipy.ndimage import sum as ndsum +from scipy.ndimage import generate_binary_structure as ndi_structure import pandas as pd from epicure.laptrack_centroids import LaptrackCentroids +import tifffile as tif # type: ignore +from napari.utils import progress # type: ignore try: from skimage.graph import RAG @@ -44,6 +50,10 @@ def show_documentation_page(page): webbrowser.open_new_tab("https://gitlab.pasteur.fr/gletort/epicure/-/wikis/"+page) return +def show_progress( viewer, show ): + """ Show.hide the napari activity bar to see processing progress """ + viewer.window._status_bar._toggle_activity_dock( show ) + def get_directory(imagepath): return os.path.dirname(imagepath) @@ -70,6 +80,10 @@ def suggest_segfile(out, imgname): return segfile else: return None + +def found_segfile( filepath ): + """ Check if the segmentation file exists """ + return os.path.exists( filepath ) def get_filename(outdir, imgname): return os.path.join( outdir, imgname ) @@ -162,7 +176,6 @@ def remove_widget(viewer, widname): wid.destroyOnClose() def opentif(imagepath, verbose=True): - import tifffile as tif img = tif.TiffFile(imagepath) metadata = img.imagej_metadata #print(metadata) @@ -188,7 +201,7 @@ def opentif(imagepath, verbose=True): scalet = float(metadata['finterval']) if 'unit' in metadata: if metadata['unit'] is not None: - unitxy = metatdata['unit'] + unitxy = metadata['unit'] except: metadatas = None #print(info) @@ -202,17 +215,19 @@ def opentif(imagepath, verbose=True): return image, nchan, scale, unitxy, scalet, unitt def writeTif(img, imgname, scale, imtype, what=""): - import tifffile + """ Write image in tif format """ if len(img.shape) == 2: - tifffile.imwrite(imgname, np.array(img, dtype=imtype), imagej=True, resolution=[1./scale, 1./scale], metadata={'unit': 'um', 'axes': 'YX'}) + tif.imwrite(imgname, np.array(img, dtype=imtype), imagej=True, resolution=[1./scale, 1./scale], metadata={'unit': 'um', 'axes': 'YX'}) else: - tifffile.imwrite(imgname, np.array(img, dtype=imtype), imagej=True, resolution=[1./scale, 1./scale], metadata={'unit': 'um', 'axes': 'TYX'}, compression="zstd") + try: + tif.imwrite(imgname, np.array(img, dtype=imtype), imagej=True, resolution=[1./scale, 1./scale], metadata={'unit': 'um', 'axes': 'TYX'}, compression="zstd") + except: + tif.imwrite(imgname, np.array(img, dtype=imtype), imagej=True, resolution=[1./scale, 1./scale], metadata={'unit': 'um', 'axes': 'TYX'}) show_info(what+" saved in "+imgname) def appendToTif(img, imgname): """ Append to RGB tif the current image """ - import tifffile - tifffile.imwrite(imgname, img, photometric="rgb", append=True) + tif.imwrite(imgname, img, photometric="rgb", append=True) def getCellValue(label_layer, event): vis = label_layer.visible @@ -247,15 +262,19 @@ def setCellValue(layer, label_layer, event, newvalue, layer_frame=None, label_fr layer.refresh() return label -def get_skeleton( seg, verbose=0 ) : +def get_skeleton( seg, viewer=None, verbose=0 ) : """ convert labels movie to skeleton (thin boundaries) """ startt = start_time() + if viewer is not None: + show_progress( viewer, show=True ) skel = np.zeros(seg.shape, dtype="uint8") - for z in range(seg.shape[0]): + for z in progress(range(seg.shape[0])): expz = expand_labels( seg[z], distance=1 ) skel[z][(seg[z] == 0) *(expz > 0)] = 1 if verbose > 0: show_duration(startt, header="Skeleton calculted in ") + if viewer is not None: + show_progress( viewer, show=False ) return skel @@ -294,14 +313,16 @@ def changeLabel( label_layer, old_value, new_value ): index = np.argwhere( label_layer.data==old_value ).tolist() setNewLabel( label_layer, index, new_value ) -def setNewLabel(label_layer, indices, newvalue, add_frame=None): +def setNewLabel(label_layer, indices, newvalue, add_frame=None, return_old=True): """ Change the label of all the pixels indicated by indices """ indexs = np.array(indices).T if add_frame is not None: indexs = np.vstack((np.repeat(add_frame, indexs.shape[1]), indexs)) changed_indices = label_layer.data[tuple(indexs)] != newvalue inds = tuple(x[changed_indices] for x in indexs) - oldvalues = label_layer.data[inds] + oldvalues = None + if return_old: + oldvalues = label_layer.data[inds] if isinstance(newvalue, list): newvalue = np.array(newvalue)[np.where(changed_indices)[0]] label_layer.data_setitem( inds, newvalue ) @@ -394,8 +415,9 @@ def getBBoxFromPts(pts, extend, imshape, outdim=None, frame=None): bbox = [int(np.min(arr[:,0])), int(np.min(arr[:,1])), int(np.min(arr[:,2])), int(np.max(arr[:,0]))+1, int(np.max(arr[:,1]))+1, int(np.max(arr[:,2]))+1] if extend > 0: for i in range(outdim): - bbox[(outdim==3)+i] = max( bbox[(outdim==3)+i] - extend, 0) - bbox[(outdim==3)+i+outdim] = min(bbox[(outdim==3)+i+outdim] + extend, imshape[(outdim==3)+i] ) + if i < 2: + bbox[(outdim==3)+i] = max( bbox[(outdim==3)+i] - extend, 0) + bbox[(outdim==3)+i+outdim] = min(bbox[(outdim==3)+i+outdim] + extend, imshape[(outdim==3)+i] ) return bbox def inside_bounds( pt, imshape ): @@ -489,7 +511,16 @@ def addFrameIndices( indices, frame ): return [ [frame, ind[0], ind[1]] for ind in indices ] def shiftFrameIndices( indices, add_frame ): - return [ [add_frame, ind[0], ind[1]] for ind in indices ] + if isinstance( indices, list ): + indices = np.array(indices) + indices[:, 0] += add_frame + return indices.tolist() + +def shiftFrames( indices, frames ): + if isinstance( indices, list ): + indices = np.array(indices) + indices[:, 0] = frames[indices[:, 0]] + return indices.tolist() def toFullMoviePos( indices, bbox, frame=None ): """ Replace indexes inside bounding box to full movie indexes """ @@ -497,6 +528,8 @@ def toFullMoviePos( indices, bbox, frame=None ): if frame is not None: frame_arr = np.full(len(indices), frame) return np.column_stack((frame_arr, indices[:, 0] + bbox[0], indices[:, 1] + bbox[1])) + if len(bbox) == 6: + return np.column_stack((indices[:, 0] + bbox[0], indices[:, 1] + bbox[1], indices[:, 2] + bbox[2])) return np.column_stack((indices[:, 0], indices[:, 1] + bbox[0], indices[:, 2] + bbox[1])) def cropBBox(img, bbox): @@ -567,6 +600,20 @@ def inv_visibility(viewer, layername): layer.visible = not layer.visible ######## Measure labels +def average_area( seg ): + """ Average area of labels (cells) """ + # Label the input image + labeled_array, num_features = ndlabel(seg) + + if num_features == 0: + return 0.0 + + # Calculate the area of each label + areas = ndsum(seg > 0, labeled_array, index=np.arange(1, num_features + 1)) + # Calculate the average area + avg_area = np.mean(areas) + return avg_area + def summary_labels( seg ): """ Summary of labels (cells) measurements """ @@ -624,6 +671,49 @@ def labels_table( labimg, intensity_image=None, properties=None, extra_propertie return regionprops_table( labimg, intensity_image=intensity_image, properties=properties, extra_properties=extra_properties ) return regionprops_table( labimg, properties=properties, extra_properties=extra_properties ) +def labels_to_table( labimg, frame ): + """ Get label and centroid """ + labels = np.unique(labimg) + labels = labels[labels != 0] + centroids = center_of_mass(labimg, labels=labimg, index=labels) + table = np.column_stack((labels, np.full(len(labels), frame), centroids)) + return table.astype(int) + +def non_unique_labels( labimg ): + """ Check if contains only unique labels """ + relab, nlabels = ndlabel( labimg ) + return nlabels > (len( np.unique(labimg) )-1) + +def reset_labels( labimg ): + """ Relabel in 3D all labels (unique labels) """ + s = ndi_structure(3,1) + ## ignore 3D connectivity (unique labels in all frames) + s[0,:,:] = 0 + s[2,:,:] = 0 + lab = ndlabel( labimg, structure=s )[0] + return lab + + +def skeleton_to_label( skel, labelled ): + """ Transform a skeleton to label image with numbers from labelled image """ + labels = ndlabel( np.invert(skel) )[0] + new_labels = find_objects( labels ) + newlab = np.zeros( skel.shape, np.uint32 ) + for i, obj_slice in enumerate(new_labels): + if obj_slice is not None: + label_mask = labels[obj_slice] == (i+1) + label_values = labelled[obj_slice][label_mask] + labvals, counts = np.unique(label_values, return_counts=True ) + labval = labvals[ np.argmax(counts) ] + newlab[obj_slice][label_mask] = labval + return newlab + +def get_most_frequent( labimg, img, label ): + """ Returns which label is the most frequent in mask """ + mask = labimg == label + vals, counts = np.unique( img[mask], return_counts=True ) + return vals[ np.argmax(counts) ] + def labels_properties( labimg ): """ Returns basic label properties """ return regionprops( labimg ) @@ -702,15 +792,73 @@ def connectivity_graph( img, distance ): """ Returns the region adjancy graph of labels """ touchlab = touching_labels( img, expand=distance ) return RAG( touchlab, connectivity=2 ) + +def get_boundary_cells( img ): + """ Return cells on tissu boundary in current image """ + dilated = binary_dilation( img > 0, disk(3) ) + zero = np.invert( dilated ) + zero = binary_dilation( zero, disk(5) ) + touching = np.unique( img[ zero ] ).tolist() + if 0 in touching: + touching.remove(0) + return touching + +def get_border_cells( img ): + """ Return cells on border in current image """ + height = img.shape[1] + width = img.shape[0] + labels = list( np.unique( img[ :, 0:2 ] ) ) ## top border + labels += list( np.unique( img[ :, (height-2): ] ) ) ## bottom border + labels += list( np.unique( img[ 0:2,] ) ) ## left border + labels += list( np.unique( img[ (width-2):,] ) ) ## right border + return labels + +def count_neighbors( label_img, label ): + """ Get the number of neighboring labels of given label """ + ## much slower than using the RAG graph + # Dilate the labeled image + dilated_mask = binary_dilation( label_img==label, disk(1) ) + nonzero = np.nonzero( dilated_mask) + + # Find the unique labels in the dilated region, excluding the current label and background + neighboring_labels = np.unique( label_img[nonzero] ).tolist() + + # Add the number of unique neighboring labels + return len(neighboring_labels) - 1 - 1*(0 in neighboring_labels) ## don't count itself or 0 + +def get_cell_radius( label, labimg ): + """ Get the radius of the cell label in labimg (2D) """ + area = np.sum( labimg == label ) + return math.sqrt( area / math.pi ) + ####### Distance measures +def consecutive_distances( pts_pos ): + """ Distance travelled by the cell between each frame """ + diff = np.diff( pts_pos, axis=0 ) + disp = np.linalg.norm(diff, axis=1) + return disp + +def velocities( pts_pos ): + """ Velocity of the cell between each frame (average between previous and next) """ + diff = np.diff( pts_pos, axis=0 ).astype(float) + diff = np.vstack( (diff[0], diff) ) + diff = np.vstack( (diff, diff[-1]) ) + kernel=np.array([0.5,0.5]) + adiff = np.zeros( (len(diff)+1, 3) ) + for i in range(3): + adiff[:,i] = np.convolve( diff[:,i], kernel ) + adiff = adiff[1:-1] + disp = np.linalg.norm(adiff[:,1:3], axis=1) + dt = adiff[:,0] + return disp/dt + def total_distance( pts_pos ): """ Total distance travelled by point with coordinates xpos and ypos """ - after = pts_pos[1:] - before = pts_pos[:(-1)] - disp = after - before - return np.sum( np.sqrt( np.square(disp[:,0]) + np.square(disp[:,1]) ) ) + diff = np.diff( pts_pos, axis=0 ) + disp = np.linalg.norm(diff, axis=1) + return np.sum(disp) def net_distance( pts_pos ): """ Net distance travelled by point with coordinates xpos and ypos """ diff --git a/src/epicure/displaying.py b/src/epicure/displaying.py index 76ecc40726a3dc20cfbcd2061d3100c560e253df..c4c932da28c48e75556672cd4e4280d77402f77e 100644 --- a/src/epicure/displaying.py +++ b/src/epicure/displaying.py @@ -347,7 +347,7 @@ class Displaying(QWidget): rect = np.array([[x*wid, y*hei], [(x+1)*wid, (y+1)*hei]]) rects.append(rect) rects_names.append(chr(65+x)+"_"+str(y)) - self.viewer.add_shapes(rects, name="EpicGrid", text=rects_names, face_color=[1,0,0,0], edge_color=self.grid_color, edge_width=gwidth, opacity=0.7) + self.viewer.add_shapes(rects, name="EpicGrid", text=rects_names, face_color=[1,0,0,0], edge_color=self.grid_color, edge_width=gwidth, opacity=0.7, scale=self.viewer.layers["Segmentation"].scale[1:]) self.viewer.layers["EpicGrid"].text.visible = False ut.set_active_layer( self.viewer, "Segmentation" ) diff --git a/src/epicure/editing.py b/src/epicure/editing.py index 977b472d5de7991850b7fc90c3dc62c4725f22bf..42b75a739e58f1770b3997c4842a4cd93d9df115 100644 --- a/src/epicure/editing.py +++ b/src/epicure/editing.py @@ -1,17 +1,17 @@ import numpy as np -import edt +import edt # type: ignore from skimage.segmentation import watershed, clear_border, find_boundaries, random_walker from skimage.measure import label, points_in_poly from skimage.morphology import binary_closing, binary_opening, binary_dilation, binary_erosion, disk -from qtpy.QtWidgets import QVBoxLayout, QWidget -from napari.layers.labels._labels_utils import interpolate_coordinates +from qtpy.QtWidgets import QVBoxLayout, QWidget # type: ignore from scipy.ndimage import binary_fill_holes, distance_transform_edt, generate_binary_structure from scipy.ndimage import label as ndlabel -from napari.layers.labels._labels_utils import sphere_indices -from napari.utils import progress +from napari.layers.labels._labels_utils import sphere_indices # type: ignore +from napari.layers.labels._labels_utils import interpolate_coordinates # type: ignore +from napari.utils import progress # type: ignore +from napari.qt.threading import thread_worker # type: ignore import epicure.Utils as ut import epicure.epiwidgets as wid -from napari.qt.threading import thread_worker class Editing( QWidget ): """ Handle user interaction to edit the segmentation """ @@ -28,6 +28,10 @@ class Editing( QWidget ): layout = QVBoxLayout() + ## Option to use default napari painting options + #self.napari_painting = wid.add_check( "Default Napari painting tools (no checks)", checked=False, check_func=self.painting_tools, descr="Use the label painting of Napari instead of customized EpiCure ones (will not perform any sanity check)" ) + #layout.addWidget( self.napari_painting ) + ## Option to remove all border cells clean_line, self.clean_vis, self.gCleaned = wid.checkgroup_help( name="Cleaning options", checked=False, descr="Show/hide options to clean the segmentation", help_link="Edit#cleaning-options", display_settings=self.epicure.display_colors, groupnb="group" ) layout.addLayout(clean_line) @@ -75,6 +79,15 @@ class Editing( QWidget ): self.paint_radius() self.disk_one = disk(radius=1) + def painting_tools( self ): + """ Choose which painting tools should be activated """ + if self.napari_painting.isChecked(): + self.epicure.seglayer.fill = self.napari_fill + self.epicure.seglayer.paint = self.napari_paint + else: + self.epicure.seglayer.fill = self.epicure_fill + self.epicure.seglayer.paint = self.lazy + def apply_settings( self, settings ): """ Load the prefered settings for Edit panel """ @@ -318,34 +331,25 @@ class Editing( QWidget ): parents = [None]*len(labels) if tframe > 0: twoframes = ut.crop_twoframes( self.epicure.seglayer.data, bbox, tframe ) - orig = np.copy( twoframes[1] ) twoframes[1] = new_cells - ut.keep_orphans( twoframes, orig, [] ) + twoframes = self.keep_orphans( twoframes, tframe ) parents = self.get_parents( twoframes, labels ) childs = [None]*len(labels) if tframe < (self.epicure.nframes-1): twoframes = np.copy( ut.cropBBox2D(self.epicure.seglayer.data[tframe+1], bbox) ) twoframes = np.stack( (twoframes, np.copy(new_cells)) ) - orig = ut.cropBBox2D(self.epicure.seglayer.data[tframe], bbox) - ut.keep_orphans( twoframes, orig, [] ) + twoframes = self.keep_orphans( twoframes, tframe ) childs = self.get_parents( twoframes, labels ) free_labels = self.epicure.get_free_labels( nlabels ) torelink = [] for i in range( len(labels) ): - print(parents[i]) - print(childs[i]) if (parents[i] is not None) and (childs[i] is not None): - ## the two propagation agrees (gap allowed) - #if parents[i] == childs[i]: - # free_labels[i] = parents[i] - # if self.epicure.verbose > 0: - # print("Link new cell with previous/next "+str(free_labels[i])) free_labels[i] = parents[i] if self.epicure.verbose > 0: print("Link new cell with previous/next "+str(free_labels[i])) - if childs[i] != parents[i]: - torelink.append( [free_labels[i], childs[i]] ) + #if childs[i] != parents[i]: + # torelink.append( [free_labels[i], childs[i]] ) ## only one link found, take it if (parents[i] is not None) and (childs[i] is None): free_labels[i] = parents[i] @@ -356,7 +360,7 @@ class Editing( QWidget ): if self.epicure.verbose > 0: print("Link new cell with previous/next "+str(free_labels[i])) - print(free_labels) + print("Added cells "+str(free_labels)) ## get the new indices and labels to draw new_labels = [] @@ -377,8 +381,8 @@ class Editing( QWidget ): self.epicure.change_labels( indices, new_labels ) ## relink child tracks if necessary - for relink in torelink: - self.epicure.replace_label( relink[1], relink[0], tframe ) + #for relink in torelink: + # self.epicure.replace_label( relink[1], relink[0], tframe ) def touching_masks(self, maska, maskb): """ Check if the two mask touch """ @@ -410,7 +414,8 @@ class Editing( QWidget ): def set_checked(layer, event): if event.type == "mouse_press": if (event.button == 1) and (len(event.modifiers) == 0): - if layer.mode == "paint": + if layer.mode == "paint": + #and not self.napari_painting.isChecked(): ### Overwrite the painting to check that everything stays within EpiCure constraints if self.shapelayer_name not in self.viewer.layers: self.create_shapelayer() @@ -830,7 +835,6 @@ class Editing( QWidget ): def split_along_line(self, tframe, positions): """ Split a label along a line drawn manually """ - bbox = ut.getBBox2DFromPts( positions, extend=0, imshape=self.epicure.imgshape2D ) bbox = ut.extendBBox2D( bbox, extend_factor=1.25, imshape=self.epicure.imgshape2D ) @@ -915,9 +919,6 @@ class Editing( QWidget ): bbox = ut.extendBBox2D( bbox, extend_factor, self.epicure.imgshape2D ) segt_crop = ut.cropBBox2D( segt, bbox ) - if self.epicure.verbose > 1: - ut.show_duration(start_time, "Merging: bbox cropped ") - ## check that labels can be merged touch = ut.checkTouchingLabels( segt_crop, startlab, endlab ) if not touch: @@ -1049,7 +1050,7 @@ class Editing( QWidget ): def create_shapelayer( self ): """ Create the layer that handle temporary drawings """ shapes = [] - shap = self.viewer.add_shapes( shapes, name=self.shapelayer_name, ndim=3, blending="additive", opacity=1, edge_width=2 ) + shap = self.viewer.add_shapes( shapes, name=self.shapelayer_name, ndim=3, blending="additive", opacity=1, edge_width=2, scale=self.viewer.layers["Segmentation"].scale ) shap.text.visible = False shap.visible = False @@ -1086,7 +1087,7 @@ class Editing( QWidget ): def create_seedlayer(self): pts = [] - points = self.viewer.add_points( np.array(pts), face_color="blue", size = 7, edge_width=0, name="Seeds" ) + points = self.viewer.add_points( np.array(pts), face_color="blue", size = 7, edge_width=0, name="Seeds", scale=self.viewer.layers["Segmentation"].scale ) def reset_seeds(self): ut.remove_layer(self.viewer, "Seeds") @@ -1335,21 +1336,12 @@ class Editing( QWidget ): ut.show_info( "Resetting everything ") self.viewer.window._status_bar._toggle_activity_dock(True) progress_bar = progress(total=5) - progress_bar.set_description("Reset: get skeleton") - progress_bar.update(0) ## get skeleton and relabel (ensure label unicity) - #skel = np.zeros(self.epicure.seg.shape, dtype="uint8") - #skel[self.epicure.seg==0] = 1 - skel = self.epicure.get_skeleton() - skel = np.uint32( skel ) - self.epicure.seg = skel - self.epicure.seglayer.data = skel progress_bar.update(1) progress_bar.set_description("Reset: relabel") self.epicure.reset_data() self.epicure.tracking.reset() - self.epicure.junctions_to_label() - self.epicure.seglayer.data = self.epicure.seg + self.epicure.reset_labels() progress_bar.update(2) progress_bar.set_description("Reset: reinit tracks") self.epicure.tracked = 0 @@ -1515,7 +1507,7 @@ class Editing( QWidget ): ut.remove_layer(self.viewer, self.grouplayer_name) data = self.epicure.seglayer.data grouped = self.epicure.draw_groups() - grouplayer = self.viewer.add_labels(grouped, name=self.grouplayer_name, opacity=0.75, blending="additive") + grouplayer = self.viewer.add_labels(grouped, name=self.grouplayer_name, opacity=0.75, blending="additive", scale=self.viewer.layers["Segmentation"].scale) ut.set_active_layer(self.viewer, "Segmentation") else: ut.remove_layer(self.viewer, self.grouplayer_name) @@ -1600,16 +1592,16 @@ class Editing( QWidget ): twoframes[1] = crop_merge # merge of the labels and 0 outside ## keep only parent labels that stop at the previous frame - orig_frame = ut.cropBBox2D(self.epicure.seglayer.data[frame], bbox) - ut.keep_orphans(twoframes, orig_frame, []) + twoframes = self.keep_orphans(twoframes, frame-1, parent=True) ## do mini-tracking to assign most likely parent parent = self.get_parents( twoframes, [1] ) if self.epicure.verbose > 0: print( "Found parent "+str(parent[0])+" to clicked cells "+str(labela)+" and "+str(labelb) ) ## add division to graph - self.epicure.tracking.add_division( labela, labelb, parent[0] ) - ## add division to event list (if active) - self.epicure.inspecting.add_division( labela, labelb, parent[0], frame ) + if parent is not None and parent[0] is not None: + self.epicure.tracking.add_division( labela, labelb, parent[0] ) + ## add division to event list (if active) + self.epicure.inspecting.add_division( labela, labelb, parent[0], frame ) def key_tracking_binding(self): """ active key bindings for tracking options """ @@ -1803,59 +1795,71 @@ class Editing( QWidget ): """ Merge track with label a with track of label b, temporally or spatially """ + if labela == labelb: + if self.epicure.verbose > 0: + print("Already the same track" ) + return if int(posb[0]) == int(posa[0]): - self.tracks_spatial_merging( labela, posa, labelb, posb ) + self.tracks_spatial_merging( labela, posa, labelb ) else: self.tracks_temporal_merging( labela, posa, labelb, posb ) - def tracks_spatial_merging( self, labela, posa, labelb, posb ): + def tracks_spatial_merging( self, labela, posa, labelb ): """ Merge spatially two tracks: labels have to be touching all along the common frames """ + start_time = ut.start_time() ## get last common frame lasta = self.epicure.tracking.get_last_frame( labela ) lastb = self.epicure.tracking.get_last_frame( labelb ) lastcommon = min(lasta, lastb) - ## check if labels are touching from first frame (posa[0]) to the last common frame - touched = True - for frame in range( int(posa[0]), lastcommon+1 ): - bbox, mask = ut.getBBox2DMerge( self.epicure.seg[frame], labela, labelb ) - bbox = ut.extendBBox2D( bbox, 1.05, self.epicure.imgshape2D ) - segt_crop = ut.cropBBox2D( self.epicure.seg[frame], bbox ) - touched = ut.checkTouchingLabels( segt_crop, labela, labelb ) - if not touched: - print("Labels "+str(labela)+" and "+str(labelb)+" are not always touching. Refusing to merge them") - return + + ## if longer than the last common, split the label(s) that continue + if lasta > lastcommon: + if self.epicure.tracking.get_first_frame( labela ) < int(posa[0]): + self.epicure.split_track( labela, lastcommon+1 ) + if lastb > lastcommon: + if self.epicure.tracking.get_first_frame( labelb ) < int(posa[0]): + self.epicure.split_track( labelb, lastcommon+1 ) ## Looks, ok, create a new track and merge the two tracks in it new_label = self.epicure.get_free_label() new_labels = [] ind_tomodif = None + footprint = disk(radius=3) for frame in range( int(posa[0]), lastcommon+1 ): bbox, merged = ut.getBBox2DMerge( self.epicure.seg[frame], labela, labelb ) bbox = ut.extendBBox2D( bbox, 1.05, self.epicure.imgshape2D ) + + ## check if labels are touching at each frame segt_crop = ut.cropBBox2D( self.epicure.seg[frame], bbox ) + touched = ut.checkTouchingLabels( segt_crop, labela, labelb ) + if not touched: + print("Labels "+str(labela)+" and "+str(labelb)+" are not always touching. Refusing to merge them") + return ## merge the two labels together joinlab = ut.cropBBox2D( merged, bbox ) - footprint = disk(radius=3) joinlab = new_label * binary_closing(joinlab, footprint) ## get the index and new values to change indmodif = ut.ind_boundaries( joinlab ) + #indmodif = ut.toFullMoviePos( indmodif, bbox, frame ) + new_labels = new_labels + [0]*len(indmodif) + curmodif = np.transpose( np.nonzero( joinlab == new_label ) ) + new_labels = new_labels + [new_label]*len(curmodif) + indmodif = np.vstack((indmodif, curmodif)) indmodif = ut.toFullMoviePos( indmodif, bbox, frame ) if ind_tomodif is None: ind_tomodif = indmodif else: ind_tomodif = np.vstack((ind_tomodif, indmodif)) - new_labels = new_labels + np.repeat(0, len(indmodif)).tolist() - curmodif = np.argwhere( joinlab == new_label ) - new_labels = new_labels + [new_label]*curmodif.shape[0] - curmodif = ut.toFullMoviePos( curmodif, bbox, frame ) - ind_tomodif = np.vstack((ind_tomodif, curmodif)) + #ind_tomodif = np.vstack((ind_tomodif, curmodif)) ## update the labels and the tracks - self.epicure.change_labels(ind_tomodif, new_labels) + self.epicure.change_labels_frommerge( ind_tomodif, new_labels, remove_labels=[labela, labelb] ) if self.epicure.verbose > 0: ut.show_info("Merged spatially "+str(labela)+" with "+str(labelb)+" from frame "+str(int(posa[0]))+" to frame "+str(lastcommon)+"\n New track label is "+str(new_label)) + if self.epicure.verbose > 1: + ut.show_duration(start_time, "Merging spatially tracks in ") def tracks_temporal_merging( self, labela, posa, labelb, posb ): @@ -1916,7 +1920,21 @@ class Editing( QWidget ): new_labels = new_labels + ([parent_label]*curmodif.shape[0]) return indmodif, new_labels, parent_labels - def inherit_parent_labels(self, myframe, labels, bbox, frame): + def keep_orphans( self, img, frame, keep_labels=[]): + """ Keep only labels that doesn't have a follower (track is finishing at that frame) """ + ## remove the labels to track + labs = np.unique(img[0]).tolist() #np.setdiff1d( img[0], labels ).tolist() + if 0 in labs: + labs.remove(0) + ## Check that it's not present at current frame + torem = [ lab for lab in labs if (lab not in keep_labels) and (self.epicure.tracking.is_in_frame( lab, frame ) ) ] + if len(torem) == 0: + return img + mask = np.isin(img[0], torem) + img[0][mask] = 0 + return img + + def inherit_parent_labels(self, myframe, labels, bbox, frame, keep_labels): """ Get parent labels if any and indices to modify with it """ if ( self.epicure.tracked == 0 ) or (frame<=0): parent_labels = [None]*len(labels) @@ -1924,8 +1942,7 @@ class Editing( QWidget ): else: twoframes = ut.crop_twoframes( self.epicure.seglayer.data, bbox, frame ) twoframes[1] = np.copy(myframe) # merge of the labels and 0 outside - orig_frame = ut.cropBBox2D(self.epicure.seglayer.data[frame], bbox) - ut.keep_orphans(twoframes, orig_frame, labels) + twoframes = self.keep_orphans( twoframes, frame, keep_labels=keep_labels) parent_labels = self.get_parents( twoframes, labels ) @@ -1951,9 +1968,7 @@ class Editing( QWidget ): return [], [] twoframes = np.stack( (twoframes, np.copy(myframe)) ) - orig_frame = ut.cropBBox2D(self.epicure.seglayer.data[frame], bbox) - prev = np.copy(twoframes[0]) - ut.keep_orphans(twoframes, orig_frame, keep_labels) + twoframes = self.keep_orphans(twoframes, frame, keep_labels=keep_labels) child_labels = self.get_parents( twoframes, labels ) if self.epicure.verbose > 0: @@ -1966,9 +1981,19 @@ class Editing( QWidget ): new_labels = [] for i, lab in enumerate(child_labels): if lab is not None: + if lab == parent_labels[i]: + ## going to propagate to itself, no need + continue after_frame = frame+1 + last_frame = self.epicure.tracking.get_last_frame( parent_labels[i] ) + if (last_frame is not None) and (last_frame >= after_frame): + ## the label to propagate is present somewhere after the current frame + self.epicure.split_track( parent_labels[i], after_frame ) inds = self.epicure.get_label_indexes( lab, after_frame ) - indmodif = indmodif + inds + if len(indmodif) == 0: + indmodif = inds + else: + indmodif = np.vstack((indmodif, inds)) new_labels = new_labels + np.repeat(parent_labels[i], len(inds)).tolist() return indmodif, new_labels @@ -1981,16 +2006,16 @@ class Editing( QWidget ): new_labels = np.repeat(0, len(indmodif)).tolist() ## get parent labels if any for each label - indmodif2, new_labels2, parent_labels = self.inherit_parent_labels(myframe, labels, bbox, frame) + indmodif2, new_labels2, parent_labels = self.inherit_parent_labels(myframe, labels, bbox, frame, keep_labels) if indmodif2 is not None: indmodif = np.vstack((indmodif, indmodif2)) new_labels = new_labels+new_labels2 if self.epicure.verbose > 1: ut.show_duration(start_time, "Propagation, parents found, ") - ## propagate the split: get child labels if any for each label + ## propagate the change: get child labels if any for each label indmodif_child, new_labels_child = self.inherit_child_labels(myframe, labels, bbox, frame, parent_labels, keep_labels) - if indmodif_child != []: + if len(indmodif_child) > 0: indmodif = np.vstack((indmodif, indmodif_child)) new_labels = new_labels + new_labels_child if self.epicure.verbose > 1: diff --git a/src/epicure/epicuring.py b/src/epicure/epicuring.py index 9fa9c87e4e353781d71c1df0e017f2a742d2506f..ef9f2012b8daf9b53583b07b3e62c1f87f5b8c40 100644 --- a/src/epicure/epicuring.py +++ b/src/epicure/epicuring.py @@ -1,16 +1,16 @@ import numpy as np import os, time, pickle import napari +import math from qtpy.QtWidgets import QPushButton, QVBoxLayout, QTabWidget, QWidget from napari.utils import progress from magicgui.widgets import TextEdit from skimage.segmentation import find_boundaries, watershed from skimage.measure import label -from skimage.morphology import skeletonize, binary_closing, binary_opening +from skimage.morphology import skeletonize, binary_closing from skimage.util import apply_parallel -from skimage.measure import regionprops, regionprops_table +from skimage.measure import regionprops from multiprocessing.pool import ThreadPool as Pool -import pandas as pand import epicure.Utils as ut from epicure.editing import Editing @@ -76,6 +76,10 @@ class EpiCure(): self.epi_metadata["MainChannel"] = 0 self.epi_metadata["Allow gaps"] = True self.epi_metadata["Verbose"] = 1 + self.epi_metadata["Scale bar"] = True + self.epi_metadata["MovieFile"] = "" + self.epi_metadata["SegmentationFile"] = "" + self.epi_metadata["EpithelialCells"] = True ## epithelial (packed) cells def get_resetbtn_color( self ): """ Returns the color of Reset buttons if defined """ @@ -91,8 +95,8 @@ class EpiCure(): def load_movie(self, imgpath): """ Load the intensity movie, and get metadata """ - self.imgpath = imgpath - self.img, nchan, self.epi_metadata["ScaleXY"], self.epi_metadata["UnitXY"], self.epi_metadata["ScaleT"], self.epi_metadata["UnitT"] = ut.opentif( self.imgpath, verbose=self.verbose>1 ) + self.epi_metadata["MovieFile"] = imgpath + self.img, nchan, self.epi_metadata["ScaleXY"], self.epi_metadata["UnitXY"], self.epi_metadata["ScaleT"], self.epi_metadata["UnitT"] = ut.opentif( self.epi_metadata["MovieFile"], verbose=self.verbose>1 ) ## transform static image to movie (add temporal dimension) if len(self.img.shape) == 2: self.img = np.expand_dims(self.img, axis=0) @@ -134,12 +138,28 @@ class EpiCure(): self.epi_metadata["Allow gaps"] = allow_gap self.forbid_gaps = not allow_gap + def set_epithelia( self, epithelia ): + """ Set the mode for cell packing (touching or not especially) """ + self.epi_metadata["EpithelialCells"] = epithelia + + def set_scalebar( self, show_scalebar ): + """ Show or not the scale bar """ + self.epi_metadata["Scale bar"] = show_scalebar + if self.viewer is not None: + self.viewer.scale_bar.visible = show_scalebar + self.viewer.scale_bar.unit = self.epi_metadata["UnitXY"] + for lay in self.viewer.layers: + lay.scale = [1, self.epi_metadata["ScaleXY"], self.epi_metadata["ScaleXY"]] + self.viewer.reset_view() + def set_scales( self, scalexy, scalet, unitxy, unitt ): """ Set the scaling units for outputs """ self.epi_metadata["ScaleXY"] = scalexy self.epi_metadata["ScaleT"] = scalet self.epi_metadata["UnitXY"] = unitxy self.epi_metadata["UnitT"] = unitt + if self.viewer is not None: + self.viewer if self.verbose > 0: ut.show_info( "Movie scales set to "+str(self.epi_metadata["ScaleXY"])+" "+self.epi_metadata["UnitXY"]+" and "+str(self.epi_metadata["ScaleT"])+" "+self.epi_metadata["UnitT"] ) @@ -173,8 +193,8 @@ class EpiCure(): def load_segmentation(self, segpath): """ Load the segmentation file """ start_time = ut.start_time() - self.segpath = segpath - self.seg,_, _,_,_,_ = ut.opentif( self.segpath, verbose=self.verbose>1 ) + self.epi_metadata["SegmentationFile"] = segpath + self.seg,_, _,_,_,_ = ut.opentif( segpath, verbose=self.verbose>1 ) self.seg = np.uint32(self.seg) ## transform static image to movie (add temporal dimension) if len(self.seg.shape) == 2: @@ -189,21 +209,19 @@ class EpiCure(): self.tracked = 0 else: self.has_been_tracked() - self.thin_boundaries() - - if np.max(self.seg) < 50000: - self.dtype = np.uint16 - self.seg = np.uint16(self.seg) + self.prepare_labels() ## define a reference size of the movie to scale default parameters self.reference_size = np.max( self.imgshape2D ) + ## define the average cell radius + #self.cell_avg_radius = int( math.sqrt( ut.average_area( self.seg[0])/math.pi ) ) # on cell area ut.summary_labels( self.seg[0] ) - if self.verbose > 1: - print("Reference size of the movie: "+str(self.reference_size)) + #if self.verbose > 1: + # print("Reference size of the movie: "+str(self.reference_size)) # display the segmentation file movie if self.viewer is not None: - self.seglayer = self.viewer.add_labels( self.seg, name="Segmentation", blending="additive", opacity=0.5 ) + self.seglayer = self.viewer.add_labels( self.seg, name="Segmentation", blending="additive", opacity=0.5, scale=self.viewer.layers["Movie"].scale ) self.viewer.dims.set_point(0,0) self.seglayer.brush_size = 4 ## default label pencil drawing size if self.verbose > 0: @@ -221,7 +239,6 @@ class EpiCure(): self.handle_gaps( track_list=None, verbose=1 ) ut.show_info(""+str(len(self.tracking.get_track_list()))+" "+tracked+" cells loaded") - def has_been_tracked(self): """ Look if has been tracked already (some labels are in several frames) """ nb = 0 @@ -236,7 +253,9 @@ class EpiCure(): def suggest_segfile(self, outdir): """ Check if a segmentation file from EpiCure already exists """ - imgname, imgdir, out = ut.extract_names( self.imgpath, outdir, mkdir=False ) + if (self.epi_metadata["SegmentationFile"] != "") and ut.found_segfile(self.epi_metadata["SegmentationFile"]): + return self.epi_metadata["SegmentationFile"] + imgname, imgdir, out = ut.extract_names( self.epi_metadata["MovieFile"], outdir, mkdir=False ) return ut.suggest_segfile(out, imgname) def outname(self): @@ -244,7 +263,7 @@ class EpiCure(): def set_names(self, outdir): """ Extract default names from imgpath """ - self.imgname, self.imgdir, self.outdir = ut.extract_names( self.imgpath, outdir, mkdir=True ) + self.imgname, self.imgdir, self.outdir = ut.extract_names( self.epi_metadata["MovieFile"], outdir, mkdir=True ) def go_epicure( self, outdir="epics", segmentation_file=None ): """ Initialize everything and start the main widget """ @@ -256,6 +275,7 @@ class EpiCure(): progress_bar.set_description( "Reading segmented image" ) ## load the segmentation self.load_segmentation( segmentation_file ) + self.epi_metadata["SegmentationFile"] = segmentation_file progress_bar.update(1) ut.set_active_layer( self.viewer, "Segmentation" ) @@ -441,13 +461,13 @@ class EpiCure(): """ Get a summary of the infos of the movie """ summ = "----------- EpiCure summary ----------- \n" summ += "--- Image infos \n" - summ += "Movie name: "+str(self.imgpath)+"\n" + summ += "Movie name: "+str(self.epi_metadata["MovieFile"])+"\n" summ += "Movie size (x,y): "+str(self.imgshape2D)+"\n" if self.nframes is not None: summ += "Nb frames: "+str(self.nframes)+"\n" summ += "\n" summ += "--- Segmentation infos \n" - summ += "Segmentation file: "+str(self.segpath)+"\n" + summ += "Segmentation file: "+str(self.epi_metadata["SegmentationFile"])+"\n" summ += "Nb tracks: "+str(len(self.tracking.get_track_list()))+"\n" tracked = "yes" if self.tracked == 0: @@ -459,6 +479,7 @@ class EpiCure(): summ += "Average cell area: "+str(mean_area)+" pixels^2\n" summ += "Nb suspect events: "+str(self.inspecting.nb_events(only_suspect=True))+"\n" summ += "Nb divisions: "+str(self.inspecting.nb_type("division"))+"\n" + summ += "Nb extrusions: "+str(self.inspecting.nb_type("extrusion"))+"\n" summ += "\n" summ += "--- Parameter infos \n" summ += "Junction thickness: "+str(self.thickness)+"\n" @@ -482,14 +503,14 @@ class EpiCure(): if "Movie" not in self.viewer.layers: if self.verbose > 0: print("Reput movie layer") - mview = self.viewer.add_image( self.img, name="Movie", blending="additive", colormap="gray" ) + mview = self.viewer.add_image( self.img, name="Movie", blending="additive", colormap="gray", scale=[1, self.epi_metadata["ScaleXY"], self.epi_metadata["ScaleXY"]]) #mview.reset_contrast_limits() mview.contrast_limits=self.quantiles() mview.gamma=0.95 if "Segmentation" not in self.viewer.layers: if self.verbose > 0: print("Reput segmentation") - self.seglayer = self.viewer.add_labels( self.seg, name="Segmentation", blending="additive", opacity=0.5 ) + self.seglayer = self.viewer.add_labels( self.seg, name="Segmentation", blending="additive", opacity=0.5, scale=self.viewer.layers["Movie"].scale ) self.finish_update() @@ -722,7 +743,37 @@ class EpiCure(): skel[[0, -1], :] = bin[[0, -1], :] # top and bottom borders skel[:, [0, -1]] = bin[:, [0, -1]] # left and right borders return skel + + def reset_labels( self ): + """ Reset all labels, ensure unicity """ + if self.epi_metadata["EpithelialCells"]: + ### packed (contiguous cells), ensure that they are separated by one pixel only + skel = self.get_skeleton() + skel = np.uint32( skel ) + self.seg = skel + self.seglayer.data = skel + self.junctions_to_label() + self.seglayer.data = self.seg + else: + self.get_cells() + + def prepare_labels(self): + """ Process the labels to be in a correct Epicurable format """ + if self.epi_metadata["EpithelialCells"]: + ### packed (contiguous cells), ensure that they are separated by one pixel only + self.thin_boundaries() + else: + self.get_cells() + + def get_cells( self ): + """ Non jointive cells: check label unicity """ + for frame in self.seg: + if ut.non_unique_labels( frame ): + self.seg = ut.reset_labels( self.seg ) + return + + def thin_boundaries(self): """" Assure that all boundaries are only 1 pixel thick """ if self.process_parallel: @@ -733,25 +784,12 @@ class EpiCure(): for z in range(self.seg.shape[0]): self.thin_seg_one_frame(z) - def skeleton_to_label( self, skel, labelled ): - """ Transform a skeleton to label image with numbers from labelled image """ - newlab = label( skel, background=np.max(skel), connectivity=1 ) - props = regionprops( newlab ) - newlab = np.zeros( newlab.shape, np.uint32 ) - for prop in props: - if prop.label != 0: - labvals, counts = np.unique(labelled[prop.slice][prop.image], return_counts=True ) - labval = labvals[ np.argmax(counts) ] - newlab[prop.slice][prop.image] = labval - return newlab - def thin_seg_one_frame(self, tframe): """ Boundaries of the frame one pixel thick """ bin_img = binary_closing( find_boundaries(self.seg[tframe], connectivity=2, mode="outer"), footprint=np.ones((3,3)) ) skel = skeletonize( bin_img ) skel = self.copy_border( skel, bin_img ) - self.seg[tframe] = self.skeleton_to_label( skel, self.seg[tframe] ) - + self.seg[tframe] = ut.skeleton_to_label( skel, self.seg[tframe] ) def add_skeleton(self): """ add a layer containing the skeleton movie of the segmentation """ @@ -759,17 +797,17 @@ class EpiCure(): if self.viewer is not None: skel = np.zeros(self.seg.shape, dtype="uint8") skel[self.seg==0] = 1 - skel = self.get_skeleton() + skel = self.get_skeleton( viewer=self.viewer ) ut.remove_layer(self.viewer, "Skeleton") - skellayer = self.viewer.add_image( skel, name="Skeleton", blending="additive", opacity=1 ) + skellayer = self.viewer.add_image( skel, name="Skeleton", blending="additive", opacity=1, scale=self.viewer.layers["Movie"].scale ) skellayer.reset_contrast_limits() skellayer.contrast_limits = (0,1) - def get_skeleton(self): + def get_skeleton(self, viewer=None): """ convert labels movie to skeleton (thin boundaries) """ if self.seg is None: return None - return ut.get_skeleton( self.seg, verbose = self.verbose ) + return ut.get_skeleton( self.seg, viewer = viewer, verbose = self.verbose ) def to_skeleton(self, mode): """ convert labels movie to skeleton (thin boundaries) """ @@ -834,29 +872,56 @@ class EpiCure(): def get_label_indexes (self, label, start_frame=0 ): """ Returns the indexes where label is present in segmentation, starting from start_frame """ indmodif = [] - while start_frame < self.nframes: - found_label = np.argwhere( self.seg[start_frame] == label ) - if found_label is None: - break ### if lose a label, stop there (don't allow gap frame) - inds = ut.shiftFrameIndices( found_label.tolist(), start_frame ) - indmodif = indmodif + inds - start_frame = start_frame + 1 + if self.verbose > 2: + start_time = ut.start_time() + pos = self.tracking.get_track_column( track_id=label, column="fullpos" ) + pos = pos[pos[:,0]>=start_frame] + ## if nothing in pos, pb with track data + if pos is None or len(pos)==0: + ut.show_warning("Something wrong in the track data. Resetting track data (can take time)") + self.tracking.reset_tracks() + self.get_label_indexes( label, start_frame ) + + indmodif = np.argwhere( self.seg[pos[:,0]] == label ) + indmodif = ut.shiftFrames( indmodif, pos[:,0] ) + if self.verbose > 2: + ut.show_duration(start_time, header="Label indexes found in ") return indmodif def replace_label( self, label, new_label, start_frame=0 ): - """ Replace label with new_label from start_frame """ + """ Replace label with new_label from start_frame - Relabelling only """ indmodif = self.get_label_indexes( label, start_frame ) new_labels = [new_label]*len(indmodif) - self.change_labels( indmodif, new_labels ) + self.change_labels( indmodif, new_labels, replacing=True ) + + def change_labels_frommerge(self, indmodif, new_labels, remove_labels ): + """ Change the value at pixels indmodif to new_labels and update tracks/graph. Full remove of the two merged labels """ + if len(indmodif)>0: + ## get effectively changed labels + indmodif, new_labels, _ = ut.setNewLabel( self.seglayer, indmodif, new_labels, add_frame=None, return_old=False ) + if len(new_labels) > 0: + self.update_added_labels( indmodif, new_labels ) + self.update_removed_labels( indmodif, remove_labels ) + self.seglayer.refresh() - def change_labels(self, indmodif, new_labels): - """ Change the value at pixels indmodif to new_labels and update tracks/graph """ + def change_labels(self, indmodif, new_labels, replacing=False): + """ Change the value at pixels indmodif to new_labels and update tracks/graph + + Assume that only label at current frame can have its shape modified. Other changed label is only relabelling at frames > current frame (child propagation) + """ if len(indmodif)>0: - bbox = ut.getBBoxFromPts( indmodif, extend=0, imshape=self.imgshape, outdim=3 ) ## get effectively changed labels indmodif, new_labels, old_labels = ut.setNewLabel( self.seglayer, indmodif, new_labels, add_frame=None ) if len(new_labels) > 0: - self.update_changed_labels( indmodif, new_labels, old_labels ) + if replacing: + self.update_replaced_labels( indmodif, new_labels, old_labels ) + else: + ## the only label to change are the current frame (smaller one), the other are only relabelling (propagation) + cur_frame = np.min( indmodif[0] ) + to_reshape = (indmodif[0] == cur_frame) + self.update_changed_labels( (indmodif[0][to_reshape], indmodif[1][to_reshape], indmodif[2][to_reshape]), new_labels[to_reshape], old_labels[to_reshape] ) + to_relab = np.invert(to_reshape) + self.update_replaced_labels( (indmodif[0][to_relab], indmodif[1][to_relab], indmodif[2][to_relab]), new_labels[to_relab], old_labels[to_relab] ) self.seglayer.refresh() def get_mask( self, label, start=None, end=None ): @@ -952,13 +1017,13 @@ class EpiCure(): def update_label(self, label, frame): """ Update the given label at given frame """ self.tracking.update_track_on_frame([label], frame) - - def update_changed_labels( self, indmodif, new_labels, old_labels ): - """ Check what had been modified, and update tracks from it """ + + def update_changed_labels( self, indmodif, new_labels, old_labels, full=False ): + """ Check what had been modified, and update tracks from it, looking frame by frame """ ## check all the old_labels if still present or not if self.verbose > 1: start_time = time.time() - frames = np.sort( np.unique( indmodif[0] ) ) + frames = np.unique( indmodif[0] ) all_deleted = [] debug_verb = self.verbose > 2 if debug_verb: @@ -967,23 +1032,70 @@ class EpiCure(): keep = indmodif[0] == frame ## check old labels if totally removed or not deleted = np.setdiff1d( old_labels[keep], self.seg[frame] ) + left = np.setdiff1d( old_labels[keep], deleted ) if deleted.shape[0] > 0: self.tracking.remove_one_frame( deleted, frame, handle_gaps=False, refresh=False ) - all_deleted = all_deleted + list(set(deleted) - set(all_deleted)) + if self.forbid_gaps: + all_deleted = all_deleted + list(set(deleted) - set(all_deleted)) + if left.shape[0] > 0: + self.tracking.update_track_on_frame( left, frame ) ## now check new labels nlabels = np.unique( new_labels[keep] ) if nlabels.shape[0] > 0: self.tracking.update_track_on_frame( nlabels, frame ) if debug_verb: print("Labels deleted at frame "+str(frame)+" "+str(deleted)+" or added "+str(nlabels)) + + def update_added_labels( self, indmodif, new_labels ): + """ Update tracks of labels that have been fully added """ + if self.verbose > 1: + start_time = time.time() + + ## Deleted labels + frames = np.unique( indmodif[0] ) + self.tracking.add_tracks_fromindices( indmodif, new_labels ) + if self.forbid_gaps: + ## Check if some gaps has been created in tracks (remove middle(s) frame(s)) + added = list(set(new_labels)) + if len(added) > 0: + self.handle_gaps( added, verbose=0 ) - ## Check if some gaps has been created in tracks (remove middle(s) frame(s)) + if self.verbose > 1: + ut.show_duration(start_time, "updated added tracks in ") + + def update_removed_labels( self, indmodif, old_labels ): + """ Update tracks of labels that have been fully removed """ + if self.verbose > 1: + start_time = time.time() + + ## Deleted labels + frames = np.unique( indmodif[0] ) + self.tracking.remove_on_frames( np.unique(old_labels), frames ) + if self.forbid_gaps: + ## Check if some gaps has been created in tracks (remove middle(s) frame(s)) + deleted = list(set(old_labels)) + if len(deleted) > 0: + self.handle_gaps( deleted, verbose=0 ) + + if self.verbose > 1: + ut.show_duration(start_time, "updated removed tracks in ") + + def update_replaced_labels( self, indmodif, new_labels, old_labels ): + """ Old_labels were fully replaced by new_labels on some frames, update tracks from it """ + if self.verbose > 1: + start_time = time.time() + + ## Deleted labels + frames = np.unique( indmodif[0] ) + self.tracking.replace_on_frames( np.unique(old_labels), np.unique(new_labels), frames ) if self.forbid_gaps: - if len(all_deleted) > 0: - self.handle_gaps( all_deleted, verbose=0 ) + ## Check if some gaps has been created in tracks (remove middle(s) frame(s)) + deleted = list(set(old_labels)) + if len(deleted) > 0: + self.handle_gaps( deleted, verbose=0 ) if self.verbose > 1: - ut.show_duration(start_time, "updated tracks in ") + ut.show_duration(start_time, "updated replaced tracks in ") def handle_gaps(self, track_list, verbose=None): """ Check and fix gaps in tracks """ diff --git a/src/epicure/epiwidgets.py b/src/epicure/epiwidgets.py index 0ff8bc8c5840cf667dc0dc70ff86b430091a1c0c..aa65867fc9ea52341b68b54b7c61f7f3f1eba47e 100644 --- a/src/epicure/epiwidgets.py +++ b/src/epicure/epiwidgets.py @@ -1,6 +1,6 @@ import epicure.Utils as ut -from qtpy.QtWidgets import QPushButton, QCheckBox, QHBoxLayout, QVBoxLayout, QLabel, QLineEdit, QComboBox, QSpinBox, QSlider, QGroupBox -from qtpy.QtCore import Qt +from qtpy.QtWidgets import QPushButton, QCheckBox, QHBoxLayout, QVBoxLayout, QLabel, QLineEdit, QComboBox, QSpinBox, QSlider, QGroupBox # type: ignore +from qtpy.QtCore import Qt # type: ignore def help_button( link, description="", display_settings=None ): """ Create a new Help button with given parameter """ @@ -88,6 +88,10 @@ def hlayout(): """ Return a horizontal layout """ return QHBoxLayout() +def vlayout(): + """ Return a vertical layout """ + return QVBoxLayout() + def add_check_tolayout( layout, check, checked, check_func=None, descr="" ): """ Add a checkbox with set parameters """ cbox = add_check( check, checked, check_func, descr ) diff --git a/src/epicure/inspecting.py b/src/epicure/inspecting.py index 2c8ea1caee7389b869c198dd63b75e2f23af1fb3..654555ab9855e7f48eaab1e3b5b19adfd8d8f363 100644 --- a/src/epicure/inspecting.py +++ b/src/epicure/inspecting.py @@ -160,7 +160,7 @@ class Inspecting(QWidget): """ Create a point layer that contains the events """ features = {} pts = [] - self.events = self.viewer.add_points( np.array(pts), properties=features, face_color="red", size = 10, symbol='x', name=self.eventlayer_name, ) + self.events = self.viewer.add_points( np.array(pts), properties=features, face_color="red", size = 10, symbol='x', name=self.eventlayer_name, scale=self.viewer.layers["Segmentation"].scale ) self.event_types = {} self.update_nevents_display() self.epicure.finish_update() @@ -173,7 +173,7 @@ class Inspecting(QWidget): if colors is None: colors = "red" symb = symbols - self.events = self.viewer.add_points( np.array(pts), properties=features, face_color=colors, size = 10, symbol=symbols, name=self.eventlayer_name, ) + self.events = self.viewer.add_points( np.array(pts), properties=features, face_color=colors, size = 10, symbol=symbols, name=self.eventlayer_name, scale=self.viewer.layers["Segmentation"].scale ) self.event_types = event_types self.update_nevents_display() self.epicure.finish_update() @@ -235,6 +235,7 @@ class Inspecting(QWidget): self.color_choice.addItem("track-1-2-*") self.color_choice.addItem("track-length") self.color_choice.addItem("track-gap") + self.color_choice.addItem("track-jump") self.color_choice.addItem("division") self.color_choice.addItem("area") self.color_choice.addItem("solidity") @@ -242,7 +243,7 @@ class Inspecting(QWidget): self.color_choice.addItem("tubeness") disp_layout.addLayout(colorlay) - esize = int(self.epicure.reference_size/75+10) + esize = int(self.epicure.reference_size/70+10) msize = 100 if esize > 70: msize = 200 @@ -393,12 +394,16 @@ class Inspecting(QWidget): for i, eclass in enumerate(self.event_class): disp["Show "+eclass] = self.show_class[i].isChecked() disp["Ignore border"] = self.ignore_borders.isChecked() + disp["Ignore boundaries"] = self.ignore_boundaries.isChecked() disp["Flag length"] = self.check_length.isChecked() + disp["Flag jump"] = self.check_jump.isChecked() disp["length"] = self.min_length.text() disp["Check size"] = self.check_size.isChecked() disp["Check shape"] = self.check_shape.isChecked() + disp["Get merging"] = self.get_merge.isChecked() disp["Get apparitions"] = self.get_apparition.isChecked() disp["Get disparitions"] = self.get_disparition.isChecked() + disp["Get extrusions"] = self.get_extrusions.isChecked() disp["Get gaps"] = self.get_gaps.isChecked() disp["threshold disparition"] = self.threshold_disparition.text() disp["Min gap"] = self.min_gaps.text() @@ -425,18 +430,26 @@ class Inspecting(QWidget): self.show_hide_events() if setting == "Ignore border": self.ignore_borders.setChecked( val ) + if setting == "Ignore boundaries": + self.ignore_boundaries.setChecked( val ) if setting == "Flag length": self.check_length.setChecked( val ) + if setting == "Flag jump": + self.check_jump.setChecked( val ) if setting == "length": self.min_length.setText( val ) if setting == "Check size": self.check_size.setChecked( val ) if setting == "Check shape": self.check_shape.setChecked( val ) + if setting == "Get merging": + self.get_merge.setChecked( val ) if setting == "Get apparitions": self.get_apparition.setChecked( val ) if setting == "Get disparitions": self.get_disparition.setChecked( val ) + if setting == "Get extrusions": + self.get_extrusions.setChecked( val ) if setting == "Get gaps": self.get_gaps.setChecked( val ) if setting == "Threshold disparition": @@ -484,7 +497,7 @@ class Inspecting(QWidget): features["label"] = np.array([label], dtype=self.epicure.dtype) features["score"] = np.array([0], dtype="uint8") pts = [pos] - self.events = self.viewer.add_points( np.array(pts), properties=features, face_color="score", size = int( self.event_size.value() ), symbol="x", name="Events", ) + self.events = self.viewer.add_points( np.array(pts), properties=features, face_color="score", size = int( self.event_size.value() ), symbol="x", name="Events", scale=self.viewer.layers["Segmentation"].scale ) self.add_event_type(0, sid, featurename) self.events.refresh() self.update_nevents_display() @@ -495,6 +508,11 @@ class Inspecting(QWidget): tframe = int(pos[0]) if label in self.border_cells[tframe]: return + + if (not force) and (self.ignore_boundaries.isChecked()) and (self.boundary_cells is not None): + tframe = int(pos[0]) + if label in self.boundary_cells[tframe]: + return ## initialise if necessary if len(self.events.data) <= 0: @@ -533,6 +551,7 @@ class Inspecting(QWidget): self.events.refresh() self.update_nevents_display() self.reset_event_range() + self.epicure.finish_update() def new_event_id(self): """ Find the first unused id """ @@ -548,10 +567,11 @@ class Inspecting(QWidget): features = {} pts = [] ut.remove_layer(self.viewer, "Events") - self.events = self.viewer.add_points( np.array(pts), properties=features, face_color="red", size = 10, symbol='x', name="Events", ) + self.events = self.viewer.add_points( np.array(pts), properties=features, face_color="red", size = 10, symbol='x', name="Events", scale=self.viewer.layers["Segmentation"].scale ) self.event_types = {} self.update_nevents_display() #self.update_nevents_display() + self.epicure.finish_update() def reset_event_type(self, feature, frame): """ Remove all event_types of given feature, for current frame or all if frame is None """ @@ -844,19 +864,37 @@ class Inspecting(QWidget): ''' Block interface of functions for error suggestions based on tracks ''' track_layout = QVBoxLayout() - self.ignore_borders = wid.add_check( "Ignore cells on border", False, None, "When adding suspect, don't add it if the cell is touching the border" ) - track_layout.addWidget(self.ignore_borders) + hign = wid.hlayout() + ignore_label = wid.label_line( "Ignore cells on:") + hign.addWidget( ignore_label ) + vign = wid.vlayout() + self.ignore_borders = wid.add_check( "border (of image)", False, None, "When adding suspect, don't add it if the cell is touching the border of the image" ) + vign.addWidget(self.ignore_borders) + + self.ignore_boundaries = wid.add_check( "tissue boundaries", False, None, "When adding suspect, don't add it if the cell is on the tissu boundaries (no neighbor in one side)" ) + vign.addWidget(self.ignore_boundaries) + hign.addLayout( vign ) + track_layout.addLayout( hign ) + + ## Look for merging tracks + self.get_merge = wid.add_check( "Flag track merging", True, None, "Add a suspect if two track merge in one" ) + track_layout.addWidget(self.get_merge) ## Look for sudden appearance of tracks self.get_apparition = wid.add_check( "Flag track apparition", True, None, "Add a suspect if a track appears in the middle of the movie (not on border)" ) track_layout.addWidget(self.get_apparition) ## Look for sudden disappearance of tracks - disp_line, self.get_disparition, self.threshold_disparition = wid.check_value( check="Flag track disparition", checkfunc=None, checked=True, value="200", descr="Add a suspect if a track disappears (not last frame, not border) and the cell area is above threshold", label="if cell area above" ) - track_layout.addLayout( disp_line ) - + dsp_layout = wid.hlayout() + self.get_disparition = wid.add_check( check="Flag track disparition", checked=True, check_func=None, descr="Add a suspect if a track disappears (not last frame, not border)" ) + disp_line, self.threshold_disparition = wid.value_line( label="cell area threshold", default_value="200", descr="Flag cell if cell area is above threshold" ) self.get_extrusions = wid.add_check( "Get extrusions", True, None, "Add extrusions events when a track is disappearing normally (below cell area threshold)" ) - track_layout.addWidget( self.get_extrusions ) + vlay = wid.vlayout() + vlay.addWidget( self.get_disparition ) + vlay.addWidget( self.get_extrusions ) + dsp_layout.addLayout( vlay ) + dsp_layout.addLayout( disp_line ) + track_layout.addLayout( dsp_layout ) ## Look for temporal gaps in tracks gap_line, self.get_gaps, self.min_gaps = wid.check_value( check="Flag track gaps", checkfunc=None, checked=True, value="1", descr="Add a suspect if a track has gaps longer than threshold (in nb of frames)", label="if gap above" ) @@ -866,10 +904,14 @@ class Inspecting(QWidget): ilengthlay, self.check_length, self.min_length = wid.check_value( check="Flag tracks smaller than", checkfunc=None, checked=True, value="1", descr="Add a suspect event for each track smaller than chosen value (in number of frames)" ) track_layout.addLayout(ilengthlay) + ## track sudden jump in position + ijumplay, self.check_jump, self.jump_factor = wid.check_value( check="Flag jump in track position", checkfunc=None, checked=True, value="3.0", descr="Add a suspect event for when the position of cell centroid moves suddenly a lot compared to the rest of the track" ) + track_layout.addLayout(ijumplay) + ## Variability in feature event_type - sizevar_line, self.check_size, self.size_variability = wid.check_value( check="Size variation", checkfunc=None, checked=False, value="1", descr="Add a suspect if the size of the cell varies suddenly in the track" ) + sizevar_line, self.check_size, self.size_variability = wid.check_value( check="Size variation", checkfunc=None, checked=False, value="3", descr="Add a suspect if the size of the cell varies suddenly in the track" ) track_layout.addLayout( sizevar_line ) - shapevar_line, self.check_shape, self.shape_variability = wid.check_value( check="Shape variation", checkfunc=None, checked=False, value="2.0", descr="Add a suspect if the shape of the cell varies suddenly in the track" ) + shapevar_line, self.check_shape, self.shape_variability = wid.check_value( check="Shape variation", checkfunc=None, checked=False, value="3.0", descr="Add a suspect if the shape of the cell varies suddenly in the track" ) track_layout.addLayout( shapevar_line ) ## merge/split combinaisons @@ -885,6 +927,7 @@ class Inspecting(QWidget): self.reset_event_type("track-1-2-*", None) self.reset_event_type("track-2->1", None) self.reset_event_type("track-length", None) + self.reset_event_type("track-jump", None) self.reset_event_type("track-size", None) self.reset_event_type("track-shape", None) self.reset_event_type("track-apparition", None) @@ -896,7 +939,13 @@ class Inspecting(QWidget): """ Find all cells that are only in one frame """ max_len = int(self.min_length.text()) labels, lengths, positions = self.epicure.tracking.get_small_tracks( max_len ) + ## remove track from first and last frame + first_tracks = self.epicure.tracking.get_tracks_on_frame( 0 ) + last_tracks = self.epicure.tracking.get_tracks_on_frame( self.epicure.nframes-1 ) for label, nframe, pos in zip(labels, lengths, positions): + if label in first_tracks or label in last_tracks: + ## present in the first or last track, don't check it + continue if self.epicure.verbose > 2: print("event track length "+str(nframe)+": "+str(label)+" frame "+str(pos[0]) ) self.add_event(pos, label, "track-length", refresh=False) @@ -906,21 +955,22 @@ class Inspecting(QWidget): """ Look for suspicious tracks """ self.viewer.window._status_bar._toggle_activity_dock(True) ut.set_visibility( self.viewer, "Events", True ) - progress_bar = progress( total=9 ) + progress_bar = progress( total=10 ) progress_bar.update(0) self.reset_tracking_event() progress_bar.update(1) - if self.ignore_borders.isChecked(): - progress_bar.set_description("Identifying border cells") - self.get_border_cells() + if self.ignore_borders.isChecked() or self.ignore_boundaries.isChecked(): + progress_bar.set_description("Identifying border and/or boundaries cells") + self.get_outside_cells() progress_bar.update(2) tracks = self.epicure.tracking.get_track_list() if self.check_length.isChecked(): progress_bar.set_description("Identifying too small tracks") self.track_length() progress_bar.update(3) - progress_bar.set_description("Inspect tracks 2->1") - self.track_21() + if self.get_merge.isChecked(): + progress_bar.set_description("Inspect tracks 2->1") + self.track_21() progress_bar.update(4) if (self.check_size.isChecked()) or self.check_shape.isChecked(): progress_bar.set_description("Inspect track features") @@ -930,14 +980,18 @@ class Inspecting(QWidget): progress_bar.set_description("Check new track apparition") self.track_apparition( tracks ) progress_bar.update(6) - if self.get_disparition.isChecked(): - progress_bar.set_description("Check track disparition") + if self.get_disparition.isChecked() or self.get_extrusions.isChecked(): + progress_bar.set_description("Check track disparition and/or extrusion") self.track_disparition( tracks, progress_bar ) progress_bar.update(7) if self.get_gaps.isChecked(): progress_bar.set_description("Check temporal gaps in tracks") self.track_gaps( tracks, progress_bar ) progress_bar.update(8) + if self.check_jump.isChecked(): + progress_bar.set_description("Check position jump in tracks") + self.track_position_jump( tracks, progress_bar ) + progress_bar.update(9) progress_bar.close() self.viewer.window._status_bar._toggle_activity_dock(False) ut.set_active_layer( self.viewer, "Segmentation" ) @@ -951,9 +1005,9 @@ class Inspecting(QWidget): for i, track_id in enumerate( ctracks) : fframe = self.epicure.tracking.get_first_frame( track_id ) ## If on the border, ignore - outside = self.epicure.cell_on_border( track_id, fframe ) - if outside: - continue + #outside = self.epicure.cell_on_border( track_id, fframe ) + #if outside: + # continue ## Not on border, check if potential division if (graph is not None) and (track_id in graph.keys()): continue @@ -980,11 +1034,11 @@ class Inspecting(QWidget): sub_bar.update( i ) lframe = self.epicure.tracking.get_last_frame( track_id ) ## If on the border, ignore - outside = self.epicure.cell_on_border( track_id, lframe ) - if outside: - continue + #outside = self.epicure.cell_on_border( track_id, lframe ) + #if outside: + # continue ## Not on border, check if potential division - if self.epicure.tracking.is_parent( track_id ): + if self.epicure.tracking.is_single_parent( track_id ): continue ## check if the cell area is below the threshold, then considered as ok (likely extrusion) @@ -1000,6 +1054,8 @@ class Inspecting(QWidget): print("Add extrusion: "+str(track_id)+" at frame "+str(lframe) ) self.add_event( pos, track_id, "extrusion", symb="diamond", color="red", refresh=False ) continue + if not self.get_disparition.isChecked(): + continue ## event disparition posxy = self.epicure.tracking.get_position( track_id, lframe ) if posxy is not None: @@ -1061,7 +1117,7 @@ class Inspecting(QWidget): if pos is not None: if self.epicure.verbose > 1: print("event 1->2->1 track: "+str(graph[parent[0]][0])+"-"+str(parent)+"-"+str(child)+" frame "+str(pos[0]) ) - self.add_event(pos, parent[0], "track-1-2-*") + self.add_event(pos, parent[0], "track-1-2-*", refresh=False) onetwoone = True if not onetwoone: @@ -1074,30 +1130,33 @@ class Inspecting(QWidget): if self.epicure.verbose > 1: print("Something weird, "+str(child)+" mean position") - self.epicure.finish_update() self.refresh_events() + def get_outside_cells( self ): + """ Get list of cells on tissu boundaries and/or on border of the movie """ + self.boundary_cells = dict() + self.border_cells = dict() + check_border = self.ignore_borders.isChecked() + check_bound = self.ignore_boundaries.isChecked() + for tframe in range(self.epicure.nframes): + img = self.epicure.seg[tframe] + if check_bound: + self.boundary_cells[tframe] = ut.get_boundary_cells( img ) + if check_border: + self.border_cells[tframe] = ut.get_border_cells( img ) + + def get_boundaries_cells(self): + """ Return list of cells that are at the tissu boundaries (touching background) """ + self.boundary_cells = dict() + for tframe in range(self.epicure.nframes): + self.boundary_cells[tframe] = ut.get_boundary_cells( self.epicure.seg[tframe] ) + def get_border_cells(self): - """ Return list of cells that are at the border (touching background) """ + """ Return list of cells that are at the border of the movie """ self.border_cells = dict() for tframe in range(self.epicure.nframes): img = self.epicure.seg[tframe] - self.border_cells[tframe] = self.get_border_cells_frame(img) - - def get_border_cells_frame(self, imframe): - """ Return cells on border in current image """ - height = self.epicure.imgshape2D[1] - width = self.epicure.imgshape2D[0] - labels = list( np.unique( imframe[ :, 0:2 ] ) ) ## top border - labels += list( np.unique( imframe[ :, (height-2): ] ) ) ## bottom border - labels += list( np.unique( imframe[ 0:2,] ) ) ## left border - labels += list( np.unique( imframe[ (width-2):,] ) ) ## right border - #graph = ut.connectivity_graph( imframe, distance=3 ) - #adj_bg = [] - #if 0 in graph.nodes: - # adj_bg = list( graph.adj[0 ]) - #return adj_bg - return labels + self.border_cells[tframe] = ut.get_border_cells(img) def get_divisions( self ): """ Get and add divisions from the tracking graph """ @@ -1126,10 +1185,10 @@ class Inspecting(QWidget): continue ## get the average first position of the childs just after division pos = self.epicure.tracking.mean_position(indexes, only_first=True) - self.add_event(pos, parent, "division", symb="o", color="#0055ffff", force=True) + self.add_event(pos, parent, "division", symb="o", color="#0055ffff", force=True, refresh=False) ## Update display to show/hide the divisions self.show_hide_divisions() - self.epicure.finish_update() + self.refresh_events() def show_hide_events( self ): """ Update which type of events to show or hide """ @@ -1218,6 +1277,25 @@ class Inspecting(QWidget): if sid in self.event_types[event]: return True return False + + def track_position_jump( self, track_ids, progress_bar ): + """ Look at jump in the track position """ + factor = float( self.jump_factor.text() ) + sub_bar = progress( total = len( track_ids ), desc="Check position jump in tracks", nest_under = progress_bar ) + for i, tid in enumerate(track_ids): + sub_bar.update(i) + track_indexes = self.epicure.tracking.get_track_indexes( tid ) + ## track should be long enough to make sense to look for outlier + if len(track_indexes) > 3: + track_velo = self.epicure.tracking.measure_speed( tid ) + jumps = self.find_jump( track_velo, factor=factor, min_value=5 ) + for tind in jumps: + tdata = self.epicure.tracking.get_frame_data( tid, tind ) + if self.epicure.verbose > 1: + print("event track jump: "+str(tdata[0])+" "+" frame "+str(tdata[1]) ) + self.add_event( tdata[1:4], tid, "track-jump", refresh=False ) + self.refresh_events() + def track_features(self): """ Look at outliers in track features """ @@ -1249,30 +1327,31 @@ class Inspecting(QWidget): tdata = self.epicure.tracking.get_frame_data( tid, out ) if self.epicure.verbose > 1: print("event track "+feature+": "+str(tdata[0])+" "+" frame "+str(tdata[1]) ) - self.add_event(tdata[1:4], tid, "track_"+featType[feature]) + self.add_event(tdata[1:4], tid, "track_"+featType[feature], refresh=False) + self.refresh_events() - def find_jump( self, tab, factor=1 ): + def find_jump( self, tab, factor=1, min_value=None ): """ Detect brutal jump in the values """ jumps = [] tab = np.array(tab) - diff = tab[:(len(tab)-2)] - 2*tab[1:(len(tab)-1)] + tab[2:] - diff = [(tab[1]-tab[0])] + diff.tolist() + [tab[len(tab)-1]-tab[len(tab)-2]] - avg = (tab[:(len(tab)-2)] + tab[2:])/2 - avg = [(tab[1]+tab[0])/2] + avg.tolist() + [(tab[len(tab)-1]+tab[len(tab)-2])/2] - eps = 0.000000001 + diff = np.diff( tab, n=2, prepend=tab[0], append=tab[-1] ) + ## get local average + if len(tab) <= 10: + avg = np.mean( tab ) + else: + kernel = np.repeat (1.0/10.0, 10 ) + avg = np.convolve( tab, kernel, mode="same") + ## normalize the difference by the average value + eps = 0.0000001 diff = np.array(diff, dtype=np.float32) avg = np.array(avg, dtype=np.float32) diff = abs(diff+eps)/(avg+eps) ## keep only local max above threshold - for i, diffy in enumerate(diff): - if (i>0) and (i<len(diff)-1): - if diffy > factor: - if (diffy > diff[i-1]) and (diffy > diff[i+1]): - jumps.append(i) - else: - if diffy > factor: - jumps.append(i) - #jumps = (np.where( diff > factor )[0]).tolist() + local_max = (np.diff( np.sign(np.diff(diff)) )<0).nonzero()[0] + 1 + if min_value is None: + jumps = [i for i in local_max if diff[i] > factor] + else: + jumps = [ i for i in local_max if (diff[i] > factor) and (tab[i] > min_value) ] return jumps def find_outliers_tuk( self, tab, factor=3, below=True, above=True ): diff --git a/src/epicure/laptrack_overlaps.py b/src/epicure/laptrack_overlaps.py index 8d372e1dba7e7c89e1ead4a50d07e3d9cb743aa1..bbf9da35d2551eadd97bd90ef6bc1391b31de308 100644 --- a/src/epicure/laptrack_overlaps.py +++ b/src/epicure/laptrack_overlaps.py @@ -4,18 +4,13 @@ Tracking with laptrack package Inspired from example https://github.com/yfukai/laptrack/blob/main/docs/examples/cell_segmentation.ipynb """ -import napari import numpy as np import pandas as pd from skimage.measure import regionprops_table -import time - from laptrack import OverLapTrack from laptrack import datasets import epicure.Utils as ut -from napari.utils import progress -from napari.utils.notifications import show_info #from multiprocessing.pool import ThreadPool as Pool #from functools import partial diff --git a/src/epicure/outputing.py b/src/epicure/outputing.py index 3c5cbf6d894f150df9ccc370b038f42ce20370ec..20247664e67786caec527f600cb5f160b06d3b98 100644 --- a/src/epicure/outputing.py +++ b/src/epicure/outputing.py @@ -1,6 +1,3 @@ -from qtpy.QtWidgets import QHBoxLayout, QVBoxLayout, QWidget, QTableWidget, QTableWidgetItem, QGridLayout, QListWidget -from qtpy.QtWidgets import QAbstractItemView as aiv -from qtpy.QtCore import Qt import pandas as pand import numpy as np import roifile @@ -10,12 +7,15 @@ import napari from napari.utils import progress import epicure.Utils as ut import epicure.epiwidgets as wid -import matplotlib as mpl -from matplotlib.backends.backend_qt5agg import FigureCanvas -from matplotlib.figure import Figure -import matplotlib.pyplot as plt +import plotly.express as px +from qtpy import QtCore +from qtpy.QtCore import Qt +QtCore.QCoreApplication.setAttribute(QtCore.Qt.AA_ShareOpenGLContexts, True) ## for QtWebEngine import to work on some computers +from qtpy.QtWebEngineWidgets import QWebEngineView +from qtpy.QtWidgets import QHBoxLayout, QVBoxLayout, QWidget, QTableWidget, QTableWidgetItem, QGridLayout, QListWidget +from qtpy.QtWidgets import QAbstractItemView as aiv +from random import sample - class Outputing(QWidget): @@ -308,7 +308,7 @@ class Outputing(QWidget): ## save filled image (for label or skeleton) to file outname = os.path.join( self.epicure.outdir, self.epicure.imgname+endname ) if self.save_choice.currentText() == "skeleton": - tosave = ut.get_skeleton( tosave, verbose=self.epicure.verbose ) + tosave = ut.get_skeleton( tosave, viewer=self.viewer, verbose=self.epicure.verbose ) ut.writeTif( tosave, outname, self.epicure.epi_metadata["ScaleXY"], 'uint8', what="Skeleton" ) else: ut.writeTif(tosave, outname, self.epicure.epi_metadata["ScaleXY"], 'float32', what="Segmentation") @@ -375,7 +375,7 @@ class Outputing(QWidget): self.measure_one_frame( frame, properties, extra_properties, other_features, do_channels, int_feat, extra_prop, iframe, labgroups) self.table = pand.DataFrame.from_dict( self.table )[self.table.keys()] - if "intensity_junction_cytoplasm" in self.table.keys(): + if "intensity_junction_cytoplasm-0" in self.table.keys(): self.table = self.table.rename(columns={"intensity_junction_cytoplasm-0": "intensity_cytoplasm", "intensity_junction_cytoplasm-1":"intensity_junction"}) self.table_selection = self.selection_choices.index(self.output_mode.currentText()) self.viewer.window._status_bar._toggle_activity_dock(False) @@ -438,33 +438,56 @@ class Outputing(QWidget): ## add features of neighbors relationship with graph do_neighbor = "NbNeighbors" in other_features - do_border = "Border" in other_features - if do_neighbor or do_border: - graph = ut.connectivity_graph( img, distance=3 ) ## be sure that labels touch and get graph + if do_neighbor: + if frame is not None: + nimg = self.epicure.seg[frame] + else: + nimg = self.epicure.seg + #start_time = ut.start_time() + graph = ut.connectivity_graph( nimg, distance=3 ) ## be sure that labels touch and get graph adj_bg = list(graph.adj[0]) if 0 in graph.nodes else [] graph.remove_node(0) if 0 in graph.nodes else None - if do_neighbor: - if first: - self.table["NbNeighbors"] = [] - self.table["NbNeighbors"].extend( [-1]*ndata ) - if do_border: - if first: - self.table["Border"] = [] - self.table["Border"].extend( [-1]*ndata ) - - for label in graph.nodes: - rlabel = np.where( (frame_table["label"] == label) )[0] - if do_neighbor: + if first: + self.table["NbNeighbors"] = [] + self.table["NbNeighbors"].extend( [-1]*ndata ) + + for label in np.unique(frame_table["label"]): + if label in graph.nodes: + rlabel = np.where( (frame_table["label"] == label) )[0] nneighbor = len(graph.adj[label]) for ind in rlabel: self.table["NbNeighbors"][ind+nrows] = nneighbor - if do_border: - outer = int(label in adj_bg) + #ut.show_duration( start_time, "Neighborhoods measured" ) + + ## measure cells on boundary + if "Boundary" in other_features: + if frame is not None: + boundimg = self.epicure.seg[frame] + else: + boundimg = self.epicure.seg + bounds = ut.get_boundary_cells( boundimg ) + if first: + self.table["Boundary"] = [] + self.table["Boundary"].extend( [0]*ndata ) + for label in np.unique(frame_table["label"]): + if label in bounds: + rlabel = np.where( (frame_table["label"] == label) )[0] for ind in rlabel: - self.table["Border"][ind+nrows] = outer + self.table["Boundary"][ind+nrows] = 1 + ## measure cells on border + if "Border" in other_features: + bounds = ut.get_border_cells( img ) + if first: + self.table["Border"] = [] + self.table["Border"].extend( [0]*ndata ) + for label in bounds: + rlabel = np.where( (frame_table["label"] == label) )[0] + for ind in rlabel: + self.table["Border"][ind+nrows] = 1 + def selection_changed(self): if self.table_selection is None: return True @@ -539,7 +562,7 @@ class Outputing(QWidget): cell = self.seglayer.data==lab mapfeat[cell] = val ut.remove_layer(self.viewer, "Map_"+featname) - self.viewer.add_image(mapfeat, name="Map_"+featname) + self.viewer.add_image(mapfeat, name="Map_"+featname, scale=self.viewer.layers["Segmentation"].scale ) self.viewer.window._status_bar._toggle_activity_dock(False) def draw_orientation( self ): @@ -585,7 +608,7 @@ class Outputing(QWidget): oriens[ (frames, xas, yas) ] = 255 oriens[ (frames, xbs, ybs) ] = 255 - self.viewer.add_image( oriens, name="CellOrientation", blending="additive", opacity=1 ) + self.viewer.add_image( oriens, name="CellOrientation", blending="additive", opacity=1, scale=self.viewer.layers["Segmentation"].scale ) self.viewer.window._status_bar._toggle_activity_dock(False) ################### Export to other plugins @@ -712,14 +735,16 @@ class Outputing(QWidget): properties = ["label", "area", "centroid"] self.table = None - if type(track_ids) == int: + if type(track_ids) == np.ndarray or type(track_ids)==np.array: + track_ids = track_ids.tolist() + if not type(track_ids) == list: track_ids = [track_ids] for itrack, track_id in progress(enumerate(track_ids)): - track_frame = self.measure_one_track( track_id ) - if self.table is None: - self.table = pand.DataFrame([track_frame]) - else: - self.table = pand.concat([self.table, pand.DataFrame([track_frame])]) + track_frame = self.measure_one_track( track_id ) + if self.table is None: + self.table = pand.DataFrame([track_frame]) + else: + self.table = pand.concat([self.table, pand.DataFrame([track_frame])]) self.table_selection = self.selection_choices.index(self.output_mode.currentText()) self.viewer.window._status_bar._toggle_activity_dock(False) @@ -838,14 +863,27 @@ class CellFeatures(QWidget): self.features = {} self.chan_list = None - other_list = ["group", "NbNeighbors", "Border"] + other_list = ["group", "NbNeighbors", "Boundary", "Border"] feat_layout = self.add_feature_group( other_list, "other" ) layout.addLayout( feat_layout ) + sel_all_b = wid.add_button( "Select spatial features", lambda: self.select_all("other"), "Select all spatial features" ) + desel_all_b = wid.add_button( "Deselect spatial features", lambda: self.deselect_all("other"), "Deselect all spatial features" ) + sel_line_b = wid.hlayout() + sel_line_b.addWidget( sel_all_b ) + sel_line_b.addWidget( desel_all_b ) + layout.addLayout( sel_line_b ) + ## Add shape features shape_list = ["centroid", "area", "area_convex", "axis_major_length", "axis_minor_length", "feret_diameter_max", "equivalent_diameter_area", "eccentricity", "orientation", "perimeter", "solidity"] feat_layout = self.add_feature_group( shape_list, "prop" ) layout.addLayout( feat_layout ) + sel_all = wid.add_button( "Select morphology features", lambda: self.select_all("prop"), "Select all morphology features" ) + desel_all = wid.add_button( "Deselect morphology features", lambda: self.deselect_all("prop"), "Deselect all morphology features" ) + sel_line = wid.hlayout() + sel_line.addWidget( sel_all ) + sel_line.addWidget( desel_all ) + layout.addLayout( sel_line ) int_lab = wid.label_line( "Intensity features:") layout.addWidget( int_lab ) @@ -863,11 +901,39 @@ class CellFeatures(QWidget): self.chan_list.setSelectionMode(aiv.MultiSelection) self.chan_list.item(0).setSelected(True) layout.addWidget( self.chan_list ) + + sel_all_int = wid.add_button( "Select intensity features", lambda: self.select_all("intensity"), "Select all spatial features" ) + desel_all_int = wid.add_button( "Deselect intensity features", lambda: self.deselect_all("intensity"), "Deselect all spatial features" ) + sel_line_int = wid.hlayout() + sel_line_int.addWidget( sel_all_int ) + sel_line_int.addWidget( desel_all_int ) + layout.addLayout( sel_line_int ) bye = wid.add_button( "Ok", self.close, "Close the window" ) layout.addWidget( bye ) self.setLayout( layout ) + def select_all( self, feat ): + """ Select all features of type feat """ + if feat == "intensity": + self.select_all( "intensity_prop" ) + self.select_all( "intensity_extra" ) + return + for featy, feat_val in self.features.items(): + if feat_val[1] == feat: + feat_val[0].setChecked( True ) + + def deselect_all( self, feat ): + """ Deselect all features of type feat """ + if feat == "intensity": + self.deselect_all( "intensity_prop" ) + self.deselect_all( "intensity_extra" ) + return + for featy, feat_val in self.features.items(): + if feat_val[1] == feat: + feat_val[0].setChecked( False ) + + def add_feature_group( self, feat_list, feat_type ): """ Add features to the GUI """ layout = QVBoxLayout() @@ -1089,8 +1155,11 @@ class TemporalPlots(QWidget): layout.addWidget(self.plot_wid) ## save plot or save data of the plot line = wid.double_button( "Save plot image", self.save_plot_image, "Save the grapic in a PNG file", "Save plot data", self.save_plot_data, "Save the value used for the plot in .csv file" ) + self.by_label = wid.add_check( "Arranged data by label", False, None, "Save the data with one column by label" ) + line.addWidget( self.by_label ) layout.addLayout( line ) self.setLayout(layout) + self.resize(1000,800) def setTable(self, table): """ Data table to plot """ @@ -1122,7 +1191,6 @@ class TemporalPlots(QWidget): return if feat == "": return - self.ax.cla() if "group" in self.table: tab = list(zip(self.table["frame"], self.table[feat], self.table["label"], self.table["group"])) @@ -1130,37 +1198,44 @@ class TemporalPlots(QWidget): else: tab = list(zip(self.table["frame"], self.table[feat], self.table["label"])) self.df = pand.DataFrame( tab, columns=["frame", feat, "label"] ) - self.df.set_index('frame', inplace=True) - if self.smooth.isChecked(): - rollsize = 20 - ## average on a smaller scale if only few frames - if np.max( self.table["frame"] ) <= 20: - rollsize = 5 - if feat+"_smooth" in self.df.columns: - feat = feat+"_smooth" - else: - self.df[feat+"_smooth"] = self.df[feat].rolling(rollsize).mean() - feat = feat+"_smooth" + shape = "linear" + if self.smooth.isChecked(): + shape = "spline" if "group" in self.table and self.avg_group.isChecked(): self.dfmean = self.df.groupby(['group', 'frame'])[feat].mean().reset_index() - self.dfmean.set_index('frame', inplace=True) self.df.columns.name = 'group' - self.dfmean.groupby('group')[feat].plot(legend=False, ax=self.ax) - self.ax.legend(np.unique(self.dfmean['group'])) + self.fig = px.line( self.dfmean, x='frame', y=feat, color='group', labels={'frame': 'Time (frame)', feat: feat}, line_shape=shape, render_mode="svg" ) + else: + if len( np.unique(self.df["label"]) ) > 1000: + ut.show_warning( "Too many lines to plot; Using a random subset instead" ) + subset = sample( np.unique(self.df["label"]).tolist(), 1000) + subdf = self.df[self.df["label"].isin(subset)] + self.fig = px.line( subdf, x=subdf["frame"], y=feat, color='label', labels={'frame': 'Time (frame)', feat: feat}, line_shape = shape, render_mode="svg") + else: + self.fig = px.line(self.df, x=self.df["frame"], y=feat, color='label', labels={'frame': 'Time (frame)', feat: feat}, line_shape = shape, render_mode="svg") + + self.browser.setHtml( self.fig.to_html(include_plotlyjs='cdn')) + + def smooth_df( self, df ): + """ Smooth temporally the dataframe by label or by group """ + rollsize = 20 + ## average on a smaller scale if only few frames + if np.max( self.table["frame"] ) <= 20: + rollsize = 5 + if feat+"_smooth" in self.df.columns: + feat = feat+"_smooth" else: - self.df.groupby('label')[feat].plot(legend=False, ax=self.ax) - self.ax.set_ylabel(''+feat) - self.ax.set_xlabel('Time (frame)') - self.fig.canvas.draw_idle() - self.ymin, self.ymax = self.ax.get_ylim() + self.df[feat+"_smooth"] = self.df[feat].rolling(rollsize, center=True).mean() + print(self.df) + feat = feat+"_smooth" def save_plot_image( self ): """ Save current plot graphic to PNG image """ feat = self.feature_choice.currentText() outfile = self.epicure.outname()+"_plot_"+feat+".png" if self.fig is not None: - self.fig.savefig( outfile ) + self.fig.write_image( outfile ) if self.epicure.verbose > 0: ut.show_info("Measures saved in "+outfile) @@ -1170,29 +1245,30 @@ class TemporalPlots(QWidget): outfile = self.epicure.outname()+"_time_"+feat+".csv" if self.avg_group.isChecked(): data = self.dfmean.reset_index()[["frame", "group", feat]] - data[["frame", "group", feat]].to_csv( outfile, sep='\t', header=True, index=False ) + if self.by_label.isChecked(): + df = pand.pivot_table( data, columns="label", index="frame", values=feat ) + df.to_csv( outfile, sep='\t', header=True, index=True ) + else: + data[["frame", "group", feat]].to_csv( outfile, sep='\t', header=True, index=False ) else: data = self.df.reset_index()[["frame", "label", feat]] - data[["frame", "label", feat]].to_csv( outfile, sep='\t', header=True, index=False ) + if self.by_label.isChecked(): + df = pand.pivot_table( data, columns="label", index="frame", values=feat ) + df.to_csv( outfile, sep='\t', header=True, index=True ) + else: + data[["frame", "label", feat]].to_csv( outfile, sep='\t', header=True, index=False ) def move_framepos(self, frame): """ Move the vertical line showing the current frame position in the main window """ - if self.ax is not None: - if self.ymin is None: - self.ymin, self.ymax = self.ax.get_ylim() - if self.vline is not None: - self.vline.remove() - ymin = float(self.ymin*1.01) - ymax = float(self.ymax*0.99) - self.vline = self.ax.vlines(x=frame, ymin=ymin, ymax=ymax, ls=':', color="0.6") - self.fig.canvas.draw_idle() + return + #if self.fig is not None: + # self.fig.add_vline( x=frame, line_dash="dash", line_color="gray" ) + # self.browser.setHtml( self.fig.to_html(include_plotlyjs='cdn')) def create_plotwidget(self): """ Create plot window """ - mpl_widget = FigureCanvas( Figure(figsize=(6,6) ) ) - self.fig = mpl_widget.figure - self.ax = mpl_widget.figure.subplots() - return mpl_widget + self.browser = QWebEngineView(self) + return self.browser diff --git a/src/epicure/preferences.py b/src/epicure/preferences.py index f5751bca21d71d5484843199840a352a26beb7c5..cb7847e21a9f4902dc617ede66aaa905e008163e 100644 --- a/src/epicure/preferences.py +++ b/src/epicure/preferences.py @@ -155,7 +155,7 @@ class Preferences(): self.add_click_shortcut( "Labels", shortname="split accross", fulltext="drag-click in the cell to split into 2 cells ", button="Right", modifiers=["Control"] ) self.add_click_shortcut( "Labels", shortname="split draw", fulltext="drag-click draw a junction to split in 2 cells", button="Right", modifiers=["Alt"] ) self.add_click_shortcut( "Labels", shortname="redraw junction", fulltext="drag-click draw a junction to correct it", button="Left", modifiers=["Alt"] ) - self.add_key_shortcut( "Labels", shortname="draw junction mode", fulltext="<key shortcut> then Right click drawing connected junction(s) to create new cell", key="j" ) + self.add_key_shortcut( "Labels", shortname="draw junction mode", fulltext="Drawd junction(s) mode ON", key="j" ) self.add_click_shortcut( "Labels", shortname="drawing junction", fulltext="Draw junction mode ON. Drag-click draw a junction to create new cell(s)", button="Left", modifiers=["Control"] ) ## Seeds (manual segmentation) shortcuts diff --git a/src/epicure/start_epicuring.py b/src/epicure/start_epicuring.py index 0bd728386e871a9d3560c596bb5d3fc627406983..894b191a86533bd964fcd376a1ad2b417bd65af5 100644 --- a/src/epicure/start_epicuring.py +++ b/src/epicure/start_epicuring.py @@ -21,6 +21,7 @@ def start_epicure(): ncpus = int(multiprocessing.cpu_count()*0.5) def set_visibility(): + """ Handle the visibility of the advanced parameters """ get_files.output_dirname.visible = get_files.advanced_parameters.value get_files.show_other_chanels.visible = get_files.advanced_parameters.value get_files.process_frames_parallel.visible = get_files.advanced_parameters.value @@ -28,13 +29,17 @@ def start_epicure(): get_files.junction_half_thickness.visible = get_files.advanced_parameters.value get_files.verbose_level.visible = get_files.advanced_parameters.value get_files.allow_gaps.visible = get_files.advanced_parameters.value + get_files.show_scale_bar.visible = get_files.advanced_parameters.value + get_files.epithelial_cells.visible = get_files.advanced_parameters.value def load_movie(): + """ Load and display the selected movie """ start_time = ut.start_time() nonlocal caxis, cval image_file = get_files.image_file.value caxis, cval = Epic.load_movie(image_file) imgdir = ut.get_directory(image_file) + get_files.segmentation_file.visible = True get_files.segmentation_file.value = pathlib.Path(imgdir) labname = Epic.suggest_segfile( get_files.output_dirname.value ) Epic.set_names( get_files.output_dirname.value ) @@ -49,6 +54,11 @@ def start_epicure(): get_files.timeframe.value = Epic.epi_metadata["ScaleT"] get_files.unit_xy.value = Epic.epi_metadata["UnitXY"] get_files.unit_t.value = Epic.epi_metadata["UnitT"] + get_files.scale_xy.visible = True + get_files.unit_xy.visible = True + get_files.timeframe.visible = True + get_files.unit_t.visible = True + get_files.segment_with_epyseg.visible = True get_files.allow_gaps.value = bool(Epic.epi_metadata["Allow gaps"]) get_files.verbose_level.value = int(Epic.epi_metadata["Verbose"]) ut.show_duration(start_time, header="Movie loaded in ") @@ -67,9 +77,25 @@ def start_epicure(): show_others() ut.show_duration(start_time, header="Movie chanel loaded in ") + def launch_napari_epyseg(): + """ Open napari-epyseg plugin to segment the intensity channel movie """ + try: + import napari_epyseg + from napari_epyseg.start_epyseg import run_epyseg + except: + ut.show_error("This option requires the plugin napari-epyseg that is missing.\nInstall it and restart") + print("Running EpySeg with default parameters on the movie. To change the settings, use the napari-epyseg plugin outside of EpiCure or EpySeg module directly") + parameters = {"tile_width":32, "tile_height":32, "model":"epyseg default(v2)"} + segres = run_epyseg( Epic.img, parameters ) + segname = str(get_files.image_file.value)+"_epyseg.tif" + ut.writeTif( segres, segname, 1.0, "uint8", what="Epyseg results saved in " ) + get_files.segmentation_file.value = segname + get_files.segment_with_epyseg.visible = False + @magicgui(call_button="Start cure", junction_chanel={"widget_type": "Slider", "min":0, "max": 0}, + segment_with_epyseg = {"widget_type": "PushButton", "label": "Segment now with EpySeg"}, scale_xy = {"widget_type": "LiteralEvalLineEdit"}, timeframe = {"widget_type": "LiteralEvalLineEdit"}, junction_half_thickness={"widget_type": "LiteralEvalLineEdit"}, @@ -79,13 +105,16 @@ def start_epicure(): def get_files( image_file = pathlib.Path(cdir), junction_chanel = 0, segmentation_file = pathlib.Path(cdir), + segment_with_epyseg = False, scale_xy = 1, unit_xy = "um", timeframe = 1, unit_t = "min", advanced_parameters = False, show_other_chanels = True, + show_scale_bar = True, allow_gaps = True, + epithelial_cells = True, process_frames_parallel = False, nbparallel_threads = ncpus, junction_half_thickness = 1, @@ -101,16 +130,25 @@ def start_epicure(): Epic.set_verbose( verbose_level ) Epic.nparallel = nbparallel_threads #Epic.load_segmentation(segmentation_file) - Epic.set_thickness(junction_half_thickness) + Epic.set_thickness( junction_half_thickness ) Epic.set_scales(scale_xy, timeframe, unit_xy, unit_t) + Epic.set_scalebar( show_scale_bar ) Epic.set_gaps_option( allow_gaps ) + Epic.set_epithelia( epithelial_cells ) Epic.go_epicure(outdir, segmentation_file) set_visibility() + get_files.segmentation_file.visible = False + get_files.segment_with_epyseg.visible = False + get_files.scale_xy.visible = False + get_files.unit_xy.visible = False + get_files.timeframe.visible = False + get_files.unit_t.visible = False get_files.junction_chanel.visible = False get_files.advanced_parameters.clicked.connect(set_visibility) get_files.show_other_chanels.clicked.connect(show_others) get_files.image_file.changed.connect(load_movie) get_files.junction_chanel.changed.connect(set_chanel) + get_files.segment_with_epyseg.clicked.connect( launch_napari_epyseg ) return get_files diff --git a/src/epicure/tracking.py b/src/epicure/tracking.py index 360a14aebc9ecfda19ecfa1f8e440040bf5abdb2..a289ce060078b852ddac70e3782eba037671e52f 100644 --- a/src/epicure/tracking.py +++ b/src/epicure/tracking.py @@ -1,11 +1,10 @@ -from qtpy.QtWidgets import QVBoxLayout, QWidget +from qtpy.QtWidgets import QVBoxLayout, QWidget # type: ignore from epicure.laptrack_centroids import LaptrackCentroids from epicure.laptrack_overlaps import LaptrackOverlaps -from laptrack.data_conversion import convert_split_merge_df_to_napari_graph -from napari.utils import progress +from laptrack.data_conversion import convert_split_merge_df_to_napari_graph # type: ignore +from napari.utils import progress # type: ignore from skimage.transform import warp from skimage.registration import optical_flow_ilk -import time import pandas as pd import numpy as np from multiprocessing.pool import ThreadPool as Pool @@ -33,6 +32,9 @@ class Tracking(QWidget): self.track_update = wid.add_button( "Update tracks display", self.update_track_layer, "Update the Track layer with the changements made since the last update" ) layout.addWidget(self.track_update) + ## Correct track button + #track_reset = wid.add_button( "Correct track data", self.reset_tracks, "Correct the track data after some track was lost" ) + #layout.addWidget(track_reset) ## Method specific track_method, self.track_choice = wid.list_line( "Tracking method", "Choose the tracking method to use and display its parameter", func=None ) @@ -128,11 +130,14 @@ class Tracking(QWidget): ## plot tracks if len(track_table) > 0: + self.clear_graph() self.viewer.add_tracks( track_table, graph=self.graph, name=self.tracklayer_name, - properties = track_prop,) + properties = track_prop, + scale = self.viewer.layers["Segmentation"].scale, + ) self.viewer.layers[self.tracklayer_name].visible=True self.viewer.layers[self.tracklayer_name].color_by="track_id" ut.set_active_layer(self.viewer, "Segmentation") @@ -196,26 +201,38 @@ class Tracking(QWidget): track = self.get_track_data( track_id ) if track.shape[0] == 0: return features + track = track[track[:,1].argsort()] start = int(np.min(track[:,1])) end = int(np.max(track[:,1])) features["Label"] = track_id features["TrackDuration"] = end - start + 1 features["TrackStart"] = start features["TrackEnd"] = end + features["NbGaps"] = end - start + 1 - len(track) if (end-start) == 0: ## only one frame features["TotalDisplacement"] = None features["NetDisplacement"] = None features["Straightness"] = None + features["MeanVelocity"] = None else: features["TotalDisplacement"] = ut.total_distance( track[:,2:4] ) features["NetDisplacement"] = ut.net_distance( track[:,2:4] ) + features["MeanVelocity"] = np.mean( ut.velocities( track[:,1:4] ) ) if features["TotalDisplacement"] > 0: features["Straightness"] = features["NetDisplacement"]/features["TotalDisplacement"] else: features["Straightness"] = None return features + def measure_speed( self, track_id ): + """ Returns the velocities of the track """ + track = self.get_track_data( track_id ) + if track.shape[0] == 0: + return None + track = track[track[:,1].argsort()] + return ut.velocities( track[:,1:4] ) + def measure_features( self, track_id, features ): """ Measure features along all the track """ mask = self.epicure.get_mask( track_id ) @@ -254,7 +271,7 @@ class Tracking(QWidget): """ Get the dataframe of the labels in the segmented image """ resdf = None for iframe, frame in progress(enumerate(segimg)): - frame_table = self.get_one_frame( frame, iframe ) + frame_table = ut.labels_to_table( frame, iframe ) if resdf is None: resdf = pd.DataFrame(frame_table) else: @@ -331,6 +348,14 @@ class Tracking(QWidget): if isinstance( track_id, int ): return (np.flatnonzero( self.track_data[:,0] == track_id ) ) return (np.flatnonzero( np.isin( self.track_data[:,0], track_id ) ) ) + + def get_track_indexes_onframes( self, track_id, frames ): + """ Get indexes of track_id tracks position in the arrays """ + if isinstance( frames, int ): + frames = [frames] + if isinstance( track_id, int ): + return (np.flatnonzero( (self.track_data[:,0] == track_id) * np.isin( self.track_data[:,1], frames) ) ) + return (np.flatnonzero( np.isin( self.track_data[:,0], track_id ) * np.isin( self.track_data[:,1], frames) ) ) def get_track_indexes_from_frame(self, track_id, frame): """ Get indexes of track_id tracks position in the arrays from the given frame """ @@ -369,12 +394,14 @@ class Tracking(QWidget): """ Return the chosen column (frame, x, y, label) of track track_id """ indexes = self.get_track_indexes( track_id ) if column == "frame": - i = 1 + return self.track_data[indexes, 1] if column == "label": - i = 0 + return self.track_data[indexes, 0] if column == "pos": - i = [2,3] - track = self.track_data[indexes, i] + return self.track_data[indexes, 2:4] + if column == "fullpos": + return self.track_data[indexes, 1:4] + track = self.track_data[indexes] return track def get_frame_data( self, track_id, ind ): @@ -411,11 +438,19 @@ class Tracking(QWidget): return None return int( np.min(track[:,1]) ) + def is_in_frame( self, track_id, frame ): + """ Returns if track_id is present at given frame """ + track = self.get_track_data( track_id ) + if len(track) > 0: + return frame in track[:,1] + return False def get_last_frame(self, track_id): """ Returns last frame where track_id is present """ track = self.get_track_data( track_id ) - return int(np.max(track[:,1])) + if len(track) > 0: + return int(np.max(track[:,1])) + return None def get_extreme_frames(self, track_id): """ Returns the first and last frames where track_id is present """ @@ -440,6 +475,16 @@ class Tracking(QWidget): self.track_data[ind, 2] = cx self.track_data[ind, 3] = cy + def replace_on_frames( self, tids, new_tids, frames ): + """ Replace the id tid by new_tid in all given frames """ + ind = self.get_track_indexes_onframes( tids, frames ) + cur_track = np.copy(self.track_data[ind]) + new_ids = np.repeat(-1, len(ind)) + for tid, new_tid in zip(tids, new_tids): + self.update_graph_frames( tid, cur_track[cur_track[:,0]==tid,1] ) + new_ids[cur_track[:,0]==tid] = new_tid + self.track_data[ind, 0] = new_ids + def swap_frame_id(self, tid, otid, frame): """ Swap the ids of two tracks at frame """ ind = int(self.get_index(tid, frame)) @@ -461,6 +506,17 @@ class Tracking(QWidget): else: cur_cell = np.array( [[tid, frame, int(x), int(y)]] ) self.track_data = np.append(self.track_data, cur_cell, axis=0) + + def add_tracks_fromindices( self, indices, track_ids ): + """ Add tracks of given track ids from the indices""" + new_data = np.empty( (0,4), int ) + for tid in np.unique(track_ids): + keep = track_ids == tid + for frame in np.unique( indices[0][keep] ): + cent0 = np.mean( indices[1][keep] ) + cent1 = np.mean( indices[2][keep] ) + new_data = np.append( new_data, np.array([[tid, frame, int(cent0), int(cent1)]]), axis=0 ) + self.track_data = np.append( self.track_data, new_data, axis=0) def add_one_frame(self, track_ids, frame, refresh=True): """ Add one frame from track """ @@ -499,9 +555,34 @@ class Tracking(QWidget): centx, centy = self.track_data[ind, 2:4].astype(int).flatten() return self.epicure.seg[frame, centx, centy] + def clear_graph( self ): + """ Check the state of the graph and removes non existing keys or values """ + if self.graph is None: + return + keys = list(self.graph.keys()) + for key in keys: + if key not in self.track_data[:,0]: + del self.graph[key] + else: + vals = self.graph[key] + if isinstance(vals, list): + for val in vals: + if val not in self.track_data[:,0]: + del self.graph[key] + break + else: + if vals not in self.track_data[:,0]: + del self.graph[key] + + def update_graph_frames( self, track_id, frames ): + """ Update graph when one label was deleted at given frames """ + fframe = np.min(frames) + lframe = np.max(frames) + self.update_graph( track_id, fframe ) + self.update_graph( track_id, lframe ) + def update_graph(self, track_id, frame): """ Update graph if deleted label was linked at that frame, assume keys are unique """ - start_time = time.time() if self.graph is not None: ## handles current node is last of his branch parents = self.last_in_graph( track_id, frame ) @@ -517,12 +598,14 @@ class Tracking(QWidget): del self.graph[track_id] else: self.update_key( track_id, current_label ) - if self.epicure.verbose > 1: - ut.show_duration( start_time, header="Graph updated in " ) def update_child(self, parent, prev_key, new_key): """ Change the value of a key in the graph """ - self.graph[parent] = [new_key if val == prev_key else val for val in self.graph[parent]] + if isinstance(self.graph[parent], list): + self.graph[parent] = [new_key if val == prev_key else val for val in self.graph[parent]] + else: + if self.graph[parent] == prev_key: + self.graph[parent] = new_key def update_key(self, prev_key, new_key): """ Change the value of a key in the graph """ @@ -532,7 +615,7 @@ class Tracking(QWidget): """ Return if the current id is in the graph (as a parent, so in values) """ if self.graph is None: return False - return any( cur_id in vals for vals in self.graph.values() ) + return any( cur_id in vals if isinstance(vals, list) else cur_id in [vals] for vals in self.graph.values() ) def add_division( self, childa, childb, parent ): """ Add info of a division to the graph of divisions/merges """ @@ -544,13 +627,27 @@ class Tracking(QWidget): """ Remove a division event from the graph """ self.graph = {key: vals for key, vals in self.graph.items() if not (isinstance(vals, list) and vals[0] == parent) and vals != parent} - def last_in_graph(self, track_id, frame): + def last_in_graph(self, track_id, frame=None, check_last=True): """ Check if given label and frame is the last of a branch, in the graph """ - return [key for key, vals in self.graph.items() if track_id in (vals if isinstance(vals, list) else [vals]) and self.get_last_frame(track_id) == frame] + if check_last: + return [key for key, vals in self.graph.items() if track_id in (vals if isinstance(vals, list) else [vals]) and self.get_last_frame(track_id) == frame] + return [key for key, vals in self.graph.items() if track_id in (vals if isinstance(vals, list) else [vals])] - def first_in_graph(self, track_id, frame): + def first_in_graph(self, track_id, frame=None, check_first=True): """ Check if the given label and frame is the first in the branch so the node in the graph """ - return track_id in self.graph and self.get_first_frame(track_id) == frame + if check_first: + return track_id in self.graph and self.get_first_frame(track_id) == frame + return track_id in self.graph + + def remove_on_frames( self, track_ids, frames ): + """ Remove tracks with given id on specified frames """ + track_ids = track_ids.tolist() + if 0 in track_ids: + track_ids.remove(0) + inds = self.get_track_indexes_onframes( track_ids, frames ) + for tid in track_ids: + self.update_graph_frames( tid, frames ) + self.track_data = np.delete( self.track_data, inds, axis=0 ) def remove_tracks(self, track_ids): """ Remove track with given id """ @@ -565,20 +662,13 @@ class Tracking(QWidget): key: vals for key, vals in self.graph.items() if (key not in track_ids_set) and ( not any( val in track_ids_set for val in (vals if isinstance(vals, list) else [vals])) ) } - """ - for tid in track_ids: - self.graph.pop( tid, None ) - - keys = self.graph.keys() - for key in keys: - kgr = self.graph[ key ] - if not isinstance( kgr, list ): - kgr = [kgr] - for val in kgr: - if val in track_ids: - del self.graph[key] - continue - """ + + def is_single_parent( self, cur_id ): + """ Return if the current id is in the graph (as a single parent, not a merge) """ + if self.graph is None: + return False + return any( cur_id in [vals] if not isinstance(vals, list) else (cur_id in vals and len(vals)==1) for vals in self.graph.values() ) + def build_tracks(self, track_df): """ Create tracks from dataframe (after tracking) """ @@ -588,9 +678,11 @@ class Tracking(QWidget): def create_tracks(self): """ Create tracks from labels (without tracking) """ - track_table = np.empty( (0,4), int ) + #track_table = np.empty( (0,4), int ) + track_tables = [] for iframe, frame in progress(enumerate(self.epicure.seg), total=self.epicure.nframes): - track_table = self.get_one_frame(frame, iframe, track_table) + track_tables.append( ut.labels_to_table( frame, iframe ) ) + track_table = np.concatenate( track_tables, axis=0 ) return track_table, None # track_prop def add_track_features(self, labels): @@ -603,14 +695,6 @@ class Tracking(QWidget): nframes[ list(cur_track) ] = len(cur_track) return nframes - def get_one_frame(self, seg, frame, track_table): - """ Get the regionprops results into the dataframe """ - frame_table = ut.labels_table(seg, properties=self.properties) - labels = frame_table["label"] - frame_data = np.column_stack((labels.astype(int), np.full(len(labels), frame), frame_table["centroid-0"].astype(int), frame_table["centroid-1"].astype(int))) - track_table = np.vstack( (track_table, frame_data) ) - #return frame_data.astype(int) - return track_table ########################################## #### Tracking functions @@ -628,10 +712,10 @@ class Tracking(QWidget): """ Find in the first frame the parents of labels from second frame """ if self.track_choice.currentText() == "Laptrack-Centroids": - return self.laptrack_centroids_twoframes(labels, twoframes) + return self.laptrack_centroids_twoframes(labels, twoframes, loose=True) if self.track_choice.currentText() == "Laptrack-Overlaps": - return self.laptrack_overlaps_twoframes(labels, twoframes) + return self.laptrack_overlaps_twoframes(labels, twoframes, loose=True) def do_tracking(self): @@ -642,7 +726,7 @@ class Tracking(QWidget): else: start = 0 end = self.nframes-1 - start_time = time.time() + start_time = ut.start_time() self.viewer.window._status_bar._toggle_activity_dock(True) self.epicure.inspecting.reset_all_events() @@ -708,7 +792,7 @@ class Tracking(QWidget): def relabel_trackids(self, track_df, splitdf, mergedf): """ Change the trackids to take the first label of each track """ - start_time = time.time() + start_time = ut.start_time() new_trackids = track_df['track_id'].copy() new_splitdf = splitdf.copy() new_mergedf = mergedf.copy() @@ -914,17 +998,20 @@ class Tracking(QWidget): pen_layout.addLayout(solidCost) self.bpenalties.setLayout(pen_layout) - def laptrack_centroids_twoframes(self, labels, twoframes): - """ Perform tracking of two frames with current parameters """ + def laptrack_centroids_twoframes(self, labels, twoframes, loose=False): + """ Perform tracking of two frames with strict parameters """ laptrack = LaptrackCentroids(self, self.epicure) - laptrack.max_distance = float(self.max_dist.text()) + laptrack.max_distance = float(self.max_dist.text()) + if loose: + laptrack.max_distance = min(50, laptrack.max_distance) ## more probable to find a parent self.region_properties = ["label", "centroid"] - if self.check_penalties.isChecked(): - self.region_properties.append("area") - self.region_properties.append("solidity") - laptrack.penal_area = float(self.area_cost.text()) - laptrack.penal_solidity = float(self.solidity_cost.text()) - laptrack.set_region_properties(with_extra=self.check_penalties.isChecked()) + #if self.check_penalties.isChecked(): + # self.region_properties.append("area") + # self.region_properties.append("solidity") + # laptrack.penal_area = float(self.area_cost.text()) + # laptrack.penal_solidity = float(self.solidity_cost.text()) + #laptrack.set_region_properties(with_extra=self.check_penalties.isChecked()) + laptrack.set_region_properties(with_extra=False) df = self.twoframes_centroid(twoframes) if set(np.unique(df["label"])) == set(labels): @@ -1025,11 +1112,13 @@ class Tracking(QWidget): progress_bar.update(6) progress_bar.close() - def laptrack_overlaps_twoframes(self, labels, twoframes): - """ Perform tracking of two frames with current parameters """ + def laptrack_overlaps_twoframes(self, labels, twoframes, loose=False): + """ Perform tracking of two frames with strict parameters """ laptrack = LaptrackOverlaps(self, self.epicure) - miniou = max( float(self.min_iou.text()), 1.0 ) + miniou = min( float(self.min_iou.text()), 0.9999 ) ## ensure that miniou is < 1 laptrack.cost_cutoff = 1.0 - miniou + if loose: + laptrack.cost_cutoff = 0.95 ## more probable to find a parent/child self.region_properties = ["label", "centroid"] laptrack.splitting_cost = False ## disable splitting option diff --git a/src/tests/test_suspects.py b/src/tests/test_suspects.py index 267cacfa40eec31807742ea208dbf73337451f7b..03491169e729999834cfd3e468b208098116f1d5 100644 --- a/src/tests/test_suspects.py +++ b/src/tests/test_suspects.py @@ -76,4 +76,23 @@ def test_suspect_track(): susp.inspect_tracks() assert susp.nb_events() > nmin +def test_boundaries(): + test_img = os.path.join(".", "data_test", "003_crop.tif") + test_seg = os.path.join(".", "data_test", "003_crop_epyseg.tif") + + ## load and initialize + epic = epi.EpiCure() + epic.load_movie(test_img) + epic.go_epicure("test_epics", test_seg) + + ## check that doesn't find any boundary cells (touching background) + epic.inspecting.get_boundaries_cells() + assert len(epic.inspecting.boundary_cells[0]) == 0 + + ## remove the border cells, so now should find boundary cells + epic.editing.remove_border() + epic.inspecting.get_boundaries_cells() + assert len(epic.inspecting.boundary_cells[0]) > 20 + #test_suspect_track() +#test_boundaries()