Mercurial > repos > mb2013 > nepenthes_3dpca
changeset 13:fc4bc4ab91b7 draft
Uploaded
author | mb2013 |
---|---|
date | Tue, 20 May 2014 03:28:13 -0400 |
parents | f69fc89bb392 |
children | 70fc3c02ca14 |
files | ReConstructor.py |
diffstat | 1 files changed, 597 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ReConstructor.py Tue May 20 03:28:13 2014 -0400 @@ -0,0 +1,597 @@ +# Reconstruction of faces of a .ply file, +# based on the symmetry plane +# MB +import math +from math import * +import sys +import numpy +from time import gmtime, strftime + +# Function main, +def main(): + var_col = 0 + name_file_ply = sys.argv[1] + output = open('new_coordinates2.ply', 'w') + file_ply = open(name_file_ply) + out_log = open(str(sys.argv[15]), 'w') + + out_log.write('Start time: %s\n\nInput user:\n'%(time())) + x_1 = sys.argv[9] + x_2 = sys.argv[10] + y_1 = sys.argv[11] + y_2 = sys.argv[12] + z_1 = sys.argv[13] + z_2 = sys.argv[14] + # to function 'user input side' + side = user_input_side(out_log) + # to function 'user input sections' + number_of_boxes_hor,number_of_boxes_ver,number_of_boxes_z = user_input_sections(out_log) + # to function 'user input factor stdev' + factor = user_input_factor_stdev(out_log) + # to function user input colors + col_r, col_g, col_b = user_input_colors(out_log) + # to function 'extracting header' + var_header, var_vertex_nm,var_face_nm = extracting_header(file_ply) + # to funtion 'array coordinates' + matrix, matrix2,matrix3 = array_coordinates(name_file_ply, var_vertex_nm, var_header, var_face_nm) + # to function 'calc min max' + amin,amax = calc_min_max(out_log, matrix, x_1,x_2,y_1,y_2,z_1,z_2) + + # to function 'range sections' + steps_x,steps_y,steps_z,x_window_min_x,x_window_max_x = range_sections(amin,amax, + number_of_boxes_hor, number_of_boxes_ver, + number_of_boxes_z) + # to function 'create list ranges' + newlist = create_list_ranges(amin,amax,number_of_boxes_hor, number_of_boxes_ver, + number_of_boxes_z, x_window_min_x, x_window_max_x,steps_x,steps_y,steps_z) + # to function 'create list coordinates' + newlist2, indexlist = create_list_coordinates(newlist) + # to function 'fill sections coordinates' + newlist, indexlist = fill_sections_coordinates(newlist, newlist2, matrix, indexlist) + # to function 'calc sections' + difference, differencemean = calc_sections(number_of_boxes_ver,number_of_boxes_hor,number_of_boxes_z,newlist) + # to function 'calc mean stdev' + mean_percentage, std_percentage = calc_mean_stdev(differencemean,factor) + # to function 'calc range' + left_range, right_range = calc_range(mean_percentage,std_percentage) + # to function 'collect mirror values' + total, listindexcount, listindex, listindextotalcount, totalcount = collect_mirror_values(out_log, side, output, + right_range, left_range, + mean_percentage, std_percentage, + newlist, indexlist, difference, + number_of_boxes_hor,number_of_boxes_ver,number_of_boxes_z, + col_r, col_g, col_b) + + # to function 'find original faces' + m_face,m_face_square = find_original_faces(matrix2,matrix3, listindex) + + # to function 'replace index values' + out_faces = replace_index_values(m_face,listindex,listindextotalcount) + out_faces_2 = [] + # if there are square faces, do function replace index values again + if len(m_face_square)!= 0: + out_faces_2 = replace_index_values(m_face_square,listindex,listindextotalcount) + output.close() + + # to function 'write output' + write_output(name_file_ply, out_log, totalcount, var_header, var_vertex_nm, var_face_nm,out_faces, out_faces_2) + #a = strftime("%a, %d %b %Y %H:%M:%S", gmtime()) + + out_log.write('End time: %s'%(time())) + out_log.close() # close log + +# Function time +def time(): + a = strftime("%a, %d %b %Y %H:%M:%S", gmtime()) + return a + +# Function user input side +def user_input_side(out_log): # user input side + side = int(sys.argv[3]) + if side == 0: + out_log.write('Left side of object is correct') + else: + out_log.write('Right side of object is correct') + return side + +# Function user input sections +def user_input_sections(out_log): + # number of boxes on x axis, this will be multiplied. + number_of_boxes_hor = int(sys.argv[4]) + # number of boxes on y axis + number_of_boxes_ver = int(sys.argv[5]) + # number of boxes on z axis + number_of_boxes_z = int(sys.argv[6]) + out_log.write("number of boxes x axis:\t%s\nnumber of boxes y axis:\t%s\nnumber of boxes z axis:\t%s\n\n"% + ((number_of_boxes_hor *2),(number_of_boxes_ver), (number_of_boxes_z))) + + return number_of_boxes_hor, number_of_boxes_ver, number_of_boxes_z + +# Function user input colors of reconstructed surface +def user_input_colors(out_log): + col = sys.argv[8] + # four colors can be chosen + if int(col) == 0: + out_log.write('color is pink\n') + col_r = 236 + col_g = 149 + col_b = 221 + elif int(col) == 1: + out_log.write('color is blue\n') + col_r = 192 + col_g = 245 + col_b = 250 + elif int(col) == 2: + out_log.write('color is green\n') + col_r = 0 + col_g = 255 + col_b = 0 + else: + out_log.write('color is red\n') + col_r = 255 + col_g = 0 + col_b = 0 + + return col_r, col_g, col_b + +# Function user input factor standard deviation +def user_input_factor_stdev(out_log): + factor = int(sys.argv[7]) + out_log.write("Factor standard deviation:\t%s\n\n"%(factor)) + + return factor + +# Function extracting values of header ply file +def extracting_header(file_ply): + for x in range(0, 20): + readheader = file_ply.readline().strip().split() + if readheader[0] == 'end_header': # when the words 'end_header' are found + var_header = x # var_header is line number + if readheader[0] == "element" and readheader[1] == "vertex": # when 'element vertex' found + var_vertex_nm = readheader[2] # amount of vertexen + if readheader[0] == "element" and readheader[1] == "face": # when 'element face' found + var_face_nm = readheader[2] # amount of faces + file_ply.close() + + return var_header, var_vertex_nm,var_face_nm + +# Function extracting coordintates of ply file +def array_coordinates(name_file_ply, var_vertex_nm, var_header,var_face_nm): + file1 = open(name_file_ply) + matrix = numpy.zeros((int(var_vertex_nm),3)) # creating empty numpy array of amount of vertexen + matrix2 = numpy.empty((int(var_face_nm),3)) # creating empty numpy array of amount of triangle faces + matrix3 = numpy.empty((int(var_face_nm),4)) # creating empty numpy array of amounf of square faces 4 + count = 0 # counter for vertexen + count2 = 0 # counter for faces + # Putting all the vertexen and faces in a matrix + for a in range(0, (int(var_header) + int(var_vertex_nm) + int(var_face_nm) + 1)): + line = file1.readline().strip().split() # reading every line + if int(var_header) < a < (int(var_vertex_nm)+ int(var_header)+1): #vertexen + matrix[count][0] = float(line[0]) # x coordinate + matrix[count][1] = float(line[1]) # y coordinate + matrix[count][2] = float(line[2]) # z coordinate + count += 1 # counter + 1 for next line in matrix + if a > (int(var_vertex_nm)+ int(var_header)): # faces + if len(line) == 4: # if triangle faces exists + matrix2[count2][0] = int(line[1]) # first coordinate for face + matrix2[count2][1] = int(line[2]) # second coordinate for face + matrix2[count2][2] = int(line[3]) # third coordinate for face + #check1 = True + + + if len(line) == 5: # if square faces exists + matrix3[count2][0] = int(line[1]) # first coordinate for face + matrix3[count2][1] = int(line[2]) # second coordinate for face + matrix3[count2][2] = int(line[3]) # third coordinate for face + matrix3[count2][3] = int(line[4]) # fourth coordinate for face + #check2 = True + count2 += 1 # counter + 1 for next line in matrix + file1.close() + + # to function calc_histogram + calc_histogram(matrix) + + return matrix,matrix2,matrix3 + +# Function histogram +def calc_histogram(matrix): # histogram of x,y, z coordinates + out_histo = open(str(sys.argv[16]), 'w') # open output + x_co = numpy.array(matrix[:,0]) + histo_x = numpy.histogram(x_co, bins = 50) # calculate x coordinate distribution + y_co = numpy.array(matrix[:,1]) + histo_y = numpy.histogram(y_co, bins = 50) # calculate y coordinate distribution + z_co = numpy.array(matrix[:,2]) + histo_z = numpy.histogram(z_co, bins = 50) # calculate z coordinate distribution + + # write histogram output + out_histo.write("x nm point\tx coordinate\ty nm point\ty coordinate\tz nm point\tz coordinate\n") + for item in range(0,len(histo_x[0])): + out_histo.write("%s\t%s\t\t%s\t%s\t\t%s\t%s\n"%(histo_x[0][item], histo_x[1][item], + histo_y[0][item], histo_y[1][item], + histo_z[0][item], histo_z[1][item], )) + + + out_histo.close() #close histogram output + +# Function calculation minimum and maximum of x,y and z +def calc_min_max(out_log, matrix, x_1,x_2,y_1,y_2,z_1,z_2): + amin = numpy.amin(matrix, axis = 0) #minima of x y and z + amax = numpy.amax(matrix, axis = 0) #maxima of x y and z + # if user input max and min values, change the min and max + if x_1 != 'standard': + amin[0] = int(x_1) + if x_2 != 'standard': + amax[0] = int(x_2) + if y_1 != 'standard': + amin[1] = int(y_1) + if y_2 != 'standard': + amax[1] = int(y_2) + if z_1 != 'standard': + amin[2] = int(z_1) + if z_2 != 'standard': + amax[2] = int(z_2) + + + # write to log + out_log.write("lowest x:\t%s\nhighest x:\t%s\nlowest y:\t%s\nhighest y:\t%s\nlowest z:\t%s\nhighest z:\t%s\n\n" + %(amin[0],amax[0],amin[1],amax[1],amin[2],amax[2])) + + return amin, amax + +# Function calculating range for sections +def range_sections(amin,amax,number_of_boxes_hor, number_of_boxes_ver,number_of_boxes_z): + # Find the highest absolute x because the boxes on the x axis has to be mirrored in the symmetry plane + if abs(amax[0]) < abs(amin[0]): # if the highest x is absolute smaller than the lowest x: + x_window_min_x = amin[0] # lowest x is lowest windows x + x_window_max_x= amin[0] *-1 # lowest x * -1 is highest windows x + steps_x = abs(amin[0]) / number_of_boxes_hor # length of x boxes + else: + x_window_min_x = amax[0] *-1 # highest x *-1 is lowest windows x + x_window_max_x = amin[0] # highest x is highest window x + steps_x = abs(amax[0]) / number_of_boxes_hor # length of x boxes + + #Calculating range of y boxes + steps_y = (abs(amax[1]) + abs(amin[1])) / number_of_boxes_ver + + #Calculating range of z boxes + steps_z = (abs(amax[2]) + abs(amin[2])) / number_of_boxes_z + + return steps_x, steps_y, steps_z, x_window_min_x,x_window_max_x + +# Function create list ranges of each section +def create_list_ranges(amin,amax,number_of_boxes_hor, number_of_boxes_ver,number_of_boxes_z, + x_window_min_x, x_window_max_x,steps_x,steps_y,steps_z): # creat list with ranges of each section + #creating the list with the ranges for every section + y_window_min = amin[1] + y_window_max = amin[1] + z_window_min = amin[2] + newlist = [] + sublist = [] + y_window_min2 = y_window_min + z_window_min2 = z_window_min + + # From front to back, from floor to ceiling + for a in range(0,int(number_of_boxes_ver)): #number of rows + x_window_min = x_window_min_x + x_window_min2 = x_window_min_x + steps_x + y_window_min = y_window_min2 + y_window_min2+= steps_y + z_window_min = amin[2] + z_window_min2 = z_window_min + for c in range(0,int(number_of_boxes_z)): # in z direction + z_window_min = z_window_min2 + z_window_min2 += steps_z + x_window_min = x_window_min_x + x_window_min2 = x_window_min_x + steps_x + for b in range(0, int((2*number_of_boxes_hor))): #number of columns + sublist.append(x_window_min) + sublist.append(x_window_min2) + sublist.append(y_window_min) + sublist.append(y_window_min2) + sublist.append(z_window_min) + sublist.append(z_window_min2) + newlist.append(sublist) + sublist = [] + x_window_min = x_window_min2 + x_window_min2 += steps_x + + return newlist + +# Function create list coordinates in sections +def create_list_coordinates(newlist): # create list coordinates in sections + #empty the list but maintaining the structure + newlist2=[] + newlist2sub = [] + newlistsub = [] + for x in range(0,len(newlist)): + for y in range(5,-1,-1): + newlist2sub.append(newlist[x].pop(y)) + newlist2.append(newlist2sub) + newlist2sub = [] + + # creating empty index list for storing line index of coordinates + indexlist = [] + for x in range(0,len(newlist)): + indexlist.append([]) + + return newlist2, indexlist + +# Function fill list sections with coordinates +def fill_sections_coordinates(newlist,newlist2,matrix,indexlist): # fill list sections with coordinates + #Filling the sections with the coordinates + for y in range(0,len(newlist2)): #notice the ranges are in reverse order of each section + rows_vertex = (numpy.where((matrix[:,0]>= float(newlist2[y][5])) + & (matrix[:,0] < float(newlist2[y][4])) + & (matrix[:,1] >= float(newlist2[y][3])) + & (matrix[:,1] < float(newlist2[y][2])) + & (matrix[:,2] >= float(newlist2[y][1])) + & (matrix[:,2] < float(newlist2[y][0])))) + + vertex_new = matrix[list(rows_vertex)] + newlist[y].append((vertex_new)) + indexlist[y].append((rows_vertex)) + + return newlist, indexlist + +# Function calculate mean and number of coordinates sections +def calc_sections(number_of_boxes_ver,number_of_boxes_hor,number_of_boxes_z,newlist): + #Selecting the sections on the symmetry axis and calculating the mean of every section + counter = 0 + counter2 = 0 + counter3 = 0 + row_count = 0 + difference = [] #whit sections which are empty + differencemean = [] #without sections which are empty + total = 0 + for x in range(0, (int(number_of_boxes_ver)*int(number_of_boxes_hor) * int(number_of_boxes_z))): # through all the lists/sections + if counter % (int (number_of_boxes_hor)*2) == 0: + counter = 0 + if counter3 == (int(number_of_boxes_hor)): + counter2 += (int(number_of_boxes_hor)) + counter3 = 0 + + first_pos = (newlist[counter2][0]) # coordinates of left section + last_pos = (newlist[counter2 + (int(number_of_boxes_hor)*2 -1)- counter][0]) # coordinates of right section + if len(first_pos) == 0: # If there are no coordinates in a section + number_coordinates_box1 = len(first_pos) + mean1 = 2000 + mean3 = '' + + else: + number_coordinates_box1 = len(first_pos) + for e in range(0,len(first_pos)): # If there are coordinates in a section + total += float(first_pos[e][0]) + mean1 = total/len(first_pos) + mean3 = mean1 + total = 0 + + if len(last_pos) == 0: # If there are no coordinates in a section + number_coordinates_box2 = len(last_pos) + mean2 = 1000 + mean4 = '' + + else: + number_coordinates_box2 = len(last_pos) + + for e in range(0,len(last_pos)): # If there are coordinates in a section + total += float(last_pos[e][0]) + mean2 = total/len(last_pos) + mean4 = mean2 + total = 0 + number_coordinates_box1 = len(first_pos) + number_coordinates_box2 = len(last_pos) + + # if number of coordinates in sections deviate to much, + if (abs(number_coordinates_box1 - number_coordinates_box2) > (number_coordinates_box1 / 2.0) + ) and(abs(number_coordinates_box1 - number_coordinates_box2) > (number_coordinates_box2 / 2.0)): + mean1 = 3000 + mean2 = 4000 + + difference.append(float(abs(mean2 + mean1))) #the means of every sections + counter += 2 + counter2 += 1 + counter3 += 1 + + try: + differencemean.append(abs(mean4 + mean3)) # the mean without empty sections + except: + continue + + return difference, differencemean + +# Function calculating the mean and the standard deviation +def calc_mean_stdev(differencemean, factor): + mean_percentage = numpy.mean(differencemean) #calculating the mean without empty sections + std_percentage = numpy.std(differencemean) # calculating the standard deviation of the list without the empty sections + std_percentage = std_percentage * factor + + return mean_percentage, std_percentage + +# Function calculating range of accepted differences +def calc_range(mean_percentage,std_percentage): + left_range = float(mean_percentage) - float(std_percentage) #left range mean minus one standard deviation + right_range = float(mean_percentage) + (float(std_percentage) )# right range mean plus one standard deviation + if left_range > 0: + left_range = 0 + return left_range, right_range + +# Function collecting incorrect coordinates +def collect_mirror_values(out_log, side, output, right_range, left_range,mean_percentage, + std_percentage, newlist, indexlist,difference, + number_of_boxes_hor,number_of_boxes_ver,number_of_boxes_z, + col_r, col_g, col_b): + + # Collecting the values of the sections which have to be mirrored + out_log.write('mean x difference:\t%s\nstandarddeviation mean x difference:\t%s\nrange:\t%s - %s\n\n'% + ((mean_percentage),(std_percentage),(left_range),(right_range))) + y = 0 + counter = 0 + counter2 = 0 + sub = [] + total = [] + left_range = "%.10f"%(left_range) + indexcount = 0 + facecount = 0 + totalcount = 0 + listindex = [] + listindexsub = [] + listindexcount = [] + listindextotalcount = [] + + for x in range(0, len(difference)): + difference_x = "%.10f"%(difference[x]) + if counter %(int(number_of_boxes_hor)) == 0 and x != 0: + counter = 0 + counter2 += (int(number_of_boxes_hor)) + counter += 1 + counter2 += 1 #left side counter + + # if left side of the object is correct + if side == 0: + if (float(difference_x) < float(left_range)) or (float(difference_x) > float(right_range)): + if len(newlist[counter2 -1][0]) != 0: + sub = (newlist[counter2 -1][0]) #coordinates in list + for y in range(0, len(sub)): + i = sub[y] + index_a = indexlist[counter2 -1][0][0][y] + try: + listindex.append(int(index_a)) # index number + listindextotalcount.append(int(totalcount)) + listindexsub.append(int(index_a)) # index number + listindexsub.append(totalcount) # index number of new coordinate + + except: + continue + listindexcount.append(listindexsub) # index number original + listindexsub = [] + number= float(sub[y][0]) * -1 # x coordinate * -1 for mirroring at symmetry plane + sub[y][0] = number + total.append(sub[y]) + totalcount += 1 + # writing new coordinates to output + output.write('%s %s %s %s %s %s\n'%(sub[y][0], float(sub[y][1]), float(sub[y][2]),(col_r), (col_g), (col_b))) + + + # if right side of the object is correct + if side == 1: + dif = (int(number_of_boxes_hor)) - counter + counter3 = counter2 + (dif*2) + 1 # right side counter + + if (float(difference[x] )> float(right_range)) or (float(difference[x]) < float(left_range)): + if len(newlist[counter3 -1][0]) != 0: + sub = (newlist[counter3 -1][0]) # coordinates in list + for y in range(0, len(sub)): + i = sub[y] + index_a = indexlist[counter3 -1][0][0][y] + try: + listindex.append(int(index_a))# index number + listindextotalcount.append(int(totalcount)) + listindexsub.append(int(index_a)) # index number + listindexsub.append(totalcount)# index number of new coordinate + except: + continue + listindexcount.append(listindexsub)# index number original + listindexsub = [] + number= float(sub[y][0]) * -1 # x coordinate * -1 for mirroring at symmetry plane + sub[y][0] = number + total.append(sub[y]) + totalcount += 1 + # writing new coordinates to output + output.write('%s %s %s %s %s %s\n'%(sub[y][0], float(sub[y][1]), float(sub[y][2]),(col_r), (col_g), (col_b))) + + return total, listindexcount, listindex, listindextotalcount,totalcount + +# Function extract original faces for corrected coordinates +def find_original_faces(matrix2, matrix3, listindex): + # finding the original faces of kopied points + ix = numpy.in1d(matrix2.ravel(),listindex).reshape(matrix2.shape) # finding where the listindex is the same as the matrix + rows, cols = numpy.where(ix) # finding index of the faces + m_face = matrix2[list(set(rows))] # extracting the rows with those faces + m_face_square = [] + # if square faces present + if matrix3[-1][0] != 0: + ix_square = numpy.in1d(matrix3.ravel(),listindex).reshape(matrix3.shape) # finding where the listindex is the same as the matrix + rows_square, cols3 = numpy.where(ix_square) # finding index of the faces + m_face_square = matrix3[list(set(rows_square))] # extracting the rows with those faces + return m_face, m_face_square + +# Function replacing index values +def replace_index_values(m_face,listindex,listindextotalcount): + c = [] + array_oi = numpy.array(m_face) + maxn = numpy.amax(array_oi) + palette = list(range(int(maxn))) + palette = numpy.array(palette) + key = numpy.array(listindex) + listindextotalcount = numpy.array(listindextotalcount) + + # sorting the lists which have to be replaced + order= numpy.argsort(key) # sorting of the array key + litc_sorted = listindextotalcount[order] # sorting the palette + key_sorted = key[order] # sorting the key + + palette2 = list(range(int(maxn))) # creating list with all possible values in it with max the highest number in the faces + key2 = [] + + counter = 0 + for item in range(0,int(maxn)): # make array containing what have to be changed + if (counter < len(key_sorted)) and item == (key_sorted[counter]): + key2.append(litc_sorted[counter]) + counter += 1 + else: + key2.append(numpy.nan) + + key2 = numpy.array(key2) # converten of list to array + index_f = numpy.digitize(array_oi.reshape(-1,), palette2)-1 # indexing which numbers has to be changed + out_sub =(key2[index_f].reshape(array_oi.shape)) # new array creating with changed values + out_faces = out_sub[~numpy.isnan(out_sub).any(axis=1)] # extracting the faces with not changed values, (nan) + + return out_faces + +# Function write output +def write_output(name_file_ply, out_log, totalcount, var_header, var_vertex_nm, var_face_nm,out_faces,out_faces_2): # write output + outfile2 = open(str(sys.argv[2]), 'w')#writing the output file + file1 = open(name_file_ply) + g = 0 + for d in range(0,(int(var_header) + int(var_vertex_nm) + int(var_face_nm) + 1)): + line2 = file1.readline().strip() + readline2 = line2.strip().split() + if readline2[0] == "element" and readline2[1] == "vertex": + outfile2.write("element vertex %s\n"%(int(var_vertex_nm) + totalcount)) #number of vertexen changing + + elif readline2[0] == "element" and readline2[1] == "face": + outfile2.write("element face %s\n"%(int(var_face_nm) + len(out_faces) + len(out_faces_2)))#number of faces changing + + elif (int(var_header) < d < (int(var_header) + int(var_vertex_nm))): #rotated vertexen, with original color code + outfile2.write('%s\n'%(line2)) + g += 1 + + elif (int(var_header) + int(var_vertex_nm)) == d: #for new vertexen, with color code 0,0,0 + outfile2.write('%s\n'%(line2)) #the last original vertex + g += 1 + with open('new_coordinates2.ply') as infile: + for line in infile: + outfile2.write(line) + + else: #everything left + outfile2.write('%s\n'%(line2)) + + #writing new faces to output + for z in range(0,len(out_faces)): + outfile2.write("3 %s %s %s\n"%((int(out_faces[z][0])+ int(var_vertex_nm)), + (int(out_faces[z][1])+ int(var_vertex_nm)), (int(out_faces[z][2])+ int(var_vertex_nm)))) + #writing new square faces to output + if len(out_faces_2) != 0: + for x in range(0,len(out_faces_2)): + outfile2.write("4 %s %s %s %s\n"%((int(out_faces_2[z][0])+ int(var_vertex_nm)), + (int(out_faces_2[z][1])+ int(var_vertex_nm)), (int(out_faces_2[z][2])+ int(var_vertex_nm)), + (int(out_faces_2[z][3] + int(var_vertex_nm))))) + + out_log.write('number of reconstructed faces:\t%s\n\n'%(len(out_faces) + len(out_faces_2))) + outfile2.close() + + +main() + + +