Optimization of Breast Tomosynthesis Visualization through 3D Volume Rendering

3D volume rendering may represent a complementary option in the visualization of Digital Breast Tomosynthesis (DBT) examinations by providing an understanding of the underlying data at once. Rendering parameters directly influence the quality of rendered images. The purpose of this work is to study the influence of two of these parameters (voxel dimension in z direction and sampling distance) on DBT rendered data. Both parameters were studied with a real phantom and one clinical DBT data set. The voxel size was changed from 0.085 × 0.085 × 1.0 mm3 to 0.085 × 0.085 × 0.085 mm3 using ten interpolation functions available in the Visualization Toolkit library (VTK) and several sampling distance values were evaluated. The results were investigated at 90º using volume rendering visualization with composite technique. For phantom quantitative analysis, degree of smoothness, contrast-to-noise ratio, and full width at half maximum of a Gaussian curve fitted to the profile of one disk were used. Additionally, the time required for each visualization was also recorded. Hamming interpolation function presented the best compromise in image quality. The sampling distance values that showed a better balance between time and image quality were 0.025 mm and 0.05 mm. With the appropriate rendering parameters, a significant improvement in rendered images was achieved.


Introduction
Breast cancer remains the most common cancer diagnosed among women and a leading cause of death. However, in the last three decades, there has been a decrease of approximately 40% in the death rate from this disease [1,2]. This fact is a direct result of the scientific advances in early detection and treatment. Early detection is mostly done through screenings [3,4]. Until recently, these screenings and breast cancer diagnosis in general were mainly performed by Digital Mammography (DM). However, DM consists of a two-dimensional (2D) acquisition of the three-dimensional (3D) breast causing tissue superposition. This often results in malignant lesions hidden between healthy tissues or normal regions considered as pathological, leading to unnecessary second examinations or biopsies, with additional costs and anxiety for patients [5,6].
Digital Breast Tomosynthesis (DBT) has consolidated its position as a technique to replace DM in both screening and clinical environments [7][8][9][10][11]. DBT has an acquisition geometry very similar to DM but it acquires a set of projection images, allowing a 3D reconstruction of the breast, reducing the tissue overlap observed with DM [12]. In this way, DBT improves the perception of the location and shape of lesions, without increasing the radiation dose to the patient, when compared to DM [13].
The 3D visualization of DBT is one of the most important and crucial aspects to correctly extract the information provided by this technique. Currently, DBT images are displayed through a 2D slice-by-slice visualization [14], with the analysis done one slice at a time or sequentially as a continuous cine loop, leading to a time-consuming process. Two other new approaches have emerged for the visualization of DBT data: synthetic mammography based on DBT data [15,16] and thicker slabs obtained by combining several slices [17,18]. Although synthetic mammography is very useful for comparison with previous DM examinations, it still presents the disadvantages of a 2D visualization (tissue overlapping) [19]. On the other hand, the thicker slabs have revealed good results in reducing time and false positives but have a lower sensitivity [18]. Additionally, a system very recently approved by the Food and Drug Administration, uses artificial intelligence to reduce the number of slices that need to be analyzed by the radiologist. It is based on thicker 6 mm slices combined with synthesized 2D images information [20,21].
A different type of visualization may play an important complementary role in breast cancer diagnosis [22]. 3D volume rendering is the process of creating realistic computer-generated images of volumetric data, yielding a true depth perception [23]. This type of supplementary visualization is already used in other tomographic medical imaging modalities such as Computed Tomography, Magnetic Resonance Imaging, or Positron Emission Tomography [24][25][26]. There are some works mentioning a few aspects of 3D volume rendering for DBT [27][28][29] and its importance to detect clusters of microcalcifications [30].
The classic volume rendering algorithm is ray casting, where rays are cast from the eye or other viewpoints and traverse a scene containing a volumetric data set. This process encompasses several parameters that decide the appearance of the final rendered image. One of these parameters is the sampling distance, which corresponds to the distance between neighboring samples taken along the ray. The value of the sampling distance should be studied and carefully selected accordingly to data set grid resolution [31]. If the distance is too large, our sampling might miss important features in the data and generate major aliasing artifacts. Yet, if we select a very small distance (the number of samples collected along the ray is increased), the amount of time required to render the image will significantly increase [32,33].
On the other hand, reconstructed DBT data typically has voxel sizes of 0.085 × 0.085 × 1.0 mm 3 . The anisotropic nature of the reconstructed DBT data is also responsible for serious quality problems in visualization techniques, namely in the direction orthogonal to the detector plane (z-direction). In this way, one solution is to make the grid isotropic through suitable interpolation functions before the rendering process [31]. The smaller the voxel size, the higher the image definition. However, there are more voxels compounding the data set and therefore the processing time of each data set is longer. The time issue is very important in medical image analysis because a large amount of data needs to be displayed and analyzed in real time. For this study, we consider the total time allocated to the visualization process as the sum of interpolation time with rendering time.
In this paper, two parameters that directly affect the quality of the rendered image-sampling distance and reconstructed voxel size-were considered. The main objective was to improve the quality of the rendered images in the z direction and to determine which options allow a better balance between quality and time. In order to transform data to an isotropic grid, several interpolation functions and their corresponding parameters were tested. Additionally, different sampling distance values were introduced in the rendering process. Qualitative and quantitative analyzes of the results were done through visualization by volume rendering of real DBT images of a phantom. Finally, some values were selected and the rendered images obtained for one clinical DBT data set were also analyzed. To the best of our knowledge, this is the first study about optimization of rendering parameters in visualization of DBT data.

