import numpy as np
import matplotlib.pyplot as plt
[docs]
class ParallelepipedProjector:
def __init__(self, parametric_coordinates, pixelSizeXY=0.5, pixelSizeXYZ=0.5, crystal_width=2, crystal_height=2,
reflector_xy=0.35, reflector_z=0.28, FoV=45, crystal_depth=30, crystals_geometry=[32, 2],
bool_consider_reflector_in_z_projection=False, bool_consider_reflector_in_xy_projection=False,
distance_crystals=60, max_center_position=None):
FoV = FoV / pixelSizeXY
if bool_consider_reflector_in_z_projection:
crystal_height = crystal_height + reflector_z
else:
crystal_height = crystal_height
if bool_consider_reflector_in_xy_projection:
crystal_width = crystal_width + reflector_xy
else:
crystal_width = crystal_width
print("ParallelepipedProjector: Init")
self.x_min_f = None
self.x_max_f = None
self.y_min_f = None
self.y_max_f = None
self.z_min_f = None
self.z_max_f = None
self.pixelSizeXY = pixelSizeXY
self.pixelSizeXYZ = pixelSizeXYZ
self.half_crystal_pitch_xy = np.float32(0.5 * crystal_width / pixelSizeXY)
self.half_crystal_pitch_z = np.float32(0.5 * crystal_height / pixelSizeXYZ)
self.distance_between_array_pixel = np.float32((distance_crystals) / pixelSizeXY)
# self.distance_between_array_pixel = np.float32(40 / pixelSizeXYZ)
# self.crystal_planes = CrystalPlanes(parametric_coordinates_object=parametric_coordinates,
# pixelSizeXY=pixelSizeXY, pixelSizeXYZ=pixelSizeXYZ)
# self.crystal_planes.multiply_by_pixel_size()
# self.crystal_planes.crystals_central_planes()
# Convert coordinates to postive values inside te image
v1, v1_normal, v2, v3, v4, p1_list, p4_list, p5_list, p9_list = self.calculateVectors(parametric_coordinates, pixelSizeXY, pixelSizeXYZ)
self.calculateVolume(FoV, crystal_height, pixelSizeXYZ)
self.planeParameters(v1, v1_normal, v2, v3, v4, p1_list, p4_list, p5_list, p9_list)
self.number_of_events = len(self.a)
print("Number of events: {}".format(self.number_of_events))
[docs]
def calculateVolume(self, FoV, crystal_height, pixelSizeXYZ):
print("Calculate Volume")
self.number_of_pixels_x = int(np.ceil(FoV) + 1)
self.number_of_pixels_y = int(np.ceil(FoV) + 1)
self.number_of_pixels_z = int(np.ceil(self.max_z + crystal_height / pixelSizeXYZ))
print("Max_x fov: {}".format(self.max_x))
print("Max_y fov: {}".format(self.max_y))
print("Min_z fov: {}".format(self.min_z))
print("Max_z fov: {}".format(self.max_z))
# Create a int range array from 0 to number of pixels)
self.x_range_lim = [np.floor((self.max_x - FoV) / 2), np.ceil((self.max_x - FoV) / 2 + np.ceil(FoV))]
self.y_range_lim = [np.floor((self.max_y - FoV) / 2), np.ceil((self.max_y - FoV) / 2 + np.ceil(FoV))]
self.z_range_lim = [np.floor(self.min_z - 0.5 * crystal_height/pixelSizeXYZ), np.ceil(self.max_z + 0.5*crystal_height/pixelSizeXYZ)+1]
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)
self.number_of_pixels_x = int(self.x_range_lim[1] - self.x_range_lim[0])
self.number_of_pixels_y = int(self.y_range_lim[1] - self.y_range_lim[0])
self.number_of_pixels_z = int(self.z_range_lim[1] - self.z_range_lim[0])
# Create 3 EMPTY arrays with images size.
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))
# 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, ...]
print(self.im_index_x.shape)
[docs]
def calculateVectors(self, parametric_coordinates, pixelSizeXY, pixelSizeXYZ):
print("Calculate Vectors")
# Stack all coordinates in one array
# 4 points are needed to make two orthogonal planes.
# p1 and p2 are the crystals centers
# p3 is a point that belongs to the crystal with the same Z value as p1
# p4 is a point that belongs to the crystal it belongs to the normal plane of plane formed by p1, p2 and p3
xi = (parametric_coordinates.xi / pixelSizeXY).astype(np.float32)
xf = (parametric_coordinates.xf / pixelSizeXY).astype(np.float32)
yi = (parametric_coordinates.yi / pixelSizeXY).astype(np.float32)
yf = (parametric_coordinates.yf / pixelSizeXY).astype(np.float32)
zi = (parametric_coordinates.zi / pixelSizeXYZ).astype(np.float32)
zf = (parametric_coordinates.zf / pixelSizeXYZ).astype(np.float32)
self.max_x = np.array([xi.max(), xf.max()]).max()
self.max_y = np.array([yi.max(), yf.max()]).max()
self.max_z = np.array([zi.max(), zf.max()]).max()
self.min_z = np.array([zi.min(), zf.min()]).min()
p1_list = np.column_stack((xi, yi, zi))
self.p1_list = p1_list
p2_list = np.column_stack((xf, yf, zf))
p3_list = np.column_stack((((parametric_coordinates.midpoint[0]) / pixelSizeXY).astype(np.float32),
((parametric_coordinates.midpoint[1]) / pixelSizeXY).astype(np.float32),
((parametric_coordinates.midpoint[2]) / pixelSizeXY).astype(np.float32)))
p9_list = np.column_stack(((parametric_coordinates.farest_vertex[0] / pixelSizeXY).astype(np.float32),
(parametric_coordinates.farest_vertex[1] / pixelSizeXY).astype(np.float32),
(parametric_coordinates.farest_vertex[2] / pixelSizeXYZ).astype(np.float32)))
p4_list = p1_list.copy()
# p5_list = p1_list.copy()
p4_list[:, 2] = p4_list[:, 2] + self.half_crystal_pitch_z
p5_list = (p1_list + p2_list) / 2
p7_list = p5_list.copy()
p6_list = p5_list + p3_list - p1_list
p7_x = p7_list[:, 2]
p7_z = -p7_list[:, 0]
p7_list[:, 0] = p7_x
p7_list[:, 2] = p7_z
# p7_list[:, 2] = p7_list[:, 2] + self.half_crystal_pitch_z
self.x_min_f = np.ascontiguousarray(np.min(np.vstack([xi, xf]), axis=0) - 2, dtype=np.short)
self.z_min_f = np.ascontiguousarray(np.min(np.vstack([zi, zf]), axis=0) - 2, dtype=np.short)
self.z_max_f = np.ascontiguousarray(np.max(np.vstack([zi, zf]), axis=0) + 2, dtype=np.short)
v1, v1_normal, v2, v3, v4 = self.calculateNormalizedVectors(p1_list, p2_list, p3_list, p4_list, p5_list,
p6_list, p7_list)
return v1, v1_normal, v2, v3, v4, p1_list, p4_list, p5_list, p9_list
[docs]
def calculateNormalizedVectors(self, p1_list, p2_list, p3_list, p4_list, p5_list, p6_list, p7_list):
# Calculate the vectors
print("Calculate Normalized Vectors")
v1_normal = p4_list - p1_list
v1 = p2_list - p1_list
v2 = p3_list - p1_list
v3 = p5_list - p6_list
v4 = p7_list - p5_list
nf_v1 = ParallelepipedProjector.normVector(v1)
nf_v1_normal = ParallelepipedProjector.normVector(v1_normal)
nf_v2 = ParallelepipedProjector.normVector(v2)
nf_v3 = ParallelepipedProjector.normVector(v3)
nf_v4 = ParallelepipedProjector.normVector(v4)
for coor in range(3):
v1[:, coor] = v1[:, coor] / nf_v1
v1_normal[:, coor] = v1_normal[:, coor] / nf_v1_normal
v2[:, coor] = v2[:, coor] / nf_v2
v3[:, coor] = v3[:, coor] / nf_v3
v4[:, coor] = v4[:, coor] / nf_v4
return v1, v1_normal, v2, v3, v4
[docs]
def calcuteVectorsImagePlane(self):
p1 = np.array([self.x_range_lim[0], self.y_range_lim[0], self.z_range_lim[0]])
p2 = np.array([self.x_range_lim[1], self.y_range_lim[0], self.z_range_lim[0]])
p3 = np.array([self.x_range_lim[0], self.y_range_lim[1], self.z_range_lim[0]])
p4 = np.array([self.x_range_lim[1], self.y_range_lim[1], self.z_range_lim[0]])
p5 = np.array([self.x_range_lim[0], self.y_range_lim[0], self.z_range_lim[1]])
p6 = np.array([self.x_range_lim[1], self.y_range_lim[0], self.z_range_lim[1]])
p7 = np.array([self.x_range_lim[0], self.y_range_lim[1], self.z_range_lim[1]])
p8 = np.array([self.x_range_lim[1], self.y_range_lim[1], self.z_range_lim[1]])
points_face_image = np.array([p1, p2, p3, p4, p5, p6, p7, p8])
vector_image_xx = p2 - p1
vector_image_yy = p3 - p1
vector_image_zz = p5 - p1
vector_image_xx_zmax = p6 - p8
vector_image_yy_zmax = p7 - p8
vector_image_zz_zmax = p4 - p8
nf_vector_image_xx = ParallelepipedProjector.normVector(vector_image_xx)
nf_vector_image_yy = ParallelepipedProjector.normVector(vector_image_yy)
nf_vector_image_zz = ParallelepipedProjector.normVector(vector_image_zz)
nf_vector_image_xx_zmax = ParallelepipedProjector.normVector(vector_image_xx_zmax)
nf_vector_image_yy_zmax = ParallelepipedProjector.normVector(vector_image_yy_zmax)
nf_vector_image_zz_zmax = ParallelepipedProjector.normVector(vector_image_zz_zmax)
for coor in range(3):
vector_image_xx[coor] /= nf_vector_image_xx
vector_image_yy[coor] /= nf_vector_image_yy
vector_image_zz[coor] /= nf_vector_image_zz
vector_image_xx_zmax[coor] /= nf_vector_image_xx_zmax
vector_image_yy_zmax[coor] /= nf_vector_image_yy_zmax
vector_image_zz_zmax[coor] /= nf_vector_image_zz_zmax
return vector_image_xx, vector_image_yy, vector_image_zz, vector_image_xx_zmax, vector_image_yy_zmax, vector_image_zz_zmax, points_face_image
[docs]
def planesImage(self):
vector_image_xx, vector_image_yy, vector_image_zz, vector_image_xx_zmax, vector_image_yy_zmax, vector_image_zz_zmax, points_face_image = self.calcuteVectorsImagePlane()
face_image_A_plane = ParallelepipedProjector.planeEquation(vector_image_xx, vector_image_yy, points_face_image[0])
face_image_B_plane = ParallelepipedProjector.planeEquation(vector_image_xx, vector_image_zz, points_face_image[0])
face_image_C_plane = ParallelepipedProjector.planeEquation(vector_image_yy, vector_image_zz, points_face_image[0])
face_image_D_plane = ParallelepipedProjector.planeEquation(vector_image_xx_zmax, vector_image_yy_zmax, points_face_image[7])
face_image_E_plane = ParallelepipedProjector.planeEquation(vector_image_xx_zmax, vector_image_zz_zmax, points_face_image[7])
face_image_F_plane = ParallelepipedProjector.planeEquation(vector_image_yy_zmax, vector_image_zz_zmax, points_face_image[7])
planes_face_image = np.array([face_image_A_plane, face_image_B_plane, face_image_C_plane, face_image_D_plane, face_image_E_plane, face_image_F_plane])
return planes_face_image
[docs]
@staticmethod
def normVector(v1):
try:
return (np.sqrt(v1[:, 0] ** 2 + v1[:, 1] ** 2 + v1[:, 2] ** 2)).astype(np.float32)
except IndexError:
return (np.sqrt(v1[0] ** 2 + v1[1] ** 2 + v1[2] ** 2)).astype(np.float32)
[docs]
def planeParameters(self, v1, v1_normal, v2, v3, v4, p1_list, p4_list, p5_list, p9_list):
print("Calculate Plane Parameters")
self.v1 = v1
self.min_max_position_of_the_event_in_image()
# These operations gives the equations of planes
self.a, self.b, self.c, self.d = ParallelepipedProjector.planeEquation(v1, v2, p1_list)
# self.distance_to_central_plane = np.ascontiguousarray(np.abs(
# self.a * p9_list[:, 0] + self.b * p9_list[:, 1] + self.c * p9_list[:, 2] - self.d) / np.sqrt(
# self.a ** 2 + self.b ** 2 + self.c ** 2) + self.half_crystal_pitch_z, dtype=np.float32)
self.distance_to_central_plane = np.ascontiguousarray(np.abs(
self.a * p9_list[:, 0] + self.b * p9_list[:, 1] + self.c * p9_list[:, 2] - self.d) / np.sqrt(
self.a ** 2 + self.b ** 2 + self.c ** 2), dtype=np.float32)
self.a_normal, self.b_normal, self.c_normal, self.d_normal = ParallelepipedProjector.planeEquation(v1_normal,
v1,
p4_list)
# self.distance_to_central_plane_normal = np.ascontiguousarray(np.abs(
# self.a_normal * p9_list[:, 0] + self.b_normal * p9_list[:, 1] + self.c_normal * p9_list[:,
# 2] - self.d_normal) / np.sqrt(
# self.a_normal ** 2 + self.b_normal ** 2 + self.c_normal ** 2) + self.half_crystal_pitch_xy,
# dtype=np.float32)
self.distance_to_central_plane_normal = np.ascontiguousarray(np.abs(
self.a_normal * p9_list[:, 0] + self.b_normal * p9_list[:, 1] + self.c_normal * p9_list[:,
2] - self.d_normal) / np.sqrt(
self.a_normal ** 2 + self.b_normal ** 2 + self.c_normal ** 2),
dtype=np.float32)
self.a_cf, self.b_cf, self.c_cf, self.d_cf = ParallelepipedProjector.planeEquation(v4, v3, p5_list)
# self.distance_to_central_plane_cf = np.ascontiguousarray(np.abs(
# self.a_cf * p9_list[:, 0] + self.b_cf * p9_list[:, 1] + self.c_cf * p9_list[:,
# 2] - self.d_cf) / np.sqrt(
# self.a_cf ** 2 + self.b_cf ** 2 + self.c_cf ** 2), dtype=np.float32)
self.distance_to_central_plane_cf = np.ascontiguousarray(np.abs(
self.a_cf * p9_list[:, 0] + self.b_cf * p9_list[:, 1] + self.c_cf * p9_list[:,
2] - self.d_cf) / np.sqrt(
self.a_cf ** 2 + self.b_cf ** 2 + self.c_cf ** 2), dtype=np.float32)
print("END- Equation Planes")
[docs]
@staticmethod
def planeEquation(v1, v2, p1):
try:
cp = np.cross(v2, v1).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)
d = np.ascontiguousarray(np.array(cp[:, 0] * p1[:, 0] + cp[:, 1] * p1[:, 1] + cp[:, 2] * p1[:, 2],
dtype=np.float32))
except IndexError:
cp = np.cross(v2, v1).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)
d = np.ascontiguousarray(np.array(cp[0] * p1[0] + cp[1] * p1[1] + cp[2] * p1[2], dtype=np.float32))
return a, b, c, d
[docs]
def cut_current_frame(self, frame_start, frame_end):
self.distance_to_central_plane = self.distance_to_central_plane[frame_start:frame_end]
self.distance_to_central_plane_normal = self.distance_to_central_plane_normal[frame_start:frame_end]
self.distance_to_central_plane_cf = self.distance_to_central_plane_cf[frame_start:frame_end]
# self.p1_list = self.p1_list[frame_start:frame_end]
# self.p2_list = self.p2_list[frame_start:frame_end]
# self.p3_list = self.p3_list[frame_start:frame_end]
# self.p4_list = self.p4_list[frame_start:frame_end]
self.a = self.a[frame_start:frame_end]
self.b = self.b[frame_start:frame_end]
self.c = self.c[frame_start:frame_end]
self.d = self.d[frame_start:frame_end]
self.a_normal = self.a_normal[frame_start:frame_end]
self.b_normal = self.b_normal[frame_start:frame_end]
self.c_normal = self.c_normal[frame_start:frame_end]
self.d_normal = self.d_normal[frame_start:frame_end]
self.a_cf = self.a_cf[frame_start:frame_end]
self.b_cf = self.b_cf[frame_start:frame_end]
self.c_cf = self.c_cf[frame_start:frame_end]
self.d_cf = self.d_cf[frame_start:frame_end]
# self.p1_list = self.p1_list[frame_start:frame_end]
# self.p2_list = self.p2_list[frame_start:frame_end]
# sle.p3_list = self.p3_list[frame_start:frame_end]
# self.p4_list = self.p4_list[frame_start:frame_end]
# self.p5_list = self.p5_list[frame_start:frame_end]
# self.p6_list = self.p6_list[frame_start:frame_end]
# self.p7_list = self.p7_list[frame_start:frame_end]
# self.p9_list = self.p9_list[frame_start:frame_end]
[docs]
def min_max_position_of_the_event_in_image(self):
# calculate the intersection of the coincidence vector v1 with image planes
planes_face_image = self.planesImage()
# calculate the determinant
i = 0
intersections = np.zeros((planes_face_image.shape[0],self.p1_list.shape[0], 3), dtype=np.float32)
for plane in planes_face_image:
# t = -(plane[3]+plane[0]*self.p1_list[:,0]+plane[1]*self.p1_list[:,1]+plane[2]*self.p1_list[:,2])/(plane[0]*plane[0]+plane[1]*plane[1]+plane[2]*plane[2])
t = -(-plane[3] + self.p1_list[:, 0] * plane[0] + self.p1_list[:, 1] * plane[1] + self.p1_list[:, 2] * plane[2])/(plane[0]*self.v1[:,0]+plane[1]*self.v1[:,1]+plane[2]*self.v1[:,2])
intersections[i] = self.p1_list + self.v1 * t[:, np.newaxis]
# for each plane set intersection the value of np.nan if the intersection is out of the image in the 3 directions, x_rangelim, y_rangelim, z_rangelim
mask_x = (intersections[i, :, 0] < self.x_range_lim[0]) | (intersections[i, :, 0] > self.x_range_lim[1])
mask_y = (intersections[i, :, 1] < self.y_range_lim[0]) | (intersections[i, :, 1] > self.y_range_lim[1])
mask_z = (intersections[i, :, 2] < self.z_range_lim[0]) | (intersections[i, :, 2] > self.z_range_lim[1])
mask = mask_x | mask_y | mask_z
intersections[i, mask] = np.nan
# t = -(plane[3] + np.dot(plane[:3], self.v1)) / np.dot(plane[:3], plane[:3])
# intersections[i] = self.v1.T + plane[:3] * t
i += 1
# if the value is bigger than range in each direction of the image, the intersection is set no np.nan
# intersections[np.where((intersections[:, :, 0] <= self.x_range_lim[0]) | (intersections[:, :, 0] >= self.x_range_lim[1]))] = np.nan
# intersections[np.where((intersections[:, :, 1] <= self.y_range_lim[0]) | (intersections[:, :, 1] >= self.y_range_lim[1]))] = np.nan
# intersections[np.where((intersections[:, :, 2] <= self.z_range_lim[0]) | (intersections[:, :, 2] >= self.z_range_lim[1]))] = np.nan
# calculate max and min of the intersection not nan
x_min = np.nanmin(intersections[:, :, 0], axis=0)
x_max = np.nanmax(intersections[:, :, 0], axis=0)
y_min = np.nanmin(intersections[:, :, 1], axis=0)
y_max = np.nanmax(intersections[:, :, 1], axis=0)
z_min = np.nanmin(intersections[:, :, 2], axis=0)
z_max = np.nanmax(intersections[:, :, 2], axis=0)
x_min[np.argwhere(np.isnan(x_min))] = self.x_range_lim[0]
x_max[np.argwhere(np.isnan(x_max))] = self.x_range_lim[1]
y_min[np.argwhere(np.isnan(y_min))] = self.y_range_lim[0]
y_max[np.argwhere(np.isnan(y_max))] = self.y_range_lim[1]
z_min[np.argwhere(np.isnan(z_min))] = self.z_range_lim[0]
z_max[np.argwhere(np.isnan(z_max))] = self.z_range_lim[1]
self.x_min_f = np.ascontiguousarray(x_min, dtype=np.short)
self.x_max_f = np.ascontiguousarray(x_max, dtype=np.short)
self.y_min_f = np.ascontiguousarray(y_min, dtype=np.short)
self.y_max_f = np.ascontiguousarray(y_max, dtype=np.short)
self.z_min_f = np.ascontiguousarray(z_min, dtype=np.short)
self.z_max_f = np.ascontiguousarray(z_max, dtype=np.short)
[docs]
def plots(self):
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
colors = ['red', 'red', 'green', 'blue', 'yellow', 'black', 'brown']
points = [self.p1_list, self.p2_list, self.p3_list, self.p4_list, self.p5_list, self.p6_list, self.p7_list]
for point, color in zip(points, colors):
ax.scatter(point[:, 0], point[:, 1], point[:, 2], color=color)
plt.show()