12162144775

Committed 04 Dec 2024 02:45PM UTC coverage: 0.77% (-0.003%) from 0.773%

Build # 12162144775

Build Type

push

github

Committed by

Pablo1990

Commit Message

Improving remodelling

Run Details

0 of 1050 branches covered (0.0%)

Branch coverage included in aggregate %.

0 of 53 new or added lines in 3 files covered. (0.0%)

1 existing line in 1 file now uncovered.

51 of 5577 relevant lines covered (0.91%)

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

Pablo1990 / pyVertexModel / 12162144775

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