11836362134

Committed 14 Nov 2024 11:19AM UTC coverage: 0.789% (-0.002%) from 0.791%

Build # 11836362134

Build Type

push

github

Committed by

Pablo1990

Commit Message

Displaying edges and different perspectives

Also, better colours

Run Details

0 of 1010 branches covered (0.0%)

Branch coverage included in aggregate %.

0 of 19 new or added lines in 2 files covered. (0.0%)

164 existing lines in 2 files now uncovered.

51 of 5454 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.axes_lengths = None
        self.Y = None
        self.globalIds = np.array([], dtype='int')
        self.Faces = []
        self.Area = None
        self.Vol = None
        self.AliveStatus = None
        self.vertices_and_faces_to_remodel = np.array([], dtype='int')

        ## Individual mechanical parameters
        # Surface area
        self.lambda_s1_perc = 1
        self.lambda_s2_perc = 1
        self.lambda_s3_perc = 1
        # Volume
        self.lambda_v_perc = 1
        # Aspect ratio/elongation
        self.lambda_r_perc = 1
        # Contractility
        self.c_line_tension_perc = 1
        # Substrate k
        self.k_substrate_perc = 1
        # Area Energy Barrier
        self.lambda_b_perc = 1

        ## Reference values
        # Aspect ratio
        self.barrier_tri0_top = None
        self.barrier_tri0_bottom = None
        # Area
        self.Area0 = None
        # Volume
        self.Vol0 = None



        # In case the geometry is not provided, the cell is empty
        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 create_pyvista_edges(self):
        """
        Create a PyVista mesh
        :return:
        """
        mesh = pv.PolyData(self.create_vtk_edges())

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

Pablo1990 / pyVertexModel / 11836362134

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