Skip to content
Snippets Groups Projects
Commit 212bd598 authored by David Johansen's avatar David Johansen
Browse files

remove old pre rename notebook

parent 8a696428
No related branches found
No related tags found
No related merge requests found
%% Cell type:markdown id: tags:
**Author:** David D. W. Johansen | **Email:** s214743@dtu.dk
%% Cell type:markdown id: tags:
### **CT padding pipeline**
This notebook provides a pipeline for reconstruction using the padding method to mitigate the ROI problem, which can lead to circle and cupping artifacts. An ordinary reconstruction is first done to see whether such artifacts are present and for comparison later in the notebook. There is also export functionality at the end.
%% Cell type:markdown id: tags:
### **Python module imports**
%% Cell type:code id: tags:
```
import glob, os, pathlib
import psutil
import qim3d
import cil
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import clear_output
from ipywidgets import interact, interactive, IntSlider
import math
import h5py
```
%% Cell type:code id: tags:
```
from cil.io import ZEISSDataReader, NikonDataReader
from cil.utilities.display import show_geometry, show2D
from cil.processors import TransmissionAbsorptionConverter
# TIGRE is faster on large-scale CT according to https://tomopedia.github.io/software/tigre/
from cil.recon import FDK # Uses TIGRE
from cil.processors import Slicer
from cil.utilities.jupyter import islicer
from cil.utilities.display import show2D
```
%% Cell type:code id: tags:
```
import viz_utils
```
%% Cell type:code id: tags:
```
import importlib
import sys
importlib.reload(sys.modules['viz_utils'])
```
%% Output
<module 'viz_utils' from '/dtu/3d-imaging-center/projects/2024_QIM_platform/analysis/3dim-reconstruction/viz_utils.py'>
%% Cell type:markdown id: tags:
### **User parameters**
%% Cell type:markdown id: tags:
All the parameters required for the pipeline are set here in the beginning of the notebook. They are also mentioned where they are used. Here they are described in detail:
- `ct_path`: the path to the CT reconstruction metadata file. Should be from **Nikon** (`.xtekct`) or **ZEISS** (`.txrm`).
- `center_height`: determines how many center-most slices to use from the projections. Setting it to `'full'` uses the whole projection volume. `center_height` cannot exceed the height of the projections.
- `pad_factor`: the amount of padding to add at each side of the sinogram in the padding method. Thus setting `pad_factor = 0.25` will make the sinogram 50% bigger. From previous experiments `pad_factor` should be set to at least `0.25` to make the method work properly (pushing the artifacts outside of the ROI). There might be benefits from choosing an even higher `pad_factor`, but this is a case-by-case basis and depends on the specific dataset.
- `save_to_disk`: boolean variable for toggling saving the padded reconstruction. **If `save_to_disk` is set to `True`, then an existing file may be overwritten.**
- `clip_range`: the intensity range that will be used to clip the volume when exporting.
- `base_filename`: the base filename will be used for exporting. The file is saved under the HDF5 format and the extension `.h5` will be appended.
%% Cell type:code id: tags:
```
ct_path = '/dtu/3d-imaging-center/projects/2021_DANFIX_Casper/raw_data_3DIM/Casper_top_3_2 [2021-03-17 16.54.39]/Casper_top_3_2_recon.xtekct'
center_height = 40
pad_factor = 0.25
save_to_disk = True
clip_range = (0.0, 0.07)
base_filename = 'top_3_2_padded'
```
%% Cell type:markdown id: tags:
### **Data reading and processing**
%% Cell type:markdown id: tags:
For faster processing, we may work on a subset of slices. `center_height` determines how many of the center-most slices should be used. Set to `'full'` to use full volume.
%% Cell type:code id: tags:
```
def create_reader(file_name, roi=None):
# Wrapper for
if file_name.endswith('txrm'):
DataReader = ZEISSDataReader
elif file_name.endswith('xtekct'):
DataReader = NikonDataReader
else:
raise ValueError("Unrecognizable CT metadata file. File extension should either be '.txrm' or '.xtekct'")
if roi is None:
return DataReader(file_name=file_name)
else:
return DataReader(file_name=file_name, roi=roi)
def get_pixel_nums(ct_path):
reader = create_reader(file_name=ct_path)
num_pixels_h = reader.get_geometry().pixel_num_h
num_pixels_v = reader.get_geometry().pixel_num_v
return num_pixels_h, num_pixels_v
num_pixels_h, num_pixels_v = get_pixel_nums(ct_path)
if center_height == 'full':
reader = create_reader(file_name=ct_path)
else:
slice_dict = {'vertical': (
num_pixels_v // 2 - center_height // 2,
num_pixels_v // 2 + center_height // 2,
1
)}
reader = create_reader(file_name=ct_path, roi=slice_dict)
```
%% Cell type:markdown id: tags:
Reading the data might take some time.
%% Cell type:code id: tags:
```
data = reader.read()
```
%% Cell type:markdown id: tags:
Usually the data is given in the transmission domain and thus we convert to the absorption domain.
%% Cell type:code id: tags:
```
data = TransmissionAbsorptionConverter()(data)
```
%% Cell type:markdown id: tags:
### **Ordinary reconstruction**
%% Cell type:code id: tags:
```
data.reorder(order='tigre')
recon = FDK(data).run(verbose=0)
```
%% Cell type:markdown id: tags:
Here we inspect the reconstruction for artifacts to see if padding is necessary.
%% Cell type:code id: tags:
```
# islicer(recon)
qim3d.viz.line_profile(recon.as_array())
```
%% Output
VBox(children=(HBox(children=(VBox(children=(HTML(value="<div style='text-align:center; font-size:16px; font-w…
%% Cell type:markdown id: tags:
### **Reconstruction with padding**
%% Cell type:code id: tags:
```
from cil.processors import Padder
dim_horizontal = data.get_data_axes_order().index('horizontal')
pad_width = round(pad_factor * data.shape[dim_horizontal])
data_padded = Padder.edge(pad_width={'horizontal': pad_width})(data)
```
%% Cell type:code id: tags:
```
data.shape, data_padded.shape
```
%% Output
((1571, 40, 1000), (1571, 40, 1500))
%% Cell type:markdown id: tags:
In the following, the padded data is passed as a sinogram to FDK, while the original image geometry is used for the volume to reconstruct. This is more efficient than the default behaviour of reconstructing on a correspondingly padded volume that we would have to crop to the original size anyway.
%% Cell type:code id: tags:
```
data_padded.reorder(order='tigre')
ig = data.geometry.get_ImageGeometry()
recon_padded = FDK(data_padded, ig).run(verbose=0)
```
%% Cell type:code id: tags:
```
print(recon_padded)
```
%% Output
Number of dimensions: 3
Shape: (40, 1000, 1000)
Axis labels: ('vertical', 'horizontal_y', 'horizontal_x')
%% Cell type:code id: tags:
```
# islicer(recon_padded, axis_labels=list(recon.dimension_labels))
qim3d.viz.line_profile(recon_padded.as_array())
```
%% Output
VBox(children=(HBox(children=(VBox(children=(HTML(value="<div style='text-align:center; font-size:16px; font-w…
%% Cell type:markdown id: tags:
### **Visual comparison**
%% Cell type:markdown id: tags:
The parameter `vrange` specifies the ranges to use. Not specifying it or setting it to `None` uses each of the volumes full intensity range. Experimenting with `vrange` can be used to finetune `clip_range`.
%% Cell type:code id: tags:
```
vrange = clip_range
viz_utils.comparator_widget(recon.as_array(), recon_padded.as_array(), titles=['No padding', 'Padding'], vrange=vrange)
```
%% Output
interactive(children=(IntSlider(value=0, description='Axis', max=2), IntSlider(value=20, description='Index', …
%% Cell type:markdown id: tags:
### **Exporting**
%% Cell type:markdown id: tags:
Choose the interval `clip_range` to clip the data to and the `base_filename` to save to (the extension is added automatically). The `clip_range` and `pad_factor` will be stored as metadata. The reconstruction volume will be converted to uint16 for exporting.
%% Cell type:code id: tags:
```
h5_filename = f'{base_filename}.h5'
```
%% Cell type:code id: tags:
```
if save_to_disk:
vol = np.clip(recon_padded.as_array(), a_min=clip_range[0], a_max=clip_range[1])
vol = (vol - clip_range[0]) / (clip_range[1] - clip_range[0])
vol = (vol * (2**16 - 1)).astype(np.uint16)
with h5py.File(h5_filename, 'w') as f:
f.create_dataset('recon_vol', data=vol)
f.attrs['clip_range'] = clip_range
f.attrs['pad_factor'] = pad_factor
```
%% Cell type:markdown id: tags:
Below we show how to read the exported file back into memory.
%% Cell type:code id: tags:
```
with h5py.File(h5_filename, 'r') as f:
vol = f['recon_vol'][...]
```
%% Cell type:code id: tags:
```
qim3d.viz.line_profile(vol)
```
%% Output
VBox(children=(HBox(children=(VBox(children=(HTML(value="<div style='text-align:center; font-size:16px; font-w…
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment