A mathematical model and iterative inversion for fluorescent optical projection tomography

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.

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 andSlaney 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 andDe 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.

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 ⊂ R 2 (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 µ : R 2 → 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 : R 2 → 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 ⊂ R 2 , that is f ∈ L 1 (C).
A common approach to modelling light propagation in tomography relies on the concept of a beam, which is a line C s,θ = {x ∈ R 2 | x · θ = s} ⊂ R 2 with constant parameters s ∈ R and θ ∈ S 1 . Here S 1 = {θ ∈ R 2 | θ = 1} is the set of unit direction vectors in R 2 , 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 is in the direction of light propagation. These vectors used in beam characterisation are illustrated in figure 1.

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 for light intensity I along the beam C s,θ gives the solution as where I s,θ 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 An assumption that the emission is isotropic, or similar in all directions, is made here. Integrating over the line C s,θ gives the solution as 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) Here the propagation direction is away from the point x ∈ R 2 , 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 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 −ϕ = P α θ with P α = cos(α) sin(α) − sin(α) cos(α) .
The total weight as a function of position x and emission direction normal θ is then as the product of the component weights. Applying it to (5) gives the following corrected light propagation model in fluorescent OPT: The integral operator R w 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.

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: Here k ∈ N 0 , G is a functional that maps the kth iterate onto an updated estimate and 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 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 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.

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.
with a suitably chosen relaxation weight r k ∈ R. Here e (k) is the iterative correction term The relaxation expression reduces to normal iteration when r k = 0. The cases r k < 0 and r k > 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 r k can be fixed over all iterations, 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 (x k ) ∞ k=0 , x k+1 = h(x k ), by the first-order Taylor expansion 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 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: 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 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.

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 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  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.

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)  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.

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, because the visually apparent quality of the images decreased quickly with the number of iterations. The noise seems to amplify in the reconstruction procedure, which leads to a trade-off between the convergence of the MSE values and the noise content of the final image. The preceding analysis of the overall reconstructions has a relatively definite outcome. When focusing on the details shown in the regions of interest in the reconstructions, however, the results are not as straightforwardly in favour of the relaxed iteration. In the Shepp-Logan region of interest in figures 5(e)-(h) it is seen that the noise in the synthetic data is best handled by the FBP algorithm, whereas the relaxed iteration shows relatively noisy objects. On the other hand, the more detailed views on the cell reconstructions in figures 5(m)-(o) show somewhat blurred shapes in contrast to the relaxed iteration in figure 5. Additionally, the plain iteration seems to have a slight edge over the relaxed iteration in terms of the distinctiveness of low-intensity sources. The SNR values, on the other hand, suggest that the relaxed iteration is still the best, especially because the high ratios are preserved when zooming in to the regions of interest. The most significant decreases in the SNR occur when comparing the FBP, the simple correction and the plain iteration results for the Shepp-Logan regions of interest to the overall reconstruction.

Zebrafish embryo reconstructions
The coronal planes of four different zebrafish embryo tail reconstructions, consisting of the four previously mentioned methods, are visualised in figure 7. The details shown in the images represent the different blood vessels found in the tail of the embryo. The thick one is seen as brighter near the centre and next to it the thin vessels form a dimmer ladder-like structure.