Data Acquisition and Reconstruction
A phantom made by us was used ( Figure 1). This phantom consists of an acrylic background to mimic breast tissue, with two embedded columns of aluminum disks to simulate high density lesions (with different diameters and 1.0 mm thick). For this study, the first column of the phantom (Figure 1b) was considered.
J. Imaging 2020, 6, x FOR PEER REVIEW 3 of 16 A phantom made by us was used ( Figure 1). This phantom consists of an acrylic background to mimic breast tissue, with two embedded columns of aluminum disks to simulate high density lesions (with different diameters and 1.0 mm thick). For this study, the first column of the phantom ( Figure  1b) was considered. Additionally, one clinical DBT data set from an anonymous patient was selected from the clinical facility (Hospital da Luz, Luz Saúde S.A, Lisbon, Portugal) database. Both phantom and clinical data set were acquired with a Siemens MAMMOMAT Inspiration system (Siemens AG, Healthcare Sector, Erlangen, Germany) and reconstructed with the manufacturer algorithm, which uses Filtered Back Projection [34]. The reconstructions have voxel sizes of 0.085 × 0.085 × 1.0 mm 3 .

Data Visualization
In volume rendering, changing the azimuth of a camera rotates its position around the focal point [32]. In this way, it is possible to have an immediate notion of the entire data volume, according to several angles. The methods under study in this work have a particular effect on the z-direction. Therefore, we were particularly interested in the visualization perpendicular to the detector plate (with the projection made on the xz planes, along y), i.e., with the camera at 90º (see Figure 2). In addition, the visualization with the camera at 0°, that is, parallel to the detector plate (projection is made on the xy planes, along z) was also considered.  [32,35] was used to develop 3D specific software in order to visualize DBT data through volume rendering. The methodologies in study were analyzed using 3D volume rendering visualization with composite technique. An Intel ® Core ™ i5-5200U CPU (2.20 GHz) @ 8 GB of memory computer was used.

Image Analysis
For phantom quantitative analysis, the profile of the 5.0 mm disk (Figure 1b in red) was obtained and three figures of merit were used: full width at half maximum (FWHM) of a Gaussian curve fitted Additionally, one clinical DBT data set from an anonymous patient was selected from the clinical facility (Hospital da Luz, Luz Saúde S.A, Lisbon, Portugal) database. Both phantom and clinical data set were acquired with a Siemens MAMMOMAT Inspiration system (Siemens AG, Healthcare Sector, Erlangen, Germany) and reconstructed with the manufacturer algorithm, which uses Filtered Back Projection [34]. The reconstructions have voxel sizes of 0.085 × 0.085 × 1.0 mm 3 .

