-
Notifications
You must be signed in to change notification settings - Fork 1
/
PrincipleComponentAnalysis.py
120 lines (102 loc) · 3.99 KB
/
PrincipleComponentAnalysis.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
import vtk
from sklearn.decomposition import PCA
import numpy as np
class PRINCIPLE_COMPONENT_ANALYSIS():
def __init__(self,volume):
self.volume=volume
#Get the outer surface of the volume
self.surface=vtk.vtkDataSetSurfaceFilter()
self.surface.SetInputData(volume)
self.surface.Update()
self.volume=self.surface
def main(self):
volume_data =self.volume.GetOutput()
surface_data=self.surface.GetOutput()
#Get the point cloud data
Centroid,Coords=self.get_point_data(volume_data)
#Perform PCA analysis
Norm1,Norm2=self.major_axis(Coords)
Norm1=Norm1*-1
#Get the Apex Point of the left ventricle
Apex=self.find_apex(Coords,Centroid,Norm1,surface_data)
#Get the size of the Myocardium (Apex to valve)
#Size1=np.max(np.power(np.power(Coords[:,0]-Apex[0],2)+np.power(Coords[:,1]-Apex[1],2)+np.power(Coords[:,2]-Apex[2],2),0.5))
Size=abs(self.get_size(Coords,Apex,Norm1,Norm2))
#print ("\n")
#print ("--- LV Morpholocial Analysis")
print ("----------The Centroid of LV: ",Centroid)
print ("----------The Major Axis Norm: ",Norm1)
print ("----------The Minor Axis Norm: ",Norm2)
print ("----------The Location of Apex: ",Apex)
print ("----------The Size of Myocardium:",Size)
del self.surface,self.volume,volume_data,surface_data,
return np.array(Centroid),Norm1,Norm2,np.array(Apex),Size
def get_size(self,Coords,Apex,Norm1,Norm2):
#Get the projection of the max distance onto the myocardial centerline
distance=np.power(np.power(Coords[:,0]-Apex[0],2)+np.power(Coords[:,1]-Apex[1],2)+np.power(Coords[:,2]-Apex[2],2),0.5)
max_id=np.argmax(distance)
coord_max=Coords[max_id,:]
vector_max=[coord_max[0]-Apex[0],coord_max[1]-Apex[1],coord_max[2]-Apex[2]]
vector_max_mag=np.linalg.norm(vector_max)
vector_max_norm=vector_max/vector_max_mag
#Now get the theta
cos_theta=np.dot(vector_max_norm,Norm1)
Size=vector_max_mag*cos_theta
return Size
def find_apex(self,Coords,Cent,Norm1,surface_data):
#Get the min and max and scale the norm
xmin=min(Coords[:,0]);xmax=max(Coords[:,0])
ymin=min(Coords[:,1]);ymax=max(Coords[:,1])
zmin=min(Coords[:,2]);zmax=max(Coords[:,2])
Length=max(xmax-xmin,ymax-ymin,zmax-zmin)
Norm1s=Norm1*Length
#Get the source point
pSource=[Cent[0]-Norm1s[0],Cent[1]-Norm1s[1],Cent[2]-Norm1s[2]]
pTarget=[Cent[0]+Norm1s[0],Cent[1]+Norm1s[1],Cent[2]+Norm1s[2]]
#Create a bounding box of the surface mesh
obbTree = vtk.vtkOBBTree()
obbTree.SetDataSet(surface_data)
obbTree.BuildLocator()
#Create a point array to store intersection points
pointsVTKintersection = vtk.vtkPoints()
code=obbTree.IntersectWithLine(pSource, pTarget, pointsVTKintersection, None)
pointsVTKIntersectionData = pointsVTKintersection.GetData()
noPointsVTKIntersection = pointsVTKIntersectionData.GetNumberOfTuples()
pointsIntersection = []
distanceCentroid=[]
for idx in range(noPointsVTKIntersection):
tup = pointsVTKIntersectionData.GetTuple3(idx)
dx=(Cent[0]-tup[0])**2
dy=(Cent[1]-tup[1])**2
dz=(Cent[2]-tup[2])**2
distanceCentroid.append((dx+dy+dz)**0.5)
pointsIntersection.append(tup)
#Now find the furthest point from Centroid as Apex
max_idx=np.argmax(distanceCentroid)
apex_coord=pointsIntersection[max_idx]
del obbTree,pointsVTKintersection,code,pointsVTKIntersectionData,noPointsVTKIntersection
return apex_coord
def major_axis(self,Coords):
pca=PCA(n_components=2)
pca.fit(Coords)
Norm1=pca.components_[0]
Norm2=pca.components_[1]
return Norm1,Norm2
def get_point_data(self,surface_data):
#Number of points
N_points=surface_data.GetNumberOfPoints()
#Get Array Points
Coords=np.zeros(shape=(N_points,3))
Centroid=[0,0,0]
for i in range(N_points):
coord_=surface_data.GetPoint(i)
Coords[i,0]=coord_[0]
Coords[i,1]=coord_[1]
Coords[i,2]=coord_[2]
Centroid[0]+=coord_[0]
Centroid[1]+=coord_[1]
Centroid[2]+=coord_[2]
Centroid[0]=Centroid[0]/N_points
Centroid[1]=Centroid[1]/N_points
Centroid[2]=Centroid[2]/N_points
return Centroid,Coords