Among various three-dimension (3D) display techniques, computer-generated holography is an ideal and more flexible technology^[1–4]. Computer-generated holograms (CGHs) make it possible to record virtual objects and store holograms digitally. However, a huge computing load and the expression of diverse types of virtual objects are problems that need to be addressed in CGH^[5–10]. Recently, the improvement of programmability of the graphics processing unit (GPU) decreases the masses of computation in CGH^[11–14], but it is not easy to implement different types of virtual objects like traditional computer graphics (CG) techniques by conventional methods for calculating CGH. Besides, some researches consider the physical phenomena effects (e.g.,  occlusions, refraction, and reflection) in holographic display^[15,16], but it is still inadequate for multiple requirements of reconstructing virtual objects. Traditional computer-generated holography used surface models, lacking the ability of expressing the interior of an object. In contrast, volume data contains the raw information of the 3D data field, including the interior detail^[17–19]. In recent researches, several methods have been proposed to generate CGH with volume data. The maximum intensity projection (MIP)-based method only focuses on the local maximum intensity between the volume data, thus limiting the application of scenarios. In the polygon-based method, volume data is transformed into polygonal models before CGH. However, the process is time-consuming, and a lot of raw information is missing^[20]. The approximate volume rendering method regards each voxel as a point light source, so it can generate holograms for transparent objects and multi-layer objects, but the computation usually requires hours^[21].

In this work, we propose a CGH generation method for volume data to flexibly display abundant information of the volumetric field without additional polygonal model transformation. The proposed method realizes optional holographic visualization of the interior and surface of objects using the adjustable ray to optical-wave conversion method. Ray-casting rendering and the optical-wave superimposed method are used in the holographic recording process. Moreover, our method is implemented on GPU in parallel to reduce the CGH generation time. The specific scheme is shown as follows.

There are two stages in generating CGH of volume data: rays are cast from the hologram plane to collect and process efficient voxels in 3D virtual space, and then composited data acquired from the rays are directly converted into back propagating optical waves to be superimposed on the holographic plane. The adjustable ray to optical-wave conversion method is illustrated in Fig. 1. A bundle of rays is emitted from the center of the hologram element after determining the direction of each ray. Equal interval path sampling is carried out if the rays intersect with the volume data. Suppose the origin of the ray is $\mathit{O}$, the normalized direction vector of the ray is $\mathit{d}$, and the ray can be expressed by parametric function $\mathit{r}(l)=\mathit{O}+l\mathit{d}$, $l\ge 0$, where $l$ represents the distance between the sampling point and the origin of the ray. If the ray enters the volume data at $l_{\mathrm{begin}}$ and exits from the volume data at $l_{\mathrm{end}}$, then the domain of $l$ is determined as $[l_{\mathrm{begin}},l_{\mathrm{end}}]$. The coordinate of $l_{\mathrm{begin}}$ can be converted by ${\mathit{r}}_{\mathrm{begin}}=\mathit{O}+{\mathit{l}}_{\mathrm{begin}}\mathit{d}$, and the next sampling position is $\mathit{r}={\mathit{r}}_{\mathrm{begin}}+\mathrm{\Delta }l\mathit{d}$, where $\mathrm{\Delta }l$ is the sampling interval. The sampling value can be derived by trilinear interpolation of eight voxel values adjacent to the position, as shown in Fig. 2(a).

The raw voxel value is a scalar value. Considering that objects of the same type have some similarity in scalar values, the scalar value of volume data can be used as a data feature for the classification basis of different regions of objects. When a scalar sample value is obtained, it is mapped to a group of four-dimensional values with optical properties (RGB color values and opacity) by transfer function, as shown in Fig. 2(b). The transfer function can be regarded as the control function of the object light wave. By adjusting the control function, the presentation or suppression of the information of the sample can be controlled to affect the overall complex amplitude of the virtual object. The mapping process is expressed as $$I_n(C_n,O_n)=I_n(R_n,G_n,B_n,O_n)=\mathrm{TF}[v(s_n)],$$(1)where $v(s_n)$ is the scalar value of the $n$th voxel $s_n$, $\mathrm{TF}(v)$ is the transfer function, and $I_n$ is the mapped value containing color component $C_n$ and the opacity $O_n$. The values of the RGB components are, respectively, represented by $R_n$, $G_n$, and $B_n$. A transfer function can be generated through some specific discrete points in the function. Once the mapping relation of specific points is determined, the transfer function can be determined by fitting and interpolation. By adjusting the number and mapping relation of these control points, $\mathrm{TF}(v)$ can be flexibly changed.

After the mapping process, information composition is analyzed with the absorption-emission model^[22]. Each voxel can absorb incident light and emit light, but there is no scattering between the voxels. The discretization representation of this process can be described as Eqs. (2) and (3): $$C_{{\text{sum}}_n}=C_{{\text{sum}}_{n-1}}+(1-O_{{\text{sum}}_{n-1}})C_n,$$(2)$$O_{{\text{sum}}_n}=O_{{\text{sum}}_{n-1}}+(1-O_{{\text{sum}}_{n-1}})O_n,$$(3)where $C_{{\text{sum}}_n}$ and $O_{{\text{sum}}_n}$ represent the cumulative values of color and opacity when the ray traverses to the $n$th voxel. $C_{{\text{sum}}_{n-1}}$ and $O_{{\text{sum}}_{n-1}}$ are the color and opacity at the last sampling position. Then, the final composite color is converted to the gray value as the amplitude of the virtual synthetic object point of this ray, and the location of the first sampling point is taken with the cumulative opacity greater than zero as the spatial coordinate of the virtual synthetic object point.

