GBIR: A Novel Gaussian Iterative Method for Medical Image Reconstruction

Thank you for your time and effort in reviewing our paper. Below are our responses to your comments:

Q1: The method introduces additional computational complexity compared to traditional deep learning methods.
Answer: As indicated in Table 3 of our initial submission, we have studied the computational cost of different methods. Notably, the proposed GBIR method is more computationally efficient than advanced diffusion models while maintaining high accuracy. We believe this demonstrates the practical advantages of GBIR in terms of computational efficiency.

Q2: The paper briefly addresses diffusion models.
Answer: The focus of this paper is not on diffusion models; rather, they serve as baseline methods for comparison. In contrast, we introduce a novel Gaussian-based iterative reconstruction method (GBIR). We have also discussed and analyzed diffusion models in detail in Section 5.4, titled "Benchmark and Findings on MORE Dataset," to provide a comprehensive comparison of their performance against our method.

Q3: How does the performance of GBIR compare on a prospective dataset?
Answer: We have thoroughly analyzed the performance of GBIR on the MORE dataset, as detailed in over 20 tables within the manuscript. Additionally, we benchmarked various methods on this dataset to provide a comprehensive evaluation. This extensive analysis illustrates the effectiveness and robustness of GBIR.

Q4: The paper lacks a figure illustrating the framework.
Answer: Figure 2 is the manuscript provides an illustration of the proposed GBIR framework. We hope this visual clarification helps in understanding the structure and operation of the method.

We would like to express our sincere gratitude for your valuable feedback. Below are our responses to your comments:

Q1: The authors are unaware of physics-driven methods that are considered the state-of-the-art in MRI reconstruction.
Answer: Thank you for your suggestion. We will revise the manuscript to include a more comprehensive discussion of related work, including physics-driven methods, to provide a broader context for our work.

Q2: The authors treat their Gaussians as isotropic, yet each Gaussian is specified to be non-isotropic with a covariance function $\Sigma_k$ in (4).
Answer: Sorry for the confusion. Our Gaussian representation is indeed isotropic, and the covariance matrix $\Sigma_k$ is diagonal. We will clarify this point in the revised manuscript.

Q3: I cannot find where the authors discuss how they learn $\Sigma$, though I can see they talk about learning $\mu$.
Answer: The covariance matrix $\Sigma_k$ is involved in the computation of the Mahalanobis distance tensor $D^2_{n,c,h,w}$ in Equation (8), which plays a role in the reconstruction process. By minimizing the difference between predicted and ground-truth measurements, the network iteratively learns the covariance matrix $\Sigma_k$. We will clarify this in the revised manuscript.

Q4: Eq. (11) for MRI suggests the authors are unaware that virtually all MRI raw data is multi-coil.
Answer: The MRI measurements used in our experiments are simulated, following the practice in our baseline method DiffusionMBIR [1]. This is different from real MRI data. The design choice is motivated by two factors: (1) to ensure a fair comparison with the baseline method, and (2) to facilitate the general framework of our proposed approach.

Q5: The authors test their MRI methods on the BRATS dataset. Unfortunately, this dataset contains no raw data.
Answer: You are correct that the BRATS dataset does not contain raw k-space data. All methods, including ours, use simulated k-space data under identical settings, following the baseline DiffusionMBIR [1], to ensure a fair comparison. We will clarify this in the revised manuscript.

Q6: MRI reconstruction and CT reconstruction do not even have the same output size.
Answer: We apologize for the confusion. In our MRI reconstruction experiments, we only reconstruct the magnitude images, which are real-valued. We will clarify this in the revised manuscript to avoid further confusion.

Q7: The MORE dataset is interesting. On a minor note, details are missing on demographics. On a more major note, it is unclear if there are any matched MRI/CT data on the same subjects or if this is a mere concatenation of an MRI dataset and a CT dataset. Finally, does the MRI data actually contain raw data or just DICOMs?
Answer: Thank you for your interest in the MORE dataset. The MRI data in the current MORE dataset only contains DICOMs, and all experiments are conducted on simulated k-space data. The MRI and CT data are not matched on the same subjects; instead, the dataset is a concatenation of separate MRI and CT datasets, aimed at providing a comprehensive benchmark for medical image reconstruction. We will add more details about demographics and the dataset in the revised manuscript.

Q8: Long reconstruction time makes this method clinically unusable.
Answer: This is a valid limitation, as acknowledged in our manuscript. However, our method is faster than many advanced diffusion models, as shown in Tables 1 and 3. We recognize that further improvements in efficiency are needed, and we plan to address this in future work.

