Radiography registration for mosaic tomography

Vescovi, R.F.C.; Cardoso, M.B.; Miqueles, E.X.

doi:10.1107/S1600577517001953

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 24| Part 3| May 2017| Pages 686-694

https://doi.org/10.1107/S1600577517001953

Radiography registration for mosaic tomography

R. F. C. Vescovi,^a,^b M. B. Cardoso ^b and E. X. Miqueles ^b ^*

^aInstitute Física Gleb-Wataghin (IFGW), University of Campinas, SP, Brazil, and ^bBrazilian Synchrotron Light Laboratory/CNPEM, Campinas, SP, Brazil
^*Correspondence e-mail: eduardo.miqueles@lnls.br

Edited by P. A. Pianetta, SLAC National Accelerator Laboratory, USA (Received 15 February 2016; accepted 6 February 2017; online 7 April 2017)

A hybrid method of stitching X-ray computed tomography (CT) datasets is proposed and the feasibility to apply the scheme in a synchrotron tomography beamline with micrometre resolution is shown. The proposed method enables the field of view of the system to be extended while spatial resolution and experimental setup remain unchanged. The approach relies on taking full tomographic datasets at different positions in a mosaic array and registering the frames using Fourier phase correlation and a residue-based correlation. To ensure correlation correctness, the limits for the shifts are determined from the experimental motor position readouts. The masked correlation image is then minimized to obtain the correct shift. The partial datasets are blended in the sinogram space to be compatible with common CT reconstructors. The feasibility to use the algorithm to blend the partial datasets in projection space is also shown, creating a new single dataset, and standard reconstruction algorithms are used to restore high-resolution slices even with a small number of projections.

Keywords: computed tomography; mosaic tomography; microtomography; non-limited field of view.

1. Introduction

Since the introduction of imaging detectors, from everyday cellphone pictures to high-resolution satellite topographies, many applications use information panoramas to increase the field of view (FOV) of their equipment. This is especially true when it is needed to image large objects compared with the experiment resolution (Miller, 2006 ; Ma et al., 2007 ).

One field of study that relies on area detectors is X-ray imaging. In conventional X-ray images the apparent pixel size is directly proportional to the FOV and the geometry of the experiment. Contrast is generated by X-rays propagating through the sample, and the intensity measured at the detector is proportional to the integral of the sample refraction index (Als-Nielsen & McMorrow, 2011 ). This makes X-ray images suitable for volumetric reconstruction through an inverse transform. The solution is exactly given by the Radon transform and has been discussed extensively by Deans (2007 ).

A technique that can acquire a dataset for the inverse Radon reconstruction is called computed tomography (CT) and relies on acquiring sample projections at different angles. One of the main difficulties with the FOV limitation is that conventional CT reconstruction algorithms rely on having datasets with finite support (equivalent to having the whole sample contained within the FOV in every projection). This may not be a problem when the sample is sufficiently small or can be cut to fit the FOV. However, in many cases the whole sample presents useful information but cannot be physically damaged. Solving this problem can be difficult since increasing the size of the initial data makes the inversion more complex and inaccurate (Hansen, 1992 ).

Many methods have been proposed to extend the FOV in tomography without losing pixel resolution. A comparison of several methods has been given by Kyrieleis et al. (2009 ). Every variation of extended FOV methods relies on accurate knowledge of the sample position in every frame of the complete dataset; this can be calculated after the experiment or ensured in the data acquisition.

The simplest approach is to move the sample along the rotation axis. Since the reconstruction is done along the axis perpendicular to the rotation, this approach does not interfere with the final reconstruction quality. Nevertheless, this solution requires a stretched sample form. Since the sample translation does not affect the reconstruction algorithm, one could find the overlap area after reconstruction.

A more complicated approach is moving the rotation axis in relation to the camera axis. This translation perpendicular to the camera makes the reconstruction sensitive to the accuracy of the registration and the quality of frame merging. With motor control, one can have a rough approximation of the shift between two images. If the final image pixel size is of the order of the translation motor resolution, this approximation is satisfactory and experiment misalignments will not interfere. When pixel size is smaller than motor resolution it is necessary to correlate the dataset's position in the mosaic grid to reconstruct the image.

Image panorama relies on finding the relative geometric transformation (shift, rotation and magnification) between two images and using the result to merge the images into a new image. Other artifacts such as different illumination and object motion cannot be easily corrected, as discussed by Brown (1992 ). Such a procedure is often called image registration.

