1. Introduction

1.1. Rendering Equation

There are two methods to render a volume: indirect and direct volume rendering (see Figure 1). Indirect rendering is used to render the isosurface of a volume [16]. First the generation of an intermediate representation of the data set is required (e.g. a polygonal representation of an isosurface generated by Marching Cubes [17]). The second step is the rendering of this representation. In a polygonal representation only a limited number of the 3D data set directly contributes to the 2D output. The interior of the volume is not considered.

Figure 1. This illustration shows a CT scan of a bonsai. The leafs are visualized by direct volume rendering whereas the trunk and the branches are visualized using an isosurface rendering. Obviously in the direct volume rendering each voxel contributes to the resulting 2D image. (image source: [15]).

Figure 1. This illustration shows a CT scan of a bonsai. The leafs are visualized by direct volume rendering whereas the trunk and the branches are visualized using an isosurface rendering. Obviously in the direct volume rendering each voxel contributes to the resulting 2D image. (image source: [15]).

Since in biological or medical application the interior is of utmost importance, for our application we use direct rendering [18]. Direct rendering techniques do not need the generation of an intermediate representation of the data set, since every voxel contributes to the resulting 2D image. To abide to core physical properties of the rendering process, we use a model that takes emission and absorption into account. In this model every particle emits and absorbs light. The combination of all these properties is maintained in the following rendering equation, a simplified version of the equation introduced by [19]:

$$\begin{array}{cc}\hfill I\left(D\right)& {\displaystyle ={\int }_0^{D}g\left(s\right)·e^{\left({\displaystyle -{\int }_s^D\tau \left(t\right)\mathrm{d}t}\right)}\mathrm{d}s}\hfill \\ \hfill & {\displaystyle ={\int }_0^{D}g\left(s\right)·t\left(s\right)\mathrm{d}s}\hfill \end{array}$$

(1)

This equation integrates from 0, the edge of the volume, to D at the eye. $I\left(D\right)$ describes the value of accumulated light intensity of one ray traversing the volume. $\tau \left(t\right)$ is the absorption at t and $g\left(s\right)$ is the sampled value at s.

In fact, the equation consists of two main terms: $g\left(s\right)$ denotes emission and $t\left(s\right)$ the transparency at point s.

1.1.1. Emission

In general, the emission part of the used model described in Eq. 1 assumes that each particle of the volume denotes a tiny light source. Thus, even if no external light source is applied, the volume can emit a certain amount of light. This behavior is motivated by the following applications:

• Microscopy

  (a)

  Self-luminous particles are modelled. This occurs during fluorescence microscopy, as certain proteins glow.

  (b)

  The particle is illuminated by a light source. This leads to a better spatial impression.

• MRI

  (a)

  Magnetic signals are processed. Varying values are transformed into visual signals.

  (b)

  No external illumination is applied/applicable while recording data. However, in a voxel grid, the signals can be compared to self-illuminated particles.

In the used model this light is not scattered. What is more, we neglect the attenuation of light on the path from the light source to an enlightened particle and we neglect the attenuation of light from the volume to the camera as well. This is not completely physically correct but enables a simple model and thus a high computational performance. Figure 2 illustrates this simplification.

1.1.2. Absorption

In our case absorption describes the attenuation of light. Eq. 2 describes the transparency between point s and the eye. If we assume a transparency of $t\in [0..1]$ we can simply reformulate this transparency such that it becomes opacity by $\alpha \left(s\right)=1-\tau \left(s\right)$.

$$t\left(s\right)=e^{\left({\displaystyle -{\int }_s^D\tau \left(t\right)\mathrm{d}t}\right)}\mathrm{d}s$$

(2)

1.1.3. Discretization

Since Eq. 1 cannot be solved analytically for all interpolation methods, to obtain a general discretization one has to solve it numerically. This can be easily done by the Riemann sum and a fixed step size $\Delta s$. The step size scales the current sampling position along the ray. More precisely, the ray gets divided into n equal segments, each of size $\Delta s$ [19].