Q9: The MRI comparisons are outdated. Nobody would use AUTOMAP for MRI reconstruction.
Answer: We acknowledge that AUTOMAP may be outdated for MRI reconstruction. However, our large benchmark on the MORE dataset includes both traditional CNN methods and advanced diffusion models, providing a comprehensive evaluation of various methods, including our proposed approach. Our goal is not only to demonstrate the superiority of our method but also to offer a comprehensive benchmark for the community. Furthermore, the comparison with Score-MRI [2] and DiffusionMBIR [1] illustrates the effectiveness of our approach.

Q10: Relevant methods to the proposed method.
Answer: Thank you for your recommendation. We have read the papers you mentioned and will cite them in the revised manuscript. Our method is distinct from these approaches in that we propose a Gaussian-based reconstruction framework for CT and MRI. This approach is similar (but distinct) to 3D Gaussian Splatting (3DGS), a 3D reconstruction method. We will provide further clarification in the revised manuscript.

Thank you sincerely for your time and valuable feedback. Below is our response to the reviewers' comments:

Q1. The paper lacks qualitative comparisons.
Answer: We provided the visualization in Appendix Figure 5. We will include more qualitative comparisons in the revised manuscript to provide a more comprehensive evaluation of our proposed method.

Q2: The proposed method may lack novelty.
Answer: We respectfully disagree with this comment. Our proposed GBIR is not based on 3D Gaussian Splatting, which is different from the X-Gaussian [1], $R^2$ Gaussian [2], and 3DGR-CAR [3]. Instead of adapting 3DGS for medical image reconstruction, we propose a different 3D Gaussian representation and as well as an efficient iterative reconstruction framework. We will further clarify this in the revised manuscript.

Q3: Lack of details regarding the MRI task setup.
Answer: We uniformly sample the k-space data following previous work DiffusionMBIR [4], with a simple visualization shown in Figure 1 (IV part) the Mask image. We will revise the manuscript to provide more details on the MRI task setup.

[1] Cai, Yuanhao, et al. "Radiative gaussian splatting for efficient x-ray novel view synthesis." European Conference on Computer Vision. Springer, Cham, 2025.
[2] Zha, Ruyi, et al. "$R^2$-Gaussian: Rectifying Radiative Gaussian Splatting for Tomographic Reconstruction." arXiv preprint arXiv:2405.20693 (2024).
[3] Fu, Xueming, et al. "3DGR-CAR: Coronary artery reconstruction from ultra-sparse 2D X-ray views with a 3D Gaussians representation." International Conference on Medical Image Computing and Computer-Assisted Intervention. Cham: Springer Nature Switzerland, 2024.
[4] Hyungjin Chung, et al. "Solving 3d inverse problems using pre-trained 2d diffusion models." Conference on Computer Vision and Pattern Recognition, 2023.

