{'fraction_non_alphanumeric': 0.027699575320895165, 'fraction_numerical': 0.007563105406308155, 'mean_word_length': 4.930249009620827, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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 two approaches for geometry estimation. Our first approach is a calibration object correlation method in which 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 second approach is projection trajectory simulation, where we use a set of known intersection points and a sequential Monte Carlo method for estimating the posterior density of the parameters. We show numerically that a large set of parameters can be used for artefact-free reconstruction. We deploy Bayesian inversion with Cauchy priors for synthetic and real sawmill data for detection of knots with a very low number of measurements and uncertain measurement geometry.


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.
In noisy industrial settings, however, X-ray scanners do not always allow for direct measurement of the scanning geometry parameters describing the location and orientation of the source-detector pair [3]. There is therefore a need for estimation of the parameters in order to obtain the reconstructions. In this work, we show numerically that the parameter set providing good reconstruction is not unique, and a large set of parameters can thus be used for reconstructions of adequate quality, in the sense that we see all the features of logs needed for sawmill process decision making. In our specific industrial setting, we have sparse measurement geometry, for which we deploy Bayesian inversion with edge-preserving Cauchy priors, which enables detection of sharp features and material interfaces with a very low number of measurements.
CT image reconstruction can be described mathematically as an inverse problem. The aim when solving such problems is to find quantity of interest (causes) given its indirect measurement (observed effects) based on some forward model that establishes the connection between the former and latter quantities [4][5][6]. A two-dimensional X-ray CT reconstruction is an example of an inverse problem where the cross-sectional image of a scanned object is estimated from the measured intensity loss of multiple X-ray beams that penetrate the object from different angles (projections). Projections taken at different angles are collected to sinogram. Mathematically, the intensity loss or attenuation of X-ray beams can be represented by line integrals over the attenuation coefficient function describing the internal structure of the scanned object. The transform that maps a function on R 2 into the set of its line integrals is known as the Radon transform, which corresponds to the forward CT model. The objective of the inverse problem is to reconstruct the image of the object, that is, the attenuation coefficient function, using a combination of the forward model and noisy sinogram data [7,8].
A CT setting where the number of projections is extremely limited (sparse or limited-angle tomography) is known to be a severely ill-posed problem. Sparse data lead to severe artefacts when the image is reconstructed by conventional methods, for example, filtered backprojection [8]. The features on the boundaries between different regions are of particular interest in reconstructions. Therefore, various edge-preserving regularisation techniques have been proposed to avoid smoothing of edges in CT [9][10][11].
CT scanner geometry describes the location and orientation of the source-detector pair. Geometric inaccuracies are known to cause errors and severe artefacts in reconstructions [12,13]. In practical industrial settings in sawmills, the structure of an X-ray scanning device might not allow for the direct measurement of the physical dimensions of the source and detector and important geometry parameters such as the distance between the source and detector or the detector tilt. Moreover, geometrical displacements may be added to the system due to factors such as vibration during motion of either the scanner or scanned object [14]. There is, therefore, a need to address the issue of unknown geometry parameters. In this work, we demonstrate how the parameters can be selected and describe a number of misspecified geometry parameters, including first scanning angle, distance between the centre of rotation and the detector, source shift, detector shift, and detector tilt.
Many methods can be found in the literature for estimating geometry parameters, that is, geometry calibration, in CT systems. However, most research is concentrated on the medical field, where a single source-detector pair rotating around the scanning object is considered [15][16][17]. In the CT scanner used in this study, the measurement setting has the opposite configuration: the source and detector are fixed, and the object is rotated while scanning to get 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) self-calibration methods (i.e., methods without a calibration object) [18].
Calibration phantom-based methods, which is 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 [19].
A pioneering study [15] on imaging system geometry estimation utilised single photon emission computed tomography. In the study, the Levenberg-Marquardt optimisation algorithm was 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 linear parameter estimation approach is proposed in [17]. 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 [18], 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 a 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 [20,21]. In [22], data-driven estimation of unknown fan-beam 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 [23]. 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 Xray 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; and 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. To further investigate the identifiability and uniqueness of the geometry parameters, we run simulations on tracking of the projection trajectories on the detector to estimate the geometry parameter distribution by the adaptive Sequential Monte Carlo sampler.
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 recently developed class of priors, firstorder isotropic Cauchy difference priors [24]. 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 tomography reconstruction 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 and presents suggestions for future research.
2 Tomography reconstruction 2.1 Models X-ray tomography is a non-invasive imaging technique that allows reconstruction of the inner structure of an object to be obtained based on the observed attenuation of X-rays penetrating the object (the Beer-Lambert law). A single X-ray beam can be considered as a straight line passing through the object from the source to the detector, and the scanning process, therefore, results in a set of line integrals of an object, represented by a real-valued function x ∈ R 2 , along all the lines l in some manifold The transform (1) mapping function x to the set of its line integrals is known as the Radon transform, which was introduced by Johann Radon in 1917. In Radon's work (see [25] for its English translation), it is shown that the inner structure of an object can be reconstructed uniquely if the projection data from all possible lines passing through the object is available.
Since only a finite number of X-ray projections is available in many applications and the sinogram measurements are in practice done in pixel detectors, a discretised model for the X-ray imaging is used to find an approximate solution computationally. The model can be written as follows where measurements y ∈ R n and forward theory matrix A ∈ R n×m are known exactly, and we know also the statistical properties of noise e, which we assume to be a zero-mean Gaussian multivariate e ∼ N (0, Σ). The objective is to estimate the unknown attenuation coefficients x ∈ R m . This can be done either in Bayesian estimation or a classical deterministic framework [5].
In the Bayesian formulation of inverse problems, the observed data y, the unknown parameter x and noise e are treated as random variables. According to Bayes' theorem, the posterior distribution of unknown x can be formulated as where 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 [5,6]. The posterior distribution π post (x) is used to obtain point-estimates including a maximum a posteriori probability (MAP) estimate The MAP estimator can be computed in closed form for the linear Gaussian models and Gaussian priors, or, for example, via quasi-Newton optimisation methods in more complicated models.

