Geometry parameter estimation for sparse X-ray log imaging

We consider geometry parameter estimation in industrial sawmill fan-beam X-ray tomography. In such industrial settings, scanners do not always allow identification of the location of the source-detector pair, which creates the issue of unknown geometry. This work considers an approach for geometry estimation based on the calibration object. We parametrise the geometry using a set of 5 parameters. To estimate the geometry parameters, we calculate the maximum cross-correlation between a known-sized calibration object image and its filtered backprojection reconstruction and use differential evolution as an optimiser. The approach allows estimating geometry parameters from full-angle measurements as well as from sparse measurements. We show numerically that different sets of parameters can be used for artefact-free reconstruction. We deploy Bayesian inversion with first-order isotropic Cauchy difference priors for reconstruction of synthetic and real sawmill data with a very low number of measurements.


Introduction
A particular case of sparse X-ray computed tomography (CT) problems is found in the sawmill industry, where X-ray scanners are employed for non-invasive detection of log features and defects such as knots, rotten parts, and foreign objects.Furthermore, X-ray scanners allow timber tracking from the raw material to the resulting product (wood fingerprinting) [1,2].Non-destructive feature detection and timber tracking provide information for selection of optimal sawing strategies and efficient process control.
CT image reconstruction can be described mathematically as an inverse problem [3][4][5].The objective of two-dimensional X-ray CT is to reconstruct the cross-sectional image of a scanned object from a collection of one-dimensional line projections.The projections are obtained by measuring the intensity loss of the beam of X-rays that penetrates the object from different view angles.The projections are collected to a sinogram.
In the sawmill industry, highly limited (sparse) tomographic data is used for log imaging.The sparse tomography is known to be a severely ill-posed problem.Sparse data lead to severe artefacts when the image is reconstructed by conventional methods like filtered backprojection [6].
To solve ill-posed inverse problems, a Bayesian approach is often applied [5].The Bayesian framework naturally allows for incorporating a prior distribution which encodes the possible internal structure of the object.First-order isotropic Cauchy difference priors are shown to promote sharp edges, which enables detection of sharp features and material interfaces with a very low number of measurements in the sawmill applications [7].
The CT scanner geometry describes the location and orientation of the source-detector pair for each projection of the CT scan with respect to the reference system of the scanned object.Geometric inaccuracies are known to cause errors and severe artefacts in reconstructions [8,9].In practical industrial settings in sawmills, the structure of an X-ray scanning device does not always allow for the direct measurement of the important geometry parameters such as the distance between the source and detector or the detector tilt.Moreover, geometric errors can arise due to factors such as vibration during motion of either the scanner or scanned object [10].There is, therefore, a need to address the issue of unknown geometry parameters.
In this work, we estimate the fan-beam geometry defined by 5 parameters: the first scanning angle, the distance between the centre of rotation and the detector, the source shift, the detector shift, and the detector tilt.In the particular application to sawmill log imaging discussed in the paper, a single source-detector pair is fixed while the scanned object is rotated to obtain its projections at different angles.

Related Work
Existing studies on geometry calibration can be roughly divided into two groups: 1) methods using calibration phantoms, and 2) selfcalibration (marker-free) methods.Calibration phantom-based methods, which are the more common approach, use one or more calibration objects (calibration phantoms) incorporating a known arrangement of radiopaque markers such as steel ball bearings.An extensive review of such approaches is given in [11].An overview of markerfree methods for estimating the crucial geometry parameter, the axis of rotation, is given in [12].Fig. 1 A scheme of the X-ray measurement geometry.The point source S and the flat-line detector D are placed around the log.The log is rotated around the centre of rotation COR to obtain different projections.The distance between the source and COR is r S and the distance between COR and the detector is r D .In the geometry, there are possible shifts of source h S and detector h D in tangent directions, and the detector tilt α D A pioneering study [13] on imaging system geometry estimation utilises single photon emission computed tomography.In the study, the Levenberg-Marquardt optimisation algorithm is used to calibrate a fan-beam geometry from projection measurements of reference objects consisting of point markers.The same principles are also appropriate for calibrating X-ray CT systems.
A calibration method for cone-beam geometry that employs a calibration phantom to estimate the geometry parameters is proposed in [14].The method uses a transparent cylindrical acrylic tube with 24 steel ball bearings as a calibration phantom, but the approach can be adapted also to other phantom configurations.
A self-calibration method is introduced in [15], where calibration is done in the cone-beam setting without calibration objects.The method performs calibration by registering two-dimensional projection data to a previously acquired threedimensional image of the scanned object.
Joint estimation of uncertain view angles determining the geometry of the forward CT model and the attenuation coefficient function of a scanned object is discussed in [16,17].
In [18], the geometry estimation, in particular, the center-of-rotation offset, is performed jointly with reconstruction based on a Bayesian approach that also provides information on the uncertainty of geometry parameters.
In [19], data-driven estimation of unknown fanbeam geometry is done by employing a neural network that learns the unknown forward operator from training data consisting of sinogram-image pairs.Another learned approach to tackle parametric uncertainties in the measurement operator is introduced in [20].The approach allows for joint automatic geometry parameter calibration along with reconstruction of the unknown CT image.

