Commit fbf9b946 authored by Pierre Paleo's avatar Pierre Paleo
Browse files

Merge branch 'half_tomo' into 'master'

Half tomo

See merge request paleo/nabu!5
parents d000fd72 1afaede5
......@@ -10,9 +10,9 @@
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
#
# import os
# import sys
# sys.path.insert(0, os.path.abspath('.'))
import os
import sys
sys.path.insert(0, os.path.abspath('../'))
# -- Project information -----------------------------------------------------
......@@ -32,13 +32,19 @@ release = '2019.1.0'
# ones.
extensions = [
'recommonmark',
'numpydoc',
# 'numpydoc',
'sphinx.ext.napoleon',
'sphinx.ext.autodoc',
'sphinx.ext.mathjax',
'sphinx.ext.viewcode',
'nbsphinx'
'nbsphinx',
# 'sphinx.ext.autosummary',
# 'sphinx.ext.doctest',
# 'sphinx.ext.inheritance_diagram',
]
#autosummary_generate = True
autodoc_member_order = 'bysource'
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
......@@ -73,3 +79,5 @@ def setup(app):
'auto_toc_tree_section': 'Contents',
}, True)
app.add_transform(AutoStructify)
......@@ -11,7 +11,7 @@ A **radiograph** (or **radio**) is the raw data at the output of a X-ray detecto
In this documentation, we distinguish between radios and sinograms.
As stated before, a radio is a raw projection data, i.e what comes out of the detector. "Corrected radios" refer to radios after some CCD corrections and flat-field normalization.
As stated before, a radio is a raw projection data, i.e what comes out of the detector. "Corrected radios" refer to radios after some CCD corrections.
Each radio has *Nx* pixels horizontally, and *Nz* pixels vertically. During a scan, a collections of *Na* radios are acquired, where *Na* is the number of projection angles.
......
......@@ -7,13 +7,20 @@ Welcome to Nabu's documentation!
================================
.. toctree::
:maxdepth: 2
:maxdepth: 1
:caption: Contents:
definitions.md
ccdprocessing.md
nabu_config_file.md
validators.md
modules/index.rst
..
apidoc/nabu.rst
..
Indices and tables
==================
......
......@@ -145,8 +145,8 @@ class MedianFilter(CudaProcessing):
else:
self._allocate_array("d_output", self.shape)
self.cuda_kernel_2d(
self.d_input.gpudata,
self.d_output.gpudata,
self.d_input,
self.d_output,
self.Nx,
self.Ny,
self.Nz,
......
......@@ -71,12 +71,3 @@ __global__ void medfilt2d(
#endif
}
......@@ -286,34 +286,6 @@ class FlatField(CCDProcessing):
return radios_corrected
"""
# It should generate a mask of valid/invalid pixels.
# This does not work yet !
def sigma_clip_mask(data, nsigma_low=4, nsigma_up=None, axis=None, maxiter=5):
NAN = np.float32(np.nan)
nsigma_up = nsigma_up or nsigma_low
res = np.zeros(data.shape, np.float32)
res[:] = data[:]
iteration = 0
outliers = 1
invalids = np.zeros(data.shape, np.bool)
valid_mask = np.ones_like(invalids)
while outliers > 0 and iteration < maxiter:
m = np.nanmean(res, axis=axis)
std = np.nanstd(res, axis=axis)
lower_bound = m - std * nsigma_low
upper_bound = m + std * nsigma_up
invalid_mask = (res[valid_mask] < lower_bound) + (res[valid_mask] > upper_bound)
valid_mask = np.logical_not(invalid_mask)
invalids = invalids + invalid_mask
res[invalid_mask] = NAN
outliers = invalid_mask.sum()
iteration += 1
return invalids
"""
class CCDCorrection(CCDProcessing):
"""
......
......@@ -6,7 +6,6 @@ import pycuda.driver as cuda
from pycuda import gpuarray as garray
from ..utils import updiv, get_cuda_srcfile, _sizeof
#~ from ..cuda.utils import copy_array, get_cuda_context
from ..cuda.processing import CudaProcessing
from ..cuda.kernel import CudaKernel
......
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