DDGS-CT: Direction-Disentangled Gaussian Splatting for Realistic Volume Rendering

Beer-Lambert law:

We would like to kindly correct the reviewer. The exponential transmittance model $T(t)$ used in NeRF/GS, to describe the light attenuation as it travels through a medium from point $r(t_0)$ to $r(t)$, is also based on Beer-Lambert law [17, r1-4]:

$T(t)=\exp(-\int_{t_0}^t a(r(s))\delta s)$,

where $r(t)=o+td$ is a 3D point at step $t$ on a ray of origin $o$ and direction $d$, and $a(p)$ is the medium's absorption/extinction (linked to its density) at point $p$. Assuming a piecewise homogeneous medium (with $a_i$ the absorption of homogeneous segment $l$ of length $\delta_l$), this transmittance can be rewritten:

$T(t)=\exp(-\sum_{l=1}^{l_t} a_l\delta_l)=\prod_{l=1}^{l_t}{(1-\sigma_l)}$,

as in our Eq. 1; with $\sigma_l=1-\exp(-a_l\delta_n)$ and $l_t$ the end segment containing $t$. See [r4] for proof.

X-ray attenuation from photoelectric effect (PE) follows the same model: at each step through the medium (e.g., voxel), the probability of photon extinction (i.e., averaged attenuation) is proportional to the atomic density. E.g., if a ray travels through only 2 voxels, and if it is attenuated by $a_1$=50% at step 1 then $a_2$=50% again at step 2, then its total attenuation is 75% (c.f. $(100-(100-50)^2)$%), in accordance to Beer-Lambert.

Hence, Eq. 1 (3DGS rendering) could be directly applied to DDR synthesis (assuming standard neglog scaling of absorption values into pixel ones [12, 33]), if we were to adopt a simplified model of X-ray imaging which only considers isotropic absorption (i.e., following the Beer-Lambert law), though the term $c_j\in\mathbb{R}^{K}$ (view-dependent radiance) which could be omitted. This is the simplified model proposed in prior work [7], as well as concurrent work [r5] (published online after our submission). In our paper, we propose to go beyond their isotropic simplification and further account for anisotropic scattering of X-rays. This is why we preserve $c_j$ and decompose it into two terms: an isotropic term $c^{\text{iso}}$ and direction-dependent non-linear term $c^{\text{dir}}$ (Eq. 3). Note, we also fix the dimension $K$ of these variables (= number of output channels) to 1 (monochromatic DRR) instead of 3 in usual 3DGS (RGB). We will clarify this in the paper.

Contributions w.r.t. 3DGS:

We would argue that while Eq. 3 is compact, it is not straightforward; nor it covers all our contributions (e.g., radiodensity-aware initialization). It may be easily overseen as it is formalized via subscripts ($i$ vs. $j$) but a key contribution is the decoupling of the isotropic and view-dependent terms, modeled by distinct Gaussians ($g^{\text{iso}}_i$ and $g^{\text{dir}}_j$) c.f. L163-164 and 171-177. This differs from prior 3DGS and DRR work [7, 24], using 1 set of Gaussians. Our dual-set idea may seem unintuitive but demonstrates better compactness and accuracy. We believe this insight benefits the community.

In comparison, X-Gaussian [7] (ECCV 2024) adapts 3DGS to DRR rendering by: (a) simplifying view-dependent radiance into an isotropic term (similar to our $c^{\text{iso}}$); and (b) replacing SfM initialization with uniform 3D space sampling (ignoring CT content). Our contributions are comparatively significant: (a) a 2nd anisotropic term $c^{\text{dir}}$ with its own Gaussians $g^{\text{dir}}$ (a more accurate physics adaptation of 3DGS to X-rays); and (b) an advanced initialization strategy considering CT data distribution; overall resulting in a +1.04db PSNR increase and a 10% decrease in needed Gaussians compared to [7] over CTPelvic1K (Tab. 1).

MC/real comparison:

Please refer to our global response. We cannot compare MC and analytical solutions against real images but:

• Tab. 1 already provides a relevant comparison between DRR tools (as candidates) and MC simulation (TIGRE [5], as GT). The high metrics for DDGS show that it approximates MC simulation well.
• Similarly, the registration experiments (e.g., Tab. 7) implicitly provide insight into the similarity between real and synthetic data, as pose optimization is done by comparing DRRs to real target scans.
• For a more direct comparison, we used our DDGS-based registration method to refine GT poses in Ljubljana data (c.f. aforementioned GT inaccuracies), then used the refined poses to render and compare DRRs to real scans. Results can be found in Tab. R1 and Fig. R1 (rebuttal PDF). E.g., when zooming, we notice that DiffDRR suffers from pixelization, unlike GS methods. Refined GT data will be made public.

