|
15
|
1 #Program for correcting landmarks.
|
|
|
2 #The object will be rotated, based on 3 symmetry points.
|
|
|
3 #The differences between the mirror points will be calculated and
|
|
|
4 # the coordinates will be corrected where needed.
|
|
|
5
|
|
|
6 #MB
|
|
|
7
|
|
|
8 import math
|
|
|
9 from math import *
|
|
|
10 import sys
|
|
|
11 import numpy
|
|
|
12 import os
|
|
|
13
|
|
|
14 # Main function
|
|
|
15 def main():
|
|
|
16 file_name = sys.argv[1] # dta file
|
|
|
17 file_connect = sys.argv[2] # connect file
|
|
|
18 side = sys.argv[3] # correct side of object
|
|
|
19
|
|
|
20 file_outputname5 = sys.argv[5] # output dta file
|
|
|
21 output5 = open(file_outputname5, 'w')
|
|
|
22
|
|
|
23 # to function 'open connect file'
|
|
|
24 sym1,sym2,sym3, listmirror, listspec = open_connect_file(file_connect)
|
|
|
25
|
|
|
26 # to function 'open dta file'
|
|
|
27 listdata = open_dta_file(file_name, output5)
|
|
|
28
|
|
|
29 # to function 'shift to zero'
|
|
|
30 shift_to_zero(listdata,sym1,sym2,sym3,listmirror,listspec, side, output5)
|
|
|
31
|
|
|
32 # Function for extracting values of connect file
|
|
|
33 def open_connect_file(file_connect):
|
|
|
34 mirrorpoints = open(file_connect)
|
|
|
35 lines= mirrorpoints.readlines()
|
|
|
36 listmirror = [] #with mirror points
|
|
|
37 listsym = [] #with symmetry points
|
|
|
38 listspec = [] #with spec points
|
|
|
39
|
|
|
40 #extracting symmetry landmarks and mirror landmarks of connect file
|
|
|
41 for x in range(0, len(lines)):
|
|
|
42 line = lines[x].strip()
|
|
|
43 if line == 'sym': # extract points on symmetry plane
|
|
|
44 for a in range(1,4):
|
|
|
45 listsym.append(lines[a].strip().split())
|
|
|
46 if line == 'mirror': # extract mirror points
|
|
|
47 for y in range(x + 1,len(lines)):
|
|
|
48 try:
|
|
|
49 if lines[y].strip() != 'spec':
|
|
|
50 if lines[y].strip() == '':
|
|
|
51 break
|
|
|
52 else:
|
|
|
53 listmirror.append(lines[y].strip().split()) # add mirror points
|
|
|
54
|
|
|
55 else:
|
|
|
56 break
|
|
|
57
|
|
|
58 except:
|
|
|
59 break
|
|
|
60 if line == 'spec':# if specs are present
|
|
|
61 for b in range(x+1, len(lines)):
|
|
|
62 listspec.append(lines[b].strip().split()) # add points in specs
|
|
|
63 #Symmetry points
|
|
|
64 sym1 = int(listsym[0][0]) -1
|
|
|
65 sym2 = int(listsym[1][0]) -1
|
|
|
66 sym3 = int(listsym[2][0]) -1
|
|
|
67
|
|
|
68 mirrorpoints.close()
|
|
|
69 return sym1, sym2, sym3, listmirror, listspec
|
|
|
70
|
|
|
71 # Function open .dta file
|
|
|
72 def open_dta_file(file_name, output5):
|
|
|
73 datapoints = open(file_name)#open file
|
|
|
74 datalines = datapoints.readlines()
|
|
|
75 listdata = [] # landmarklist
|
|
|
76 listdta = [] # coordinates
|
|
|
77 sublist = [] # sublist
|
|
|
78 for a in range(0,len(datalines)):
|
|
|
79 aline = datalines[a].strip()
|
|
|
80 listdta.append(datalines[a].strip().split(' ')) # notice space
|
|
|
81
|
|
|
82 #numbers in dta file to list
|
|
|
83 if len(listdta[a]) == 3:
|
|
|
84 sublist.append(float(listdta[a][0]))
|
|
|
85 sublist.append(float(listdta[a][1]))
|
|
|
86 sublist.append(float(listdta[a][2]))
|
|
|
87 listdata.append(sublist)
|
|
|
88 sublist = []
|
|
|
89
|
|
|
90 #write header of dta file to outputfile
|
|
|
91 else:
|
|
|
92 output5.write("%s\n"%(aline))
|
|
|
93
|
|
|
94 datapoints.close()
|
|
|
95 return listdata
|
|
|
96
|
|
|
97 # Function to calculate the values in the Z rotation matrix
|
|
|
98 def Rz_matrix(z_angle): # Rz rotation matrix
|
|
|
99 return [[cos(math.radians(z_angle)), -sin(math.radians(z_angle)), 0.0],[sin(math.radians(z_angle)),cos(math.radians(z_angle)),0.0],[0.0, 0.0, 1.0]]
|
|
|
100
|
|
|
101 # Function to calculate the new coordinates rotated with Z rotation matrix
|
|
|
102 def Z_rotation(point2, z_angle): # multiplication rotation matrix and coordinates
|
|
|
103 r_z = Rz_matrix(z_angle)
|
|
|
104 rotated_z = []
|
|
|
105 for i in range(3):
|
|
|
106 rotated_z.append((sum([r_z[i][j] * point2[j] for j in range(3)])))
|
|
|
107 return rotated_z
|
|
|
108
|
|
|
109 # Function to calculate the values in the X rotation matrix
|
|
|
110 def Rx_matrix(x_angle): #rotation matrix x-axis
|
|
|
111 return [[1, 0, 0],[0,cos(math.radians(x_angle)),-sin(math.radians(x_angle))],[0,sin(math.radians(x_angle)),cos(math.radians(x_angle))]]
|
|
|
112
|
|
|
113 # Function to calculate the new coordinates rotated with X rotation matrix
|
|
|
114 def X_rotation(point3, x_angle): #multiplication rotation matrix and coordinates
|
|
|
115 r_x = Rx_matrix(x_angle)
|
|
|
116 rotated_x = []
|
|
|
117 for i in range(3):
|
|
|
118 rotated_x.append((sum([r_x[i][j] * point3[j] for j in range(3)])))
|
|
|
119 return rotated_x
|
|
|
120
|
|
|
121 # Function to calculate the values in the Y rotation matrix
|
|
|
122 def Ry_matrix(y_angle): # Ry rotation matrix
|
|
|
123 return [[cos(math.radians(y_angle)), 0.0, sin(math.radians(y_angle))],[0.0, 1.0, 0.0],[-sin(math.radians(y_angle)),0.0, cos(math.radians(y_angle))]]
|
|
|
124
|
|
|
125 # Function to calculate the new coordinates rotated with Y rotation matrix
|
|
|
126 def Y_rotation(point4, y_angle): #multiplication rotation matrix and coordinates
|
|
|
127 r_y = Ry_matrix(y_angle)
|
|
|
128 rotated_y = []
|
|
|
129 for i in range(3):
|
|
|
130 rotated_y.append((sum([r_y[i][j] * point4[j] for j in range(3)])))
|
|
|
131 return rotated_y
|
|
|
132
|
|
|
133 # Function shift object to zeropoint
|
|
|
134 def shift_to_zero(listdata, sym1,sym2,sym3, listmirror, listspec, side,output5):
|
|
|
135 # shifting object that sym3 point will be the zeropoint #
|
|
|
136 zeropoint = int(sym3)
|
|
|
137 co_zeropoint = listdata[zeropoint] #coordinates of sym3
|
|
|
138
|
|
|
139 listdata2 = []
|
|
|
140 listdata3 = []
|
|
|
141
|
|
|
142 # minus x,y,z of point shift to zeropoint of every point
|
|
|
143 for x in range(0,len(listdata)):
|
|
|
144 listdata2.append(listdata[x][0] - co_zeropoint[0])
|
|
|
145 listdata2.append(listdata[x][1] - co_zeropoint[1])
|
|
|
146 listdata2.append(listdata[x][2] - co_zeropoint[2])
|
|
|
147 listdata3.append(listdata2)
|
|
|
148 listdata2 = []
|
|
|
149 rotate_z_axis(listdata3,sym1,sym2,sym3, listmirror, listspec, side, output5)
|
|
|
150
|
|
|
151
|
|
|
152 # Function for rotating the object around z axis
|
|
|
153 def rotate_z_axis(listdata3, sym1,sym2,sym3, listmirror, listspec, side, output5):
|
|
|
154 # If object is upside down, 180 degrees rotation around the Z axis is needed.#
|
|
|
155 listdatatemp = []
|
|
|
156
|
|
|
157 if listdata3[sym1][1] < listdata3[sym3][1]:
|
|
|
158 angle = 180
|
|
|
159 for coordinates in range(0,len(listdata3)):
|
|
|
160 listdatatemp.append(Z_rotation(listdata3[coordinates], angle))
|
|
|
161 listdata3 = listdatatemp
|
|
|
162
|
|
|
163
|
|
|
164 # calculate angle rotation z
|
|
|
165 len_z_a = listdata3[int(sym1)][0] - listdata3[int(sym3)][0]
|
|
|
166 len_z_b = listdata3[int(sym1)][1] - listdata3[int(sym3)][1]
|
|
|
167 z_angle = (math.degrees(math.atan(len_z_a/len_z_b)))
|
|
|
168
|
|
|
169 # calculate new coordinates with rotation matrix of Z
|
|
|
170 listdata4 = []
|
|
|
171 for coordinates in range(0, len(listdata3)):
|
|
|
172 listdata4.append(Z_rotation(listdata3[coordinates], z_angle))
|
|
|
173
|
|
|
174 rotate_x_axis(listdata4,sym1,sym2,sym3, listmirror, listspec, side, output5)
|
|
|
175
|
|
|
176 # Function for rotating the object around x axis
|
|
|
177 def rotate_x_axis(listdata4,sym1,sym2,sym3, listmirror, listspec, side, output5):
|
|
|
178 #calculate angle rotation x
|
|
|
179 len_x_a = listdata4[int(sym1)][2] - listdata4[int(sym3)][0]
|
|
|
180 len_x_b = listdata4[int(sym1)][1] - listdata4[int(sym3)][0]
|
|
|
181 x_angle = -(math.degrees(math.atan(len_x_a/len_x_b)))
|
|
|
182
|
|
|
183 # calculate new coordinates with rotation matrix of X
|
|
|
184 listdata5 = []
|
|
|
185 for d in range(0, len(listdata4)):
|
|
|
186 listdata5.append(X_rotation(listdata4[d], x_angle))
|
|
|
187
|
|
|
188 rotate_y_axis(listdata5,sym1,sym2,sym3, listmirror, listspec, side, output5)
|
|
|
189
|
|
|
190 # Function for rotating the object around y axis
|
|
|
191 def rotate_y_axis(listdata5,sym1,sym2,sym3,listmirror, listspec, side, output5):
|
|
|
192 #calculate angle rotation y
|
|
|
193 len_y_a = (listdata5[int(sym2)][0] - listdata5[int(sym3)][0])
|
|
|
194 len_y_b = (listdata5[int(sym2)][2] - listdata5[int(sym3)][2])
|
|
|
195 y_angle = -(math.degrees(math.atan(len_y_a/len_y_b)))
|
|
|
196
|
|
|
197 # calculate new coordinates with rotation matrix of Y
|
|
|
198 listdata6 = []
|
|
|
199 for d in range(0, len(listdata5)):
|
|
|
200 listdata6.append(Y_rotation(listdata5[d], y_angle))
|
|
|
201
|
|
|
202
|
|
|
203 #Rotate 180 degrees around y axis when object is backwards.#
|
|
|
204 listdatatemp = []
|
|
|
205 if listdata6[sym2][0] < listdata6[(int(listmirror[0][0]))][0]: #point sym2_x < point mirror 1_x
|
|
|
206 angle = 180
|
|
|
207 for coordinates in range(0,len(listdata6)):
|
|
|
208 listdatatemp.append(Y_rotation(listdata6[coordinates], angle))
|
|
|
209 listdata6 = listdatatemp
|
|
|
210 write_rotate_co(listdata6)
|
|
|
211
|
|
|
212 # to correcting the landmarks
|
|
|
213 if len(listspec) == 0:
|
|
|
214 correct_landmarks(listdata6, listmirror, side, output5)
|
|
|
215 else:
|
|
|
216 correct_specs(listdata6, listmirror, listspec, output5)
|
|
|
217
|
|
|
218 # Function writing rotated coordinates to file
|
|
|
219 def write_rotate_co(listdata6):
|
|
|
220 file_outputname4 = sys.argv[4]
|
|
|
221 output4 = open(file_outputname4, 'w')
|
|
|
222 #writing new coordinates to outputfile
|
|
|
223 for x in range(0,len(listdata6)):
|
|
|
224 output4.write("%.7f\t%.7f\t%.7f\n"%(listdata6[x][0], listdata6[x][1], listdata6[x][2]))
|
|
|
225
|
|
|
226 #Function for correcting the landmarks
|
|
|
227 def correct_landmarks(listdata6, listmirror, side,output5):
|
|
|
228
|
|
|
229 percentage_distance = []
|
|
|
230 #calculate marge of distance
|
|
|
231 for x in range(0,len(listmirror)):
|
|
|
232 left = listmirror[x][0] #left landmark
|
|
|
233 right = listmirror[x][1] #right landmark
|
|
|
234 co_left = listdata6[int(left)-1] #left landmark coordinates
|
|
|
235 co_left_x = co_left[0] # left landmark coordinate x
|
|
|
236 co_right = listdata6[int(right)-1]#right landmark coordinates
|
|
|
237 co_right_x = co_right[0] #right landmark coordinate x
|
|
|
238 distance_x = float(co_left[0]) + float(co_right[0])
|
|
|
239
|
|
|
240 # if left is correct
|
|
|
241 if int(side) == 0:
|
|
|
242 percentage_distance.append(abs(distance_x) / abs(co_left_x)) #percentage
|
|
|
243
|
|
|
244 # if right is correct
|
|
|
245 elif int(side) == 1:
|
|
|
246 percentage_distance.append(abs(distance_x) / abs(co_right_x)) #percentage
|
|
|
247 else:
|
|
|
248 print 'something wrong'
|
|
|
249
|
|
|
250
|
|
|
251 mean_percentage = numpy.mean(percentage_distance) #mean of percentages
|
|
|
252 std_percentage = numpy.std(percentage_distance)
|
|
|
253 #range of correct values
|
|
|
254 left_range = mean_percentage - std_percentage #left range mean minus one standard deviation
|
|
|
255 right_range = mean_percentage + std_percentage # rigth range mean plus one standard deviation
|
|
|
256
|
|
|
257 #correcting the landmarks coordinates#
|
|
|
258 for x in range(0,len(listmirror)):
|
|
|
259 left = listmirror[x][0]
|
|
|
260 right = listmirror[x][1]
|
|
|
261 co_left = listdata6[int(left)-1]
|
|
|
262 co_right = listdata6[int(right)-1]
|
|
|
263 #if there is to much deviation take correct coordinate and project it to the other side.
|
|
|
264 if (percentage_distance[x] > right_range) or (percentage_distance[x] < left_range):
|
|
|
265
|
|
|
266 #if left side is correct
|
|
|
267 if int(side) == 0:
|
|
|
268 listdata6[int(right)-1][0] = float(co_left[0]) * -1
|
|
|
269 listdata6[int(right)-1][1] = float(co_left[1])
|
|
|
270 listdata6[int(right)-1][2] = float(co_left[2])
|
|
|
271 #if right side is correct
|
|
|
272 if int(side) == 1:
|
|
|
273 listdata6[int(left)-1][0] = float(co_right[0]) * -1
|
|
|
274 listdata6[int(left)-1][1] = float(co_right[1])
|
|
|
275 listdata6[int(left)-1][2] = float(co_right[2])
|
|
|
276 write_output(listdata6,output5)
|
|
|
277
|
|
|
278 #Function for correcting landmarks, defined in specs
|
|
|
279 def correct_specs(listdata6, listmirror, listspec,output5):
|
|
|
280 for x in range(0,len(listspec)):
|
|
|
281 number = listspec[x][0]
|
|
|
282 for spec in range(0, len(listmirror)):
|
|
|
283 try:
|
|
|
284 pos = listmirror[spec].index(number) #find number in mirror list
|
|
|
285 #extracting the opposite landmark number
|
|
|
286 if pos == 0:
|
|
|
287 number2 = listmirror[spec][1]
|
|
|
288 elif pos == 1 :
|
|
|
289 number2 = listmirror[spec][0]
|
|
|
290 else:
|
|
|
291 print 'wrong'
|
|
|
292 #replace the wrong coordinates
|
|
|
293 listdata6[int(number)-1][0] = float(listdata6[int(number2)-1][0])*-1
|
|
|
294 listdata6[int(number)-1][1] = float(listdata6[int(number2)-1][1])
|
|
|
295 listdata6[int(number)-1][2] = float(listdata6[int(number2)-1][2])
|
|
|
296 except:
|
|
|
297 continue
|
|
|
298 write_output(listdata6,output5)
|
|
|
299
|
|
|
300 #Function for writing the corrected landmarks to outputfile
|
|
|
301 def write_output(listdata6,output5):
|
|
|
302
|
|
|
303 for x in range(0,len(listdata6)):
|
|
|
304 output5.write("%.7f %.7f %.7f\n"%(listdata6[x][0], listdata6[x][1], listdata6[x][2]))
|
|
|
305 output5.close()
|
|
|
306 main()
|