Our Contributions
We propose a new approach for identifying unknown geometry parameters in fan-beam X-ray log tomography by using a calibration object of known size that is scanned jointly with a log.
Compared to other existing works on geometry calibration, the difference in this approach is that: 1) it employs an easy-to-produce phantom for calibration of X-ray CT geometry with a fixed source-detector pair and rotating scanned object; 2) the approach itself is based on calculating the maximum cross-correlation between the ground truth image of the calibration object and its filtered backprojection (FBP) reconstruction; and 3) the approach allows estimating geometry parameters even for a limited number of projection angles.
Additionally, we demonstrate a number of possible artefacts that may appear in the reconstructions in the eventuality of misspecified geometry parameters when using a fan-beam configuration with a flat line detector.In extreme cases, misspecifications can lead to full unidentifiability, and in less severe cases, the reconstructions can contain duplicated features, halos, shifts, and blurring.We aim to demonstrate these features and present a way to overcome the issue with suitable parameter choices.
Finally, we consider the sparse-angle tomography and employ maximum a posteriori estimates with a recently developed class of priors, firstorder isotropic Cauchy difference priors [21].We demonstrate the amenability, for our specific use case, of this method over conventional methods (filtered backprojection, Tikhonov regularisation) when the number of projection data is extremely limited.

Outline of the Paper
The rest of this paper is organised as follows: In Section 2, we review the tomography reconstruction model and methods.The new geometry parameter estimation algorithm is presented in Section 3. Synthetic and real data CT examples are demonstrated in Section 4. Section 5 concludes the study.

Tomography reconstruction 2.1 Fan-beam geometry setting
The schematic plot of the X-ray scanner geometry used in the numerical experiments is shown in Figure 1.An X-ray source, placed around the log, emits fan-shaped X-ray beams that are recorded by the corresponding flat line detector.
The geometry of the fan-beam model is determined by seven parameters in total.There are two known parameters: • r S = 859.46(mm) -distance between the X-ray source and centre of rotation (COR), and • n D = 768 (px) -number of detector elements, the pixel size is 2mm, and the other five parameters are treated as unknown: • α 0 (rad) -the first scanning angle; the projection angles are defined w.r.t. the first scanning angle, and are assumed to be uniformly distributed between the first and last angles, • r D (mm) -distance between the COR and detector, • h S (mm)-the source shift (in the tangent direction), • h D (mm)-the detector shift (in the tangent direction), and • α D -defines the detector tilt.In the default configuration, the initial source-to-detector vector is [0, 1], and the initial detector axis is perpendicular to it (i.e.[1, α D = 0]).However, if the detector is tilted then component α D should be adjusted accordingly.
We further denote those five unknown parameters as a 5-dimensional vector

