key: cord-0788037-t1m4xg31
authors: Lukić, Tibor; Balázs, Péter
title: Limited-view binary tomography reconstruction assisted by shape centroid
date: 2021-01-12
journal: Vis Comput
DOI: 10.1007/s00371-020-02044-8
sha: c379041b151ee661e372bbaa32fc129249524d95
doc_id: 788037
cord_uid: t1m4xg31

In this paper, the binary tomographic reconstruction problem for very limited projection data availability is considered. Being this inverse problem highly ill-posed, we propose a new reconstruction model that uses a shape centroid-based regularization term, i.e., we assume that the center of gravity of the object of interest is known, at least approximately, in advance. Motivation for this regularization is found in the close connection between the projection data and the object centroid, as we will show. Experimental evaluation underpins that reasonable results can be obtained from practically minimal amount of projection data, gathered from just one projection direction.

Tomography [10] is a field of image processing which deals with reconstruction of images from a number of projections. It has a particularly important role in investigating the non-accessible or non-visible interior of objects, in a non-invasive way. Tomography methods allow, among others, non-destructive industrial testing [7, 16, 34] and also a big variety of diagnostic approaches in medicine [1, 6, 13, 17, 28, 29, 32, 35] . In the standard protocol for COVID-19 case ascertainment, computer tomography (CT) chest scan of the patient plays also an essential role, which makes the tomographic reconstruction additionally topical.

The aim of this paper is to study the shape centroid as a priori information in the context of binary tomography. The centroid point of an object can be calculated in terms of discrete geometric moments. This property makes this feature more suitable for use in discrete (digital) spaces (see [25] ). Based on the knowledge of the centroid, we introduce B Tibor Lukić tibor@uns.ac.rs Péter Balázs pbalazs@inf.u-szeged.hu 1 Faculty of Technical Sciences, University of Novi Sad, Novi Sad, Serbia 2 Department of Image Processing and Computer Graphics, University of Szeged, Szeged, Hungary a regularization term into the reconstruction model, with the purpose to make reconstructions possible from very limited projection data, practically from a minimum quantity, taken from a single projection direction.

The main contribution of our work include:

-We reveal the connection between horizontal and vertical projections and the centroid of a binary image. -We construct an energy function incorporating prior information on the centroid of the image to reconstruct from its projections. -We derive the gradient of this energy function in an analytical way to optimize the function by the spectral projected gradient algorithm. -We underpin by experiments that a single direction projection data can be augmented by the prior information on the centroid of the image in order to increase reconstruction quality.

The paper begins with a brief overview of related works and the results achieved in the field (Sect. 2). Then, we give the formalization of the reconstruction problem (Sect. 3) and recall the concept of geometrical moments (Sect. 4). In Sect. 5, we show that reconstructions of objects from two orthogonal projection directions completely preserve the centroid. Knowing this fact, we analyze the opposite possibility: can the centroid, as a priori information, "replace" projection data? Hoping for a positive answer, the object centroid as prior knowledge is incorporated into the reconstruction process, with a new regularization term (Sect. 6). Our aim is to develop a method capable of providing acceptable reconstructions based on a single projection direction. The "missing projection data" obtainable from another (orthogonal) direction is compensated by the information about the object centroid. Experimental results, presented in Sect. 7, confirm the validity of this idea. Finally, in Sect. 8 we conclude our results.

From mathematical point of view, the tomographic projection process can be defined by integrals or sums. For a continuous image function u(x, y), the calculation of projection values can be modeled by the Radon transform [31] 

where α is the angle of the projection ray relative to the given coordinate axis, while s denotes the distance from the ray to the origin of the coordinate system. The tomographic reconstruction problem belongs to the class of inverse problems and deals with the reconstruction of the original function u(x, y) based on the given projection data R(α, s), for some fixed α. The filtered backprojection and its variants [26, 27] are the most often used methods for image reconstruction. However, for a good image quality, they require a large amount of projection data. An important class of tomographic reconstruction problems using much fewer projections belongs to the field of Discrete Tomography (DT) [11, 12] . In DT problems, the range of the image function is a finite and discrete set, known in advance. In practical applications, DT deals with reconstructions of digital images that have only a few different gray levels. Binary tomography (BT) is a special case of DT and its scope is the reconstruction of binary images. Although BT is a restricted subfield of DT, it has various applications, as in many real situations, the object being studied consists of a single known material, and therefore, the absence or presence of this material can be represented by 0 and 1 in the images, respectively.

