Reconstruction Clipart to Build Again Clipart
This is part two of a tutorial that describes how to use the ASTRA Toolbox to create a 3D reconstruction from second projection images that were taken with a cone-axle CT scanner. In part ane, I've created a synthetic dataset. In this office, I'll use that dataset to create a reconstruction. Note that, if you start from a real dataset instead of from this simulated i, you might nevertheless have to convert your 10-ray images into projection images. This is explained in part 1 of my previous series of manufactures on tomography. However, apart from that, these projection images are completely equivalent to real ones, admitting that they look every bit if taken with a mechanically perfect (but non noiseless!) scanner.
With the take chances of repeating a few things from role one, I'll list all steps to take to go from the projection images to the concluding reconstruction.
Scan Parameters
The first thing to do is to assemble all relevant information about the scan. This data will be in some sort of logfile that comes with the scan, since the parameters can be adjusted for each scan, even if the same scanner is used. For example, the position of the rotating phase with the sample can exist shifted closer to or farther abroad from the detector, to modify the magnification. For the purpose of this article, I'll use the Python fragment from part one of this tutorial, which is repeated beneath.
distance_source_origin = 300 # [mm] distance_origin_detector = 100 # [mm] detector_pixel_size = one.05 # [mm] detector_rows = 200 # Vertical size of detector [pixels]. detector_cols = 200 # Horizontal size of detector [pixels]. num_of_projections = 180 angles = np.linspace ( 0 , 2 * np.pi , num=num_of_projections, endpoint= Faux )
This "logfile" provides enough data to be able to create a reconstruction from the project images. Of course, you also take to be sure that the projections were actually taken with a cone-beam scanner.
The origin of the coordinate system of the scanner is the center of the object (assuming that the object is in the middle of the rotating stage, of course). From the Python code, you can run across that the altitude betwixt the source and the origin, allow's phone call that \(d_\mathrm{so}\), is 300 mm, and that the distance between the origin and the detector, \(d_\mathrm{od}\), is 100 mm. This means that the magnification is \((d_\mathrm{and so}+d_\mathrm{sd})/d_\mathrm{so}=4/three\) in this instance. The detector has square pixels with size \(s_\mathrm{det}=i.05\,\mathrm{mm}\), and has 200 rows in the vertical direction and 200 columns in the horizontal direction. The angles are equally spaced in 180 steps over 360 degrees, which means that a project was made every two degrees.
You should be able to adapt these values by looking at the logfiles of your own scans. Sometimes the distances are given as source–origin and source–detector, so you might to have to adapt for that. In the Python code that follows, I assume that your projections are in a directory called dataset and that your files are named proj####.tif .
Loading the Dataset
Loading the dataset is done in two steps. The get-go step is reading the dataset from disk and putting information technology in a 3D NumPy array, equally follows.
projections = np.zeros ( (detector_rows, num_of_projections, detector_cols) ) for i in range (num_of_projections): im = imread(bring together(input_dir, 'proj%04d.tif' % i) ).astype ( float ) im /= 65535 projections[:, i, :] = im
This is straighforward, just notation that information technology is axis one that increases with each projection. The code assumes 16-scrap TIFF files. This is quite common, just you might have to accommodate this if your projection images are in a different format. I assume hither that the projections use the full range of 0 to 65535. This wil normally never exist the case with actual data, so y'all volition have to rescale your data. The actual maximum is not important, because we are working with relative numbers here. I have rescaled for a maximum of one in this example. The minimum, however, should be naught later you have converted the X-ray images into projection images (over again, see my existing series of articles on tomography). This is because, in places where there is no object in between the source and the detector, the measured density of the (non-existant) material should be exactly nada.
Creating the Projection Geometry
The 2d step is then to put the NumPy variable projections into the toolbox past using the function astra.data3d.create() . To do that, you lot beginning have to depict the geometry of your projection data by creating a projection geometry through calling astra.create_proj_geom() .
The cosmos of the projection geometry is ane of the crucial steps that you need to understand well to successfully utilize the ASTRA Toolbox. There is ane of import holding of the geometries of the ASTRA Toolbox, which is that the voxels of the reconstruction area have a size of 1 (so 1×1×1 in the case of the 3D reconstruction that we are working on). This implies that you have to scale the rest of your parameters to make this so. A little "play a joke on" to make this easier to handle is to put the detector about in the center of the object. This is not mandatory, merely it simplifies thinking about information technology. Or, at least, that's what I think.
The parameters of the scanner as given in the Python lawmaking fragment to a higher place are the following.
\[\begin{align}
d_\mathrm{so}&=300\,\mathrm{mm}\\
d_\mathrm{od}&=100\,\mathrm{mm}\\
s_\mathrm{det}&=i.05\,\mathrm{mm}
\stop{align}\]
Shifting the detector to the origin amounts to setting the origin–detector distance to cipher and rescaling the size of the detector pixels to compensate for that. This results in a new fix of parameters:
\[\brainstorm{marshal}
d_\mathrm{so}'&=d_\mathrm{then}=300\,\mathrm{mm}\\
d_\mathrm{od}'&=0\,\mathrm{mm}\\
s_\mathrm{det}'&=\frac{d_\mathrm{and so}}{d_\mathrm{so}+d_\mathrm{od}}s_\mathrm{det}\approx 0.79\,\mathrm{mm}
\stop{align}\]
Afterward shifting the detector to the origin, the parameter \(s_\mathrm{det}'\) gives you the bodily size of the reconstruction voxels in real-globe units. In the instance scanner, each voxel in the reconstructed object volition correspond to a cube of 0.79×0.79×0.79 mm. This is of import information if y'all need to know or measure the actual dimensions of the scanned object in real-world units.
The next pace is then to scale the parameters so that the reconstruction voxels get a size of i×1×1. This results in the terminal gear up of parameters:
\[\begin{marshal}
d_\mathrm{so}''&=\frac{d_\mathrm{so}'}{s_\mathrm{det}'}=\frac{d_\mathrm{so}+d_\mathrm{od}}{s_\mathrm{det}}\approx 381\\
d_\mathrm{od}''&=0\\
s_\mathrm{det}''&=one
\end{align}\]
I've left out the units in this terminal set of expressions, since it seems confusing to keep including them after this scaling. In Python, this results in the post-obit code.
proj_geom = \ astra.create_proj_geom ( 'cone' , i , one , detector_rows, detector_cols, angles, (distance_source_origin + distance_origin_detector) / detector_pixel_size, 0 ) projections_id = astra.data3d.create ( '-sino' , proj_geom, projections)
The result of this footstep is that yous now have a projections_id handle to the projection data in the Toolbox that you need to provide to subsequent calls to the Toolbox.
Reconstruction
The reconstruction itself is commonly not the hard part. You create a volume geometry and then create a 3D information object inside of the Toolbox to hold the reconstruction. This provides you lot with a reconstruction_id handle. You then choose an algorithm and specify its parameters. For a cone-beam dataset, the FDK_CUDA algorithm is the obvious i to start with. The resulting Python code is shown beneath.
vol_geom = astra.creators.create_vol_geom (detector_cols, detector_cols, detector_rows) reconstruction_id = astra.data3d.create ( '-vol' , vol_geom, data= 0 ) alg_cfg = astra.astra_dict ( 'FDK_CUDA' ) alg_cfg[ 'ProjectionDataId' ] = projections_id alg_cfg[ 'ReconstructionDataId' ] = reconstruction_id algorithm_id = astra.algorithm.create (alg_cfg) astra.algorithm.run (algorithm_id) reconstruction = astra.data3d.get (reconstruction_id)
Afterwards this role of the code has run, the NumPy array reconstruction contains the 3D reconstructed object. Yous could then cutting off negative values, which should be due to dissonance if the scaling of your dataset was done correctly. To salvage the reconstruction as a series of individual slices, you lot tin can convert the reconstruction book to an 8-bit (or 16-bit) unsigned int and salve the slices equally PNG or TIFF files:
reconstruction[reconstruction < 0 ] = 0 reconstruction /= np.max (reconstruction) reconstruction = np.circular (reconstruction * 255 ).astype (np.uint8 ) for i in range (detector_rows): im = reconstruction[i, :, :] im = np.flipud (im) imwrite(join(output_dir, 'reco%04d.png' % i) , im)
Note the essential call to flipud() . This necessary because of the inverted Y-axis that you get when you convert from a NumPy array to an image. The resulting PNG for the and then-chosen central slice, which is piece 100 in this case, is shown in Figure 1.
Figure ane. Central piece of the reconstructed object.
I hope that this extensive tutorial commodity helps you in applying the ASTRA Toolbox to your own cone-axle CT datasets!
Consummate Python Plan
The complete Python program follows beneath. Note that yous need an NVIDIA GPU to run this code, since the 3D algorithms of the ASTRA Toolbox are implemented equally CUDA GPU code.
from __future__ import partitioning import numpy as np from os import mkdir from os.path import join, isdir from imageio import imread, imwrite import astra # Configuration. distance_source_origin = 300 # [mm] distance_origin_detector = 100 # [mm] detector_pixel_size = one.05 # [mm] detector_rows = 200 # Vertical size of detector [pixels]. detector_cols = 200 # Horizontal size of detector [pixels]. num_of_projections = 180 angles = np.linspace ( 0 , two * np.pi , num=num_of_projections, endpoint= False ) input_dir = 'dataset' output_dir = 'reconstruction' # Load projections. projections = np.zeros ( (detector_rows, num_of_projections, detector_cols) ) for i in range (num_of_projections): im = imread(bring together(input_dir, 'proj%04d.tif' % i) ).astype ( float ) im /= 65535 projections[:, i, :] = im # Copy projection images into ASTRA Toolbox. proj_geom = \ astra.create_proj_geom ( 'cone' , 1 , ane , detector_rows, detector_cols, angles, (distance_source_origin + distance_origin_detector) / detector_pixel_size, 0 ) projections_id = astra.data3d.create ( '-sino' , proj_geom, projections) # Create reconstruction. vol_geom = astra.creators.create_vol_geom (detector_cols, detector_cols, detector_rows) reconstruction_id = astra.data3d.create ( '-vol' , vol_geom, data= 0 ) alg_cfg = astra.astra_dict ( 'FDK_CUDA' ) alg_cfg[ 'ProjectionDataId' ] = projections_id alg_cfg[ 'ReconstructionDataId' ] = reconstruction_id algorithm_id = astra.algorithm.create (alg_cfg) astra.algorithm.run (algorithm_id) reconstruction = astra.data3d.go (reconstruction_id) # Limit and scale reconstruction. reconstruction[reconstruction < 0 ] = 0 reconstruction /= np.max (reconstruction) reconstruction = np.round (reconstruction * 255 ).astype (np.uint8 ) # Salvage reconstruction. if not isdir(output_dir): mkdir(output_dir) for i in range (detector_rows): im = reconstruction[i, :, :] im = np.flipud (im) imwrite(join(output_dir, 'reco%04d.png' % i) , im) # Cleanup. astra.algorithm.delete (algorithm_id) astra.data3d.delete (reconstruction_id) astra.data3d.delete (projections_id)
Source: https://tomroelandts.com/articles/astra-toolbox-tutorial-reconstruction-from-projection-images-part-2
0 Response to "Reconstruction Clipart to Build Again Clipart"
Post a Comment