Discrete CT imaging model
Using the Beer-Lambert law [6], the intensity loss of X-ray beams passing through a scanned object can be represented as where x is the attenuation coefficient describing the internal structure of the object, I 0 and I φ,τ are the initial and final intensities, respectively, and l φ,τ denotes the line corresponding to the acquisition angle φ ∈ [0, 2π) and detector location τ ∈ R.
Taking logarithms in Equation ( 1) we obtain the following where the left-hand side y := − ln(I φ,τ /I 0 ) is measured data known as projection.The collection of projections is commonly referred to as a sinogram.
The right-hand side of ( 2) is known as the Radon transform x(s)ds, that maps the unknown attenuation coefficient function x into the set of its line integrals [22].In other words, the known measured data is represented by the Radon transform, which corresponds to the forward CT model.In practice, only a finite number of X-ray projections is available and the sinogram measurements are recorded by the pixel detectors.If the number of acquisition angles is K, the number of detector pixels is M , and a scanned object is discretised into a N × N grid, we obtain the discrete model for the X-ray CT imaging where x ∈ R N 2 ×1 is the unknown vector of atten- uation coefficients, A θ ∈ R KM ×N 2 is a so-called system matrix (forward matrix) describing intersection lengths for each of KM lines.The system matrix is assumed constant throughout the reconstruction and depends on the scanning geometry parameters θ that are not known exactly.Vector y ∈ R KM ×1 specifies the measurement data (projections) taken during the scanning process, and e denotes measurement noise.
The corresponding inverse problem aims to reconstruct x from the available measured sinogram data y.This can be done either in Bayesian estimation or a classical deterministic framework [4].In the paper, we use and compare three different reconstruction methods (filtered backprojection, Tikhonov regularisation and MAP estimate) discussed in Section 2.3.

Reconstruction methods
The traditional method to reconstruct the crosssectional image of an object from the measurement data is filtered backprojection (FBP).If y denotes the measurement data and θ is the geometry parameter vector, then the discretised form of the FBP formula is defined as where ϕ is a convolutional filter in the projection domain, B θ is discrete backprojection, and * denotes the discrete convolution.
Several frequency domain filters are commonly used in FBP, including the most basic Ram-Lak filter [23] and its windowed modifications such as Shepp-Logan, Cosine, Hamming and Hann filters.The filter function significantly affects the robustness of the reconstruction and should be carefully selected based on the particular reconstruction task.
In sparse CT where the number of projections is extremely limited, the FBP method leads to severe artefacts in reconstructions due to illposedness of the problem [6].Different regularisation techniques have been used to solve ill-posed problems.One of the most common regularisation techniques, Tikhonov regularisation, is defined as a solution of the following minimisation problem The regularisation parameter α > 0 is chosen depending on the noise strength to find a tradeoff between the solution stability and its ability to explain the measurements.
In the Bayesian formulation of inverse problem (3), the observed data y, the unknown parameter x and noise e are treated as random variables.
The noise e ∼ N (0, Σ) is assumed to be a zero-mean Gaussian multivariate with covariance matrix Σ.
According to Bayes' theorem, the posterior distribution of unknown x can be formulated as The prior distribution π pr (x) incorporates a priori information about the unknown x before the measurements were obtained, and π(y|x) is the data-likelihood that contains information about the forward operator A θ and assumptions about the statistical distribution of noise e [4,5].
The features on the boundaries between different regions are of particular interest in reconstructions.In sawmill applications, sharp-interface reconstructions are needed for measuring the size and shape of tree knots.Various edge-preserving regularisation techniques have been proposed to avoid smoothing of edges in CT [24][25][26].
As an edge-preserving method, we employ the maximum a posteriori (MAP) estimate of the posterior distribution with the first-order isotropic Cauchy difference prior [7,21,27] and the following data likelihood term Unlike the standard Tikhonov regularisation (5), which assumes no special structure for x except that it should be close to zero, the first-order Cauchy difference prior incorporates prior information on how the first-order derivatives of x will be distributed.Employing Cauchy distribution for the distribution of the differences effectively allows the existence of discontinuities If x is assumed to be discretised in a twodimensional uniform grid indexed by (i, j) with size N × N , then the isotropic Cauchy difference prior distribution can be expressed as follows [21] where β 2 is a parameter that controls the strength of the prior in terms of how close to zero the firstorder differences should be.

Geometry parameter estimation
For geometry calibration, we use two different calibration phantoms of known size and shape: an L-shaped calibration phantom and a calibration phantom with a hole (see Figure 2).

