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/vertexModel.py

UNCOV 1	import logging	×
UNCOV 2	import os	×
UNCOV 3	from abc import abstractmethod	×
UNCOV 4	from itertools import combinations	×
5
UNCOV 6	import imageio	×
UNCOV 7	import numpy as np	×
UNCOV 8	import pandas as pd	×
UNCOV 9	import pyvista as pv	×
UNCOV 10	from scipy.stats import zscore	×
UNCOV 11	from skimage.measure import regionprops	×
12
UNCOV 13	from src.pyVertexModel.algorithm import newtonRaphson	×
UNCOV 14	from src.pyVertexModel.geometry import degreesOfFreedom	×
UNCOV 15	from src.pyVertexModel.geometry.geo import Geo	×
UNCOV 16	from src.pyVertexModel.mesh_remodelling.remodelling import Remodelling	×
UNCOV 17	from src.pyVertexModel.parameters.set import Set	×
UNCOV 18	from src.pyVertexModel.util.utils import save_state, save_backup_vars, load_backup_vars, copy_non_mutable_attributes	×
19
UNCOV 20	logger = logging.getLogger("pyVertexModel")	×
21
22
UNCOV 23	def generate_tetrahedra_from_information(X, cell_edges, cell_height, cell_centroids, main_cells,	×
24	neighbours_network, selected_planes, triangles_connectivity,
25	vertices_of_cell_pos, geo):
26	"""
27	Generate tetrahedra from the information of the cells.
28	:param X:
29	:param cell_edges:
30	:param cell_height:
31	:param cell_centroids:
32	:param main_cells:
33	:param neighbours_network:
34	:param selected_planes:
35	:param triangles_connectivity:
36	:param vertices_of_cell_pos:
37	:param geo:
38	:return:
39	"""
40	bottom_plane = 0	×
41	top_plane = 1	×
42	if bottom_plane == 0:	×
43	z_coordinate = [-cell_height, cell_height]	×
44	else:
45	z_coordinate = [cell_height, -cell_height]	×
46	Twg = []	×
47	for idPlane, numPlane in enumerate(selected_planes):	×
48	# Using the centroids and vertices of the cells of each 2D image as ghost nodes
49	X, Xg_faceIds, Xg_ids, Xg_verticesIds = (	×
50	add_faces_and_vertices_to_x(X,
51	np.hstack((cell_centroids[numPlane][:, 1:3], np.tile(z_coordinate[idPlane],
52	(len(cell_centroids[
53	numPlane][:,
54	1:3]), 1)))),
55	np.hstack((np.fliplr(vertices_of_cell_pos[numPlane]),
56	np.tile(z_coordinate[idPlane],
57	(len(vertices_of_cell_pos[numPlane]), 1))))))
58
59	# Fill Geo info
60	if idPlane == bottom_plane:	×
61	geo.XgBottom = Xg_ids	×
62	elif idPlane == top_plane:	×
63	geo.XgTop = Xg_ids	×
64
65	# Create tetrahedra
66	Twg_numPlane = create_tetrahedra(triangles_connectivity[numPlane], neighbours_network[numPlane],	×
67	cell_edges[numPlane], main_cells, Xg_faceIds, Xg_verticesIds)
68
69	Twg.append(Twg_numPlane)	×
70	Twg = np.vstack(Twg)	×
71	return Twg, X	×
72
73
UNCOV 74	def add_faces_and_vertices_to_x(X, Xg_faceCentres2D, Xg_vertices2D):	×
75	"""
76	Add faces and vertices to the X matrix.
77	:param X:
78	:param Xg_faceCentres2D:
79	:param Xg_vertices2D:
80	:return:
81	"""
82	Xg_nodes = np.vstack((Xg_faceCentres2D, Xg_vertices2D))	×
83	Xg_ids = np.arange(X.shape[0] + 1, X.shape[0] + Xg_nodes.shape[0] + 1)	×
84	Xg_faceIds = Xg_ids[0:Xg_faceCentres2D.shape[0]]	×
85	Xg_verticesIds = Xg_ids[Xg_faceCentres2D.shape[0]:]	×
86	X = np.vstack((X, Xg_nodes))	×
87	return X, Xg_faceIds, Xg_ids, Xg_verticesIds	×
88
89
UNCOV 90	def create_tetrahedra(triangles_connectivity, neighbours_network, edges_of_vertices, x_internal, x_face_ids,	×
91	x_vertices_ids):
92	"""
93	Add connections between real nodes and ghost cells to create tetrahedra.
94
95	:param triangles_connectivity: A 2D array where each row represents a triangle connectivity.
96	:param neighbours_network: A 2D array where each row represents a pair of neighboring nodes.
97	:param edges_of_vertices: A list of lists where each sublist represents the edges of a vertex.
98	:param x_internal: A 1D array representing the internal nodes.
99	:param x_face_ids: A 1D array representing the face ids.
100	:param x_vertices_ids: A 1D array representing the vertices ids.
101	:return: A 2D array representing the tetrahedra.
102	"""
103	x_ids = np.concatenate([x_face_ids, x_vertices_ids])	×
104
105	# Relationships: 1 ghost node, three cell nodes
106	twg = np.hstack([triangles_connectivity, x_vertices_ids[:, None]])	×
107
108	# Relationships: 1 cell node and 3 ghost nodes
109	new_additions = []	×
110	for id_cell, num_cell in enumerate(x_internal):	×
111	face_id = x_face_ids[id_cell]	×
112	vertices_to_connect = edges_of_vertices[id_cell]	×
113	new_additions.extend(np.hstack([np.repeat(np.array([[num_cell, face_id]]), len(vertices_to_connect), axis=0),	×
114	x_vertices_ids[vertices_to_connect]]))
115	twg = np.vstack([twg, new_additions])	×
116
117	# Relationships: 2 ghost nodes, two cell nodes
118	twg_sorted = np.sort(twg[np.any(np.isin(twg, x_ids), axis=1)], axis=1)	×
119	internal_neighbour_network = [neighbour for neighbour in neighbours_network if	×
120	np.any(np.isin(neighbour, x_internal))]
121	internal_neighbour_network = np.unique(np.sort(internal_neighbour_network, axis=1), axis=0)	×
122
123	new_additions = []	×
124	for num_pair in range(internal_neighbour_network.shape[0]):	×
125	found = np.isin(twg_sorted, internal_neighbour_network[num_pair])	×
126	new_connections = np.unique(twg_sorted[np.sum(found, axis=1) == 2, 3])	×
127	if len(new_connections) > 1:	×
128	new_connections_pairs = np.array(list(combinations(new_connections, 2)))	×
129	new_additions.extend([np.hstack([internal_neighbour_network[num_pair], new_connections_pair])	×
130	for new_connections_pair in new_connections_pairs])
131	else:
132	raise ValueError('Error while creating the connections and initial topology')	×
133	twg = np.vstack([twg, new_additions])	×
134
135	return twg	×
136
137
UNCOV 138	def calculate_cell_height_on_model(img2DLabelled, main_cells, c_set):	×
139	"""
140	Calculate the cell height on the model regarding the diameter of the cells.
141	:param img2DLabelled:
142	:param main_cells:
143	:return:
144	"""
145	properties = regionprops(img2DLabelled)	×
146	# Extract major axis lengths
147	avg_diameter = np.mean([prop.major_axis_length for prop in properties if prop.label in main_cells])	×
148	cell_height = avg_diameter * c_set.CellHeight	×
149	return cell_height	×
150
151
UNCOV 152	class VertexModel:	×
153	"""
154	The main class for the vertex model simulation. It contains the methods for initializing the model,
155	iterating over time, applying Brownian motion, and checking the integrity of the model.
156	"""
157
UNCOV 158	def __init__(self, c_set=None, create_output_folder=True, update_derived_parameters=True):	×
159	"""
160	Vertex Model class.
161	:param c_set:
162	"""
163	self.colormap_lim = None	×
164	self.OutputFolder = None	×
165	self.numStep = None	×
166	self.backupVars = None	×
167	self.geo_n = None	×
168	self.geo_0 = None	×
169	self.tr = None	×
170	self.t = None	×
171	self.X = None	×
172	self.didNotConverge = False	×
173	self.geo = Geo()	×
174
175	# Set definition
176	if c_set is not None:	×
177	self.set = c_set	×
178	else:
179	# TODO Create a menu to select the set
180	self.set = Set()	×
181	# self.set.cyst()
182	self.set.wing_disc()	×
183	if self.set.ablation:	×
184	self.set.wound_default()	×
185
186	if update_derived_parameters:	×
187	self.set.update_derived_parameters()	×
188
189	# Redirect output
190	if self.set.OutputFolder is not None and create_output_folder:	×
191	self.set.redirect_output()	×
192
193	# Degrees of freedom definition
194	self.Dofs = degreesOfFreedom.DegreesOfFreedom()	×
195
196	self.relaxingNu = False	×
197	self.EnergiesPerTimeStep = []	×
198	self.t = 0	×
199	self.tr = 0	×
200	self.numStep = 1	×
201
UNCOV 202	@abstractmethod	×
UNCOV 203	def initialize(self):	×
204	pass	×
205
UNCOV 206	def brownian_motion(self, scale):	×
207	"""
208	Applies Brownian motion to the vertices of cells in the Geo structure.
209	Displacements are generated with a normal distribution in each dimension.
210	:param scale:
211	:return:
212	"""
213
214	# Concatenate and sort all tetrahedron vertices
215	all_tets = np.sort(np.vstack([cell.T for cell in self.geo.Cells]), axis=1)	×
216	all_tets_unique = np.unique(all_tets, axis=0)	×
217
218	# Generate random displacements with a normal distribution for each dimension
219	displacements = scale * (self.geo.Cells[0].X - self.geo.Cells[1].X) * np.random.randn(all_tets_unique.shape[0],	×
220	3)
221
222	# Update vertex positions based on 3D Brownian motion displacements
223	for cell in [c for c in self.geo.Cells if c.AliveStatus is not None and c.ID not in self.geo.BorderCells]:	×
224	_, corresponding_ids = np.where(np.all(np.sort(cell.T, axis=1)[:, None] == all_tets_unique, axis=2))	×
225	cell.Y += displacements[corresponding_ids, :]	×
226
UNCOV 227	def iterate_over_time(self):	×
228	"""
229	Iterate the model over time. This includes updating the degrees of freedom, applying boundary conditions,
230	updating measures, and checking for convergence.
231	:return:
232	"""
233	temp_dir = os.path.join(self.set.OutputFolder, 'images')	×
234	if not os.path.exists(temp_dir):	×
235	os.makedirs(temp_dir)	×
236
237	if self.set.Substrate == 1:	×
238	self.Dofs.GetDOFsSubstrate(self.geo, self.set)	×
239	else:
240	self.Dofs.get_dofs(self.geo, self.set)	×
241
242	self.geo.remodelling = False	×
243	if self.geo_0 is None:	×
244	self.geo_0 = self.geo.copy(update_measurements=False)	×
245
246	if self.geo_n is None:	×
247	self.geo_n = self.geo.copy(update_measurements=False)	×
248
249	# Count the number of faces in average has a cell per domain
250	self.geo.update_barrier_tri0_based_on_number_of_faces()	×
251	self.backupVars = save_backup_vars(self.geo, self.geo_n, self.geo_0, self.tr, self.Dofs)	×
252
253	print("File: ", self.set.OutputFolder)	×
254	self.save_v_model_state()	×
255
256	while self.t <= self.set.tend and not self.didNotConverge:	×
257	gr = self.single_iteration()	×
258
259	if np.isnan(gr):	×
260	break	×
261
262	return self.didNotConverge	×
263
UNCOV 264	def single_iteration(self, post_operations=True):	×
265	"""
266	Perform a single iteration of the model.
267	:return:
268	"""
269	self.set.currentT = self.t	×
270	logger.info("Time: " + str(self.t))	×
271	if not self.relaxingNu:	×
272	self.set.i_incr = self.numStep	×
273
274	# Ablate cells if needed
275	if self.set.ablation:	×
276	if self.set.ablation and self.set.TInitAblation <= self.t and self.geo.cellsToAblate is not None:	×
277	self.save_v_model_state(file_name='before_ablation')	×
278	self.geo.ablate_cells(self.set, self.t)	×
279	self.geo_n = self.geo.copy()	×
280	# Update the degrees of freedom
281	self.Dofs.get_dofs(self.geo, self.set)	×
282
283	self.Dofs.ApplyBoundaryCondition(self.t, self.geo, self.set)	×
284	# IMPORTANT: Here it updates: Areas, Volumes, etc... Should be
285	# up-to-date
286	self.geo.update_measures()	×
287	if self.set.implicit_method is True:	×
288	g, K, _, energies = newtonRaphson.KgGlobal(self.geo_0, self.geo_n, self.geo, self.set,	×
289	self.set.implicit_method)
290	else:
291	K = 0	×
292	g, energies = newtonRaphson.gGlobal(self.geo_0, self.geo_n, self.geo, self.set,	×
293	self.set.implicit_method)
294	for key, energy in energies.items():	×
295	logger.info(f"{key}: {energy}")	×
296	self.geo, g, __, __, self.set, gr, dyr, dy = newtonRaphson.newton_raphson(self.geo_0, self.geo_n, self.geo,	×
297	self.Dofs, self.set, K, g,
298	self.numStep, self.t,
299	self.set.implicit_method)
300	if not np.isnan(gr) and post_operations:	×
301	self.post_newton_raphson(dy, dyr, g, gr)	×
302	return gr	×
303
UNCOV 304	def post_newton_raphson(self, dy, dyr, g, gr):	×
305	"""
306	Post Newton Raphson operations.
307	:param dy:
308	:param dyr:
309	:param g:
310	:param gr:
311	:return:
312	"""
313	if (gr < self.set.tol and dyr < self.set.tol and np.all(~np.isnan(g[self.Dofs.Free])) and	×
314	np.all(~np.isnan(dy[self.Dofs.Free]))):
315	self.iteration_converged()	×
316	# if self.set.implicit_method is False:
317	# self.set.tol = gr
318	# if self.set.tol < self.set.tol0:
319	# self.set.tol = self.set.tol0
320	else:
321	self.iteration_did_not_converged()	×
322
323	self.Dofs.get_dofs(self.geo, self.set)	×
324
UNCOV 325	def iteration_did_not_converged(self):	×
326	"""
327	If the iteration did not converge, the algorithm will try to relax the value of nu and dt.
328	:return:
329	"""
330	self.geo, self.geo_n, self.geo_0, self.tr, self.Dofs = load_backup_vars(self.backupVars)	×
331	self.relaxingNu = False	×
332	if self.set.iter == self.set.MaxIter0 and self.set.implicit_method:	×
333	self.set.MaxIter = self.set.MaxIter0 * 3	×
334	self.set.nu = 10 * self.set.nu0	×
335	else:
336	if (self.set.iter >= self.set.MaxIter and	×
337	(self.set.dt / self.set.dt0) > 1e-6):
338	self.set.MaxIter = self.set.MaxIter0	×
339	self.set.nu = self.set.nu0	×
340	self.set.dt = self.set.dt / 2	×
341	self.t = self.set.last_t_converged + self.set.dt	×
342	else:
343	self.didNotConverge = True	×
344
UNCOV 345	def iteration_converged(self):	×
346	"""
347	If the iteration converged, the algorithm will update the values of the variables and proceed to the next step.
348	:return:
349	"""
350	if self.set.nu / self.set.nu0 == 1:	×
351	# STEP has converged
352	logger.info(f"STEP {str(self.set.i_incr)} has converged ...")	×
353
354	# Remodelling
355	if abs(self.t - self.tr) >= self.set.RemodelingFrequency:	×
356	if self.set.Remodelling:	×
357	save_state(self,	×
358	os.path.join(self.set.OutputFolder,
359	'data_step_before_remodelling_' + str(self.numStep) + '.pkl'))
360	remodel_obj = Remodelling(self.geo, self.geo_n, self.geo_0, self.set, self.Dofs)	×
361	self.geo, self.geo_n = remodel_obj.remodel_mesh(self.numStep)	×
362	# Update tolerance if remodelling was performed to the current one
363	if self.set.implicit_method is False:	×
364	g, energies = newtonRaphson.gGlobal(self.geo_0, self.geo_n, self.geo, self.set,	×
365	self.set.implicit_method)
366	self.Dofs.get_dofs(self.geo, self.set)	×
367	gr = np.linalg.norm(g[self.Dofs.Free])	×
368
369	# Append Energies
370	# energies_per_time_step.append(energies)
371
372	# Build X From Y
373	self.geo.build_x_from_y(self.geo_n)	×
374
375	# Update last time converged
376	self.set.last_t_converged = self.t	×
377
378	# Test Geo
379	# TODO: CHECK
380	# self.check_integrity()
381
382	if abs(self.t - self.tr) >= self.set.RemodelingFrequency:	×
383	self.save_v_model_state()	×
384
385	# Reset noise to be comparable between simulations
386	self.reset_noisy_parameters()	×
387	self.tr = self.t	×
388
389	# Brownian Motion
390	if self.set.brownian_motion is True:	×
391	self.brownian_motion(self.set.brownian_motion_scale)	×
392
393	self.t = self.t + self.set.dt	×
394	self.set.dt = np.min([self.set.dt + self.set.dt * 0.5, self.set.dt0])	×
395	self.set.MaxIter = self.set.MaxIter0	×
396	self.numStep = self.numStep + 1	×
397	self.backupVars = {	×
398	'Geo_b': self.geo.copy(),
399	'Geo_n_b': self.geo_n.copy(),
400	'Geo_0_b': self.geo_0.copy(),
401	'tr_b': self.tr,
402	'Dofs': self.Dofs.copy()
403	}
404	self.geo_n = self.geo.copy()	×
405	self.relaxingNu = False	×
406	else:
407	self.set.nu = np.max([self.set.nu / 2, self.set.nu0])	×
408	self.relaxingNu = True	×
409
UNCOV 410	def save_v_model_state(self, file_name=None):	×
411	"""
412	Save the state of the vertex model.
413	:param file_name:
414	:return:
415	"""
416	# Create VTK files for the current state
417	self.geo.create_vtk_cell(self.set, self.numStep, 'Edges')	×
418	self.geo.create_vtk_cell(self.set, self.numStep, 'Cells')	×
419	temp_dir = os.path.join(self.set.OutputFolder, 'images')	×
420	self.screenshot(temp_dir)	×
421	# Save Data of the current step
422	if file_name is None:	×
423	save_state(self, os.path.join(self.set.OutputFolder, 'data_step_' + str(self.numStep) + '.pkl'))	×
424	else:
425	save_state(self, os.path.join(self.set.OutputFolder, file_name + '.pkl'))	×
426
UNCOV 427	def reset_noisy_parameters(self):	×
428	"""
429	Reset noisy parameters.
430	:return:
431	"""
432	for num_cell in range(len(self.geo.Cells)):	×
433	c_cell = self.geo.Cells[num_cell]	×
434	self.geo.Cells[num_cell].contractlity_noise = None	×
435	self.geo.Cells[num_cell].lambda_s1_noise = None	×
436	self.geo.Cells[num_cell].lambda_s2_noise = None	×
437	self.geo.Cells[num_cell].lambda_s3_noise = None	×
438	self.geo.Cells[num_cell].lambda_v_noise = None	×
439	for n_face in range(len(c_cell.Faces)):	×
440	face = c_cell.Faces[n_face]	×
441	for n_tri in range(len(face.Tris)):	×
442	tri = face.Tris[n_tri]	×
443	tri.ContractilityValue = None	×
444	tri.lambda_r_noise = None	×
445	tri.lambda_b_noise = None	×
446	tri.k_substrate_noise = None	×
447	# tri.edge_length_time.append([self.t, tri.edge_length])
448	self.geo.Cells[num_cell].Faces[n_face].Tris[n_tri] = tri	×
449
UNCOV 450	def check_integrity(self):	×
451	"""
452	Performs tests on the properties of cells, faces, and triangles (tris) within the Geo structure.
453	Ensures that certain geometrical properties are above minimal threshold values.
454	"""
455
456	# Define minimum error thresholds for edge length, area, and volume
457	min_error_edge = 1e-5	×
458	min_error_area = min_error_edge ** 2	×
459	min_error_volume = min_error_edge ** 3	×
460
461	# Test Cells properties:
462	# Conditions checked:
463	# - Volume > minimum error volume
464	# - Initial Volume > minimum error volume
465	# - Area > minimum error area
466	# - Initial Area > minimum error area
467	for c_cell in self.geo.Cells:	×
468	if c_cell.AliveStatus:	×
469	assert c_cell.Vol > min_error_volume, "Cell volume is too low"	×
470	assert c_cell.Vol0 > min_error_volume, "Cell initial volume is too low"	×
471	assert c_cell.Area > min_error_area, "Cell area is too low"	×
472	assert c_cell.Area0 > min_error_area, "Cell initial area is too low"	×
473
474	# Test Faces properties:
475	# Conditions checked:
476	# - Area > minimum error area
477	# - Initial Area > minimum error area
478	for c_cell in self.geo.Cells:	×
479	if c_cell.AliveStatus:	×
480	for face in c_cell.Faces:	×
481	assert face.Area > min_error_area, "Face area is too low"	×
482	assert face.Area0 > min_error_area, "Face initial area is too low"	×
483
484	# Test Tris properties:
485	# Conditions checked:
486	# - Edge length > minimum error edge length
487	# - Any Lengths to Centre > minimum error edge length
488	# - Area > minimum error area
489	for c_cell in self.geo.Cells:	×
490	if c_cell.AliveStatus:	×
491	for face in c_cell.Faces:	×
492	for tris in face.Tris:	×
493	assert tris.EdgeLength > min_error_edge, "Triangle edge length is too low"	×
494	assert any(length > min_error_edge for length in	×
495	tris.LengthsToCentre), "Triangle lengths to centre are too low"
496	assert tris.Area > min_error_area, "Triangle area is too low"	×
497
UNCOV 498	def analyse_vertex_model(self):	×
499	"""
500	Analyse the vertex model.
501	:return:
502	"""
503	# Initialize average cell properties
504	cell_features = []	×
505	debris_features = []	×
506
507	wound_centre, debris_cells = self.geo.compute_wound_centre()	×
508	list_of_cell_distances = self.geo.compute_cell_distance_to_wound(debris_cells, location_filter=None)	×
509	list_of_cell_distances_top = self.geo.compute_cell_distance_to_wound(debris_cells, location_filter=0)	×
510	list_of_cell_distances_bottom = self.geo.compute_cell_distance_to_wound(debris_cells, location_filter=2)	×
511
512	# Analyse the alive cells
513	for cell_id, cell in enumerate(self.geo.Cells):	×
514	if cell.AliveStatus:	×
515	cell_features.append(cell.compute_features(wound_centre))	×
516	elif cell.AliveStatus is not None:	×
517	debris_features.append(cell.compute_features())	×
518
519	# Calculate average of cell features
520	all_cell_features = pd.DataFrame(cell_features)	×
521	all_cell_features["cell_distance_to_wound"] = list_of_cell_distances	×
522	all_cell_features["cell_distance_to_wound_top"] = list_of_cell_distances_top	×
523	all_cell_features["cell_distance_to_wound_bottom"] = list_of_cell_distances_bottom	×
524	all_cell_features["time"] = self.t	×
525	avg_cell_features = all_cell_features.mean()	×
526
527	# Compute wound features
528	try:	×
529	wound_features = self.compute_wound_features()	×
530	avg_cell_features = pd.concat([avg_cell_features, pd.Series(wound_features)])	×
531	except Exception as e:	×
532	pass	×
533
534	return all_cell_features, avg_cell_features	×
535
UNCOV 536	def compute_wound_features(self):	×
537	"""
538	Compute wound features.
539	:return:
540	"""
541	wound_features = {	×
542	'num_cells_wound_edge': len(self.geo.compute_cells_wound_edge()),
543	'num_cells_wound_edge_top': len(self.geo.compute_cells_wound_edge(location_filter="Top")),
544	'num_cells_wound_edge_bottom': len(self.geo.compute_cells_wound_edge(location_filter="Bottom")),
545	'wound_area_top': self.geo.compute_wound_area(location_filter="Top"),
546	'wound_area_bottom': self.geo.compute_wound_area(location_filter="Bottom"),
547	'wound_volume': self.geo.compute_wound_volume(),
548	'wound_height': self.geo.compute_wound_height(),
549	'wound_aspect_ratio_top': self.geo.compute_wound_aspect_ratio(location_filter="Top"),
550	'wound_aspect_ratio_bottom': self.geo.compute_wound_aspect_ratio(location_filter="Bottom"),
551	'wound_perimeter_top': self.geo.compute_wound_perimeter(location_filter="Top"),
552	'wound_perimeter_bottom': self.geo.compute_wound_perimeter(location_filter="Bottom")
553	}
554
555	return wound_features	×
556
UNCOV 557	def screenshot(self, temp_dir, selected_cells=None):	×
558	"""
559	Create a screenshot of the current state of the model.
560	:param selected_cells:
561	:param temp_dir:
562	:return:
563	"""
564	# if exists variable export_images in set
565	if hasattr(self.set, 'export_images'):	×
566	if self.set.export_images is False:	×
567	return	×
568
569	total_real_cells = len([cell.ID for cell in self.geo.Cells if cell.AliveStatus is not None])	×
570
571	# Create a colormap_lim
572	if self.colormap_lim is None:	×
573	self.colormap_lim = [0, total_real_cells]	×
574
575	# Create a plotter
576	if selected_cells is None:	×
577	selected_cells = []	×
578
579	plotter = pv.Plotter(off_screen=True)	×
580	for _, cell in enumerate(self.geo.Cells):	×
581	if cell.AliveStatus == 1 and (cell.ID in selected_cells or selected_cells is not []):	×
582	# Load the VTK file as a pyvista mesh
583	mesh = cell.create_pyvista_mesh()	×
584
585	# Add the mesh to the plotter
586	plotter.add_mesh(mesh, scalars='ID', lighting=True, cmap="prism", clim=self.colormap_lim,	×
587	show_edges=True, edge_opacity=0.5, edge_color='grey')
588	# Set a fixed camera zoom level
589	fixed_zoom_level = 1	×
590	plotter.camera.zoom(fixed_zoom_level)	×
591
592	# Add text to the plotter
593	if self.set.ablation:	×
594	timeAfterAblation = float(self.t) - float(self.set.TInitAblation)	×
595	text_content = f"Ablation time: {timeAfterAblation:.2f}"	×
596	plotter.add_text(text_content, position='upper_right', font_size=12, color='black')	×
597	else:
598	text_content = f"Time: {self.t:.2f}"	×
599	plotter.add_text(text_content, position='upper_right', font_size=12, color='black')	×
600
601	# Render the scene and capture a screenshot
602	img = plotter.screenshot()	×
603	# Save the image to a temporary file
604	temp_file = os.path.join(temp_dir, f'vModel_perspective_{self.numStep}.png')	×
605	imageio.imwrite(temp_file, img)	×
606
607	# True 2D
608	plotter.enable_parallel_projection()	×
609	plotter.enable_image_style()	×
610
611	# Set the camera to the top view
612	plotter.view_xy()	×
613
614	img = plotter.screenshot()	×
615	temp_file = os.path.join(temp_dir, f'vModel_top_{self.numStep}.png')	×
616	imageio.imwrite(temp_file, img)	×
617
618	# Set the camera to the front view
619	plotter.view_xz()	×
620
621	img = plotter.screenshot()	×
622	temp_file = os.path.join(temp_dir, f'vModel_front_{self.numStep}.png')	×
623	imageio.imwrite(temp_file, img)	×
624
625	# Set the camera to the bottom view
626	plotter.view_xy(negative=True)	×
627
628	img = plotter.screenshot()	×
629	temp_file = os.path.join(temp_dir, f'vModel_bottom_{self.numStep}.png')	×
630	imageio.imwrite(temp_file, img)	×
631
632	# Close the plotter
633	plotter.close()	×
634
UNCOV 635	def copy(self):	×
636	"""
637	Copy the VertexModel object.
638	:return:
639	"""
640	new_v_model = VertexModel()	×
641	copy_non_mutable_attributes(self, '', new_v_model)	×
642
643	return new_v_model	×
644
UNCOV 645	def calculate_error(self, K, initial_recoil, error_type=None):	×
646	"""
647	Calculate the error of the model.
648	:return:
649	"""
650	# The error consist on:
651	# - There shouldn't be any cells with very small area in the top or bottom domain.
652	# - It should get until the end of the simulation (tend).
653	# - When ablating, it should get to: 165 percentage of area at 35.8 minutes.
654	error = 0	×
655
656	# Check if the simulation reached the end
657	if self.t < self.set.tend:	×
658	error += (self.t - self.set.tend) ** 2	×
659
660	# # Check how many cells have a very small area
661	if error_type == 'None' or 'SmallArea' in error_type:	×
662	std_area_top = np.std([cell.compute_area(location_filter=0) for cell in self.geo.Cells if cell.AliveStatus == 1])	×
663	std_area_bottom = np.std([cell.compute_area(location_filter=2) for cell in self.geo.Cells if cell.AliveStatus == 1])	×
664	mean_area_top = np.mean([cell.compute_area(location_filter=0) for cell in self.geo.Cells if cell.AliveStatus == 1])	×
665	mean_area_bottom = np.mean([cell.compute_area(location_filter=2) for cell in self.geo.Cells if cell.AliveStatus == 1])	×
666	zscore_area_top = std_area_top / mean_area_top	×
667	zscore_area_bottom = std_area_bottom / mean_area_bottom	×
668	error += zscore_area_top ** 2	×
669	error += zscore_area_bottom ** 2	×
670
671	# Check how similar the recoil from in vivo is to the initial recoil and K value
672	correct_K = 0.126	×
673	correct_initial_recoil = 0.213	×
674	if error_type == 'None' or 'K' in error_type:	×
675	try:	×
676	error += np.abs(K[0] - correct_K) * 100	×
677	except IndexError:	×
678	error += np.abs(K - correct_K) * 100	×
679
680	if error_type == 'None' or 'InitialRecoil' in error_type:	×
681	try:	×
682	error += np.abs(initial_recoil[0] - correct_initial_recoil) * 100	×
683	except IndexError:	×
684	error += np.abs(initial_recoil - correct_initial_recoil) * 100	×
685
686	return error	×

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