# HG changeset patch # User mb2013 # Date 1400570939 14400 # Node ID 60ed96f5706e6401480fe103ea342d5999e9b0aa # Parent 70fc3c02ca14444e9808f4d7ff1e81158e9176d0 Uploaded diff -r 70fc3c02ca14 -r 60ed96f5706e ReMarker.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ReMarker.py Tue May 20 03:28:59 2014 -0400 @@ -0,0 +1,306 @@ +#Program for correcting landmarks. +#The object will be rotated, based on 3 symmetry points. +#The differences between the mirror points will be calculated and +# the coordinates will be corrected where needed. + +#MB + +import math +from math import * +import sys +import numpy +import os + +# Main function +def main(): + file_name = sys.argv[1] # dta file + file_connect = sys.argv[2] # connect file + side = sys.argv[3] # correct side of object + + file_outputname5 = sys.argv[5] # output dta file + output5 = open(file_outputname5, 'w') + + # to function 'open connect file' + sym1,sym2,sym3, listmirror, listspec = open_connect_file(file_connect) + + # to function 'open dta file' + listdata = open_dta_file(file_name, output5) + + # to function 'shift to zero' + shift_to_zero(listdata,sym1,sym2,sym3,listmirror,listspec, side, output5) + +# Function for extracting values of connect file +def open_connect_file(file_connect): + mirrorpoints = open(file_connect) + lines= mirrorpoints.readlines() + listmirror = [] #with mirror points + listsym = [] #with symmetry points + listspec = [] #with spec points + + #extracting symmetry landmarks and mirror landmarks of connect file + for x in range(0, len(lines)): + line = lines[x].strip() + if line == 'sym': # extract points on symmetry plane + for a in range(1,4): + listsym.append(lines[a].strip().split()) + if line == 'mirror': # extract mirror points + for y in range(x + 1,len(lines)): + try: + if lines[y].strip() != 'spec': + if lines[y].strip() == '': + break + else: + listmirror.append(lines[y].strip().split()) # add mirror points + + else: + break + + except: + break + if line == 'spec':# if specs are present + for b in range(x+1, len(lines)): + listspec.append(lines[b].strip().split()) # add points in specs + #Symmetry points + sym1 = int(listsym[0][0]) -1 + sym2 = int(listsym[1][0]) -1 + sym3 = int(listsym[2][0]) -1 + + mirrorpoints.close() + return sym1, sym2, sym3, listmirror, listspec + +# Function open .dta file +def open_dta_file(file_name, output5): + datapoints = open(file_name)#open file + datalines = datapoints.readlines() + listdata = [] # landmarklist + listdta = [] # coordinates + sublist = [] # sublist + for a in range(0,len(datalines)): + aline = datalines[a].strip() + listdta.append(datalines[a].strip().split(' ')) # notice space + + #numbers in dta file to list + if len(listdta[a]) == 3: + sublist.append(float(listdta[a][0])) + sublist.append(float(listdta[a][1])) + sublist.append(float(listdta[a][2])) + listdata.append(sublist) + sublist = [] + + #write header of dta file to outputfile + else: + output5.write("%s\n"%(aline)) + + datapoints.close() + return listdata + +# Function to calculate the values in the Z rotation matrix +def Rz_matrix(z_angle): # Rz rotation matrix + 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]] + +# Function to calculate the new coordinates rotated with Z rotation matrix +def Z_rotation(point2, z_angle): # multiplication rotation matrix and coordinates + r_z = Rz_matrix(z_angle) + rotated_z = [] + for i in range(3): + rotated_z.append((sum([r_z[i][j] * point2[j] for j in range(3)]))) + return rotated_z + +# Function to calculate the values in the X rotation matrix +def Rx_matrix(x_angle): #rotation matrix x-axis + 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))]] + +# Function to calculate the new coordinates rotated with X rotation matrix +def X_rotation(point3, x_angle): #multiplication rotation matrix and coordinates + r_x = Rx_matrix(x_angle) + rotated_x = [] + for i in range(3): + rotated_x.append((sum([r_x[i][j] * point3[j] for j in range(3)]))) + return rotated_x + +# Function to calculate the values in the Y rotation matrix +def Ry_matrix(y_angle): # Ry rotation matrix + 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))]] + +# Function to calculate the new coordinates rotated with Y rotation matrix +def Y_rotation(point4, y_angle): #multiplication rotation matrix and coordinates + r_y = Ry_matrix(y_angle) + rotated_y = [] + for i in range(3): + rotated_y.append((sum([r_y[i][j] * point4[j] for j in range(3)]))) + return rotated_y + +# Function shift object to zeropoint +def shift_to_zero(listdata, sym1,sym2,sym3, listmirror, listspec, side,output5): + # shifting object that sym3 point will be the zeropoint # + zeropoint = int(sym3) + co_zeropoint = listdata[zeropoint] #coordinates of sym3 + + listdata2 = [] + listdata3 = [] + + # minus x,y,z of point shift to zeropoint of every point + for x in range(0,len(listdata)): + listdata2.append(listdata[x][0] - co_zeropoint[0]) + listdata2.append(listdata[x][1] - co_zeropoint[1]) + listdata2.append(listdata[x][2] - co_zeropoint[2]) + listdata3.append(listdata2) + listdata2 = [] + rotate_z_axis(listdata3,sym1,sym2,sym3, listmirror, listspec, side, output5) + + +# Function for rotating the object around z axis +def rotate_z_axis(listdata3, sym1,sym2,sym3, listmirror, listspec, side, output5): + # If object is upside down, 180 degrees rotation around the Z axis is needed.# + listdatatemp = [] + + if listdata3[sym1][1] < listdata3[sym3][1]: + angle = 180 + for coordinates in range(0,len(listdata3)): + listdatatemp.append(Z_rotation(listdata3[coordinates], angle)) + listdata3 = listdatatemp + + + # calculate angle rotation z + len_z_a = listdata3[int(sym1)][0] - listdata3[int(sym3)][0] + len_z_b = listdata3[int(sym1)][1] - listdata3[int(sym3)][1] + z_angle = (math.degrees(math.atan(len_z_a/len_z_b))) + + # calculate new coordinates with rotation matrix of Z + listdata4 = [] + for coordinates in range(0, len(listdata3)): + listdata4.append(Z_rotation(listdata3[coordinates], z_angle)) + + rotate_x_axis(listdata4,sym1,sym2,sym3, listmirror, listspec, side, output5) + +# Function for rotating the object around x axis +def rotate_x_axis(listdata4,sym1,sym2,sym3, listmirror, listspec, side, output5): + #calculate angle rotation x + len_x_a = listdata4[int(sym1)][2] - listdata4[int(sym3)][0] + len_x_b = listdata4[int(sym1)][1] - listdata4[int(sym3)][0] + x_angle = -(math.degrees(math.atan(len_x_a/len_x_b))) + + # calculate new coordinates with rotation matrix of X + listdata5 = [] + for d in range(0, len(listdata4)): + listdata5.append(X_rotation(listdata4[d], x_angle)) + + rotate_y_axis(listdata5,sym1,sym2,sym3, listmirror, listspec, side, output5) + +# Function for rotating the object around y axis +def rotate_y_axis(listdata5,sym1,sym2,sym3,listmirror, listspec, side, output5): + #calculate angle rotation y + len_y_a = (listdata5[int(sym2)][0] - listdata5[int(sym3)][0]) + len_y_b = (listdata5[int(sym2)][2] - listdata5[int(sym3)][2]) + y_angle = -(math.degrees(math.atan(len_y_a/len_y_b))) + + # calculate new coordinates with rotation matrix of Y + listdata6 = [] + for d in range(0, len(listdata5)): + listdata6.append(Y_rotation(listdata5[d], y_angle)) + + + #Rotate 180 degrees around y axis when object is backwards.# + listdatatemp = [] + if listdata6[sym2][0] < listdata6[(int(listmirror[0][0]))][0]: #point sym2_x < point mirror 1_x + angle = 180 + for coordinates in range(0,len(listdata6)): + listdatatemp.append(Y_rotation(listdata6[coordinates], angle)) + listdata6 = listdatatemp + write_rotate_co(listdata6) + + # to correcting the landmarks + if len(listspec) == 0: + correct_landmarks(listdata6, listmirror, side, output5) + else: + correct_specs(listdata6, listmirror, listspec, output5) + +# Function writing rotated coordinates to file +def write_rotate_co(listdata6): + file_outputname4 = sys.argv[4] + output4 = open(file_outputname4, 'w') + #writing new coordinates to outputfile + for x in range(0,len(listdata6)): + output4.write("%.7f\t%.7f\t%.7f\n"%(listdata6[x][0], listdata6[x][1], listdata6[x][2])) + +#Function for correcting the landmarks +def correct_landmarks(listdata6, listmirror, side,output5): + + percentage_distance = [] + #calculate marge of distance + for x in range(0,len(listmirror)): + left = listmirror[x][0] #left landmark + right = listmirror[x][1] #right landmark + co_left = listdata6[int(left)-1] #left landmark coordinates + co_left_x = co_left[0] # left landmark coordinate x + co_right = listdata6[int(right)-1]#right landmark coordinates + co_right_x = co_right[0] #right landmark coordinate x + distance_x = float(co_left[0]) + float(co_right[0]) + + # if left is correct + if int(side) == 0: + percentage_distance.append(abs(distance_x) / abs(co_left_x)) #percentage + + # if right is correct + elif int(side) == 1: + percentage_distance.append(abs(distance_x) / abs(co_right_x)) #percentage + else: + print 'something wrong' + + + mean_percentage = numpy.mean(percentage_distance) #mean of percentages + std_percentage = numpy.std(percentage_distance) + #range of correct values + left_range = mean_percentage - std_percentage #left range mean minus one standard deviation + right_range = mean_percentage + std_percentage # rigth range mean plus one standard deviation + + #correcting the landmarks coordinates# + for x in range(0,len(listmirror)): + left = listmirror[x][0] + right = listmirror[x][1] + co_left = listdata6[int(left)-1] + co_right = listdata6[int(right)-1] + #if there is to much deviation take correct coordinate and project it to the other side. + if (percentage_distance[x] > right_range) or (percentage_distance[x] < left_range): + + #if left side is correct + if int(side) == 0: + listdata6[int(right)-1][0] = float(co_left[0]) * -1 + listdata6[int(right)-1][1] = float(co_left[1]) + listdata6[int(right)-1][2] = float(co_left[2]) + #if right side is correct + if int(side) == 1: + listdata6[int(left)-1][0] = float(co_right[0]) * -1 + listdata6[int(left)-1][1] = float(co_right[1]) + listdata6[int(left)-1][2] = float(co_right[2]) + write_output(listdata6,output5) + +#Function for correcting landmarks, defined in specs +def correct_specs(listdata6, listmirror, listspec,output5): + for x in range(0,len(listspec)): + number = listspec[x][0] + for spec in range(0, len(listmirror)): + try: + pos = listmirror[spec].index(number) #find number in mirror list + #extracting the opposite landmark number + if pos == 0: + number2 = listmirror[spec][1] + elif pos == 1 : + number2 = listmirror[spec][0] + else: + print 'wrong' + #replace the wrong coordinates + listdata6[int(number)-1][0] = float(listdata6[int(number2)-1][0])*-1 + listdata6[int(number)-1][1] = float(listdata6[int(number2)-1][1]) + listdata6[int(number)-1][2] = float(listdata6[int(number2)-1][2]) + except: + continue + write_output(listdata6,output5) + +#Function for writing the corrected landmarks to outputfile +def write_output(listdata6,output5): + + for x in range(0,len(listdata6)): + output5.write("%.7f %.7f %.7f\n"%(listdata6[x][0], listdata6[x][1], listdata6[x][2])) + output5.close() +main()