Among the possible applications, we mention the field of human X-ray angiography, where the aim is to reconstruct images representing blood vessels or heart chambers, using X-ray tomography methods. Injecting a contrast agent with high linear attenuation coefficient into the part of the body being examined, the problem can be solved with BT: one can seek for the presence or absence of the contrast agent in certain positions [8, 30] .

Another field of the applications of BT belongs to industrial non-destructive testing, where the goal is to examine the quality or degree of degradation of the interior of objects made of homogeneous materials, without damaging the objects themselves. Such objects include, for example, jet engine turbine blades or cylinder heads in internal combustion engines [4, 14] .

There are a number of excellent reconstruction methods suggested for discrete or binary tomography, such as the discrete algebraic reconstruction technique [3, 42] , convexconcave regularized algorithms [33] , simulated annealingbased approaches [40] and multi-well potential-based techniques [18] . However, they require data from at least two or even more projection directions. In BT, the projection data is often amended by prior information of the image to be reconstructed, such as, e.g., smoothness, convexity, texture [2, 33, 38, 40] .

To our knowledge, the only BT reconstruction method suggested for single projection availability is based on shape orientation [19] . Our intention is to go a step further and give a new contribution in this direction of research.

We follow the strategy of formalizing the binary tomographic reconstruction problem as a system of linear equations

where N and M denote the number of image pixels and projection rays, respectively. Furthermore, u is the vector representation of the image to reconstruct, and A is the socalled projection matrix whose row entries a i,· represent the length of the intersection of the pixels and the i-th projection ray passing through them (see Fig. 1b ). The projection vector b contains the acquired projection data. Coordinates of this vector b i are calculated as it is explained in Fig. 1b . The side length of each pixel is assumed to be one unit. Moreover, we assume parallel beam projection geometry, i.e., for each projection direction a number of equidistantly placed parallel projection rays are taken (being the distance equal to the side length of the pixels), as it is shown in Fig. 1a .

Under the problem of reconstruction, we understand the determination of the solution image u of (1). (The matrix A and the projection vector b are given.) It is important to notice that this system has a binary constraint, and it is also often underdetermined (N > M), i.e., there are more pixels (unknowns) than projection rays (equations). Transforming (1) into an optimization problem and using regularization to incorporate prior information about the solution is a commonly used technique to treat this issue (see, e.g., [19, 33] ). 

Given a 2-dimensional continuous image u (i.e., a function u : Ω → R, Ω ⊆ R 2 ), the geometrical moment (or shortly, the moment) of order

When u is a digital image (i.e., a function u : Γ → Z, Γ ⊆ Z 2 ), we can calculate the digitized moments as

It is worth noting that this discrete version approximates very well the original continuous definition (2) (for more details we refer to [15] ). In addition, moments are desirable operators in discrete spaces, and many regularization procedures in image inverse problems, such as tomographic reconstruction or denoising, are based on discrete moments (see, e.g., [19, 25] ). Especially, we focus now on the center of gravity (or centroid) of the image which can be expressed in the terms of first-order geometrical moments. The centroid of image u is defined by

Let u(i, j) be a digital image of size m × n (i = 1, . . . , m; j = 1, . . . , n). The vertical projection data corresponding to the k-th column, v k , is calculated in the following way:

In the same manner, the horizontal projection data is calculated as

Now, we can easily show that the vertical and horizontal projection data completely determines the centroid. Indeed, moments in (3) can be expressed as functions of h k and v k :

We are interested in the opposite issue, that is, how useful it is to utilize the centroid as a priori information, to achieve a successful regularization? In other words, can the centroid be a good substitute for projection data, especially for vertical or horizontal? We believe in a positive answer and therefore propose a new reconstruction method with centroid based regularization. Our focus is on the reconstructions obtained from projection data acquired in a single direction. Since in this case only a minimal amount of data is available, it is normal to expect reconstructions being far from perfect. Nevertheless, our goal is to reach acceptable reconstructions, that already allow us to get important information about the original objects.

We propose an energy-minimization reconstruction model, where the information about the centroid of the original object is incorporated in the reconstruction process. Formally, the proposed model is given by

where u is a vector representation of the image u(i, j), u(1, n) , . . . , u(m, n)] T , while the energy or objective function is given by