Reconstruction methods
In order to reconstruct the internal structure of the object from its X-ray measurements, we use three reconstruction methods: Tikhonov regularisation, filtered backprojection and Bayesian inversion with Cauchy priors.

Tikhonov regularisation
Since the image reconstruction is an ill-posed inverse problem, different regularisation techniques are used to obtain a stable solution. One of the most common regularisation techniques, Tikhonov regularisation, is defined as a solution of the following minimisation problem where 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.

Filtered backprojection
The traditional method to reconstruct attenuation function x from the measurement data y is filtered backprojection (FBP), which is defined by the following equation where φ is a filter function, B is the backprojection operator, and * denotes the convolution operation [8]. Analogously, the discretised form of the filtered backprojection formula is defined as where φ is a filter, 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 [26] 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 this paper, FBP reconstructions with the Ram-Lak filter [26] were obtained using the Operator Discretization Library (ODL) [27]. ODL is a Python library developed at KTH Royal Institute of Technology, Stockholm. ODL focuses on inverse problems including CT reconstruction.

Bayesian inversion with Cauchy prior
As the third image reconstruction method, we employ the maximum a posteriori (MAP) estimate of the posterior distribution with the first-order isotropic Cauchy difference prior [24,28,29] 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 or the first-order differences when the function is discretised in a uniformly spaced lattice will be distributed.
Employing Cauchy distribution for the distribution of the differences effectively allows the existence of discontinuities in x by favouring features that are piecewise close to constant. These characteristics of the Cauchy difference prior are notably different compared to Gaussian difference priors and most other non-deep Gaussian random field priors due to the infinite variance of the Cauchy distribution [24]. We do not claim that the Cauchy difference prior is superior to Gaussian or other random field priors in log tomography in terms of reconstruction capabilities, but instead we use the Cauchy difference prior as an optional reconstruction tool due to its simplicity and unique features.
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 where β 2 is a parameter that controls the strength of the prior in terms of how close to zero the firstorder differences should be.
While the prior itself is improper, since we do not specify the distribution of the boundary nodes of x to make the integral of the prior probability density function to integrate to a positive constant, the likelihood function of the tomography setting virtually makes the posterior distribution proper. We find the MAP estimate by using the low-memory variant of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm.

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

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. Let us denote by θ the vector of the unknown geometry parameters. We seek to find vector of parameters θ based on the calibration phantom image x ref of N × N pixels and its mirrored image x ref . Let FBP(θ; y) be the corresponding FBP reconstruction (N × N pixels) 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 FBP(θ; y) at the shift index N 2 , N 2 , and the correlation between x ref and FBP(θ; y) at the same shift where denotes correlation. In other words, the objective function evaluates two possible calibration object candidates and returns the negated maximum of their correlations with the FBP reconstruction with given geometry parameters.
We also experimented with an objective function based on sinogram data directly rather than FBP, namelỹ where A θ is the forward CT operator matrix that depends on the geometry parameters θ. However, the objective function (10) demonstrated worse performance in comparison with the objective function (9) since the reconstructions with corresponding parametrisations contained severe artefacts. This can be explained by the regularity of the reconstruction space, which is more suitable for global optimisation than the sinogram space. From a numerical viewpoint, simultaneous evaluation of two calibration object scenarios in the objective function (9) 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 (9) for the unified optimisation workflow.
We seek parameters θ * that minimise the objective function J(θ) We follow the convention that the output of cross-correlation has the same dimensions as the first argument. Since the dimensions of images x ref and FBP(θ; y) are the same, their crosscorrelation at shift N 2 , N 2 is equal to the Frobenius inner product of the images when they are interpreted as matrices. If the filtered backprojection is denoted as f := FBP (θ; y), then where tr means the trace of a matrix, and the superscript T means transpose.

