SIViP DOI 10.1007/s11760-014-0693-9

ORIGINAL PAPER

3D image reconstruction using Radon transform Mario D’Acunto · Antonio Benassi · Davide Moroni · Ovidio Salvetti

Received: 7 March 2014 / Revised: 29 August 2014 / Accepted: 31 August 2014 © Springer-Verlag London 2014

Abstract Computed tomography as well as magnetic resonance or positron electron tomography is currently the most commonly used medical imaging modalities for the analysis of human body complex structures and organs, where diseases must be recognized and identified. The image reconstruction process used in these tomography techniques is usually based on the Radon transform (RT). In this paper, an algorithm based on Markov chain Monte Carlo methods for reconstruction of 2D–3D structures is introduced, including correction of noise directly from the RT sinograms. The 3D reconstruction is independent by the instrumental technique and can be applied to any RT-based medical imaging technologies. Keywords Radon transform · Image reconstruction · Markov chain Monte Carlo

M. D’Acunto (B) · A. Benassi · D. Moroni · O. Salvetti Institute of Information Science and Technologies, National Research Council of Italy (ISTI-CNR), Via G. Moruzzi 1, 56124 Pisa, Italy e-mail: [email protected] A. Benassi e-mail: [email protected] D. Moroni e-mail: [email protected] O. Salvetti e-mail: [email protected] M. D’Acunto Institute of Structure of Matter, National Research Council of Italy (ISM-CNR), Via Fosso del Cavaliere, 100, 00133 Tor Vergata, Rome, Italy

1 Introduction An important problem in image processing is to reconstruct a cross section of an object from several images of its projections. A projection is a shadowgram obtained by illuminating an object by penetrating radiation. Each horizontal line shown in this figure is a one-dimensional projection of a horizontal slice of the object. Each pixel of the projected image represents the total absorption of the penetrating radiation along its path from the source to the detector. By rotating the source-detector assembly around the object, projections for several different angles can be obtained. The goal of image reconstruction is to obtain an image of a cross section of the object from a collection of projections. Imaging systems that generate such slice views are called computerized tomography (CT) scanners. The RT is the underlying fundamental concept used for CT scanning, as well for a wide range of other disciplines, including radar imaging, geophysical imaging, nondestructive testing and medical imaging [1–8]. In the context of biomedical imaging, the reconstruction of complete 3D models of dendritic structures can be an important tool in many imaging applications. In general, a dendritic model can accurately quantify both 3D topology (connectivity) and morphology (radius, length, branching angles, etc) of the dendritic structure. Recently, several algorithms were produced where the main emphasis is put on topology or morphology, or both [9]. Inferential methods have been largely used for the identification of vascular dendritic structures. A powerful way is to use a Bayesian method, such as maximum a posteriori estimation (MAP). Nevertheless, a principal difficulty with Bayesian techniques is a computational one. In fact, such methods require the use of Markov chain Monte Carlo (MCMC) algorithms, which may run for hundreds of thousands of iterations to yield reliable results. A partial solution is given by the usage of methods that combine local

123

SIViP

likelihood maximization with global structure determination using a general model of structures within a collection of specific samples. In this paper, we introduce an algorithm that uses RT for the reconstruction of 3D structures recorded in conditions of low signal-to-noise ratio (SNR). The reconstruction is applied to a nucleolus of an osteoblast cell recorded using a confocal microscope. Poissonian and Gaussian noise sources, as well as some instrumental artifacts, are treated before the 3D reconstruction [10]. The results should improve the resolution of the medical images, making possible to quantify spatial sources observed by indirect measurements as in the functional magnetic resonance (fMRI) applied to brain study, where the measured signal is due to a changes of hemoglobin, a paramagnetic molecule, due to increased brain activity as a response of an external stimulus.

where the unit vector η using spherical coordinates is now written as: η = (cos φ sin θ, sin φ sin θ, cos θ ). Limiting ourselves to bi-dimensional case, since P(−η, ρ + π ) = P(η, ρ), we will mostly use a single-valued parameterization assuming that −∞ < η < ∞, 0 ≤ ρ ≤ π . In addition, we can consider a class of function { f } to which to apply the projection theorem prescribing that: P(χ , ρ) = f (χ cos ρ, χ sin ρ) where



