Unsupervised Trajectory Optimization for 3D Registration in Serial Section Electron Microscopy using Neural ODEs

Weaknesses 1

Thank you for your suggestion. We have explored the speed–accuracy trade-off for large-scale EM data. Our optimization significantly reduces runtime (from 2.46s to 0.95s), while preserving SOTA SSIM accuracy (see our responses to Cons 3 from Reviewer L5DT).

Weaknesses 2

We agree that trajectory tracking plays an important role in the overall pipeline performance. Our current design adopts a lightweight 2D U-Net to enable fast and memory-efficient tracking between adjacent slices, but we acknowledge that its limited temporal context may cause errors in ambiguous regions. To address this, we explored two strategies to improve performance under more challenging scenarios (see our responses to Q.2 from Reviewer Luv5). But we acknowledge the importance of modeling longer temporal context, and we plan to explore trajectory models with extended time horizons to handle challenging scenarios with low texture or strong noise.

Moreover, to evaluate robustness, we tested our method on two additional real datasets (see responses to Cons 1 from Reviewer L5DT). All corresponding results will be included in the revision.

Weaknesses 3

We apologize for the typo in Eq. (10) of the main text (implementation and experiments were correct). The original equation:

where $\mathcal{F}_\theta(\mathcal{P}^{(t_0)})$ denotes a shorthand for $\text{ODESolver}$

should be corrected to:

Define

1. Gauss-Seidel:
   Converges under standard conditions for SPD linear systems, solving the quadratic objective

via iterations

with spectral radius $\rho(T) < 1$.

2. Neural ODE:
   Then the ODE has a unique solution $\mathcal{P}(t)$ on $[t_0, t_N]$, and numerical solvers yield approximations with bounded error:

3. Adaptive Blending:

• Both modules converge individually.
• Blending weight $\alpha$ increases monotonically as error decreases via a sigmoid function.
• If $\|\mathcal{P}^* - \mathcal{P}^\dagger\|$ is sufficiently small, the combined update is stable and converges.

Under these conditions, the adaptive blending update $\mathcal{P}^{(t_N)}$ converges reliably.

Weaknesses 4

We thank the reviewer for emphasizing the importance of parameter $\beta$ in Eq. (11). It controls the sigmoid steepness and thus the transition speed between Gauss-Seidel and Neural ODE modules. A larger $\beta$ makes $\alpha$ shift faster from 0 to 1, enabling a quicker switch, while a smaller $\beta$ yields a smoother transition. We will add convergence curves and performance comparisons for different $\beta$ values in the revised version.

Ablation studies were also repeated on four additional datasets (see responses to Weaknesses 2 from Reviewer zfxm).

Weaknesses 5

We employ Gaussian-smoothed random displacement fields to ensure fair comparison. Specifically, we strictly follow the deformation protocols established in widely used baselines such as GaussReg, EFSR, and EMReg, which introduce slice-wise perturbations that closely mimic real-world misalignments observed in ssEM data.

We agree with the reviewer that this model may not fully represent the diverse range of real-world artifacts resulting from physical sectioning. Therefore, in addition to the four real FemFlyBrain blocks included in the supplementary material, we further conduct experiments on six blocks from two widely used real datasets (see our responses to Cons 1 from Reviewer L5DT), providing stronger validation under real-world scenarios. Additional results will be included in the revision.

Weaknesses 6

Thank you for the thoughtful question. Our trajectory inversion step aims to obtain displacement vectors defined on a pixel-aligned grid, which are then used to construct the full deformation field.

It is important to note that this inversion operates on displacements rather than pixel intensities. But it does not distribute a single point’s value across multiple pixels, and thus does not introduce the kind of smoothing effect typically associated with image-based interpolation. The resulting deformation field is applied once to warp the input stack using STN [1]. Therefore, bilinear splatting is only used as a technical means to discretize continuous displacements, and it does not additionally degrade spatial precision.

[1] Spatial transformer networks. NeurIPS 2015

Q.1