Data Visualization
In volume rendering, changing the azimuth of a camera rotates its position around the focal point [32]. In this way, it is possible to have an immediate notion of the entire data volume, according to several angles. The methods under study in this work have a particular effect on the z-direction. Therefore, we were particularly interested in the visualization perpendicular to the detector plate (with the projection made on the xz planes, along y), i.e., with the camera at 90º (see Figure 2). In addition, the visualization with the camera at 0 • , that is, parallel to the detector plate (projection is made on the xy planes, along z) was also considered.
J. Imaging 2020, 6, x FOR PEER REVIEW 3 of 16 A phantom made by us was used ( Figure 1). This phantom consists of an acrylic background to mimic breast tissue, with two embedded columns of aluminum disks to simulate high density lesions (with different diameters and 1.0 mm thick). For this study, the first column of the phantom ( Figure  1b) was considered. Additionally, one clinical DBT data set from an anonymous patient was selected from the clinical facility (Hospital da Luz, Luz Saúde S.A, Lisbon, Portugal) database. Both phantom and clinical data set were acquired with a Siemens MAMMOMAT Inspiration system (Siemens AG, Healthcare Sector, Erlangen, Germany) and reconstructed with the manufacturer algorithm, which uses Filtered Back Projection [34]. The reconstructions have voxel sizes of 0.085 × 0.085 × 1.0 mm 3 .

Data Visualization
In volume rendering, changing the azimuth of a camera rotates its position around the focal point [32]. In this way, it is possible to have an immediate notion of the entire data volume, according to several angles. The methods under study in this work have a particular effect on the z-direction. Therefore, we were particularly interested in the visualization perpendicular to the detector plate (with the projection made on the xz planes, along y), i.e., with the camera at 90º (see Figure 2). In addition, the visualization with the camera at 0°, that is, parallel to the detector plate (projection is made on the xy planes, along z) was also considered.  [32,35] was used to develop 3D specific software in order to visualize DBT data through volume rendering. The methodologies in study were analyzed using 3D volume rendering visualization with composite technique. An Intel ® Core ™ i5-5200U CPU (2.20 GHz) @ 8 GB of memory computer was used.

Image Analysis
For phantom quantitative analysis, the profile of the 5.0 mm disk (Figure 1b in red) was obtained and three figures of merit were used: full width at half maximum (FWHM) of a Gaussian curve fitted  [32,35] was used to develop 3D specific software in order to visualize DBT data through volume rendering. The methodologies in study were analyzed using 3D volume rendering visualization with composite technique. An Intel ® Core ™ i5-5200U CPU (2.20 GHz) @ 8 GB of memory computer was used.

Image Analysis
For phantom quantitative analysis, the profile of the 5.0 mm disk (Figure 1b in red) was obtained and three figures of merit were used: full width at half maximum (FWHM) of a Gaussian curve fitted to the 5.0 mm disk's profile, contrast to noise ratio (CNR), and a measure of profile smoothness. FWHM 90º was considered as an indicator of the disk's definition at 90 • . As already mentioned, DBT presents a lower quality in z than in xy and this is exposed in the spreading of structures of high intensity in z (such as these disks). Although the disks have a thickness of 1.0 mm, the observed values will be higher. In this way, a lower value of FWHM 90º will mean a greater definition of the feature at 90 • . For CNR, a region of interest (ROI) over the 5.0 mm disk and other two ROIs over the surrounding background were drawn. CNR was calculated using Equation (1): µ 5.0mm and µ BG stand for mean pixel values in ROI over the 5.0 mm disk and background, respectively; and σ BG stands for the mean of standard deviations in background ROIs. In order to obtain a measure for the profile smoothness (important for the interpolation quality analysis), the STEYX Microsoft Excel ® (Microsoft Office 2013) function was used [36] and its inverse was calculated (Equation (2)).The STEYX function gives a measure of the variability of the data in a given range. The degree of smoothness was calculated considering intensity levels between z-distance (16, 24) mm since it corresponded to an area with high variation in intensity: The quantitative analysis was performed at 90 • , where the methods under study have the greatest effect. For qualitative purposes, displays of the phantom and one clinical data set obtained with volume rendering visualization at 0 • and 90 • are presented.

