Commit 46d7b6d9 authored by Pierre Paleo's avatar Pierre Paleo
Browse files

Start SinoProcessing

parent d000fd72
import numpy as np
class SinoProcessing(object):
"""
A base class for processing sinograms.
"""
def __init__(self, shape=None, radios=None):
"""
Initialize a SinoProcessing instance.
Parameters
----------
shape: tuple of int
Shape of the sinogram or stack of sinograms, in the form
(n_z, n_angles, n_x).
Note: in the radio domain, the shape is (n_angles, n_z, n_x).
If the tuple has two elements, it is assumed to be only one sinogram,
i.e the tuple is parsed as (n_angles, n_x).
radios: array-like
Initialize the SinoProcessing instance from a chunk of radios.
In this case, the parameter "shape" is ignored.
"""
if (shape is None) ^ (radios is None) == 0:
raise ValueError("Please provide either shape= or radios=")
self.radios = radios
if shape is not None:
if len(shape) == 2:
n_z = 1
n_angles, n_x = shape
elif len(shape) == 3:
n_z, n_angles, n_x = shape
else:
raise ValueError("Expected a 2-tuple or 3-tuple for shape")
if radios is not None:
if radios.ndim != 3:
raise ValueError("Expected a 3D array of radios")
n_angles, n_z, n_x = radios.shape
self.shape = (n_z, n_angles, n_x)
self.n_angles = n_angles
self.n_x = n_x
self.n_z = n_z
self.sinos = None
def radios_to_sinos(self, radios=None, output=None):
"""
Transpose a chunk of radios to a stack of sinograms.
Parameters
-----------
radios: array, optional
Radios to transpose. If not provided, the radios provided at class
instantiation are taken.
output: array, optional
Output sinograms array, pre-allocated
"""
if radios is None:
if self.radios is None:
raise ValueError("Please provide a radios array")
radios = self.radios
# todo check radios shape
sinos = np.rollaxis(radios, 1, 0) # view
if output is None:
res = np.ascontiguousarray(sinos) # copy
else:
# todo check output shape
for i in range(output.shape[0]):
output[i] = sinos[i]
res = output
self.sinos = res
return res
def _set_sinos(self, sino):
if sino is None:
if self.sinos is None:
raise ValueError("Need to provide sino array")
sino = self.sinos
if sino.ndim == 2:
sino = sino.reshape((1,) + sino.shape)
return sino
# TODO support drift in rotation axis
def convert_to_halftomo(self, rot_center, sinos=None, output=None):
"""
Convert a (stack of) sinogram acquired over 360 degrees, to a sinogram
with extended FOV with the "half tomography" setting.
Parameters
----------
rot_center: float or array
Position of the rotation axis. A single float indicates that the rotation
axis position is the same for all projections.
An array indicates that the axis position varies as a function of the
projections.
sinos: array, optional
Sinogram or stack of sinograms.
output: array, optional
Output
"""
if not (np.isscalar(rot_center)):
raise ValueError("rot_center = <array> is not supported yet")
sinos = self._set_sinos(sinos)
n_z, n_a, n_x = sinos.shape
if (n_a % 2) != 0:
raise ValueError("Expected an even number of angles")
n_a2 = n_a // 2
r = int(round(rot_center))
out_dwidth = 2 * r
if output is not None:
if output.shape[-1] != out_dwidth:
raise ValueError(
"Expected output sinogram last dimension to have %d elements"
% out_dwidth
)
if output.shape[-2] != n_a2:
raise ValueError("Expected output sinogram to have %d angles" % n_a2)
if output.ndim == 2:
output = output.reshape((1,) + output.shape)
else:
output = np.zeros((n_z, n_a2, out_dwidth), dtype=np.float32)
d = n_x - r
for i in range(n_z):
output[i][:] = convert_halftomo(sinos[i], r)
return output
def convert_halftomo(sino, rotation_axis_position):
"""
Converts a sinogram into a sinogram with extended FOV with the "half tomography"
setting.
"""
assert sino.ndim == 2
assert (sino.shape[0] % 2) == 0
na, nx = sino.shape
na2 = na // 2
r = rotation_axis_position
d = nx - r
res = np.zeros((na2, 2 * r), dtype="f")
sino1 = sino[:na2, :]
sino2 = sino[na2:, ::-1]
res[:, : nx - d] = sino1[:, : nx - d]
#
w1 = np.linspace(0, 1, d, endpoint=True)
res[:, nx - d : nx] = (1 - w1) * sino1[:, nx - d :] + w1 * sino2[:, d : 2 * d]
#
res[:, nx:] = sino2[:, 2 * d :]
return res
import numpy as np
import pycuda.gpuarray as garray
from .sinogram import SinoProcessing
class CudaSinoProcessing(SinoProcessing):
def __init__(self, shape=None, radios=None):
super().__init__(shape=shape, radios=radios)
def radios_to_sinos(self, radios=None, output=None):
if radios is None:
if self.radios is None:
raise ValueError("Please provide a radios array")
radios = self.radios
# todo check radios shape
if output is None:
output = garray.zeros(self.shape, np.float32)
else:
# todo check output shape
pass
for i in range(output.shape[0]):
output[i][:] = radios[:, i, :] # memcpy D2D
return output
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment