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

UNCOV 1	import logging	×
UNCOV 2	import math	×
UNCOV 3	import statistics	×
4
UNCOV 5	import numpy as np	×
UNCOV 6	from scipy.optimize import minimize	×
UNCOV 7	from scipy.spatial import Delaunay	×
8
UNCOV 9	from src.pyVertexModel.algorithm.vertexModel import VertexModel	×
10
UNCOV 11	logger = logging.getLogger("pyVertexModel")	×
12
13
UNCOV 14	def AreTri(p1, p2, p3):	×
15	return 0.5 * np.linalg.norm(np.cross(p2 - p1, p3 - p1))	×
16
17
UNCOV 18	def check_replicateed_nodes(X, nX, h):	×
19	ToBeRemoved = np.zeros(nX.shape[0], dtype=bool)	×
20	for jj in range(nX.shape[0]):	×
21	m = np.linalg.norm(X - nX[jj], axis=1)	×
22	m = np.min(m)	×
23	if m < 1e-2 * h:	×
24	ToBeRemoved[jj] = True	×
25	nX = nX[~ToBeRemoved]	×
26	return nX	×
27
28
UNCOV 29	def SeedNodeTet(X, XgID, Twgi, h):	×
30	XTet = X[Twgi, :]	×
31	Center = np.mean(XTet, axis=0)	×
32	nX = np.zeros((4, 3))	×
33	for i in range(4):	×
34	vc = Center - XTet[i, :]	×
35	dis = np.linalg.norm(vc)	×
36	dir = vc / dis	×
37	offset = h * dir	×
38	if dis > np.linalg.norm(offset):	×
39	nX[i, :] = XTet[i, :] + offset	×
40	else:
41	nX[i, :] = XTet[i, :] + vc	×
42
43	mask = np.isin(Twgi, XgID)	×
44	nX = nX[~mask, :]	×
45	nX = np.unique(nX, axis=0)	×
46	nX = check_replicateed_nodes(X, nX, h)	×
47	nXgID = np.arange(X.shape[0], X.shape[0] + nX.shape[0])	×
48	X = np.vstack((X, nX))	×
49	XgID = np.concatenate((XgID, nXgID))	×
50	return X, XgID	×
51
52
UNCOV 53	def SeedNodeTri(X, XgID, Tri, h):	×
54	XTri = X[Tri, :]	×
55	Center = np.mean(XTri, axis=0)	×
56	nX = np.zeros((3, 3))	×
57	for i in range(3):	×
58	vc = Center - XTri[i, :]	×
59	dis = np.linalg.norm(vc)	×
60	dir = vc / dis	×
61	offset = h * dir	×
62	if dis > np.linalg.norm(offset):	×
63	nX[i, :] = XTri[i, :] + offset	×
64	else:
65	nX[i, :] = XTri[i, :] + vc	×
66
67	mask = np.isin(Tri, XgID)	×
68	nX = nX[~mask, :]	×
69	nX = np.unique(nX, axis=0)	×
70	nX = check_replicateed_nodes(X, nX, h)	×
71	nXgID = np.arange(X.shape[0], X.shape[0] + nX.shape[0])	×
72	X = np.vstack((X, nX))	×
73	XgID = np.concatenate((XgID, nXgID))	×
74	return X, XgID	×
75
76
UNCOV 77	def delaunay_compute_entities(tris, X, XgID, XgIDBB, nCells, s):	×
78	# Initialize variables
79	Side = np.array([[0, 1, 2], [0, 1, 3], [1, 2, 3], [0, 2, 3]])	×
80	Edges = np.array([[0, 1], [1, 2], [0, 2], [0, 3], [1, 3], [2, 3]])	×
81	Vol = np.zeros(tris.shape[0])	×
82	AreaFaces = np.zeros((tris.shape[0] * 3, 4))	×
83	LengthEdges = np.zeros((tris.shape[0] * 3, 6))	×
84	Arc = 0	×
85	Lnc = 0	×
86
87	# Compute volume, area and length of each tetrahedron
88	for i in range(tris.shape[0]):	×
89	for j in range(4):	×
90	if np.sum(np.isin(tris[i, Side[j]], XgID)) == 0:	×
91	p1, p2, p3 = X[tris[i, Side[j]]]	×
92	AreaFaces[i, j] = AreTri(p1, p2, p3)	×
93	Arc += 1	×
94
95	for j in range(6):	×
96	if np.sum(np.isin(tris[i, Edges[j]], XgID)) == 0:	×
97	p1, p2 = X[tris[i, Edges[j]]]	×
98	LengthEdges[i, j] = np.linalg.norm(p1 - p2)	×
99	Lnc += 1	×
100
101	# Seed nodes in big entities (based on characteristic Length h)
102	for i in range(tris.shape[0]):	×
103	for j in range(4):	×
104	if np.sum(np.isin(tris[i, Side[j]], XgID)) == 0 and AreaFaces[i, j] > s ** 2:	×
105	X, XgID = SeedNodeTri(X, XgID, tris[i, Side[j]], s)	×
106
107	for j in range(6):	×
108	if np.sum(np.isin(tris[i, Edges[j]], XgID)) == 0 and LengthEdges[i, j] > 2 * s:	×
109	X, XgID = SeedNodeTet(X, XgID, tris[i], s)	×
110	break	×
111
112	# Seed on ghost tetrahedra
113	for i in range(len(Vol)):	×
114	if np.sum(np.isin(tris[i], XgID)) > 0:	×
115	X, XgID = SeedNodeTet(X, XgID, tris[i], s)	×
116
117	X = np.delete(X, XgIDBB, axis=0)	×
118	XgID = np.arange(nCells, X.shape[0])	×
119
120	return X, XgID	×
121
122
UNCOV 123	def generate_points_in_sphere(total_cells):	×
124	"""
125	Generate points in a sphere
126	:param total_cells: The total number of cells
127	:return: The X, Y, Z coordinates of the points
128	"""
129	r_unit = 1	×
130
131	# Calculating area, distance, and increments for theta and phi
132	Area = 4 * math.pi * r_unit ** 2 / total_cells	×
133	Distance = math.sqrt(Area)	×
134	M_theta = round(math.pi / Distance)	×
135	d_theta = math.pi / M_theta	×
136	d_phi = Area / d_theta	×
137
138	# Initializing lists for X, Y, Z coordinates
139	X, Y, Z = [], [], []	×
140	N_new = 0	×
141
142	for m in range(M_theta):	×
143	Theta = math.pi * (m + 0.5) / M_theta	×
144	M_phi = round(2 * math.pi * math.sin(Theta) / d_phi)	×
145
146	for n in range(M_phi):	×
147	Phi = 2 * math.pi * n / M_phi	×
148
149	# Updating node count
150	N_new += 1	×
151
152	# Calculating and appending coordinates
153	X.append(math.sin(Theta) * math.cos(Phi))	×
154	Y.append(math.sin(Theta) * math.sin(Phi))	×
155	Z.append(math.cos(Theta))	×
156
157	return X, Y, Z, N_new	×
158
159
UNCOV 160	def generate_first_ghost_nodes(X):	×
161	# Bounding Box 1
162	nCells = X.shape[0]	×
163	r0 = np.average(X, axis=0)	×
164	r0[0] = statistics.mean(X[:, 0])	×
165	r0[1] = statistics.mean(X[:, 1])	×
166	r0[2] = statistics.mean(X[:, 2])	×
167
168	r = 5 * np.max(np.abs(X - r0))	×
169	# Define bounding nodes: bounding sphere
170	theta = np.linspace(0, 2 * np.pi, 5)	×
171	phi = np.linspace(0, np.pi, 5)	×
172	theta, phi = np.meshgrid(theta, phi, indexing='ij') # Ensure the order matches MATLAB	×
173	# Phi and Theta should be transpose as it is in Matlab
174	phi = phi.T	×
175	theta = theta.T	×
176
177	# Convert to Cartesian coordinates
178	x = r * np.sin(phi) * np.cos(theta)	×
179	y = r * np.sin(phi) * np.sin(theta)	×
180	z = r * np.cos(phi)	×
181	# Reshape to column vectors, ensuring the same order as MATLAB
182	x = x.flatten('F')	×
183	y = y.flatten('F')	×
184	z = z.flatten('F')	×
185	# Offset the points by r0 and combine into a single array
186	Xg = np.column_stack((x, y, z)) + r0	×
187	# Find unique values considering the tolerance
188	tolerance = 1e-6	×
189	_, idx = np.unique(Xg.round(decimals=int(-np.log10(tolerance))), axis=0, return_index=True)	×
190	Xg = Xg[idx]	×
191
192	# Add new bounding nodes to X
193	XgID = np.arange(nCells, nCells + Xg.shape[0])	×
194	XgIDBB = XgID.copy()	×
195	X = np.vstack((X, Xg))	×
196	return X, XgID, XgIDBB, nCells	×
197
198
UNCOV 199	def build_topo(c_set, nx=None, ny=None, nz=None, columnar_cells=False):	×
200	"""
201	This function builds the topology of the mesh.
202	:param nx: Number of nodes in x direction
203	:param ny: Number of nodes in y direction
204	:param nz: Number of nodes in z direction
205	:param c_set: Set class
206	:param columnar_cells: Boolean to indicate if the cells are columnar
207	:return: X: Nodal positions
208	X_Ids: Nodal IDs
209	"""
210	X = np.empty((0, 3))	×
211	X_Ids = []	×
212	if c_set.InputGeo == 'Bubbles':	×
213	for numZ in range(nz):	×
214	x = np.arange(nx)	×
215	y = np.arange(ny)	×
216	x, y = np.meshgrid(x, y, indexing='ij')	×
217	x = x.flatten('F')	×
218	y = y.flatten('F')	×
219	z = np.ones_like(x) * numZ	×
220	X = np.vstack((X, np.column_stack((x, y, z))))	×
221
222	if columnar_cells:	×
223	X_Ids.append(np.arange(len(x)))	×
224	else:
225	X_Ids = np.arange(X.shape[0])	×
226
227	elif c_set.InputGeo == 'Bubbles_Cyst':	×
228	X, Y, Z, _ = generate_points_in_sphere(c_set.TotalCells)	×
229
230	X = np.array([X, Y, Z]).T * 10	×
231
232	# Lumen as the first cell
233	lumenCell = np.mean(X, axis=0)	×
234	X = np.vstack([lumenCell, X])	×
235	c_set.TotalCells = X.shape[0]	×
236
237	return X, X_Ids	×
238
239
UNCOV 240	def SeedWithBoundingBox(X, s):	×
241	"""
242	This function seeds nodes in desired entities (edges, faces and tetrahedrons) while cell-centers are bounded
243	by ghost nodes.
244	:param X:
245	:param s:
246	:return:
247	"""
248
249	X, XgID, XgIDBB, nCells = generate_first_ghost_nodes(X)	×
250
251	N = 3 # The dimensions of our points	×
252	options = 'Qt Qbb Qc' if N <= 3 else 'Qt Qbb Qc Qx' # Set the QHull options	×
253	Tri = Delaunay(X, qhull_options=options)	×
254
255	# first Delaunay with ghost nodes
256	X, XgID = delaunay_compute_entities(Tri.simplices, X, XgID, XgIDBB, nCells, s)	×
257	return XgID, X	×
258
259
UNCOV 260	def fit_ellipsoid_to_points(points):	×
261	"""
262	Fit an ellipsoid to a set of points using the least-squares method
263	:param points:
264	:return:
265	"""
266	# Extract coordinates from the input array
267	x, y, z = points[:, 0], points[:, 1], points[:, 2]	×
268
269	# Define the objective function for ellipsoid fitting
270	def ellipsoidError(c_points):	×
271	"""
272	Calculate the sum of squared distances from the ellipsoid surface to the input points
273	:param c_points: The input points
274	:return: The sum of squared distances from the ellipsoid surface to the input points
275	"""
276	a, b, c = c_points	×
277	distances = (x 2 / a 2) + (y 2 / b 2) + (z 2 / c 2) - 1	×
278	error = np.sum(distances ** 2)	×
279	return error	×
280
281	# Initial guess for the semi-axis lengths
282	initialGuess = np.array([np.std(x), np.std(y), np.std(z)])	×
283
284	# Perform optimization to find the best-fitting ellipsoid parameters
285	result = minimize(ellipsoidError, x0=initialGuess, method='BFGS')	×
286
287	# Extract optimized parameters and normalize
288	paramsOptimized = result.x	×
289	a, b, c = paramsOptimized / np.max(paramsOptimized)	×
290
291	return abs(a), abs(b), abs(c), abs(paramsOptimized)	×
292
293
UNCOV 294	def extrapolate_ys_faces_ellipsoid(geo, c_set):	×
295	"""
296	Extrapolate the vertices of the cells to the ellipsoid
297	:param geo:
298	:param c_set:
299	:return:
300	"""
301	# Original axis values
302	Ys_top = np.concatenate([cell.Y for cell in geo.Cells[1:c_set.TotalCells]])	×
303
304	a, b, c, paramsOptimized_top = fit_ellipsoid_to_points(Ys_top)	×
305	a, b, c, paramsOptimized_bottom = fit_ellipsoid_to_points(geo.Cells[0].Y)	×
306
307	# Normalised based on those
308	ellipsoid_axis_normalised1 = c_set.ellipsoid_axis1 / paramsOptimized_top[0]	×
309	ellipsoid_axis_normalised2 = c_set.ellipsoid_axis2 / paramsOptimized_top[1]	×
310	ellipsoid_axis_normalised3 = c_set.ellipsoid_axis3 / paramsOptimized_top[2]	×
311	lumen_axis_normalised1 = c_set.lumen_axis1 / paramsOptimized_bottom[0]	×
312	lumen_axis_normalised2 = c_set.lumen_axis2 / paramsOptimized_bottom[1]	×
313	lumen_axis_normalised3 = c_set.lumen_axis3 / paramsOptimized_bottom[2]	×
314
315	# Extrapolate top layer as the outer ellipsoid, the bottom layer as the lumen, and lateral is rebuilt.
316	allTs = np.unique(np.sort(np.concatenate([cell.T for cell in geo.Cells[:c_set.TotalCells]]), axis=1), axis=0)	×
317	topTs = allTs[np.any(np.isin(allTs, geo.XgTop), axis=1)]	×
318	bottomsTs = allTs[np.any(np.isin(allTs, geo.XgBottom), axis=1)]	×
319
320	# Changes vertices of other cells
321	for tetToCheck in topTs:	×
322	for nodeInTet in tetToCheck:	×
323	if (nodeInTet not in geo.XgTop and geo.Cells[nodeInTet] is not None and	×
324	geo.Cells[nodeInTet].Y is not None):
325	newPoint = geo.Cells[nodeInTet].Y[	×
326	np.all(np.isin(geo.Cells[nodeInTet].T, tetToCheck), axis=1)]
327	newPoint_extrapolated = extrapolate_points_to_ellipsoid(newPoint, ellipsoid_axis_normalised1,	×
328	ellipsoid_axis_normalised2,
329	ellipsoid_axis_normalised3)
330	geo.Cells[nodeInTet].Y[	×
331	np.all(np.isin(geo.Cells[nodeInTet].T, tetToCheck), axis=1)] = newPoint_extrapolated
332
333	for tetToCheck in bottomsTs:	×
334	for nodeInTet in tetToCheck:	×
335	if (nodeInTet not in geo.XgTop and geo.Cells[nodeInTet] is not None and	×
336	geo.Cells[nodeInTet].Y is not None):
337	newPoint = geo.Cells[nodeInTet].Y[	×
338	np.all(np.isin(geo.Cells[nodeInTet].T, tetToCheck), axis=1)]
339	newPoint_extrapolated = extrapolate_points_to_ellipsoid(newPoint, lumen_axis_normalised1,	×
340	lumen_axis_normalised2,
341	lumen_axis_normalised3)
342	geo.Cells[nodeInTet].Y[	×
343	np.all(np.isin(geo.Cells[nodeInTet].T, tetToCheck), axis=1)] = newPoint_extrapolated
344
345	# Recalculating face centres here based on the previous changes
346	geo.rebuild(geo.copy(), c_set)	×
347	geo.build_global_ids()	×
348	geo.update_measures()	×
349	for cell in geo.Cells:	×
350	cell.Area0 = c_set.cell_A0	×
351	cell.Vol0 = c_set.cell_V0	×
352	geo.Cells[0].Area0 = c_set.lumen_V0 * (c_set.cell_A0 / c_set.cell_V0)	×
353	geo.Cells[0].Vol0 = c_set.lumen_V0	×
354
355	# Calculate the mean volume excluding the first cell
356	meanVolume = np.mean([cell.Vol for cell in geo.Cells[1:c_set.TotalCells]])	×
357	logger.info(f'Average Cell Volume: {meanVolume}')	×
358	# Calculate the standard deviation of volumes excluding the first cell
359	stdVolume = np.std([cell.Vol for cell in geo.Cells[1:c_set.TotalCells]])	×
360	logger.info(f'Standard Deviation of Cell Volumes: {stdVolume}')	×
361	# Display the volume of the first cell
362	firstCellVolume = geo.Cells[0].Vol	×
363	logger.info(f'Volume of Lumen: {firstCellVolume}')	×
364	# Calculate the sum of volumes excluding the first cell
365	sumVolumes = np.sum([cell.Vol for cell in geo.Cells[1:c_set.TotalCells]])	×
366	logger.info(f'Tissue Volume: {sumVolumes}')	×
367
368	return geo	×
369
370
UNCOV 371	def extrapolate_points_to_ellipsoid(points, ellipsoid_axis_normalised1, ellipsoid_axis_normalised2,	×
372	ellipsoid_axis_normalised3):
373	points[:, 0] = points[:, 0] * ellipsoid_axis_normalised1	×
374	points[:, 1] = points[:, 1] * ellipsoid_axis_normalised2	×
375	points[:, 2] = points[:, 2] * ellipsoid_axis_normalised3	×
376
377	return points	×
378
379
UNCOV 380	class VertexModelBubbles(VertexModel):	×
UNCOV 381	def __init__(self, set_test=None):	×
382	super().__init__(set_test)	×
383
UNCOV 384	def initialize(self):	×
385	"""
386	Initialize the geometry and the topology of the model.
387	:return:
388	"""
389	# Build nodal mesh
390	self.generate_Xs(self.geo.nx, self.geo.ny, self.geo.nz)	×
391
392	# This code is to match matlab's output and python's
393	# N = 3 # The dimensions of our points
394	# options = 'Qt Qbb Qc' if N <= 3 else 'Qt Qbb Qc Qx' # Set the QHull options
395	Twg = Delaunay(self.X).simplices	×
396
397	# Remove tetrahedras formed only by ghost nodes
398	Twg = Twg[~np.all(np.isin(Twg, self.geo.XgID), axis=1)]	×
399	# Remove weird IDs
400
401	# Re-number the surviving tets
402	uniqueTets, indices = np.unique(Twg, return_inverse=True)	×
403	self.geo.XgID = np.arange(self.geo.nCells, len(uniqueTets))	×
404	self.X = self.X[uniqueTets]	×
405	Twg = indices.reshape(Twg.shape)	×
406
407	if self.set.InputGeo == 'Bubbles_Cyst':	×
408	self.geo.XgBottom = [0]	×
409	self.geo.XgTop = self.geo.XgID	×
410	self.geo.XgID = np.append(self.geo.XgID, 0)	×
411	else:
412	Xg = self.X[self.geo.XgID]	×
413	self.geo.XgBottom = self.geo.XgID[Xg[:, 2] < np.mean(self.X[:, 2])]	×
414	self.geo.XgTop = self.geo.XgID[Xg[:, 2] > np.mean(self.X[:, 2])]	×
415
416	self.geo.Main_cells = range(len(self.geo.nCells))	×
417	self.geo.build_cells(self.set, self.X, Twg)	×
418
419	if self.set.InputGeo == 'Bubbles_Cyst':	×
420	# Extrapolate Face centres and Ys to the ellipsoid
421	self.geo = extrapolate_ys_faces_ellipsoid(self.geo, self.set)	×
422
423	# Define upper and lower area threshold for remodelling
424	self.geo.init_reference_cell_values(self.set)	×
425
UNCOV 426	def generate_Xs(self, nx=None, ny=None, nz=None):	×
427	"""
428	Generate the nodal positions of the mesh based on the input geometry
429	:return:
430	"""
431	self.X, X_IDs = build_topo(self.set, nx, ny, nz)	×
432	self.geo.nCells = self.X.shape[0]	×
433	# Centre Nodal position at (0,0)
434	self.X[:, 0] = self.X[:, 0] - np.mean(self.X[:, 0])	×
435	self.X[:, 1] = self.X[:, 1] - np.mean(self.X[:, 1])	×
436	self.X[:, 2] = self.X[:, 2] - np.mean(self.X[:, 2])	×
437
438	if self.set.InputGeo == 'Bubbles_Cyst':	×
439	a, b, c, paramsOptimized = fit_ellipsoid_to_points(self.X)	×
440
441	ellipsoid_axis_normalised1 = np.mean([self.set.ellipsoid_axis1, self.set.lumen_axis1]) / paramsOptimized[0]	×
442	ellipsoid_axis_normalised2 = np.mean([self.set.ellipsoid_axis2, self.set.lumen_axis2]) / paramsOptimized[1]	×
443	ellipsoid_axis_normalised3 = np.mean([self.set.ellipsoid_axis3, self.set.lumen_axis3]) / paramsOptimized[2]	×
444
445	# Extrapolate Xs
446	self.X = extrapolate_points_to_ellipsoid(self.X, ellipsoid_axis_normalised1, ellipsoid_axis_normalised2,	×
447	ellipsoid_axis_normalised3)
448	# Perform Delaunay
449	self.geo.XgID, self.X = SeedWithBoundingBox(self.X, self.set.s)	×
450	if self.set.Substrate == 1:	×
451	Xg = self.X[self.geo.XgID, :]	×
452	self.X = np.delete(self.X, self.geo.XgID, 0)	×
453	Xg = Xg[Xg[:, 2] > np.mean(self.X[:, 2]), :]	×
454	self.geo.XgID = np.arange(self.X.shape[0], self.X.shape[0] + Xg.shape[0] + 2)	×
455	self.X = np.concatenate((self.X, Xg, [np.mean(self.X[:, 0]), np.mean(self.X[:, 1]), -50]), axis=0)	×
456
UNCOV 457	def copy(self):	×
458	"""
459	Copy the object
460	:return:
461	"""
462	return super().copy()	×

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