Study of Interpolation Functions
Due to the acquisition process, DBT data has a finer resolution within slices (in xy planes) and a coarser resolution between slices (in the z-direction). This leads to an anisotropic grid, greatly reducing the quality of rendering. To deal with the anisotropic grid spacing problem, one hypothesis is to change the size of each voxel, considering the smallest dimension, so it can correspond to a perfect cube [37]. Therefore, in this study, the z-resolution has been modified to match the resolution within the slice (which has not changed). This results in a homogeneous resolution, improving data quality after reconstruction and before rendering. At the same time, the number of voxels compounding the data set is also increased and thus more computational memory is required. We can refer to a z-interpolation, since the voxel size went from 0.085 × 0.085 × 1.0 mm 3 to 0.085 × 0.085 × 0.085 mm 3 , that is, the change occurred only in the third dimension.
As this work is based on C++ software developed with the VTK library, the volume data was resampled into an isotropic grid using appropriate interpolation functions available in this library [38]. In VTK, data interpolation is done internally by several classes. The two main classes considered here were vtkImageInterpolator and vtkImageSincInterpolator. The first one is the default interpolator and provides linear, cubic, and nearest-neighbor interpolation. The second is responsible for an approximation to sinc interpolation by multiplying one of the available window functions, in order to limit the kernel size. The window functions studied in this work were: Lanczos, Kaiser (adjusmet parameter off), Cosine, Hann, Hamming, Blackman, and Nuttall [38,39]. In Fourier space, the resolution of spectral window is reduced to the order of the half-width of the sinc function [40]. In this way, the quality of sinc interpolation is related not only to the window function, but also to the window width. For this reason, the window half-width (WHW) was also addressed in this study with values ranging between 1 and 16. vtkImageSincInterpolator class has an option to set blur factors in x, y, and z directions in order to blur the data while interpolating [38]. The size of the kernel is automatically increased by the blur factor (BF), increasing the computation time. As here the interpolation is in the z-direction, some values for BF in z were tested. A summary of the functions and parameters analyzed during the interpolation process is shown in Table 1. For the interpolation study, the default value of sampling distance (1.0 mm) was used.

Study of Sampling Distance
In the ray casting volume rendering algorithm, a ray traverses a volume data set. Along the ray, contributions (samples) based on the intensity values weighted by transparency or opacity (transfer functions) are accumulated at discrete locations of the ray, separated by a certain distance-sampling distance. The process of selecting these locations and, therefore, the distance between them, is subject to the sampling theorem. In this case, it is translated on the condition that the distance between two accumulations must be less than or equal to twice the respective smallest voxel spacing [31].
Taking into account the best results obtained with the interpolation functions and their parameters, several sampling distance values were studied. Considering the sampling theorem, values smaller than 2 × 0.085 mm (0.170 mm) were tested. As previously mentioned, the lower the sampling distance value, the higher the rendering quality should be. As the computation time increases exponentially for lower sampling distance values, our aim was to test several values smaller than 0.170 mm and try to understand where the best balance between quality and time could be achieved. In addition, five other values outside this range were tested. The default value for sampling distance in VTK is 1.0 mm and it was the maximum value included in our study. Included in the class corresponding to volume ray casting, VTK has an option which automatically computes the sampling distance from the data spacing. Using this option, an automatic sampling distance of 0.195 mm was generated. This value was also included in this work (from a brief previous study, it was concluded that, for our DBT data, this automatically generated value was related to the number of voxels according to an approximation of the expression autoSD = 3580.5 × (nr voxels) (−0.621) ). On the other hand, in order to contextualize the results between 0.195 mm and 1.0 mm, some intermediate values (0.40 mm, 0.60 mm, and 0.80 mm) were also considered. In summary, the sampling distance values analyzed were: 0.010 mm, 0.025 mm, 0.050 mm, 0.075 mm, 0.100 mm, 0.145 mm, 0.170 mm (2 × 0.085 mm), 0.195 mm (automatic Sampling distance), 0.4 mm, 0.6 mm, 0.8 mm, and 1.0 mm (default value in VTK).

Results
As mentioned, the time involved to render data is very important in order to make visualization through volume rendering useful. Here, we have separated total time spent in visualization as interpolation time plus render time (from the moment the original data is opened until it reaches the screen). The first depends on voxel size and interpolation functions used in the rescaling. The second is related with ray casting process, namely the sampling distance value. Render time values recorded for rendering the original data, as well as data after interpolation (with sampling distance 1.0 mm) are shown in Figure 3. As the interpolation was done only in the third dimension, all the results presented in this section were measured considering phantom visualization at 90° (i.e., through the z direction).

Linear, Cubic, and Nearest-Neighbor Interpolation
The 5.0 mm disk's profiles obtained after rescaling with the linear, cubic, and nearest-neighbor interpolators are presented in Figure 4. Additionally, it is also shown the profile obtained before rescaling (with voxel size of 0.085 × 0.085 × 1.0 mm 3 ). FWHM and smoothness values measured for each visualization at 90° are presented in Table 2. The total time required in the process is also shown and it stands for the rendered time ( Figure 3) plus interpolation time.

