A mathematical model and iterative inversion for fluorescent optical projection tomography

Physics in Medicine & Biology PAPER • OPEN ACCESS A mathematical model and iterative inversion for fluorescent optical projection tomography To cite this article: Ville Koljonen et al 2019 Phys. Med. Biol. 64 045017 View the article online for updates and enhancements. This content was downloaded from IP address 130.230.187.150 on 02/09/2019 at 08:18 © 2019 Institute of Physics and Engineering in Medicine 1. Introduction The classical x-ray computed tomography (CT) has become a standard of tomographic imaging. The method allows for high accuracy results in applications ranging from the industrial to the clinical setting. However, the high energy x-rays are harmful to biological targets. In this and other regards, optical projection tomography (OPT) (Sharpe et al 2002) provides a suitable alternative in the mesoscopic scale with samples of 1 – 10 mm in diameter (Figueiras et al 2014, Soto et al 2016, Belay et al 2017). Optical imaging provides naturally better contrast with soft tissues and, at lower intensities, light as non-ionizing radiation should have a low, even negligible, impact on living organisms in imaging phase. This allows for long-term monitoring of the samples with OPT if the sample preparation and handling is properly designed. On the other hand, OPT differs from x-ray CT with respect to signal quality and optical properties. Firstly, the relatively low energy visible light used in OPT allows much less signal penetration capability than the high energy x-ray CT. Secondly, a light beam does not propagate exactly linearly, but rather it forms a beam with varying intensity and width along its path (Trull et al 2017). Thirdly, because of isotropic emission from within the sample, in the fluorescent mode a single set of parallel beams is not enough to model the light propagation between the transmitter and the sensor. In this study, we focus on the third point and strive to enhance the inversion of the fluorescent projection images. As in transmission OPT, the sample is imaged from multiple different directions using a high-resolution camera. The light visible in the images is emitted by fluorescent molecules, which act like point sources activated using excitation light. The inverse problem in fluorescent OPT is to estimate an emission function describing the sources inside the target from the projection images. This is different from transmission tomography, where V Koljonen et al A mathematical model and iterative inversion for fluorescent optical projection tomography Printed in the UK 045017 PHMBA7 © 2019 Institute of Physics and Engineering in Medicine 64 Phys. Med. Biol. PMB 1361-6560 10.1088/1361-6560/aafd63 4 1 15 Physics in Medicine & Biology IOP 18 February 2019 A mathematical model and iterative inversion for fluorescent optical projection tomography Ville Koljonen1,2,4, Olli Koskela2, Toni Montonen2, Atena Rezaei1, Birhanu Belay2, Edite Figueiras3, Jari Hyttinen2 and Sampsa Pursiainen1 1 Faculty of Information Technology and Communication Sciences, Tampere University, Tampere, Finland 2 BioMediTech Institute and Faculty of Medicine and Health Technology, Tampere University, Tampere, Finland 3 Champalimaud Research, Champalimaud Foundation, Lisbon, Portugal 4 Author to whom any correspondence should be addressed. E-mail: v[email protected] Keywords: optical projection tomography, fluorescence tomography, beam modelling, weighted Radon transform, iterative reconstruction Abstract Solving the fluorophore distribution in a tomographic setting has been difficult because of the lack of physically meaningful and computationally applicable propagation models. This study concentrates on the direct modelling of fluorescence signals in optical projection tomography (OPT), and on the corresponding inverse problem. The reconstruction problem is solved using emission projections corresponding to a series of rotational imaging positions of the sample. Similarly to the bright field OPT bearing resemblance with the transmission x-ray computed tomography, the fluorescent mode OPT is analogous to x-ray fluorescence tomography (XFCT). As an improved direct model for the fluorescent OPT, we derive a weighted Radon transform based on the XFCT literature. Moreover, we propose a simple and fast iteration scheme for the slice-wise reconstruction of the sample. The developed methods are applied in both numerical experiments and inversion of fluorescent OPT data from a zebrafish embryo. The results demonstrate the importance of propagation modelling and our analysis provides a flexible modelling framework for fluorescent OPT that can easily be modified to adapt to different imaging setups. PAPER 2019 Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. RECEIVED 20 September 2018 REVISED 21 December 2018 ACCEPTED FOR PUBLICATION 10 January 2019 PUBLISHED 18 February 2019 OPEN ACCESS https://doi.org/10.1088/1361-6560/aafd63 Phys. Med. Biol. 64 (2019) 045017 (15pp) 2 V Koljonen et al the objective function describes the attenuation of the transmitted light via absorption and scattering instead. Attenuation, however, is also present in fluorescence measurements with both the excitation and emission light. At the moment, the filtered backprojection (Natterer 1986, Kak and Slaney 1988), or FBP, is the predominant inversion technique applied in fluorescent OPT studies. For its relatively low computational cost, it is practically the only method of choice in studies, such as Bassi et al (2015) and Belay et al (2017), involving full-resolution (for example 2048-by-2048) projection images. In its discrete form, the FBP can achieve highly accurate slice reconstructions with the main limiting factor being the number of available views on the sample. In the FBP method, the projection rays are assumed to be parallel and to penetrate through the whole sample for each angle. These assumptions do not hold in the fluorescent mode of the OPT because of the distinction between excitation and emission light. Previous work by Darrell et al (2008) and Trull et al (2017) consider the defocus aberration and the beam shape, respectively, in an attempt to solve the issue. Here, we derive a weighted Radon transform, similar to those used under the far-field approximation in single photon emission computed tomography (SPECT), positron emission tomography (PET) (Chang 1978, Natterer 1986) and x-ray fluorescence tomography (XFCT) (Hogan et al 1991, Miqueles and De Pierro 2011), as an improved direct model for the fluorescent OPT. Additionally, we extend the XFCT iterative reconstruction method introduced in Miqueles and De Pierro (2011) to fluorescent OPT. The iteration makes use of the FBP to achieve sufficiently short computation times. In this paper, we aim to improve the capability of fluorescent OPT in applications demanding high accuracy. The measurement setup and the relevant physical concepts are reviewed briefly and interpreted in a mathematical way. The direct model and the iterative reconstruction procedure are introduced along with specific issues and their solutions relating to the performed measurements. We benchmark the iteration against the FBP and related procedures by comparing the reconstructions of both synthetic and zebrafish embryo data to demonstrate the achieved improvements. 2. Methods 2.1. The fluorescent OPT measurement Slice-based tomographic reconstruction, which we will adopt, imposes several assumptions on the shape and other properties of the OPT sample. Firstly, the slice itself is necessarily a bounded and convex set S⊂R2 (see figure 1) to enable the use of the 2D transforms related to tomographic modelling. This is an idealisation of the discrete measurement. Secondly, the material in the sample interacts with incident light via absorptive and scattering processes. This decreases the intensity I of the light and the overall effect is considered to be attenuation. For beams of light, the attenuation phenomenon is governed by the Beer–Lambert law that expresses changes in the intensity in terms of an attenuation coefficient µ . It has the value µ=0 when there is no attenuation, and a value µ>0 when there is. For a combination of these possibilities, we define the continuous attenuation function µ:R2→R with compact support contained in S. It is worth noticing that the attenuation coefficient is also a function of the energy carried by the propagating light. In fluorescent OPT measurements, however, both the excitation and emission light are nearly monochromatic, which resolves the issue. For emission tomography, it is necessary to define an emission function f that describes the additive increase of the light intensity at each point in the sample slice. The emission function has properties similar to the attenuation function and it is also defined as a continuous function f: R2→R with compact support contained in S. It can be shown now that f (and µ as well) is Lebesgue integrable over a continuous curve C⊂R2 , that is f∈L1 ( C) . A common approach to modelling light propagation in tomography relies on the concept of a beam, which is a line Cs,θ = {x∈R2|x· θ= s}⊂R2 with constant parameters s∈R and θ∈S1 . Here S1 = { θ ∈R2| θ  =1 } is the set of unit direction vectors in R2 , s is the signed perpendicular distance of the line from the origin and θ is its unit normal vector. Beams have a direction as well, and here the normal vector is always chosen so that θ ⊥=Tθwith T= 01 −10  (1) is in the direction of light propagation. These vectors used in beam characterisation are illustrated in figure 1. 2.2. Direct model construction In transmission OPT the objective is to recover the attenuation function µ from projection images. It is assumed that the light incident on the sample consists of parallel beams, meaning that they share the same direction vector. Integrating the differential Beer–Lambert law θ⊥·∇ I(x)= − µ(x)I(x ) (2) for light intensity I along the beam Cs,θ gives the solution as Phys. Med. Biol. 64 (2019) 045017 (15pp) 3 V Koljonen et al g(s,θ)=−ln  Is,θ 1 I0  =  x · θ=s µ(x)dx=(Rµ)(s,θ) , (3) where Is,θ 1 is the intensity after passing through the sample and g represents the projection data arranged into a sinogram, that is, the Radon transform R of f . Notice the logarithmic intensity dependence of the sinogram values. The same model can be derived analogously in fluorescent OPT. Interpreting the value of the emission function f as the added intensity of light in a beam through a point, the directional derivative of the light intensity I along a beam with normal θ is simply θ⊥·∇ I(x)=f(x) . (4) An assumption that the emission is isotropic, or similar in all directions, is made here. Integrating over the line Cs,θ gives the solution as g(s,θ)=Is,θ 1= x·θ=s f(x)dx=( R f)(s,θ) . (5) Here the intensity dependence is linear, rather than logarithmic as in (3). Equation (5) is the most extensively used starting point for direct models in fluorescent OPT methodology and often applied as such (e.g. Belay et al (2017)). However, it does not take into account all of the physical phenomena observed in measurements: most notably it neglects the attenuation of both the excitation and emission light. The modelling errors related to applying (5) to the fluorescent OPT measurements can be corrected by introducing a weight for the emission function. The weight should comprise of distinct parts for the excitation and emission phases and take into account that the propagation in both has an endpoint. The Radon transform modified for rays of light is the divergent beam transform (Natterer 1986) (D µ)(x,θ)= ∞ 0 µ(x+tθ⊥)dt . (6) Here the propagation direction is away from the point x∈R2 , which is true for the fluorescent emission. By the linearity of the propagation, the only change required for light travelling towards the point is to change the sign of the direction vector. Now, assuming ray propagation for both the emission and excitation phases, and a linear intensity dependence for the weight functions, the weight expressions can be derived as w µ (x,ϕ) = exp (−(Dµ)(x,−ϕ)) w λ (x,θ) = exp ( − ( D λ)(x,θ)) . (7) Here µ and ϕ are the attenuation function and normal vector for the excitation rays, and λ and θ are those for the emission rays. Knowing the directed angle α from the emission ray to the excitation ray (see figure 1), it is easily seen that Figure 1. Emission tomography beam geometry. The slice S is penetrated by several excitation beams with normal ϕ . This causes fluorescent emission towards the measurement plane, with normal θ . Phys. Med. Biol. 64 (2019) 045017 (15pp) 4 V Koljonen et al − ϕ=Pαθwith Pα=  cos(α) sin(α) − sin(α) cos(α) . (8) The total weight as a function of position x and emission direction normal θ is then w(x,θ)=exp(−(Dµ)(x,Pαθ)−(Dλ)(x,θ)) (9) as the product of the component weights. Applying it to (5) gives the following corrected light propagation model in fluorescent OPT: g (s,θ)= x·θ=s f(x)w(x,θ)dx=(Rwf)(s,θ) . (10) The integral operator Rw in (10) is known as the weighted Radon transform, and also as the attenuated Radon transform when w depends exponentially on an attenuation function. It has been used successfully in, for example, SPECT and PET, and even more importantly in XFCT, which is the closest analogous emission tomography method for fluorescent OPT. Also Darrell et al (2008) implicitly use this type of a direct model with a different choice of the weight function. This underlines the flexibility and adaptability of such a direct model to different views on fluorescent OPT phenomena. 2.3. Reconstruction methods Let R−1 denote the FBP based inverse Radon transform. Now, the iterative approach introduced in Miqueles and De Pierro (2011) to reconstruct XFCT data can be expressed as follows: f(k+1)=G(f(k))=f(k)+R −1 (d−R wf (k) ) a ,f(0) ≡0. (11) Here k∈N0 , G is a functional that maps the kth iterate onto an updated estimate and a (x)= 1 2πS 1 w(x,θ)d θ (12) is the directional average of the weight function w in the weighted Radon transform. The division by a compensates for the loss of contrast resulting from the Radon inversion by FBP, whereas the difference between the measurement data and the weighted Radon transform of the current estimate represents a correction needed to arrive at a better reconstruction. This method benefits from not relying on information about the weight function, making it applicable whenever (10) is used as the direct model. Conditions for the convergence of this iteration are investigated by Kunyansky (1992) and Miqueles and De Pierro (2013). The iteration (11) can be reduced back to two simpler inversion methods in a straightforward manner. Firstly, restricting k⩽1 and choosing a≡1 gives simply the FBP. Secondly, letting a instead be as in (12) yields f(1)=R −1 d a, (13) for the reconstructed emission function, which is, in fact, the SPECT simple correction method proposed by Chang (1978) applied in light-based tomography. This was already noticed in Miqueles and De Pierro (2011). A similar approach to fluorescent OPT reconstruction was also indirectly taken by Darrell et al (2008). On the other hand, releasing the restrictions on k gives a simpler version of the iteration (11) as f (k+1) =f (k) + R−1 (d −R w f (k) ),f (0)≡0. (14) The attenuation compensated iteration and these three methods are chosen for comparison to assess the effects of both the attenuation compensation and iterating the solution. 2.4. Relaxation In order to enhance the performance of (11), we apply relaxation: the actual iterate is obtained as a weighted mean of the previous one and its updated value, i.e. f (k+1) =(1 − r k )G(f (k) )+r k f (k) =f (k) +(1 − r k )e (k) (15) with a suitably chosen relaxation weight rk∈R . Here e(k) is the iterative correction term e (k)=G(f(k)) − f(k)=R −1 (d−R wf (k) ) a. (16) The relaxation expression reduces to normal iteration when rk = 0. The cases rk < 0 and rk > 0 are sometimes referred to as over-relaxation and under-relaxation, respectively. In the former one, the updated value is favoured more, whereas, in the latter one, the previous iterate is emphasized. The choice of rk can be fixed over all iterations, Phys. Med. Biol. 64 (2019) 045017 (15pp) 5 V Koljonen et al or it can vary between them. Employing relaxation adds control over the emission function reconstruction process, potentially helping it to arrive at a better solution. One possibility for systematically choosing the relaxation weights arises from the Aitken ∆2 method (Aitken 1927) for iteratively computing polynomial roots. A modern interpretation is as follows: approximate the next element in a scalar-valued sequence (xk)∞ k=0 , xk+1=h(xk) , by the first-order Taylor expansion xk+1=h(xk)=x+β(xk−x), (17) where β∈R and x=h(x) is the limit point of the sequence. Applying this at two consecutive known members of the sequence gives an estimate to the limit point as x =xk+( x k+1− x k) 2 x k+2− x k+1− (x k+1− x k )=xk+(∆ x k) 2 ∆2x k , (18) where ∆ is the forward difference operator for the sequence. In our case it has the value of the iterative correction e(k). Irons and Tuck (1969) reformulate this method as relaxation and generalise it for non-scalar sequences as the following iteration: r k+1=rk+(rk−1)  e (k) ,e (k−1) −e (k)   e(k − 1) − e(k)   2,r0= 1. (19) The explicit dependence on not only the previous relaxation factor but also the two previous correction values means that all the preceding iterations will affect the final result. This helps to make the iteration more robust with very simple and quick calculations. These properties have been observed in practice for example in fluidstructure interaction solver research (Küttler and Wall 2008). Algorithm 1. Simultaneous fluorescent OPT reconstruction and centre of rotation correction. The subscripts refer to position in a vector, not in an iterative sequence. Input: projection data p , weighted Radon operator Rw , mean weight function a, FBP operator R−1 Output: emission function estimate f // prepare for iteration 1 g← vector of sinograms from p 2 f← vector of slices set to zero 3 r← vector of relaxation weights set to one 4 ∆← vector of differences set to a large constant 5 g←R wf 6 p← vector of projections from g 7 for i←1 to number of iterations do 8 for j←1 to number of sinograms do // calculate the residual and differences 9 e←R −1 (gj − g  j )/a 10 ∆2←∆j−e 11 ∆j←e // calculate the relaxation weight 12 r j ←r j+( r j −1 )  ∆j , ∆ 2 /  ∆ 22 // iterate and reproject 13 fj←fj+(1−rj)∆j 14 g j←R w f j 15 end 16 end Algorithm 1 describes the steps necessary for reconstructing a 3D volume from the projection data in different directions. After initialising the variables, an iterate of the emission function sequence (15) is calculated, and then it is reprojected using the proposed direct model for fluorescent OPT. Such a process is iterated for a desired number of times. 2.5. Simulation and experimental setup Both numerical experiments and reconstruction of real data were performed to test the direct model in practice. Artificial measurement data were generated using (10) with both the Shepp–Logan phantom (Shepp and Logan 1974) and a more sparse phantom as the emission function. The latter represents a cross-section of a sample of Phys. Med. Biol. 64 (2019) 045017 (15pp) 6 V Koljonen et al individual cells, motivated by a promising application for fluorescent OPT (Belay et al 2017). See appendix and the supplementary material for synthesis details for the cell phantom. The two 512-by-512 emission functions are illustrated in figure 2. For computational simplicity, the attenuation functions were chosen to be µ=λ=1 within the inscribing circle of the phantom images, and zero elsewhere. An equispaced sequence of 400 projection angles between 0.0◦ and 359.1◦ , inclusive, was assumed as the sample for the projection angle θ . The angle from the emission beam to the excitation beam was chosen to be α=−10◦ . Finally, Gaussian noise with a standard deviation of 1.0% of the maximum value of the emission function was introduced in the simulated sinograms. The experimental data were obtained from a zebrafish embryo using an in-house OPT measurement system (Figueiras et al 2014, Soto et al 2016). Simplistically, the system consists of a light source, a sample chamber, the detection optics and a camera (see figure 3). The sample itself is submerged in water in the chamber to minimise the effects of a changing refractive index, and it is held and rotated with the help of a rod from above. The details on the light source and detection machinery are described below. The projection images were acquired using a sCMOS camera with an infinity-corrected, long working distance, 10 × objective. The fluorescent marker (eGFP, excitation and emission maxima 480 nm and 510 nm ) was excited using a 480 nm LED with bandwidth narrowed to 488 ±10 nm using a bandpass filter. The excitation light was filtered out from the acquired emission signal with a 520 ±36 nm filter. The instrumentation also included an iris diaphragmthat was not used, that is, the iris was fully open. The projection data of both samples were acquired at 0.9◦ intervals from 0.0◦ to 359.1◦ , 400 images in total. Before imaging, the sample was manually centered in the horizontal direction using symmetric rotation and along the projection direction using the tube walls as a visual aid. The bright field projection data of the zebrafish was imaged with the same setup, except for the light source and filters: a white light LED was used for transillumination and no filters were applied. For the purposes of reconstruction, the same assumptions regarding the emission-to-excitation angle α and the attenuation functions µ and λ were applied. This is unrealistic, but necessary to enable the fastest reconstruction. Reconstruction from the bright field data was beyond the scope of this work, but a better estimate to the values taken on by the attenuation functions could be taken from this mode. This is because attenuation is, in fact, the objective function in the transmission tomography problem. Thus it is possible to arrive at an experimental estimate of the weight function of the direct model and the attenuation compensation factor of the reconstruction as well. The zebrafish embryo sample used (Tg(fli:1a:eGFP)) was two days post fertilization (2dpf) and expressed eGFP in endothelium. The embryo was prepared in 1.5 % agarose (Sigma Aldrich, Finland) hydrogel with added Tricaine (Sigma Aldrich, Finland) for anesthetizing the fish. The sample was inserted in a fluorinated propylene ethylene (FEP) tube for the imaging. The FEP tube had an inside diameter of 1 nm . The sample was imaged Table 1. Reconstruction SNR values. The values’ location in the table corresponds to the image-method pairs shown in figure 5. SNR (dB) FBP Simple correction Plain iteration Relaxed iteration Shepp–Logan 0.0245 0.0830 0.0671 0.0861 Shepp–Logan ROI 0.0093 0.0651 0.0394 0.0879 Cell 0.0012 0.0051 0.0040 0.0064 Cell ROI 0.0013 0.0055 0.0043 0.0066 (a) (b) Figure 2. The Shepp–Logan (a) and cell phantoms (b) used in the numerical experiments. Phys. Med. Biol. 64 (2019) 045017 (15pp) 7 V Koljonen et al with the OPT system using a bright field exposure time of 13 ms and a fluorescence exposure time of 3000 ms , respectively. Only one focal plane, placed manually in the sample center, was used for all image sets. Examples of both bright field and fluorescence projection images are shown in figure 4. Figures 4(a)–(c) show the head of the zebrafish embryo and figures 4(d) and (e) show its tail. The measured fluorescent projection images of the embryo tail were cropped from 2048-by-2048 dimensions to 512-by-512 to achieve shorter reconstruction times. The full width of tail fit into the cropped data. For the head projections, 1024-by-1024 resolution cropped versions of the original images were used, and also average down-scaled version of 1024-by-1024 to size 512-by-512 for comparison. 2.6. Reconstruction and implementation The slice-wise reconstructions were computed with four different methods: the filtered backprojection, the correction method (13), the plain iteration (14) and the relaxed iteration (15) with relaxation factors calculated using (19). The number of iterations was five (5) in the last two methods for numerical experiments and three (3) for experimental data. Every FBP operation was computed using the iradon function from the Image Processing Toolbox. The Hamming filter compressed to 50% of the full frequency spectrum was applied. Before the actual reconstructions, the projection data was manually corrected for centre of rotation misalignment. This is a result of the non-ideal rotation of the sample, where its vertical centre axis is not perfectly aligned with the axis of rotation. This causes the centre of rotation to deviate from the geometric centre of at least some of the slices. The inversion algorithms, however, assume that these two points coincide, and this causes artifacts in the reconstructed slices. Walls et al (2005) Different ways of achieving this correction were experimented with, including the variance (Walls et al 2005, Dong et al 2013, Figueiras et al 2014), the centre of mass (Azevedo et al 1990, Donath et al 2006, Dong et al 2013, Jun and Yoon 2017) and the cross-correlation (Yang et al 2013) based methods. All of these depend on optimising an image metric value between different reconstructions, and then deducing the correct shift needed to counter the centre of rotation error. It was also tested, whether this process could be embedded within the iterative inversion procedure, following the work of Gürsoy et al (2017). Finally, the actual method used was based on the ideas presented by Figueiras et al (2014): the centre of rotation deviation was manually assessed for the top and bottom slices of the data set and interpolated for the remaining ones. The slice-wise sinograms were corrected by shifting them circularly in the coordinate direction by the amount of the deviation. The 512-by-512 2D image results of the numerical experiments are illustrated using the MATLAB imagesc function with a greyscale colourmap. The visualisation of real data is done in two ways. The simpler of these is to choose a vertically oriented cross-section (coronal plane) of the 3D reconstruction volume, and to show this in two dimensions. The more complex method is the isosurfacing method based on the assumption of related areas on the zebrafish embryo having similar emission function values. This enables the surface of the sample to be shown at least reasonably well by meshing an appropriate isosurface of the 3D reconstruction volume. The mesh was generated and illustrated respectively using the isosurface and patch functions from the Image Processing Toolbox. The emission function values were normalised to their maximum, and a range of illustrations with different isosurface values was tested manually. This visualisation approach is only practical at a relatively low resolution, and the results between reconstruction algorithms are not directly comparable, but it gives a more complete view on the 3D structure of the sample. All computations were performed using a high-end Lenovo P910 Workstation equipped with two Intel Xeon E5-2697 processors and 256 GB of RAM. The MATLAB version R2016b 64-bit (The MathWorks, Inc.) was employed as the computation platform. The functions and scripts written for the purposes of this paper are published as supplementary material. See appendix for an overview of this content. Figure 3. Visualisation of the OPT measurement setup (Figueiras et al 2014, Soto et al 2016). The imaging system has a fixed position while the sample rotates around a vertical axis. Phys. Med. Biol. 64 (2019) 045017 (15pp) 8 V Koljonen et al 3. Results 3.1. Numerical experiments The results of the numerical experiments for both phantoms are shown in figure 5. The 128-by-128 regions of interest (ROI) were selected to see the quality of the reconstruction at a smaller scale as well. For the Shepp– Logan phantom this is the central area containing the two smaller circles and for the cell phantom the bottomleft corner with the four squares of differing brightness. Table 1 shows the signal-to-noise ratio (SNR) values corresponding to each of the reconstructed phantoms. Additionally, the mean-squared error (MSE) of each of the reconstructions are shown by iteration in figure 6 to show the rates of convergence for each of the iterative methods. Differences in the reconstructions are visible between both the iterated and the non-iterated methods, and the methods with and without attenuation compensation. Firstly, the non-compensated FBP and plain iteration show significantly darker central areas in the Shepp–Logan slice (figures 5(a) and (c)). A similar phenomenon can be detected with the corresponding cell reconstructions in figures 5(i) and (k) as well but to a much lesser degree. There a more striking similar feature is the overly bright rightmost square. Moreover, in the plain iteration reconstructions, these properties seem to be weaker than in the FBP. The attenuation compensated simple correction (figures 5(b) and (j)) and relaxed iteration methods (figures 5(d) and (l)), on the other hand, show an overall level of contrast much closer to the original phantom. The above visual assessment of the phantom reconstructions is supported by the SNR values computed between the reconstructed image and the original. Here the deviation of the reconstruction from the original is included in the noise. The relaxed iteration gives the highest SNR values in all cases, whereas the FBP performs the worst. As expected the cell phantom reconstruction values are much lower than those of the Shepp–Logan reconstructions, because of the higher sparsity of the emission function. A similar trend is apparent in the MSE comparison of the reconstruction methods in figure 6. The attenuation compensated relaxed iteration shows slightly lower overall deviation from the original phantom in both cases. The same results give insight into the differences between the iterated and non-iterated methods as well. The error values of the FBP and simple correction methods coincide with the first iterate ones for the plain and relaxed iteration, respectively. All graphs in figure 6 show decreasing MSE values over the five iterations, meaning that in this sense the iteration can improve the reconstruction. The relaxed iterative solution is found quicker than the plain iterated one, as the MSE shows better convergence. The iteration can and needs to be stopped here, (a) (b) (c) (d) (e) Figure 4. The head (top row) and the tail (bottom row) of the zebrafish embryo sample. (a) and (d) Show the bright field OPT projections, whereas (b)–(c) and (e) are fluorescent images. The bright regions in (b) and (c) correspond to high emitted intensity from the yolk sac of the embryo. Contrast in all images has been enhanced for print. The scalebar corresponds to 260 µ m. Phys. Med. Biol. 64 (2019) 045017 (15pp) 15 V Koljonen et al Jari Hyttinen https://orcid.org/0000-0003-1850-3055 Sampsa Pursiainen https://orcid.org/0000-0002-9131-9070 References Aitken A C 1927 XXV.–On Bernoulli’s numerical solution of algebraic equations Proc. R. Soc. Edinburgh 46 289–305 Azevedo S G, Schneberk D J, Fitch J P and Martz H E 1990 Calculation of the rotational centers in computed tomography sinograms IEEE Trans. Nucl. Sci. 37 1525–40 Bassi A, Schmid B and Huisken J 2015 Optical tomography complements light sheet microscopy for in toto imaging of zebrafish development Development 142 1016–20 Belay B, Koivisto J T, Vuornos K, Montonen T, Koskela O, Lehti-Polojärvi M, Miettinen S, Kellomäki M, Figueiras E and Hyttinen J 2017 Optical projection tomography imaging of single cells in 3D gellan gum hydrogel EMBEC & NBC 2017 (Berlin: Springer) pp 996–9 Chang L T 1978 A method for attenuation correction in radionuclide computed tomography IEEE Trans. Nucl. Sci. 25 638–43 Darrell A, Marias K, Garofalakis A, Meyer H, Brady M and Ripoll J 2006 Accounting for Point Source Propagation Properties in 3D Fluorescence OPT. EMBS’06 (IEEE) pp 6513–6 Darrell A, Meyer H, Marias K, Brady M and Ripoll J 2008 Weighted filtered backprojection for quantitative fluorescence optical projection tomography Phys. Med. Biol. 53 3863–81 Donath T, Beckmann F and Schreyer A 2006 Automated determination of the center of rotation in tomography data JOSA A 23 1048–57 Dong D et al 2013 Automated recovery of the center of rotation in optical projection tomography in the presence of scattering IEEE J. Biomed. Health Inform. 17 198–204 Figueiras E et al 2014 Optical projection tomography as a tool for 3D imaging of hydrogels Biomed. Opt. Express 5 3443–9 Gürsoy D et al 2017 Rapid alignment of nanotomography data using joint iterative reconstruction and reprojection Sci. Rep. 7 11818 Hogan J P, Gonsalves R A and Krieger A S 1991 Fluorescent computer tomography: a model for correction of x-ray absorption IEEE Trans. Nucl. Sci. 38 1721–7 Irons B M and Tuck R C 1969 A version of the Aitken accelerator for computer iteration Int. J. Numer. Methods Eng. 1 275–7 Jun K and Yoon S 2017 Alignment solution for CT image reconstruction using fixed point and virtual rotation axis Sci. Rep. 7 41218 Kak A C and Slaney M 1988 Principles of Computerized Tomographic Imaging (Piscataway, NJ: IEEE) Koljonen V, Koskela O and Montonen T 2018 fOPT Tools (https://doi.org/10.5281/zenodo.1422207) Kunyansky L A 1992 Generalized and attenuated Radon transforms: restorative approach to the numerical inversion Inverse Problems 8 809–19 Küttler U and Wall W A 2008 Fixed-point fluid-structure interaction solvers with dynamic relaxation Comput. Mech. 43 61–72 Landweber L 1951 An iteration formula for Fredholm integral equations of the first kind Am. J. Math. 73 615–24 Miqueles E X and De Pierro A R 2011 Iterative reconstruction in x-ray fluorescence tomography based on Radon inversion IEEE Trans. Med. Imaging 30 438–50 Miqueles E X and De Pierro A R 2013 On the iterative inversion of generalized attenuated Radon transforms J. Inverse Ill-Posed Problems 21 695–712 Natterer F 1986 The Mathematics of Computerized Tomography (Philadelphia: SIAM) Piana M and Bertero M 1997 Projected Landweber method and preconditioning Inverse Problems 13 441 Sharpe J, Ahlgren U, Perry P, Hill B, Ross A, Hecksher-Sørensen J, Baldock R and Davidson D 2002 Optical projection tomography as a tool for 3D microscopy and gene expression studies Science 296 541–5 Shepp L A and Logan B F 1974 The Fourier reconstruction of a head section IEEE Trans. Nucl. Sci. 21 21–43 Soto A M, Koivisto J T, Parraga J E, Silva-Correia J, Oliveira J M, Reis R L, Kellomäki M, Hyttinen J and Figueiras E 2016 Optical projection tomography technique for image texture and mass transport studies in hydrogels based on gellan gum Langmuir 32 5173–82 Strand O N 1974 Theory and methods related to the singular-function expansion and Landweber’s iteration for integral equations of the first kind SIAM J. Numer. Anal. 11 798–825 Trull A K, van der Horst J, Palenstijn W J, van Vliet L J, van Leeuwen T and Kalkman J 2017 Point spread function based image reconstruction in optical projection tomography Phys. Med. Biol. 62 7784–97 Trull A K, van der Horst J, van Vliet L J and Kalkman J 2018 Comparison of image reconstruction techniques for optical projection tomography Appl. Opt. 57 1874–82 Walls J R, Sled J G, Sharpe J and Henkelman R M 2005 Correction of artefacts in optical projection tomography Phys. Med. Biol. 50 4645 Yang M, Pan J, Zhang J, Song S J, Meng F, Li X and Wei D 2013 Center of rotation automatic measurement for fan-beam CT system based on sinogram image features Neurocomputing 120 250–7 Phys. Med. Biol. 64 (2019) 045017 (15pp)