#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Mon Aug 5 15:56:47 2019 @author: pavel """ """ Class that extends the vispy SceneCanvas to draw 3D tubes """ from vispy import gloo, scene, app from vispy.gloo import set_viewport, set_state, clear, set_blend_color, context from vispy.gloo import gl from vispy.util.transforms import perspective, translate, rotate, scale import vispy.gloo.gl as glcore from vispy.util.quaternion import Quaternion import numpy as np import math import network_dep as nwt from tube_shaders import FRAG_SHADER, VERT_SHADER DEBUG = False if DEBUG: from mpl_toolkits.mplot3d import Axes3D import matplotlib import matplotlib.pyplot as plt class TubeDraw(scene.SceneCanvas): #sigUpdate = QtCore.pyqtSignal(float, float, float) #Initiates the canvas. def __init__(self, **kwargs): #Initialte the class by calling the superclass scene.SceneCanvas.__init__(self, size=(512,512), keys='interactive', **kwargs) #unfreeze the drawing area to allow for dynamic drawing and interaction self.unfreeze() self.edge_dict = {} self.select = False #generate dummy buffers for the meshes self.program = gloo.Program(VERT_SHADER, FRAG_SHADER) self.cylinder_data = np.zeros(5*5, dtype=[('a_position', np.float32, 3), ('a_normal', np.float32, 3), ('a_fg_color', np.float32, 4), ('a_selection', np.float32, 1), #('a_linewidth', np.float32, 1), ]) self.triangle_data = np.random.randint(size=(5, 3), low=0, high=(4-1)).astype(np.uint32) self.vbo = gloo.VertexBuffer(self.cylinder_data) self.triangles = gloo.IndexBuffer(self.triangle_data) #generate an index buffer for the caps. self.cap_data = np.random.randint(size=(5, 3), low=0, high=(4-1)).astype(np.uint32) self.caps = gloo.IndexBuffer(self.cap_data) self.program.bind(self.vbo) self.scale = [1,1,1] self.r1 = np.eye(4, dtype=np.float32) self.r2 = np.eye(4, dtype=np.float32) set_viewport(0,0,*self.physical_size) #set_state(clear_color='white', depth_test=True, blend=True, # blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('less'), cull_face='back') set_state(clear_color='white', depth_test=True, blend=True, blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal')) #set_blend_color(color='black') #set_state('translucent') #self.program['u_LightPos'] = [0., 0., -1000.] #self.camera = self.central_widget.add_view() #self.camera.camera = 'turntable' self.down = False self.camera = np.asarray([0.0, 0.0, 200.0], dtype=np.float32) self.up = np.asarray([0., 1., 0.], dtype=np.float32) #self.init_camera = [0.,0.,1000.] ##### prototype ##### #Set the visualization matrices self.program['u_selection'] = 0.0 self.program['u_eye'] = self.camera self.program['u_up'] = self.up self.program['u_target'] = np.asarray([0., 0., 0.], dtype=np.float32) # self.bb = np.ones((26, 3), dtype=np.float32) # self.program['u_bb'] = self.bb #Load the data necessary to draw all of the microvessels def set_data(self, G, bbu, bbl, num_sides): self.G = G self.num_sides = num_sides self.bbu = bbu self.bbl = bbl bb = nwt.AABB(G).resample_sides(3) print(self.bbu, self.bbl, self.camera) #create program self.gen_cylinder_vbo(self.G, self.num_sides) self.vbo = gloo.VertexBuffer(self.cylinder_data) #self.triangle_data = np.append(self.triangle_data, self.cap_data, axis=0) self.triangles = gloo.IndexBuffer(self.triangle_data) self.caps = gloo.IndexBuffer(self.cap_data) #self.view = np.eye(4, dtype=np.float32) self.model = np.eye(4, dtype=np.float32) self.projection = np.eye(4, dtype=np.float32) self.projection = perspective(90.0, self.physical_size[0]/self.physical_size[1], 1.0, 1000.0) #self.projection = perspective(90.0, 1.0, -1.0, 1.0) self.program['u_model'] = self.model #self.program['u_LightPos'] = [0., 0., 1000.] #self.program['u_view'] = self.view self.program['u_projection'] = self.projection self.program.bind(self.vbo) print(self.bbu, self.bbl, self.camera) gloo.set_clear_color('white') self.center = (bbu-bbl)/2.0 print(self.bbu, self.bbl, self.camera) self.translate = [-self.center[0], -self.center[1], -self.center[2]] self.bb = np.ones((26, 3), dtype=np.float32) for i in range(26): for j in range(3): self.bb[i,j] = bb[i][j] self.program['u_bb'] = self.bb if DEBUG: print('bb is ', self.bb) # for i in range(len(self.translate)): # self.camera[i] += self.translate[i] ##### prototype ##### self.camera = self.camera - self.translate print(self.camera) self.program['u_eye'] = self.camera self.up = np.cross((np.asarray(self.center, dtype=np.float32)-np.asarray(self.camera, dtype=np.float32)), np.asarray(self.up)) self.program['u_up'] = self.up self.program['u_target'] = self.translate #self.show() #Called during resize of the window in order to redraw the same image in the #larger/smaller area. def on_resize(self, event): width, height = event.physical_size gloo.set_viewport(0, 0, width, height) if DEBUG: print(self.physical_size) #overloaded function called during the self.update() call to update the current #frame using the GLSL frag/vert shaders def on_draw(self, event): clear(color='white', depth=True) gloo.set_clear_color('white') if self.select: set_state(clear_color='white', depth_test=True, blend=True, blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('always'), cull_face='back') else: set_state(clear_color='white', depth_test=True, blend=True, blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal'), cull_face='back') self.program.draw('triangles', self.triangles) if self.select: set_state(clear_color='white', depth_test=True, blend=True, blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('always'), cull_face=False) else: set_state(clear_color='white', depth_test=True, blend=True, blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal'), cull_face=False) self.program.draw('triangles', self.caps) self.projection = perspective(90.0, self.physical_size[0]/self.physical_size[1], 1.0, 1000.0) self.program['u_projection'] = self.projection def refresh(self): self.update() app.Canvas.update(self) #Creates a cylinder around ever segment in the microvascular network. def gen_cylinder_vbo(self, G, num_sides = 32): i = 0 num_pts = 0 num_tri = 0 num_tri_caps = 0 for e in G.edges(): num_pts += len(self.G.edge_properties["x"][e]) num_tri += (len(self.G.edge_properties["x"][e])-1)*num_sides*2 num_tri_caps += (2*(num_sides-2)) self.cylinder_data = np.zeros(num_pts*num_sides, dtype=[('a_position', np.float32, 3), ('a_normal', np.float32, 3), ('a_fg_color', np.float32, 4), ('a_selection', np.float32, 1), ]) self.triangle_data = np.random.randint(size=(num_tri, 3), low=0, high=(G.num_edges()-1)).astype(np.uint32) self.cap_data = np.random.randint(size = (num_tri_caps, 3), low = 0, high=(G.num_edges()-2)).astype(np.uint32) index = 0 t_index = 0 c_index = 0 #for each edge generate a cylinder. for e in G.edges(): #print("done") #for each fiber get all the points and the radii X = self.G.edge_properties["x"][e] Y = self.G.edge_properties["y"][e] Z = self.G.edge_properties["z"][e] R = self.G.edge_properties["r"][e] color = G.edge_properties["RGBA"][e] pts = np.array([X,Y,Z]).T circle_pts = np.zeros((pts.shape[0], num_sides, 3), dtype = np.float32) step = 2*np.pi/num_sides # U = np.zeros(pts.shape, dtype=np.float32) # V = np.zeros(pts.shape, dtype=np.float32) # direction = np.zeros(pts.shape, dtype=np.float32) #for every point in the edge for p in range(pts.shape[0]): #if first point, generate the circles. #In this case we want to generate a cap if we see the first circle or the last. if(p == 0): #get the direction direction = (pts[p+1] - pts[p]) #normalize direction direction = direction/np.sqrt(direction[0]*direction[0] + direction[1]*direction[1] + direction[2]*direction[2]) #generate a vector to use as an element of the cross product Y = np.zeros((3,), dtype = np.float32) Y[0] = 1. #if direction and Y are parallel, change Y if(np.cos(np.dot(Y, direction)) < 0.087): Y[0] = 0.0 Y[2] = 1.0 #generate first plane vector U = np.cross(direction, Y) U = U/np.sqrt(U[0]*U[0] + U[1]*U[1] + U[2]*U[2]) #generate second plane vector V = np.cross(direction, U) V = V/np.sqrt(V[0]*V[0] + V[1]*V[1] + V[2]*V[2]) #print(R[p],pts[p]) #generate circle. for s in range(num_sides): circle_pts[p][s][0] = R[p]*np.cos(s*step)*V[0]*0.5 + R[p]*np.sin(s*step)*U[0]*0.5 circle_pts[p][s][1] = R[p]*np.cos(s*step)*V[1]*0.5 + R[p]*np.sin(s*step)*U[1]*0.5 circle_pts[p][s][2] = R[p]*np.cos(s*step)*V[2]*0.5 + R[p]*np.sin(s*step)*U[2]*0.5 #if last point, copy the previous circle. elif(p == pts.shape[0]-1): for s in range(num_sides): circle_pts[p][s] = circle_pts[p-1][s] for v in range(pts.shape[0]): circle_pts[v,:,0] += pts[v][0] circle_pts[v,:,1] += pts[v][1] circle_pts[v,:,2] += pts[v][2] #otherwise, rotate the circle else: #print(R[p], pts[p]) #generate a new direction vector. dir_new = (pts[p+1]-pts[p]) dir_new = dir_new/np.sqrt(dir_new[0]*dir_new[0] + dir_new[1]*dir_new[1] + dir_new[2]*dir_new[2]) dir_new2 = (pts[p]-pts[p-1]) dir_new2 = dir_new2/np.sqrt(dir_new2[0]*dir_new2[0] + dir_new2[1]*dir_new2[1] + dir_new2[2]*dir_new2[2]) dir_new = dir_new+dir_new2 dir_new = dir_new/np.sqrt(dir_new[0]*dir_new[0] + dir_new[1]*dir_new[1] + dir_new[2]*dir_new[2]) #print(dir_new, direction) #generate the quaternion rotation vector for the shortest arc k = 1.0 + np.dot(direction, dir_new) s = 1/np.sqrt(k+k) r = s*np.cross(direction, dir_new) theta = k*s #r = np.cross(direction, dir_new) #r = r/np.sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]) #calculate the degree of quaternion rotation. #cos_theta = np.sqrt(np.sqrt(np.dot(dir_new, dir_new)) * np.sqrt(np.dot(direction, direction))) + np.dot(dir_new, direction) #cos_theta = np.dot(direction, dir_new) #theta = np.arccos(cos_theta)/2.0 #print(r, cos_theta, theta) #quat = np.append(theta, r) q = Quaternion(w=theta, x = r[0], y = r[1], z = r[2]).normalize() #print(quat) #q = np.quaternion(quat[0], quat[1], quat[2], quat[3]) #rot = Rotation.from_quat(quat, normalized=False) #rot.as_quat() for s in range(num_sides): circle_pts[p][s] = q.rotate_point(circle_pts[p-1][s]) #circle_pts[p][s] = rot.apply(circle_pts[p-1][s]) #circle_pts[p][s] = q.rotate(circle_pts[p-1][s]) #circle_pts[p][s] = np.quaternion.rotate_vectors(q, circle_pts[p][s]) #circle_pts[p][s] = q.rotate_vectors(q, circle_pts[p][s]) #circle_pts[p][s] = circle_pts[p][s] direction = dir_new #generate the triangles triangles = np.random.randint(size=((pts.shape[0]-1)*2*(num_sides), 3), low=0, high=(G.num_edges()-1)).astype(np.uint32) cap_triangles = np.random.randint(size=((num_sides-2)*2 , 3), low=0, high=(num_sides-2)).astype(np.uint32) t_index_temp = 0 c_index_temp = 0 #for every ring in the fiber for ring in range(0, pts.shape[0], pts.shape[0]-1): #if we have the first ring or the last ring add a cap. idx_ori = index+ring*num_sides for side in range(0, num_sides-2): idx_current_point = index+ring*num_sides + side if(side < num_sides-2): cap_tri = [idx_ori, idx_current_point+1, idx_current_point+2] else: cap_tri = [idx_ori, idx_current_point+1, idx_ori] cap_triangles[c_index_temp] = cap_tri self.cap_data[c_index] = cap_tri c_index_temp += 1 c_index += 1 for ring in range(pts.shape[0]-1): #otherwise generate the sides. for side in range(num_sides): if(side < num_sides-1): idx_current_point = index+ring*num_sides + side idx_next_ring = index + (ring+1) * num_sides + side triangle1 = [idx_current_point, idx_next_ring, idx_next_ring+1] triangle2 = [idx_next_ring+1, idx_current_point+1, idx_current_point] triangles[t_index_temp] = [idx_current_point, idx_next_ring, idx_next_ring+1] triangles[t_index_temp+1] = [idx_next_ring+1, idx_current_point+1, idx_current_point] self.triangle_data[t_index] = triangle1 self.triangle_data[t_index+1] = triangle2 t_index += 2 t_index_temp += 2 else: idx_current_point = index + ring*num_sides + side idx_next_ring = index + (ring+1)*num_sides + side triangle1 = [idx_current_point, idx_next_ring, idx_next_ring-num_sides+1] triangle2 = [idx_next_ring-num_sides+1, idx_current_point-num_sides+1, idx_current_point] triangles[t_index_temp] = [idx_current_point, idx_next_ring-num_sides, idx_next_ring-num_sides+1] triangles[t_index_temp+1] = [idx_next_ring-num_sides+1, idx_current_point-num_sides+1, idx_current_point] self.triangle_data[t_index] = triangle1 self.triangle_data[t_index+1] = triangle2 t_index += 2 t_index_temp += 2 #generate the points data structure circle_pts_data = circle_pts.reshape((pts.shape[0]*num_sides, 3)) self.cylinder_data['a_position'][index:(pts.shape[0]*num_sides+index)] = circle_pts_data self.cylinder_data['a_fg_color'][index:(pts.shape[0]*num_sides+index)] = color self.cylinder_data['a_selection'][index:(pts.shape[0]*num_sides+index)] = 0.0 #generate the normals data structure pts_normals = circle_pts.copy() for p in range(pts.shape[0]): pts_normals[p][:] = pts_normals[p][:] - pts[p] for s in range(num_sides): pts_normals[p][s] = pts_normals[p][s]/np.sqrt(pts_normals[p][s][0]*pts_normals[p][s][0] \ + pts_normals[p][s][1]*pts_normals[p][s][1] + pts_normals[p][s][2]*pts_normals[p][s][2]) self.cylinder_data['a_normal'][index:(pts.shape[0]*num_sides+index)] = \ pts_normals.reshape((pts.shape[0]*num_sides, 3)) index += pts.shape[0]*num_sides self.edge_dict[(int(e.source()), int(e.target()))] = (index-pts.shape[0]*num_sides, index) #Add the caps for each of the endpoints. if DEBUG: if(i == 2): fig = plt.figure() ax = fig.add_subplot(111, projection='3d') #ax.scatter(circle_pts[:,:,0], circle_pts[:,:,1], circle_pts[:,:,2]) ax.plot(pts[:,0], pts[:,1], pts[:,2]) for j in range(pts.shape[0]): ax.plot(circle_pts[j,:,0], circle_pts[j,:,1], circle_pts[j,:,2]) for j in range(cap_triangles.shape[0]): tri = np.zeros((3,4)) tri[:,0] = self.cylinder_data['a_position'][cap_triangles[j][0]] tri[:,1] = self.cylinder_data['a_position'][cap_triangles[j][1]] tri[:,2] = self.cylinder_data['a_position'][cap_triangles[j][2]] tri[:,3] = self.cylinder_data['a_position'][cap_triangles[j][0]] ax.plot(tri[0,:], tri[1,:], tri[2,:], c='b') for j in range(triangles.shape[0]): tri = np.zeros((3,4)) tri[:,0] = self.cylinder_data['a_position'][triangles[j][0]] tri[:,1] = self.cylinder_data['a_position'][triangles[j][1]] tri[:,2] = self.cylinder_data['a_position'][triangles[j][2]] tri[:,3] = self.cylinder_data['a_position'][triangles[j][0]] ax.plot(tri[0,:], tri[1,:], tri[2,:], c='b') for j in range(triangles.shape[0]): tri = np.zeros((3,3)) tri[:,0] = self.cylinder_data['a_position'][triangles[j][0]] tri[:,1] = self.cylinder_data['a_position'][triangles[j][1]] tri[:,2] = self.cylinder_data['a_position'][triangles[j][2]] norm = np.zeros((3,3)) norm[:,0] = self.cylinder_data['a_normal'][triangles[j][0]] norm[:,1] = self.cylinder_data['a_normal'][triangles[j][1]] norm[:,2] = self.cylinder_data['a_normal'][triangles[j][2]] ax.quiver(tri[0,:], tri[1,:], tri[2,:], norm[0,:], norm[1,:], norm[2,:], colors = 'r') plt.show() i+=1 #create the data. def select_edges(self, edges): #gloo.context.set_current_canvas(gloo.context.get_current_canvas()) if len(edges) > 0: self.select = True #gloo.set_depth_func('always') #clear(color='white', depth=True) #gl.GL_DEPTH_FUNC('always') #print(gloo.wrappers.get_gl_configuration()) self.program['u_selection'] = 1.0 for e in edges: if e in self.edge_dict: idx = self.edge_dict[e] #print(self.canvas.edge_dict[e]) elif (e[1], e[0]) in self.edge_dict: idx = self.edge_dict[(e[1],e[0])] else: print("WHAT THE FUCK HAPPENED") self.cylinder_data["a_selection"][idx[0]:idx[1]] = 1.0 self.vbo = gloo.VertexBuffer(self.cylinder_data) self.program.bind(self.vbo) #self.update() else: self.program['u_selection'] = 0.0 self.select = False #set_state(clear_color='white', depth_test=True, blend=True, # blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal')) #gloo.set_depth_func('lequal') #clear(color='white', depth=True) #print(gloo.wrappers.get_gl_configuration()) self.refresh() #print(self.canvas.edge_dict[(e[1], e[0])]) #Handles the mouse wheel event, i.e., zoom def on_mouse_wheel(self, event): # self.scale[0] = self.scale[0] + self.scale[0]*event.delta[1]*0.05 # self.scale[1] = self.scale[1] + self.scale[1]*event.delta[1]*0.05 # self.scale[2] = self.scale[2] + self.scale[2]*event.delta[1]*0.05 ## self.view[0][0] = self.scale[0] ## self.view[1][1] = self.scale[1] ## self.view[2][2] = self.scale[2] # print("in mouse wheel ", self.r1, self.r2) # self.view = np.matmul(translate((self.translate[0], self.translate[1], self.translate[2])), self.r1) # self.view = np.matmul(self.view, self.r2) # self.view = np.matmul(self.view, scale((self.scale[0], self.scale[1], self.scale[2]))) # #self.view = np.matmul(self.view, self.r1) # #self.view = np.matmul(self.view, self.r2) # #self.view = np.matmul(translate((self.translate[0], self.translate[1], self.translate[2])), scale((self.scale[0], self.scale[1], self.scale[2]))) # #self.rotate = [0., 0.] # #self.camera = # # #self.view[1][1] = self.view[1][1]+self.view[1][1]*event.delta[1]*0.05 # #print(self.view[0][0], " ",self.view[1][1]) # #print(self.view) # self.camera = [0.0, 0.0, -100.0, 1.0] ## for i in range(len(self.translate)): ## self.camera[i] += self.translate[i] # self.program['u_view'] = self.view # if(np.asarray(self.camera).all() == np.asarray([0., 0., 0.]).all()): # self.camera = np.asarray([0., 0., 0.]) # else: direction = np.asarray(self.translate) - np.asarray(self.camera) direction = direction/np.sqrt(np.dot(direction, direction)) for i in range(3): self.camera[i] = self.camera[i] + direction[i]*event.delta[1]*2.0 self.program['u_eye'] = self.camera #print(self.view) #print(event.delta[1]) self.refresh() def update_view(self, event): self.location = event.pos #self.program['u_view'] = self.view self.down = True #Handles the mouse press event to rotate the camera if the left mouse button #if clicked down. def on_mouse_press(self, event): if(event.button == 1): self.update_view(event) #Handles the rotation of the camera using a quaternion centered around the #focus point. def on_mouse_move(self, event): if(self.down == True): coord = self.transforms.get_transform('canvas', 'render').map(self.location)[:2] coord2 = self.transforms.get_transform('canvas', 'render').map(event.pos)[:2] #coord[1] = 0 #coord2[1] = 0 phi = (coord[0]-coord2[0])*360.0/2.0/np.pi theta = (coord[1]-coord2[1])*360.0/2.0/np.pi if DEBUG: print(theta*360.0/2.0/np.pi, -phi*360.0/2.0/np.pi) self.camera = self.camera - self.translate q1 = Quaternion.create_from_axis_angle(angle=phi, ax=0.0, ay=1.0, az=0.0, degrees=True) q2 = Quaternion.create_from_axis_angle(angle=theta, ax=1.0, ay=0.0, az=0.0, degrees=True) #q1 = Quaternion(w=theta, x = 0, y = 1, z = 0).inverse().normalize() #q2 = Quaternion(w=-phi, x = 1, y = 0, z = 0).inverse().normalize() q = q1*q2 # #print("Angle in Degrees is ", theta, " ", phi, coord[0] - coord2[0]) # self.r1 = rotate(theta, (0, 1, 0)) # self.r2 = rotate(-phi, (1, 0, 0)) # # print("in on mouse move ", self.r1, self.r2) # # self.view = np.matmul(self.view, self.r1) # #print("r1", self.view) # self.view = np.matmul(self.view, self.r2) # #print("r2", self.view) # ## self.view = np.matmul(self.view, q1.get_matrix().T) ## self.view = np.matmul(self.view, q2.get_matrix().T) # # self.program['u_view'] = self.view # self.location = event.pos # #print("Angle in Degrees is ", theta, " ", phi) # #print(self.camera) self.camera = np.asarray(q.rotate_point(self.camera), dtype=np.float) self.camera = self.camera + self.translate self.up = np.asarray(q.rotate_point(self.up), dtype=np.float) self.up = self.up/np.sqrt(np.dot(self.up, self.up)) #self.rotate[0] = self.rotate[0] + theta #self.rotate[1] = self.rotate[1] + phi #print(self.rotate)f #radius = np.sqrt(np.dot(self.center, self.center))*2 #test = np.linalg.inv(self.view).T #print(test) #self.camera = sph2cart(self.rotate[0]/360.0*2.0*np.pi+np.pi, self.rotate[1]/360.0*2.0*np.pi+np.pi, 1000.0) #self.camera[0] = self.camera[0] + self.center[0] #self.camera[1] = self.camera[1] + self.center[1] #self.camera[2] = self.camera[2] - self.center[2] if DEBUG: print("current position ", self.camera, " and up vector ", self.up) self.program['u_eye'] = self.camera self.program['u_up'] = self.up #self.program['u_LightPos'] = [self.camera[0], self.camera[1], self.camera[2]] self.refresh() #reverts the mouse state during release. def on_mouse_release(self, event): self.down = False