Geometry estimation model
The parameters of the partially unknown geometry of the measurement equipment are estimated using an objective function that evaluates the correlation of a reference image to an FBP reconstruction with given parameters.More precisely, we compare both the reference image and its mirrored image to the reconstruction, because the ambiguity of the direction of the detector plates and the rotation direction of the object inside the scanning device may render the reconstructions mirrored as well.The ambiguity of the handedness of the reconstructions can be taken into account in the objective function by calculating the maximum of correlations of both handednesses.We seek to find a vector of unknown geometry parameters θ based on the calibration phantom image X ref of N ×N pixels and its mirrored image X ref .To separate the reference image from its flattened vector form x ref , we denote the image (N × N pixels) with the capitalised X ref .
Let f θ = FBP(θ; y) be the corresponding FBP reconstruction that depends on the input X-ray (sinogram) data y and geometry parameters θ.Then the objective function J(θ) is defined as the negative maximum of the correlation between X ref and f θ and the correlation between where Tr means the trace of a matrix, and the superscript T means transpose.In other words, the objective function evaluates the Frobenius inner product of the FBP reconstruction with given geometry parameters and two possible calibration object candidates.It returns the negated maximum of the values.In this study, we do not employ the normalised correlation in the method, because the unnormalised correlation proved effective and robust.For the second, usage of the normalised correlation requires addressing the exceptional case when the standard deviation of the reconstruction is zero.This might be the case when the differential evolution proposes a set of geometry parameters that results in an empty reconstruction in the domain of the reference object, although such a case is very unlikely.We also experimented with an objective function based on sinogram data directly rather than FBP, namely where A θ is the forward CT operator matrix that depends on the geometry parameters θ.However, the objective function (9) demonstrated worse performance in comparison with the objective function (8).More specifically, FBP reconstructions using the set of found geometry parameters contained severe artifacts and were often misaligned with the reference object.We argue this is because the optimisation schemes cannot propose and find good sets of geometry parameters to align the reference and measured sinograms when the correlation is performed in the projection space, where the differences between the perturbed and unperturbed geometry parameters appear relatively subtle.Instead, the effect of perturbation of the geometry parameters is typically more substantial in the reconstruction space, which helps the optimiser to find the potential global optimum of the cross-correlation, which, in turn, ensures the proper alignment of the geometry.From a numerical viewpoint, simultaneous evaluation of two calibration object scenarios in the objective function (8) may pose additional challenges to the optimisation algorithms as the likelihood of getting trapped to a local minimum might increase.Evaluating the two calibration object cases separately is plausible.However, even an objective function that evaluates the correlation with one calibration object at times is non-differentiable, so we use the objective function presented in Equation ( 8) for the unified optimisation workflow.

Differential Evolution
The aim of geometry parameter estimation task is to find parameters θ * that minimise the objective function J(θ) To find a solution to the optimisation problem (10), we use the differential evolution (DE) method [28].The method belongs to evolutionary algorithms that apply natural mechanisms of evolution theory such as mutation and selection to evolve a solution to an optimisation problem.Let S ⊂ R D be the search space of the optimisation problem to be solved.Then the DE algorithm keeps a population of N pop individuals (D-dimensional vectors) each of these individuals stands for a possible solution of the problem.At every pass through the population, the algorithm mutates each candidate solution by other candidate solutions to create a mutant candidate.Different mutation strategies are used in DE [29,30].In this work, we adhere to the DE/best/1/bin strategy and use its implementation in the SciPy Python library [31] (the pseudocode is listed in Algorithm 1).In this strategy, the best member of the population, x best , is mutated (the Mutation stage) by the difference of two other randomly selected members x k1 and x k2 , (k 1 ̸ = k 2 ̸ = best) to form a mutant vector m i as follows where µ ∈ [0, 2] is the mutation constant that controls the amplification of the differential variation (x k1 − x k2 ).
Next, at the Crossover stage, the trial vector is formed as follows u i,j = m i,j , if r i,j < P cross , x i,j , otherwise, for j = 1, . . ., D − 1, where r i,j ∼ U (0, 1) is a random number drawn from the uniform distribution on the interval [0, 1], and P cross ∈ [0, 1] is a crossover (recombination) probability [28].In the numerical experiments, we use P cross = 0.7.The final entity u i,D is always calculated as Eventually, at the Selection stage, the fitness of the trial candidate is assessed: if the trial candidate is better than the original candidate then it takes its place.
The choice of the DE method for our problem of geometry parameter optimisation stems from the fact that DE is able to cope with nonlinear and non-differentiable objective functions and requires only a few control variables.It is also a global for i = 1..N pop do 5: Generation: Take the best individual 6: x best in the current population; 7: Choose 2 random integers 8: Mutation: Form a mutant vector 11: Crossover: 13: for j = 1..D − 1 do 14: Generate r i,j ∼ U (0, 1); Selection: x best ← u i ;

