"""Tools to manage array square cell (index <=> position)
"""
import logging
import numpy as np
import matplotlib.pyplot as plt
from numba import njit
from ecpi.common.instru.array_convention import yzegp_to_ijdet_array
#pylint:disable=C0326
[docs]@njit
def intersection_rect(p_r1, p_r2): # pragma: no cover
"""Compute the rectangle formed by the
intersection of two rectangles.
Each rectangle is formed by an array-like
structure of the form [[ll_x, ll_y], [ur_x, ur_y]]
where the lower left (ll) corner has coordinates (ll_x, ll_y)
and the upper right (ur) corner has coordinates (ur_x, ur_y).
If the two rectangles do not overlap then
the function returns an array filled with zeros.
:param p_r1: first rectangle
:type p_r1: array-like structured as [[ll_x, ll_y], [ur_x, ur_y]]
:param p_r2: second rectangle
:type p_r2: array-like structured as [[ll_x, ll_y], [ur_x, ur_y]]
:return: intersection of the two input rectangles
:rtype: array-like structured as [[ll_x, ll_y], [ur_x, ur_y]]
"""
rec_inter_0 = np.zeros((2, 2), dtype=np.float32)
rec_inter = np.empty_like(rec_inter_0)
corner_lb_p = p_r1[0,:]
corner_ru_p = p_r1[1,:]
corner_lb_pp = p_r2[0,:]
corner_ru_pp = p_r2[1,:]
# corner lb intersection
if corner_lb_pp[0] < corner_lb_p[0]:
rec_inter[0, 0] = corner_lb_p[0]
elif corner_lb_pp[0] < corner_ru_p[0]:
rec_inter[0, 0] = corner_lb_pp[0]
else:
return rec_inter_0
if corner_lb_pp[1] < corner_lb_p[1]:
rec_inter[0, 1] = corner_lb_p[1]
elif corner_lb_pp[1] < corner_ru_p[1]:
rec_inter[0, 1] = corner_lb_pp[1]
else:
return rec_inter_0
if corner_ru_pp[0] < corner_lb_p[0]:
return rec_inter_0
elif corner_ru_pp[0] < corner_ru_p[0]:
rec_inter[1, 0] = corner_ru_pp[0]
else:
rec_inter[1, 0] = corner_ru_p[0]
if corner_ru_pp[1] < corner_lb_p[1]:
return rec_inter_0
elif corner_ru_pp[1] < corner_ru_p[1]:
rec_inter[1, 1] = corner_ru_pp[1]
else:
rec_inter[1, 1] = corner_ru_p[1]
return rec_inter
[docs]def surf_intersection_rect(p_r1, p_r2):
"""
:param p_r1: first rectangle
:type p_r1: array-like structured as [[ll_x, ll_y], [ur_x, ur_y]]
:param p_r2: second rectangle
:type p_r2: array-like structured as [[ll_x, ll_y], [ur_x, ur_y]]
:return: intersection of the two input rectangles
"""
inter_rec = intersection_rect(p_r1, p_r2)
surface = (inter_rec[1, 0] - inter_rec[0, 0]) * (inter_rec[1, 1] - inter_rec[0, 1])
return surface
[docs]class ArraySquareCell():
"""Features between index and physical position for a square cell array
In cartesian convention, ie index (0,0) is the corner left bottom of array.
::
Example :
- ArraySquareCell[1, 2] = 1
- ArraySquareCell[0, 0] = 5
Second axis
^
|
|
3 +-+-+
|2|1|
2 +-+-+
|3|4|
1 +-+-+
|5|6|
0 +-+-+-+----> First axis
0 1 2 3
"""
def __init__(self, ncel_x=1, ncel_y=1, size_cm=1.0):
"""*Constructor*
:param ncel_x: number of cells along x axis.
:type ncel_x: int>0. Default is 1.
:param ncel_y: number of cells along y axis.
:type ncel_y: int>0. Default is 1.
:param size_cm: area size of the square array in cm2.
:type size_cm: float>0. DEfault is 1.0
"""
self.max_ix = ncel_x - 1
self.max_iy = ncel_y - 1
self.size_cm = size_cm * 1.0
self.surf_cell = self.size_cm ** 2
self.ar_val = np.zeros((ncel_x, ncel_y), dtype=np.float32)
# vector to pass to origin index to origin referential (in cm)
self.to_origin_cm = np.array([0.0, 0.0], dtype=np.float32)
self.l_verb = 0
# declare logger
self.logger = logging.getLogger(__name__)
[docs] def zeros_like(self, p_asc):
"""Initiation of current ArraySquareCell object
with parameter of ArraySquareCell p_asc.
:param p_asc: array model
:type p_asc: ArraySquareCell
"""
assert isinstance(p_asc, ArraySquareCell)
self.max_ix = p_asc.max_ix
self.max_iy = p_asc.max_iy
self.size_cm = p_asc.size_cm
self.surf_cell = self.size_cm ** 2
self.ar_val = np.zeros((self.max_ix + 1, self.max_iy + 1), dtype=np.float32)
self.to_origin_cm = p_asc.to_origin_cm
[docs] def cast_array(self, p_type):
'''Cast array with p_type
:param p_type: type number
'''
self.ar_val = self.ar_val.astype(p_type)
[docs] def set_array(self, p_ar):
"""Set the current ArraySquareCell object
number of cells with those the p_ar ndarray.
:param p_ar: array model
:type p_ar: array
"""
assert isinstance(p_ar, np.ndarray)
if (self.get_nb_cell_y() != p_ar.shape[1]):
return None
if (self.get_nb_cell_x() != p_ar.shape[0]):
return None
self.ar_val[:,:] = p_ar[:,:]
[docs] def add_hit(self, apos):
"""Add one on pixel associated to position
:param apos: array plan position of photon
:type apos: numpy array (n,2)
"""
if apos.size == 0:
return
idx_pix = self.pos2idx_inf(apos)
it = np.nditer([idx_pix[:,0], idx_pix[:,1]], [], [['readonly'], ['readonly']])
for il, ic in it:
self.ar_val[il, ic] += 1
[docs] def set_origin_center(self):
'''Put the origin of the positions in the center of the array
'''
self.to_origin_cm = self.get_sizecm() / 2
[docs] def set_pos_corner_lb(self, p_pos_ori: np.array):
"""
vector OP where P=pos_array_lb and O: origin
lb = left bottom
:param p_pos_ori: position of the origin
:type p_pos_ori: array(float)
"""
self.to_origin_cm = -p_pos_ori.copy()
[docs] def set_to_origin(self, p_to_ori: np.array):
"""
vector PO where P=pos_array_lb and O: origin
:param p_to_ori: vector lb to origin
:type p_to_ori: array(float)
"""
self.to_origin_cm = p_to_ori.copy()
[docs] def get_nb_cell_x(self):
'''
Return number of cell in first axis
'''
return self.ar_val.shape[0]
[docs] def get_nb_cell_y(self):
'''
Return number of cell in second axis
'''
return self.ar_val.shape[1]
[docs] def get_nb_cell(self):
'''
return total cell in array
'''
return self.ar_val.size
[docs] def get_sizecm(self):
'''
Return size array in cm
'''
array_size = np.array([0, 0], dtype=np.float32)
array_size[0] = np.float32(self.ar_val.shape[0]) * self.size_cm
array_size[1] = np.float32(self.ar_val.shape[1]) * self.size_cm
return array_size
[docs] def get_surface(self):
"""
return : [cm]^2
"""
array_size = self.get_sizecm()
return array_size[0] * array_size[1]
[docs] def pos_rect_array(self):
'''
Return current rectangle [[x_lb, y_lb], [x_ru, y_ru]] where
lb: left bottom, ru: right upper in cm
:return: rectangle [[x_lb, y_lb], [x_ru, y_ru]]
:rtype: array 2x2
'''
rect = np.empty((2, 2), dtype=np.float32)
rect[0,:] = self.pos_corner_lb()
rect[1,:] = rect[0,:] + self.get_sizecm()
return rect
[docs] def pos_rect_cell(self, idx):
"""
for one cell
"""
rect = np.empty((2, 2), dtype=np.float32)
a_idx = np.array([idx])
rect[0,:] = self.idx2pos_inf(a_idx)
rect[1,:] = rect[0,:] + self.size_cm
return rect
[docs] def pos_rect_cells(self, aidx):
"""Returns physical position of the lower left and
upper right corners.
:param aidx: list or numpy array-like containing the indices
of the lower left corner
:type aidx: [[ll_x_idx1, ll_y_idx1], [ll_x_idx2, ll_y_idx2], ...]
:return: list of the physical coord. of the ll and ur corners.
:rtype: list of numpy-like arrays
"""
n_rect = aidx.shape[0]
arect = np.empty((n_rect, 2, 2), dtype=np.float32)
arect[:,0,:] = self.idx2pos_inf(aidx)
arect[:,1,:] = self.idx2pos_sup(aidx)
return arect
[docs] def pos_corner_lb(self):
"""Return position corner 'l'eft 'b'ottom in absolute referential
:return: [cm] position corner left bottom
:rtype: array (2,)
"""
return -self.to_origin_cm
[docs] def pos_corner_ru(self):
"""Return position corner 'r'ight 'u'pper in absolute referential
:return: [cm] position corner right upper
:rtype: array (2,)
"""
return self.get_sizecm() - self.to_origin_cm
[docs] def pos_center_cells(self):
"""Return 2 array x and y with center cells
.. warning: to be finished.
:return: 2 array x and y with center cells
:rtype: numpy array
"""
x_center = np.arange(self.get_nb_cell_x()) * self.size_cm + self.size_cm / 2
x_center -= self.to_origin_cm[0]
y_center = np.arange(self.get_nb_cell_y()) * self.size_cm + self.size_cm / 2
y_center -= self.to_origin_cm[1]
return (x_center, y_center)
[docs] def array_pos_center_cells(self):
"""Return 2 array x and y with center cells
.. warning: to be finished.
:return: 2 array x and y with center cells
:rtype: numpy array
"""
x_center, y_center = self.pos_center_cells()
a_center = np.empty((self.get_nb_cell_x(), self.get_nb_cell_y(), 2), dtype=np.float)
for i_x in range(self.get_nb_cell_x()):
a_center[i_x,:,0] = x_center[i_x]
a_center[i_x,:,1] = y_center
return a_center
[docs] def value_at_pos_check(self, pos):
'''
Value in cell associated to pos, with limit check
:param pos: x,y in array
'''
idx = self.pos2idx_inf(pos)
idx0 = idx[0]
if self.l_verb > 0:
self.logger.debug(f'max: {self.get_nb_cell_x()}, {self.get_nb_cell_y()}')
if idx0 < 0 or idx0 >= self.get_nb_cell_x():
self.logger.error('index out of range.')
return None
idx1 = idx[1]
if idx1 < 0 or idx1 >= self.get_nb_cell_y():
self.logger.error('index out of range.')
return None
return self.value(idx)
[docs] def value_at_pos(self, pos):
"""Return the value at a given
physical position.
:param pos: physical position
:type pos: [ll_x, ll_y]
:return: value
:rtype: float
"""
idx = self.pos2idx_inf(pos)
return self.value(idx)
[docs] def value(self, idx):
"""Return value of a given pixel index.
:param idx: pixel index
:type idx: [idx_x, idx_y]
:return: value
:rtype: float
"""
return self.ar_val[idx[0], idx[1]]
[docs] def value_array(self, aidx: np.array):
"""Return value at a given set of pixel indices.
Generalisation of value() for several indices.
:param idx: pixel indices
:type idx: [[idx_x, idx_y], ... ]
:return: values
:rtype: array(float)
"""
return self.ar_val[aidx[:,0], aidx[:,1]]
[docs] def get_shape(self):
return self.ar_val.shape
[docs] def set_val(self, p_val):
'''
Set p_val for all cells of array
:param p_val:
:type p_val: number
'''
self.ar_val[:,:] = p_val
[docs] def get_ij_cell_surrounded_no_0(self):
a_ij = np.empty((0, 2), dtype=np.int16)
p_ar = self.ar_val
for p_i in range(1, p_ar.shape[0] - 1):
for p_j in range(1, p_ar.shape[1] - 1):
self.logger.debug(f'p_i: {p_i} / p_j: {p_j}')
if (p_ar[p_i + 1, p_j - 1] * p_ar[p_i + 1, p_j] * p_ar[p_i + 1, p_j + 1] == 0.0):
continue
if (p_ar[p_i, p_j - 1] * p_ar[p_i, p_j + 1] == 0.0):
continue
if (p_ar[p_i - 1, p_j - 1] * p_ar[p_i - 1, p_j] * p_ar[p_i - 1, p_j + 1] == 0.0):
continue
a_ij = np.vstack((a_ij, np.array([p_i, p_j], dtype=np.int16)))
return a_ij
[docs] def pos2idx_rect(self, p_rect, p_boundary=False):
"""
return index rectangle associated to pos rectangle p_rect
:param p_rect: rectangle
:type p_rect:
:param p_boundary:
:type p_boundary:
"""
r_idx = np.empty((2, 2), dtype=np.int16)
r_idx[0] = self.pos2idx_inf(p_rect[0])
r_idx[1] = self.pos2idx_inf(p_rect[1])
if p_boundary:
r_idx = np.where(r_idx < 0, 0, r_idx)
r_idx[:,0] = np.where(r_idx[:,0] > self.max_ix, self.max_ix, r_idx[:,0])
r_idx[:,1] = np.where(r_idx[:,1] > self.max_iy, self.max_iy, r_idx[:,1])
return r_idx
[docs] def idx2pos(self, array_idx, pos="CTR"):
"""return array position
:param array_idx: array index
:type array_idx: int
:param pos: pos = {Center: CTR, left bottom :LB, right up :RU}
:type pos: string
"""
array_idx_c = array_idx.copy().astype(np.float32)
if pos == "CTR":
array_idx_c += np.array([0.5, 0.5], dtype=np.float32)
elif pos == "RU":
array_idx_c += np.array([1, 1])
pos = array_idx_c * self.size_cm - self.to_origin_cm
return pos
[docs] def idx2pos_inf(self, array_idx):
return self.idx2pos(array_idx, "LB")
[docs] def idx2pos_sup(self, array_idx):
return self.idx2pos(array_idx, "RU")
[docs] def pos2idx_noround(self, pos_cm):
return (pos_cm + self.to_origin_cm) / self.size_cm
[docs] def pos2idx(self, pos_cm, round_func=None):
if round_func == None:
raise
array_idx = (pos_cm + self.to_origin_cm) / self.size_cm
idx = round_func(array_idx).astype(np.int16)
return idx
[docs] def pos2idx_inf(self, pos_cm):
"""Return index x,y lower left corner of cell containing pos_cm
"""
ret = self.pos2idx(pos_cm, np.floor)
return ret
[docs] def pos2idx_sup(self, pos_cm):
return self.pos2idx(pos_cm, np.ceil)
[docs] def random_value(self, pmin, pmax):
aa = np.random.randint(pmin, pmax + 1, self.ar_val.size).reshape((self.get_nb_cell_x(), -1))
self.set_array(aa)
[docs] def random_pos(self, nb_pos):
"""
random position in all array
:param nb_pos: number position to random
:type nb_pos: int
"""
msize = self.get_sizecm()
apos = np.empty((nb_pos, 2), dtype=np.float32)
epsilon32 = np.finfo(np.float32).eps
apos[:,0] = np.random.uniform(0, msize[0] * (1 - epsilon32), nb_pos)
apos[:,1] = np.random.uniform(0, msize[1] * (1 - epsilon32), nb_pos)
retpos = apos - self.to_origin_cm
return retpos
[docs] def random_pos_rect(self, nb_pos, rec):
"""
random position in rectangle given by parameter
return : nb_pos random position in rec
"""
# n = len(nb_pos)
# epsilon = np.finfo(np.float32).eps
# apos = np.empty((n,2), dtype=np.float32)
# apos[:,0] = np.random.uniform(rec[0,0]*(1+epsilon), rec[1,0]*(1-epsilon), n)
# apos[:,1] = np.random.uniform(rec[0,1]*(1+epsilon), rec[1,1]*(1-epsilon), n)
# return apos
epsilon = np.finfo(np.float32).eps
apos = np.empty((nb_pos, 2), dtype=np.float32)
apos[:,0] = np.random.uniform(rec[0,0] * (1 + epsilon), rec[1,0] * (1 - epsilon), nb_pos)
apos[:,1] = np.random.uniform(rec[0,1] * (1 + epsilon), rec[1,1] * (1 - epsilon), nb_pos)
return apos
[docs] def intersection(self, o_asc, vec_trans):
"""
Return the intersection rectangle between current array and o_asc translated of vec_trans
:param o_asc: array
:type o_asc : ArraySquareCell
:param vec_trans: translation to apply at o_asc
:type vec_trans : vector
:return: intersection rectangle between current array and o_asc translated of vec_trans
:rtype: rectangle [[x_lb, y_lb], [x_ru, y_ru]] lb: left bottom, ru: right upper in cm
"""
r1 = self.pos_rect_array()
r2 = o_asc.pos_rect_array()
r2[0, :] += vec_trans
r2[1, :] += vec_trans
return intersection_rect(r1, r2)
[docs] def sample_up(self, pfact):
"""
pfact is int
result: repeat pfactxpfact each element
"""
o_asc = ArraySquareCell(self.get_nb_cell_x() * pfact, self.get_nb_cell_y() * pfact, \
self.size_cm / pfact)
su_ar = self.ar_val.repeat(pfact).reshape((self.get_nb_cell_x(), \
self.get_nb_cell_y() * pfact)).T
su_ar = su_ar.repeat(pfact).reshape((self.get_nb_cell_x() * pfact, \
self.get_nb_cell_y() * pfact)).T
o_asc.set_array(su_ar)
return o_asc
[docs] def plot_ijdet(self, p_mes=""): # pragma: no cover
plt.figure()
plt.title(p_mes)
plt.imshow(yzegp_to_ijdet_array(self.ar_val))
plt.colorbar()
[docs] def plot(self, p_mes=""): # pragma: no cover
plt.figure()
plt.title(p_mes)
plt.imshow(np.transpose(self.ar_val), origin='lower')
plt.colorbar()