Mercurial > repos > mb2013 > nepenthes_3dpca

comparison ReMarker.py @ 15:60ed96f5706e draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mb2013
date Tue, 20 May 2014 03:28:59 -0400
parents
children
comparison
equal deleted inserted replaced
14:70fc3c02ca14 15:60ed96f5706e
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()