Sinc Interpolation with Different Window Functions
The sinc function was multiplied by different window functions (Lanczos, Kaiser, Cosine, Hann, Hamming, Blackman, and Nuttall). The 5.0 mm disk's profiles obtained with WHW values ranging from 1 to 16 are presented in Figure 5. As the interpolation was done only in the third dimension, all the results presented in this section were measured considering phantom visualization at 90 • (i.e., through the z direction).

Linear, Cubic, and Nearest-Neighbor Interpolation
The 5.0 mm disk's profiles obtained after rescaling with the linear, cubic, and nearest-neighbor interpolators are presented in Figure 4. Additionally, it is also shown the profile obtained before rescaling (with voxel size of 0.085 × 0.085 × 1.0 mm 3 ). FWHM and smoothness values measured for each visualization at 90 • are presented in Table 2. The total time required in the process is also shown and it stands for the rendered time ( As the interpolation was done only in the third dimension, all the results presented in this section were measured considering phantom visualization at 90° (i.e., through the z direction).

Linear, Cubic, and Nearest-Neighbor Interpolation
The 5.0 mm disk's profiles obtained after rescaling with the linear, cubic, and nearest-neighbor interpolators are presented in Figure 4. Additionally, it is also shown the profile obtained before rescaling (with voxel size of 0.085 × 0.085 × 1.0 mm 3 ). FWHM and smoothness values measured for each visualization at 90° are presented in Table 2. The total time required in the process is also shown and it stands for the rendered time ( Figure 3) plus interpolation time.   Figure 5.

Sinc Interpolation with Different Window Functions
The sinc function was multiplied by different window functions (Lanczos, Kaiser, Cosine, Hann, Hamming, Blackman, and Nuttall). The 5.0 mm disk's profiles obtained with WHW values ranging from 1 to 16 are presented in Figure 5. Based on profiles of Figure 5 and for each window function and WHW, the corresponding smoothness and FWHM values were determined and the results are displayed in Figure 6.  Figure 7 from (a) to (f), respectively. As for the WHW, based on profiles of Figure 7, the corresponding smoothness and FWHM values for each BF (z) and window function were determined and the results are presented in Figure 8. Based on profiles of Figure 5 and for each window function and WHW, the corresponding smoothness and FWHM values were determined and the results are displayed in Figure 6.  Figure 7 from (a) to (f), respectively. As for the WHW, based on profiles of Figure 7, the corresponding smoothness and FWHM values for each BF (z) and window function were determined and the results are presented in Figure 8. The profiles calculated with BF (z) values ranging from 1.0 to 4.0 are shown in Figure 7 from (a) to (f), respectively. As for the WHW, based on profiles of Figure 7, the corresponding smoothness and FWHM values for each BF (z) and window function were determined and the results are presented in Figure 8.

Selection of Functions and Parameters for the Sampling Distance Study
The study of interpolation functions and their parameters is, by itself, quite extensive. Thus, to simplify the analysis of the results obtained with different sampling distance values, some selections were made to proceed.

Selection of Functions and Parameters for the Sampling Distance Study
The study of interpolation functions and their parameters is, by itself, quite extensive. Thus, to simplify the analysis of the results obtained with different sampling distance values, some selections were made to proceed.

Selection of Functions and Parameters for the Sampling Distance Study
The study of interpolation functions and their parameters is, by itself, quite extensive. Thus, to simplify the analysis of the results obtained with different sampling distance values, some selections were made to proceed. In summary, the next section results are for cubic interpolation and Hamming window function (sinc) with WHW = 5 or/and BF(z) = 2.

Sampling Distance Study
For each sampling distance, the time required to render the original data (z: 1.0 mm, black in Figure 9) was recorded. In addition, the average rendering times of each selected interpolator (cubic, Hamming with WHW = 5, Hamming with BF (z) = 2 and Hamming with WHW = 5 and BF(z) = 2) were also calculated (gray in Figure 9). J. Imaging 2020, 6, x FOR PEER REVIEW 9 of 16 > 2, these values become very similar. Thus, we choose BF (z) = 2 as the value that represents a better compromise between smoothness, FWHM, and time.

•
Interpolator: Nearest and linear interpolators were excluded since the corresponding profiles showed low smoothness when compared to the others. The cubic interpolator was selected to proceed as it presented smoothness and FWHM values comparable to the other interpolations with a similar interpolation time. For the sinc interpolator, considering the results with WHW = 5 and BF (z) = 2, window functions were sorted according to the profile smoothness (in decreasing order). The Hamming window function presented the best correspondent result between the two options (for WHW = 5: Kaiser > Nuttall > Hamming; for BF (z) = 2: Blackman > Lanczos > Hann > Hamming). Since, among them, window functions presented very close results, this selection of a single function to proceed was done only to simplify and concentrate the next results.
In summary, the next section results are for cubic interpolation and Hamming window function (sinc) with WHW = 5 or/and BF(z) = 2.