High-dim. residual contribution:

Thank you, we will clarify. X-ray scattering is a complex phenomenon causing imaging noise (residual impact). The noise distribution is influenced by so many parameters (source direction and energy, medium properties, bounces, etc.), only a costly high-dimensional function could accurately model it (see Sec. 2.1 of [8] for details). Therefore, we only claim to approximate the residual impact of scattering.

CT noise:

We agree that the quality of DRRs is bound to the quality of the input CT scans. This issue affects all DRR and MC methods. E.g., MC-based tools segment CT scans into materials and use the segmented data as simulation input. If a 3D scan is noisy, so will the material segmentation and simulation results (see discussion in [33]). When our method uses simulators to get target images, it will face similar noise issues. However, our experiments show that DDGS barely introduces additional noise (highest PSNR in Tabs. 1-2). Compensating for CT noise is an interesting, under-explored research direction that could benefit all DRR methods; but, we believe, beyond the scope of our study.

Clinical usage:

When clinicians use DRRs for quick visualization, they know the tools' limitations. E.g., DRRs are typically not used for diagnosis. We agree with the reviewer that image quality is still important, hence our experiments showing DDGS accuracy over similar DRR works.

We appreciate g5Vp's feedback and are glad to hear that "the author's answer solved some of my doubts." Given that some of the previously mentioned weaknesses have been addressed, may we respectfully ask g5Vp to kindly elaborate on their decision to nevertheless opt for borderline reject? We are eager to address any remaining concerns and to further improve our submission based on the feedback.

In our view, we have responded to the key concerns raised by g5Vp, specifically:

• Weakness #1 = theoretical misunderstanding by g5Vp. We hope that we have clarified the application of the Beer-Lambert Law to correct the reviewer's understanding.

• Weakness #2 = significant contributions compared to SOTA. We acknowledge that the novelty and impact of a paper are subjective, so we have highlighted the theoretical and empirical significance of our contributions relative to the SOTA [7, 11, 12].

• Weakness #3 = additional results provided. We exceeded the standard analytical DRR evaluation [5, 7, 11, 12, 33] by comparing our solution to real projections, as requested by the reviewer.

We have also thoroughly addressed the questions of other reviewers, particularly those raised by ZGY5 regarding the SOTA status of [7, 11, 12] and the paper’s organization. We note that overall, the reviewers' concerns/questions do not have major overlaps, as opposed to their positive feedback ("novelty" [Zg8U, YFhr, ZGY5], "substantial improvement on target tasks" [all], "clear motivation/methodology" [all], etc.).

Finally, we would like to emphasize that we submitted this paper to NeurIPS's "Machine learning for healthcare" primary area. We believe our solution is well-suited to this application area, given its ML contributions and focus on clinical tasks. Thus, we respectfully disagree with the reviewer's suggestion that our submission might be better suited to a different venue.

We kindly ask g5Vp to consider the arguments above; and if g5Vp still believes that our submission is not appropriate for NeurIPS (Machine learning for healthcare), we would be thankful for an explanation.

Sincerely,

The Authors

Effectiveness of the radiodensity-aware dual sampling:

We believe that the accuracy increase compared to the SOTA [7] (i.e., performing uniform sampling) brought by this single contribution (novel radiodensity-based sampling) is still significant.

We should further highlight that our novel CT-specific initialization scheme significantly improves the convergence of 3DGS models for DRR rendering, as shown in the new Table R3 (see rebuttal PDF). E.g., after 500 iterations, our model with the proposed sampling strategy outperforms the one with SOTA even-sampling [7] by +2.06dB, and the one with random-sampling [24] by +1.31dB. The gap slowly reduces as models converge and reach saturation in terms of image quality; but this shows that our contribution facilitates convergence. We will add this insight to the manuscript.

CT volume visualization:

We thank the reviewer for their insightful suggestion. We will add some volumetric visualizations of input CT scans to our qualitative figures. Please refer to the rebuttal Figure R1 (in the attached PDF) to see how they would look. Note that we used the standard medical imaging platform Slicer [r6] to generate these CT visualizations.

Error map visualization:

We thank the reviewer for another valid suggestion. We will add signed error maps to Figure 3 (we find signed error maps to be more insightful than error maps for DRRs, as they highlight where the methods may underestimate or overestimate the overall photon absorption). An example of how these error maps would look can be found in the new Figure R1 (in the attached PDF), where we compared DRR outputs to real X-ray scans (c.f. discussion with Reviewer g5Vp).