26:
end for 27: end while optimisation algorithm, which is important for our geometry parameter search.There are several local optima in the parameter space when evaluating the correlation between a reference object and the parametrised FBP reconstructions.Consequently, local optimisation algorithms that do not require the objective function to be differentiable, such as the Nelder-Mead method [32], are not suitable for this task.It must be recalled that the objective function likely does not even have a unique optimum: there might be multiple parameter vectors that give almost the same cross-correlation with the reference object when used in the FBP reconstruction function.However, the uniqueness of the geometry parameters In addition, we claim that the nature of the objective function (8) favors DE over simulated annealing (SA), which is another common global optimisation algorithm [33,34].Differential evolution utilises a set of particles that mimic the natural selection, evolution and mutation of the best candidates for the global maximum of the objective function.When some of the geometry parameter candidates have been mutated closer to a possible optimum, the evolution steps help to find the precise maximum in a more robust manner than SA.

Results and discussion
In the numerical experiments, we used the Operator Discretization Library (ODL) [35].ODL is a Python library developed at KTH Royal Institute of Technology, Stockholm.The library focuses on inverse problems including CT reconstruction.
We simulated the scanning geometry, computed forward projections, Tikhonov and FBP reconstructions using ODL.To compute MAP estimates we used the low-memory variant of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm implemented in Julia.
To assess the quality of reconstruction, we report relative reconstruction error defined as where image reconstruction x and reference image x true are in the vector form.As the other image quality metric, we employ the structural similarity index measure (SSIM) since it is shown to mimics well the perceived quality of an image by the human visual system [36,37].In numerical experiments, the full-angle FBP reconstruction is used as a reference.

Synthetic data analysis
Here, we illustrate the proposed method for geometry calibration with an application to a digital log phantom (Figure 3).The phantom consists of concentric rings modelling tree growth rings, two ellipses representing knots inside a tree trunk, and one small ellipse with higher density representing a foreign object, for example, a metallic bullet.1: Relative error ε rel and SSIM for full-angle FBP reconstructions of log phantom obtained with different parametrisations.The parametrisations were found using the L-shaped calibration phantom and the calibration phantom with a hole.True parameters are provided as a reference.Geometry parameter vectors for the L-shaped phantom correspond to those in Figure 5 Fig. 6 Pairwise plot of geometry parameters θ = [α 0 , r D , h S , h D , α D ] T and the objective function J(θ).Optimal geometry parameter values (in blue) were obtained by running the calibration method with full-angle (360 projections) measurements of the L-shaped calibration phantom 100 times.True values of θ and J(θ) are marked as red dots