FBP
simple correction plain iteration r elaxed iteration Comparing the illustrations of the different reconstructions reveals only small differences between them. All four images show an adequately sharp image along the coronal plane with the two vertical vessels easily recognisable. The most significant differences relate to the brightness of the thicker vessel and the distinctiveness of the horizontal ones. The relaxed iteration solution (figure 7) shows overall somewhat lower intensity than especially  the non-compensated reconstructions in figures 7(a) and (c). The simply corrected result in figure 7 is located between the others in this respect. As a result, it is more difficult to distinguish the horizontal vessels in the relaxed iterative solution. It was not expected, however, that the coronal slices would show them clearly, because the plane runs between the two laddered vessel structures found in zebrafish embryos. In this sense, the relaxed iteration seems to yield the best outcome. While the iterated solution can seemingly be more accurate with respect to the fine details in the measurement data, it clearly lacks robustness to noise. Out of the four images in figure 7, the FBP and the simple correction show slightly better contrast, suggesting noise amplification in the iteration process. This observation is consistent with the numerical experiments. It was, on the other hand, also observed in the numerical experiments on the cell phantom that the attenuation compensated solutions have better balance between the interior and the exterior parts of the slices, a property that could explain the different distinctiveness of the relaxed iteration in figure 7(d).
The central vertical slices of the FBP and the relaxed iteration reconstructions of the zebrafish embryo head are shown in figure 8. Here the most striking features are the dimly recognisable curvature of the head and the bright yolk sac. Some internal features can also be seen. The reconstructions are calculated here at two resolutions. For the higher 1024-by-1024 resolution the original projection images were cropped to contain the head only, and the 512-by-512 one was down-sampled and interpolated from this.
Structurally, the reconstructions share similar properties, with the curvature at the back of the head and shapes in its centre, but they differ vastly in the brightness. The FBP reconstruction in figure 8 shows the expectedly high intensity emission from the yolk sac. Looking at the relaxed iteration solution in figure 8, there is a clear loss of brightness of the yolk sac and contrast with the background compared to the FBP. Overall, the differences between the two methods are more difficult detect here, possibly due to the higher emitter density of the head part of the zebrafish embryo. The higher resolution images in figures 8(c) and (d) are essentially lower-intensity versions of the two other images, and thus the differences are also similar.
Finally, the reconstruction time from the 512-by-512 images is just over 9 min with three iteration steps. At 1024-by-1024 resolution the time required increases to around 100 min. For comparison, the FBP reconstruction from 2048-by-2048 images takes approximately 30 min. The MSE difference of second and third iterate in relaxed iteration for the tail was 0.21 and for the head 0.21 with 512-by-512 resolution and 0.23 with 1024-by-1024 resolution, suggesting good convergence of the solution.
The isosurfacing approach to visualising the zebrafish embryo data is included in the supplementary material. Videos show the evolution of the isosurface structure as the segmentation parameter is adjusted. The isosurface videos further show the similarities in detected vessel boundaries. An additional video is showing the sequence of coronal slices through the zebrafish embryo head reconstructions. See appendix for a more detailed description.

