diff --git a/Tests/test_vertexModel.py b/Tests/test_vertexModel.py index 24fe329..662911b 100644 --- a/Tests/test_vertexModel.py +++ b/Tests/test_vertexModel.py @@ -12,9 +12,9 @@ from src.pyVertexModel.algorithm.vertexModel import create_tetrahedra from src.pyVertexModel.algorithm.vertexModelBubbles import build_topo, SeedWithBoundingBox, generate_first_ghost_nodes, \ delaunay_compute_entities, VertexModelBubbles -from src.pyVertexModel.algorithm.vertexModelVoronoiFromTimeImage import build_triplets_of_neighs, calculate_neighbours, \ - VertexModelVoronoiFromTimeImage, add_tetrahedral_intercalations, build_2d_voronoi_from_image, \ - populate_vertices_info, calculate_vertices, get_four_fold_vertices, divide_quartets_neighbours, process_image +from src.pyVertexModel.algorithm.vertexModelVoronoiFromTimeImage import build_triplets_of_neighs, \ + VertexModelVoronoiFromTimeImage, add_tetrahedral_intercalations, \ + get_four_fold_vertices, divide_quartets_neighbours, process_image from src.pyVertexModel.geometry.degreesOfFreedom import DegreesOfFreedom from src.pyVertexModel.util.utils import save_backup_vars @@ -496,3 +496,84 @@ def test_initialize_voronoi_from_time_image(self): assert_array1D(g_test, mat_info['g']) assert_matrix(K_test, mat_info['K']) + def test_initialize_cells_with_numpy_array(self): + """ + Test that initialize_cells can accept a numpy array directly. + :return: + """ + import scipy.io + from src.pyVertexModel.parameters.set import Set + + # Load an existing image as a numpy array + mat_data = scipy.io.loadmat('resources/LblImg_imageSequence.mat') + img_array = mat_data['imgStackLabelled'] + + # Create a simple 2D labeled image for testing + # Using a subset to make it faster + img_2d = img_array[:, :, 0] + + # Create settings + set_test = Set(set_option='voronoi_from_image') + set_test.TotalCells = 50 # Use fewer cells for faster testing + set_test.CellHeight = 10 + + # Test with numpy array input + vModel_test = VertexModelVoronoiFromTimeImage(set_option='voronoi_from_image', set_test=set_test, + create_output_folder=False) + vModel_test.initialize_cells(img_2d) + + # Verify that the geometry was created + assert vModel_test.geo is not None, "Geometry should be initialized" + assert vModel_test.geo.nCells > 0, "Should have cells" + assert len(vModel_test.geo.Cells) > 0, "Should have Cell objects" + + def test_process_image_with_numpy_array(self): + """ + Test that process_image can handle numpy array input. + :return: + """ + # Create a simple labeled image + test_img = np.zeros((100, 100), dtype=np.uint16) + # Create some labeled regions + test_img[10:30, 10:30] = 1 + test_img[40:60, 40:60] = 2 + test_img[70:90, 70:90] = 3 + + # Test process_image with numpy array + img2d, img_stack = process_image(test_img) + + # Verify the output + assert img2d is not None, "2D image should be returned" + assert img_stack is not None, "Image stack should be returned" + assert img2d.shape == test_img.shape, "2D image should have same shape as input" + + def test_initialize_with_numpy_array(self): + """ + Test that initialize method can accept a numpy array. + :return: + """ + import scipy.io + from src.pyVertexModel.parameters.set import Set + + # Load an existing image as a numpy array + mat_data = scipy.io.loadmat('resources/LblImg_imageSequence.mat') + img_array = mat_data['imgStackLabelled'] + + # Use a 2D slice for faster testing + img_2d = img_array[:, :, 0] + + # Create settings + set_test = Set(set_option='voronoi_from_image') + set_test.TotalCells = 50 + set_test.CellHeight = 10 + + # Test initialize with numpy array input + vModel_test = VertexModelVoronoiFromTimeImage(set_option='voronoi_from_image', set_test=set_test, + create_output_folder=False) + vModel_test.initialize(img_2d) + + # Verify that the geometry was created + assert vModel_test.geo is not None, "Geometry should be initialized" + assert vModel_test.geo.nCells > 0, "Should have cells" + + diff --git a/src/pyVertexModel/algorithm/vertexModel.py b/src/pyVertexModel/algorithm/vertexModel.py index 0eb1a43..176c616 100644 --- a/src/pyVertexModel/algorithm/vertexModel.py +++ b/src/pyVertexModel/algorithm/vertexModel.py @@ -251,9 +251,11 @@ def __init__(self, set_option='wing_disc', c_set=None, create_output_folder=True self.tr = 0 self.numStep = 1 - def initialize(self): + def initialize(self, img_input=None): """ Initialize the geometry and the topology of the model. + :param img_input: Optional. Either a filename (str) or a numpy array containing the image. + If None, uses the filename from settings. """ filename = os.path.join(PROJECT_DIRECTORY, self.set.initial_filename_state) @@ -285,7 +287,10 @@ def initialize(self): self.geo = Geo(mat_info['Geo']) self.geo.update_measures() else: - self.initialize_cells(filename) + if img_input is None: + self.initialize_cells(filename) + else: + self.initialize_cells(img_input) # Resize the geometry to a given cell volume average self.geo.resize_tissue() diff --git a/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py b/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py index f8dfe6d..750c6f5 100644 --- a/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py +++ b/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py @@ -197,81 +197,157 @@ def divide_quartets_neighbours(img_neighbours_all, labelled_img, quartets): return img_neighbours_all -def process_image(img_filename="src/pyVertexModel/resources/LblImg_imageSequence.tif", redo=False): +def process_image(img_input="src/pyVertexModel/resources/LblImg_imageSequence.tif", redo=False): """ Process the image and return the 2D labeled image and the 3D labeled image. - :param redo: - :param img_filename: - :return: + :param redo: Whether to redo the processing even if cached version exists + :param img_input: Either a filename (str) or a numpy array containing the image + :return: img2DLabelled, imgStackLabelled """ - # Load the tif file from resources if exists - if os.path.exists(img_filename): - if img_filename.endswith('.tif'): - if os.path.exists(img_filename.replace('.tif', '.xz')) and not redo: - imgStackLabelled = pickle.load(lzma.open(img_filename.replace('.tif', '.xz'), "rb")) - - imgStackLabelled = imgStackLabelled['imgStackLabelled'] - if imgStackLabelled.ndim == 3: - img2DLabelled = imgStackLabelled[0, :, :] - else: - img2DLabelled = imgStackLabelled + # Check if input is a numpy array or a filename + if isinstance(img_input, np.ndarray): + # Input is a numpy array, process it directly + imgStackLabelled = img_input.copy() + + # Reordering cells based on the centre of the image + if imgStackLabelled.ndim == 3: + img2DLabelled = imgStackLabelled[0, :, :] + else: + img2DLabelled = imgStackLabelled + + # Check if the image is already segmented (has labels > 1) or needs segmentation + if np.max(imgStackLabelled) <= 1: + # Image is binary (0 and 1 only), needs segmentation + # Invert so cells are 1 and background is 0, then label + imgStackLabelled = (imgStackLabelled > 0).astype(np.uint16) + from scipy.ndimage import label as scipy_label + imgStackLabelled, num_features = scipy_label(imgStackLabelled) + if imgStackLabelled.ndim == 3: + img2DLabelled = imgStackLabelled[0, :, :] else: - imgStackLabelled = io.imread(img_filename) - - # Reordering cells based on the centre of the image - if imgStackLabelled.ndim == 3: - img2DLabelled = imgStackLabelled[0, :, :] + img2DLabelled = imgStackLabelled + + # For already segmented images (max > 1), we still need to label background regions + # This is consistent with the file-loading behavior + # Use appropriate structure based on dimensionality + if imgStackLabelled.ndim == 3: + structure_element = np.array([ + [[0, 0, 0], [0, 1, 0], [0, 0, 0]], + [[0, 1, 0], [1, 1, 1], [0, 1, 0]], + [[0, 0, 0], [0, 1, 0], [0, 0, 0]] + ]) + else: + structure_element = [[0, 1, 0], [1, 1, 1], [0, 1, 0]] + + imgStackLabelled, num_features = label(imgStackLabelled == 0, + structure=structure_element) + + props = regionprops_table(imgStackLabelled, properties=('centroid', 'label',), ) + + # The centroids are now stored in 'props' as separate arrays 'centroid-0', 'centroid-1', etc. + centroids = np.column_stack([props['centroid-0'], props['centroid-1']]) + centre_of_image = np.array([img2DLabelled.shape[0] / 2, img2DLabelled.shape[1] / 2]) + + # Sorting cells based on distance to the middle of the image + distanceToMiddle = cdist([centre_of_image], centroids) + distanceToMiddle = distanceToMiddle[0] + sortedId = np.argsort(distanceToMiddle) + sorted_ids = np.array(props['label'])[sortedId] + + oldImgStackLabelled = copy.deepcopy(imgStackLabelled) + newCont = 1 + for numCell in sorted_ids: + if numCell != 0: + imgStackLabelled[oldImgStackLabelled == numCell] = newCont + newCont += 1 + + # Remaining cells that are not in the image + for numCell in np.arange(newCont, np.max(img2DLabelled) + 1): + imgStackLabelled[oldImgStackLabelled == numCell] = newCont + newCont += 1 + + if imgStackLabelled.ndim == 3: + img2DLabelled = imgStackLabelled[0, :, :] + else: + img2DLabelled = imgStackLabelled + + if imgStackLabelled.ndim == 3: + # Transpose the image stack to have the shape (Z, Y, X) + imgStackLabelled = np.transpose(imgStackLabelled, (2, 0, 1)) + + elif isinstance(img_input, str): + # Input is a filename, load from file + img_filename = img_input + # Load the tif file from resources if exists + if os.path.exists(img_filename): + if img_filename.endswith('.tif'): + if os.path.exists(img_filename.replace('.tif', '.xz')) and not redo: + imgStackLabelled = pickle.load(lzma.open(img_filename.replace('.tif', '.xz'), "rb")) + + imgStackLabelled = imgStackLabelled['imgStackLabelled'] + if imgStackLabelled.ndim == 3: + img2DLabelled = imgStackLabelled[0, :, :] + else: + img2DLabelled = imgStackLabelled else: - img2DLabelled = imgStackLabelled - - imgStackLabelled, num_features = label(imgStackLabelled == 0, - structure=[[0, 1, 0], [1, 1, 1], [0, 1, 0]]) - props = regionprops_table(imgStackLabelled, properties=('centroid', 'label',), ) - - # The centroids are now stored in 'props' as separate arrays 'centroid-0', 'centroid-1', etc. - centroids = np.column_stack([props['centroid-0'], props['centroid-1']]) - centre_of_image = np.array([img2DLabelled.shape[0] / 2, img2DLabelled.shape[1] / 2]) - - # Sorting cells based on distance to the middle of the image - distanceToMiddle = cdist([centre_of_image], centroids) - distanceToMiddle = distanceToMiddle[0] - sortedId = np.argsort(distanceToMiddle) - sorted_ids = np.array(props['label'])[sortedId] - - oldImgStackLabelled = copy.deepcopy(imgStackLabelled) - # imgStackLabelled = np.zeros_like(imgStackLabelled) - newCont = 1 - for numCell in sorted_ids: - if numCell != 0: + imgStackLabelled = io.imread(img_filename) + + # Reordering cells based on the centre of the image + if imgStackLabelled.ndim == 3: + img2DLabelled = imgStackLabelled[0, :, :] + else: + img2DLabelled = imgStackLabelled + + imgStackLabelled, num_features = label(imgStackLabelled == 0, + structure=[[0, 1, 0], [1, 1, 1], [0, 1, 0]]) + props = regionprops_table(imgStackLabelled, properties=('centroid', 'label',), ) + + # The centroids are now stored in 'props' as separate arrays 'centroid-0', 'centroid-1', etc. + centroids = np.column_stack([props['centroid-0'], props['centroid-1']]) + centre_of_image = np.array([img2DLabelled.shape[0] / 2, img2DLabelled.shape[1] / 2]) + + # Sorting cells based on distance to the middle of the image + distanceToMiddle = cdist([centre_of_image], centroids) + distanceToMiddle = distanceToMiddle[0] + sortedId = np.argsort(distanceToMiddle) + sorted_ids = np.array(props['label'])[sortedId] + + oldImgStackLabelled = copy.deepcopy(imgStackLabelled) + # imgStackLabelled = np.zeros_like(imgStackLabelled) + newCont = 1 + for numCell in sorted_ids: + if numCell != 0: + imgStackLabelled[oldImgStackLabelled == numCell] = newCont + newCont += 1 + + # Remaining cells that are not in the image + for numCell in np.arange(newCont, np.max(img2DLabelled) + 1): imgStackLabelled[oldImgStackLabelled == numCell] = newCont newCont += 1 - # Remaining cells that are not in the image - for numCell in np.arange(newCont, np.max(img2DLabelled) + 1): - imgStackLabelled[oldImgStackLabelled == numCell] = newCont - newCont += 1 + if imgStackLabelled.ndim == 3: + img2DLabelled = imgStackLabelled[0, :, :] + else: + img2DLabelled = imgStackLabelled + + # Save the labeled image stack as tif + io.imsave(img_filename.replace('.tif', '_labelled.tif'), imgStackLabelled.astype(np.uint16)) + save_variables({'imgStackLabelled': imgStackLabelled}, + img_filename.replace('.tif', '.xz')) if imgStackLabelled.ndim == 3: - img2DLabelled = imgStackLabelled[0, :, :] + # Transpose the image stack to have the shape (Z, Y, X) + imgStackLabelled = np.transpose(imgStackLabelled, (2, 0, 1)) + elif img_filename.endswith('.mat'): + imgStackLabelled = scipy.io.loadmat(img_filename)['imgStackLabelled'] + if imgStackLabelled.ndim == 3: + img2DLabelled = imgStackLabelled[:, :, 0] else: img2DLabelled = imgStackLabelled - - # Save the labeled image stack as tif - io.imsave(img_filename.replace('.tif', '_labelled.tif'), imgStackLabelled.astype(np.uint16)) - - save_variables({'imgStackLabelled': imgStackLabelled}, - img_filename.replace('.tif', '.xz')) - if imgStackLabelled.ndim == 3: - # Transpose the image stack to have the shape (Z, Y, X) - imgStackLabelled = np.transpose(imgStackLabelled, (2, 0, 1)) - elif img_filename.endswith('.mat'): - imgStackLabelled = scipy.io.loadmat(img_filename)['imgStackLabelled'] - if imgStackLabelled.ndim == 3: - img2DLabelled = imgStackLabelled[:, :, 0] - else: - img2DLabelled = imgStackLabelled + else: + raise ValueError('Image file not found %s' % img_filename) else: - raise ValueError('Image file not found %s' % img_filename) + raise TypeError(f'img_input must be a string (filename) or numpy.ndarray, got {type(img_input)}') return img2DLabelled, imgStackLabelled @@ -319,13 +395,25 @@ def __init__(self, set_option='wing_disc', set_test=None, update_derived_paramet create_output_folder=create_output_folder) self.dilated_cells = None - def initialize_cells(self, filename): + def initialize_cells(self, img_input): """ Initialize the cells from the image. - :param filename: + :param img_input: Either a filename (str) or a numpy array containing the image :return: """ - output_filename = filename.replace('.tif', f'_{self.set.TotalCells}cells_{self.set.CellHeight}.pkl') + # Generate a unique output filename based on input type + if isinstance(img_input, str): + # Input is a filename + output_filename = img_input.replace('.tif', f'_{self.set.TotalCells}cells_{self.set.CellHeight}.pkl') + elif isinstance(img_input, np.ndarray): + # Input is a numpy array - create a generic filename based on array shape and settings + output_filename = f'vertex_model_array_{img_input.shape}_{self.set.TotalCells}cells_{self.set.CellHeight}.pkl' + # Save it in the output folder if available, otherwise in current directory + if hasattr(self.set, 'OutputFolder') and self.set.OutputFolder: + output_filename = os.path.join(self.set.OutputFolder, output_filename) + else: + raise TypeError(f'img_input must be a string (filename) or numpy.ndarray, got {type(img_input)}') + if exists(output_filename): # Check date of the output_filename and if it is older than 1 day from today, redo the file # if os.path.getmtime(output_filename) < (time.time() - 24 * 60 * 60): @@ -337,12 +425,18 @@ def initialize_cells(self, filename): return # Load the image and obtain the initial X and tetrahedra - Twg, X = self.obtain_initial_x_and_tetrahedra() + Twg, X = self.obtain_initial_x_and_tetrahedra(img_input) # Build cells self.geo.build_cells(self.set, X, Twg) - # Save screenshot of the initial state - image_file = '/'+ os.path.join(*filename.split('/')[:-1]) - screenshot_(self.geo, self.set, 0, output_filename.split('/')[-1], image_file) + + # Save screenshot of the initial state (only if we have a filename) + if isinstance(img_input, str): + image_file = '/'+ os.path.join(*img_input.split('/')[:-1]) + screenshot_(self.geo, self.set, 0, output_filename.split('/')[-1], image_file) + else: + # For numpy array input, try to save screenshot in output folder if available + if hasattr(self.set, 'OutputFolder') and self.set.OutputFolder: + screenshot_(self.geo, self.set, 0, output_filename.split('/')[-1], self.set.OutputFolder) # Save state with filename using the number of cells save_state(self.geo, output_filename) @@ -532,17 +626,18 @@ def calculate_vertices(self, labelled_img, neighbours): return vertices_info - def obtain_initial_x_and_tetrahedra(self, img_filename=None): + def obtain_initial_x_and_tetrahedra(self, img_input=None): """ Obtain the initial X and tetrahedra for the model. - :return: + :param img_input: Either a filename (str) or a numpy array containing the image. If None, uses default from settings. + :return: Twg, X """ self.geo = Geo() - if img_filename is None: - img_filename = PROJECT_DIRECTORY + '/' + self.set.initial_filename_state + if img_input is None: + img_input = PROJECT_DIRECTORY + '/' + self.set.initial_filename_state selectedPlanes = [0, 0] - img2DLabelled, imgStackLabelled = process_image(img_filename) + img2DLabelled, imgStackLabelled = process_image(img_input) # Building the topology of each plane trianglesConnectivity = {} diff --git a/src/pyVertexModel/geometry/geo.py b/src/pyVertexModel/geometry/geo.py index b00d4b8..801d2d1 100644 --- a/src/pyVertexModel/geometry/geo.py +++ b/src/pyVertexModel/geometry/geo.py @@ -324,11 +324,12 @@ def build_cells(self, c_set, X, twg): self.Cells[c].Area = self.Cells[c].compute_area() self.Cells[c].Vol = self.Cells[c].compute_volume() + # Unique Ids for each point (vertex, node or face center) used in K + self.build_global_ids() + # Initialize reference values self.init_reference_cell_values(c_set) - # Unique Ids for each point (vertex, node or face center) used in K - self.build_global_ids() if c_set.Substrate == 1: for c, c_cell in enumerate(self.Cells): @@ -356,6 +357,7 @@ def init_reference_cell_values(self, c_set): :return: None """ logger.info('Initializing reference cell values') + # Assemble nodes from all cells that are not None self.AssembleNodes = [i for i, cell in enumerate(self.Cells) if cell.AliveStatus is not None]