$$\begin{array}{cc}\hfill t\left(s\right)=e^{\left({\displaystyle -{\int }_s^D\tau \left(t\right)\mathrm{d}t}\right)}\mathrm{d}s& \approx e^{\left({\sum }_{i=i+1}^n\tau (k·\Delta t)\Delta t\right)}\hfill \\ \hfill & {\displaystyle =\prod _{j=i+1}^n{e}^{\left(\tau (k·\Delta t)\Delta t\right)}}\hfill \\ \hfill & {\displaystyle =\prod _{j=i+1}^n{t}_j}\hfill \end{array}$$

(3)

Thus, we end up with the discretized rendering equation that approximates Eq. 1:

$$\begin{array}{cc}\hfill I\left(D\right)& {\displaystyle ={\int }_0^{D}g\left(s\right)·e^{\left({\displaystyle -{\int }_s^D\tau \left(t\right)\mathrm{d}t}\right)}\mathrm{d}s}\hfill \\ \hfill & {\displaystyle \approx \sum _{i=1}^{n}g(i·\Delta s)\Delta s·\prod _{j=i+1}^n{t}_j}\hfill \\ \hfill & {\displaystyle =\sum _{i=1}^n{g}_i·\prod _{j=i+1}^n{t}_j}\hfill \end{array}$$

(4)

2. Material and Methods

2.2. Stereo Rendering

There are many possibilities for a stereo rendering setup. Some of these approaches do not produce correct results.

Toe-in is easy to implement. Since in this projection the two stereo cameras are pointing to the same focal point, a rotation of the object suffices to render the two images. Unfortunately, one suffers from vertical parallax. The higher the distance to the center of the projection plane, the more this effect occurs. Thus this approach does not result in correct stereo images. Figure 4 illustrates this rendering setup.

2.2.1. Off-Axis

Off-axis is the correct projection. No vertical parallax is introduced. Both cameras have a different focal point (see Figure 5). The viewing directions are parallel. Thus, it is necessary to render two images with different camera view frustums. However, the two extended camera frustums do not overlap completely. To obtain the off-axis projection plane of both images, one has to trim the projection of the extended frustums. Therefore, one has to compute the non-overlapping amount of pixels $\delta$ [26]:

$$\delta =\frac{b·w}{2·f_o·tan\left(\frac{\alpha }{2}\right)},$$

(8)

where w is the image width in pixel and $f_o$ is the focal length. b is the stereo base, i.e. half of camera x-offset. The angle of view $\alpha$ can be computed as follows:

$$\begin{array}{cc}\hfill \alpha & =2·arctan\left(\frac{d}{2·f_o}\right)\hfill \\ \hfill & =2·arctan\left(\frac{1}{f_o}\right)\hfill \end{array}$$

(9)

where d is the width of the normalized image plane. In our case $d=2$ because the range of the normalized image plain goes from -1 to 1.

If now a stereo image of resolution $w\mathrm{x}h$ is rendered, first this resolution will be extended by $\delta$. Then the left and right images with a resolution of $(w+\delta )\mathrm{x}h$ are rendered. Finally, to obtain the off-axis projection plane both images are trimmed respectively to $w\mathrm{x}h$ again.

Figure 5. The correct off-axis stereo camera setup. The extended frustums are depicted. To obtain the off-axis projection plane one has to trim the projection plane of each extended frustum (image source: [26]).

Figure 5. The correct off-axis stereo camera setup. The extended frustums are depicted. To obtain the off-axis projection plane one has to trim the projection plane of each extended frustum (image source: [26]).

2.4. MATLAB ^® Interface

The compiled MATLAB^® command of the volume renderer exposes several input parameters. Thus, we decided to write complementary wrapper classes to further simplify the developer experience.

Figure 9. To provide call by reference instead of call by value Volume and VolumeRender inherit handle. Additionally, the assignment of members does not return a deep copy of the object. Since, instances of LightSource consume low memory it does not inherit handle.