Discussion
In this study, we discussed a simple and fast inversion technique for fluorescent OPT inspired by the work of Miqueles and De Pierro (2011) on x-ray fluorescence CT. The two-phase attenuated light propagation in the fluorescent OPT measurement is taken into account by adding a weight function to the simple Radon transform forward model. The reconstruction is achieved via an iterative use of the FBP, the classical inverse Radon transform.
In the simulated and real data reconstructions it was found that the discussed forward model and its iterative inversion can enhance the fluorescent OPT imaging outcome as compared to the plain FBP. In the numerical experiments, the performance of the simple correction, the plain iteration and the relaxed iteration were observed to vary in terms of the contrast, noise and sharpness of shapes. Similarly, the distinctiveness of the blood vessels of the tail and the overall intensity level of both the tail and the head were affected by the choice of method in the zebrafish embryo reconstructions. The discussed relaxed iterative methods presented some improvement over the FBP in the experiments. Moreover, they are simple to implement, and for these reasons we propose to replace the conventional FBP by a weighted Radon transform model and iterative inversion for a wide range of fluorescent OPT imaging situations. Along with the relaxed iterations, the plain iteration might also be preferable for low-intensity image details, such as a weakly stained cell in a sparse distribution. In general, the proposed approaches were found to allow finding at least as good reconstruction as the FBP without remarkably extending the computation time. The value of the proposed methodology is in its potential for future development as outlined below.
It was also seen that the iterative methods did not fare very well with the additive noise in the numerical experiments with the Shepp-Logan phantom. On the other hand, in the case of the cell phantom especially the relaxed iteration was observed to produce a sharper reconstruction than the others. It is possible that there is variation in the applicability of a single method between different types of samples. A similar distinction between sparse and non-sparse samples was observed by Trull et al (2018) in their comparison of non-iterated deconvolution and iterated point spread function methods.
Different, or at least further developed methods are needed to still increase the resolution of the reconstructions. Noise is an apparent issue, especially when experimenting with dimensions closer to the full 2048-by-2048 projection images. Down-sampling the measurement data uses interpolation, which helps to reduce the random noise translated into the final reconstructions. This effectively regularises the inverse problem, making it easier to estimate its solution. On the other hand, it is easy to lose structural information with too heavily applied interpolation.
The inverse iteration was observed to produce the optimal imaging outcome with a relatively low number of steps. Using more than five steps led to noisy and deteriorated results unless a very small relaxation parameter was used. The observed behaviour is an example of semi-convergence (Natterer 1986), where the iterates first seem to enhance the reconstruction quality, but then become increasingly noisy and lose structural information. While Kunyansky (1992) and later Miqueles and De Pierro (2013) find conditions for the convergence of the proposed iteration, their analyses do not consider the effects of noise in the measurement data. Thus, even though the discussed iterations may converge in the sense that the function sequence ( f (k) ) ∞ k=0 has a limit, the final iterate is not necessarily a suitable reconstruction. Therefore we suggest the use of 3-5 iterative steps.
To further improve the imaging outcome might necessitate using more advanced inversion techniques. The presented method for weighted Radon inversion is, in fact, a modified Landweber iteration, the simplest known iterative regularisation technique for linear inverse problems. The original method (Landweber 1951, Strand 1974, formulated for tomography, would make use of only the plain backprojection instead of the FBP. This is, however, known to lead to serious blurring in the reconstructed image, and far more iterative steps would be required to find acceptable results. Further, theoretically developed options specifically for the Landweber-type iterations include the use of projection methods and preconditioning, which allows incorporating a priori information in the inversion procedure (Piana and Bertero 1997). In fact, the attenuation compensation used here is already a crude form of the latter. The possibilities for improvement of the iterative reconstruction method seem promising. Another strategy to improve the inversion results would be to develop the proposed forward mapping for fluorescent OPT, that is, the weighted Radon transform. A direct model of the same form has been considered indirectly by Darrell et al (2008), but with a different choice of the weight function. Their objective is to contain the effects of defocus aberration within that function, whereas our approach considers the attenuation of the light propagating towards the fluorophores and the camera. Also, their weight function is determined exper imentally, which would also have be done with our model in applied studies. The similarity, however, suggests that the direct model in fluorescent OPT could be improved simply by developing the weight function expression. Apart from the defocus issue, features currently not taken into account include the beam width and intensity which, in OPT, are not constant but vary along the path of the beam obeying approximately the Gaussian beam model. Moreover, the attenuation function was assumed constant, which is equivalent to assuming a completely homogeneous propagation medium. The discussed method includes the inhomogeneity, but this could not be implemented because of the high computational cost. It is nevertheless likely that the weighted Radon transform, with the choice of a suitable weight function, is applicable as a direct model to a range of OPT imaging scenarios extending even beyond just the fluorescent mode. More work needs to be done on understanding the effects of a given weight function, but our simplified example acts as a starting point.
To model fluorescent OPT with the presented weighted Radon transform, where the emitted light is also assumed to consist of parallel beams, can first seem erroneous. The fluorophores inside the sample are effectively point sources, and thus a physically accurate model would describe the light propagation as a spherical wave. However, the well known far-field approximation can be applied to simplify the model. Far away from the point source, the wavefront shape approaches that of the plane wave, which is the wave-optical counterpart of the beam of light. Under this idealisation, the fluorescent source emission can be reduced back to the beam representation. At the same time, far away from the sources the mesoscopic differences in the macroscopic distance to the detector become negligible, and thus the inverse-square effect on light intensity discussed in Darrell et al (2006) can be omitted as well.
In future work, the goal is to use the presented direct-inverse model pair in applied studies, which will provide more knowledge on choosing and optimising the model parameters in practice. We will also investigate ways to improve the noise-robustness of the iterative algorithms, as well as develop a direct model for fluorescent OPT incorporating non-linear optical beam properties, such as the Gaussian shape.

Conclusion
This study focused on direct modelling and iterative inversion in fluorescent OPT. As an improved direct model for the fluorescent OPT, we presented a weighted Radon transform. This was introduced in the most general form, and an example of a simple weight function model was provided. It is not, however, necessary to restrict the modelling to our method, but rather the weighted Radon transform could be used as a modelling framework in future studies. Additionally, we discussed iterative inversion techniques to enhance the quality of the reconstruction compared to the filtered backprojection. The need for such methods arises from the different, more accurate direct model. It was found that especially the relaxed iteration among the investigated methods enhances the fluorescent OPT imaging outcome over the plain FBP. The performance of the compared methods was observed to vary with respect to the noise in the measurement data and the distinctiveness and connectedness of the blood vessels in the tail reconstructions. Due to its simplicity and decent overall performance, we therefore suggest the weighted Radon transform coupled with the relaxed iteration as the method of choice to replace FBP in most of the fluorescent OPT imaging scenarios.
Protection of Animals Used for Scientific or Educational Purposes (497/2013) and the Government Decree on the Protection of Animals Used for Scientific or Educational Purposes (564/2013).

Appendix. An overview of the supplementary material
The supplementary material fopt-tools consists of MATLAB functions and scripts for fluorescent OPT inversion using the discussed methods. The zebrafish embryo dataset along with video material of the final reconstructions is also provided. In its published form, the package is ready for computing the reconstructions and illustrations shown in this paper, with no changes required. Concise instructions are attached to the package, and a more thorough overview of the workflow is given here.
The functions were written for use with MATLAB 64-bit version R2016b or later, and they rely on the Image Processing and the Parallel Computing Toolboxes. A working installation of each of these components is required. The measurement data should be placed into a subdirectory, such as projections, as projection images. Moreover, these files should be named systematically, and the names should contain an identifier between 1 and the number of images. Additionally, the working directory has to contain sinograms and reconstructions as subdirectories.
The actual reconstruction is done by running the script in the file reconstruct.m. Before doing this the model parameters, that is the projection angles, attenuation coefficients and the angle from the emission ray to the excitation ray, need to be set to appropriate values. The image dimension, the number of iterations and options for attenuation compensation, relaxation and centre of rotation correction can also be varied between reconstructions. All of these changes are to be made by editing the script file directly.
The last step in the workflow is to visualise the reconstructions. The script file show3Dreco.m is meant for this purpose. Again, before executing the script, the file needs to be edited to show the correct reconstructions with correct image dimension and desired isosurface values.
The source files for all of the functions written for this material are located in src. The files are well-documented, and the reader is encouraged to refer to them for more information. The most crucial functions are contained in expradon.m, iterinversion.m and fluorphantom.m, which are used to calculate the weighted Radon transform in the case of constant attenuation functions, to handle the iterative reconstruction process and to construct the cell phantom used in the numerical experiments. All files are published separately Koljonen et al (2018) under the Open Source Initiative MIT license and are available at https://doi.org/10.5281/ zenodo.1422207.
The isosurfacing method of visualising the zebrafish embryo reconstructions is illustrated in the video files isosurface-tail-fbp-128.avi and isosurface-tail-iter-128.avi. Additionally, figure A1 shows a single frame of the relaxed iteration reconstruction video as an example. Here the resolution of the reconstructed 3D volume is 128-by-128-by-128. Despite this, the hollow structure of the vessels is shown quite well, but some of the horizontal vessels are not connected as they are in reality. Coronal slice sequences of the head are also provided as a video (coronal-head-both-1024.avi) in the supplementary material. Jari Hyttinen https://orcid.org/0000-0003-1850-3055 Sampsa Pursiainen https://orcid.org/0000-0002-9131-9070