The first step is to know how to load data to workspace. Different
communities use different file formats: hdf/nxs, mrc, txrm, xrm, tif, dicom,
raw,… and they need specific Python libraries to work with. The common format
used by synchrotron/neutron community is hdf/nxs. To load a dataset from a hdf/nxs
file we need to know the key or path to the data. This can be done using Hdfviewer
or Algotom’s function as shown in section 1.2. For a tomographic
hdf file, some basic metadata we need to know:
Keys to projections images, flat-field images, and dark-field images.
Key to rotation angles corresponding to projection images. This information may be
not needed if the data was acquired in the ange range of [0; 180-degree].
If the data is from a helical scan, extra information such as pixel size,
pitch, and translation positions is needed.
Information which is used by specific data processing methods such as pixel
size, sample-detector distance, or X-ray energy.
Fig. 1.4.1 shows the content of a hdf file where users can find key
to datasets used for tomographic reconstruction.
importalgotom.io.loadersaveraslosafile="C:/data/tomo_00064.h5"proj_img=losa.load_hdf(file,key_path="exchange/data")# This is an hdf object, no data being loaded yet.flat_img=losa.load_hdf(file,key_path="exchange/data_white")dark_img=losa.load_hdf(file,key_path="exchange/data_dark")angles=losa.load_hdf(file,key_path="exchange/theta")
Another tomographic hdf/nxs file acquired at beamline I12, Diamond Light Source (DLS) where
the file structure and keys are different to the one before. In this data,
dark-field images and flat-field images are at the same dataset as projection images
where there is a dataset named “image_key” used to distinguish these images.
Fig. 1.4.2 Datasets of a tomographic hdf file acquired at DLS. Image-key value of 2 is
for dark-field, 1 is for flat-field, and 0 is for projection image.¶
We can load the data, extract a projection image, and save it to tif.
importalgotom.io.loadersaveraslosafile="E:/Tomo_data/68067.nxs"data_img=losa.load_hdf(file,key_path="entry1/tomo_entry/data/data")# This is an hdf object.# Extract a subset of data and save to tifprint(data_img.shape)# 1861 x 2160 x2560losa.save_image("E:/output/image_00061.tif",data_img[61])
There are many Algotom’s functions in the IO module
to handle different tasks such as converting tif images to hdf, displaying the
hdf tree, or saving output to a hdf file.
The flat-field correction process is based on the Beer-Lambert’s law
\[\frac{I}{I_0} = \int_{}e^{-\alpha (x,y,z) dx}\]
in practice, it is done using the following formula
\[\frac{P_{\theta}-D}{F-D}\]
where \(P_{\theta}\) is a projection image of a sample at a rotation
angle of \(\theta\), \(D\) is a dark-field image (camera’s dark noise)
taken with a photon source off, and \(F\) is a flat-field image taken without
the sample. This can be done using Algotom as follows; data used in this
demonstration can be download from here
importnumpyasnpimportalgotom.io.loadersaveraslosafile="E:/Tomo_data/68067.nxs"data_img=losa.load_hdf(file,key_path="entry1/tomo_entry/data/data")# This is an hdf object.# Get image keyikey=losa.load_hdf(file,key_path="entry1/tomo_entry/instrument/detector/image_key")ikey=np.squeeze(np.asarray(ikey[:]))# Load data and convert to numpy 1d-array.# Use image_key to load flat-field images and average themdark_field=np.mean(np.asarray(data_img[np.squeeze(np.where(ikey==2.0)),:,:]),axis=0)flat_field=np.mean(np.asarray(data_img[np.squeeze(np.where(ikey==1.0)),:,:]),axis=0)# Get indices of projection imagesproj_idx=np.squeeze(np.where(ikey==0))# Apply flat-field correction to the first projection image.proj_img=data_img[proj_idx[0]]flat_dark=flat_field-dark_fieldnmean=np.mean(flat_dark)flat_dark[flat_dark==0.0]=nmean# Handle zero divisionproj_norm=(proj_img-dark_field)/flat_dark# Save imageslosa.save_image("E:/output/proj_before.tif",proj_img)losa.save_image("E:/output/proj_after.tif",proj_norm)
Running the code gives the output images
We can apply the process to a sinogram.
# Generate sinogram at the middle of an image height(depth,height,width)=data_img.shapesino_idx=height//2start=proj_idx[0]stop=proj_idx[-1]+1sinogram=data_img[start:stop,sino_idx,:]# Apply flat-field correction the sinogramsino_norm=(sinogram-dark_field[sino_idx])/flat_dark[sino_idx]# Save imageslosa.save_image("E:/output/sino_before.tif",sinogram)losa.save_image("E:/output/sino_after.tif",sino_norm)
Zingers are prominent bright dots in images caused by scattered X-rays hitting
the detector system CCD or CMOS chip (Fig. 1.4.3 (a,b)). They produce
line artifacts across a reconstructed image (Fig. 1.4.3 (c)).
Fig. 1.4.3 Artifacts caused by zingers. (a) Zingers in the sinogram space. (b) Zingers
in the projection space. (c) Line artifacts caused by the zingers.¶
Zingers are easily removed by using a method
in Algotom
Causes of ring artifacts and methods for removing them [R19] have been documented in detailed
here. There are many methods to choose from
in Algotom. However the combination of methods
has been proven to be the most effective way to clean most of ring artifact types.
Note that in the sinogram space, ring artifacts appear as stripe artifacts. Example
of how to use the methods
There is a method for enhancing the contrast of an image, known as the Paganin filter
which is commonly used at synchrotron facilities. Algotom implements a
simplified version of this filter, named
the Fresnel filter as it is based on the Fresnel propagator. There is
a widespread misunderstanding in the community that the resulting image of
the Paganin filter is a phase-contrast image. It is not. Because the filter acts
as a low-pass filter, it reduces noise and the dynamic range of an image. This helps
to enhance the contrast between low-contrast features which can be confused if
this enhancement comes from the phase effect. Detailed demonstration for the
argument is at here.
Note that a denoising filter or smoothing filter should not be used before the above
pre-processing methods (zinger removal, ring artifact removal, center calculation).
Blurring an image will impact the performance of these methods.
There are many choices for reconstruction methods and open-source software. In
the current version (<=1.1), Algotom implements two FFT-based methods which is fast enough
for a 2k x 2k x 2k dataset. Algotom also provides wrappers for other reconstruction
methods available in Tomopy (gridrec)
and Astra Toolbox (FBP, SIRT, SART, CGLS,…).
Examples of comparing reconstructed images before and after artifacts removal:
importalgotom.rec.reconstructionasrec# No need to pass angles if it's a 180-degree sinogramrec_img1=rec.dfi_reconstruction(sino_norm,center,angles=None)rec_img2=rec.dfi_reconstruction(sino_rem2,center,angles=None)losa.save_image("E:/output/rec_with_artifacts.tif",rec_img1)losa.save_image("E:/output/rec_artifacts_removed.tif",rec_img2)
Examples of comparing reconstructed images after using the Fresnel filter with
different strengths:
Calculate distortion coefficients [R18] using the Discorpy package. The output is a text file.
Use the calculated coefficients for correction.
importnumpyasnpimportalgotom.io.loadersaveraslosaimportalgotom.prep.correctionascorrimportalgotom.prep.removalasremoimportalgotom.prep.calculationascalcimportalgotom.prep.filteringasfiltimportalgotom.rec.reconstructionasreco# Paths to data. Download at: https://doi.org/10.5281/zenodo.3339629proj_path="E:/data/tomographic_projections.hdf"flat_path="E:/data/flats.hdf"dark_path="E:/data/darks.hdf"coef_path="E:/data/coefficients_bw.txt"key_path="/entry/data/data"# Where to save the outputsoutput_base="E:/output/"# Load data of projection images as an hdf objectproj_data=losa.load_hdf(proj_path,key_path)(depth,height,width)=proj_data.shape# Load flat-field images and dark-field images, average each of themflat_field=np.mean(losa.load_hdf(flat_path,key_path)[:],axis=0)dark_field=np.mean(losa.load_hdf(dark_path,key_path)[:],axis=0)# Load distortion coefficientsxcenter,ycenter,list_fact=losa.load_distortion_coefficient(coef_path)# Apply distortion correction to dark- and flat-field image.flat_discor=corr.unwarp_projection(flat_field,xcenter,ycenter,list_fact)dark_discor=corr.unwarp_projection(dark_field,xcenter,ycenter,list_fact)# Generate a sinogram with distortion correction.index=800sinogram=corr.unwarp_sinogram(proj_data,index,xcenter,ycenter,list_fact)sinogram=corr.flat_field_correction(sinogram,flat_discor[index],dark_discor[index])sinogram=remo.remove_all_stripe(sinogram,3.0,51,17)center=calc.find_center_vo(sinogram,width//2-50,width//2+50)# Reconstruct image from the sinogramrec_img=reco.dfi_reconstruction(sinogram,center,angles=None,apply_log=True)losa.save_image(output_base+"/rec_00800.tif",rec_img)
1.4.8.2. Sinogram stitching for a half-acquisition scan¶
Half-acquisition scanning technique are being used more often at synchrotron facilities.
It is a simple technique to double the field-of-view (FOV) of a tomography system
by shifting the rotation axis to a side of the FOV then acquiring data in the
angle range of [0, 360-degree]. To process the data, a 360-degree sinogram is
converted to an equivalent 180-degree sinogram by stitching two halves of
the 360-degree sinogram, before reconstruction. For stitching, we need to know
either the center-of-rotation, or the overlap-area and overlap-side between
two halves of the sinogram. Algotom provides methods[C1] for automatically finding these parameters.
importnumpyasnpimportalgotom.io.loadersaveraslosaimportalgotom.prep.correctionascorrimportalgotom.prep.removalasremoimportalgotom.prep.calculationascalcimportalgotom.prep.conversionasconvimportalgotom.rec.reconstructionasrecoinput_base="E:/data/"output_base="E:/output/"# Data at: https://doi.org/10.5281/zenodo.4386983proj_path=input_base+"/scan_00008/projections_00000.hdf"flat_path=input_base+"/scan_00009/flats_00000.hdf"dark_path=input_base+"/scan_00009/darks_00000.hdf"meta_path=input_base+"/scan_00008/scan_00008.nxs"key_path="/entry/data/data"angle_key="/entry1/tomo_entry/data/rotation_angle"data=losa.load_hdf(proj_path,key_path)(depth,height,width)=data.shapeangles=np.squeeze(np.asarray(losa.load_hdf(meta_path,angle_key)[:]))# Load dark-field images and flat-field images, averaging each result.flat_field=np.mean(losa.load_hdf(flat_path,key_path)[:],axis=0)dark_field=np.mean(losa.load_hdf(dark_path,key_path)[:],axis=0)# Generate a sinogram and perform flat-field correction.index=height//2sino_360=corr.flat_field_correction(data[:,index,:],flat_field[index],dark_field[index])# Calculate the center-of-rotation, the overlap-side and overlap-area used for stitching(center0,overlap,side,_)=calc.find_center_360(sino_360,100)# Remove zingerssino_360=remo.remove_zinger(sino_360,0.08)# Remove ring artifactssino_360=remo.remove_all_stripe(sino_360,3,51,17)# Convert the 360-degree sinogram to the 180-degree sinogram.sino_180,center1=conv.convert_sinogram_360_to_180(sino_360,center0)losa.save_image(output_base+"/sino_360.tif",sino_360)losa.save_image(output_base+"/sino_180.tif",sino_180)# Perform reconstructionrec_img=reco.dfi_reconstruction(sino_180,center1,apply_log=True)losa.save_image(output_base+"/rec_img_1.tif",rec_img)# 2nd way: extend the 360-degree sinogram. It's useful for tomography fly-scans# where the two halves of a 360-degree sinogram are mismatch due to the angle# step is not divisible.(sino_ext,center2)=conv.extend_sinogram(sino_360,center0)# Perform reconstruction# Using fbp-method for angle range > 180 degreeimg_rec=reco.fbp_reconstruction(sino_ext,center2,angles=angles*np.pi/180.0,apply_log=False,gpu=True)losa.save_image(output_base+"/rec_img_2.tif",img_rec)
