11363842044

Committed 16 Oct 2024 10:30AM UTC coverage: 0.794% (-8.0%) from 8.812%

Build # 11363842044

Build Type

push

github

Committed by

Pablo1990

Commit Message

Big refactor on interface type

might affect performance, but minimise errors

Run Details

0 of 1002 branches covered (0.0%)

Branch coverage included in aggregate %.

1 of 64 new or added lines in 10 files covered. (1.56%)

524 existing lines in 33 files now uncovered.

51 of 5423 relevant lines covered (0.94%)

0.01 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/src/pyVertexModel/algorithm/vertexModelVoronoiFromTimeImage.py

UNCOV 1	import copy	×
UNCOV 2	import logging	×
UNCOV 3	import lzma	×
UNCOV 4	import os	×
UNCOV 5	import pickle	×
UNCOV 6	from itertools import combinations	×
7
UNCOV 8	import numpy as np	×
UNCOV 9	import scipy	×
UNCOV 10	from numpy import arange	×
UNCOV 11	from scipy.spatial.distance import squareform, pdist, cdist	×
UNCOV 12	from skimage import io	×
UNCOV 13	from skimage.measure import regionprops_table	×
UNCOV 14	from skimage.morphology import dilation, disk, square	×
UNCOV 15	from skimage.segmentation import find_boundaries	×
16
UNCOV 17	from src import PROJECT_DIRECTORY	×
UNCOV 18	from src.pyVertexModel.algorithm.vertexModel import VertexModel, generate_tetrahedra_from_information, \	×
19	calculate_cell_height_on_model
UNCOV 20	from src.pyVertexModel.geometry.geo import Geo	×
UNCOV 21	from src.pyVertexModel.util.utils import ismember_rows, save_variables, load_state	×
22
23
UNCOV 24	def calculate_neighbours(labelled_img, ratio_strel):	×
25	"""
26	Calculate the neighbours of each cell
27	:param labelled_img:
28	:param ratio_strel:
29	:return:
30	"""
31	se = disk(ratio_strel)	×
32
33	cells = np.sort(np.unique(labelled_img))	×
34	if np.sum(labelled_img == 0) > 0:	×
35	# Deleting cell 0 from range
36	cells = cells[1:]	×
37
38	img_neighbours = [None] * (np.max(cells) + 1)	×
39
40	for idx, cell in enumerate(cells, start=1):	×
41	BW = find_boundaries(labelled_img == cell, mode='inner')	×
42	BW_dilate = dilation(BW, se)	×
43	neighs = np.unique(labelled_img[BW_dilate == 1])	×
44	img_neighbours[idx] = neighs[(neighs != 0) & (neighs != cell)]	×
45
46	return img_neighbours	×
47
48
UNCOV 49	def build_quartets_of_neighs_2d(neighbours):	×
50	"""
51	Build quartets of neighboring cells.
52
53	:param neighbours: A list of lists where each sublist represents the neighbors of a cell.
54	:return: A 2D numpy array where each row represents a quartet of neighboring cells.
55	"""
56	quartets_of_neighs = []	×
57
58	for n_cell, neigh_cell in enumerate(neighbours):	×
59	if neigh_cell is not None:	×
60	intercept_cells = [None] * len(neigh_cell)	×
61
62	for cell_j in range(len(neigh_cell)):	×
63	if neighbours[neigh_cell[cell_j]] is not None and neigh_cell is not None:	×
64	common_cells = list(set(neigh_cell).intersection(neighbours[neigh_cell[cell_j]]))	×
65	if len(common_cells) > 2:	×
66	intercept_cells[cell_j] = common_cells + [neigh_cell[cell_j], n_cell]	×
67
68	intercept_cells = [cell for cell in intercept_cells if cell is not None]	×
69
70	if intercept_cells:	×
71	for index_a in range(len(intercept_cells) - 1):	×
72	for index_b in range(index_a + 1, len(intercept_cells)):	×
73	intersection_cells = list(set(intercept_cells[index_a]).intersection(intercept_cells[index_b]))	×
74	if len(intersection_cells) >= 4:	×
75	quartets_of_neighs.extend(list(combinations(intersection_cells, 4)))	×
76
77	if len(quartets_of_neighs) > 0:	×
78	quartets_of_neighs = np.unique(np.sort(quartets_of_neighs, axis=1), axis=0)	×
79
80	return quartets_of_neighs	×
81
82
UNCOV 83	def get_four_fold_vertices(img_neighbours):	×
84	"""
85	Get the four-fold vertices of the cells.
86	:param img_neighbours:
87	:return:
88	"""
89	quartets = build_quartets_of_neighs_2d(img_neighbours)	×
90	if len(quartets) == 0:	×
91	return None, 0	×
92
93	percQuartets = quartets.shape[0] / len(img_neighbours)	×
94
95	return quartets, percQuartets	×
96
97
UNCOV 98	def build_triplets_of_neighs(neighbours):	×
99	"""
100	Build triplets of neighboring cells.
101
102	:param neighbours: A list of lists where each sublist represents the neighbors of a cell.
103	:return: A 2D numpy array where each row represents a triplet of neighboring cells.
104	"""
105	triplets_of_neighs = []	×
106
107	for i, neigh_i in enumerate(neighbours):	×
108	if neigh_i is not None:	×
109	for j in neigh_i:	×
110	if j > i:	×
111	neigh_j = neighbours[j]	×
112	if neigh_j is not None:	×
113	for k in neigh_j:	×
114	if k > j and neighbours[k] is not None:	×
115	common_cell = {i}.intersection(neigh_j, neighbours[k])	×
116	if common_cell:	×
117	triangle_seed = sorted([i, j, k])	×
118	triplets_of_neighs.append(triangle_seed)	×
119
120	if len(triplets_of_neighs) > 0:	×
121	triplets_of_neighs = np.unique(np.sort(triplets_of_neighs, axis=1), axis=0)	×
122
123	return triplets_of_neighs	×
124
125
UNCOV 126	def calculate_vertices(labelled_img, neighbours, ratio):	×
127	"""
128	Calculate the vertices for each cell in a labeled image.
129
130	:param labelled_img: A 2D array representing a labeled image.
131	:param neighbours: A list of lists where each sublist represents the neighbors of a cell.
132	:param ratio: The radius of the disk used for morphological dilation.
133	:return: A dictionary containing the location of each vertex and the cells connected to each vertex.
134	"""
135	se = disk(ratio)	×
136	neighbours_vertices = build_triplets_of_neighs(neighbours)	×
137	# Initialize vertices
138	vertices = [None] * len(neighbours_vertices)	×
139
140	# Calculate the perimeter of each cell for efficiency
141	dilated_cells = [None] * (np.max(labelled_img) + 1)	×
142
143	for i in arange(1, np.max(labelled_img) + 1):	×
144	BW = np.zeros_like(labelled_img)	×
145	BW[labelled_img == i] = 1	×
146	BW_dilated = dilation(find_boundaries(BW, mode='inner'), se)	×
147	dilated_cells[i] = BW_dilated	×
148
149	# The overlap between cells in the labeled image will be the vertices
150	border_img = np.zeros_like(labelled_img)	×
151	border_img[labelled_img > -1] = 1	×
152
153	for num_triplet in range(len(neighbours_vertices)):	×
154	BW1_dilate = dilated_cells[neighbours_vertices[num_triplet][0]]	×
155	BW2_dilate = dilated_cells[neighbours_vertices[num_triplet][1]]	×
156	BW3_dilate = dilated_cells[neighbours_vertices[num_triplet][2]]	×
157
158	row, col = np.where((BW1_dilate * BW2_dilate * BW3_dilate * border_img) == 1)	×
159
160	if len(row) > 1:	×
161	if round(np.mean(col)) not in col:	×
162	vertices[num_triplet] = [round(np.mean([row[col > np.mean(col)], col[col > np.mean(col)]]))]	×
163	vertices.append([round(np.mean([row[col < np.mean(col)], col[col < np.mean(col)]]))])	×
164	else:
165	vertices[num_triplet] = [round(np.mean(row)), round(np.mean(col))]	×
166	elif len(row) == 0:	×
167	vertices[num_triplet] = [None, None]	×
168	else:
169	vertices[num_triplet] = [[row, col]]	×
170
171	# Store vertices and remove artifacts
172	vertices_info = {'location': vertices, 'connectedCells': neighbours_vertices}	×
173
174	not_empty_cells = [v[0] is not None for v in vertices_info['location']]	×
175	if len(vertices_info['location'][0]) == 2:	×
176	vertices_info['location'] = [vertices_info['location'][i] for i in range(len(not_empty_cells)) if	×
177	not_empty_cells[i]]
178	vertices_info['connectedCells'] = [vertices_info['connectedCells'][i] for i in range(len(not_empty_cells)) if	×
179	not_empty_cells[i]]
180	else:
181	vertices_info['location'] = [vertices_info['location'][i] for i in range(len(not_empty_cells)) if	×
182	not_empty_cells[i]]
183	vertices_info['connectedCells'] = [vertices_info['connectedCells'][i] for i in range(len(not_empty_cells)) if	×
184	not_empty_cells[i]]
185
186	return vertices_info	×
187
188
UNCOV 189	def boundary_of_cell(vertices_of_cell, neighbours=None):	×
190	"""
191	Determine the order of vertices that form the boundary of a cell.
192
193	:param vertices_of_cell: A 2D array where each row represents the coordinates of a vertex.
194	:param neighbours: A 2D array where each row represents a pair of neighboring vertices.
195	:return: A 1D array representing the order of vertices.
196	"""
197	# If neighbours are provided, try to order the vertices based on their neighbors
198	if neighbours is not None:	×
199	initial_neighbours = neighbours	×
200	neighbours_order = neighbours[0]	×
201	next_neighbour = neighbours[0][1]	×
202	next_neighbour_prev = next_neighbour	×
203	neighbours = np.delete(neighbours, 0, axis=0)	×
204
205	# Loop until all neighbours are ordered
206	while neighbours.size > 0:	×
207	match_next_vertex = np.any(np.isin(neighbours, next_neighbour), axis=1)	×
208
209	neighbours_order = np.vstack((neighbours_order, neighbours[match_next_vertex]))	×
210
211	next_neighbour = neighbours[match_next_vertex][0]	×
212	next_neighbour[next_neighbour == next_neighbour_prev] = 0	×
213	neighbours = np.delete(neighbours, match_next_vertex, axis=0)	×
214
215	next_neighbour_prev = next_neighbour	×
216
217	_, vert_order = ismember_rows(neighbours_order, np.array(initial_neighbours))	×
218
219	new_vert_order = np.vstack((vert_order, np.hstack((vert_order[1:], vert_order[0])))).T	×
220
221	return new_vert_order	×
222
223	# If ordering based on neighbours failed or no neighbours were provided,
224	# order the vertices based on their angular position relative to the centroid of the cell
225	imaginary_centroid_mean_vert = np.mean(vertices_of_cell, axis=0)	×
226	vector_for_ang_mean = vertices_of_cell - imaginary_centroid_mean_vert	×
227	th_mean = np.arctan2(vector_for_ang_mean[:, 1], vector_for_ang_mean[:, 0])	×
228	vert_order = np.argsort(th_mean)	×
229
230	new_vert_order = np.vstack((vert_order, np.hstack((vert_order[1:], vert_order[0])))).T	×
231
232	return new_vert_order	×
233
234
UNCOV 235	def build_2d_voronoi_from_image(labelled_img, watershed_img, total_cells):	×
236	"""
237	Build a 2D Voronoi diagram from an image
238	:param labelled_img:
239	:param watershed_img:
240	:param total_cells:
241	:return:
242	"""
243
244	ratio = 2	×
245
246	labelled_img[watershed_img == 0] = 0	×
247
248	# Create a mask for the edges with ID 0
249	edge_mask = labelled_img == 0	×
250
251	# Get the closest labeled polygon for each edge pixel
252	closest_id = dilation(labelled_img, square(5))	×
253
254	filled_image = closest_id	×
255	filled_image[~edge_mask] = labelled_img[~edge_mask]	×
256
257	labelled_img = copy.deepcopy(filled_image)	×
258
259	img_neighbours = calculate_neighbours(labelled_img, ratio)	×
260
261	# Calculate the network of neighbours from the cell with ID 1
262	if not isinstance(total_cells, np.ndarray):	×
263	main_cells = img_neighbours[1]	×
264	while len(main_cells) < total_cells:	×
265	for c_cell in main_cells:	×
266	for c_neighbour in img_neighbours[c_cell]:	×
267	if len(main_cells) >= total_cells:	×
268	break	×
269
270	if c_neighbour not in main_cells:	×
271	main_cells = np.append(main_cells, c_neighbour)	×
272
273	main_cells = np.sort(main_cells)	×
274	else:
275	main_cells = total_cells	×
276
277	# Calculate the border cells from main_cells
278	border_cells_and_main_cells = np.unique(np.block([img_neighbours[i] for i in main_cells]))	×
279	border_ghost_cells = np.setdiff1d(border_cells_and_main_cells, main_cells)	×
280	border_cells = np.intersect1d(main_cells, np.unique(np.block([img_neighbours[i] for i in border_ghost_cells])))	×
281
282	border_of_border_cells_and_main_cells = np.unique(	×
283	np.concatenate([img_neighbours[i] for i in border_cells_and_main_cells]))
284
285	quartets, _ = get_four_fold_vertices(img_neighbours)	×
286	if quartets is not None:	×
287	img_neighbours = divide_quartets_neighbours(img_neighbours, labelled_img, quartets)	×
288
289	# labelled_img[~np.isin(labelled_img, border_of_border_cells_and_main_cells)] = 0
290	vertices_info = populate_vertices_info(border_cells_and_main_cells, img_neighbours,	×
291	labelled_img, main_cells, ratio)
292
293	neighbours_network = generate_neighbours_network(img_neighbours, main_cells)	×
294
295	triangles_connectivity = np.array(vertices_info['connectedCells'])	×
296	cell_edges = vertices_info['edges']	×
297	vertices_location = vertices_info['location']	×
298
299	# Remove Nones from the vertices location
300	cell_edges = [cell_edges[i] for i in range(len(cell_edges)) if cell_edges[i] is not None]	×
301
302	return (triangles_connectivity, neighbours_network, cell_edges, vertices_location, border_cells,	×
303	border_of_border_cells_and_main_cells, main_cells)
304
305
UNCOV 306	def generate_neighbours_network(neighbours, main_cells):	×
307	neighbours_network = []	×
308	for num_cell in main_cells:	×
309	current_neighbours = np.array(neighbours[num_cell])	×
310	current_cell_neighbours = np.vstack(	×
311	[np.ones(len(current_neighbours), dtype=int) * num_cell, current_neighbours]).T
312
313	neighbours_network.extend(current_cell_neighbours)	×
314	return neighbours_network	×
315
316
UNCOV 317	def divide_quartets_neighbours(img_neighbours_all, labelled_img, quartets):	×
318	"""
319	Divide the quartets of neighboring cells into two pairs of cells.
320	:param img_neighbours_all:
321	:param labelled_img:
322	:param quartets:
323	:return:
324	"""
325	props = regionprops_table(labelled_img, properties=('centroid', 'label',))	×
326	# The centroids are now stored in 'props' as separate arrays 'centroid-0', 'centroid-1', etc.
327	# You can combine them into a single array like this:
328	face_centres_vertices = np.column_stack([props['centroid-0'], props['centroid-1']])	×
329	# Loop through the quartets and split them into two pairs of cells
330	for num_quartets in range(quartets.shape[0]):	×
331	# Get the face centres of the current quartet whose ids correspond to props['label']
332	quartets_ids = [np.where(props['label'] == i)[0][0] for i in quartets[num_quartets]]	×
333	current_centroids = face_centres_vertices[quartets_ids, :]	×
334
335	# Get the distance between the centroids of the current quartet
336	distance_between_centroids = squareform(pdist(current_centroids))	×
337
338	# Get the maximum distance between the centroids
339	max_distance = np.max(distance_between_centroids)	×
340	row, col = np.where(distance_between_centroids == max_distance)	×
341
342	# Split the quartet into two pairs of cells
343	current_neighs = img_neighbours_all[quartets[num_quartets][col[0]]]	×
344	current_neighs = current_neighs[current_neighs != quartets[num_quartets][row[0]]]	×
345	img_neighbours_all[quartets[num_quartets][col[0]]] = current_neighs	×
346
347	current_neighs = img_neighbours_all[quartets[num_quartets][row[0]]]	×
348	current_neighs = current_neighs[current_neighs != quartets[num_quartets][col[0]]]	×
349	img_neighbours_all[quartets[num_quartets][row[0]]] = current_neighs	×
350
351	return img_neighbours_all	×
352
353
UNCOV 354	def populate_vertices_info(border_cells_and_main_cells, img_neighbours_all, labelled_img, main_cells,	×
355	ratio):
356	"""
357	Populate the vertices information.
358	:param border_cells_and_main_cells:
359	:param img_neighbours_all:
360	:param labelled_img:
361	:param main_cells:
362	:param ratio:
363	:return:
364	"""
365	vertices_info = calculate_vertices(labelled_img, img_neighbours_all, ratio)	×
366	total_cells = np.max(border_cells_and_main_cells) + 1	×
367	vertices_info['PerCell'] = [None] * total_cells	×
368	vertices_info['edges'] = [None] * total_cells	×
369	for idx, num_cell in enumerate(main_cells):	×
370	vertices_of_cell = np.where(np.any(np.isin(vertices_info['connectedCells'], num_cell), axis=1))[0]	×
371	vertices_info['PerCell'][idx] = vertices_of_cell	×
372	current_vertices = [vertices_info['location'][i] for i in vertices_of_cell]	×
373	current_connected_cells = [vertices_info['connectedCells'][i] for i in vertices_of_cell]	×
374
375	# Remove the current cell 'num_cell' from the connected cells
376	current_connected_cells = [np.delete(cell, np.where(cell == num_cell)) for cell in current_connected_cells]	×
377
378	vertices_info['edges'][idx] = vertices_of_cell[	×
379	boundary_of_cell(current_vertices, current_connected_cells)]
380	assert (len(vertices_info['edges'][idx]) ==	×
381	len(img_neighbours_all[num_cell])), 'Error missing vertices of neighbours'
382	return vertices_info	×
383
384
UNCOV 385	def process_image(img_filename="src/pyVertexModel/resources/LblImg_imageSequence.tif", redo=False):	×
386	"""
387	Process the image and return the 2D labeled image and the 3D labeled image.
388	:param redo:
389	:param img_filename:
390	:return:
391	"""
392	# Load the tif file from resources if exists
393	if os.path.exists(img_filename):	×
394	if img_filename.endswith('.tif'):	×
395	if os.path.exists(img_filename.replace('.tif', '.xz')) and not redo:	×
396	imgStackLabelled = pickle.load(lzma.open(img_filename.replace('.tif', '.xz'), "rb"))	×
397
398	imgStackLabelled = imgStackLabelled['imgStackLabelled']	×
399	img2DLabelled = imgStackLabelled[0, :, :]	×
400	else:
401	imgStackLabelled = io.imread(img_filename)	×
402
403	# Reordering cells based on the centre of the image
404	img2DLabelled = imgStackLabelled[0, :, :]	×
405	props = regionprops_table(img2DLabelled, properties=('centroid', 'label',), )	×
406
407	# The centroids are now stored in 'props' as separate arrays 'centroid-0', 'centroid-1', etc.
408	centroids = np.column_stack([props['centroid-0'], props['centroid-1']])	×
409	centre_of_image = np.array([img2DLabelled.shape[0] / 2, img2DLabelled.shape[1] / 2])	×
410
411	# Sorting cells based on distance to the middle of the image
412	distanceToMiddle = cdist([centre_of_image], centroids)	×
413	distanceToMiddle = distanceToMiddle[0]	×
414	sortedId = np.argsort(distanceToMiddle)	×
415	sorted_ids = np.array(props['label'])[sortedId]	×
416
417	oldImgStackLabelled = copy.deepcopy(imgStackLabelled)	×
418	# imgStackLabelled = np.zeros_like(imgStackLabelled)
419	newCont = 1	×
420	for numCell in sorted_ids:	×
421	if numCell != 0:	×
422	imgStackLabelled[oldImgStackLabelled == numCell] = newCont	×
423	newCont += 1	×
424
425	# Remaining cells that are not in the image
426	for numCell in np.arange(newCont, np.max(img2DLabelled) + 1):	×
427	imgStackLabelled[oldImgStackLabelled == numCell] = newCont	×
428	newCont += 1	×
429
430	img2DLabelled = imgStackLabelled[0, :, :]	×
431
432	save_variables({'imgStackLabelled': imgStackLabelled},	×
433	img_filename.replace('.tif', '.xz'))
434
435	imgStackLabelled = np.transpose(imgStackLabelled, (2, 0, 1))	×
436	elif img_filename.endswith('.mat'):	×
437	imgStackLabelled = scipy.io.loadmat(img_filename)['imgStackLabelled']	×
438	img2DLabelled = imgStackLabelled[:, :, 0]	×
439	else:
440	raise ValueError('Image file not found %s' % img_filename)	×
441
442	return img2DLabelled, imgStackLabelled	×
443
444
UNCOV 445	def add_tetrahedral_intercalations(Twg, xInternal, XgBottom, XgTop, XgLateral):	×
446	allCellIds = np.concatenate([xInternal, XgLateral])	×
447	neighboursMissing = {}	×
448
449	for numCell in xInternal:	×
450	Twg_cCell = Twg[np.any(np.isin(Twg, numCell), axis=1)]	×
451
452	Twg_cCell_bottom = Twg_cCell[np.any(np.isin(Twg_cCell, XgBottom), axis=1), :]	×
453	neighbours_bottom = allCellIds[np.isin(allCellIds, Twg_cCell_bottom)]	×
454
455	Twg_cCell_top = Twg_cCell[np.any(np.isin(Twg_cCell, XgTop), axis=1), :]	×
456	neighbours_top = allCellIds[np.isin(allCellIds, Twg_cCell_top)]	×
457
458	neighboursMissing[numCell] = np.setxor1d(neighbours_bottom, neighbours_top)	×
459	for missingCell in neighboursMissing[numCell]:	×
460	tetsToAdd = np.sort(allCellIds[	×
461	np.isin(allCellIds, Twg_cCell[np.any(np.isin(Twg_cCell, missingCell), axis=1), :])])
462	assert len(tetsToAdd) == 4, f'Missing 4-fold at Cell {numCell}'	×
463	if not np.any(np.all(tetsToAdd == Twg, axis=1)):	×
464	Twg = np.vstack((Twg, tetsToAdd))	×
465	return Twg	×
466
467
UNCOV 468	class VertexModelVoronoiFromTimeImage(VertexModel):	×
UNCOV 469	def __init__(self, set_test=None, update_derived_parameters=True, create_output_folder=True):	×
470	super().__init__(set_test, update_derived_parameters, create_output_folder)	×
471
UNCOV 472	def initialize(self):	×
473	"""
474	Initialize the geometry and the topology of the model.
475	"""
476	filename = os.path.join(PROJECT_DIRECTORY, self.set.initial_filename_state)	×
477
478	if not os.path.exists(filename):	×
479	logging.error(f'File {filename} not found')	×
480
481	if filename.endswith('.pkl'):	×
482	output_folder = self.set.OutputFolder	×
483	load_state(self, filename, ['geo', 'geo_0', 'geo_n'])	×
484	self.set.OutputFolder = output_folder	×
485	elif filename.endswith('.mat'):	×
486	mat_info = scipy.io.loadmat(filename)	×
487	self.geo = Geo(mat_info['Geo'])	×
488	self.geo.update_measures()	×
489	else:
490	# Load the image and obtain the initial X and tetrahedra
491	Twg, X = self.obtain_initial_x_and_tetrahedra()	×
492	# Build cells
493	self.geo.build_cells(self.set, X, Twg)	×
494
495	# Save state with filename using the number of cells
496	filename = filename.replace('.tif', f'_{self.set.TotalCells}cells.pkl')	×
497	# save_state(self.geo, 'voronoi_40cells.pkl')
498
499	# Add border cells to the shared cells
500	for cell in self.geo.Cells:	×
501	if cell.ID in self.geo.BorderCells:	×
502	for face in cell.Faces:	×
503	for tris in face.Tris:	×
504	tets_1 = cell.T[tris.Edge[0]]	×
505	tets_2 = cell.T[tris.Edge[1]]	×
506	shared_cells = np.intersect1d(tets_1, tets_2)	×
507	if np.any(np.isin(self.geo.BorderGhostNodes, shared_cells)):	×
508	shared_cells_list = list(tris.SharedByCells)	×
509	shared_cells_list.append(shared_cells[np.isin(shared_cells, self.geo.BorderGhostNodes)][0])	×
510	tris.SharedByCells = np.array(shared_cells_list)	×
511
512	# Create periodic boundary conditions
513	self.geo.apply_periodic_boundary_conditions(self.set)	×
514
515	if self.set.ablation:	×
516	self.geo.cellsToAblate = self.set.cellsToAblate	×
517
518	# Define upper and lower area threshold for remodelling
519	if self.geo.lmin0 is None:	×
520	self.geo.init_reference_cell_values(self.set)	×
521
UNCOV 522	def obtain_initial_x_and_tetrahedra(self, img_filename=None):	×
523	"""
524	Obtain the initial X and tetrahedra for the model.
525	:return:
526	"""
527	if img_filename is None:	×
528	img_filename = PROJECT_DIRECTORY + '/src/pyVertexModel/resources/LblImg_imageSequence.tif'	×
529
530	selectedPlanes = [0, 99]	×
531	img2DLabelled, imgStackLabelled = process_image(img_filename)	×
532
533	# Building the topology of each plane
534	trianglesConnectivity = {}	×
535	neighboursNetwork = {}	×
536	cellEdges = {}	×
537	verticesOfCell_pos = {}	×
538	borderCells = {}	×
539	borderOfborderCellsAndMainCells = {}	×
540	cell_centroids = {}	×
541	main_cells = self.set.TotalCells	×
542	for numPlane in selectedPlanes:	×
543	current_img = imgStackLabelled[:, numPlane, :]	×
544	(triangles_connectivity, neighbours_network,	×
545	cell_edges, vertices_location, border_cells,
546	border_of_border_cells_and_main_cells,
547	main_cells) = build_2d_voronoi_from_image(current_img, current_img, main_cells)
548
549	trianglesConnectivity[numPlane] = triangles_connectivity	×
550	neighboursNetwork[numPlane] = neighbours_network	×
551	cellEdges[numPlane] = cell_edges	×
552	verticesOfCell_pos[numPlane] = vertices_location	×
553	borderCells[numPlane] = border_cells	×
554	borderOfborderCellsAndMainCells[numPlane] = border_of_border_cells_and_main_cells	×
555
556	props = regionprops_table(current_img, properties=('centroid', 'label',))	×
557	zeros_column = np.zeros((props['centroid-1'].shape[0], 1))	×
558	cell_centroids[numPlane] = np.column_stack([props['label'], props['centroid-0'], props['centroid-1'],	×
559	zeros_column])
560
561	# Select nodes from images
562	all_main_cells = np.unique(	×
563	np.concatenate([borderOfborderCellsAndMainCells[numPlane] for numPlane in selectedPlanes]))
564
565	# Obtain the average centroid between the selected planes
566	max_label = int(np.max([np.max(cell_centroids[numPlane][:, 0]) for numPlane in selectedPlanes]))	×
567	X = np.zeros((max_label + 1, 3))	×
568	for label in arange(1, max_label + 1):	×
569	list_centroids = []	×
570	for numPlane in selectedPlanes:	×
571	found_centroid = np.where(cell_centroids[numPlane][:, 0] == label)[0]	×
572	if len(found_centroid) > 0:	×
573	list_centroids.append(cell_centroids[numPlane][found_centroid, 1:3])	×
574
575	if len(list_centroids) > 0:	×
576	X[label, 0:2] = np.mean(list_centroids, axis=0)	×
577	else:
578	print(f'Cell {label} not found in all planes')	×
579
580	# Basic features
581	cell_height = calculate_cell_height_on_model(img2DLabelled, main_cells, self.set)	×
582
583	# Generate tetrahedra from information of images
584	Twg, X = generate_tetrahedra_from_information(X, cellEdges, cell_height, cell_centroids, main_cells,	×
585	neighboursNetwork, selectedPlanes, trianglesConnectivity,
586	verticesOfCell_pos, self.geo)
587
588	# Fill Geo info
589	self.geo.nCells = len(main_cells)	×
590	self.geo.Main_cells = main_cells	×
591	self.geo.XgLateral = np.setdiff1d(all_main_cells, main_cells)	×
592	self.geo.XgID = np.setdiff1d(np.arange(1, X.shape[0] + 1), main_cells)	×
593	# Define border cells
594	self.geo.BorderCells = np.unique(np.concatenate([borderCells[numPlane] for numPlane in selectedPlanes]))	×
595	self.geo.BorderGhostNodes = self.geo.XgLateral	×
596
597	# Create new tetrahedra based on intercalations
598	Twg = add_tetrahedral_intercalations(Twg, main_cells, self.geo.XgBottom, self.geo.XgTop, self.geo.XgLateral)	×
599
600	# After removing ghost tetrahedras, some nodes become disconnected,
601	# that is, not a part of any tetrahedra. Therefore, they should be
602	# removed from X
603	Twg = Twg[~np.all(np.isin(Twg, self.geo.XgID), axis=1)]	×
604
605	# Remove tetrahedra with cells that are not in all_main_cells
606	# cells_to_remove = np.setdiff1d(range(1, np.max(all_main_cells) + 1), all_main_cells)
607	# Twg = Twg[~np.any(np.isin(Twg, cells_to_remove), axis=1)]
608
609	# Re-number the surviving tets
610	Twg, X = self.renumber_tets_xs(Twg, X)	×
611	# Normalise Xs
612	X = X / img2DLabelled.shape[0]	×
613
614	return Twg, X	×
615
UNCOV 616	def renumber_tets_xs(self, Twg, X):	×
617	"""
618	Renumber the tetrahedra and the coordinates.
619
620	This function renumbers the tetrahedra and the coordinates based on the unique values in the tetrahedra array.
621	It also updates the ghost node indices in the geometry object.
622
623	Parameters:
624	Twg (numpy.ndarray): A 2D array where each row represents a tetrahedron and each column represents a node of the tetrahedron.
625	X (numpy.ndarray): A 2D array where each row represents a node and each column represents a coordinate of the node.
626
627	Returns:
628	Twg (numpy.ndarray): The renumbered tetrahedra array.
629	X (numpy.ndarray): The renumbered coordinates array.
630	"""
631	# Get the unique values in the tetrahedra array and their inverse mapping
632	oldIds, oldTwgNewIds = np.unique(Twg, return_inverse=True)	×
633	# Create a new array of indices
634	newIds = np.arange(len(oldIds))	×
635	# Update the coordinates array based on the old indices
636	X = X[oldIds - 1, :]	×
637	# Reshape the inverse mapping to match the shape of the tetrahedra array
638	Twg = oldTwgNewIds.reshape(Twg.shape)	×
639	# Update the ghost node indices in the geometry object
640	self.geo.XgBottom = newIds[np.isin(oldIds, self.geo.XgBottom)]	×
641	self.geo.XgTop = newIds[np.isin(oldIds, self.geo.XgTop)]	×
642	self.geo.XgLateral = newIds[np.isin(oldIds, self.geo.XgLateral)]	×
643	self.geo.XgID = newIds[np.isin(oldIds, self.geo.XgID)]	×
644	self.geo.BorderGhostNodes = self.geo.XgLateral	×
645	self.geo.Main_cells = newIds[np.isin(oldIds, self.geo.Main_cells)]	×
646	# Return the renumbered tetrahedra and coordinates arrays
647	return Twg, X	×
648
UNCOV 649	def copy(self):	×
650	"""
651	Copy the object.
652	"""
653	return super().copy()	×
654
UNCOV 655	def calculate_error(self, K, initial_recoil, error_type=None):	×
656	"""
657	Calculate the error.
658	"""
659	return super().calculate_error(K, initial_recoil, error_type)	×

Pablo1990 / pyVertexModel / 11363842044

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous