vamtoolbox.projector.pyTorchRayTrace#
Module Contents#
Classes#
pyTorch implementation of custom ray tracing operation |
|
ODE Solver to iterate over ray operations |
|
State of a set of rays. The set can be all rays, rays emitted per angle, per z-slice, or per angle-z-slice |
- class vamtoolbox.projector.pyTorchRayTrace.PyTorchRayTracingPropagator(target_geo, proj_geo, output_torch_tensor=False)#
pyTorch implementation of custom ray tracing operation This operator is the A in Ax = b, or the P in Pf = g, where x or f is the real space object (flattened into 1D array) and b or g is the sinogram (flattened into 1D array). Here naming convention Pf = g is used to avoid collision with x which commonly denote position of ray in ray tracing literature.
Since storing all curved ray path is prohibitively memory-intensive, this tracing operation is designed to have as low memory requirement as possible. Compute on GPU by default but fallback to CPU computation if GPU is not found. Currently occulsion is not supported yet.
Internally, 2D and 3D problem is handled identically. Conditional handling only happens in I/O. #Implement ray_setup, coordinate systems
- forward(f)#
- backward(ray_energy_0)#
- inverseBackward(f, filter_option='ram-lak', **kwargs)#
This function compute the inverse of backpropagation. Similar to forward propagation, this function maps from a spatial quantity to a sinogram quantity.
- buildPropagationMatrix()#
This function builds the forward propagation matrix that corresponds to the forward propagation operation in discrete form. This operator is the A in Ax = b, or the P in Pf = g, where x or f is the real space object (flattened into 1D array) and b or g is the sinogram (flattened into 1D array). This propagation matrix is returned as a pytorch COO tensor if “output_torch_tensor” attribute is True, and as a scipy sparse COO array otherwise.
- perAngleTrace()#
- perZTrace()#
- perRayTrace()#
- class vamtoolbox.projector.pyTorchRayTrace.RayTraceSolver(device, index_model, attenuation_model, absorption_model, ray_trace_method, eikonal_parametrization, ode_solver, max_num_step=750, num_step_per_exit_check=10)#
ODE Solver to iterate over ray operations
Store defined step size, selected stepping scheme
index and its gradient field
- solveUntilExit(ray_state, step_size, callback=None, tracker_on=False, track_every=5)#
- depositEnergyUntilExit(ray_state, step_size, callback=None, tracker_on=False, track_every=5)#
- integrateEnergyUntilExit(ray_state, step_size, real_space_distribution, callback=None, tracker_on=False, track_every=5)#
- recordEnergyUntilExit(ray_state, step_size, tracker_on=False, track_every=5)#
This function record the energy deposited to voxels by rays of unity energy. The result is recorded as a pyTorch sparse COO tensor.
- _dx_dsigma(_, v, *arg)#
- _dv_dsigma(x, _, known_n=None)#
- _dx_ds(x, v, known_n=None)#
- _dv_ds(x, _, *arg)#
- _dx_dopl(x, v, known_n=None)#
- _dv_dopl(x, _, known_n=None)#
- _forwardSymplecticEuler(ray_state, step_size)#
- _forwardEuler(ray_state, step_size)#
- _leapfrog_init(ray_state, step_size)#
- _no_step_init(ray_state, _)#
- deposit(ray_state, step_size, desposition_grid)#
- integrate(ray_state, step_size, real_space_distribution)#
- record(ray_state, step_size, propagation_matrix_shape)#
- depositEnergyOnAdjacentVoxel(ray_energy, step_size_idx_unit, x_arr_idx, voxel_idx_x, adjacent_voxel_number: int, desposition_grid)#
In 3D problems, this function handles the dose deposition of one particular voxel (lower/upper corner) out of 8 voxels adjacent to positions at x. (In 2D problems, this function handles one out of 4 adjacent voxels.)
The adjacent_voxel_number is integer from 0-7 for 3D problem and 0-3 for 2D problem
This function works for any number of rays, as long as the shape[0] of all ray-numbered inputs are consistent. Internally, it only process rays that have a valid adjacent voxel position.
- integrateEnergyFromAdjacentVoxel(ray_energy, step_size_idx_unit, x_arr_idx, voxel_idx_x, adjacent_voxel_number: int, real_space_distribution, ray_integral)#
In 3D problems, this function handles the dose deposition of one particular voxel (lower/upper corner) out of 8 voxels adjacent to positions at x. (In 2D problems, this function handles one out of 4 adjacent voxels.)
The adjacent_voxel_number is integer from 0-7 for 3D problem and 0-3 for 2D problem
This function works for any number of rays, as long as the shape[0] of all ray-numbered inputs are consistent. Internally, it only process rays that have a valid adjacent voxel position.
- recordEnergyOnAdjacentVoxel(ray_energy, step_size_idx_unit, x_arr_idx, voxel_idx_x, adjacent_voxel_number: int, propagation_matrix_shape, active_ray_idx)#
In 3D problems, this function handles the dose deposition of one particular voxel (lower/upper corner) out of 8 voxels adjacent to positions at x. (In 2D problems, this function handles one out of 4 adjacent voxels.)
The adjacent_voxel_number is integer from 0-7 for 3D problem and 0-3 for 2D problem
This function works for any number of rays, as long as the shape[0] of all ray-numbered inputs are consistent. Internally, it only process rays that have a valid adjacent voxel position.
voxel_idx_x is multi-dimensional index in real space domain active_ray_idx is 1-dimensional index of rays
- expressPositionInArrayIndices(x)#
Express position x in grid indices
- getAdjacentVoxelIndicesAtLocation(x_arr_idx: torch.Tensor)#
Get adjacent voxel indices. Output indices are presented by binary permutation of lower and upper indices in every dimension Note: Not all rays are required. Able to only process the activate ray segments.
- linearCombinationRayEnergyToGrid(desposition_grid_shape, voxel_idx_select_at_valid_idx, energy_interaction_at_valid_idx)#
This function builds a sparse tensor for depositing ray energy onto the grid. The tensor adds up the contribution of multiple rays towards the same elements. This avoid the racing condition of grid elements in direct insertion. The first dimension matches the numel() of the deposition grid and the second dimension matches the number of rays that have a valid adjacent voxel index.
- coalesceRayEnergyToGrid(desposition_grid_shape, voxel_idx_select_at_valid_idx, energy_interaction_at_valid_idx)#
- discreteSurfaceIntersectionCheck()#
- exitCheck(ray_state)#
When the rays are outside the bounding box AND pointing away. As long as this condition is satisfied in one of the dimension, the rest of the straight ray will not intersect the domain.
#for (position < min-tol & direction < 0) or (position > max+tol & direction > 0)
#exit = exit_x | exit_y | exit_z
- class vamtoolbox.projector.pyTorchRayTrace.RayState(device, tensor_dtype=torch.float32, x_0: torch.Tensor = None, v_0: torch.Tensor = None, sino_shape: tuple = None, sino_coord: list = None)#
State of a set of rays. The set can be all rays, rays emitted per angle, per z-slice, or per angle-z-slice This class is used as a dataclass to store current state of a set of rays, and where these rays originally belong to (x,theta,z) or (phi, theta, psi). For parallelization and generality, x_0 and v_0 has shape of (n_rays,n_dim).
- class RaySelector(ray_state)#
- _dtype_to_tensor_type#
- setupRays(ray_trace_ray_config, target_coord_vec_list, azimuthal_angles_deg, inclination_angle_deg=0, ray_density=1, ray_width_fun=None)#
By default, ray_density = 1. It means per azimuthal projection angle, there are max(nX, nY)*nZ number of rays, where nX,nY,nZ = target.shape. When ray_density =/= 1, number of rays = max(nX, nY)*nZ*(ray_density^2) because the density increase/decrease per dimension.
- setupRaysParallel(target_coord_vec_list)#
- resetRaysIterateToInitial()#
- allocateSinogramEnergyToRay(sinogram)#
- allocateRayEnergyToSinogram()#
- reshape(ray_quantity)#
- plot_ray_init_position(angles_deg, color='black')#
- selectCoord(angles, mode='and')#
- selectInverse()#
- _selectSino0(sino0_coord)#
- _selectSino1(sino0_coord)#
- _selectSino2(sino0_coord)#