Thank you for the suggestion. We have added runtime comparisons of Euler, Dopri5, and RK4 solvers on datasets beyond EFPL (see Table 13). The results show that while Dopri5 and RK4 provide marginal improvements in performance compared to Euler (e.g., SSIM increases from 0.872 to 0.877 and 0.876 on the mouse heart dataset), they come at a significant computational cost (runtime increases from 2.45s to 8.33s and 7.68s, respectively). Given this trade-off, the Euler solver offers a better balance between accuracy and efficiency in our pipeline. For the training of Neural ODEs, we already employ the adjoint sensitivity method to enable memory-efficient and scalable gradient computation.

Table13

Q.2

Thank you for your suggestion. We explored replacing the original U-Net with more advanced architectures while keeping the same training pipeline (see our responses to Cons 2 from Reviewer L5DT). The results show that incorporating more advanced models can further improve performance, yielding 0.3–0.5 higher SSIM scores.

Different dense prediction modules don’t require retraining the trajectory smoothing module. However, sparse prediction are unsuitable because our trajectory inversion uses bilinear splatting to get displacement vectors on the pixel grid. Sparse trajectories can cause missing pixel values and reduce registration quality. Since accurate registration needs a dense, continuous deformation field, dense prediction modules are preferable.

We have added ablation studies on trajectory sampling intervals (see our responses to Cons 3 from Reviewer L5DT). The results demonstrate that overly sparse trajectory sampling leads to performance degradation (e.g., SSIM drops from 0.874 to 0.831), especially for EM images where subpixel-level alignment is critical.

Q.3

Thank you for your insightful question. We have provided a theoretical analysis regarding the convergence of the combined Gauss-Seidel + Neural ODE solver in our above response to Weaknesses.3. As for the behavior of the weighting coefficient $\alpha$, it is dynamically modulated by the relative residual through a sigmoid-shaped function controlled by the parameter $\beta$. Specifically, a larger $\beta$ causes $\alpha$ to transition more rapidly from 0 to 1 once the relative error has sufficiently decreased, thereby accelerating the switch from the linear Gauss-Seidel module to the nonlinear Neural ODE module. In contrast, a smaller $\beta$ leads to a more gradual and smoother transition (see our above response to Weaknesses 4).

Q.4

We thank the reviewer for raising this important question. The $(\lambda, \mu)$ search ranges were initially determined through a grid search on a small subset of the EFPL dataset during preliminary experiments. To assess generalization, we further conducted ablation studies on additional datasets (see responses to Weaknesses 2 from Reviewer zfxm) and consistently observed clear performance peaks around the same parameter values. These results suggest that our choice of $(\lambda, \mu)$ offers good generalization across datasets, reducing the need for extensive re-tuning when applying our method to new domains.

Q.5

We thank the reviewer for the suggestion. We have added new results on four datasets with varied inter-slice intervals to simulate missing sections (see Table 14). Due to time constraints, we were unable to retrain the trajectory tracking module, resulting in reduced performance as the slice interval increases. However, based on our ablation studies on the tracking module (see our response to Cons 2 from Reviewer L5DT), we believe that training on anisotropic data and replacing the tracker with more advanced networks can effectively mitigate this issue.

As for folded sections, handling them is particularly challenging because it requires not only detecting the folds but also unfolding the distorted tissue, which is a complex nonlinear operation [2]. In practice, it is common to treat folded slices as missing sections and exclude them from the registration process to avoid introducing errors [3]. We will include these clarifications and relevant discussions in the revision.

[2] A unified deep learning framework for ssTEM image restoration. IEEE TMI, 2022.

[3] A complete electron microscopy volume of the brain of adult Drosophila melanogaster. Cell, 2018.

Table14

I appreciate the authors’ detailed rebuttal. The authors have addressed key weaknesses such as adding more data validation and improving runtime trade-offs, but shared concerns remain about deployment practicality, scalability, and clarity of writing.

Weaknesses 1

Thank you for the question. The “previous works” in Figure 1 refer to pairwise registration methods. The demo in Figure 1 uses the EFPL dataset, and the color-coded represents deformation field. We apologize for the missing reference to Figure 3(a), which is intended to contrast grid interpolation and splatting. We will clarify these figures with more explicit explanations in the revision.

Weaknesses 2

