11368168500

Committed 16 Oct 2024 02:48PM UTC coverage: 0.794% (+0.001%) from 0.793%

Build # 11368168500

Build Type

push

github

Committed by

Pablo1990

Commit Message

refactor to obtain tris_areas

Run Details

0 of 1003 branches covered (0.0%)

Branch coverage included in aggregate %.

0 of 2 new or added lines in 1 file covered. (0.0%)

113 existing lines in 3 files now uncovered.

51 of 5422 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/geometry/cell.py

import numpy as np
import pyvista as pv
import vtk
from sklearn.decomposition import PCA

from src.pyVertexModel.geometry import face
from src.pyVertexModel.geometry.face import get_interface
from src.pyVertexModel.util.utils import copy_non_mutable_attributes


def face_centres_to_middle_of_neighbours_vertices(Geo, c_cell):
    """
    Move the face centres to the middle of the neighbours vertices.
    :param Geo:
    :param c_cell:
    :return:
    """
    for num_face, _ in enumerate(Geo.Cells[c_cell].Faces):
        all_edges = []
        for tri in Geo.Cells[c_cell].Faces[num_face].Tris:
            all_edges.append(tri.Edge)

        all_edges = np.unique(np.concatenate(all_edges))
        Geo.Cells[c_cell].Faces[num_face].Centre = np.mean(
            Geo.Cells[c_cell].Y[all_edges, :], axis=0)


def compute_2d_circularity(area, perimeter):
    """
    Compute the 2D circularity of the cell
    :return:
    """
    if perimeter == 0:
        return 0
    else:
        return 4 * np.pi * area / perimeter ** 2


def compute_y(geo, T, cellCentre, Set):
    """
    Compute the new Y of the cell

    :param geo:
    :param T:
    :param cellCentre:
    :param Set:
    :return:
    """
    x = [geo.Cells[i].X for i in T]
    newY = np.mean(x, axis=0)
    if sum([geo.Cells[i].AliveStatus is not None for i in T]) == 1 and "Bubbles" in Set.InputGeo:
        vc = newY - cellCentre
        dir = vc / np.linalg.norm(vc)
        offset = Set.f * dir
        newY = cellCentre + offset

    if "Bubbles" not in Set.InputGeo:
        if any(i in geo.XgTop for i in T):
            newY[2] /= sum(i in geo.XgTop for i in T) / 2
        elif any(i in geo.XgBottom for i in T):
            newY[2] /= sum(i in geo.XgBottom for i in T) / 2
    return newY