(3)

P(η, ρ)eiχ η dη

P(χ , ρ) =

(4a)

and f (χ cos ρ, χ sin ρ)   = f (ρ cos θ − l sin θ, ρ sin θ + l cos θ ) eiχ η d2 r

2 Radon transform: theoretical background

(4b)

The Radon transform (RT) of a two-dimensional function f (x, y) is defined as follows. Let A be a 2D model in 2 and let f (x) be the bi-dimensional binary function of A, defined as f (x) = 1, when x lies within the 2D model’s surface, and f (x) = 0 otherwise. Let, also, η = (cos θ, sin θ) be a unit vector in 2 and ρ a real number. The 2D Radon transformation R f (η, ρ) of f (x) is a function which associates with each pair (η, ρ) the integral of f (x) on the plane P(η, ρ) = {x | xt · η = ρ}, where the superscript t indicates vector transposition. This plane is normal to the direction η and a distance ρ to the origin [11] +∞ +∞ f (x, y) δ (x cos θ + y sin θ −ρ) dxdy

R f (η, ρ) ≡ −∞ −∞ +∞

=

f (ρ cos θ − l sin θ, ρ sin θ + l cos θ )dl −∞

(1) where the δ-function converts the two-dimensional integral to a line integral dl along the line x cos θ + y sin θ = ρ. The transformed function R f (η, ρ) is often referred to as the sinogram of f (x, y), because a point δ-function in f transforms to a sinusoidal line δ-function in P. The generalization of Eq. (1), in the case of f (x, y, z), is straightforward defined as: +∞ +∞ +∞ R f (η, ρ) =

f (x, y, z)δ −∞ −∞ −∞

× (x cos ϕ sin θ + y sin ϕ sin θ + z cos θ − ρ) dxdydz (2)

123

The projection theorem (3) implies that in the class { f } we can use the inversion formula: 1 f (r ) = 4π 2



∞ dρ

|λ|e−iχ r,η dλ

 

eiχ η P(η, ρ)dη

0

(5) Equation (5) describes the reconstruction algorithm of a mathematical model where { f } is the class of analyzed objects and P(η, ρ) is the array of experimental data. Now let us discuss the physical assumptions that form the basis for the model described above. Let us have a point radiation source M and a point detector D. The measurements of P(η, ρ) corresponds to the case when both the source and the detector lie on the line (η, ρ) with the radiation propagating along this line. The above model corresponds to an infinite number of measurements (all (η, ρ)) performed without error by a point detector. It is clear that these conditions are not satisfied in practice, since there are no error-free point detectors and many measurements should require an infinite time to become errors irrelevant. We can consider noise or instrumental induced artifacts assuming the quantities [12]:  (6) Ii, j = Θ(ρi − ρ)P(η j , ρ)dρ where Θ(ρ) is defined as the instrumental function (IF), which at the moment we do not specify. Generally, the IF possesses the following property Θ(ρ) = Θ(−ρ). In addition, it is defined locally, i.e., Θ(ρ) = 0 for |ρ| > ρΘ , where ρΘ is a threshold value. Analogously, the noise can be considered taking the quantities: Ii, j = Ii, j + Ii, j

(7)

SIViP

where Ii, j are random variables. The relation (7) can be easily extended to a 3D case. For a more detailed description on how to treat the noise sources, see “Appendix”. 2.1 Inverse problem associated to RT Let us consider the bi-dimensional case of reconstruction of a real image. The unavailable image that is being solved for, f (x, y) is an N × N image that can be represented by an N 2 × 1 vector, Φ: Φ = vec( f (x, y)) where ⎛ f (x1 , y1 ) . . . ⎜ .. .. f (x, y) = ⎝ . .

⎞ f (x N , y1 ) ⎟ .. ⎠ . f (x1 , y N ) . . . f (x N , y N )

(8)