Thank you for the suggestion. We have extended our ablation studies to four additional datasets, as presented in Tables 8–11 below. The consistent performance patterns observed across all datasets demonstrate the generalizability and robustness of our method beyond a single test set. In addition, we have included several other ablation experiments in response to suggestions from other reviewers (see our responses to Cons 2 from Reviewer L5DT). We will include all results and discussion in the revised version.

Table 8

Table9

Table10

Table 11

Weaknesses 3

Thank you for your question. As reported in Table 2 of the Supplementary Materials, our method achieves the best accuracy with an SSIM of 0.87, though the runtime is 2.46s. To further improve efficiency, we have now added new experiments focused on optimization strategies. We significantly reduce the runtime (from 2.46s to 0.95s) while maintaining SOTA-level accuracy. This makes our approach more applicable for large-scale EM use cases. Please refer to our response to Cons3 from Reviewer L5DT for detailed experimental results and analysis.

Weaknesses 4

Thank you for your comments. We addressed challenging scenarios like high noise and nonlinear deformations by: 1) data augmentation by adding noise during training and 2) upgrading the baseline tracker with advanced optical flow networks (see our response to Q.2 from Reviewer Luv5). The results, e.g., on the mouse heart dataset with 10% noise, demonstrate that both strategies significantly improve performance: data augmentation improves SSIM from 0.832 to 0.852, and upgrading the tracking network further boosts it to 0.883.

For varying nonlinear deformations, Table 12 shows our method consistently outperforms GaussReg (e.g., SSIM on mouse heart drops from 0.832 to 0.742 for GaussReg, but only from 0.872 to 0.803 for ours).

Our biological prior—based on the smooth manifold nature of structures—introduces temporal smoothness and spatial consistency losses (see Line 140 of the main text) to ensure plausible deformations. Ablation studies in Table 2 of the main text and additional datasets beyond EFPL (see our above response to Weaknesses 2) confirm the effectiveness and robustness of this design.

We agree that using a stronger tracking network improves results, with advanced optical flow models boosting SSIM by 0.3–0.5 (see our response to Cons 2 from Reviewer L5DT). However, our work presents a novel paradigm, not just an incremental upgrade (as noted by Reviewer TWfx). A simple network already demonstrates great potential of our framework, though it leaves room for further improvements with better tracking components.

We optimized the speed–accuracy trade-off for large-scale EM data, reducing runtime from 2.46s to 0.95s while maintaining SOTA SSIM accuracy (see our response to Cons 3 from Reviewer L5DT). Scalability and memory for TB-scale datasets are also addressed (see our below response to Q.4). All results will be included in the revised version.

Table 12

Q.1

Thank you for raising this important question. For the synthetic experiments, we adopt the high-quality volumes from the OpenOrganelle dataset as the ground-truth references. These volumes are acquired with isotropic resolution and minimal distortion, making them suitable for constructing controlled benchmarks.

To simulate realistic z-axis deformations, we synthetically warp the volume along the axial direction using smooth deformation fields generated by randomly sampling Gaussian noise fields (see Eq. (1) in the Supplementary Materials). Independent perturbations are applied to each slice (except for the first) to introduce realistic inter-slice misalignments (see Section A2.1 for details in the Supplementary Materials).

To ensure the realism of the synthetic distortions, we strictly follow the deformation protocols established in widely used baselines such as GaussReg, EFSR, and EMReg, which introduce slice-wise perturbations that closely mimic real-world misalignments observed in ssEM data. This benchmark design enables controlled evaluation while preserving biological plausibility.

Q.2

Thank you for the question. To the best of our knowledge, EMReg, as presented in the original paper, focuses solely on pairwise registration, and neither the method nor the reported results extend to full-stack or global 3D registration. Therefore, EMReg has not been adapted to optimize over larger blocks with global regularization. If the reviewer is aware of subsequent work where EMReg has been applied for global regularization across larger volumes, we would be happy to analyze and compare accordingly.

Regarding baseline implementations, we used the official EMReg implementation and extended it to 3D registration via sequential pairwise registration. For TrakEM2, we employed the FIJI “Elastically Align Stack” plugin. For SEAMLeSS, we used the official Corgie package for 3D registration. For EFSR and GaussReg, we followed the settings of the original papers. All baselines were trained on the same data as our method, with hyperparameters set according to their respective papers. We will clarify all baseline configurations and resolve citation inconsistencies in the revised version.