class Cell:
    """
    Class that contains the information of a cell.
    """

    def __init__(self, mat_file=None):
        """
        Constructor of the class
        :param mat_file:
        """
        self.opposite_cell = None
        self.lambda_s1_noise = None
        self.lambda_s2_noise = None
        self.lambda_s3_noise = None
        self.lambda_v_noise = None
        self.barrier_tri0_top = None
        self.barrier_tri0_bottom = None
        self.contractility_noise = None
        self.SubstrateLambda = None
        self.InternalLambda = None
        self.ExternalLambda = None
        self.axes_lengths = None
        self.Y = None
        self.globalIds = np.array([], dtype='int')
        self.Faces = []
        self.Area = None
        self.Area0 = None
        self.Vol = None
        self.Vol0 = None
        self.AliveStatus = None
        self.vertices_and_faces_to_remodel = np.array([], dtype='int')
        # TODO: Save contractile forces (g) to output
        self.substrate_g = None
        self.lambdaB_perc = 1

        if mat_file is None:
            self.ID = None
            self.X = None
            self.T = None
        else:
            self.ID = mat_file[0][0][0] - 1
            self.X = np.array(mat_file[1][0], dtype=np.float64)
            self.T = mat_file[2] - 1

            if len(mat_file[4]) > 0:
                self.Y = np.array(mat_file[3], dtype=np.float64)
                for c_face in mat_file[4][0]:
                    self.Faces.append(face.Face(c_face))
                if len(mat_file[5]) > 0:
                    self.Vol = mat_file[5][0][0]
                    self.Vol0 = mat_file[6][0][0]
                    self.Area = mat_file[7][0][0]
                    self.Area0 = mat_file[8][0][0]
                    self.globalIds = np.concatenate(mat_file[9]) - 1
                    self.cglobalids = mat_file[10][0][0] - 1
                    self.AliveStatus = mat_file[11][0][0]

    def copy(self):
        """
        Copy the cell
        :return:
        """
        new_cell = Cell()

        copy_non_mutable_attributes(self, 'Faces', new_cell)

        new_cell.Faces = [f.copy() for f in self.Faces]

        return new_cell

    def compute_area(self, location_filter=None):
        """
        Compute the area of the cell
        :param location_filter:
        :return:
        """
        total_area = 0.0
        for f in range(len(self.Faces)):
            if location_filter is not None:
                if get_interface(self.Faces[f].InterfaceType) == get_interface(location_filter):
                    total_area = total_area + self.Faces[f].Area
            else:
                total_area = total_area + self.Faces[f].Area

        return total_area

    def compute_volume(self):
        """
        Compute the volume of the cell
        :return:
        """
        v = 0.0
        for f in range(len(self.Faces)):
            c_face = self.Faces[f]
            for t in range(len(c_face.Tris)):
                y1 = self.Y[c_face.Tris[t].Edge[0], :] - self.X
                y2 = self.Y[c_face.Tris[t].Edge[1], :] - self.X
                y3 = c_face.Centre - self.X
                ytri = np.array([y1, y2, y3])

                current_v = np.linalg.det(ytri) / 6
                # If the volume is negative, switch two the other option
                if current_v < 0:
                    ytri = np.array([y2, y1, y3])
                    current_v = np.linalg.det(ytri) / 6

                v += current_v

        self.Vol = v
        return v

    def create_vtk(self):
        """
        Create a vtk cell
        :return:
        """
        points = vtk.vtkPoints()
        points.SetNumberOfPoints(len(self.Y) + len(self.Faces))
        for i in range(len(self.Y)):
            points.SetPoint(i, self.Y[i, 0], self.Y[i, 1], self.Y[i, 2])

        cell = vtk.vtkCellArray()
        # Go through all the faces and create the triangles for the VTK cell
        total_tris = 0
        for f in range(len(self.Faces)):
            c_face = self.Faces[f]
            points.SetPoint(len(self.Y) + f, c_face.Centre[0], c_face.Centre[1], c_face.Centre[2])
            for t in range(len(c_face.Tris)):
                cell.InsertNextCell(3)
                cell.InsertCellPoint(c_face.Tris[t].Edge[0])
                cell.InsertCellPoint(c_face.Tris[t].Edge[1])
                cell.InsertCellPoint(len(self.Y) + f)

            total_tris += len(c_face.Tris)

        vpoly = vtk.vtkPolyData()
        vpoly.SetPoints(points)

        vpoly.SetPolys(cell)

        # Get all the different properties of a cell
        properties = self.compute_features()

        # Go through the different properties of the dictionary and add them to the vtkFloatArray
        for key, value in properties.items():
            # Create a vtkFloatArray to store the properties of the cell
            property_array = vtk.vtkFloatArray()
            property_array.SetName(key)

            # Add as many values as tris the cell has
            for i in range(total_tris):
                property_array.InsertNextValue(value)

            # Add the property array to the cell data
            vpoly.GetCellData().AddArray(property_array)

            if key == 'ID':
                # Default parameter
                vpoly.GetCellData().SetScalars(property_array)

        return vpoly

    def create_pyvista_mesh(self):
        """
        Create a PyVista mesh
        :return:
        """
        mesh = pv.PolyData(self.create_vtk())

        return mesh

    def compute_features(self, centre_wound=None):
        """
        Compute the features of the cell and create a dictionary with them
        :return: a dictionary with the features of the cell
        """
        features = {'ID': self.ID,
                    'Area': self.compute_area(),
                    'Area_top': self.compute_area(location_filter=0),
                    'Area_bottom': self.compute_area(location_filter=2),
                    'Area_cellcell': self.compute_area(location_filter=1),
                    'Volume': self.compute_volume(),
                    'Height': self.compute_height(),
                    'Width': self.compute_width(),
                    'Length': self.compute_length(),
                    'Perimeter': self.compute_perimeter(),
                    '2D_circularity_top': compute_2d_circularity(self.compute_area(location_filter=0),
                                                                 self.compute_perimeter(filter_location=0)),
                    '2d_circularity_bottom': compute_2d_circularity(self.compute_area(location_filter=2),
                                                                    self.compute_perimeter(filter_location=2)),
                    '2D_aspect_ratio_top': self.compute_2d_aspect_ratio(filter_location=0),
                    '2D_aspect_ratio_bottom': self.compute_2d_aspect_ratio(filter_location=2),
                    '2D_aspect_ratio_cellcell': self.compute_2d_aspect_ratio(filter_location=1),
                    '3D_aspect_ratio_0_1': self.compute_3d_aspect_ratio(),
                    '3D_aspect_ratio_0_2': self.compute_3d_aspect_ratio(axis=(0, 2)),
                    '3D_aspect_ratio_1_2': self.compute_3d_aspect_ratio(axis=(1, 2)),
                    'Sphericity': self.compute_sphericity(),
                    'Elongation': self.compute_elongation(),
                    'Ellipticity': self.compute_ellipticity(),
                    'Neighbours': len(self.compute_neighbours()),
                    'Neighbours_top': len(self.compute_neighbours(location_filter=0)),
                    'Neighbours_bottom': len(self.compute_neighbours(location_filter=2)),
                    'Tilting': self.compute_tilting(),
                    'Perimeter_top': self.compute_perimeter(filter_location=0),
                    'Perimeter_bottom': self.compute_perimeter(filter_location=2),
                    'Perimeter_cellcell': self.compute_perimeter(filter_location=1),
                    }

        if centre_wound is not None:
            features['Distance_to_wound'] = self.compute_distance_to_wound(centre_wound)

        return features

    def create_vtk_edges(self):
        """
        Create a vtk with only the information on the edges of the cell
        :return:
        """
        points = vtk.vtkPoints()
        points.SetNumberOfPoints(len(self.Y))
        for i in range(len(self.Y)):
            points.SetPoint(i, self.Y[i, 0], self.Y[i, 1], self.Y[i, 2])

        cell = vtk.vtkCellArray()
        # Go through all the faces and create the triangles for the VTK cell
        total_edges = 0
        for f in range(len(self.Faces)):
            c_face = self.Faces[f]
            for t in range(len(c_face.Tris)):
                cell.InsertNextCell(2)
                cell.InsertCellPoint(c_face.Tris[t].Edge[0])
                cell.InsertCellPoint(c_face.Tris[t].Edge[1])

            total_edges += len(c_face.Tris)

        vpoly = vtk.vtkPolyData()
        vpoly.SetPoints(points)

        vpoly.SetLines(cell)

        # Get all the different properties of a cell
        properties = self.compute_edge_features()

        # Go through the different properties of the dictionary and add them to the vtkFloatArray
        for key, value in properties[0].items():
            # Create a vtkFloatArray to store the properties of the cell
            property_array = vtk.vtkFloatArray()
            property_array.SetName(key)

            # Add as many values as tris the cell has
            for i in range(total_edges):
                if properties[i][key] is None:
                    property_array.InsertNextValue(0)
                else:
                    property_array.InsertNextValue(properties[i][key])

            # Add the property array to the cell data
            vpoly.GetCellData().AddArray(property_array)

            if key == 'ContractilityValue':
                # Default parameter
                vpoly.GetCellData().SetScalars(property_array)

        return vpoly

    def compute_edge_features(self):
        """
        Compute the features of the edges of a cell and create a dictionary with them
        :return:
        """
        cell_features = np.array([])
        for f, c_face in enumerate(self.Faces):
            for t, c_tris in enumerate(c_face.Tris):
                cell_features = np.append(cell_features, c_tris.compute_features())

        return cell_features

    def compute_height(self):
        """
        Compute the height of the cell regardless of the orientation
        :return:
        """
        return np.max(self.Y[:, 2]) - np.min(self.Y[:, 2])

    def compute_principal_axis_length(self):
        """
        Compute the principal axis length of the cell
        :return:
        """
        # Perform PCA to find the principal axes
        pca = PCA(n_components=3)
        pca.fit(self.Y)

        # Project points onto the principal axes
        projected_points = pca.transform(self.Y)

        # Calculate the lengths of the ellipsoid axes
        min_values = projected_points.min(axis=0)
        max_values = projected_points.max(axis=0)
        self.axes_lengths = max_values - min_values

        return self.axes_lengths

    def compute_width(self):
        """
        Compute the width of the cell
        :return:
        """
        return np.max(self.Y[:, 0]) - np.min(self.Y[:, 0])

    def compute_length(self):
        """
        Compute the length of the cell
        :return:
        """
        return np.max(self.Y[:, 1]) - np.min(self.Y[:, 1])

    def compute_neighbours(self, location_filter=None):
        """
        Compute the neighbours of the cell
        :return:
        """
        neighbours = []
        for f in range(len(self.Faces)):
            if location_filter is not None:
                if get_interface(self.Faces[f].InterfaceType) == get_interface(location_filter):
                    for t in range(len(self.Faces[f].Tris)):
                        neighbours.append(self.Faces[f].Tris[t].SharedByCells)
            else:
                for t in range(len(self.Faces[f].Tris)):
                    neighbours.append(self.Faces[f].Tris[t].SharedByCells)

        # Flatten the list of lists into a single 1-D array
        if len(neighbours) > 0:
            neighbours_flat = np.concatenate(neighbours)
            neighbours_unique = np.unique(neighbours_flat)
            neighbours_unique = neighbours_unique[neighbours_unique != self.ID]
        else:
            neighbours_unique = []

        return neighbours_unique

    def compute_tilting(self):
        """
        Compute the tilting of the cell
        :return:
        """
        return np.arctan(self.compute_height() / self.compute_width())

    def compute_sphericity(self):
        """
        Compute the sphericity of the cell
        :return:
        """
        return 36 * np.pi * self.Vol ** 2 / self.Area ** 3

    def compute_perimeter(self, filter_location=None):
        """
        Compute the perimeter of the cell
        :return:
        """
        perimeter = 0.0
        for f in range(len(self.Faces)):
            if filter_location is not None:
                if get_interface(self.Faces[f].InterfaceType) == get_interface(filter_location):
                    perimeter = perimeter + self.Faces[f].compute_perimeter()
            else:
                perimeter = perimeter + self.Faces[f].compute_perimeter()

        return perimeter

    def compute_ellipticity(self):
        """
        Compute the ellipticity of the cell
        :return:
        """
        return self.compute_width() / self.compute_length()

    def compute_elongation(self):
        """
        Compute the elongation of the cell
        :return:
        """
        return self.compute_length() / self.compute_height()

    def compute_3d_aspect_ratio(self, axis=(0, 1)):
        """
        Compute the 3D aspect ratio of the cell
        :return:
        """
        return self.compute_principal_axis_length()[axis[0]] / self.compute_principal_axis_length()[axis[1]]

    def compute_2d_aspect_ratio(self, filter_location=None):
        """
        Compute the 2D aspect ratio of the cell
        :return:
        """
        perimeter = self.compute_perimeter(filter_location)
        if perimeter == 0:
            return 0
        else:
            return self.compute_area(filter_location) / perimeter ** 2

    def build_y_from_x(self, geo, c_set):
        Tets = self.T
        dim = self.X.shape[0]
        Y = np.zeros((len(Tets), dim))
        for i in range(len(Tets)):
            Y[i] = compute_y(geo, Tets[i], self.X, c_set)
        return Y

    def compute_distance_to_wound(self, centre_wound):
        """
        Compute the distance from the centre of the cell to the centre of the wound
        :return:
        """
        return np.linalg.norm(self.Y - centre_wound)

    def kill_cell(self):
        """
        Kill the cell
        :return:
        """
        self.AliveStatus = None
        self.Y = None
        self.Vol = None
        self.Vol0 = None
        self.Area = None
        self.Area0 = None
        self.globalIds = np.array([], dtype='int')
        self.Faces = []
        self.T = None
        self.X = None
        self.lambda_s1_noise = None
        self.lambda_s2_noise = None
        self.lambda_s3_noise = None
        self.lambda_v_noise = None
        self.barrier_tri0_top = None
        self.barrier_tri0_bottom = None
        self.contractility_noise = None
        self.SubstrateLambda = None
        self.InternalLambda = None
        self.ExternalLambda = None
        self.axes_lengths = None

    def compute_distance_to_centre(self, centre_of_tissue):
        """
        Compute the distance from the centre of the cell to the centre of the tissue
        :return:
        """
        return np.linalg.norm(self.Y - centre_of_tissue)

