Source code for toor.Corrections.CT.Projector.pyramidalprojector

import numpy as np
# from Projector import GeneralProjector


[docs] class PyramidalProjector: def __init__(self, voxelSize=None, FoVAxial=45, FoVRadial=25, FoVTangencial=30, FovRadialStart=None, FovRadialEnd=None, fov=45, only_fov=False): """ Pyramidal Projector :param voxelSize: list with the voxel dimensions in mm :param FoVAxial: Field of view in axial direction in mm :param FoVRadial: Field of view in radial direction in mm :param FoVTangencial: Field of view in tangencial direction in mm """ if voxelSize is None: voxelSize = [0.25, 0.25, 0.25] # super().__init__(voxelSize=voxelSize,FoVAxial=45, FoVRadial=25, FoVTangencial=30, FovRadialStart=FovRadialStart, # FovRadialEnd=FovRadialEnd) self.type = "Pyramidal" self.voxelSize = voxelSize self.FoVAxial = FoVAxial self.FoVRadial = FoVRadial/voxelSize[0] self.FoVRadial_min = FovRadialStart/voxelSize[0] self.FoVRadial_max = FovRadialEnd/voxelSize[0] self.FoVTangencial = FoVTangencial self._only_fov = only_fov self.fov = fov / voxelSize[0] self.pointCenterList = None self.pointCorner1List = None self.pointCorner2List = None self.pointCorner3List = None self.pointCorner4List = None self.planes = None # Plane Left self.aLeft = None self.bLeft = None self.cLeft = None self.dLeft = None # Plane Right self.aRight = None self.bRight = None self.cRight = None self.dRight = None # Plane Front self.aFront = None self.bFront = None self.cFront = None self.dFront = None # Plane Back self.aBack = None self.bBack = None self.cBack = None self.dBack = None self.number_of_pixels_x = int(np.ceil(self.FoVRadial)) self.number_of_pixels_y = int(np.ceil(self.FoVTangencial)) self.number_of_pixels_z = int(np.ceil(self.FoVAxial)) self._absMinMin = np.zeros((3,), dtype=np.float32) self.x_range_lim = None self.y_range_lim = None self.z_range_lim = None self._countsPerPosition = None @property def countsPerPosition(self): """ Return the number of counts per position Valid for non ListMode data :return: """ return self._countsPerPosition
[docs] def setCountsPerPosition(self, countsPerPosition): """ Set the number of counts per position Valid for non ListMode data :param countsPerPosition: :return: """ self._countsPerPosition = countsPerPosition self._countsPerPosition = np.ascontiguousarray(self._countsPerPosition, dtype=np.int32)
[docs] def setAbsMinMin(self): self._absMinMin = np.zeros((3,), dtype=np.float32) for coor in range(3): self._absMinMin[coor] = np.abs(np.min([np.min(self.pointCenterList[:, coor]),np.min(self.pointCorner1List[:, coor]), np.min(self.pointCorner2List[:, coor]), np.min(self.pointCorner3List[:, coor]), np.min(self.pointCorner4List[:, coor])]))
[docs] def transformIntoPositivePoints(self): """ Transform the points into positive points """ for coor in range(3): # abs_min_min = np.abs(np.min([np.min(self.pointCenterList[:, coor]),np.min(self.pointCorner1List[:, coor]), np.min(self.pointCorner2List[:, coor]), # np.min(self.pointCorner3List[:, coor]), np.min(self.pointCorner4List[:, coor])])) abs_min_min = self._absMinMin[coor] # abs_min_min = np.abs(np.min(self.pointCenterList[:, coor])) self.pointCorner1List[:, coor] += abs_min_min self.pointCorner2List[:, coor] += abs_min_min self.pointCorner3List[:, coor] += abs_min_min self.pointCorner4List[:, coor] += abs_min_min self.pointCenterList[:, coor] += abs_min_min
[docs] def amplifyPointsToGPUCoordinateSystem(self): """ Amplify the points to the GPU coordinate system acorddingly to the voxel size :return: """ for coor in range(3): self.pointCenterList[:, coor] = (self.pointCenterList[:, coor]/self.voxelSize[coor]).astype(np.float32) self.pointCorner1List[:, coor] = (self.pointCorner1List[:, coor]/self.voxelSize[coor]).astype(np.float32) self.pointCorner2List[:, coor] = (self.pointCorner2List[:, coor]/self.voxelSize[coor]).astype(np.float32) self.pointCorner3List[:, coor] = (self.pointCorner3List[:, coor]/self.voxelSize[coor]).astype(np.float32) self.pointCorner4List[:, coor] = (self.pointCorner4List[:, coor]/self.voxelSize[coor]).astype(np.float32)
[docs] def createVectorialSpace(self): """ Create the vectorial space Needs pointCenterList, pointCorner1List, pointCorner2List, pointCorner3List, pointCorner4List to be defined first """ print("Creating vectorial space") print("Transforming points to positive points") # Create a int range array from 0 to number of pixels) # self.x_range_lim = [np.floor((np.max(x) - FoV) / 2), np.ceil((np.max(x) - FoV) / 2 + np.ceil(FoV))] # self.y_range_lim = [np.floor((np.max(y) - FoV) / 2), np.ceil((np.max(y) - FoV) / 2 + np.ceil(FoV))] # self.z_range_lim = [np.floor(np.min(z) - 1.5 * crystal_height), np.ceil(np.max(z) + 1.5 * crystal_height)] # self.x_range_lim = [-self.FoVRadial/2, self.FoVRadial/2] # self.y_range_lim = [-self.FoVTangencial/2, self.FoVTangencial/2] # self.z_range_lim = [-self.FoVAxial/2, self.FoVAxial/2] self.max_x = self.pointCorner1List[:,0].max() self.max_y = self.pointCorner1List[:,1].max() self.min_x = self.pointCorner1List[:,0].min() self.min_y = self.pointCorner1List[:,1].min() print(f"Max X: {self.max_x} Min X: {self.min_x}") print(f"Max Y: {self.max_y} Min Y: {self.min_y}") extra_pixel_x = np.abs(self.pointCorner1List[:,0].min() - self.pointCorner2List[:,0].min()) extra_pixel_y = np.abs(self.pointCorner1List[:,1].min() - self.pointCorner4List[:,1].min()) extra_pixel_z = np.abs(self.pointCorner1List[:,2].min() - self.pointCorner3List[:,2].min()) if self._only_fov: self.x_range_lim = [np.floor((self.max_x - self.fov) / 2), np.ceil((self.max_x - self.fov) / 2 + np.ceil(self.fov))] self.y_range_lim = [np.floor((self.max_y - self.fov) / 2), np.ceil((self.max_y - self.fov) / 2 + np.ceil(self.fov))] # self.x_range_lim = [0, 60] # self.y_range_lim = [0, 60] else: self.x_range_lim = [np.ceil(self.pointCorner1List[:,0].min()), np.ceil(self.pointCorner1List[:,0].max())] self.y_range_lim = [np.ceil(self.pointCorner1List[:,1].min()), np.ceil(self.pointCorner1List[:,1].max())] self.z_range_lim = [np.floor(self.pointCorner3List[:,2].min()), np.ceil(self.pointCorner1List[:,2].max())] # self.z_range_lim = [0, 73/self.voxelSize[2]] self.number_of_pixels_x = int(np.ceil(self.x_range_lim[1]-self.x_range_lim[0])) self.number_of_pixels_y = int(np.ceil(self.y_range_lim[1]-self.y_range_lim[0])) self.number_of_pixels_z = int(np.ceil(self.z_range_lim[1]-self.z_range_lim[0])) self.im_index_x = np.ascontiguousarray( np.empty((self.number_of_pixels_x, self.number_of_pixels_y, self.number_of_pixels_z), dtype=np.int32)) self.im_index_y = np.ascontiguousarray( np.empty((self.number_of_pixels_x, self.number_of_pixels_y, self.number_of_pixels_z), dtype=np.int32)) self.im_index_z = np.ascontiguousarray( np.empty((self.number_of_pixels_x, self.number_of_pixels_y, self.number_of_pixels_z), dtype=np.int32)) print(f"Image Shape{self.im_index_z.shape}") x_range = np.arange(self.x_range_lim[0], self.x_range_lim[1], dtype=np.int32) y_range = np.arange(self.y_range_lim[0], self.y_range_lim[1], dtype=np.int32) z_range = np.arange(self.z_range_lim[0], self.z_range_lim[1], dtype=np.int32) # Create 3 EMPTY arrays with images size. # Repeat values in one direction. Like x_axis only grows in x_axis (0, 1, 2 ... number of pixels) # but repeat these values on y an z axis self.im_index_x[:] = x_range[..., None, None] self.im_index_y[:] = y_range[None, ..., None] self.im_index_z[:] = z_range[None, None, ...]
[docs] def createPlanes(self): """ Create the planes of the projector Defines the planes equations ax + by + cz + d = 0 for each plane. The planes are defined by the 4 corners of the collimator and the focal point :return: """ v1 = PyramidalProjector.calcVector(self.pointCenterList, self.pointCorner1List) v2 = PyramidalProjector.calcVector(self.pointCenterList, self.pointCorner2List) v3 = PyramidalProjector.calcVector(self.pointCenterList, self.pointCorner3List) v4 = PyramidalProjector.calcVector(self.pointCenterList, self.pointCorner4List) n1 = PyramidalProjector.norm(v1) n2 = PyramidalProjector.norm(v2) n3 = PyramidalProjector.norm(v3) n4 = PyramidalProjector.norm(v4) for coor in range(3): v1[:, coor] = v1[:, coor] / n1 v2[:, coor] = v2[:, coor] / n2 v3[:, coor] = v3[:, coor] / n3 v4[:, coor] = v4[:, coor] / n4 self.aLeft, self.bLeft, self.cLeft, cp = PyramidalProjector.crossProduct(v1, v2) self.dLeft = PyramidalProjector.calcD(self.pointCenterList, cp) self.aFront, self.bFront, self.cFront, cp = PyramidalProjector.crossProduct(v2, v3) self.dFront = PyramidalProjector.calcD(self.pointCenterList, cp) self.aRight, self.bRight, self.cRight, cp = PyramidalProjector.crossProduct(v3, v4) self.dRight = PyramidalProjector.calcD(self.pointCenterList, cp) self.aBack, self.bBack, self.cBack, cp = PyramidalProjector.crossProduct(v4, v1) self.dBack = PyramidalProjector.calcD(self.pointCenterList, cp) self.planes = np.array([[self.aLeft, self.bLeft, self.cLeft, self.dLeft], [self.aRight, self.bRight, self.cRight, self.dRight], [self.aFront, self.bFront, self.cFront, self.dFront], [self.aBack, self.bBack, self.cBack, self.dBack]])
[docs] @staticmethod def norm(vector): """ Calculate the norm of a vector :param vector: :return: """ return np.sqrt(vector[:, 0] ** 2 + vector[:, 1] ** 2 + vector[:, 2] ** 2)
[docs] @staticmethod def calcVector(p1, p2): """ Calculate the vector between two points :param p1: :param p2: :return: """ return p1-p2
[docs] @staticmethod def calcD(point, crossproduct): """ Calculate the d value of the plane equation :param point: :param crossproduct: :return: """ cp = crossproduct[:, 0] * point[:, 0] + crossproduct[:, 1] * point[:, 1]+crossproduct[:, 2] * point[:, 2] cp = np.ascontiguousarray(np.array(cp, dtype=np.float32)) return cp
[docs] @staticmethod def crossProduct(v1,v2): """ Calculate the cross product of two vectors :param v1: :param v2: :return: """ cp = np.cross(v1, v2).astype(np.float32) a = np.ascontiguousarray(cp[:, 0], dtype=np.float32) b = np.ascontiguousarray(cp[:, 1], dtype=np.float32) c = np.ascontiguousarray(cp[:, 2], dtype=np.float32) return a, b, c, cp
if __name__ == "__main__": from Geometry import HeadGeometry spectGeometry = HeadGeometry() spectGeometry.calculateInitialGeometry() pointCenter = np.copy(spectGeometry.CZTModules[0].initialMatrix).T p1 = spectGeometry.collimators[0].vertex1.T p2 = spectGeometry.collimators[0].vertex2.T p3 = spectGeometry.collimators[0].vertex3.T p4 = spectGeometry.collimators[0].vertex4.T p = PyramidalProjector(FovRadialStart=0, FovRadialEnd=20, FoVRadial=20) p.pointCenterList = pointCenter p.pointCorner1List = p1 p.pointCorner2List = p2 p.pointCorner3List = p3 p.pointCorner4List = p4 p.createVectorialSpace() p.createPlanes() # gd = GeometryDesignerObject() # for i in range(spectGeometry.numberOfModules): # Czt1 = spectGeometry.CZTModules[i] # # gd.drawDetectors(geometryVector=[Czt1.initialMatrix[0], Czt1.initialMatrix[1], Czt1.initialMatrix[2]], # moduleNumber=0) # p = [ # pointCenter, # p1, # p2, # p3, # p4 # ] # gd.makePyramid(p) # gd.makeAxes() # # # gd.drawPlanes(plane_values=[p.aLeft, p.bLeft, p.cLeft]) # # gd.drawPlanes(plane_values=[p.aRight, p.bRight, p.cRight]) # gd.renderWin.Render() # gd.renderInteractor.Start()