'''------------------------------------------------------------ Title: Fibre sampling of discretized beam cross section Date: 25/11/2021 Author: J. K. Patel Version: 1.0 ################################################################# Fuctions for simulating the readout of fibres which sample a discretized plane in a beam of protons. ################################################################## -------------------------------------------------------------''' import numpy as np import pickle as pkl import os from util import * def sampled_plane(plane_size, element_size): ''' Returns sparse co-ordinates of specified size and density, corresponding to an 'image' plane at which the beam is sampled, with origin at centre. Parameters ----------- plane_size: float, int Length of dimension of square plane. element_size: float, int Size/length of sampling elements in plane. Returns -------- xo: 2D array (row) x co-ordinates in discretized sampled plane yo: 2D array (column) y co-ordinates of discretized sampled plane ''' # check that elements are smaller than plane if element_size >= plane_size: raise ValueError('elementsize ({}) must be smaller than planesize ({})'.format(element_size, plane_size)) # create co-ordinate matrices representing locations of beam elements #x, y = np.arange(0,planesize,planeLengthElements), np.arange(0,planesize,planeLengthElements) yy, xx = np.ogrid[0:plane_size:element_size, 0:plane_size:element_size] # set origin at center centreIndex = plane_size/2 xo, yo = xx - centreIndex, yy - centreIndex return xo, yo def fibre_mask_with_indices(x, y, centre, width): ''' Generates single boolean fibre along x co-ordinate axis (rotated from horizontal if x contains rotated co-ordinates), centred at y co-ordinate, 'centre', and with width in the co-ordinate space dimensions, 'width'. ''' boolFibre = np.zeros_like(x*y) # indices of fibre in x,y cordinates fibreXY = np.where((y<=centre+width/2)&(y>=centre-width/2)&(np.abs(x)<=np.max(np.abs(x)))) boolFibre[fibreXY] = 1 return boolFibre, fibreXY def round_fiber_Xsection(x, r): return normalize(np.nan_to_num( 2*np.sqrt(r**2-x**2))) def roundFiber_with_indices(x,y,center,dia): # indices of fibre in x,y cordinates fibreXY = np.where((y<=center+dia/2)&(y>=center-dia/2)&(np.abs(x)<=np.max(np.abs(x)))) return np.ones_like(x)*round_fiber_Xsection(y-center, dia/2), fibreXY def get_coarse_fibre_indices(fibre_indices, fine_pixels_per_coarse_pixel): ''' Returns list of indices into pixels in sampled plane which coincide with fibre area. ''' # get indices of coarse array where finely sampled fibre is 'True' fibreIndices = np.vstack((fibre_indices[0],fibre_indices[1])) # just getting as np.array from tuple # divide fine index by number of fine elements in a coarse pixel, then round down to integer: this is coarse index coarseIndex = np.floor(fibreIndices/fine_pixels_per_coarse_pixel) coarseIndex = np.array(np.unique(coarseIndex,axis=1), dtype = np.int32).T # get unique locations, transpose for nicer format return coarseIndex def coarse_pixel_weightings(bool_fibre, coarse_fibre_indices, n_fine_in_coarse_pixel): ''' Returns list of weightings for pixels in sampled plane which coincide with fibre area, where weighting at each pixel is the fraction of the pixel which is occupied by fibre area ''' coarseWindow = np.zeros_like(bool_fibre) # empty array we will use to pass a top-hat over coarse pixels ratioList = [] # list to add fill ratios to # helper fine index into coarse pixels, as long as fine array, but in nfineInCoarsePixel steps i = np.arange(0, bool_fibre.shape[0]+1, n_fine_in_coarse_pixel) # loop through coarse fibre pixels for Y, X in coarse_fibre_indices: # x, y range in fine pixels x0, x1 = i[X], i[X+1] y0, y1 = i[Y], i[Y+1] coarseWindow[y0:y1, x0:x1] = 1 # top-hat (1 in coarse pixel, zero everywhere else) # window finely sampled fibre coarseWindow *= bool_fibre fillRatio = np.sum(coarseWindow)/n_fine_in_coarse_pixel**2 ratioList.append(fillRatio) coarseWindow = np.zeros_like(bool_fibre) # reset to zero everywhere return ratioList def weight_coarse_pixels(ratios, coarse_fibre_indices, coarse_fibre): ''' Returns coarse fibre with pixels weighted by fraction of the pixel which is occupied by fibre area ''' #check ratios and fibre indices in coarse grid match assert len(ratios) == coarse_fibre_indices.shape[0], \ f'length of weightings ({len(ratios)}) should be the same as number of coarse fibre indices ({coarse_fibre_indices.shape[0]})' # initialise weighted coarse matrix with ones like coarse fibre boolean matrix weightedCoarseFibre = coarse_fibre.copy() # loop through ratios in list and weight corresponding coarse pixel in bool matrix for en, r in enumerate(ratios): # go to coarse fibre = 1 weightedCoarseFibre[coarse_fibre_indices[en,0], coarse_fibre_indices[en,1]] *= r return weightedCoarseFibre def rotate_coordinates(xo, yo, angle = 0, shift_origin = None): ''' Returns co-ordinates rotated clockwise by some 'angle' in degrees ''' if shift_origin != None: xo += shift_origin[1] yo += shift_origin[0] angleRad = np.deg2rad(angle) xrtd = xo*np.cos(angleRad) + yo*np.sin(angleRad) yrtd = -xo*np.sin(angleRad) + yo*np.cos(angleRad) if shift_origin != None: xrtd -= shift_origin[1] yrtd -= shift_origin[0] return xrtd, yrtd def fibre_centres(xo, yo, fibre_width, fibre_sep, fibre_num = None): ''' Returns 1D array of perpendicular distances to fibre centres, either by setting: 'fibre_sep': separation between fibres -- will return the maximum number of fibres, or; 'fibre_num': the number of fibres in a panel -- will adjust the separation of fibres accordingly Prioritises fibre_sep if passed as a value not 'None' ''' def max_coord(xo,yo): return max(np.abs(np.max(xo)), np.abs(np.min(xo)), np.abs(np.max(yo)), np.abs(np.min(yo))) # limits of fibre centre positions range = max_coord(xo, yo) lim = range - fibre_width/2 if fibre_sep != None: # step between fibre centres dy = fibre_width + fibre_sep centres = np.arange(-lim, lim, dy) elif fibre_num != None: assert fibre_num*fibre_width <= 2*lim, f'Too many fibres ({fibre_num}) of width ({fibre_width}). Use fewer, or narrower fibres' centres = np.linspace(-lim, lim, fibre_num) return centres def fibre_panel(angle, fibre_centres, fibre_width, fine_element = 0.01, coarse_element = 1, plane_size = 25): ''' Returns: -------- fibrePanel: full panel of fibres (high resolution) rList: coarse fibre pixel weighting ratios coarseIndexList: coarse fibre pixel indices ''' # generate fine co-ordinate grid, and rotate by angle xo, yo = sampled_plane(plane_size, fine_element) xr, yr = rotate_coordinates(xo,yo, angle) nfineInCoarsePixel = int(coarse_element/fine_element) # number of fine elements in coarse pixel # initialize empty lists for ratios and coarse fibre indices rList, coarseIndexList = [], [] # empty array to fill with full res panel of boolean fibres fibrePanel = np.zeros_like(xr*yr) # loop over fibres by relative position of centers for centre in fibre_centres: # get fibre and indices and add to panel boolFibre, fibreIndices = roundFiber_with_indices(xr,yr, centre, fibre_width) fibrePanel += boolFibre # index into a coarsely sampled fibre grid using no. of fine elements in a coarse pixel coarseIndices = get_coarse_fibre_indices(fibreIndices, nfineInCoarsePixel) coarseIndexList.append(coarseIndices) # get ratio of 'filled' fine pixel elements in each coarse pixel ratios = coarse_pixel_weightings(boolFibre, coarseIndices, nfineInCoarsePixel) rList.append(ratios) return fibrePanel, rList, coarseIndexList def save_ratios_and_indices(dir, r_list, index_list, angle): ''' Saves lists of ratios and coarse indices to pickle files ''' angleString = str(angle) if not os.path.exists(f'{dir}/{angleString}'): os.mkdir(f'{dir}/{angleString}') with open(os.path.abspath(f'{dir}/{angleString}/PSSI_{angleString}deg_ratioList_JP_v1.pkl'),'wb') as file: pkl.dump(r_list,file) with open(os.path.abspath(f'{dir}/{angleString}/PSSI_{angleString}deg_indexList_JP_v1.pkl'),'wb') as file2: pkl.dump(index_list,file2) def loadRatiosAndPositions(dir, angle): ''' Loads lists of ratios and coarse indices from pickle files ''' angleString = str(angle) with open(os.path.abspath(f'{dir}/{angleString}/PSSI_{angleString}deg_ratioList_JP_v1.pkl'), 'rb') as file: ratios = pkl.load(file) with open(os.path.abspath(f'{dir}/{angleString}/PSSI_{angleString}deg_indexList_JP_v1.pkl'),'rb') as file2: indexes = pkl.load(file2) return ratios,indexes def save_fibre_panel(dir, angle, fibre_width, fibre_sep, fine_element = 0.01, coarse_element = 1, plane_size = 25, fibre_num = None): ''' Generates single panel of parallel fibres, and saves minimalist information to describe their position in a sampled plane into a new folder Parameters: ----------- dir: str Directory in which to save fibre positions and weightings. angle: float Angle at which to sample plane. fibre_width: float Width of scintillating fibres. fibre_separation: float Separation of scintillating fibres. If 'None', will use 'fibre_num' parameter to determine number of fibres and positions. fine_element: float Resolution at which to generate fibres to sample more coarsely discretized plane (should be less than 'coarse_element'). coarse_element: float Resolution of sampled plane. Should be > 'fine_element'. plane_size: float Length of dimension of square plane to be sampled at 'coarse_element' resolution. fibre_num: int Number of fibres in sampling plane - to be used if 'fibre_separation' = None. Default is to use 'fibre_separation' and 'fibre_num' = None. ''' # create fine meshgrid and set fibre width and separation xo, yo = sampled_plane(plane_size, fine_element) # get list of fibre centres in a panel rotated at 'angle' to horizontal centres = fibre_centres(xo, yo, fibre_width, fibre_sep, fibre_num) # get list of ratios and corresponding indices into coarse co-ordinate grid rList, coarseIndexList = fibre_panel(angle = angle, fibre_centres = centres, fibre_width = fibre_width, fine_element = fine_element, coarse_element = coarse_element)[1:] save_ratios_and_indices(dir, rList, coarseIndexList, angle) def load_layer_of_panels(dir, angles): ''' Loads ratios and corresponding indices into coarsely sampled grid from pickle files labelled by angle of each panel. ''' # empty lists to fill with lists of ratios for each panel/angle layerRatioList, layerIndexList = [], [] for a in angles: ratioList, indexList = loadRatiosAndPositions(dir, a) layerRatioList.append(ratioList) layerIndexList.append(indexList) return layerRatioList, layerIndexList def save_layer_of_panels(parentDir, angles, fibre_width, fibre_separation, fine_element = 0.01, coarse_element = 1, plane_size = 25, fibre_num = None): ''' Generates panels of parallel fibres at series of angles, and saves minimalist information to describe their position in a sampled plane into a new folder Parameters: ----------- parentDir: str Directory in which to create new folder containing fibre positions and weightings. angles: 1D array, float Angles at which to sample plane. fibre_width: float Width of scintillating fibres. fibre_separation: float Separation of scintillating fibres. If 'None', will use 'fibre_num' parameter to determine number of fibres and positions. fine_element: float Resolution at which to generate fibres to sample more coarsely discretized plane (should be less than 'coarse_element'). coarse_element: float Resolution of sampled plane. Should be > 'fine_element'. plane_size: float Length of dimension of square plane to be sampled at 'coarse_element' resolution. fibre_num: int Number of fibres in sampling plane - to be used if 'fibre_separation' = None. Default is to use 'fibre_separation' and 'fibre_num' = None. Returns: -------- newDir: str Name of new directory into which fibre positions and weightings have been saved. ''' # create new folder in directory with name given by passed parameters if not os.path.exists(parentDir): os.mkdir(parentDir) if fibre_num == None: fibre_num = plane_size//(fibre_width + fibre_separation) newDir = os.path.abspath(f'{parentDir}fWidth{fibre_width}_nFibres{fibre_num}_res{coarse_element}mm_planeSize{plane_size}') if not os.path.exists(newDir): os.mkdir(newDir) # write .txt file with meta data with open(f'{newDir}/metadata.txt', 'w') as f: f.write(f'No. of panels = {len(angles)}\n'\ f'angles = {angles}\n'\ f'fine element [mm] = {fine_element}\n'\ f'coarse element [mm] = {coarse_element}\n'\ f'sampled plane size [mm] = {plane_size}\n'\ f'fibre width [mm] = {fibre_width}\n'\ f'fibre separation [mm] = {fibre_separation}\n' f'number of fibres = {fibre_num}') print('Metadata saved.') # create fine meshgrid and set fibre width and separation xo, yo = sampled_plane(plane_size, fine_element) # get list of fibre centres in a panel rotated at 'angle' to horizontal centres = fibre_centres(xo, yo, fibre_width, fibre_separation, fibre_num) # loop through panels at different angles for en, a in enumerate(angles): # get list of ratios and corresponding indices into coarse co-ordinate grid, then save to pickle files labelled by angle fibrePanel, rList, coarseIndexList = fibre_panel(a, fibre_centres = centres, fibre_width = fibre_width, fine_element = fine_element, coarse_element = coarse_element, plane_size = plane_size) save_ratios_and_indices(newDir, rList, coarseIndexList, a) print(f'Panel {en+1} at {a} degrees saved') return newDir def atten(xr, decay_length, dist_to_detector = 0): ''' Exponential attenuation function for signal attenuation in scintillating fibre and fibre readout to sensor. Parameters ----------- xr: 2D array x co-ordinates of sampled plane in panel frame of reference (i.e. will be rotated by panel angle) decay_length: float signal attenuation length of scintillating fibres and also readout fibres (carrying signal from sampled plane to sensor) dist_to_detector: float length of the readout fibre from edge of sampled plane to sensor Returns -------- attenuation_profile: 2D array signale decay profile along passed 'xr' co-ordinate axis ''' lenDiag = np.sqrt(2*xr.shape[0]**2)/2 # half max length of fibre across plane # add length of fibre from edge of detection plane to detector fibreLength = dist_to_detector + lenDiag attenuation_profile = np.exp(-(xr+fibreLength)/decay_length) return attenuation_profile def atten_fibre(fibre, xr, decay_length, dist_to_detector = 0, direction = 'l'): ''' Masks a fibre weight vector (i.e. image of a fibre in discretized sample plane, weighted by pixel fill fraction) by an attenuation function which depends on the direction of readout. Parameters ----------- fibre: 2D array weight vector corresponding to an image of a fibre in discretized sample plane, weighted by pixel fill fraction xr: 2D array x co-ordinates of sampled plane in panel/fibre frame of reference (i.e. will be rotated by panel/fibre angle) decay_length: float signal attenuation length of scintillating fibres and also readout fibres (carrying signal from sampled plane to sensor) dist_to_detector: float length of the readout fibre from edge of sampled plane to sensor Returns -------- attenuated_fibre: 2D array product of attenuation mask from 'atten()' and passed fibre weight vector (image) ''' # for detection on the left, signal shoud be attenuated more for pixels on the right, so flip x values if direction == 'r': xr, _ = rotate_coordinates(xr, xr, 180) attenuated_fibre = fibre*atten(xr, decay_length, dist_to_detector) return attenuated_fibre def projection(image, resolution, angle, ratio_list, index_list, attenuation = True, dist_to_detector = 0, decay_length = 25): ''' Returns projection of input image along parallel fibres in panel at specified angle Parameters ----------- image: 2D array Discretized representation of a beam cross-section in sampled plane angle: float Angle at which to sample plane. ratio_list: list, floats List of weightings corresponding to ratio of fibre area coinciding with a pixel in sampled plane to the total pixel area index_list: list, int List of indices corresponding to pixels in the sampled plane at which fibre area coincides attenuation: bool Whether to account for attenuation effects of fibre length on optical signal. If 'True' will return projections from both ends of fibres. dist_to_detector: float length of the readout fibre from edge of sampled plane to sensor decay_length: float signal attenuation length of scintillating fibres and also readout fibres (carrying signal from sampled plane to sensor) Returns -------- prj: 1D/2D array, floats Array of 'ray sums': integrals of 'image' signal along fibres. If 'attenution' = True, will have shape (2, fibre_num), where first element is projection to a detector in the negative x-direction (i.e. 'left') when angle = 0, and second element is projection to a detector on the 'right'. ''' imgSize = image.shape[0] planeSize = int(imgSize*resolution) nfibres = len(index_list) x, y = sampled_plane(planeSize, resolution) xr, yr = rotate_coordinates(x, y, angle) if attenuation == True: prj = np.zeros((2,len(index_list))) # empty array to fill with fibre integrals panel_projection_matrix = np.zeros((2, nfibres, imgSize, imgSize)) else: prj = np.zeros(len(index_list)) # only need one projection if we're not attenuating panel_projection_matrix = np.zeros((nfibres, imgSize, imgSize)) # loop through fibres for fibre in range(nfibres): # get x,y indices, and corresponding ratios of coarse elements with overlapping fibre fX, fY = index_list[fibre][:,1], index_list[fibre][:,0] # create boolean fibre image, then weight with ratios coarseFibre = np.zeros((imgSize, imgSize)) coarseFibre[fY, fX] = 1 weightedCoarseFibre = weight_coarse_pixels(ratio_list[fibre], index_list[fibre], coarseFibre) if attenuation == True: # get both left- and right-hand projections accounting for attenuation left_atten_fibre = atten_fibre(weightedCoarseFibre, xr, decay_length, dist_to_detector, direction = 'l') right_atten_fibre = atten_fibre(weightedCoarseFibre, xr, decay_length, dist_to_detector, direction = 'r') panel_projection_matrix[0,fibre,:,:] = left_atten_fibre panel_projection_matrix[1,fibre,:,:] = right_atten_fibre prj[0,fibre] = np.sum(image[fY, fX]*left_atten_fibre[fY, fX]) prj[1,fibre] = np.sum(image[fY, fX]*right_atten_fibre[fY, fX]) else: # get single projection, not accounting for attenuation prj[fibre] = np.sum(image[fY, fX]*weightedCoarseFibre[fY, fX]) panel_projection_matrix[fibre,:,:] = weightedCoarseFibre return prj, panel_projection_matrix def sinogram(image, resolution, angles, RList, IList, parent_dir = None, attenuation = True, decay_length = 25, noise = 0): ''' Get full sinograms of input image, along with weightings (ratios) and positions (indices) corresponding to fibre vectors (image in discretized sample plane) Parameters ----------- image: 2D array Discretized representation of a beam cross-section in sampled plane angles: float Angles at which to sample plane with panels. parent_dir: str Directory from which to load fibre weightings (ratios) and positions (indices) in sampled plane attenuation: bool Whether to account for attenuation effects of fibre length on optical signal. If 'True' will return projections from both ends of fibres. decay_length: float signal attenuation length of scintillating fibres and also readout fibres (carrying signal from sampled plane to sensor) Returns -------- sino: 2D array 'Ray sums' from each fibre (axis 1) at each panel angle (axis 0) RList: list, floats Nested list of fibre weightings (ratios) indexed like [panel][fibre] IList: list, floats Nested list of fibre positions (indices) indexed like [panel][fibre] ''' sino = [] projection_matrix = [] if RList == None and IList == None: RList, IList = load_layer_of_panels(parent_dir, angles) for en, angle in enumerate(angles): rList, iList = RList[en], IList[en] sino.append(projection(image, resolution, angle, rList, iList, attenuation = attenuation, decay_length = decay_length)[0]) projection_matrix.append(projection(image, resolution, angle, rList, iList, attenuation = attenuation, decay_length = decay_length)[1]) # print(f'Panel {en+1} complete.') # rearrange matrix and sinogram projection_matrix = np.swapaxes(np.array(projection_matrix), 0, 1) sino = np.swapaxes(np.array(sino), 0,1) # add gaussian noise if noise != 0: noise = np.random.normal(0, noise, sino.shape) sino += noise return sino, RList, IList, projection_matrix def back_projection(sinogram, rList, iList, panelSize, resolution, attenuation = True): ''' Returns simple back projection of fibre 'ray-sums', or returns a pair of back-projections if pair of sinograms (left- and right-hand readout) are passed. ''' # empty array with shape of image to reconstruct imgSize = int(panelSize/resolution) if attenuation == True: # if attenuation = 'True' then will get both left- and right-hand projection reconstructions backProjection = np.zeros((2, imgSize, imgSize)) else: backProjection = np.zeros((imgSize, imgSize)) # helper array to update back-projection ratioArray = np.zeros((imgSize, imgSize)) # loop through panels for panel in range(sinogram.shape[0]): ratios, indices = rList[panel], iList[panel] # loop through fibres for fibre in range(sinogram.shape[-1]): fX, fY = indices[fibre][:,1], indices[fibre][:,0] # get fibre indices ratioArray[fY,fX] = ratios[fibre] # set fibre elements to ratios # create images with fibre elements of values = fibre integral * ratio if attenuation == True: # when accounting for attenuation, get both L and R sinograms prjL = sinogram[panel,0,fibre] prjR = sinogram[panel,1,fibre] backProjection[0][fY,fX] += prjL*ratioArray[fY,fX] backProjection[1][fY,fX] += prjR*ratioArray[fY,fX] else: prj = sinogram[panel,fibre] backProjection[fY,fX] += prj*ratioArray[fY,fX] return backProjection def get_prj_matrix(rList, iList, panelSize, resolution): ''' Returns projection matrix with each fiber as a row vector with size of image ''' # empty array with shape of image to reconstruct imgSize = int(panelSize/resolution) # flatten ratios and indices into 1D arrays/vectors ratios, indices = rList, iList nFibres = len(ratios[0]) nPanels = len(ratios) # create projection matrix to fill with fibre ratio values prjMatrix = np.zeros((nFibres*nPanels, imgSize, imgSize)) panels = range(nPanels) resorted_panels = [] for a in range(nPanels//2): resorted_panels.append(panels[a::nPanels//2]) panels_rsrt = np.array(resorted_panels).ravel() f = 0 for panel in panels_rsrt: # loop through fibres for fibre in range(nFibres): fX, fY = indices[panel][fibre][:,1], indices[panel][fibre][:,0] # get fibre indices prjMatrix[f,fY,fX] = ratios[panel][fibre] # set fibre elements to ratios f += 1 return np.asmatrix(prjMatrix.reshape((nFibres*nPanels, imgSize**2)))