From 950c425b8f128b628e5304fe7f35268907fdce77 Mon Sep 17 00:00:00 2001 From: yliu88au <75292881+yliu88au@users.noreply.github.com> Date: Tue, 24 May 2022 21:45:58 +1000 Subject: [PATCH] Rewrite plot_geometry.py and demos Rewrite plot_geometry.py with many new features: ie., flexible geometry display with animations --- Python/demos/d01_CreateGeometry.py | 11 +- Python/tests/test_plot_geometry.py | 74 +++ .../utilities/visualization/plot_geometry.py | 484 ++++++++++-------- 3 files changed, 355 insertions(+), 214 deletions(-) create mode 100644 Python/tests/test_plot_geometry.py diff --git a/Python/demos/d01_CreateGeometry.py b/Python/demos/d01_CreateGeometry.py index 883b4cbf..278ea36e 100644 --- a/Python/demos/d01_CreateGeometry.py +++ b/Python/demos/d01_CreateGeometry.py @@ -24,6 +24,7 @@ #%% Import statements import tigre import numpy as np +from matplotlib import pyplot as plt #%% Geometry definition # @@ -129,4 +130,12 @@ #%% Plot your geometry geo = tigre.geometry_default() # Default cone beam geometry -tigre.plot_geometry(geo, angle=-np.pi / 6) +tigre.plot_geometry(geo) + +# animation +geo = tigre.geometry_default() # Default cone beam geometry +angles=np.linspace(0,np.pi,50) # half circle +ani = tigre.plot_geometry(geo,angles,animate=True,fname='d01_Create_Geometry') # a *.gif or *.mp4 file will be saved +plt.show() + + diff --git a/Python/tests/test_plot_geometry.py b/Python/tests/test_plot_geometry.py new file mode 100644 index 00000000..6a3e8654 --- /dev/null +++ b/Python/tests/test_plot_geometry.py @@ -0,0 +1,74 @@ + +## test + +import numpy as np +from matplotlib import pyplot as plt +import tigre + +# # plot 1, check zero angle position, scan rotation direction, detector offset + +# geo = tigre.geometry(nVoxel=np.array([256,128,256]),default=True) +# geo.sVoxel = geo.nVoxel * np.array([1,1,1]) # (z,y,x) +# geo.nDetector = np.array([256,128]) +# geo.dDetector = np.array([0.8, 0.8])*2 +# geo.sDetector = geo.dDetector * geo.nDetector +# geo.offDetector=np.array([0,0]) # viewing from S, D move (up, right) => (v, u) +# geo.rotDetector=np.array([30,0,0])/180*np.pi # [roll, pitch, yaw] viewing from S to D +# geo.offOrigin = np.array([0,0,0]) # (z,y,x) +# geo.COR=0 +# angles=np.linspace(0,np.pi,100) +# ani1=tigre.plot_geometry(geo,angles,10,animate=True) # angle=0, S is at (x=DSO, y=0, z=0) +# ani1 +# # confirm the plot with projection +# from scipy.io import loadmat +# head=loadmat('head.mat')['img'].transpose(2,1,0).copy() +# head=head[:,:128,:].copy() +# proj = tigre.Ax(head,geo,angles) +# plt.figure() +# plt.subplot(1,2,1) +# plt.imshow(head[:,:,128],origin='lower') +# plt.title('dim2=128') +# plt.ylabel('dim0 ->') +# plt.xlabel('dim1 ->') +# plt.subplot(1,2,2) +# plt.imshow(proj[0,:,:],origin='lower') +# plt.ylabel('v ->') +# plt.xlabel('u ->') + + +# plot 2, check staticDetectorGeo() for tomosymthesis setup + +geo = tigre.geometry_default() +angles=np.linspace(-60,60,31)/180*np.pi +geo = tigre.staticDetectorGeo(geo,angles,60) + +ani2=tigre.plot_geometry(geo,angles,0,animate=True,fname='Tomosynthesis') +ani2 + + +## plot 3, fixed target object and detector positions and orientations, source moving linearly + +geo = tigre.geometry_default() +df = np.linspace(-510,510,64) # source position on +geo.DSO = 750 +geo.DSD = 1000 + +geo, angles = tigre.staticDetLinearSourceGeo(geo,df,10,60) + +ani3=tigre.plot_geometry(geo, angles, 0, animate=True, fname='Linear_Tomosynthesis') +ani3 + + +## plot 4, helical CT +geo = tigre.geometry_default(high_resolution=False) + +angles = np.linspace(0, 2 * np.pi, 100) +angles = np.hstack([angles, angles, angles]) # loop 3 times + +# This makes it helical, axis order (z,y,x) for python +geo.offOrigin = np.zeros((angles.shape[0], 3)) +geo.offOrigin[:, 0] = np.linspace( + -1024 / 2 + 128, 1024 / 2 - 128, angles.shape[0]) + +ani4 = tigre.plot_geometry(geo, angles, 0, animate=True, fname='Helical_CT') +ani4 diff --git a/Python/tigre/utilities/visualization/plot_geometry.py b/Python/tigre/utilities/visualization/plot_geometry.py index 9eb65739..0838e6dd 100644 --- a/Python/tigre/utilities/visualization/plot_geometry.py +++ b/Python/tigre/utilities/visualization/plot_geometry.py @@ -1,217 +1,273 @@ -import matplotlib.patches -import numpy as np -import tigre -from mpl_toolkits.mplot3d import art3d - -# https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html - - -class Arrow3D(matplotlib.patches.FancyArrowPatch): - def __init__(self, xs, ys, zs, *args, **kwargs): - matplotlib.patches.FancyArrowPatch.__init__(self, (0, 0), (0, 0), *args, **kwargs) - self._verts3d = xs, ys, zs - - def draw(self, renderer): - from mpl_toolkits.mplot3d import proj3d - - xs3d, ys3d, zs3d = self._verts3d - xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M) - self.set_positions((xs[0], ys[0]), (xs[1], ys[1])) - matplotlib.patches.FancyArrowPatch.draw(self, renderer) - - -ROT_DEFAULT = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) -TRANS_DEFAULT = np.array([0, 0, 0]) - - -def pathpatch_2d_to_3d_affine(pathpatch, mat_rot=ROT_DEFAULT, vec_trans=TRANS_DEFAULT): - """ - Transforms a 2D Patch to a 3D patch using the affine tranform - of the given rotation matrix and translation vector. - The pathpatch is assumed to be on the plane Z = 0. - """ - path = pathpatch.get_path() # Get the path and the associated transform - trans = pathpatch.get_patch_transform() - - path = trans.transform_path(path) # Apply the transform - - pathpatch.__class__ = art3d.PathPatch3D # Change the class - pathpatch._code3d = path.codes # Copy the codes - pathpatch._facecolor3d = pathpatch.get_facecolor # Get the face color - - verts = path.vertices # Get the vertices in 2D - - M = np.array( - [ - [mat_rot[0, 0], mat_rot[0, 1], mat_rot[0, 2], vec_trans[0]], - [mat_rot[1, 0], mat_rot[1, 1], mat_rot[1, 2], vec_trans[1]], - [mat_rot[2, 0], mat_rot[2, 1], mat_rot[2, 2], vec_trans[2]], - ] - ) - - pathpatch._segment3d = np.array([np.dot(M, (x, y, 0, 1)) for x, y in verts]) - - -def plot_geometry(geo, angle=0): - """Plots the given geometry.""" - import mpl_toolkits.mplot3d.art3d as art3d - import matplotlib.pyplot as plt - - fig = plt.figure(figsize=(6, 6)) - fig.suptitle("Cone Beam Compute Tomography geometry") - ax = fig.add_subplot(111, projection="3d") - ax.set_title("Current CBCT geometry, in scale") - - limXY = max(geo.DSO, geo.DSD - geo.DSO) - limZ = geo.sVoxel[0] - - ax.set_box_aspect((limXY, limXY, limZ)) - - ax.set_xlabel("X") - ax.set_ylabel("Y") - ax.set_zlabel("Z") - ax.set_xlim3d(-limXY * 1.2, limXY * 1.2) - ax.set_ylim3d(-limXY * 1.2, limXY * 1.2) - ax.set_zlim3d(-limZ * 1.2, limZ * 1.2) - - # Trajectory of Source - # https://matplotlib.org/devdocs/api/_as_gen/matplotlib.patches.Circle.htm - circ = matplotlib.patches.Circle((0, 0), geo.DSO, color="black", fill=False, ls="-.", lw=0.5) - ax.add_patch(circ) - art3d.pathpatch_2d_to_3d(circ, z=0, zdir="z") - - # Trajectory of Detector - circ = matplotlib.patches.Circle( - (0, 0), geo.DSD - geo.DSO, color="black", fill=False, ls="-.", lw=0.5 - ) - ax.add_patch(circ) - art3d.pathpatch_2d_to_3d(circ, z=0, zdir="z") - - # Source - # ax.scatter([0], [0], [0], color="g", s=100) - sourcePos3D = [geo.DSO * np.cos(angle), geo.DSO * np.sin(angle), 0] # xyz - ax.scatter([sourcePos3D[0]], [sourcePos3D[1]], [sourcePos3D[2]], color="steelblue", s=100) - - # Axes XYZ - length_axis = geo.sVoxel[0] - x_axis = Arrow3D( - [0, length_axis], - [0, 0], - [0, 0], - mutation_scale=10, - shrinkA=0, - lw=1, - arrowstyle="-|>", - color="r", - ) - y_axis = Arrow3D( - [0, 0], - [0, length_axis], - [0, 0], - mutation_scale=10, - shrinkA=0, - lw=1, - arrowstyle="-|>", - color="b", - ) - z_axis = Arrow3D( - [0, 0], - [0, 0], - [0, length_axis], - mutation_scale=10, - shrinkA=0, - lw=1, - arrowstyle="-|>", - color="g", - ) - ax.add_artist(x_axis) - ax.add_artist(y_axis) - ax.add_artist(z_axis) - - # Detector - print("sDetector = {}".format(geo.sDetector)) - alpha_detector = 0.7 - detectorPos3D = [ - (geo.DSD - geo.DSO) * np.cos(angle + np.pi) + geo.offDetector[1] * np.sin(angle + np.pi), - (geo.DSD - geo.DSO) * np.sin(angle + np.pi) - geo.offDetector[1] * np.cos(angle + np.pi), - geo.offDetector[0], - ] # xyz - detector = matplotlib.patches.Rectangle( - (-geo.sDetector[1] / 2, -geo.sDetector[0] / 2), - geo.sDetector[1], - geo.sDetector[0], - angle=0, # angle, - color="maroon", - fill=True, - alpha=alpha_detector, - ) - print("detector={}".format(detector)) - ax.add_patch(detector) - mat_rot = np.array( - [[-np.sin(angle), 0, np.cos(angle)], [np.cos(angle), 0, np.sin(angle)], [0, 1, 0]] - ) - pathpatch_2d_to_3d_affine(detector, mat_rot, detectorPos3D) - - # Image FOV - alpha_img = 0.1 - offOrigin = np.array([geo.offOrigin[2], geo.offOrigin[1], geo.offOrigin[0]]) - mat_rot_xy = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=np.float) - for idx in range(2): - img_face_xy = matplotlib.patches.Rectangle( - (-geo.sVoxel[2] / 2, -geo.sVoxel[1] / 2), # xy order - geo.sVoxel[2], - geo.sVoxel[1], - angle=0, # angle, - color="black", - fill=True, - alpha=alpha_img, - ) - ax.add_patch(img_face_xy) - pathpatch_2d_to_3d_affine( - img_face_xy, - mat_rot_xy, - np.array([0, 0, -geo.sVoxel[0] / 2 if idx == 0 else geo.sVoxel[0] / 2]) + offOrigin, - ) - - mat_rot_yz = np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]], dtype=np.float) - for idx in range(2): - img_face_yz = matplotlib.patches.Rectangle( - (-geo.sVoxel[1] / 2, -geo.sVoxel[0] / 2), # xy order - geo.sVoxel[1], - geo.sVoxel[0], - angle=0, # angle, - color="black", - fill=True, - alpha=alpha_img, - ) - ax.add_patch(img_face_yz) - pathpatch_2d_to_3d_affine( - img_face_yz, - mat_rot_yz, - np.array([-geo.sVoxel[2] / 2 if idx == 0 else geo.sVoxel[2] / 2, 0, 0]) + offOrigin, - ) - - mat_rot_zx = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]], dtype=np.float) - for idx in range(2): - img_face_zx = matplotlib.patches.Rectangle( - (-geo.sVoxel[0] / 2, -geo.sVoxel[2] / 2), # xy order - geo.sVoxel[0], - geo.sVoxel[2], - angle=0, # angle, - color="black", - fill=True, - alpha=alpha_img, - ) - ax.add_patch(img_face_zx) - pathpatch_2d_to_3d_affine( - img_face_zx, - mat_rot_zx, - np.array([0, -geo.sVoxel[1] / 2 if idx == 0 else geo.sVoxel[1] / 2, 0]) + offOrigin, - ) - plt.show() +import numpy as np +from matplotlib import pyplot as plt +from numpy import sin, cos +import mpl_toolkits.mplot3d.art3d as art3d +from matplotlib import animation +import tigre + +def plot_geometry(geo,angles=np.linspace(0,2*np.pi,100),pos=0,animate=False,fname=None): +#PLOT_GEOMETRY(GEO,ANGLES,POS,Animate,fname) plots a simplified version of the CBCT geometry with the +# given geomerty GEO and scanning angles ANGLES at angle POS. If angles is +# not given, [0,2pi] at 100 steps will be given, and pos=0 will be chosen. If +# Animate=TRUE, an animation file will be generated. Animation file name can be +# specified by fname. +# +# h=PLOT_GEOMETRY(...) will return the figure handle +#-------------------------------------------------------------------------- +#-------------------------------------------------------------------------- +# This file is part of the TIGRE Toolbox +# +# Copyright (c) 2015, University of Bath and +# CERN-European Organization for Nuclear Research +# All rights reserved. +# +# License: Open Source under BSD. +# See the full license at +# https://github.com/CERN/TIGRE/blob/master/LICENSE +# +# Contact: tigre.toolbox@gmail.com +# Codes: https://github.com/CERN/TIGRE/ +# Coded by: Ander Biguri, modified by Yi Liu +#-------------------------------------------------------------------------- + + # check geo, which makes DSD, DSO etc all matrices + try: + geo.check_geo(angles) + except: + print('Please check that geometry and angles are consistent') + + if pos < angles[0] or pos > angles[-1]: + raise ValueError('Pos should be within the range of [{}, {}]'.format(angles[0],angles[-1])) + + pos = abs(angles-pos).argmin() + + ## Figure stuff + fig=plt.figure(figsize=(10,8)) + ax=plt.axes(projection='3d') + ax.view_init(azim=52,elev=26) + + ln = geo.sVoxel.min() + + ## source trajectory, cordinates in (z,y,x) + thz,thx,thy = geo.angles[:,0], geo.angles[:,1], geo.angles[:,2] + Rs = rotmat3d(geo.angles,order='ZYZ') # source and detector rotation are in opposite direction to that of object + stj = np.zeros_like(geo.angles) + # Note: angles = (0,0,0) <==> source at (x=DSO, y=0, z=0) + scent = np.array([geo.DSO,0*thy,0*thz]).T # source centre (x,y,z) before scan rotation + for j in range(len(thx)): + stj[j,:] = np.matmul(Rs[:,:,j], scent[j,:]) # no offset for source + # displacement in y for geo.COR + if hasattr(geo, 'COR'): + stj[:,1] += geo.COR + if np.ptp(stj,axis=0).any(): + ax.plot3D(stj[:,0],stj[:,1],stj[:,2],color='grey',ls='',marker=".",markersize=2.5,mfc='grey',mec="grey") + # source centre at pos + source=ax.scatter(stj[pos,0],stj[pos,1],stj[pos,2],color='r',s=5) + stext=ax.text(stj[pos,0],stj[pos,1],stj[pos,2]+15,'S',None) + + ## detector trajectory: scan -> offset -> rotate + # detector centre (x,y,z) scan rotation. + dcent = np.array([-geo.DSD+geo.DSO, geo.offDetector[:,1], geo.offDetector[:,0]]).T + dtj = np.zeros_like(geo.angles) + for j in range(len(thx)): + dtj[j,:] = np.matmul(Rs[:,:,j], dcent[j,:]) + # displacement in y for geo.COR + if hasattr(geo, 'COR'): + dtj[:,1] += geo.COR + # # detector offset + if np.ptp(dtj,axis=0).any(): + ax.plot3D(dtj[:,0],dtj[:,1],dtj[:,2],color='grey',ls='',marker=".",markersize=2.5,mfc='grey',mec="grey") + # detector centre at pos + det=ax.scatter(dtj[pos,0],dtj[pos,1],dtj[pos,2],color='brown',s=5) + dtext=ax.text(dtj[pos,0],dtj[pos,1],dtj[pos,2]+15,'D',None) + # after scan to pos, then rotate the detector (roll,pitch,yaw) <==> (x,y,z) + R = np.zeros_like(Rs) + for j in range(Rs.shape[2]): + R[:,:,j] = np.matmul( Rs[:,:,j], rotmat3d(geo.rotDetector[j,:],order='XYZ')) + # Detector, cp returns four cordinates of corners closest to source + ddp = 5 # detector depth + dsz = np.array([ddp,geo.sDetector[1],geo.sDetector[0]]) # at angles (0,0,0) + dverts = calCube(dtj,dsz,R) + dcube = art3d.Poly3DCollection(dverts[pos],color='brown',alpha=0.3) + ax.add_collection3d(dcube) + + ## origin trajectory. NOTE: offOrigin is in (z,y,x) + otj = np.fliplr(geo.offOrigin).astype(np.float32) + # displacement in y for geo.COR + if hasattr(geo, 'COR'): + otj[:,1] += geo.COR + if np.ptp(otj,axis=0).any(): + ax.plot3D(otj[:,0],otj[:,1],otj[:,2],color='grey',ls='',marker=".",markersize=2.5,mfc='grey',mec="grey") + # Cordinates Arrows from origin + ax.quiver(0,0,0,1,0,0,length=ln,color='r') + ax.quiver(0,0,0,0,1,0,length=ln,color='b') + ax.quiver(0,0,0,0,0,1,length=ln,color='g') + ax.text(-10,-10,-10,'O',None) + # CUBE/Image, sVoxel dims order (z,y,x) in python + overts = calCube(otj,geo.sVoxel[[2,1,0]]) + ocube = art3d.Poly3DCollection(overts[pos],color='c',alpha=0.1) + ax.add_collection3d(ocube) + + # effective beam profile + beampf = [0,0,0,0] + for i in range(4): + beampf[i]=ax.plot3D(*zip(stj[pos,:],dverts[pos][0][i]),color='y') + # cenrtal beam + cbeam=ax.plot3D(*zip(stj[pos,:],dtj[pos,:]),color='pink') + + # set tight limits and aspect + limX = ax.get_xlim3d() + limY = ax.get_ylim3d() + limZ = ax.get_zlim3d() + ax.set_box_aspect((limX[1]-limX[0], limY[1]-limY[0], limZ[1]-limZ[0])) + + # set up plot + roll, pitch, yaw = geo.rotDetector[:,0]/np.pi*180, geo.rotDetector[:,1]/np.pi*180, geo.rotDetector[:,1]/np.pi*180 + plt.title('Current CBCT geometry at angles {} (deg), in scale'.format(list(map('{:.1f}'.format,geo.angles[pos]/np.pi*180))) \ + +'\n\nDetector rotation [{:.1f} - {:.1f}, {:.1f} - {:.1f}, {:.1f} - {:.1f}] (deg)'.format(roll[0],roll[-1],pitch[0],pitch[-1],yaw[0],yaw[-1])) + ax.set_xlabel('X'); + ax.set_ylabel('Y'); + ax.set_zlabel('Z'); + + ax.view_init(azim=52,elev=26) + plt.tight_layout() + + def update(pos,stj,dtj,otj,dverts,overts): + source.set_offsets(stj[pos,:2]) + source.set_3d_properties(stj[pos,2], 'z') + stext.set_position(stj[pos,:2]+15) + det.set_offsets(dtj[pos,:2]) + det.set_3d_properties(dtj[pos,2], 'z') + dtext.set_position(dtj[pos,:2]+15) + cbeam[0].set_data(*zip(stj[pos,:2],dtj[pos,:2])) + cbeam[0].set_3d_properties((stj[pos,2],dtj[pos,2]),'z') + for i in range(4): + beampf[i][0].set_data(*zip(stj[pos,:2],dverts[pos][0][i][:2])) + beampf[i][0].set_3d_properties((stj[pos,2],dverts[pos][0][i][2]),'z') + dcube.set_verts(dverts[pos]) + ocube.set_verts(overts[pos]) + ax.set_title('Current CBCT geometry at angles {} (deg), in scale'.format(list(map('{:.1f}'.format,geo.angles[pos]/np.pi*180))) \ + +'\n\nDetector rotation [{:.1f} - {:.1f}, {:.1f} - {:.1f}, {:.1f} - {:.1f}] (deg)'.format(roll[0],roll[-1],pitch[0],pitch[-1],yaw[0],yaw[-1])) + + + if animate: + ani = animation.FuncAnimation(fig, update, len(angles), fargs=(stj,dtj,otj,dverts,overts), interval=100) + if isinstance(fname, str): + fname += '_geometry' + try: + ani.save('%s.mp4'%(fname),writer='ffmpeg',fps=30) + except ValueError: + print('Movie writer "ffmpeg" unavailable, try "Pillow"') + try: + ani.save('%s.gif'%(fname),writer='pillow',fps=30) + except ValueError: + print('Movie writer "Pillow" unavailable, try "ImageMagick"') + try: + ani.save('%s.gif'%(fname),writer='imagemagick',fps=30) + except ValueError: + print('Movie writer unavailable, animation not saved.') + return ani + else: + return ax + +# other useful functions +def calCube(centre,size,R='eye'): + # Calculate vertices of a cuboid, centred at "centre" with "size" and rotated by R, 3D rotation matrix + # centre is either (n,3) or (3,) array and R is (3,3,n) or (3,3) array + ## 3D cube corners (x,y,z) + corner = np.array([[-1, -1, -1], + [1, -1, -1], + [1, 1, -1], + [-1, 1, -1], + [-1, -1, 1], + [1, -1, 1], + [1, 1, 1], + [-1, 1, 1]]) + + if centre.ndim>1: + n = centre.shape[0] + verts = [] + for i in range(n): + # scaling to cuboid and rotation + if isinstance(R, str): + c = np.matmul( corner*size/2, np.eye(3) ) + centre[i,:] + else: + c = np.matmul( corner*size/2, R[:,:,i].T ) + centre[i,:] + + # vertices of the cube + verts.append( [[c[1], c[2], c[6], c[5]], # front + [c[0], c[1], c[2], c[3]], # bottom + [c[4], c[5], c[6], c[7]], # top + [c[0], c[1], c[5], c[4]], # left + [c[2], c[3], c[7], c[6]], # right + [c[0], c[3], c[7], c[4]]] ) # back + else: + # scaling and rotation + if isinstance(R, str): + c = np.matmul( corner*size/2, np.eye(3) ) + centre + else: + c = np.matmul( corner*size/2, R.T ) + centre + + # vertices of the cube + verts = [[c[1], c[2], c[6], c[5]], # front + [c[0], c[1], c[2], c[3]], # bottom + [c[4], c[5], c[6], c[7]], # top + [c[0], c[1], c[5], c[4]], # left + [c[2], c[3], c[7], c[6]], # right + [c[0], c[3], c[7], c[4]]] # back + + return verts + + +def rotmat3d(ang,order='ZYX'): + # Calculate 3D rotation matrices of any Eular configuation (right-hand) + # ang: can be (3,) or (n,3) angles in rads + # order: str of any combination of "X", "Y", "Z", in any order + # cordinates in (x,y,z) + + if ang.ndim == 1: + an = np.expand_dims(ang,0) + else: + an = ang.copy() + + r = np.eye(3) + R = np.stack([r for _ in range(an.shape[0])], axis=2) + + for j in range(an.shape[0]): + for i in range(len(order)): + if order[i] in 'X': + R[:,:,j] = np.matmul(rotX(an[j,i]), R[:,:,j]) + elif order[i] in 'Y': + R[:,:,j] = np.matmul(rotY(an[j,i]), R[:,:,j]) + elif order[i] in 'Z': + R[:,:,j] = np.matmul(rotZ(an[j,i]), R[:,:,j]) + return np.squeeze(R) + + +def rotZ(a): + # rotation around z axis of a rads (yaw) + Rz = np.array([[cos(a),-sin(a),0*a], + [sin(a),cos(a),0*a], + [0*a,0*a,0*a+1]]) + return Rz + +def rotY(b): + # rotation around y axis of b rads (pitch) + Ry = np.array([[cos(b),0*b,sin(b)], + [0*b,0*b+1,0*b], + [-sin(b),0*b,cos(b)]]) + return Ry + +def rotX(c): + # rotation aroud x axis of c rads (roll) + Rx = np.array([[0*c+1,0*c,0*c], + [0*c,cos(c),-sin(c)], + [0*c,sin(c),cos(c)]]) + return Rx + if __name__ == "__main__": n_voxel_z = 128 n_voxel_y = 256 @@ -233,5 +289,7 @@ def plot_geometry(geo, angle=0): geo.offDetector = np.array([off_detector_v, off_detector_u]) geo.offOrigin = np.array([off_origin_z, off_origin_y, off_origin_x]) print(geo) - angle = -np.pi / 6 - plot_geometry(geo, angle) + angles = np.linspace(0.0,np.pi,180) + pos = np.pi / 6 + ani=plot_geometry(geo, angles, pos, animate=True) + plt.show()