Sampling Distance Study
For each sampling distance, the time required to render the original data (z: 1.0 mm, black in Figure 9) was recorded. In addition, the average rendering times of each selected interpolator (cubic, Hamming with WHW = 5, Hamming with BF (z) = 2 and Hamming with WHW = 5 and BF(z) = 2) were also calculated (gray in Figure 9). The 5.0 mm disk's profiles were obtained from the rendered images with each sampling distance value and for each interpolator. The measured profiles between 20 mm and 36 mm (z distance) are shown in Figure 10. The 5.0 mm disk's profiles were obtained from the rendered images with each sampling distance value and for each interpolator. The measured profiles between 20 mm and 36 mm (z distance) are shown in Figure 10.
For a quantitative analysis of the quality of rendered images, smoothness, CNR, and FWHM values were calculated for each sampling distance and the results obtained with the original data and after interpolation are shown in Figure 11. For a quantitative analysis of the quality of rendered images, smoothness, CNR, and FWHM values were calculated for each sampling distance and the results obtained with the original data and after interpolation are shown in Figure 11. For qualitative inspection, sinc interpolation with Hamming window and BF (z) = 2 were selected and sampling distance of 0.025 mm were used. Images of the 5.0 mm disk of the phantom are presented in Figure 12. Images achieved with volume rendering at 0º and 90º are shown in the first and second column of Figure 12, respectively. The first row represents original data with default visualization (without interpolation and with sampling distance 1.0 mm), the second row presents "processed" data with Hamming BF (z) = 2 interpolation, and sampling distance 0.025 mm. The quantitative analysis corresponding to the images in Figure 12 is summarized in Table 3.  For a quantitative analysis of the quality of rendered images, smoothness, CNR, and FWHM values were calculated for each sampling distance and the results obtained with the original data and after interpolation are shown in Figure 11. For qualitative inspection, sinc interpolation with Hamming window and BF (z) = 2 were selected and sampling distance of 0.025 mm were used. Images of the 5.0 mm disk of the phantom are presented in Figure 12. Images achieved with volume rendering at 0º and 90º are shown in the first and second column of Figure 12, respectively. The first row represents original data with default visualization (without interpolation and with sampling distance 1.0 mm), the second row presents "processed" data with Hamming BF (z) = 2 interpolation, and sampling distance 0.025 mm. The quantitative analysis corresponding to the images in Figure 12 is summarized in Table 3. For qualitative inspection, sinc interpolation with Hamming window and BF (z) = 2 were selected and sampling distance of 0.025 mm were used. Images of the 5.0 mm disk of the phantom are presented in Figure 12. Images achieved with volume rendering at 0º and 90º are shown in the first and second column of Figure 12, respectively. The first row represents original data with default visualization (without interpolation and with sampling distance 1.0 mm), the second row presents "processed" data with Hamming BF (z) = 2 interpolation, and sampling distance 0.025 mm. The quantitative analysis corresponding to the images in Figure 12 is summarized in Table 3.

Figure 12.
Volume rendering images at 0º and 90º for the 5-mm disk obtained for the original data with default visualization (top row) and interpolated data with Hamming window with BF (z) = 2 and sampling distance 0.025 mm (bottom row). Table 3. Summary of the results obtained for phantom at 0° and 90°, with default visualization options (voxels with 0.085 × 0.085 × 1.0 mm 3 and sampling distance 1.0 mm) and the options selected in our study (voxels with 0.085 × 0.085 × 0.085 mm 3 after interpolation with Hamming window function and BF (z) = 2 and sampling distance 0.025 mm).

Clinical Data
To evaluate the consistency of the results, sinc interpolator with Hamming window function with BF (z) = 2 was used in the rescaling of one clinical case (0.085 × 0.085 × 1.0 mm 3 → 0.085 × 0.085 × 0.085 mm 3 ). Volume rendering of clinical data was obtained with sampling distance value of 0.025 mm since it showed good results in terms of smoothness, CNR, and FWHM. 2D displays of composite volume rendering of clinical data obtained at 0° and 90° are shown in Figure 13.