Figure 9. To provide call by reference instead of call by value Volume and VolumeRender inherit handle. Additionally, the assignment of members does not return a deep copy of the object. Since, instances of LightSource consume low memory it does not inherit handle.

2.4.1. Handle Superclass

The realization of the memory management described in Section 2.3.3 requires special techniques on the part of MATLAB^®. Usually a MATLAB^® class member work with call by value. But since we require the possibility to check the pointer of the volumes, call by reference is necessary. More precisely, on the C++ side we want to check whether the volumes are pointing to different pointer addresses or to the same one. Consequently only the unique volume data are copied to the GPU whereas the other assignments are realized by texture mapping to the assigned volume data. For us this seemed to be the most user-friendly way to implement our memory model.

Fortunately, MATLAB^® offers the possibility for call by reference. A class that inherits from the special handle superclass automatically uses call by reference instead of call by value [29]. Hence, the class Volume inherits handle. Through this, all properties stored in a Volume object are pointers.

What is more, a regular value assignment through a setter can be expensive, since the old object is replaced by the new one. This new object is finally returned.

3. Results

3.2. Volume Design

Our renderer supports the capability to define each volume, i.e. emission, absorption and reflection, allowing custom renderings to simulate different materials with different light properties. This is possible to do at a fine-granular scale for each voxel by defining distinct volumes or for an entire volume by a scalar. For our zebrafish volume, Figure 10 shows certain permutations of the scale value. The rendering setup is similar to the scene, but the absorption volume was set to a volume with value 1 for each voxel. The visual impact of the scales and therefore the different volumes becomes clearly discernible.

Figure 10. Rendered images of a zebrafish embryo average brain with different emission, absorption and reflection factors.

Figure 10. Rendered images of a zebrafish embryo average brain with different emission, absorption and reflection factors.

Obviously, without any emission, absorption and reflection information, the image is not visible. On the other hand, while increasing the absorption and adding some reflection, the opacity increases, as the light is reflected and absorbed by the object. High reflection values lead to a shiny appearance what can be counteracted by absorption, while no reflection at all makes the object disappear. Increasing the emission signals increases the effects. Although this illustration demonstrates the effect on the entire volume using the scaling values, it is important to emphasize that it is certainly possible to manipulate an object’s voxels to achieve specific effects in targeted areas. Boolean operations with masking can be used to select areas of interest.

3.3. Case Study

The rendered images demonstrate a high degree of realism, capturing the fine details of zebrafish embryo anatomy. By adjusting the emission, absorption, and reflection parameters, the visualizations effectively differentiate between various tissue types and fluorescent markers. The illumination setup further enhances the depth and contrast of the images, making the internal structures more discernible. Our case study consisted of 240 rendered image containing the zebrafish embryo average brain and some structure. The fish was rotated by 1200 °and one half of the average brain faded out in that image sequence. Figure 11 shows an excerpt of 8 images in which the object rotated 360°and the half of the average brain faded out. The entire video scene can be seen in the Supplementary Materials. The GPU rendering took 119.32 seconds on our hardware setup.

Figure 11. A rendered image of a zebrafish embryo average brain and some brain structures rendered in a separate pass, in pink color. In the sequence shown, the fish rotates 360°while one side of the average brain is faded out and finally rendered transparent, while the structure remains visible in pink. The first image is at the top left and the last at the bottom right.

Figure 11. A rendered image of a zebrafish embryo average brain and some brain structures rendered in a separate pass, in pink color. In the sequence shown, the fish rotates 360°while one side of the average brain is faded out and finally rendered transparent, while the structure remains visible in pink. The first image is at the top left and the last at the bottom right.

Specifically, the first sequence took 4.57s, the second 58.58s and the third one 6.22s. The internal structure was rendered in 49.95s. Obviously, the second scene was the most expensive one, as due to the fading effect, parts of the volume were changed for each frame, requiring a data sync with the GPU. The computational runtime per frame in that sub-scene was 0.33s whereas the other scenes in average costed 0.20s, revealing the costs of host to GPU transfers.