Although the process can be simple for extending the FOV of a single radiography, our approach proposes to merge whole partial tomogram datasets in a manner compatible with normal workstation and tomography setups. Two assumptions are made in order to make the process feasible and avoid reconstruction artifacts. The first one is that every partial dataset has images taken at the same rotation angle. The second is that for different datasets the rotation axis and image magnification are the same. Other approaches have already been in use on other beamlines for microtomography using TomoPy (Gürsoy et al., 2014 ) and transmission X-ray microscopy with TXM Wizard (Liu et al., 2012 ).

Fig. 1(a) shows the definition of the mathematical axes used in this manuscript. Radiographies have axes x₁,x₂ and reconstructions will have axes y₁,y₂. Fig. 1(b) shows that there may be two types of misalignment in our approach. One is between the partial datasets rotations axis and the other is between the final dataset axis and the camera axis. In order to ensure a good reconstruction it is necessary to align all partial datasets axes with the rotation axis. The main problem with this misalignment is defining the reconstruction paths. Normal reconstruction uses the saved data axis and does not calculate non-linear paths.

Figure 1
(a) Imaging plane I represented by axis $[{\bf{x}}]$ = (x₁, x₂) and slice plane S represented by axis $[{\bf{y}}]$ = ( y₁,y₂). (b) Representation of two images taken in a Cartesian axis (on the imaging plane I) that is different from the camera axis.

For the ideal case where all the experimental assumptions are satisfied, i.e. there is no rotation between measurements, the problem is simplified to finding the translation between two images. A brute-force algorithm can be used to test every possible transformation and minimize the residue of the overlapping area. This approach may be possible if the registration presents few degrees of freedom (such as simple translation). For registrations with more degrees of freedom (such as rotation and magnification), the computational complexity of the problem can be unsolvable. This is due to the fact that the brute-force algorithm may be computationally unfeasible. The difference between images may not be just geometrical but also include different noise and background levels, which makes the brute-force method even less accurate. For the generic case, finding the correlation between two images can be done using image marks (Pulli et al., 2012 ; Szeliski, 2006 ) and other techniques, such as the Fourier phase correlation method (Foroosh et al., 2002 ). For radiographies where there is not enough information to register images, another approach is the use of fiducial markers (Lemieux & Jagoe, 1994 ).

This article presents a novel way to correlate images in the radiography space and stitch in either radiography or sinogram space. This method expands the FOV of X-ray tomography experiments without needing to save the whole data into a new dataset. Projection correlation can be performed using the brute-force and cross-correlation algorithms presented in §2. Three challenging samples were measured in order to test our approach, as described in §3, and the results can be found in §4. In §5 the obtained results and further reconstruction ideas are discussed, and §6 contains a summary.

2. Standard correlation method

Let $[f \colon U \subseteq {\bf{R}}^2 \to {\bf{R}}]$ and $[h \colon V \subseteq {\bf{R}}^2 \to {\bf{R}}]$ be two-dimensional functions that represent distinct images as pictured in Fig. 2. Here, the domains U and V are such that $[U \cup V]$ will represent the domain of the resulting stitched area. Referring to Fig. 2, U and V are typically defined as

$[U = [-\infty, a] \times {\bf{R}}, \quad V = [a, \infty] \times {\bf{R}},]$

where x = a indicates the point where we assume that f and h correlate. We use the convention $[U\cup V \subseteq [-1,1]\times[0,1]]$ . The main idea of the process is to translate function h in such a way that at the boundary x₁ = a we obtain an optimal correlation. This means that, for every reasonably small $[\delta\,\gt\,0]$ , we look for a displacement $[{\bf{x}}]$ = $[{\bf{c}} \in {\bf{R}}^2]$ such that

$[f({\bf{x}}) = {\cal{T}}_{\!\!{\bf{c}}\,} [h]({\bf{x}}), \quad {\hbox{for all}}\,\, \left|x_1 - a\right| \,\lt\, \delta, \quad x_2 \in {\bf{R}}, \eqno(1)]$

where $[{\cal{T}}_{\!\!\bf c}]$ is the translation operator

$[{\cal{T}}_{\!\!{\bf{c}}\,}[h]({\bf{x}}) = h({\bf{x}} + {\bf{c}}). \eqno(2)]$

Here, we look for a displacement $[{\bf{c}}]$ that minimizes the quadratic residual