The energy function (5) is a sum of three nonnegative terms. The first one, also called the data fitting term, expresses consent of a solution with the given projection data. The role of the second term, also called smooth regularization term, is to ensure the smoothness or homogeneity of the solution. Here, r and b mark the index of the pixel adjacent to u i , to the right and to the bottom, respectively. The application of this regularization is based on the assumption about the compactness of the original image, in the sense that it consists of holeless regions. The third term measures the accordance of the centroid point of the current solution (C x (u), C y (u)) with the given centroid coordinates (C * x , C * y ), being this latter the centroid of the original object (or at least its proper approximation). In industrial applications, such information can be gained either by a CAD model or by studying a blueprint object. In medicine, in case of slice-by-slice imaging, the similarity of adjacent cross-sections can be exploited. Adjacency can be regarded not only in space but also in time when the same slice is reconstructed several times within a short period to investigate dynamic processes [30] . We note that in case of in-vivo experiments, binary tomography is only applicable when an imaging protocol ensures that only the organ of interest is X-ray projected. This can be achieved, e.g., by Substraction Computed Tomography (SCT) [9, 17, 37] .

The impact of different terms in the energy function (5) is controlled by parameters w P > 0, w H > 0, and w C > 0 (associated to projection data, homogenization, and centroid, respectively). In the special case, when w C = 0, the energy (5) reduces to the energy function that has been already successfully applied in several tomographic reconstruction problems [20, 33, 40] which motivates our choice to insert the second term into (5) .

Owing to the binary conditions, the minimization problem (4) belongs to the class of constrained minimization problems. The feasible set {0, 1} N is discrete and finite. Therefore, application of a gradient type optimization approach is not possible in this formulation. However, we can transform this problem into a relaxed form, with continuous feasible set:

where τ = [1, 1, . . . , 1] T . The concave regularization term u, τ − u is added to enforce the binary solution. Its minimum value (zero) is achieved when all pixel intensities of u are equal to zero or one, i.e., u is a binary image. Therefore, minimization of this term during the reconstruction process enforces the binary solution. This effect is used to ensure the final binary solution. Parameter μ controls the magnitude of the influence of the term u, τ − u in the reconstruction process. We note that this type of binary regularization for tomography problems was first introduced in [33] (where the reader can find more details), and later successfully applied in a number of similar problems, see [19, 40] . The gradient of the energy function in (6) can be determined in analytical manner:

Elements of the gradient of the smooth regularization term are calculated by

where u a and u l designate neighbors above and to the left of the pixel u i , respectively. For the gradient of the centroid based regularization term, we have

Applying elementary calculus, we can show that 1, 1, . . . , 1, 2, 2, . . . , 2, 3, 3, . . . , 3, . . . , m, m, . . . , m] T .

With this, we have fully determined the calculation of the gradient ∇ E R (u) of the proposed energy function in (6) . Therefore, the proposed reconstruction model can be solved be a gradient based minimization method. For this task, we use the spectral projected gradient iterative algorithm, originally proposed in [5] . Our choice is motivated by the fact that this algorithm has been successfully applied in similar problems [19, 23, 25] . The process is repeated, gradually increasing the value of μ, until the final (binary) solution is reached. The pseudocode of the proposed reconstruction method is shown in Algorithm 1. 