The ability to generate stereo images and animations adds another layer of depth to the visualizations, making it easier to comprehend the spatial relationships between different anatomical features. These realistic visualizations not only aid in detailed anatomical studies but also serve as valuable tools for educational purposes, providing clear and comprehensible images of complex biological structures. Figure 12 shows an image of the scene rendered in stereo using anaglyph. The entire scene in anaglyph can be seen in Supplementary Materials.

Figure 12. Rendered images of a zebrafish embryo average brain and some brain structure in anaglyph. Allures of the part which is fading out are still visible.

Figure 12. Rendered images of a zebrafish embryo average brain and some brain structure in anaglyph. Allures of the part which is fading out are still visible.

4. Discussion

5. Conclusion

References

1. Kaufman, A.E. 43 - Volume Visualization in Medicine. In Handbook of Medical Imaging; Bankman, I.N., Ed.; Biomedical Engineering, Academic Press, 2000; pp. 713–730. [CrossRef]
2. Ronneberger, O.; Liu, K.; Rath, M.; Rueß, D.; Mueller, T.; Skibbe, H.; Drayer, B.; Schmidt, T.; Filippi, A.; Nitschke, R.; et al. ViBE-Z: a framework for 3D virtual colocalization analysis in zebrafish larval brains. Nature Methods 2012, 9, 735–742. [Google Scholar] [CrossRef] [PubMed]
3. Dickie, D.A.; Shenkin, S.D.; Anblagan, D.; Lee, J.; Blesa Cabez, M.; Rodriguez, D.; Boardman, J.P.; Waldman, A.; Job, D.E.; Wardlaw, J.M. Whole Brain Magnetic Resonance Image Atlases: A Systematic Review of Existing Atlases and Caveats for Use in Population Imaging. Frontiers in Neuroinformatics, 11.
4. Kikinis, R.; Pieper, S.D.; Vosburgh, K.G. 3D Slicer: A Platform for Subject-Specific Image Analysis, Visualization, and Clinical Support. In Intraoperative Imaging and Image-Guided Therapy; Jolesz, F.A., Ed.; Springer, 2014; pp. 277–289. [CrossRef]
5. Schindelin, J.; Arganda-Carreras, I.; Frise, E.; Kaynig, V.; Longair, M.; Pietzsch, T.; Preibisch, S.; Rueden, C.; Saalfeld, S.; Schmid, B.; et al. Fiji: an open-source platform for biological-image analysis. Nature Methods, 9, 676–682. Number: 7 Publisher: Nature Publishing Group. [CrossRef]
6. Royer, L.A.; Weigert, M.; Günther, U.; Maghelli, N.; Jug, F.; Sbalzarini, I.F.; Myers, E.W. ClearVolume: open-source live 3D visualization for light-sheet microscopy. Nature Methods, 12, 480–481. Number: 6 Publisher: Nature Publishing Group. [CrossRef]
7. Dhawan, A.P. Medical image analysis; John Wiley & Sons, 2011.
8. Reyes-Aldasoro, C.C. Biomedical image analysis recipes in MATLAB: for life scientists and engineers; John Wiley & Sons, 2015.
9. Mathotaarachchi, S.; Wang, S.; Shin, M.; Pascoal, T.A.; Benedet, A.L.; Kang, M.S.; Beaudry, T.; Fonov, V.S.; Gauthier, S.; Labbe, A.; et al. VoxelStats: A MATLAB Package for Multi-Modal Voxel-Wise Brain Image Analysis. Frontiers in Neuroinformatics, 10.
10. Wait, E.; Winter, M.; Cohen, A.R. Hydra image processor: 5-D GPU image analysis library with MATLAB and python wrappers. Bioinformatics, 35, 5393–5395. [CrossRef]
11. Explore 3-D Volumetric Data with Volume Viewer App - MATLAB & Simulink - MathWorks Deutschland.
12. Kroon, D.J. Volume Render.
13. Robertson, S. Ray Tracing Volume Renderer.
14. Fujii, Y.; Azumi, T.; Nishio, N.; Kato, S.; Edahiro, M. Data Transfer Matters for GPU Computing. In Proceedings of the 2013 International Conference on Parallel and Distributed Systems, 2013, pp. 275–282. ISSN: 1521-9097. [CrossRef]
15. Röttger, S.; Kraus, M.; Ertl, T. Hardware-Accelerated Volume and Isosurface Rendering Based on Cell-Projection, 2000.
16. Wiki, O. Vertex Rendering — OpenGL Wiki,, 2022. [Online; accessed 23-May-2024].
17. Lorensen, W.E.; Cline, H.E. Marching cubes: A high resolution 3D surface construction algorithm. SIGGRAPH Comput. Graph. 1987, 21, 163–169. [Google Scholar] [CrossRef]
18. Levoy, M. Display of Surfaces from Volume Data. IEEE Computer Graphics and Applications 1988, 8, 29–37. [Google Scholar] [CrossRef]
19. Max, N.L. Optical Models for Direct Volume Rendering. IEEE Trans. Vis. Comput. Graph. 1995, 1, 99–108. [Google Scholar] [CrossRef]
20. Wikipedia. Volume ray casting — Wikipedia, The Free Encyclopedia. 2023. Available online: http://en.wikipedia.org/w/index.php?title=Volume%20ray%20casting&oldid=1146671341 (accessed on 18 April 2023).
21. Shannon, C.E. Communication in the presence of noise. In Proceedings of the Proceedings of the Institute of Radio Engineers (IRE), 1949, Vol. 37, pp. 10–21.
22. Henyey, L.; Greenstein, J. Diffuse radiation in the galaxy. The Astrophysical Journal 1941, 93, 70–83. [Google Scholar] [CrossRef]
23. Blinn, J.F. Models of light reflection for computer synthesized pictures. SIGGRAPH Comput. Graph. 1977, 11, 192–198. [Google Scholar] [CrossRef]
24. Porter, T.; Duff, T. Compositing digital images. ACM SIGGRAPH Computer Graphics 1984, 18, 253–259. [Google Scholar] [CrossRef]
25. Ikits, M.; Kniss, J.; Lefohn, A.; Hansen, C. GPU GEMS Chapter 39, Volume Rendering Techniques, 5th ed.; Addison Wesley, 2007; chapter 39.4.3 Rendering.
26. Bourke, P. Calculating Stereo Pairs, 1999.
27. Williams, A.; Barrus, S.; Morley, R.K.; Shirley, P. An Efficient and Robust Ray-Box Intersection Algorithm. journal of graphics, gpu, and game tools 2005, 10, 49–54. [Google Scholar] [CrossRef]
28. NVIDIA Corporation. NVIDIA CUDA C Programming Guide, 2011. Version 4.0.
29. The MathWorks, I. Object-Oriented Programming. Object-Oriented Programming 2011, R2011b. [Google Scholar]
30. Devkota, S.; Pattanaik, S. Deep Learning based Super-Resolution for Medical Volume Visualization with Direct Volume Rendering, 2022, [arXiv:cs.GR/2210.08080].
31. Weiss, J.; Navab, N. Deep Direct Volume Rendering: Learning Visual Feature Mappings From Exemplary Images, 2021, [arXiv:cs.GR/2106.05429].
32. Hu, J.; Yu, C.; Liu, H.; Yan, L.; Wu, Y.; Jin, X. Deep Real-time Volumetric Rendering Using Multi-feature Fusion. In Proceedings of the ACM SIGGRAPH 2023 Conference Proceedings, New York, NY, USA, July 2023; SIGGRAPH ’23, pp. 1–10. [CrossRef]
33. Ikits, M.; Kniss, J.; Lefohn, A.; Hansen, C. GPU GEMS Chapter 39, Volume Rendering Techniques, 5th ed.; Addison Wesley, 2007; chapter 39.5.1 Volumetric Lighting.