Typos:

Thank you, we will further proofread our manuscript before final submission.

For the mentioned typo example, however, we would like to highlight that $]a, b]$ and $(a, b]$ are both valid, equivalent, notations for the left-open interval $a < x \leq b$ (c.f. standard ISO 31-11 for mathematical notations). We understand that the second notation may be more prevalent in our domain and will change the paper accordingly.

SOTA status of X-Gaussian and DiffDRR/DiffPose:

We respectfully disagree with the reviewer w.r.t. where the SOTA stands. Since we started writing this paper, DiffPose [12] (an extension of DiffDRR [11] by the same authors) has been presented at CVPR 2024 (oral), and X-Gaussian [7] has been approved at ECCV 2024. We believe that this demonstrates that our analysis of the state of the field is correct and that our comparison to these methods is essential. Note that, at a fundamental level, we do agree that the current trend in machine learning to consider preprints as legitimate work can be counterproductive to science's steady progress. In this case, however, we made a decision, now validated by peers, based on the quality of these works [7, 11, 12] and our knowledge of the field.

The application of novel inverse graphics techniques (e.g., 3DGS) to DRR rendering is a fast-moving field, but we did our best to reference the relevant methods. We especially consider X-Gaussian [7] and DiffDRR/DiffPose [11, 12] to be the most important for comparison, as they each represents the state-of-the-art in their respective type of analytical DRR rendering (X-Gaussian for 3DGS-based DRR rendering, and DiffPose as a differentiable extension of Siddon's voxel-based algorithm [30] enabling optimization-based tasks), directly relevant to our work.

Manuscript organization:

We thank the reviewer for finding our motivation and contributions clear, and for the suggestion to shuffle a bit the Related Work and Preliminaries sections. We agree that our manuscript may benefit from moving some application-related sentences to the former section.

Regarding the suggestion to express our intentions as clearly as our motivation: our key goal is to define a DRR rendering solution that (a) is fast and differentiable in order to enable intraoperative tasks (e.g., CT/X-ray registration); but also that (b) matches more closely the complex X-ray imaging process without sacrificing performance. Our intuition is that a better trade-off between realism and computational performance (compared to recent 3DGS-based solutions for DRRs, which ignore anisotropic properties) could further help downstream tasks, as demonstrated in the paper (accurate MC simulation approximation and pose estimation).

Additional ablation results:

We understand that performing some of the ablation studies on data subsets is not ideal. We believed that the selected instances were statistically representative enough and wanted to preserve space and computing power for the other quantitative results. It is not ideal, but also common with rich 3D datasets to perform ablation studies on a subset [8, 17, 31, 35, 36, etc.].

We complete our ablation studies, with results shared in Tables R2-R3 (in the rebuttal PDF). More will be added to the paper and supplementary material.

Symbol explanations:

Could the reviewer specify which unexplained symbols they are referring to, please? Note that we went through our manuscript and found 1 typo, L208: we meant to write "loss weight $\lambda$" rather than "loss weight $\delta$" (referring to the symbol in Equation 2). We thank the reviewer for helping us correct this.

Regarding $g_i^{\text{iso}}$ and $c_i^{\text{iso}}$ (as well as $g_j^{\text{dir}}$ and $c_j^{\text{dir}}$), those are different concepts (L163-165). We used the notation $g$ to represent Gaussians (each defined by its 3D mean position, covariance, opacity, and color/intensity features c.f. L113). So $g_i^{\text{iso}}$ and $g_j^{\text{dir}}$ each represents a distinct set of Gaussians, with different 3D positions, covariances, opacities, and, more importantly, different feature vectors (resp. isotropic $\mathbf{f}_i^{\text{iso}}$ and direction-dependent $\mathbf{f}_j^{\text{dir}}$).

In contrast, we used the symbol $c$ in the paper to represent the color/absorption contribution of each Gaussian to the image (c.f. Equation 1). Therefore, $c_i^{\text{iso}}$ represents the contribution of the Gaussian $g_i^{\text{iso}}$, calculated according to Equation 3 (i.e., multiplying its feature vector $\mathbf{f}_i^{\text{iso}}$ to the global isotropic basis vector and then applying a sigmoid activation), and similarly for $c_j^{\text{dir}}$ being the contribution of $c_j^{\text{dir}}$ to the pixel, after applying our proposed anisotropic model (c.f. Equation 3).

We would gladly take into account any further suggestions to make our manuscript accessible to a larger audience. E.g., we will improve our pipeline figure (Fig. 1) by adding the above-mentioned notations within, in order to help the readers understand their application.

Isotropic vs. anisotropic sets — motivation, initialization, and rendering:

We hope that the following points will clarify some misunderstandings:

• We define 2 different functions for the absorption contribution of the 2 Gaussian sets: one function is isotropic to approximate average radio-absorption ($c_i^{\text{iso}}$), and one function is direction-dependent to approximate the contribution of anisotropic Compton scattering to the image ($c_j^{\text{dir}}$).
• Since we have 2 distinct Gaussian sets, they could indeed be rasterized separately, as done in Figure 4 purely for qualitative ablation. However, the two sets need to be rendered together to ensure correct simulation. Otherwise, Gaussians occluded by others belonging to the other set could incorrectly contribute to the ray absorption. Please see Figure R2 in rebuttal PDF for an example: there, point $g_2^{\text{dir}}$ should not contribute to the final pixel value, as it is occluded by $g_2^{\text{iso}}$ (i.e., it absorbs all the remaining energy before the ray can reach $g_2^{\text{dir}}$). This could not be properly simulated if each Gaussian set is rasterized separately.
• By having 2 sets of Gaussians with distinct yet complementary contribution functions (isotropic vs. direction-dependent functions), our solution can better model DRR imaging. During optimization, each Gaussian set will approximate its respective imaging effect, conditioned by its respective function.
• We randomly split $\mathbf{P}^{\text{mc}}$ (3D points sampled based on the CT scan's distribution) into 2 sets, each subset contributing to the initialization positions of $g_i^{\text{iso}}$ and $g_j^{\text{dir}}$ (half the 3D points in $\mathbf{P}^{\text{mc}}$ will serve as initialization for $g_i^{\text{iso}}$, along with the entire $\mathbf{P}^{\text{mc}}$ set; and the other half of $\mathbf{P}^{\text{mc}}$ will serve as initialization for $g_j^{\text{dir}}$). After initialization, during 3DGS-based model optimization, the points in each set of Gaussians can move/split/drop independently (along with the optimization of their respective covariance/opacity/feature values). Since the optimization of each set (isotropic and anisotropic) is conditioned on their respective differentiable absorption function, they will acquire distinct properties.

We will add more details to the paper and will provide the new figure in the annex of the camera-ready version.

Comparison to "DRR-based methods":

Since our method belongs to the family of analytical DRR rendering methods [7, 11, 29, 30, 24] and to the more recent sub-family of 3DGS-based solutions [7, 24], we evaluated it accordingly. As standard in this domain [5, 7, 11, 12, 33] (see discussion with Reviewer g5Vp), we compared on meaningful downstream tasks (CT/X-ray pose estimation and 3D/2D keypoint registration) with corresponding SOTA methods [11, 12, 7] and baselines [17] (note that we did not consider [24], as its contributions w.r.t. [17] are marginal). Additionally, though less common in the domain of analytical DRR rendering (e.g., [7, 11] only compare to a few methods), we also provided an image-quality evaluation over several datasets, comparing to Monte-Carlo (MC) simulation results (used as GT, which is standard in DRR evaluation [5, 7, 11, 33]). Based on a suggestion by Reviewer g5Vp, we further provide some quantitative comparison to real projections (see rebuttal PDF).

Therefore, by combining image-quality evaluation (comparing to MC and real data) and downstream-task evaluation on multiple datasets, against several relevant SOTA analytical methods, we believe that we have been more thorough than most DRR papers.

Voxel vs. Gaussian representations:

The reviewer is correct that both the original CT scan representation (voxel grid) and the proposed 3DGS-based one are explicit. However, as shown in our experiments, our 3DGS representation is much more compact. E.g., for the CTPelvic1K data, our model can represent a CT scan with only ~42,625 Gaussians (c.f. Table 1), defined by 19 float values each (3D position/rotation + 3D covariance + opacity + feature vector) and a shared 32-dimensional basis vector $\mathbf{b}$, so $809,907$ float values in total; whereas the original CT scans are each composed of $512\times 512\times 500=93,363,200$ values. E.g., the voxel data can contain large homogeneous regions which can be approximated with few Gaussians, resulting in a significant compression rate.

Moreover, the original voxel grids contain raw absorption values, which need further heavy processing to be projected into DRRs (volume traversal and rendering). In comparison, the 3DGS representation is already meant for fast rendering, i.e., facilitating the rasterization process. For the sake of DRR generation, we demonstrate (e.g., in Table 7, against DiffDRR [11, 12] which relies on voxel representation) that 3DGS models are much lighter, both in terms of forward rendering and backward propagation (less operations), making them better suited for integration into larger automation systems.