Differential Evolution
To find a solution to the optimisation problem (11), we use the differential evolution (DE) method [30]. 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 [31,32]. In this work, we adhere to the DE/best/1/bin strategy and use its implementation in the SciPy Python library [33] (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 for j = 1, . . . , D − 1, where r i,j ∼ U (0, 1), and P cross ∈ [0, 1] is a crossover (recombination) constant [30]. 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. Generation: Take the best individual 6: x best in the current population;

22:
Selection: 23: 24: x best ← u i ; 25: end if 26: end for 27: end while 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 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 parameterised FBP reconstructions. Consequently, local optimisation algorithms that do not require the objective function to be differentiable, such as the Nelder-Mead method [34], 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 of the measurement equipment is not interesting as such, since we are only interested in the quality of reconstructions.
In addition, we claim that the nature of the objective function (9) favors DE over simulated annealing (SA), which is another common global optimisation algorithm [35,36]. 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.
In addition, the population methodology of DE helps to explore the search space. The DE algorithm also seems to be rather insensitive to the choice of hyperparameters such as mutation rate and crossover probability, which also favours its selection for the given optimisation problem.

Measurement 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: Fig. 1 A scheme of the X-ray measurement geometry • 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, and the other five parameters are treated as unknown: • α 0 -the first scanning angle, • r D -distance between the COR and detector, • h S -the source shift (in the tangent direction), • h D -the detector shift (in the tangent direction), and • α D -the detector tilt, that is, the angle between the source-detector axis and the axis defining the detector orientation.

Synthetic data
Here, we illustrate the proposed method for geometry calibration with an application to a digital log phantom consisting 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, see Figure 2. We also simulate full-angle tomography measurements of a log slice phantom and images of two calibration phantoms using ODL.
To demonstrate the effect of misspecified geometry parameters on the quality of full-angle FBP reconstructions of log phantom, we tested various combinations of geometry parameters (see Figure 2): (a) the source and detector radii, r S and r D , are specified correctly; no source and detector shifts, h S and h D , are specified, that is, the X-ray source, COR and centre of the detector are assumed to lie on the same line; and no detector tilt, α D , is specified, that is, the detector is assumed to be perpendicular to the source-detector axis; (b) the source and detector radii, r S and r D , source shift, h S , are specified correctly; no detector shift, h D , and no detector tilt, α D ; (c) the source and detector radii, r S and r D , and source and detector shifts, h S and h D , are specified correctly; no detector tilt, α D ; (d) the source and detector radii, r S and r D , and detector tilt, α D , are specified correctly; no source and detector shifts, h S and h D ; (e) the detector radius, r D , is slightly perturbed, and all other parameters are specified correctly; and (f) all the parameters are specified correctly. The described test cases showed that misspecified geometry parameters lead to severe artefacts in the reconstructions.
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 seek to find optimal values for five unknown geometry parameters, that is, θ = [α 0 , r D , h S , h D , α D ] T , based on simulated X-ray data of two calibration phantoms with true parameter values that we set to [α 0 , r D , h S , h D , α D ] T = [2.55, 715, 320, 44, 0.28] T . Figures 3 and 4 show the results of the optimal parameter search for the 5-parameter case using the calibration phantom with a hole and the L-shaped calibration phantom, respectively. Although there is no unique parameter set and uncertainty in the parameter estimates is pretty high, all the different parametrisations provide satisfactory quality in the FBP reconstructions.
In additional experiments, we tested the role of the first scanning angle, α 0 , in the geometry estimation. We ran two sets of experiments: in the first set, we considered a fixed but incorrect value for the first angle, and in the second set, we considered the correct value for the parameter. In both cases, we ran optimisation for the other four parameters, that is, The results of the first set of experiments for the calibration phantom with a hole and the L-shaped calibration phantom are shown in Figures 5 and 6, respectively. In the case of the L-shaped phantom, the FBP reconstructions obtained with different parametrisations contain severe artefacts such as duplicated knots and halos around the log phantom, whereas the reconstructions are satisfactory with the calibration phantom with a hole. In both cases, however, the obtained optimal parameter sets are far from the true parameter values. Figures 7 and 8 show the results of the second set of experiments when the first scanning angle is specified correctly. The quality of the FBP reconstructions is good for both calibration phantoms, and the optimal parameter values found by the algorithm are quite close to the true values of the geometry parameters with much less uncertainty compared to the 5-parameter test case.
Based on the results of this study, it can be concluded that there is no unique parameter set providing good reconstructions when optimising five parameters, and uncertainty in the parameters is very high. If we fix the first scanning angle at its true value and optimise the remaining parameters, we reduce the uncertainty significantly and obtain values that approach the true parameter values while keeping the reconstruction quality high.
Although the scanning geometry is identified correctly, 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, that is, in the case of the sparse angle tomography normally found in industrial applications. To tackle the issue of artefacts, we employed Bayesian inversion with edge-preserving Cauchy priors as described in Section 2.2.3. As shown in Figure 9, in a very sparse case (10 projection angles), the FBP reconstruction contains severe star 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.

Real data
In this section, the proposed method for geometry calibration is applied to real X-ray data of a log with attached calibration disks. The X-ray measurements were obtained with an X-ray log scanning system provided by Finnos Oy.

Image registration
As a result of the X-ray log scanning, we obtained a number of projection images acquired from 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 mechanical instability of the measurement system and vibrations, there is an issue of missing pixel rows at either the top or bottom of the projection image, and there therefore is a need to overlay images corresponding to projections of the same log taken from different angles, which is a technique known as image registration [37].
In the context of our task, we adhere to the classical image registration method consisting of four consecutive steps [37]: feature detection, feature matching, transform model estimation, image resampling and transformation. The algorithm was implemented using the OpenCV Python library. Figure 10 shows artefacts arising from poorly defined geometry parameters. The simulation is done analogously to that in Section 4.2 but applied to the real log X-ray data. We estimated five unknown geometry parameters using the reference images of the calibration phantoms and the corresponding X-ray data obtained from the measurement device. The results of the geometry parameter search using the calibration phantom with a hole and the L-shaped calibration phantom are shown in Figures 11  and 12, respectively. The FBP reconstructions obtained with different parametrisations are good enough for detection of log features such as knots.

Optimal geometry parameter estimation
Ground truth (log phantom) Fig. 2 Ground truth image of a log phantom and its FBP reconstructions with different geometry parameters: (a) only the source and detector radii are specified (no source and detector shifts, no detector tilt); (b) the source and detector radii and source shift are specified (no detector shift, no detector tilt); (c) the source and detector radii, and source and detector shifts are specified (no detector tilt); (d) the source and detector radii and detector tilt are specified (no detector and source shifts); (e) the detector radius is slightly perturbed; (f) all the parameters are specified correctly After identifying geometry parameters based on the calibration phantoms, a set of optimal geometry parameters was used to obtain reconstructions of the scanned log using different reconstruction methods (FBP, Tikhonov regularisation, Bayesian inversion with Cauchy prior) as well as different numbers of projection angles (see Figure 13). It can be noted that the reconstruction quality remains good even when only 20 projection angles are considered. Moreover, the artefacts on the reconstructions arise from the reconstruction methods themselves (see, for example, the explicit star artefacts on FBP reconstructions with 45 and 20 projection angles) rather than from the geometry misspecifications.

Identifiability of the geometry parameters
In addition to the calibration object experiments with real and synthetic sinogram data, we also studied the identifiability and uniqueness of the geometry parameters in X-ray tomography through additional simulations. Instead of employing calibration objects, we calculated the intersection points of a detector plate and straight lines passing the X-ray source p 0 and randomly selected six points [p 1 , p 2 , . . . , p 6 ] T = P, with ∀i : p i ∈ [0, 1] 2 . First, we seek to estimate the distribution of five parameters θ = [α 0 , r D , h S , h D , α D ] T using      Reconstructions with different parametrisations Fig. 8 Optimal parameter search (4-parameter scenario, the first angle is fixed at its true value) for the synthetic data case using the L-shaped calibration phantom    the simulated measurement data set Y of projections at the detector plate when the distance r S between the X-ray source and COR is fixed.
In the second case, we treat r S as unknown as well, and append it to the end of θ. In total, Y consists of 20 × 6 scalar values, as we use 20 equispaced rotation angles between α 0 and α 0 + 2π. As the ground truth, we set the geometry parameters to [α 0 , r D , h S , h D , α D , r S ] T = [π/3, 8.0, −0.5, 2.8, 0.9, 4.5] T to generate the measurement data Y . We do not add any noise to Y , and we assume that the detector plate can be infinitely long. The simulation principle is shown in Figure 14, which illustrates the computation of Y at one rotation angle of the virtual fan-beam system. Variable y t,i denotes the one-dimensional coordinate at the detector plate of the projection of point p i when the system has been rotated by radians from α 0 . We model the posterior distribution π(θ|Y ) of the parameters through a Gaussian likelihood with a standard deviation of σ s = 10 −5 and a uniform prior where h(θ; P) is the function that calculates the intersections of the given points P at the detector plate using the geometry parameters θ, and π(θ) = U (l θ , u θ ). Function h returns a vector, so Y is a flattened vector as well.
Unlike in the experiments with the real and synthetic X-ray CT data, we use arbitrary geometry units in the trajectory simulations, and set l θ = [0, 0, −3, −3, −2] T and u θ = [2π, 10, 3, 3, 2] T . Likewise, when we have r S as one of the parameters, we set a uniform prior U (0, 10) for it.
We approximate the posterior distribution of the parameters with an adaptive Sequential Monte Carlo (SMC) sampler which leverages an adaptive Markov Chain Monte Carlo (MCMC) kernel to rejuvenate its particles [38]. The SMC algorithm samples the particles from their prior distribution and weights them by a tempered likelihood before resampling the particles if their effective sample size (ESS) is below a certain threshold. The algorithm then proceeds to the MCMC rejuvenation. The algorithm iterates through tempered distributions whose likelihood is powered to 0 < T j ≤ 1. The temperature T j must be greater than the previous temperature T j−1 , until the final temperature of 1.0 is achieved. The temperature schedule is selected so that the ESS decreases by a factor of 0.9 between SMC iterations [39,40].
We use 15000 particles within the sampler and use 200 MCMC steps at each temperature per every particle. SMC excels at sampling lowdimensional distributions that are multimodal or whose components are correlated in an ill-posed manner that poses problems for MCMC samplers. The components of the distributions we sample turned out to be very correlated, so preferring SMC to a plain MCMC turned out to be crucial. For more information regarding SMC, please see [41]. A pseudocode for the implementation used is given in Algorithm 2.
The distributions of the parameters of both experiments are plotted in Figure 15. To illustrate the identifiability issue of the geometry parameters, two projection trajectories of one of the six points were plotted. The red trajectory corresponds to the trajectory that is included in the simulated data set Y , and it was obtained with the ground truth geometry parameters. The blue trajectories are the projections of the same point, but the geometry parameters were drawn from the distribution given by the SMC sampler. The difference of the trajectory calculated with the sampled parameters to the ground truth trajectory is plotted in black.

5:
Tempering: Select T j > T j−1 so that 6: ESS of the particles after the next step 7: reduces by a factor 0.9, unless T j = 1 8: decreases the ESS less.

24:
end if 25: end while 26: procedure ESS(w) 27: The obtained parameter distributions agree with our previous experiments using real and synthetic X-ray data, even though in the previous experiments, we apply global optimisation to minimise a loss function rather than the Bayesian methodology. The parameters of the 5-parameter scenario are highly linearly correlated, as seen in Figure 15. SMC is able to find the ground truth values of the geometry parameters as the mean of the sampled parameter distributions. The low variance of the distributions and the small difference between the reference and a sample trajectory may indicate that a trajectory which is exactly the same as a reference one is obtainable only with the same five geometry parameters. In addition, the results would mean that the parameters are uniquely determinable if there is a way to precisely track the trajectories of known pointlike objects at the detector. We recall that our real and synthetic data experiments are based on projections of known and discretised objects at a detector plate which is also discretised into detector cells, so it is not possible in practice to track the trajectories of the objects with arbitrary precision.
Additionally, distributions having a fixed initial angle α 0 and r S as the fifth variable are yet to be studied. The effect of knowing the α 0 might be beneficial in terms of parameter identifiability. Such low uncertainty of the parameters is absent in the 6-parameter scenario in Figure 15. The inclusion of the additional geometry parameter (source radius) seems to introduce a way to compensate a change in one parameter with a nonlinear change to the other parameter while keeping the trajectories of the projections of the points at the detector almost the same. The geometry of the distribution is notably more difficult for the SMC, since the narrow marginal distributions of the parameters should be perfectly contiguous, which is not the case.
The 6-parameter scenario raises concerns whether the parameters are actually unidentifiable or should even smaller variance in the Gaussian likelihood be used along with a better SMC sampler implementation to get the parameter distributions to converge towards the Dirac distribution. As the prior distributions of the parameters could be set for even further range, we are careful with claims regarding the general identifiability of the parameters.
Conceptually, both the reference object correlation method and the projection trajectory simulations are similar, so for practical X-ray measurement equipment geometry parameter estimations, some of the parameters should be known if the others are to be determined uniquely up to certain uncertainty that arises from the measurement noise. If the objective is to estimate the parameters just to enable artefact-free CT imaging, non-uniqueness of the parameters can be disregarded.

Conclusions
In this paper, we have proposed an algorithm for 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. In addition, we studied parameter estimation via projection trajectory simulation and sequential Monte Carlo.
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, there is a large set of optimal parameters that delivers 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 the MAP estimate of the posterior distribution with first-order isotropic 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.
Future work should include segmentation of the cross-sectional log reconstructions obtained after the optimal geometry parameter search proposed in this paper. Segmentation allows for detection and location of knots, possible defects, and foreign objects inside the log, which is crucial for choosing the optimal sawing strategy for improved product quality and minimised waste.