The vec( f (x, y) is a column by column concatenation of f (x, y), where the total pixels in the image is represented by the number N × N . Similarly, the sinogram data R f (η, ρ) are represented by Γ : Γ = vec(R f (η, ρ)) where ⎛ f (η1 , ρ1 ) . . . ⎜ . .. R f (η, ρ) = ⎝ .. .

⎞ f (η M , ρ1 ) ⎟ .. ⎠ .

erally, we need to combine a statistical model for image degradation with a priori information. A Markov chain Monte Carlo (MCMC) method for computing the optimal estimate of the underlying real structure to be reconstructed is adopted. The key problem in removing the instrumental artifacts and noisy sources during image reconstruction is summarized in Eqs. (8) and (9). Many authors have contributed to develop useful image restoration methods using a maximum likelihood approach [9,13–19]. Our basic original contribution is that we apply a double MCMC approach, firstly, to remove the noise and instrumental aberrations, and secondly, to reconstruct the image, after artifacts have been removed. In this second case, we can demonstrate that, by a general point of view, the maximum likelihood can be computed incorporating the prior information in the estimation process of the reconstructed image. The reconstructed image Φ given a set of structure parameters (ω = ω1 , ω2 , ω3 ) is the maximum a posteriori (MAP) estimate:

    arg min P Φ, Γ, ω a = arg min p Γ Φ p Φ ω p ω a (Φ,ω)

(9)

(Φ,ω)

   = arg min log p Γ Φ + log p Φ ω + log p ω a (Φ,ω)

(10)

f (η1 , ρ L ) . . . f (η M , ρ L )

where the sinogram data are an M × L image, where each element in the vector represents one pixel and the value associated with that pixel is the sum of f (x, y) along the line ρ = x cos θ + y sin θ . In absence of noise, the reconstruction problem involves solving the following inverse problem: Γ = H Φ, where H is the forward transformation operator. The problem can be solved by inverting the transformation Φ = H −1 Γ . Such inversion presents various difficulties. The forward transformation operator is often singular or nearly singular, hence, for a given set of projections data there may be many or even no solutions, and the inverse of H may not exist. Also in the case where H −1 exists, the matrix can be too large to compute. For the reconstruction algorithms, a Bayesian approach can be used. In the Bayesian approach, an important role is played by the prior knowledge of some physiological parameters, for example in the CT tomography, the attenuation coefficients for the various organs [9,13–17]. A reconstruction algorithm based on Bayesian approach is detailed in Sect. 3.

3 Image reconstruction under low SNR conditions using a Bayesian approach In this section, we describe a Bayesian approach for improving a low SNR image-based reconstruction method. Gen-

where a = (a1 , a2 , a3 ) is a set of adjustable parameters and the joint probability can be written as [10]: p (Γ |Φ) =

N  i=1

√

1 2π σΦ

(Γi −(H Φ)i )2

e

2σ 2 Φ

where σΦ is the SD of the image Φ,   N  3 (Γi −μω )2  1 2 p (Φ|ω) = e 2σ ω √ 2π σω i=1 ω=1

(11)

(12)

and finally: p (ω|a) =

N =3 

p (ωi |ai )

(13)

i=1

Similarly, in the first approach of MCMC method, we have to consider that we have collected a signal I as given by Eq. (8), and ψin , ψout denote, respectively, the Poisson counts of photon fluorescence inside and outside the object to be measured, Π is an unidentified function taking in account the instrumental aberrations, so that the problem can be addressed to find the posterior distribution p(I |Iobs , ψin , ψout , Π ). Note that the SNR coefficient is √ expressed as SNR = (ψin − ψout )/ ψout . The Bayes’ theorem declares that the posterior distribution can be expressed as a normalized product of two terms:

123

SIViP Fig. 1 Example of 3D reconstruction from a stack sequence using a double MCMC approach. It is relevant the identification of components in the right figure

p (I |Iobs , ψin , ψout , Π ) ∝ p (Iobs |I, ψin , ψout , Π ) p(I ) (14) where the first term on the right side of Eq. (14) is the likelihood of observing Iobs given I and the second term is the prior probability of I . Because the spin counts in any voxel are constitute by Poisson random variables, the likelihood of the observed array of counts can be written as: p (Iobs |I, ψin , ψout , Π ) =

 e−ψi ψ n i i

ni !

i

(15)

which log-likelihood expression is:  log p (Iobs |I, ψin , ψout , Π ) = (n i logψi − ψi ) + const i

(16) Following Fudenberg and Paninski [9], we can introduce a limited maximum likelihood to simulate possible contribution of the Π unknown function, so that we have to compute an estimator as  (17) Eˆ = arg max (n i logψi − ψi ) − λ f (Π ) i

where λ and f (Π ) are a parameter and a function, respectively, taking in account the incidence of instrumental aberrations. Because λ and f (Π ) represents unknown contributions, we can quantify the uncertainty in the estimator Eq. (14) choosing suitable values for such quantities. One standard way to do this is to draw samples from the posterior p(I | Iobs ) and then to quantify the variability of these samples. A Gibbs sampler can be used for doing this quantification. The idea is that, given the observed data Iobs , we iteratively choose a pixel i randomly and a set of states of pixel reproducing the incidence of artifacts, for example, discontinuities for Shepp–Logan phantoms or suitable different extinction factor for real medical images [9,13,19]. The global procedure for image reconstruction from stacks is displayed in Fig. 1, and the case study is a nucleolus of a osteoblast cell imaged with confocal microscopy. A sinogram of single stacks is made both for an as-recorded image, both for Poisson-like or Gaussian noise removed image, and

123

then the IRT procedure is applied to the sinograms, for further details, see “Appendix”. The final result of the 3D reconstruction is summarized in Fig. 1: on the left, we have an image reconstructed from a stack sequence of an optical noisy low SNR nucleolus recorded using a confocal microscopy. Its sinogram is simulated applying a RT function and then is denoised using the algorithm developed in “Appendix”, [10]. Finally, 3D image reconstruction is made using the MCMC approach as described in this section (central image). Special features corresponding to inner nucleoli are enhanced with the same approach removing some possible discontinuities. No instrumental artifacts were taken in consideration, because specific information on the standard performance of the instrument used was not available. In Fig. 1, on the left, we have an optical low SNR nucleolus image from confocal microscopy. Its sinogram is simulated applying a RT function and then is denoised using the algorithm developed in “Appendix”. Finally, 3D image reconstruction is made using the MCMC approach as described in this section, (central image). Special features corresponding to inner nucleoli are enhanced with the same approach removing some possible discontinuities (right image). One general request for the post-processing 3D reconstruction in biological systems is the possibility of a suitable identification of sub-organs and components. The MCMC approach adopted in this paper, based on the Gibbs sampling, is rather feasible and meets this general request. In addition, the most relevant aspect of our approach is that, almost in principle, the 3D reconstruction can be applied to any tomographic technique involving a RT for a 3D reconstruction. The next efforts will be addressed to use our 3D reconstruction method to fMRI, as well as X-ray tomography. In the next section, we compare the performances of our procedure with the FDK method. 4 Experiments In this section, we compare two methods for the 3D Radon reconstruction: the FDK method and our method (OM) based on the MCMC algorithm as described before. The FDK method [20] describes an approximate reconstruction algo-

SIViP Fig. 2 Comparison of the two reconstruction methods for 64 slices. From the left to right on the top row, the FDK application to reconstruction (PSNR = 25.6 dB), background subtraction and edge identification. The bottom row shows the analogous results with OM, in this case, the reconstruction image presents PSNR = 26.2 dB

rithm for circular cone-beam tomography. It is an approximate method in the sense that the reconstruction result will deviate somewhat from the measured object regardless of the measurement resolution. For moderate cone-angles, these differences are, however, small and often acceptable. The simplicity of the FDK method has made it the most used algorithm for cone-beam reconstruction. Another important advantage, is that FDK, in contrast to the exact methods, handles data truncated in the longitudinal direction, so making such useful for a reconstruction from a stack imaging [21,22]. A confocal microscope 64 slices image of size 523×471 pixels of a nucleolus (Fig. 1) has been considered for training and evaluation of the two methods for the 3D Radon reconstruction. All simulations are run by Matlab on Intel® CoreTM i72600 CPU(@) 3.40 GHz machine with 16 GB RAM running 64 Bit Ubuntu Linux 13.04. The Matlab code for the FDK algorithm has been modified starting from Faridani [23]. Precision of the reconstruction was assessed by comparing the original image with the reconstructed result, using a conventional peak SNR (PSNR) measurement. The reconstruction process was performed with the idea that the series of stack recorded with confocal microscope represents a cross section of the object that is equivalent by taking a series of photographs at a number of angles around the nucleolus. Or analogously, the light reflected from the object acts as a source of a tomographic process. So, the series of the corresponding profile of the image taken at a number of angles around the object can be used to construct the cross-sectional image of the nucleolus. The volumetric data are then constructed from a stack of cross-sectional image of all horizontal line. Finally, the nucleolus is then extracted and 3D rendered. Further

Table 1 Reconstruction quality in decibels for the two tested methods at different resolution levels denoted by the number of slices used for the reconstruction Slice

16

32

64

Method

FDK

OM

FDK

OM

FDK

OM

PSNR

16.8

16.5

20.3

20.0

25.6

26.2

Table 2 Reconstruction times in seconds for the two tested methods at different resolution levels denoted by the number of slices used for the reconstruction Slice

16

32

64

Method

FDK

OM

FDK

OM

FDK

OM

Time (s)

2.2

2.2

5.3

4.4

10.1

9.0

operations can be made on the object of interest, such as background subtraction and edges identifications. Denoising procedures of Gaussian and Poissonian noise sources are performed on the single slice and such operations are detailed in “Appendix”. In Fig. 2, we present the reconstruction results obtained with FDK and our method (OM), respectively. The results are rather equivalent in terms of the PSNR parameter, but OM is slightly less time expensive. In turn, we resume the PSNR data correspondent to the reconstruction procedure for the two tested methods at different reconstruction levels, i.e., using 16, 32, 64 slices, Table 1. And the correspondent reconstruction time (s), Table 2. It is remarkable that OM requires less time to compute the reconstruction. The main advantage in time computation of OM is due to avoidance of MCMC algorithms by taking advantage of the conditional conjugated structure of

123

SIViP

Gaussian distributions given the partitions. In addition, the Gaussian distributions offer also the possibility to match the subsequent slices so that the convergence of chains during the reconstruction procedure is more robust. But the basic advantage is that our procedure is that can be applied to any tomographic technique as well to confocal microscope imaging making equivalent the tomography reconstruction to 3D reconstruction from a stack recorded imaging technique, such as confocal microscopy. 5 Conclusions Radon Transform is the basic tool for the image reconstruction in many medical diagnostic devices, such as CT, fMRI or SPECT, etc. In any of such applications, the capability to discriminate the noise or instrumental known artifacts is a main target to be reached. In this paper, we applied some algorithms derived by a double Markov chain Monte Carlo approach for a 3D reconstruction of a biological sample recorded with a confocal microscope. Poissonian and Gaussian noises were treated before applying 3D reconstruction procedure. The reconstruction method is rather general, and it can be applied to all tomographic methods using Radon Transform for a further improvement of such image-based diagnostic techniques. The reconstruction procedure of our method requires less time to be computed if compared to the FDK method. The main advantage in time computation is due to avoidance of MCMC algorithms by taking advantage of the conditional conjugated structure of Gaussian distributions given the partitions. In addition, the Gaussian distributions offer also the possibility to match the subsequent slices so that the convergence of chains during the reconstruction procedure is more robust. But the basic advantage is that, almost in principle, our procedure is that can be applied to any tomographic technique as well to confocal microscope imaging making equivalent the tomography reconstruction to 3D reconstruction from a stack recorded imaging technique. In turn, the reconstruction procedure takes advantage by a denoising approach both in the case of Gaussain and/or Poissonian source noise, with computational results comparable or better to Bayesian CRP algorithms. Acknowledgments The authors like to thank the NanoICT laboratory for useful support and S. Danti, from the University of Pisa, for the confocal images of the fibroblast cells.

Appendix: Identification of noise and denoising procedure The tomography images may suffer from different kinds of noise and artifacts. The most significant source of noise

123

occurs in the photon detection process. Artifacts can occur from the physics of the system, from partial volume effects or from patient artifacts [1,13]. Image noise occurs with the statistical fluctuation of photon detection. Generally, this type of noise can be modeled in two main ways, as Poissonian noise or Gaussian noise. The general form of a Poisson distribution can be written as: P(y) =

λ y e−λ y!

(18)

where P(y) is the probability of a random variable y, having mean λ. The variance of the distribution is proportional to the average number of photons since photons are not emitted or detected uniformly. Poisson distribution reveals important characteristics, the mean and the variance are identical, as a consequence, the signal-to-noise ratio (SNR) is proportional to the square root of the detector measurements counts producing a much better SNR. A Gaussian distribution is given as follows: P(y) = √

1 2π σ

e

y− y˜ 2 2σ 2

(19)

where σ is the variance. In the case of a Gaussian distribution, the number of photon counts is not proportional to variance. On the contrary, image artifacts have different sources [10,13,19]. The procedure for removing noise sources is summarized as follows. Let S, X, Y, N be, respectively, the matrices corresponding to biomedical source, observed signal, output 2D image (the approach is valid for all the tomography techniques, DT, MR, PET or SPECT) and noise. Then, the RT and the approximation of its inversion set following relationships between matrices X = RT(S) + N and Y = IRT(X ), where IRT denotes the inverse RT. The noise can be correlated or not. Generally, Y = IRT(RT(S) + N ) ≈ S + IRT(N ), and the noise for the image Y − S = IRT(N ) is strongly correlated, this can be a problem for Poissonian noise. To develop a denoising algorithm independent by the noise type, we apply a filter v to be specified so that the denoised image can be written as S = IRT(v(RT(X ))). A Shepp–Logan phantom (SLP) and a nucleolus as in Fig. 3 were used to test the algorithms for removing the noise. The signal-to-noise ratio of the noisy sinogram, in relation to the original sinogram, denoted as ξtheo , can be calculated as follows: ηmax 180 2 ρ=0 |P (η(ρ))| η=0 ξtheo = 10 log10 (20) σ2 where P(η, ρ) is the original sinogram. The sinograms of the Shepp–Logan and the nucleolus were generated using 64 projections at the angles varying from 0◦ −180◦ . Initial results show that the denoising approaches feel again the incidence of the noise and the SNR. It is evidenced

SIViP Fig. 3 On the left, 256 × 256 Shepp–Logan phantom. On the middle image, the filtered backprojection, with the variance of the noise as σ 2 = 0.001 and PSNR = 9.11. On the right, the filtered backprojection, with the variance of the noise as σ 2 = 0.05 and PSNR = 14.15

Fig. 4 On the left, the nucleolus as in Fig. 1 with Poisson noise. On the middle image, the filtered backprojection, with the variance of the noise as σ 2 = 0.01 and PSNR = 19.11. On the right, the filtered backprojection, with the variance of the noise as σ 2 = 0.05 and PSNR = 24.25

that in the Shepp–Logan phantom, Fig. 3, some residual noise effects are still present. In the case of the biological sample, the nucleolus considered in the preview section, and between the middle and the right images no relevant difference is outlined, Fig. 4. Denoising sinograms can be improved also by incorporating prior knowledge of the tissues being scanned. This is currently our next research direction for reconstruction algorithms. The comparison between our denoising approach and other methods is made considering separately Gaussian and Poisson noise sources. The performances are compared both visually and numerically. The common criteria used are as follows: the mean absolute errors (MAEs), mean squared errors (MSE) and the Hausdorff distance (HD), [24]. All simulations are run by Matlab on Intel® CoreTM i7-2600 CPU(@) 3.40 GHz machine with 16 GB RAM running 64 Bit Ubuntu Linux 13.04. Gaussian denoising existing methods considered in this appendix are the Bayesian CRP (BRCP), [25], and the nonparametric Bayesian dictionary learning, (BPFA) approach and running median (RM) [26]. Here, we report the numerical value of Gaussian denoising methods. Table 3 shows that our denoising method has performances very close with the BCRP method, with less timeconsuming. In addition, we have a supplemental indication of the goodness of our method if we compare MAE, MSE and HD for a noise level of σ = 0.1. Table 4 confirms that our method, when applied to Gaussian noise source, has comparable performances with the method BCRP as proposed by Li and Ghoshal, with

Table 3 Numerical comparison of denoising methods for Gaussian noise source for a Shepp–Logan 256 × 256 phantom, as in Fig. 3, for different PSNR in terms of MSE(×10−2 ) Method

PSNR

Time (s)

19.11

13.27

5.74

2.95

BCRP

0.09

0.16

0.34

0.38

238.1

BPFA

0.10

0.18

0.39

0.45

421.7

RM

0.14

0.20

0.49

0.61

80.0

OM

0.09

0.15

0.32

0.36

183.5

The PNSR is calculated using the Eq. (20). The running times for the various noise levels are very close one each other, so, indicative a time averaged on all the different running times has been reported. The number of running simulations was fixed indicatively in 100. OM indicated our method Table 4 Numerical comparison of Gaussian denoising methods for the Shepp–Logan 256 × 256 phantom when noise SD is σ = 0.1, and the PSNR = 19.11 Method

MAE

MSE

HD

Time (s)

BCRP

1.31

0.03

40.11

52.2

BPFA

1.34

0.04

41.23

64.5

RM

2.08

0.08

62.70

11.7

OM

1.35

0.05

37.78

43.5

less time to compute the estimates. In analogy to BCRP, the advantage in time computation of our procedure is due to avoidance of MCMC algorithm by taking advantage of the conditional conjugated structure of Gaussian distributions given the partitions.

123

SIViP Table 5 Numerical comparison of denoising methods for Poissonian noise source for a Shepp–Logan 256×256 phantom for different PSNR in terms of MSE(×10−2 ) Method

PSNR

Time (s)

19.11

13.27

5.74

2.95

PURE-LET

0.42

0.46

0.74

1.12

279.8

OM

0.48

0.45

0.89

1.41

253.0

The PNSR considered are the same as in Table 3, and they are calculated using the Eq. (20) and varying opportunely λ. The number of running simulations was fixed indicatively in 100 Table 6 Numerical comparison of Poissonian denoising methods for the Shepp–Logan 256 × 256 phantom, when PSNR = 19.11 Method

MAE

MSE

HD

Time (s)

PURE-LET

1.29

0.7

70.54

72.6

OM

1.25

0.7

77.09

63.4

Finally, we present the analogous of Tables 3 and 4 for the Poissonian noise source varying the parameter λ in Eqs. (18) and (19) substituting σ with λ. The comparison is made with the PURE-LET algorithm [26,27] (Table 5). As in the case of Gaussian noise, our denoising method has performances very close with the PURE-LET method, with less time-consuming. In addition, we have a supplemental indication of the goodness of our method if we compare MAE, MSE and HD for a noise level denoted by PSNR = 19.11. Table 6 confirms that our method, when applied to Poissonian noise source, has comparable performance with the PURE-LET method, with some advantages inherent the time to compute the estimates. Future efforts will be devoted to treat cases where the Gaussian and Poissonian noise sources are mixed one each other.

References 1. Shung, K.K., Smith, M.B., Tsui, B.: Principles of Medical Imaging. Academic Press, San Diego (1992) 2. Al-Kofahi, K.A., Lasek, S., Szarowski, D.H., Pacde, C.J., Nagy, G., Turner, J.N.: Rapid automated three-dimensional tracing of neurons from confocal image stacks. IEEE Trans. Inf. Technol. Biomed. 6(2), 171–187 (2002) 3. He, W., Hamilton, T.A., Cohen, A.R., Holmes, T.J., Pace, C., Szarowski, D.H.: Automated tracing of neurons in confocal and brightfield images. Microsc. Microanal. 9(4), 296–310 (2003) 4. Colonna, F., Easley, G.: Generalized discrete Radon transforms and their use in the ridgelet transform. J. Math. Imaging Vis. 23, 145–165 (2005) 5. Yu, H.P., Liu, Y.X., Liu, Z.B.: Auto detection of wood texture orientation by Radon transform. J. For. Res. 16, 1–4 (2005) 6. Jadhav, D.V., Holambe, R.S.: Feature extraction using Radon and wavelet transforms with applications to face recognition. Neurocomputing 72, 1951–1959 (2009)

123

7. Raffei, A.F.M., Asmuni, H., Hassan, R., Othman, R.M.: Feature extraction for different distances of visible reflection iris using multiscale sparse representation of local Radon transform. Pattern Recogn. 46, 2622–2633 (2013) 8. Hasegawa, M., Tabbone, S.: Amplitude-only log Radon transform for geometric invariant shape descriptor. Pattern Recogn. 47, 643– 658 (2014) 9. Fudenberg, G., Paninski, L.: Bayesian image recovery for dendritic structures under low signal-to-noise conditions. IEEE Trans. Image Process. 18(3), 471–482 (2009) 10. D’Acunto, M., Benassi, A., Moroni, D., Salvetti, O.: Radon transform: image reconstruction and identification of noise and instrumental artifacts. In: Proceedings of IEEE 22nd Signal Processing and Communications Applications Conference, SIU 2014, pp. 2280–2284 (2014) 11. Olafsson, G., Quinto, E.T.: The Radon transform, inverse problems and topography. In: Proceedings Symposia in Applied Mathematics, vol. 63 (2005) 12. Popov, D.A., Sokolova, E.B., Sushko, D.: Mathematical models in Radon tomography. In: Gindikin, S. (ed.) Applied Problems of Radon Transform, vol. 162, no. 2, pp. 129–207. AMS Translations (2005) 13. Buzug, T.M: Computed Tomography: From Photon Statistics to Modern Cone-Beam CT. Springer, Leipzig (2008) 14. Conchello, J., McNally, J.: Fast regularization technique for expectation maximization algorithm for optical sectioning microscopy. SPIE Proc. 2655, 199–208 (1996) 15. Richardson, W.H.: Bayesian-based iterative method for image restoration. J. Opt. Soc. Am. A 2(1), 55–59 (1972) 16. Markham, J., Conchello, J.: Fast maximum-likelihood image restoring algorithms for three-dimensional fluorescence microscopy. J. Opt. Soc. Am. A 18, 1062–1071 (2001) 17. Chantas, G., Galatsanos, N.P., Molina, R., Katsaggelos, A.K.: Variational Bayesian image restoration with a product of spatially weighted total variation image priors. IEEE Trans. Image Proc. 19(2), 351–362 (2010) 18. Coupé, P., Menz, M., Manjón, J.V., Ruthazer, E.S., Collins, D.L.: A CANDLE for a deeper in vivo insight. Med. Image Anal. 16(4), 849–864 (2012) 19. Yu, H., Zhao, S., Wang, G.: A differentiable Shepp–Logan phantom and its applications in exact cone-beam CT. Phys. Med. Biol. 50, 5583–5595 (2005) 20. Feldkamp, L.A., Davis, L.C., Kress, J.W.: Practical cone-beam algorithm. J. Opt. Soc. Am. 1, 612–619 (1984) 21. Bharkhada, D., Yu, H., Zeng, K., Bai, E.-W., Wang, G.: A comparative study on interpolation methods for controlled cardiac CT. Int. J. Imaging Syst. Technol. 17, 91–98 (2007) 22. Gomez, F., Santa Marta, C., Romero, E.: Multiresolution reconstruction for cone-beam tomography from raw data projections using 3D ridgelets. J. Digital Imaging 24, 1087–1095 (2011) 23. http://www.mathworks.com/matlabcentral/fileexchange/ 6032-iconebeam; see also http://people.oregonstate.edu/ ~faridana/preprints/preprints.html 24. Li, M., Ghosal, S.: Bayesian multiscale smoothing of gaussian noised images. Bayesian Anal. 9, 1–26 (2014) 25. Zhou, M., Chen, H., Paisley, J., Ren, L., Sapiro, G., Carin, L.: Non parametric Bayesin dictionary learning for sparse image representation. IEEE Trans. Image Process. 21, 13–144 (2012) 26. Luisier, F., Blu, T., Unser, M.: Image denoising in mixed poissongaussian noise. IEEE Trans. Image Process. 20, 696–708 (2011) 27. Le Montagner, Y., Angelini, E.D., Olivo-Marin, J.C.: An unbiased risk estimator for image denoising in the presence of mixed Poisson–Gaussian noise. IEEE Trans. Image Process. 23, 1255– 1268 (2014)