Geometry parameters
We simulated full-angle tomography measurements of the log phantom and images of two calibration phantoms using ODL with true geometry parameter values [α 0 , r D , h S , h D , α D ] T = [2.55,715, 320, 44, 0.28] T .To avoid the inverse crime [4], we created the synthetic measurement data on a grid with N = 1013, then interpolated to a coarser grid with N = 256 and added relative error noise of 2%.
Figure 4 shows the effect of misspecified geometry parameters on the quality of full-angle FBP reconstructions of log phantom.We tested various combinations of geometry parameters (zero value means that the corresponding geometry parameter is misspecified, otherwise -specified correctly): In the figure, reconstructions corresponding to cases (b), (c) and (e) contain severe artefacts such as duplicated features in the log phantom To tackle the issue of unknown geometry parameters, we ran the parameter search according to the method described in Section 3. In the parameter search, we aimed to find optimal values for five unknown geometry parameters, that is, θ = [α 0 , r D , h S , h D , α D ] T , based on Figure 5 shows the full-angle FPB reconstructions with different geometry parameters found using the proposed calibration method.The geometry parameter search was done using the simulated full-angle measurements of the L-shaped calibration phantom, the reconstruction with the true geometry parameters is provided in the figure as a reference.It can be noted that all the different parametrisations provide satisfactory quality in the FBP reconstructions although there is no unique optimal parameter set.
We also ran the geometry parameter search using the other calibration phantom -phantom with a hole.The obtained results were similar to those of L-shaped calibration phantom in Figure 5, so we omit the plot.Instead, in Table 1, we report relative reconstruction error ε rel and SSIM for FBP reconstructions with different geometry parameters found by the geometry parameter search using two calibration phantoms.It can be concluded that both calibration phantoms perform well in the geometry parameter search algorithm and the obtained estimates for the geometry parameters result in the reconstructions of good quality.
To provide a view on the uncertainty of the found parameters, we ran the optimisation algorithm with L-shaped phantom 100 times and created the pairwise plot of optimal parameters together with the true parameter values (see Figure 6).It can be noted that the parameters are highly correlated and their distribution is very narrow.When using the calibration phantom with a hole we observed similar behaviour of the parameters so we omitted the corresponding plot.
In addition, we tested the performance of the calibration method when the sparse tomography measurements of the calibration phantoms were used.In Figure 7, we show the results of geometry parameter search using sparse measurements (20 projection angles) of L-shaped phantom.The results are similar to those in experiment in Figure 6, where the full-angle measurements of the calibration phantom were used.In Figure 8, two full-angle FBP reconstructions of log phantom Fig. 10 Two-dimensional reconstructions of the log slice obtained by different reconstruction methods with various numbers of projection angles (image resolution 256 × 256 pixels).The geometry parameters were obtained by the calibration method using the full-angle measurements of the L-shaped calibration phantom are shown for comparison: geometry parameters for the left reconstruction were obtained using the full-angle measurements of L-shaped phantom while for the right reconstruction, sparse measurements of the phantom were used.Red contour lines highlighting the phantom features are used to show that there is a small angular misalignment between two reconstructions.Despite the misalignment both reconstructions are of good quality, therefore, the proposed calibration method works in the sparse-angle scenario as well.
Although the scanning geometry is estimated correctly using the calibration method, severe artefacts still appear in the reconstructions if classical reconstruction approaches such as FBP and Tikhonov regularisation are used and the number of projection angles is very low.To tackle the issue of artefacts, we employed Bayesian inversion with edge-preserving Cauchy priors from Section 2.3.As shown in Figure 9, in the sparse case (20 projection angles), the FBP reconstruction contains severe streaking artefacts and the Tikhonov regularised reconstruction is extremely blurry, while the quality of the MAP estimate with a Cauchy prior is still good enough and allows for detection of knots and other features in the log phantom.
For MAP reconstructions in the figure, the relative error ε rel grows representing the expected quality deterioration as the number of projection angles decreases.However, SSIM rises suggesting the quality improvement, which contradicts the visual quality assessment and ε rel .This phenomena can be explained by the reduced ability of SSIM to represent the visual quality for CT reconstructions [38] and the images containing salt-and-pepper (impulse) noise [39].

Real data analysis
In this section, the proposed method for geometry calibration is applied to real X-ray data of a log.The X-ray measurements were obtained with an X-ray log scanning system provided by Finnos Oy.A tree log was scanned jointly with the wooden calibration phantoms (Figure 2) attached to the saw cut of the one side of the log.The calibration phantom plane was, therefore, perpendicular to the log axis of rotation.As a result of the sequential ("step and shoot") X-ray scanning, we obtained projection images capturing the log with attached calibration object from 360 different angles.The width of the projection image is fixed and corresponds to the number of detector elements n D .The height of the projection image depends on the number of slices that have been scanned (the width of 1 slice is 1 px, or equivalently, 2 mm), that is, the length of the log scanned.Due to the peculiarities of the measurement system, there is an issue of missing pixel rows at either the top or bottom of the projection image, and there is a need to overlay images corresponding to projections of the same log taken from different angles using an image registration technique [40].In the context of our task, we adhere to the following image registration method consisting of four consecutive steps: feature detection, feature matching, transform model estimation, image resampling and transformation.The algorithm was implemented using the OpenCV Python library.
We estimated five unknown geometry parameters (the first scanning angle, distance between the COR and detector, the source shift, the detector shift, and the detector tilt) using the reference images of the calibration phantoms and the corresponding X-ray data obtained from the measurement device.In Figure 10, we show the cross-sectional reconstructions of one slice of the log obtained after geometry parameter estimation using the full-angle measurements of the L-shaped calibration phantom.The figure compares different reconstruction methods (FBP, Tikhonov regularisation, MAP estimate) for various numbers of projection angles.In sparse measurements scenario (45 and 20 projection angles), the quality of FBP reconstructions is compromised by streaking artefacts, while Tikhonov and MAP reconstructions show relatively better quality.
To show the strength of the calibration method applied to sparse data, we provide two full-angle reconstructions of the same cross-sectional slice of the log obtained with different parametrisations (see Figure 11).Geometry parameters of the left reconstruction were obtained using the full-angle measurements of the the L-shaped calibration phantom, whereas parameters of the right reconstruction were estimated from sparse measurements (20 projection angles).Both reconstructions are of good quality and can be used to locate the tree knots and possible defects in the log.