Discussion
This type of visualization is an alternative and complementary approach to the standard time consuming slice-by-slice visualization. Here, we optimized the volume rendering visualization for DBT data and our analysis was focused on two main parameters: interpolation methods used before

Clinical Data
To evaluate the consistency of the results, sinc interpolator with Hamming window function with BF (z) = 2 was used in the rescaling of one clinical case (0.085 × 0.085 × 1.0 mm 3 → 0.085 × 0.085 × 0.085 mm 3 ). Volume rendering of clinical data was obtained with sampling distance value of 0.025 mm since it showed good results in terms of smoothness, CNR, and FWHM. 2D displays of composite volume rendering of clinical data obtained at 0 • and 90 • are shown in Figure 13.

Clinical Data
To evaluate the consistency of the results, sinc interpolator with Hamming window function with BF (z) = 2 was used in the rescaling of one clinical case (0.085 × 0.085 × 1.0 mm 3 → 0.085 × 0.085 × 0.085 mm 3 ). Volume rendering of clinical data was obtained with sampling distance value of 0.025 mm since it showed good results in terms of smoothness, CNR, and FWHM. 2D displays of composite volume rendering of clinical data obtained at 0° and 90° are shown in Figure 13.

Discussion
This type of visualization is an alternative and complementary approach to the standard time consuming slice-by-slice visualization. Here, we optimized the volume rendering visualization for DBT data and our analysis was focused on two main parameters: interpolation methods used before  0º and 90º ((a,b), respectively) for original data with default sampling distance (1.0 mm) (a1,b1) and interpolated data with sampling distance 0.025 mm (a2,b2).