$[q({\bf{c}}) = \int\limits_{[a-\delta,a]\times{\bf{R}}} \big[\,f({\bf{x}}) - {\cal{T}}_{\!\!{\bf{c}}\,}[h]({\bf{x}})\big]^2\,{\rm{d}} {\bf{x}}. \eqno(3)]$

Operator $[{\cal{T}}_{\!\!{\bf{c}}}]$ is easily implemented using the Fourier transform through shifting property. In the computational framework, where functions {h,f} are represented by image matrices $[{\bf{H}}, {\bf{F}} \in {\bf{R}}^{n \times n}]$ we are looking for a vector $[{\bf{c}} \in {\bf{R}}^2]$ such that the extended matrix $[{\bf{P}}]$ ,

$[{\bf{P}} = \big[\,{\bf{F}},{\bf{T}}_{\bf{c}\,}{\bf{H}}\big] \in {\bf{R}}^{n \times 2n}, \eqno(4)]$

is not discontinuous at the boundary of images $[{\bf{H}}]$ and $[{\bf{F}}]$ , respectively. Matrix $[{\bf{P}}]$ represents the stitched image. Here, $[{\bf{T}}_{{\bf{c}}}]$ indicates the translation in pixels units, in any of the four possible directions. This operator is easily computed using the fast Fourier transform. With the above notation, we search for integers ( j,k) such that

$[\|{\bf{F}}(\colon, N \colon n) - {\bf{T}}_{(\,j,k)}\,{\bf{H}} (\colon,1 \colon N) \|_{\rm F}^2 \eqno(5)]$

is minimized, where $[\|\ldots\|_{\rm F}]$ is the Frobenius norm (Golub & Van Loan, 2012 ). Constant N in the above equation is the discrete equivalent of the parameter δ defined in (1), and can be given as a user input. In fact, before running the stitching process, it is visually easy to approximately define a number of columns to search the optimal shift.

Figure 2
Representation of the brute-force registration approach on the imaging plane. For a given translation vector $[{\bf{c}} \in {\bf{R}}^2]$ , function $[h({\bf{x}}+{\bf{c}})]$ correlates with f within the square $[|x_1 - a|\,\lt\,\delta]$ ; see text for details.

2.1. Phase correlation method

The method used for the mosaic reconstruction of whole datasets must be robust and fast for large images with different noise and illumination. Then, it is logical to use a less sensitive technique to find the registration between images. Also, using fiducial markers or image characteristics is not straightforward for most of the samples. A correlation method that suits the problem is to find the actual Fourier phase correlation and extract the shift from it. Although the phase correlation method (PCM) provides satisfactory results for some images, it sometimes generates the wrong shift (Preibisch et al., 2009 ).

The Fourier transform $[{\cal{F}}\colon f({\bf{x}}) \,\mapsto {\cal{F}}[\,f\,]({\bf{w}})]$ can be used in Fourier cross-correlation image processing (f.c.c.). Indeed, the f.c.c. output, r = $[r({\bf{x}})]$ , of two images f = $[f({\bf{x}})]$ and h = $[h({\bf{x}})]$ (see Fig. 2) is given by

$[r({\bf{x}}) = {\cal{F}}^{-1}\left [{{{{\cal{F}}[h]({\bf{w}})\,{\cal{F}}[\,f^*]({\bf{w}})} \over {\left|{\cal{F}}[h]({\bf{w}})\,{\cal{F}}[\,f^*]({\bf{w}})\right|}}} \right] ({\bf{x}}), \eqno(6)]$

with ^* standing for the complex conjugate and $[|\ldots|]$ for the absolute value. The maximum of this correlation image gives the absolute value of the linear translation that maximizes the correlation of the two images, i.e.

$[{\bf{c}} = \mathop {\rm argmax}\limits_{{\bf{x}}\,\in\,{\bf{R}}^2} r({\bf{x}}). \eqno(7)]$

PCM gives the argument of the shift vector between two datasets. As described by Preibisch et al. (2009), each PCM maximum gives four possible shifts for two-dimensional images and a subsequent pixelwise comparison of the overlap sectors finds the correct shift. Due to X-ray image noise and low contrast there are many local maxima in the cross-correlation image and the global maxima may not represent the true shift vector.

2.2. Hybrid correlation method

To make the correlation problem between complex images more feasible, we propose to use the rough approximation of the experimental motor shifts to create a correlation mask. The pixelwise multiplication of the correlation image and the mask makes it easy to find the correct shift. This mask can be used both for the brute-force and PCM algorithms. The a priori knowledge of the experiment geometry withdraws the need to test every possible shift, since the shifted image's relative motion is known.