Conclusions
In this paper, we have proposed an algorithm for the estimation of unknown geometry parameters in fan-beam X-ray computed tomography (CT) for log imaging.The algorithm employs easy-toproduce calibration phantoms that are scanned jointly with the log and applies differential evolution optimisation to the objective function represented by the maximum cross-correlation between the ground truth image of the calibration phantom and its filtered backprojection (FBP) reconstruction.The algorithm is able to estimate geometry parameters not only from full-angle measurements of the calibration phantom but also from sparse measurements.
We demonstrated numerically that there might be multiple parameter vectors that give almost the same cross-correlation with the reference object when used in the FBP reconstruction, that is, different sets of optimal parameters deliver an adequate quality of reconstructions.However, the uniqueness of the geometry parameters of the measurement equipment does not pose any interest, since we are only interested in the quality of the reconstructions.
Additionally, we demonstrated the effect of misspecified geometry parameters, including the first scanning angle, the centre of rotation (COR) to detector distance, the source and detector shifts, and the detector tilt, on the quality of CT reconstructions.For reconstruction purposes, we employed Bayesian inversion (the MAP estimate of the posterior distribution with Cauchy difference priors) as opposed to classical reconstruction approaches such as FBP and Tikhonov regularisation in order to improve the reconstruction quality in the sparse CT setting.Bayesian inversion showed better performance on sparse data than FBP method that suffered from severe streaking artefacts compromising the quality of reconstructions.

Fig. 2
Fig. 2 Images of the L-shaped calibration phantom (left) and the calibration phantom with a hole (right).All dimensions are given in mm

Fig. 3 A
Fig. 3 A digital log phantom representing the cross-sectional image of a tree log (image resolution 256 × 256 pixels)

Algorithm 1 1 :
Differential Evolution (DE) Input parameters: D -problem dimensionality; N pop -population size; µ -mutation constant; P cross -crossover probability; fobjective function under consideration; 2: Initialisation: Create population with N pop individuals and initialise them with random positions in the search space; 3: while the termination condition is not satisfied do 4:

Fig. 4
Fig. 4 Example artefacts in reconstructions when geometry parameters are incorrectly specified: (a) the source and detector shifts, detector tilt are misspecified; (b) the detector shift and tilt are misspecified; (c) the detector tilt is misspecified; (d) the detector and source shifts are misspecified; (e) the detector radius is slightly perturbed; (f) all the parameters are specified correctly.FBP reconstruction (f) with true parameters is used as a reference for computing ε rel and SSIM

Fig. 5
Fig. 5 Full-angle FBP-reconstructions of log phantom (image resolution 256 × 256 pixels): true geometry parameters (framed subplot) and different parameter vectors found by the calibration method with the L-shaped phantom.FBP reconstruction with true parameters is used as a reference for computing ε rel and SSIM

Fig. 7 Fig. 8
Fig. 7 Pairwise plot of geometry parameters θ = [α 0 , r D , h S , h D , α D ] T and the objective function J(θ).Optimal geometry parameter values (in blue) were obtained by running the calibration method with sparse (20 projections) measurements of the L-shaped calibration phantom 100 times.True values of θ and J(θ) are marked as red dots

Fig. 9
Fig.9Reconstructions of the log phantom obtained by different reconstruction methods with various numbers of projection angles (image resolution 256×256 pixels).The geometry parameters were obtained by the calibration method using the full-angle measurements of the L-shaped calibration phantom

Fig. 11
Fig. 11 Full-angle reconstructions of log slice.Left: geometry parameters are estimated using full-angle measurements (360 projection angles); right: geometry parameters are estimated using sparse measurements (20 projection angles)