Q.3

We appreciate your mention of SOFIMA, whose distributed block processing strategy has been inspiring to us. After carefully reviewing its source code, we found that it employs an "elastic mesh optimizer" for registration but acknowledges that the optimization proceeds sequentially, section by section. This approach differs fundamentally from our design of global trajectory optimization.As Reviewer TWfx noted, we conceptually redefine registration as a global trajectory optimization problem—an elegant and significant contribution that directly addresses the fundamental issue of error accumulation inherent in pairwise methods. Due to the time constraint for rebuttal, we will discuss this in the revised version.

Q.4

Thank you for the suggestion. We optimized the speed–accuracy trade-off, reducing runtime from 2.46s to 0.95s while maintaining SOTA SSIM accuracy (see our response to Cons 3 from Reviewer L5DT). Our framework can scale to TB-scale ssEM datasets using a parallel, block-wise strategy like SOFIMA: 1) split volumes into chunks; 2) run parallel trajectory tracking and inversion on multi-GPUs; 3) stitch local deformation fields; 4) apply the global deformation.

As noted in Section A2.3 of Supplementary Materials, for datasets containing thousands of slices, we split them into blocks of 100 slices, which can be processed comfortably on a 40GB GPU. We will detail this scalable design in the revision.

Limitation 1

Thank you for your suggestion. We have explored the speed–accuracy trade-off for large-scale EM data. Our optimization significantly reduces runtime (from 2.46s to 0.95s), while preserving SOTA-level SSIM accuracy (see our responses to Cons 3 from Reviewer L5DT).

Limitation 2

Thank you for your suggestion. TrakEM2 is a widely used EM image processing toolkit that supports both stitching and registration. However, existing frameworks—including TrakEM2—typically adopt a sequential pipeline that performs stitching first, followed by registration. Our method is also compatible with such pipelines.

But we agree that integrating stitching and registration into a unified framework is a promising direction. A potential solution is to leverage shared encoder features between the stitching and trajectory tracking modules, enabling joint optimization.

Moreover, we agree that stitching errors can affect the trajectory tracking. To assess the impact, we added evaluations under more challenging conditions (see our responses to Q.2 from Reviewer Luv5), and further replaced the tracking module with more advanced networks to improve robustness against artifacts (see our responses to Cons 2 from Reviewer L5DT). We will incorporate all these clarifications in the revised version.

Dear Reviewer zfxm,

Thank you very much for your time and valuable comments on our work.

As the author-reviewer discussion phase will draw to its end soon, we would greatly appreciate it if you could let us know whether our rebuttal has addressed your concerns, or if any part still requires further clarification.

Best regards,

Authors

Cons 1

Thank you for the valuable suggestion. We already showed registration results on four $1024 \times 1024 \times 1299$ subvolumes from different regions of the real FemFlyBrain dataset in Figures 12–13 of our Supplementary Materials. As you can see, our method achieves improved structural continuity and better restoration of correct axial structures compared to baselines.

Following your suggestion, we have now additionally evaluated our method on two $1024 \times 1024 \times 700$ blocks from the FAFB dataset [1] and four $1024 \times 1024 \times 500$ blocks from the mouse cortical dataset [2]. Due to the lack of established quantitative metrics for this task on real data and the restriction of not being able to include figures in rebuttal, we will include corresponding results and visual comparisons in the revised version.

[1] A complete electron microscopy volume of the brain of adult Drosophila melanogaster, Cell, 2018.
[2] Saturated reconstruction of a volume of neocortex, Cell, 2015.

Cons 2

We agree with the reviewer that accurate tracking is a critical foundation for the subsequent optimization stages. In response, we have now added ablation studies on the trajectory module, including:
(i) replacing it with a stronger dense tracking network (Table 1), and
(ii) a sensitivity analysis of the key hyperparameter $\lambda$ (Table 2).

To ensure generality, these experiments are also conducted on datasets beyond EFPL. Table 1 shows that incorporating more advanced optical flow models can further improve the performance, yielding 0.3--0.5 higher SSIM scores. Table 2 demonstrates that our method remains robust and maintains SOTA-level accuracy across a wide range of $\lambda$ values (between 1.5 and 7.5).