Thank you for your prompt responses. However, my major concerns are not addressed. Firstly, this work aims to propose a novel method for medical inverse imaging, so visual comparisons between the proposed method and baselines are required. Currently, Fig. 5 in the Appendix only illustrates the reconstructions of the proposed method across iterations. Moreover, you mention that "Our proposed GBIR is not based on 3D Gaussian Splatting." However, the key formulas of the proposed GBIR (Eq. 3 and Eq. 4 in this submission) appear to align with the core concepts of R2-Gaussian based on 3D GS (Eq. 3 and Eq. 4 (https://arxiv.org/pdf/2405.20693)) in essence. Please provide further clarification on these points.

Thank you very much for your additional feedback!

Volume Visualization: We show the visualization of different methods at the Anonymous Website.

Clarification of GBIR: We would like to clarify that our work is fundamentally based on the classical Gaussian distribution, introduced by Carl Friedrich Gauss in 1809, and not derived from 3D Gaussian Splatting (3DGS).

1. Point 1: Our formula is based on Gauss's theory
   The exponent term of Gaussian distribution is given by: $G(\mathbf{x}, {\mu}, {\Sigma}) = \exp\left(-\frac{1}{2}(\mathbf{x} - {\mu})^\top {\Sigma}^{-1} (\mathbf{x} - {\mu})\right)$ where $\mathbf{x} \in \mathbb{R}^d$ is a 3D point, $\mu$ is the mean, and $\Sigma$ is the covariance matrix. This formula is the Gaussian's exponent, a standard result in probability theory, not specific to 3DGS.

2. Point 2: Our Reconstruction Target is Based on GMM
   We formulate the target volume as: $\mathbf{V} = \sum_{i=1}^n G(\mathbf{x}, {\mu}_i, {\Sigma}_i) \cdot I_i$ where $I_i$ is the intensity of the $i$-th Gaussian, which is learnable. This formulation is a variant of the Gaussian Mixture Model (GMM). While GMMs use a sum of Gaussians with normalized weights that sum to one, our formulation uses learnable intensities $I_i$​ as weights. Thus, our reconstrcution target is a learnable-weight version of the GMM, which is a standard approach in machine learning for modeling distributions, and is not specific to methods like 3DGS or $R^2$ Gaussian.

3. Point 3: Different Reconstruction and Optimization
   Reconstruction: Our method discretizes Gaussians into a 3D grid to directly reconstruct the target volume in an end-to-end manner. We maintain differentiability through contribution alignment and significantly accelerate reconstruction using 3-sigma truncation and efficient large Einstein summation decomposition.

   Optimization: We perform optimization in the measurement domain by transforming the reconstructed volume into the measurement space and minimizing the difference between the transformed volume and the measured data. This approach is fundamentally different from 3DGS, which optimizes by randomly rendering a view, then minimizing the difference between the rendered view and the ground truth image. $R^2$-Gaussian also optimizes in the image domain: and they render the X-ray projection and compute the loss with the ground truth image.

Best Regards.

Thanks sincerely for your additional discussion and question. We will include the visualization in the revised manuscript.

Innovation of scene representation:

Yes, 3DGS is the pionner of using Gaussians in representing the scene, and we also stated that we were inspired by this excellent work in our Introduction Section. In the revised manuscript, we will further emphasize our innovative reconstruction and optimization processes.

Optimizing in image domain or measurement domain?

We are confident that both 3DGS and $R^2$-Gaussian optimize in the image domain. Below, we explain why.

For 3D Gaussian Splatting:
In 3DGS, the optimization process aligns rendered images with target images by minimizing their differences. During each iteration, a viewpoint is randomly selected, an image is rendered, and the loss is computed between the rendered and ground-truth images. This confirms that the optimization occurs in the image domain.

We show how they optimize in the image domain below. (Link of Code)

Select the viewpoint, and render the image from Gaussians

render_pkg = render(viewpoint_cam, gaussians, pipe, bg, use_trained_exp=dataset.train_test_exp, separate_sh=SPARSE_ADAM_AVAILABLE)

Get the rendered image

image, viewspace_point_tensor, visibility_filter, radii = render_pkg["render"], render_pkg["viewspace_points"], render_pkg["visibility_filter"], render_pkg["radii"] 
if viewpoint_cam.alpha_mask is not None:
    alpha_mask = viewpoint_cam.alpha_mask.cuda()
    image *= alpha_mask

Then the loss is computed with the rendered image and the ground truth image, which is in the image domain.

gt_image = viewpoint_cam.original_image.cuda()
Ll1 = l1_loss(image, gt_image) # 
if FUSED_SSIM_AVAILABLE:
    ssim_value = fused_ssim(image.unsqueeze(0), gt_image.unsqueeze(0))
else:
    ssim_value = ssim(image, gt_image)
loss = (1.0 - opt.lambda_dssim) * Ll1 + opt.lambda_dssim * (1.0 - ssim_value)

For the $R^2$-Gaussian:
In $R^2$-Gaussian, they optimize the loss with rendered projection and measured projection in Formula 8. Below is their implementation and we explain how they optimize in the image domain below. (Link of Code)

Select a viewpoint and render projection data $I_r$ in this selected viewpoint

# Get one camera for training
if not viewpoint_stack:
    viewpoint_stack = scene.getTrainCameras().copy()
viewpoint_cam = viewpoint_stack.pop(randint(0, len(viewpoint_stack) - 1)) # Similar to 3DGS to select a viewpoint
render_pkg = render(viewpoint_cam, gaussians, pipe) # Render X-ray projection in this selected viewpoint

Get the rendered projection data $I_r$

image, viewspace_point_tensor, visibility_filter, radii = (
    render_pkg["render"],
    render_pkg["viewspace_points"],
    render_pkg["visibility_filter"],
    render_pkg["radii"],
)

Then Compute the loss between rendered projection $I_r$ and ground truth projection $I_m$, both are X-ray images instead of measurement. Note that there is not any Radon Transform involved in transforming the image domain into measurement domain.

gt_image = viewpoint_cam.original_image.cuda() 
loss = {"total": 0.0}
render_loss = l1_loss(image, gt_image) 

Thus, we are sure that the 3DGS and $R^2$-Gaussian are optimized in the image domain instead of measurement domain, and they do not involve any Radon Transform to project the image domain into measurement domain.

We also show how our GBIR works:
On the contrary, our GBIR directly optimize under the measurement domain, i.e., without the need of image for supervision. We show how we optimize in the measurement domain below.

Directly reconstruct the volume from the Gaussians Set，without any rendering process.

reconstruct_volume, impact_shape, std = reconstruct_3d_volume(gaussians, gaussians.volume_shape) 

Project the volume into measurement domain

projs = gaussians.Fan_ray_trafo(reconstruct_volume)

Compute the L1-Loss in the measurement domain

Ll1 = l1_loss(projs, gaussians.gt_projs)

As shown above, our Reconstruction as well as Optimization are innovative: we did not involve any viewpoint selection / rendering process in the reconstruction. Instead, we propose a highly efficient reconstruction directly from a set of Gaussians; and our optimization is completely done in the measurement domain, without any supervision signal in the image domain.

Thank you sincerely for your time and valuable feedback. Below are our responses to your comments:

Q1: The discussion on previous work on Compressed-Sensing MRI is very limited.
Answer: We appreciate your suggestion. Due to space constraints, we did not detail the compressed-sensing MRI literature extensively, as our primary focus was on the proposed method and dataset. In the revised manuscript, we will include a more comprehensive discussion of related work to address this concern.

Q2: Target volume in MRI is not real-valued but complex-valued.
Answer: We apologize for the confusion. This paper focuses solely on reconstructing the magnitude images, which are real-valued, as shown in Figure 6. We will clarify this in the revised manuscript to avoid further misunderstandings.

Q3: Measurements are complex-valued k-space measurements for MRI. How do you calculate SSIM between measurements?
Answer: In our implementation, we calculate SSIM separately for the real and imaginary parts of the complex-valued measurements:

ssim_loss = (ssim(projs.real, gt_projs.real) + ssim(projs.imag, gt_projs.imag))

This structural similarity is used in the loss function to guide the optimization process, and it works well in practice. To address this concern, we will include an ablation study in the revised manuscript to demonstrate the necessity of SSIM in the loss function.

Q4: The dataset used (BRATS) does not contain raw k-space data.
Answer: You are correct that BRATS does not include raw k-space data. All methods, including ours, use simulated k-space data under identical settings, following the baseline DiffusionMBIR [1]. We will clarify this point in the revised manuscript.

Q5: No details regarding the undersampling parameters were provided.
Answer: We followed the uniform k-space sampling strategy outlined in previous work [1]. A simple visualization of the sampling mask is shown in Figure 1 (IV part) of the manuscript. We will revise the manuscript to include more detailed explanations of the undersampling parameters and task setup.

Q6: Quantitative results in SV-CT and CS-MRI tables (Tables 4 and 5) only provide average values.
Answer: Thank you for this suggestion. While we conducted over 20 tables of experiments and utilized significant computational resources, providing standard errors for all baselines is infeasible. However, we will include the standard error for our proposed method in the revised manuscript to enhance the robustness of our results.

Q7: The baselines used are relatively weak and may not represent the state-of-the-art.
Answer: We appreciate this feedback and will incorporate additional state-of-the-art baselines in the revised manuscript to ensure a more comprehensive comparison.

Q8: How do you choose the $\lambda_s$ in Equation 12?
Answer: In our implementation, we set $\lambda_1=1$, $\lambda_2=500$, and $\lambda_3=500$. This information will be added to the revised manuscript for clarity.

Q9: DuDoRNet is not mentioned anywhere in the text but is used in the table.
Answer: Thank you for pointing this out. We will include a brief description of DuDoRNet in the revised manuscript to ensure consistency between the text and tables.

[1] Hyungjin Chung, et al. "Solving 3d inverse problems using pre-trained 2d diffusion models." Conference on Computer Vision and Pattern Recognition, 2023.

Thank you so much for your patient review and additional comments! I once thought that our setting was good since our team followed a heated published baseline. After reading your and 4ada's review, I thought maybe our adopted settings were not in line with the mainstream.

Q2: Without complex-valued $x$, one cannot go to the measurement domain by taking the Fourier transform.
Answer: We used the process below to transform the real-valued volume $x$ into complex-valued measurement $y$:

    volume = torch.rand(256, 1, 256, 256) # dtype=torch.float32  
    kspace = torch.fft.fftshift(torch.fft.fft2(volume), dim=[-1, -2]) # dtype=torch.complex64

We followed previous work to simulate the k-space data in the BRATS dataset in this way. As shown, the real-valued volume $x$ is indeed transformed into a complex-valued measurement $y$. We will find a more appropriate setting in the revised work.

Q3: $\text{1-SSIM}$ can be used as a loss function, but the formula 12 does not make sense since $M$ are k-space measurements in MRI, not images.
Answer: Yes, you are right that $\text{1-SSIM}$ is usually used as the loss in the image domain instead of the measurement domain. Our idea of using such loss to measurement domain is due to: the only supervision signal is the measurement, and we try to keep the structural similarity between measurement $y$ and recovered measurement $\hat y$. We studied this loss term and found it did improve the performance. We will further investigate the reason behind this phenomenon to justify whether it is reasonable to add such a loss term. Thank you for your questions and inquiries about this.

Q4: This is not a proper choice, as there are many publicly available k-space datasets.
Answer: We maintained the identical setting as our baseline to make the fair comparisons as possible, but only to find the sad fact that this setting is not in line with the mainstream. Your suggestion is good, and we will seek a better setting and revise this work accordingly. Thanks for your time and patience again.

Best Regards.