In the current approach the final translation vector $[{\bf{c}}]$ is chosen by calculating the translation on every pair and removing the outliers and finding the mean. This approach is specially important when the sample goes out of the partial tomogram (on the lateral edges of the mosaic). Calculating the translation for every pair in a tomogram may be too demanding, so the proposed code also allows the number of samples to be reduced (i.e. calculate the pair for some given projections).

The calculated maps can be used for cone-beam geometry if the camera is shifted instead of the sample. That way the sample position in the cone is unchanged and the geometrical corrections can be made later at the mosaic reconstruction of each frame, and CT reconstruction algorithms of the Feldkamp–Davis–Kreuss (Feldkamp et al., 1984 ) family can be used. The proposed method increases the accuracy of the registration and decreases the computational load of finding the correct PCM shift.

3. Materials and methods

The experiment was carried out at the IMX beamline at the Brazilian Synchrotron Light Source and three challenging samples were measured in order to test the approach:

(i) The first sample was a Rosary seed (Abrus precatorius) of size ∼5 mm × 5 mm. It was first measured and stitched only with seven partial radiographs to show the artifacts that appear if there is no flat- and dark-field correction. The tomography approach was carried out with two partial datasets and the displacement vector $[{\bf{c}}]$ was found with both the brute-force approach and PCM.

(ii) The algorithm was tested with vertical filling in the partial datasets and the sample was a wood-fibre cylinder of radius 1 mm. It also consisted of two partial datasets that were acquired in a full rotation (of 1000 angles) and later broken into two partial datasets of 500 angles.

(iii) The approach was tested with a very large lateral mosaic array of a fire beetle (Pyrophorus noctilucus), using six partial datasets with 1000 angles.

The main reason for the sample choices was the challenge to apply the procedure using different types of materials and applications.