We would like to emphasize that our work introduces a novel paradigm rather than merely an incremental improvement (as aptly noted by Reviewer TWfx). To highlight this, we intentionally adopted a simple baseline network architecture in the original manuscript, which leaves significant room for further performance gains by integrating more sophisticated tracking components. All the additional results will be included in the revised version.

Table 1.

RAFT: Raft: Recurrent all-pairs field transforms for optical flow. ECCV, 2020

GAFlow: GAFlow: Incorporating Gaussian Attention into Optical Flow, ICCV, 2023

SEA-RAFT: SEA-RAFT: Simple, Efficient, Accurate RAFT for Optical Flow, ECCV 2024

Table 2.

Cons 3

Thank you for the insightful suggestion regarding runtime efficiency. In response, we have added new experiments focused on efficiency optimization through:
(i) sparse trajectory sampling (every 2/3/4 pixels) (Table 3), and
(ii) resolution reduction via down-sampling followed by up-sampling (Table 4).

Table 3 shows that sparse sampling significantly reduces runtime (from 2.46s to 0.95s), while preserving SOTA-level accuracy, with SSIM only decreasing slightly from 0.874 to 0.831.
Table 4 demonstrates that reducing resolution also improves speed (from 2.46s to 0.91s), with a drop in accuracy (SSIM from 0.874 to 0.827), which still outperforms the GaussReg baseline. The observed performance drop may be attributed to the use of a simple bilinear up-sampling strategy. We believe that combining sparse trajectory sampling with resolution reduction via down-sampling followed by up-sampling can achieve a better trade-off between runtime and accuracy.We will include the above results and further discussion in the revised version.

For more situations without labels, we also include more experiments on real datasets and compare with GaussReg (see our response to Cons 1 above). We will include the complete results and extended discussion in the revised version.

Table 3.

Table 4.

Cons 4

We appreciate your valuable suggestion. As reported in Table 3 of our Supplementary Materials, we have already provided Dice scores for quantitative evaluation. In accordance with your request, we additionally report the Folds metric in Table 5 to assess the diffeomorphic property of the deformation field. The results show that our method achieves the lowest Folds values, consistently around 0.01, whereas most baseline methods report values above 0.1. This indicates significantly fewer grid foldings and better preservation of topological consistency in the deformation fields produced by our approach.

Table 5.

Q.1

Thank you for your valuable suggestion. As stated in the main text (Table 1), all six test datasets are unseen ssEM datasets. We train a single model on training data and evaluate it across these six unseen datasets to demonstrate the generalization ability of our method. Compared to traditional methods TrakEM2, our approach achieves significantly improved performance, for instance, on the mouse heart dataset, the SSIM improves from 0.57 to 0.87.

For real-world datasets, we also adopt a similar protocol: training on a small region and testing on separate, unseen regions (as described in Line 110, Section A2.3 of our Supplementary Materials). The visualizations in Figures 12–13 of Supplementary Materials demonstrate that our method better restores the correct axial structure and generalizes well to new regions that were not observed during training.

Following your suggestion, we additionally include experiments on 6 blocks from two more real-world ssEM datasets (see our response to Cons 1 above). We also provide ablation studies on datasets beyond EFPL (see our response to Weaknesses 2 from Reviewer zfxm), further validating the robustness and generalizability of our framework.

Q.2

Thank you for your helpful suggestion. We now provide a sensitivity analysis of the key hyperparameter λ in Table 2. The results show that our method remains robust and maintains SOTA-level accuracy across a wide range of λ values (from 1.5 to 7.5).

In addition, regarding the tracking module, we replace the original 2D Unet network with more advanced optical flow models: RAFT (5.3M parameters), GAFlow (10.1M), and SEA-RAFT (4.9M). Table 1 shows that incorporating these more powerful tracking networks can further boost performance, achieving 0.3–0.5 improvements in SSIM scores. Please also refer to our response to Cons 2 above for additional detailed results and discussions.

Q.3