Parameters: out = 10 −3 ; μ Δ = 0.01; u init = [0.5, 0.5, . . . 0.5] T ; μ = 0. while max i {min{u init i , 1 − u init i }} > out do /*

In this part of the paper, we deal with the presentation and evaluation of experimental results. The goal is to investigate how the newly proposed centroid prior can enhance the reconstruction quality, especially in situations of very limited projection availability.

In the experiments, nine binary test images of size 64 × 64 are used as originals (see Fig. 2 ). The first six images are synthetic (phantom) ones, and they represent arbitrarily created different shapes. It is obvious that for each shape the centroid also takes significantly different positions. The last three test images (liver, bone, and cell) are acquired from real data. Liver is based on a binary segmented CT image [41] . Bone image is obtained from a histological CT image of a bone implant, inserted in a leg of a rabbit [36] . Cell is a binary segmented fluorescence image of a calcein-stained Chinese hamster ovary cell [24, 39] . The results obtained by the proposed method, in the following denoted by CENT, based on the centroid prior are compared to other methods relevant for the case of low projection availability. Algorithms used for comparison are: the spectral projected gradient (SPG) algorithm, suggested in [21] and also applied in [22, 23] ; and the approach based on the shape orientation prior information (ORI), proposed in [19] . The SPG and ORI are energy minimization type methods, and their functioning can be shortly illustrated by the following common energy function:

where Φ(u) denotes the orientation of the shape presented in image u, and α * is the given a priori known orientation value. The weighting parameter w O ≥ 0 controls the orientation term. The first and second terms are identical as in the energy (5) . According to suggestions in [19, 21] , in the presented experiments the value of w H is set to 0.5 in both methods. The SPG method does not use a priori information, and its energy is derived from (7) when w O = 0. The ORI method, however, relies on the orientation prior information and in that case w O is set as a positive number. In our applications, w O is set as 0.1, according to the suggestion in [19] . The best values are typeset in boldface Comparison of results for phantom images obtained using three different methods: SPG, ORI and CENT

The parameters in the proposed energy function (5) are experimentally determined to provide the best possible results and are set to w P = 0.1, w H = 0.5, and w C = 0.2. In Fig. 3 , the effect of the binary solution enforcing regularization of (6) is illustrated. During the reconstruction process, the parameter μ is gradually increased by μ Δ . Its proper adjustment is important: Too large values can lead to an overemphasized binarization effect and consequently to a solution of lower quality; on the other hand, too small values of μ Δ may cause unwanted increase in execution time. In our experiments, we set μ Δ = 0.01.

We analyze reconstructions of test images obtained for projection data acquired from different projection directions. Three error measures are considered to express the quality of the reconstructions. The first one indicates the number of misclassified pixels in reconstructions, and it is denoted by P E (pixel error). This error is calculated by

where u r and u * represent the reconstructed and the original image, respectively. We note that P E can also be understood as the Hamming distance of the original and the reconstructed images. By r P E we denote the relative number of misclassified pixels, and its value is calculated as

Finally, P R E refers to the projection error, which measures how much the reconstruction matches with the projection 

In addition, the deviation in centroid (DC) is calculated by

where u r and u * are the reconstructed and the original images, respectively. As we already showed in Sect. 5, vertical and horizontal projection data completely determines the centroid. In Fig. 4 , we investigate how the centroid is preserved when the angle between the two projection directions is gradually decreasing. We deduce that the centroid is well-preserved (the DC values are small) even when the angle between the two direc-tions is very small. Major changes happen only when the angular difference is about 1 • or smaller (see greater DC values). This short experimental study also confirms that the centroid prior is worth to be applied mostly in cases when only one projection is accessible. Therefore, in the following, we consider reconstructions obtained by a single projection direction.

In Figs. 5 and 6, we focus on reconstructions relying on a single (horizontal and vertical, respectively) projection, which, in general, does not ensure sufficient information for preserving the centroid of the original object. The absolute and relative pixel errors indicate that in all of these experiments introducing the suggested centroid prior significantly contributes to improving the quality of reconstructions. Without a priori information about the centroid, the object may be slightly translated in the fixed projection direction, which can cause a high pixel error, even when the reconstruction is "visually appealing." The alternative way to apply the information about the centroid is to use it after the reconstruction process, by simply translating the obtained result such that its centroid coincides with the centroid given beforehand. The reconstruction error in that case can be expressed by the centered pixel error, first introduced in [19] , given by

where ϕ is a linear image translation operator that maps the center of gravity of u r into the center of gravity of u * . The corresponding relative centered pixel error is then defined by

In Figs. 5 and 6, the C P E value for each considered test image is also indicated. In all cases, these values are greater than P E values for the proposed method with centroid prior (CENT). This confirms that it is more beneficial to incorporate the centroid prior directly in the process of reconstruction, as we suggested, rather than just utilizing it after the reconstruction is completed. Similar results can be obtained when, instead of the vertical and horizontal directions, we use arbitrarily chosen ones. In Fig. 7 , a representative example is presented showing reconstructions of the PH5 test image using projection data acquired from angles 160 • or 20 • . The results confirm that the proposed method with the use of the centroid prior provides significant enhancement in quality.

Tables 1 and 2 collect the reconstruction errors for all test images using projections from different directions (0 • , 45 • , 90 • , 135 • ). For each projection direction, results obtained by three different reconstruction methods are presented. Regarding the relative pixel errors, in 26 out of 36 test cases (i.e., in 72.22% of the cases) the application of centroid prior (CENT) provided better results than the other two approaches, and only in one experiment (PH6, 0 • ) the method without any of both features (SPG) provided the best result. The SPG is better than CENT in the following cases also: (Cell, 0 • , 90 • ) and (PH5, 45 • ). These test images (PH5, PH6 and Cell) represent thin shaped figures, when the centroid does not have to carry enough shape information. Taking a look at the projection error we, observe that-except in the case of PH6 and projection direction of 90 • -in all cases the projection error is also higher when there is no prior information. Finally, not surprisingly, in almost all the cases the DC distance is the smallest for the centroid-based approach.

The goal of this paper was to provide a new binary tomographic approach for very limited (or minimal) projection availability. We proposed the centroid of the original object as a new regularization feature. We then designed the reconstruction as an energy minimization model and solved the optimization by a deterministic gradient-based algorithm. In the experimental evaluation, we compared the new method to another recently proposed one which uses the orientation prior [19] as well as with a classical reconstruction method that uses no prior information [21, 40] . Experiments show clear advantages of the proposed method for reconstructions using data from a single projection.

Based on the above, our general conclusion is that the proposed reconstruction method which includes the shape centroid regularization is a useful tool in reconstruction problems with very limited projection accessibility, especially in situations when just a single projection is available.

Interactive skeletonization of intensity volumes

A measure of directional convexity inspired by binary tomography

DART: A fast heuristic algebraic reconstruction algorithm for discrete tomography

Discrete tomography methods for nondestructive testing

Algorithm: 813: SPGsoftware for convex-constrained optimization

The Biomedical Engineering Handbook, Third Edition -3 Volume Set

Industrial X-ray Computed Tomography

The reconstruction of three-dimensional obejcts from two orthogonal projections and its application to cardiac cineangiography

Imaging of pulmonary perfusion using subtraction CT angiography is feasible in clinical practice

Fundamentals of Computerized Tomography: Image Reconstruction from Projection

Discrete Tomography: Foundations, Algorithms and Applications

Advances in Discrete Tomography and Its Applications

An implicit skeleton-based method for the geometry reconstruction of vasculatures

Image reconstruction and correction methods in neutron and X-ray tomography

On discrete moments of unbounded order

Image-based material characterization of complex microarchitectured additively manufactured structures

A parallelized 4D reconstruction algorithm for vascular structures and motions based on energy optimization

Discrete tomography reconstruction based on the multiwell potential

Binary tomography reconstruction based on shape orientation

Regularized image denoising based on spectral gradient optimization

A spectral projected gradient optimization for binary tomography

Deterministic discrete tomography reconstruction method for images on triangular grid

Regularized binary tomography on the hexagonal grid

Deterministic defuzzification based on spectral projected gradient optimization

A non-gradient-based energy minimization approach to the image denoising problem

Radon's inversion formulas

Why do commercial ct scanners still employ traditional, filtered backprojection for image reconstruction

Segmentation of cam-type femurs from CT scans

3D segmentation and labeling of fractured bone from CT images

Binary reconstruction of the heart chambers from biplane angiographic image sequences

On the determination of functions from their integral values along certain manifolds

Maxillofacial Cone Beam Computed Tomography: Principles, Techniques and Clinical Applications

Discrete tomography by convex-concave regularization and D.C. programming

Registration of CAD mesh models with CT volumetric model of assembly of machine parts

Automatic detection of pulmonary nodules in CT images based on 3D res-i network

Defuzzification of spatial fuzzy sets by feature distance minimization

Pretransplantation evaluation of the cirrhotic liver with explantation correlation:accuracy of CT arterioportography and digital subtraction hepatic angiography in revealing hepatocellular carcinoma

Binary tomography using variants of local binary patterns as texture priors

Algorithms for cytoplasm segmentation of fluorescence labelled cells

A benchmark evaluation of large-scale optimization approaches to binary tomography

Automatic liver segmentation method in CT images

A multi-channel dart algorithm

Acknowledgements The research of Tibor Lukić was supported by projects OI-174008 and III-44006 of the Ministry of Education and Sciences of the R. of Serbia. He also acknowledges to the DOMUS project of the Hungarian Academy of Sciences. Péter Balázs was supported by Grant TUDFO/47138-1/2019-ITM of the Ministry for Innovation and Technology, Hungary. This research was supported by the project "Integrated program for training new generation of scientists in the fields of computer science," No. EFOP-3.6.3-VEKOP-16-2017-00002. The project has been supported by the European Union and co-funded by the European Social Fund.

Conflict of interest The authors declare that they have no conflict of interest.Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.