Discussion
This type of visualization is an alternative and complementary approach to the standard time consuming slice-by-slice visualization. Here, we optimized the volume rendering visualization for DBT data and our analysis was focused on two main parameters: interpolation methods used before rendering to obtain an isotropic grid (by modifying reconstructed voxel size); and sampling distance values.
For a constant sampling distance in the interpolation study, rendering time proved to be similar between the different interpolation methods (the minimum value recorded was 0.29 seconds for nearest and the maximum was 0.37 s for Nuttall) (Figure 3). On the other hand, when the number of voxels increases from approximately 8 million (0.085 × 0.085 × 1.0 mm 3 ) to 95 million (0.085 × 0.085 × 0.085 mm 3 ), a greater difference was observed. This means that changing the number of voxels had a greater impact on rendering time than the interpolation method itself.
Taking into account the profiles in Figure 4, with a smaller voxel size, a narrower profile was observed, confirming an improvement of image resolution. Among the three interpolators from vtkImageInterpolator class, nearest presented a profile with a blocky appearance resulting from the discontinuous interpolation between neighboring voxels. This translates into a lower smoothness value (greater variability) and a higher FWHM value (Table 3). Linear and cubic interpolators showed similar results for FWHM, with cubic showing a greater smoothness of the profile.
Window functions considered here (Lanczos, Kaiser, Cosine, Hann, Hamming, Blackman, and Nuttall) are already a selection of functions that produce high quality interpolations. Thus, by expecting a similar quality resultant from the different functions, our analysis was focused on two influential parameters: WHW and BF (z). By default, these options are set to 3 and 1.0, respectively, in VTK. Starting with WHW, there was a visible difference in the variability of the profiles in Figure 5, in particular until WHW = 5, with no major differences between interpolators. This observation was corroborated by the calculation of smoothness and FWHM values shown in Figure 6. In Figure 6a, there are six distinct groups, corresponding to the six WHW values tested and, for WHW values above 5, there was no noticeable increase in the smoothness value, increasing only the interpolation time. On the other hand, the variation of this parameter did not cause large fluctuations in the FWHM values for the different interpolators, as can be seen in Figure 6b. About the BF in z, it was found that this parameter has a great influence on the intensity fluctuations existent in the images. For example, from BF (z) = 1.0 (default) to BF (z) = 1.5 (Figure 7a,b, respectively), there was a significant modification in the profile of the disk at 90 • . The results observed in the profiles were confirmed by the numerical analysis of smoothness ( Figure 8a) which increases significantly for higher BF (z) values. It would be expected that the introduction of a blur parameter would increase the smoothness by compromising dispersion in the z-direction (increasing the FWHM value). However, as seen in Figure 8b, the introduction of this factor in the interpolation (up to a certain limit) helps delineating the structures during rendering. Until BF (z) = 2.5, there is a marked decrease in FWHM values for all window functions. Taking into account the temporal information of Figure 8a on the x-axis, we can see that these improvements were achieved with an increase in the interpolation time.
In Section 3.2, the selection considered for the sampling distance study has already been explained (WHW = 5, BF (z) = 2.0 and Hamming window function). While the choice of WHW = 5 was quite simple, for BF (z), we could have opted for 2.0 or 2.5. We selected the first instead of the second value because the interpolation time was shorter. As for the window function in sinc interpolator, it was not an obvious selection, since the different functions presented very similar results. In addition, for example, the Kaiser function has a parameter (α) responsible for the balance between blurring and ringing [31] and so, the effect of α on the results should be considered in a future work. Despite this and in view of the obtained results, Hamming function was the one with the best quality considering both WHW = 5 and BF (z) = 2.
By analyzing the time results obtained for various sampling distance values, if we compare Figures 3 and 9, we found that sampling distance had much more impact on the rendering time than the number of voxels involved. Figure 9 shows an exponential curve with lower sampling distance values corresponding to much longer times. This was expected because sampling distance has a predominant role in the ray casting process.
In Figure 10, it is possible to notice once again the difference in the profiles obtained before and after rescaling (as it had already been observed before in Figures 5 and 7 with sampling distance 1.0 mm). As for the interpolated data, in contrast to the profiles obtained with Hamming, the cubic interpolator showed some oscillation in intensity between 20 mm and 26 mm for all sampling distance values (Figure 10b). On the other hand, it is evident that the smaller the sampling distance, the greater the definition of the plateau, which corresponds to a well-defined disk in the rendered image. In Figure 11, the highlight goes to sampling distance of 0.025 mm, which shows a peak in smoothness and CNR and the lowest value of FWHM. In the same figure, taking into account the different sampling distances, the results for the four interpolations showed to be similar for CNR and FWHM. For the smoothness level of the profiles, Hamming with BF (z) = 2 and Hamming with BF (z) = 2 and WHW = 5 stood out, with the first one presenting the best results in general. All four interpolations lead to a significant improvement in the quality, when compared to the original data (black in Figure 11). This improvement is visible in the rendered images of Figure 12 where, particularly at 90º, the thickness of the disk can be observed with great definition (FWHM of 4.06 mm - Table 3), while with default values what we see is a blurred disk (FWHM of 12.38 mm - Table 3). This translates into a decrease of about 67% in FWHM values. Furthermore, there was an increase in CNR and smoothness of around 500% and 127%, respectively. At 0º, there is also an increase in CNR and a slight decrease in FWHM (about 4%). On the other hand, with interpolation and sampling distance of 0.025 mm, the time required for all the visualization process increases by approximately a factor of 10. Despite remaining at an acceptable value, time continues to be a crucial parameter in rendering and should be optimized.
Another viable option for the sampling distance could be 0.05 mm as it is 1.4 times faster than 0.025 mm, despite suffering some losses in image quality. Here is a summary of the numerical results obtained for the image quality with sampling distance of 0.05 mm (for comparison with the last column of Table 3 A sampling distance of 0.025 mm and rescaling with Hamming BF (z) = 2 were applied to the visualization of a clinical data set. Rendered images of the clinical data ( Figure 13) were in accordance with the results obtained with the phantom in terms of improving image quality, namely at 90 • . There was a visual increase in contrast in the final image ( Figure 13b) and the calcification present in this case had showed better defined contours. The preliminary clinical results are intended to consolidate the results obtained with the phantom. In the future, additional accurate studies with volume rendering of DBT clinical data should be considered to perform quantitative analysis of this type of visualization.

Conclusions
DBT visualization by volume rendering is a new field of research that may support breast cancer diagnosis. With this type of visualization, there is the advantage of observing the entire volume data set at once, from different angles. This is a complex process of visualizing volumetric data set, which includes several factors crucial for the final rendered image. The optimization of these factors is extremely important so that visualization through volume rendering can, along with 2D visualization, have true clinical value. In this way, the main objective of this work was to study and optimize two of these parameters: the interpolation used in the transformation of an anisotropic into an isotropic grid and appropriate sampling distance values, taking into account the entire time required. This is a very extensive study, so the analysis of other interpolation functions and their parameters should be considered in future work.

Acknowledgments:
The authors would like to thank Luis Peralta, who was responsible for the construction of the phantom used.

Conflicts of Interest:
The authors declare no conflict of interest.

2D
Two dimensional 3D Three