Thank you for your question. For large-scale long-slice ssEM datasets, we adopt a divide-and-conquer strategy by splitting large volumes into manageable blocks. Our method estimates deformation fields in parallel for each block and then fuses them to reconstruct the full field. Please refer to our response to Q.4 from Reviewer zfxm for more details.

For handling long slices, we process them in blocks (e.g., 100 slices per batch), as described in Line 113, Section A2.3 of our Supplementary Materials.

To address the inference burden, we introduce two efficient strategies as discussed in our response to Cons 3 above: 1) sparse trajectory sampling, and 2) downsample-then-upsample resolution. Experiments show that these strategies significantly reduce inference time from 2.46s to 0.95s, while preserving SOTA accuracy.

Q.4

We appreciate your valuable suggestion. In Table 3 of our Supplementary Materials, we have already provided Dice scores for quantitative evaluation.

In accordance with your request, we additionally report the Folds metric in Table 5 above to assess the diffeomorphic property of the deformation field. The results show that our method achieves the lowest Folds values, consistently around 0.01, whereas most baseline methods report values above 0.1. This indicates significantly fewer grid foldings and better preservation of topological consistency in the deformation fields produced by our approach.

Weaknesses 1

Thank you for raising concerns on novelty. While Neural ODEs have been applied in many other domains, to the best of our knowledge, we are the first to introduce Neural ODEs for 3D registration in ssEM data.

More importantly, our key contribution lies not merely in the use of Neural ODEs, but in the conceptual reformulation of registration as an unsupervised global trajectory optimization problem—an innovation that directly addresses the long-standing issue of error accumulation in traditional pairwise methods (as succinctly summarized by Reviewer TWfx). Neural ODEs are simply a natural and effective mechanism for realizing this paradigm.

Weaknesses 2

Thank you for the valuable suggestion. We already showed registration results on four 1024 × 1024 × 1299 subvolumes from different regions of the real FemFlyBrain dataset in Figures 12–13 of our Supplementary Materials. As you can see, our method achieves improved structural continuity and better restoration of correct axial structures compared to baselines. Following your suggestion, we have now additionally evaluated our method on two blocks from the FAFB dataset and four blocks from the mouse cortical dataset (see our responses to Cons 1 from Reviewer L5DT). We will include the corresponding results in the revision.

Q.1

Thank you for the valuable question. The biological prior refers to the smooth manifold nature of biological structures, where nonlinear deformation primarily arises from gradual local curvature changes (as mentioned at the end of the third paragraph in Section 1 of the main text).

Based on this prior, we introduce two regularization terms in our trajectory optimization module: a temporal smoothness loss and a spatial consistency loss (see Line 140 in the main text). These are designed to promote biologically plausible deformations.

To validate the correctness of this biological prior modeling, we conduct ablation studies in Table 2 of the main text. In addition, we now include the same ablation experiments on two more datasets beyond EFPL (see our response to Weaknesses 2 from Reviewer zfxm). The results show consistent patterns in peak performance improvements, further supporting the robustness of our prior-based design. We will clarify this point more explicitly in the revised version.

Q.2

Thank you for the suggestion. We agree with you that noise and imaging artifacts would affect the trajectory tracking stage. To address this issue, we adopt two strategies: 1) introducing data augmentation during training and 2) replacing the baseline tracker with more advanced optical flow networks. To validate these approaches, we conducted experiments on four datasets with varying noise levels, comparing our method with GaussReg (Table 6).

The results (e.g., on the mouse heart dataset with 10% noise) demonstrate that both strategies significantly improve performance: data augmentation improves SSIM from 0.832 to 0.852, and upgrading the tracking network further boosts it to 0.883. In contrast, GaussReg achieves only 0.788. These results highlight that our method remains robust under noisy conditions. Importantly, this indicates that our work represents a novel paradigm in ssEM registration rather than a mere incremental improvement, with significant room for further gains even using more powerful networks and training tricks.

Table6

SEA-RAFT: SEA-RAFT: Simple, Efficient, Accurate RAFT for Optical Flow, ECCV 2024

Q.3

Thank you for your suggestion. We provided ablation studies in Table 5 of our Supplementary Materials by removing Neural ODEs. The results show that removing the Neural ODE component leads to performance drop, with SSIM decreasing from 0.746 (full model) to 0.721.