The synthetic object point propagates light waves to the hologram plane in the form of spherical waves. Therefore, the complex amplitude distributions $E(x,y)$ in the hologram plane are parallelly calculated by $$E(x,y)=\sum _{j=1}^M{A}_j\exp (i{\phi }_j),$$(4)$${\phi }_j=\frac{2\pi }{\lambda }{r}_j,$$(5)$$r_j=\sqrt{{(x_j-x)}^2+{(y_j-y)}^2+z_j^2},$$(6)where $x$ and $y$ represent the horizontal and vertical coordinates of every hologram sample, ${\phi }_j$ is the phase information, $\lambda$ represents the wavelength, $M$ is the number of rays that hit the object, $A_j$ is the amplitude of the synthetic object point of the specific ray $j$, and $r_j$ is the distance between the virtual synthetic object point $(x_j,y_j,z_j)$ and the hologram sample $(x,y,0)$.

Figure 3 summarizes the whole process described above, and Fig. 4 is the light path setup. There is a beam splitter for changing the light path, a polarizer for converting the polarization state of laser beam, and a $4f$ system with a filter for reducing low-frequency noise interference. A central processing unit (CPU) launches the GPU kernel after configuring parameters of the transfer function. GPU parallelly computes the value of hologram units using the information of each ray. If the scale and offset of the transfer function are changed, or there are interactive operations, the ray-casting process is restarted for re-calculating CGH. Lastly, the CGH data is loaded on the spatial light modulator (SLM) for holographic display.

The volume data is bound to texture memory of the GPU to achieve random and fast reading. We use CUDA-9.1, the parallel computing architecture of NVIDIA, to accelerate the CGH calculation process. The CPU and GPU are Intel Core i7-7700HQ and Nvidia GTX 1060. Intersection detection and early ray termination strategy are adopted to accelerate the ray-casting process. Intersection detection terminates the missing ray that does not hit the object in advance. Early ray termination strategy terminates the ray when the accumulated opacity reaches 1.0. By reducing meaningless operation, these strategies have a positive effect on acceleration. The time trend of the proposed method is shown in Fig. 5 (the size of input volume data: $256\times 256\times 256$). In the traditional polygon-based method in Ref. [20], it takes 23.9 s to generate CGHs for volume data with a GPU. In the traditional approximate volume rendering method in Ref. [21], it takes up to 3.6 h to generate CGHs for volume data. Our method only takes 0.34 s (the number of rays: $100\times 100$, the size of hologram: $1920\times 1080$). The results of Fig. 5 show good time performance in CGH generation for volume data. Besides, the time will decrease when the number of rays or the size of CGH reduces, which indicates that our work can provide dynamic interactive holographic display.

The optical experiment in Fig. 6 verifies the 3D properties of reconstructions. It uses the industrial volume data [computed tomography (CT) scan of two cylinders of an engine block, size: $256\times 256\times 128$] as the data source. The hologram is placed at the origin of the $z$ axis, the $z$ coordinate of volume data was set between 34 cm and 42 cm, and the $x$ and $y$ coordinates of volume data are both set between $-1\mathrm{cm}$ and 1 cm. Changes in the clear and blur parts of the object can be perceived when the camera is focused on different positions. It illustrates that the proposed method can represent the continuous depth information effectively. Figure 7 presents reconstructed images at different view directions. In the conventional holographic display method of volume data based on MIP in Ref. [20], the reconstructed object is symmetrical when viewing from opposite directions. The information loss is due to the limitations of the method itself, since it only takes maximum values in the path of the rays. The proposed method shows more information of the object. Observers can perceive different reconstructed images in opposite directions.

To demonstrate the adjustable reconstructions of the proposed method, we used the medical volume data (rotational $c$-arm X-ray scan of a human skull and a human foot, size: $256\times 256\times 256$) as the data source. By adjusting the transfer function, the reconstruction can display different content to meet different viewing requirements. Figure 8(a) shows the different regions of interest in the skull volume data under different transfer functions, and Fig. 8(b) shows different concerns on tissue and bone information of the foot. Both the surface and the interior information of an object are expressed. Especially, the relationship between different organizational structures is easier to be understood, and diverse information in objects can be optically reconstructed.

In summary, a CGH generation method for expressing adjustable reconstructed content of volume data is proposed. It composes optical properties along the viewing direction and directly accumulates the synthetic data into the propagation process of light waves. Different from traditional monotonous holographic reconstructions of the objects’ surfaces, the method reconstructs the interior of objects optionally. Besides, with the adjustable ray to optical-wave conversion, the proposed method can flexibly and directly calculate the complex amplitude of a 3D virtual scene instead of the conventional input of constant 2D color and depth image. The proposed method offers a whole procedure of adjustable hologram of volume data, which can dynamically reconstruct different regions of objects according to different observation requirements. It is expected to be combined with several fields such as medical and industrial fields to achieve more diverse applications in the future.