Every projection was acquired using a pco.2000 camera (https://www.pco.de/sensitive-cameras/pco2000/) coupled to a scintillator without binning (2048 × 2048 pixels with 16-bit unsigned integer values). Data processing was carried out with a 32-bit floating-point precision in order to keep the numerical error small. The calculated pixel size for the final images was 1 µm and the sample-to-detector distance was kept constant during the measurement of each dataset. This distance was optimized in order to obtain good phase-contrast conditions (Nesterets et al., 2005 ).

The proposed experimental approach relied on four separate steps: (a) acquiring datasets in a mosaic array; (b) finding the registration between datasets; (c) merging the datasets into a new single mosaic dataset; (d) reconstructing the new dataset.

Data acquisition was carried out in the same way as for a normal tomography experiment but the data were taken with only a part of the full dataset. As long as all the datasets combined ensure that the sample is contained within the new extended FOV, the tomographic reconstruction is possible. Images for correction of dark current and illumination structure (dark and flat images) were taken for every partial dataset in order to correct the frames before the final mosaic reconstruction (Als-Nielsen & McMorrow, 2011).

In the proposed scheme, image registration is performed by finding the relative phase between subsequent radiographs. Even though the sample translation stage may have a high resolution, the exact translation between frames cannot be taken with sub-pixel precision. In our experiment the vector $[{\bf{c}}]$ is constant along every mosaic reconstruction and it is possible to merge the final sinograms instead of the projections. In the case where the shift c is not constant along the partial datasets, stitching the sinograms would be a very challenging approach and projections should be merged instead.

The same approach was used in a full rotation tomography with the rotation axis shifted from the centre. The datasets for the mosaic reconstruction were obtained by separating the first and second half rotations, flipping the second one and using these as separate datasets for the proposed approach. To ensure compatibility with the already functional reconstruction algorithms and programs, the partial datasets were also rewritten into a complete dataset. This makes the partial acquisition and mosaic reconstruction invisible to the final dataset. There is no need to correct the pixel size with respect to experiment geometry in an experiment where the path of the X-rays through the sample is perpendicular to the camera axis. This means that overlapping pixels present the same information about the sample and no complex minimization is necessary to ensure the continuity of information in the final mosaic reconstruction. We found that the pixelwise mean of the overlapping area gives a good result in the final reconstruction (Als-Nielsen & McMorrow, 2011).

The reconstructions were carried out after processing the new sinograms with a centre correction and ring reduction (Miqueles et al., 2014a ) algorithms. The reconstruction was performed using a normal filtered backprojection (FBP) approach. According to the Nyquist sampling criterion, the optimal number of angles, $[N_\theta]$ , for a given number of elements in the reconstruction direction (rays), N_rays, is given by (Kak et al., 2002 )

$[N_\theta = {{\pi} \over {2}} \, N_{\rm rays}. \eqno(8)]$

To reconstruct our test samples with the appropriate number of projections according to equation (8) would make the experiment time unfeasible. Also, data storage and processing would be challenging for normal computers. One solution for this is to perform the mosaic tomography with fewer projections and reconstruct the data with an iterative algorithm (Miqueles et al., 2014b ; Miqueles & Helou, 2014 ; Sidky et al., 2010 ; Wen & Chan, 2012 ; Beck & Teboulle, 2009 ).

4. Results

Three samples were measured in order to validate and test the algorithm under challenging experimental conditions. Table 1 presents results obtained from these three experiments, which will be described next. It is important to note that, due to the nature of the experiments, it is natural to lose some slices at the top and bottom of the frame. In the third experiment shown below, we present intentionally lost slices, in order to state clearly the displacement vector $[{\bf{c}}]$ discussed in §2.

Table 1
Sample experiment description

Sample	Number of datasets	Number of angles	c	Final frame	Incomplete slices
Seed	2	500	(21231)	2050 × 3279	12
Wood	2	500	(8863)	2056 × 2911	12
Beetle	6	1000	(121845)	2111 × 11279	140

4.1. Experiment I: rosary seed

We begin by presenting the restored mosaic image without flat and dark corrections. Such an approach carries periodical artifacts, as shown in Fig. 3. The highlighted regions (marked by rectangles) appearing in Fig. 3 illustrate the presence of periodical artifacts and illumination differences due to the beam intensity distribution in the camera.

Figure 3
Mosaic reconstruction of a rosary seed without flat and dark correction. The insets shows the presence of periodical artifacts and different illumination between the images.

Using proper flat and dark corrections for each partial dataset, it is easy to notice the effect of the background correction, as shown in Fig. 4. Now, each partial frame has a background close to zero and pixel values corresponding to positive absorption information of the sample.

Figure 4
Partial radiographies of the rosary seed (with flat and dark correction).

Using the partial images of Fig. 4 and the proposed hybrid registration method, we found that the images had a translation $[{\bf{c}}]$ = (21231) pixels. Fig. 5 shows the map found using the PCM method and the brute-force residue map, described in §2. The expanded images correspond to the area after application of the mask. The results were the same for both methods but the computational time was significantly smaller for PCM since it relies on the fast Fourier transform and fewer operations. Fig. 6 shows the mosaic reconstruction of the partial projections in Fig. 4.

Figure 5
(Left) PCM correlation image and (right) correlation residue image (using the brute-force algorithm). The insets show the area of the mask and the local minima found.

Figure 6
Mosaic reconstruction of the rosary seed radiographies (depicted in Fig. 4

(with flat and dark correction). This image has dimensions 2050 × 3279 (slices × rays).

The reconstruction of the dataset was carried out with a normal FBP algorithm and one slice can be seen in Fig. 7. It is important to observe that no stitching artifacts or periodical artifacts are observed in this reconstruction. The final size of each restored projection was 2050 × 3279 (slices × rays) with 1 µm × 1 µm pixel size. The reconstructed mosaic dataset generated sinograms with 3279 × 500 (rays × angles) pixels. This number is below the Nyquist limit (Chesler et al., 1977 ) and analytical methods are not suitable for inversion of the slices.

Figure 7
Reconstruction of a single slice using the stitched 3D volume, as in Fig. 3

, for a sinogram with 3279 rays and 500 angles. The standard FBP was applied to the sinogram to obtain this reconstructed image.

The last result is sinogram registration instead of radiographs. In theory, if there are no displacements along the axis of the device, sinogram registration is equivalent to image registration of the frames. With noisy data and with several tiny displacements in the device, this is no longer true. Fig. 8 shows partial sinograms of the reconstructed slice of Fig. 7. Since the shift vector $[{\bf{c}}]$ has a component along the x₂ axis (the slice axis, see Fig. 1) it would be impossible to find the correct sinogram pair (a) and (b) without calculating $[{\bf{c}}]$ in the frames first. The partial sinograms in Fig. 8 differ in the slice number by $[{\bf{c}}_{y}]$ . For comparison, Fig. 8(c) depicts the resulting stitched sinogram of parts (a) and (b).

Figure 8
Truncated sinograms (a) and (b) giving rise to a complete sinogram (c), after a stitching process using the PCM method with pointing vector $[{\bf{c}}]$ in equation (7)

. The dashed red line shows the same ray on each sinogram.

4.2. Experiment II: wood fiber

Samples that generate quasi-periodical structures in radiographs are challenging to correlate. Fig. 9 shows the restored frame for a toothpick sample, acquired with a full rotation and later reconstructed using the mosaic approach. The final restored dataset was 2911 × 2056 × 500 (rays × slice × angle).

Figure 9
Mosaic reconstruction of a toothpick using only two datasets. This frame has dimensions 2911 × 2056 (rays × slices). The displacement vector $[{\bf{c}}]$ generates an empty space that cannot be seen in the image.

As described in Table 1, some slices were lost after mosaic reconstruction. Although imperceptible in Fig. 9, there is a blank space at the bottom that shows the displacement vector $[{\bf{c}}]$ . Removing the broken slices we finally perform the reconstruction, which is presented in Fig. 10 with resolution 2911 × 2911. No stitching marks or reconstruction artifacts can be seen in the reconstructed image.

Figure 10
Slice reconstruction of the toothpick using the restored dataset of Fig. 9

. Reconstructed image with dimensions 2911 × 2911.

4.2.1. Experiment III: beetle

A beetle (Cetonia Aurata) of size ∼1 cm was exposed to the imaging device. Measurements were made with six partial datasets to ensure that the algorithm was able to work with larger mosaic grids. Each dataset was gathered with 1000 angles and resulted in a final volume with dimensions 11279 × 2111 × 1000 pixels. Every single sinogram has 11279 × 1000 resolution (rays × angles), resulting in a final reconstruction slice of 11279 × 11279 pixels. Since this experiment was performed with a long lateral translation, the final component $[{\bf c}_{y}]$ was 60 = 5 × 12 pixels. This shift leads to several incomplete parts and a challenge to correlate the partial images. Figs. 11 and 12 show the mosaic reconstruction of one frame, using our approach and using FIJI MosaicJ (Schindelin et al., 2012 ), respectively. The black areas presented at the top and bottom of those images show the areas that do not present measured information and have to be discarded. It is important to notice that, even though the FIJI approach finds the best solution with angles between the partial datasets, it is not suitable for CT reconstruction. The angle between the datasets makes the correct reconstruction axis impossible to find, as described in the theory section.

Figure 11
Mosaic reconstruction of a beetle insect radiograph using six datasets. This image is intentionally displayed out of scale, since the horizontal axis is five times larger than the vertical one.

Figure 12
Mosaic reconstruction of a beetle insect radiograph using six datasets with FIJI MosaicJ (Schindelin et al., 2012

). The figure clearly shows a curvature, which does not favour a reconstruction scheme. The image is intentionally displayed out of scale, similar to Fig. 11

Fig. 13 shows a slice of the final reconstructed image obtained using the proposed approach. The size of each reconstruction is 11279 × 11279 pixels. An image with such dimensions presents a computational challenge in terms of reading and writing to normal hard disks. Currently, there are new methods able to reconstruct large data and the bottleneck for large image reconstructions is the backprojection, as described by Miqueles & Helou (2014). The reconstructed image in Fig. 13 presents strong streak artifacts due to a small number of angles. In fact, the FBP algorithm is not the best reconstruction scheme for this large dataset. Iterative techniques certainly provide better results, but iteration with such a large sinogram is still a challenge.

Figure 13
Tomographic reconstruction of the beetle insect with dimensions 11279 × 11279 pixels. Streak artifacts are clearly visible in this reconstruction due to the small number of angles.

5. Discussion

The proposed method presents a reliable way to extend the FOV of CT without needing to change the experimental setup. For a small increase in the FOV it was found that the method can find the shift of the sample axis in full rotation tomography acquisition. Then, this shift can be incorporated into other reconstruction routines or used to blend the dataset from two separate partial datasets.

Samples bigger than the FOV were imaged without losing pixel resolution or generating artifacts in the final reconstruction. It is not the purpose of this paper to discuss the amount of data generated, but it will be a challenge to handle large images for the computer process: storage $[\to]$ reconstruction $[\to]$ visualization. As data expand far from the rotation axis, even a very small angle within the camera can make the volume slicing obtain information from several slices. This problem would require realigning the whole data block in the memory before reconstruction and such an algorithm is not available yet.

In this article, all the slices were reconstructed using a standard FBP algorithm. As pointed out earlier, this is not the best reconstruction strategy for our tomographic setup. Indeed, since there are many missing angles, a constrained total-variation reconstructed image f^* would certainly provide better results, i.e.

$[f^{\,*} = {\hbox{argmin}} \big\{ {\hbox{TV}}(\,f)\semi \quad f \in S \big\}, \eqno(9)]$

with TV being the total variation operator (Velikina et al., 2007 ) and S the set of all two-dimensional mappings satisfying a consistency condition. A set S is determined by the Fourier slice theorem (Deans, 2007), i.e.

$[S = \big\{\, f \in {\cal{U}} \colon {\cal{F}}[\,f\,](\sigma_k \cos\theta_i, \sigma_k \sin \theta_i) = {\cal{F}}[\,g\,](\sigma_k, \theta_i)\big\}.]$

In the above equation, g is a sinogram image, typically obtained from a conventional imaging device, while $[\{(\sigma_k,\theta_i)\}]$ is a polar grid in the frequency domain. In the approach used in this manuscript, the sinogram g comes from stitching of the 3D volume in such a way that a gridding algorithm (Schomberg & Timmer, 1995 ) or conventional FBP can be used to recover the slice. For future applications of our stitching strategy, an ideal image reconstruction algorithm has to deal independently with each dataset and also with the pointing vector $[{\bf{c}}]$ defined in equation (7). With such an approach, it will not be necessary to store new datasets, and each part of the slice can be reconstructed independently, using a strategy like equation (9). Indeed, the memory needed to handle the reconstruction process grows linearly with the number of datasets used in the stitching part.

Fig. 14 illustrates a Fourier representation of two projections of a sample (at the same angle), e.g. $[S_1(\theta)]$ and $[S_2(\theta)]$ , giving rise to an incomplete frequency polar domain. Each acquisition comes from truncated sinograms, such as the ones shown in Fig. 8. In this manuscript, prior stitching is done so that $[S_1 \cup S_2]$ is a new dataset and the frequency domain is numerically dense so that a reconstruction scheme can be applied. Further reconstruction strategies as described by Miao et al. (2010 ) can also be applied to this problem.

Figure 14
Incomplete frequency domains, for a given angle θ, from two distinct unstitched datasets.

In this sense other acquisition–reconstruction methods can be designed to ensure the Fourier space density and solve the missing wedge problem (Arslan et al., 2006 ) that would arise if the sample could not achieve a full rotation in the proposed scheme.

6. Conclusion

We have shown that it is possible to find a more reliable correlation between partial tomographic datasets using existent methods and experimental information and generate a new sinogram without needing to save a new dataset. Our methodology was compared with commercial software (Schindelin et al., 2012), providing similar results. Using our approach it would be possible to extend the FOV without having to test the correctness of the mosaic reconstruction. It also does not need any change in the experimental setup for synchrotron tomography and is independent of the reconstruction algorithm. Hence, any reconstruction scheme that is already in use can take advantage of the hybrid mosaic approach to image samples larger than the camera FOV. Although increasing the size of CT data leads to several difficulties in reconstruction, it was shown that the normal FBP reconstruction can give coherent results as a first-order approximation of the solution.

Acknowledgements

This project was funded by a Brazilian National Council of Research and Development PhD scholarship. We would like to thank the staff of the IMX beamline, Frank O'Dowd and Gabriel Moreno, for help using the beamline, and also Mrs Edna Scola Klein for the rosary seeds.

References

Als-Nielsen, J. & McMorrow, D. (2011). Elements of Modern X-ray Physics. New York: John Wiley and Sons. Google Scholar
Arslan, I., Tong, J. R. & Midgley, P. A. (2006). Ultramicroscopy, 106, 994–1000. Web of Science CrossRef PubMed CAS Google Scholar
Beck, A. & Teboulle, M. (2009). IEEE Trans. Image Processing, 18, 2419–2434. Web of Science CrossRef Google Scholar
Brown, L. G. (1992). ACM Comput. Surv. 24, 325–376. CrossRef Google Scholar
Chesler, D. A., Riederer, S. J. & Pelc, N. J. (1977). J. Comput. Assist. Tomogr. 1, 64–74. CrossRef CAS PubMed Google Scholar
Deans, S. R. (2007). The Radon Transform and Some of its Applications. New Yok: Dover. Google Scholar
Feldkamp, L., Davis, L. & Kress, J. (1984). J. Opt. Soc. Am. A, 1, 612–619. CrossRef Web of Science Google Scholar
Foroosh, H., Zerubia, J. B. & Berthod, M. (2002). IEEE Trans Image Processing, 11, 188–200. Web of Science CrossRef Google Scholar
Golub, G. H. & Van Loan, C. F. (2012). Matrix Computations, Vol. 3. Baltimore: Johns Hopkins University Press. Google Scholar
Gürsoy, D., De Carlo, F., Xiao, X. & Jacobsen, C. (2014). J. Synchrotron Rad. 21, 1188–1193. Web of Science CrossRef IUCr Journals Google Scholar
Hansen, P. C. (1992). SIAM Rev. 34, 561–580. CrossRef Web of Science Google Scholar
Kak, A. C., Slaney, M. & Wang, G. (2002). Med. Phys. 29, 107. CrossRef Google Scholar
Kyrieleis, A., Ibison, M., Titarenko, V. & Withers, P. (2009). Nucl. Instrum. Methods Phys. Res. A, 607, 677–684. Web of Science CrossRef CAS Google Scholar
Lemieux, L. & Jagoe, R. (1994). Phys. Med. Biol. 39, 1915–1928. CrossRef PubMed CAS Web of Science Google Scholar
Liu, Y., Meirer, F., Williams, P. A., Wang, J., Andrews, J. C. & Pianetta, P. (2012). J. Synchrotron Rad. 19, 281–287. Web of Science CrossRef CAS IUCr Journals Google Scholar
Ma, B., Zimmermann, T., Rohde, M., Winkelbach, S., He, F., Lindenmaier, W. & Dittmar, K. E. (2007). Micron, 38, 492–499. Web of Science CrossRef PubMed Google Scholar
Mao, Y., Fahimian, B. P., Osher, S. J. & Miao, J. (2010). IEEE Trans. Image Processing, 19, 1259–1268. Web of Science CrossRef Google Scholar
Miller, C. C. (2006). Int. J. Geogr. Inf. Geovisualization, 41, 187–199. CrossRef Google Scholar
Miqueles, E. & Helou, E. (2014). European Conference on Mathematics for Industry (ECMI 2014), 9–13 June 2014, Taormina, Italy. Google Scholar
Miqueles, E., Helou, E. & De Pierro, A. R. (2014b). J. Phys. Conf. Ser. 490, 012148. CrossRef Google Scholar
Miqueles, E. X., Rinkel, J., O'Dowd, F. & Bermúdez, J. S. V. (2014a). J. Synchrotron Rad. 21, 1333–1346. Web of Science CrossRef CAS IUCr Journals Google Scholar
Nesterets, Y. I., Wilkins, S., Gureyev, T., Pogany, A. & Stevenson, A. (2005). Rev. Sci. Instrum. 76, 093706. Web of Science CrossRef Google Scholar
Preibisch, S., Saalfeld, S. & Tomancak, P. (2009). BioInformatics, 25, 1463–1465. Web of Science CrossRef PubMed CAS Google Scholar
Pulli, K., Baksheev, A., Kornyakov, K. & Eruhimov, V. (2012). Commun. ACM, 55, 61–69. Web of Science CrossRef Google Scholar
Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., Preibisch, S., Rueden, C., Saalfeld, S., Schmid, B., Tinevez, J. Y., White, D. J., Hartenstein, V., Eliceiri, K., Tomancak, P. & Cardona, A. (2012). Nat. Methods, 9, 676–682. Web of Science CrossRef CAS PubMed Google Scholar
Schomberg, H. & Timmer, J. (1995). IEEE Trans. Med. Imaging, 14, 596–607. CrossRef PubMed CAS Web of Science Google Scholar
Sidky, E. Y., Anastasio, M. A. & Pan, X. (2010). Opt. Express, 18, 10404–10422. Web of Science CrossRef PubMed Google Scholar
Szeliski, R. (2006). Found. Trends Comput. Graph. Vis. 2, 1–104. CrossRef Google Scholar
Velikina, J., Leng, S. & Chen, G.-H. (2007). Proc. SPIE, 6510, 651020. CrossRef Google Scholar
Wen, Y.-W. & Chan, R. H. (2012). IEEE Trans. Image Processing, 21, 1770–1781. Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 24| Part 3| May 2017| Pages 686-694

https://doi.org/10.1107/S1600577517001953