1	import numpy as np	×
2	import pyvista as pv	×
3	import vtk	×
4	from sklearn.decomposition import PCA	×
5
6	from src.pyVertexModel.geometry import face	×
7	from src.pyVertexModel.geometry.face import get_interface	×
8	from src.pyVertexModel.util.utils import copy_non_mutable_attributes	×
9
10
11	def face_centres_to_middle_of_neighbours_vertices(Geo, c_cell):	×
12	"""
13	Move the face centres to the middle of the neighbours vertices.
14	:param Geo:
15	:param c_cell:
16	:return:
17	"""
18	for num_face, _ in enumerate(Geo.Cells[c_cell].Faces):	×
19	all_edges = []	×
20	for tri in Geo.Cells[c_cell].Faces[num_face].Tris:	×
21	all_edges.append(tri.Edge)	×
22
23	all_edges = np.unique(np.concatenate(all_edges))	×
24	Geo.Cells[c_cell].Faces[num_face].Centre = np.mean(	×
25	Geo.Cells[c_cell].Y[all_edges, :], axis=0)
26
27
28	def compute_2d_circularity(area, perimeter):	×
29	"""
30	Compute the 2D circularity of the cell
31	:return:
32	"""
33	if perimeter == 0:	×
34	return 0	×
35	else:
36	return 4 * np.pi * area / perimeter ** 2	×
37
38
39	def compute_y(geo, T, cellCentre, Set):	×
40	"""
41	Compute the new Y of the cell
42
43	:param geo:
44	:param T:
45	:param cellCentre:
46	:param Set:
47	:return:
48	"""
49	x = [geo.Cells[i].X for i in T]	×
50	newY = np.mean(x, axis=0)	×
51	if sum([geo.Cells[i].AliveStatus is not None for i in T]) == 1 and "Bubbles" in Set.InputGeo:	×
52	vc = newY - cellCentre	×
53	dir = vc / np.linalg.norm(vc)	×
54	offset = Set.f * dir	×
55	newY = cellCentre + offset	×
56
57	if "Bubbles" not in Set.InputGeo:	×
58	if any(i in geo.XgTop for i in T):	×
59	newY[2] /= sum(i in geo.XgTop for i in T) / 2	×
60	elif any(i in geo.XgBottom for i in T):	×
61	newY[2] /= sum(i in geo.XgBottom for i in T) / 2	×
62	return newY	×
63
64
65	class Cell:	×
66	"""
67	Class that contains the information of a cell.
68	"""
69
70	def __init__(self, mat_file=None):	×
71	"""
72	Constructor of the class
73	:param mat_file:
74	"""
75	self.opposite_cell = None	×
76	self.lambda_s1_noise = None	×
77	self.lambda_s2_noise = None	×
78	self.lambda_s3_noise = None	×
79	self.lambda_v_noise = None	×
80	self.barrier_tri0_top = None	×
81	self.barrier_tri0_bottom = None	×
82	self.contractility_noise = None	×
83	self.SubstrateLambda = None	×
84	self.InternalLambda = None	×
85	self.ExternalLambda = None	×
86	self.axes_lengths = None	×
87	self.Y = None	×
88	self.globalIds = np.array([], dtype='int')	×
89	self.Faces = []	×
90	self.Area = None	×
91	self.Area0 = None	×
92	self.Vol = None	×
93	self.Vol0 = None	×
94	self.AliveStatus = None	×
95	self.vertices_and_faces_to_remodel = np.array([], dtype='int')	×
96	# TODO: Save contractile forces (g) to output
97	self.substrate_g = None	×
98	self.lambdaB_perc = 1	×
99
100	if mat_file is None:	×
101	self.ID = None	×
102	self.X = None	×
103	self.T = None	×
104	else:
105	self.ID = mat_file[0][0][0] - 1	×
106	self.X = np.array(mat_file[1][0], dtype=np.float64)	×
107	self.T = mat_file[2] - 1	×
108
109	if len(mat_file[4]) > 0:	×
110	self.Y = np.array(mat_file[3], dtype=np.float64)	×
111	for c_face in mat_file[4][0]:	×
112	self.Faces.append(face.Face(c_face))	×
113	if len(mat_file[5]) > 0:	×
114	self.Vol = mat_file[5][0][0]	×
115	self.Vol0 = mat_file[6][0][0]	×
116	self.Area = mat_file[7][0][0]	×
117	self.Area0 = mat_file[8][0][0]	×
118	self.globalIds = np.concatenate(mat_file[9]) - 1	×
119	self.cglobalids = mat_file[10][0][0] - 1	×
120	self.AliveStatus = mat_file[11][0][0]	×
121
122	def copy(self):	×
123	"""
124	Copy the cell
125	:return:
126	"""
127	new_cell = Cell()	×
128
129	copy_non_mutable_attributes(self, 'Faces', new_cell)	×
130
131	new_cell.Faces = [f.copy() for f in self.Faces]	×
132
133	return new_cell	×
134
135	def compute_area(self, location_filter=None):	×
136	"""
137	Compute the area of the cell
138	:param location_filter:
139	:return:
140	"""
141	total_area = 0.0	×
142	for f in range(len(self.Faces)):	×
143	if location_filter is not None:	×
144	if get_interface(self.Faces[f].InterfaceType) == get_interface(location_filter):	×
145	total_area = total_area + self.Faces[f].Area	×
146	else:
147	total_area = total_area + self.Faces[f].Area	×
148
149	return total_area	×
150
UNCOV 151	def compute_volume(self):	×
152	"""
153	Compute the volume of the cell
154	:return:
155	"""
UNCOV 156	v = 0.0	×
157	for f in range(len(self.Faces)):	×
158	c_face = self.Faces[f]	×
159	for t in range(len(c_face.Tris)):	×
160	y1 = self.Y[c_face.Tris[t].Edge[0], :] - self.X	×
161	y2 = self.Y[c_face.Tris[t].Edge[1], :] - self.X	×
162	y3 = c_face.Centre - self.X	×
163	ytri = np.array([y1, y2, y3])	×
164
UNCOV 165	current_v = np.linalg.det(ytri) / 6	×
166	# If the volume is negative, switch two the other option
UNCOV 167	if current_v < 0:	×
168	ytri = np.array([y2, y1, y3])	×
169	current_v = np.linalg.det(ytri) / 6	×
170
UNCOV 171	v += current_v	×
172
UNCOV 173	self.Vol = v	×
174	return v	×
175
UNCOV 176	def create_vtk(self):	×
177	"""
178	Create a vtk cell
179	:return:
180	"""
UNCOV 181	points = vtk.vtkPoints()	×
182	points.SetNumberOfPoints(len(self.Y) + len(self.Faces))	×
183	for i in range(len(self.Y)):	×
184	points.SetPoint(i, self.Y[i, 0], self.Y[i, 1], self.Y[i, 2])	×
185
UNCOV 186	cell = vtk.vtkCellArray()	×
187	# Go through all the faces and create the triangles for the VTK cell
UNCOV 188	total_tris = 0	×
189	for f in range(len(self.Faces)):	×
190	c_face = self.Faces[f]	×
191	points.SetPoint(len(self.Y) + f, c_face.Centre[0], c_face.Centre[1], c_face.Centre[2])	×
192	for t in range(len(c_face.Tris)):	×
193	cell.InsertNextCell(3)	×
194	cell.InsertCellPoint(c_face.Tris[t].Edge[0])	×
195	cell.InsertCellPoint(c_face.Tris[t].Edge[1])	×
196	cell.InsertCellPoint(len(self.Y) + f)	×
197
UNCOV 198	total_tris += len(c_face.Tris)	×
199
UNCOV 200	vpoly = vtk.vtkPolyData()	×
201	vpoly.SetPoints(points)	×
202
UNCOV 203	vpoly.SetPolys(cell)	×
204
205	# Get all the different properties of a cell
UNCOV 206	properties = self.compute_features()	×
207
208	# Go through the different properties of the dictionary and add them to the vtkFloatArray
UNCOV 209	for key, value in properties.items():	×
210	# Create a vtkFloatArray to store the properties of the cell
UNCOV 211	property_array = vtk.vtkFloatArray()	×
212	property_array.SetName(key)	×
213
214	# Add as many values as tris the cell has
UNCOV 215	for i in range(total_tris):	×
216	property_array.InsertNextValue(value)	×
217
218	# Add the property array to the cell data
UNCOV 219	vpoly.GetCellData().AddArray(property_array)	×
220
UNCOV 221	if key == 'ID':	×
222	# Default parameter
UNCOV 223	vpoly.GetCellData().SetScalars(property_array)	×
224
UNCOV 225	return vpoly	×
226
UNCOV 227	def create_pyvista_mesh(self):	×
228	"""
229	Create a PyVista mesh
230	:return:
231	"""
UNCOV 232	mesh = pv.PolyData(self.create_vtk())	×
233
UNCOV 234	return mesh	×
235
UNCOV 236	def compute_features(self, centre_wound=None):	×
237	"""
238	Compute the features of the cell and create a dictionary with them
239	:return: a dictionary with the features of the cell
240	"""
UNCOV 241	features = {'ID': self.ID,	×
242	'Area': self.compute_area(),
243	'Area_top': self.compute_area(location_filter=0),
244	'Area_bottom': self.compute_area(location_filter=2),
245	'Area_cellcell': self.compute_area(location_filter=1),
246	'Volume': self.compute_volume(),
247	'Height': self.compute_height(),
248	'Width': self.compute_width(),
249	'Length': self.compute_length(),
250	'Perimeter': self.compute_perimeter(),
251	'2D_circularity_top': compute_2d_circularity(self.compute_area(location_filter=0),
252	self.compute_perimeter(filter_location=0)),
253	'2d_circularity_bottom': compute_2d_circularity(self.compute_area(location_filter=2),
254	self.compute_perimeter(filter_location=2)),
255	'2D_aspect_ratio_top': self.compute_2d_aspect_ratio(filter_location=0),
256	'2D_aspect_ratio_bottom': self.compute_2d_aspect_ratio(filter_location=2),
257	'2D_aspect_ratio_cellcell': self.compute_2d_aspect_ratio(filter_location=1),
258	'3D_aspect_ratio_0_1': self.compute_3d_aspect_ratio(),
259	'3D_aspect_ratio_0_2': self.compute_3d_aspect_ratio(axis=(0, 2)),
260	'3D_aspect_ratio_1_2': self.compute_3d_aspect_ratio(axis=(1, 2)),
261	'Sphericity': self.compute_sphericity(),
262	'Elongation': self.compute_elongation(),
263	'Ellipticity': self.compute_ellipticity(),
264	'Neighbours': len(self.compute_neighbours()),
265	'Neighbours_top': len(self.compute_neighbours(location_filter=0)),
266	'Neighbours_bottom': len(self.compute_neighbours(location_filter=2)),
267	'Tilting': self.compute_tilting(),
268	'Perimeter_top': self.compute_perimeter(filter_location=0),
269	'Perimeter_bottom': self.compute_perimeter(filter_location=2),
270	'Perimeter_cellcell': self.compute_perimeter(filter_location=1),
271	}
272
UNCOV 273	if centre_wound is not None:	×
274	features['Distance_to_wound'] = self.compute_distance_to_wound(centre_wound)	×
275
UNCOV 276	return features	×
277
UNCOV 278	def create_vtk_edges(self):	×
279	"""
280	Create a vtk with only the information on the edges of the cell
281	:return:
282	"""
UNCOV 283	points = vtk.vtkPoints()	×
284	points.SetNumberOfPoints(len(self.Y))	×
285	for i in range(len(self.Y)):	×
286	points.SetPoint(i, self.Y[i, 0], self.Y[i, 1], self.Y[i, 2])	×
287
UNCOV 288	cell = vtk.vtkCellArray()	×
289	# Go through all the faces and create the triangles for the VTK cell
UNCOV 290	total_edges = 0	×
291	for f in range(len(self.Faces)):	×
292	c_face = self.Faces[f]	×
293	for t in range(len(c_face.Tris)):	×
294	cell.InsertNextCell(2)	×
295	cell.InsertCellPoint(c_face.Tris[t].Edge[0])	×
296	cell.InsertCellPoint(c_face.Tris[t].Edge[1])	×
297
UNCOV 298	total_edges += len(c_face.Tris)	×
299
UNCOV 300	vpoly = vtk.vtkPolyData()	×
301	vpoly.SetPoints(points)	×
302
UNCOV 303	vpoly.SetLines(cell)	×
304
305	# Get all the different properties of a cell
UNCOV 306	properties = self.compute_edge_features()	×
307
308	# Go through the different properties of the dictionary and add them to the vtkFloatArray
UNCOV 309	for key, value in properties[0].items():	×
310	# Create a vtkFloatArray to store the properties of the cell
UNCOV 311	property_array = vtk.vtkFloatArray()	×
312	property_array.SetName(key)	×
313
314	# Add as many values as tris the cell has
UNCOV 315	for i in range(total_edges):	×
316	if properties[i][key] is None:	×
317	property_array.InsertNextValue(0)	×
318	else:
UNCOV 319	property_array.InsertNextValue(properties[i][key])	×
320
321	# Add the property array to the cell data
UNCOV 322	vpoly.GetCellData().AddArray(property_array)	×
323
UNCOV 324	if key == 'ContractilityValue':	×
325	# Default parameter
UNCOV 326	vpoly.GetCellData().SetScalars(property_array)	×
327
UNCOV 328	return vpoly	×
329
UNCOV 330	def compute_edge_features(self):	×
331	"""
332	Compute the features of the edges of a cell and create a dictionary with them
333	:return:
334	"""
UNCOV 335	cell_features = np.array([])	×
336	for f, c_face in enumerate(self.Faces):	×
337	for t, c_tris in enumerate(c_face.Tris):	×
338	cell_features = np.append(cell_features, c_tris.compute_features())	×
339
UNCOV 340	return cell_features	×
341
UNCOV 342	def compute_height(self):	×
343	"""
344	Compute the height of the cell regardless of the orientation
345	:return:
346	"""
UNCOV 347	return np.max(self.Y[:, 2]) - np.min(self.Y[:, 2])	×
348
UNCOV 349	def compute_principal_axis_length(self):	×
350	"""
351	Compute the principal axis length of the cell
352	:return:
353	"""
354	# Perform PCA to find the principal axes
UNCOV 355	pca = PCA(n_components=3)	×
356	pca.fit(self.Y)	×
357
358	# Project points onto the principal axes
UNCOV 359	projected_points = pca.transform(self.Y)	×
360
361	# Calculate the lengths of the ellipsoid axes
UNCOV 362	min_values = projected_points.min(axis=0)	×
363	max_values = projected_points.max(axis=0)	×
364	self.axes_lengths = max_values - min_values	×
365
UNCOV 366	return self.axes_lengths	×
367
UNCOV 368	def compute_width(self):	×
369	"""
370	Compute the width of the cell
371	:return:
372	"""
UNCOV 373	return np.max(self.Y[:, 0]) - np.min(self.Y[:, 0])	×
374
UNCOV 375	def compute_length(self):	×
376	"""
377	Compute the length of the cell
378	:return:
379	"""
UNCOV 380	return np.max(self.Y[:, 1]) - np.min(self.Y[:, 1])	×
381
UNCOV 382	def compute_neighbours(self, location_filter=None):	×
383	"""
384	Compute the neighbours of the cell
385	:return:
386	"""
UNCOV 387	neighbours = []	×
388	for f in range(len(self.Faces)):	×
389	if location_filter is not None:	×
390	if get_interface(self.Faces[f].InterfaceType) == get_interface(location_filter):	×
391	for t in range(len(self.Faces[f].Tris)):	×
392	neighbours.append(self.Faces[f].Tris[t].SharedByCells)	×
393	else:
UNCOV 394	for t in range(len(self.Faces[f].Tris)):	×
395	neighbours.append(self.Faces[f].Tris[t].SharedByCells)	×
396
397	# Flatten the list of lists into a single 1-D array
UNCOV 398	if len(neighbours) > 0:	×
399	neighbours_flat = np.concatenate(neighbours)	×
400	neighbours_unique = np.unique(neighbours_flat)	×
401	neighbours_unique = neighbours_unique[neighbours_unique != self.ID]	×
402	else:
UNCOV 403	neighbours_unique = []	×
404
UNCOV 405	return neighbours_unique	×
406
UNCOV 407	def compute_tilting(self):	×
408	"""
409	Compute the tilting of the cell
410	:return:
411	"""
UNCOV 412	return np.arctan(self.compute_height() / self.compute_width())	×
413
UNCOV 414	def compute_sphericity(self):	×
415	"""
416	Compute the sphericity of the cell
417	:return:
418	"""
UNCOV 419	return 36 * np.pi * self.Vol 2 / self.Area 3	×
420
UNCOV 421	def compute_perimeter(self, filter_location=None):	×
422	"""
423	Compute the perimeter of the cell
424	:return:
425	"""
UNCOV 426	perimeter = 0.0	×
427	for f in range(len(self.Faces)):	×
428	if filter_location is not None:	×
429	if get_interface(self.Faces[f].InterfaceType) == get_interface(filter_location):	×
430	perimeter = perimeter + self.Faces[f].compute_perimeter()	×
431	else:
UNCOV 432	perimeter = perimeter + self.Faces[f].compute_perimeter()	×
433
UNCOV 434	return perimeter	×
435
UNCOV 436	def compute_ellipticity(self):	×
437	"""
438	Compute the ellipticity of the cell
439	:return:
440	"""
UNCOV 441	return self.compute_width() / self.compute_length()	×
442
UNCOV 443	def compute_elongation(self):	×
444	"""
445	Compute the elongation of the cell
446	:return:
447	"""
UNCOV 448	return self.compute_length() / self.compute_height()	×
449
UNCOV 450	def compute_3d_aspect_ratio(self, axis=(0, 1)):	×
451	"""
452	Compute the 3D aspect ratio of the cell
453	:return:
454	"""
UNCOV 455	return self.compute_principal_axis_length()[axis[0]] / self.compute_principal_axis_length()[axis[1]]	×
456
UNCOV 457	def compute_2d_aspect_ratio(self, filter_location=None):	×
458	"""
459	Compute the 2D aspect ratio of the cell
460	:return:
461	"""
UNCOV 462	perimeter = self.compute_perimeter(filter_location)	×
463	if perimeter == 0:	×
464	return 0	×
465	else:
UNCOV 466	return self.compute_area(filter_location) / perimeter ** 2	×
467
UNCOV 468	def build_y_from_x(self, geo, c_set):	×
469	Tets = self.T	×
470	dim = self.X.shape[0]	×
471	Y = np.zeros((len(Tets), dim))	×
472	for i in range(len(Tets)):	×
473	Y[i] = compute_y(geo, Tets[i], self.X, c_set)	×
474	return Y	×
475
UNCOV 476	def compute_distance_to_wound(self, centre_wound):	×
477	"""
478	Compute the distance from the centre of the cell to the centre of the wound
479	:return:
480	"""
UNCOV 481	return np.linalg.norm(self.Y - centre_wound)	×
482
UNCOV 483	def kill_cell(self):	×
484	"""
485	Kill the cell
486	:return:
487	"""
UNCOV 488	self.AliveStatus = None	×
489	self.Y = None	×
490	self.Vol = None	×
491	self.Vol0 = None	×
492	self.Area = None	×
493	self.Area0 = None	×
494	self.globalIds = np.array([], dtype='int')	×
495	self.Faces = []	×
496	self.T = None	×
497	self.X = None	×
498	self.lambda_s1_noise = None	×
499	self.lambda_s2_noise = None	×
500	self.lambda_s3_noise = None	×
501	self.lambda_v_noise = None	×
502	self.barrier_tri0_top = None	×
503	self.barrier_tri0_bottom = None	×
504	self.contractility_noise = None	×
505	self.SubstrateLambda = None	×
506	self.InternalLambda = None	×
507	self.ExternalLambda = None	×
508	self.axes_lengths = None	×
509
UNCOV 510	def compute_distance_to_centre(self, centre_of_tissue):	×
511	"""
512	Compute the distance from the centre of the cell to the centre of the tissue
513	:return:
514	"""
UNCOV 515	return np.linalg.norm(self.Y - centre_of_tissue)	×

Pablo1990 / pyVertexModel / 11368168500

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