To further validate the importance of this module, we additionally performed the same ablation experiments on multiple datasets beyond EFPL (see our response to Weaknesses 2 from Reviewer zfxm). Consistent patterns were observed across these datasets, confirming that the Neural ODEs are critical to the effectiveness of our method.

Q.4

Thank you for the valuable question. GaussReg performs optimization within local slice windows, which inherently limits its ability to maintain global consistency across long-range sections. In contrast, our method adopts a global trajectory-based formulation that enables coherent alignment over large spatial extents. This key advantage is quantitatively supported by the error accumulation scatter plots shown in Figure 4 of the main text and Figures 6–11 of the Supplementary materials, where our method demonstrates reduced error accumulation compared to GaussReg. We will further elaborate on this important distinction in the revision.

Q.5

Thank you for your suggestion. Regarding the effectiveness of our method on more real-world data, we already presented registration results on four 1024 × 1024 × 1299 subvolumes from different regions of the real FemFlyBrain dataset in Figures 12–13 of the Supplementary Materials. As you can see, our method demonstrates improved structural continuity and more accurate restoration of axial anatomical structures compared to the baselines. Following your suggestion, in addition to the four real FemFlyBrain blocks included in the supplementary material, we further conduct experiments on six blocks from two widely used real datasets (see our responses to Cons 1 from Reviewer L5DT), we will report the corresponding results in the revised version.

Q.6

Thank you for the valuable suggestion. We employ the Gauss-Seidel method, a classical iterative solver for linear systems $A\mathbf{x} = \mathbf{b}$ that leverages the most recent updates to accelerate convergence, particularly in sparse, tridiagonal systems.

Our smoothing objective (Eq. (2) in Section A.3 of Supplementary Materials) is:

As derived in Section A.3 of Supplementary Materials, minimizing this energy leads to the following linear system:

This system can be equivalently written as a symmetric tridiagonal linear system:

where:

• $\mathbf{p} = [\mathcal{P}(1), \ldots, \mathcal{P}(N{-}1)]^\top$ is the unknown trajectory,
• $\mathbf{b} = [\lambda \mathcal{P}^{(0)}(1), \ldots, \lambda \mathcal{P}^{(0)}(N{-}1)]^\top$,
• and $A \in \mathbb{R}^{(N{-}1) \times (N{-}1)}$ is a tridiagonal matrix with structure:

Applying the Gauss-Seidel method to this system, the general update rule for the $i$-th variable is:

For our tridiagonal matrix, each row $z$ only has nonzero entries at $z-1$, $z$, and $z+1$. Substituting these into the general Gauss-Seidel formula, the update rule for $\mathcal{P}(z)$ becomes:

where:

• $A_{zz} = 4 + \lambda$,
• $A_{z,z-1} = A_{z,z+1} = -2$,
• and we use the latest available updates for $\mathcal{P}^{(k+1)}(z{-}1)$, along with the previous values for $\mathcal{P}^{(k)}(z{+}1)$.

Substituting these in yields:

For simplicity, we divide both numerator and denominator by 2, resulting in the final explicit Gauss-Seidel update:

This completes the derivation of the iterative update rule used for trajectory smoothing, as presented in Eq.8 of the main text.

Q.7

Thank you for the valuable question. We adopt bilinear splatting (not Gaussian splatting) in the inversion step due to its favorable speed–accuracy trade-off. To validate this choice, we now provide a comparison among bilinear, trilinear, and classical radial basis function (RBF) interpolation in Table 7.

The results show that bilinear interpolation achieves an SSIM of 0.87 within only 0.008 seconds. In contrast, classical RBF interpolation, even when accelerated by GPU, takes 1.2 seconds but yields only marginal improvement (SSIM = 0.89). While learning-based interpolation methods may potentially yield further improvements, we leave this for future work due to time constraints for rebuttal.

Table7

Dear Reviewer Luv5,

Thank you very much for your time and valuable comments on our work.

As the author-reviewer discussion phase will draw to its end soon, we would greatly appreciate it if you could let us know whether our rebuttal has addressed your concerns, or if any part still requires further clarification.

Best regards,

Authors