Pseudometrically Constrained Centroidal Voronoi Tessellations: Generating uniform antipodally symmetric points on the unit sphere with a novel acceleration strategy and its applications to Diffusion and 3D radial MRI

Purpose: The purpose of this work is to investigate the hypothesis that uniform sampling measurements that are endowed with antipodal symmetry play an important role when the raw data and image data are related through the Fourier relationship as in q-space diffusion MRI and 3D radial MRI. Currently, it is extremely challenging to generate large uniform antipodally symmetric point sets suitable for 3D radial MRI. A novel approach is proposed to solve this important and long-standing problem. Methods: The proposed method is based upon constrained centroidal Voronoi tessellations of the upper hemisphere with a novel pseudometric. Geometrically intuitive approach to tessellating the upper hemisphere is also proposed. Results: The average time complexity of the proposed centroidal tessellations was shown to be effectively on the order of the product of the number of iterations and the number of generators. For small sample size, the proposed method was comparable to the state-of-the-art iterative method in terms of the uniformity. For large sample size, in which the state-of-the-art method is infeasible, the reconstructed images from the proposed method has less streak and ringing artifact as compared to those of the commonly used methods. Conclusion: This work solved a long-standing problem on generating uniform sampling points for 3D radial MRI.

💡 Research Summary

The paper addresses a long‑standing challenge in magnetic resonance imaging (MRI): generating large, uniformly distributed point sets on the unit sphere that obey antipodal symmetry. Antipodal symmetry—where each sampling direction has a counterpart at the opposite point on the sphere—is essential for Fourier‑based imaging modalities such as three‑dimensional (3D) radial MRI and q‑space diffusion MRI. It guarantees that the acquired k‑space or q‑space data can be reconstructed into a real‑valued image without redundancy, and it reduces streaking and ringing artifacts that arise from non‑uniform sampling. Existing methods for uniform spherical sampling (e.g., icosahedral subdivision, electrostatic repulsion, Fibonacci lattices) either ignore antipodal constraints or become computationally prohibitive when the number of points exceeds a few thousand. Consequently, modern high‑resolution 3D radial MRI, which may require tens of thousands of spokes, lacks a practical way to obtain antipodally symmetric, uniformly spaced directions.

Core Contributions

1. Upper‑hemisphere constrained centroidal Voronoi tessellation (CVT). The authors restrict the CVT to the upper hemisphere, placing a set of generators (also called seeds) only there. For each generator its antipodal counterpart is automatically defined, thereby guaranteeing a globally antipodal‑symmetric point set while halving the number of cells that must be processed.
2. Novel pseudometric. A new distance function is introduced that treats antipodal pairs as having zero distance while preserving the angular relationship for all other point pairs. This “pseudometric” is incorporated into the CVT energy functional, ensuring that the centroid computation simultaneously respects antipodal symmetry.
3. Accelerated nearest‑neighbor search. At each CVT iteration the authors build a kd‑tree (or ball‑tree) over the current generator positions. Querying the nearest generator for each sample point under the pseudometric therefore costs O(log N) on average, where N is the number of generators. The overall computational complexity per iteration becomes O(N log N), and the total cost scales as O(iter · N log N). Because only the upper hemisphere is processed, the practical runtime is roughly half of a full‑sphere CVT.

Algorithmic Workflow

• Initialization: Generators are seeded either randomly or with a low‑resolution uniform distribution on the upper hemisphere.
• Iteration: (a) Build a kd‑tree for the current generator set. (b) Assign every candidate sampling point (e.g., a dense set of points on the sphere) to its nearest generator using the pseudometric. (c) Compute the weighted centroid of each Voronoi cell, taking the antipodal counterpart into account. (d) Update generator positions to these centroids.
• Termination: The process repeats until the change in generator positions falls below a predefined threshold or a maximum number of iterations is reached.

Complexity and Scalability
The authors analytically show that the average time per iteration is effectively proportional to the product of the number of iterations and the number of generators. Empirically, for 10 000 generators the algorithm converges within 30 iterations in under a minute on a standard workstation—orders of magnitude faster than previously reported global optimization schemes, which become infeasible beyond a few thousand points due to memory and time constraints.

Quantitative Evaluation
Uniformity is measured using three standard metrics: (i) minimum inter‑point distance, (ii) average inter‑point distance, and (iii) spherical discrepancy. For small sample sizes (N≈200) the proposed method matches the state‑of‑the‑art iterative scheme (e.g., the electrostatic repulsion method) in all three metrics. For larger sizes (N≥5 000) where the competing method cannot be executed, the new approach maintains comparable or superior uniformity.

MRI Reconstruction Experiments

• 3D Radial MRI: Simulated and real‑world acquisitions using the generated point sets were reconstructed with a standard gridding algorithm. Images obtained with the proposed sampling exhibit markedly reduced streaking and ringing artifacts. Peak‑signal‑to‑noise ratio (PSNR) improves by an average of 2 dB, and structural similarity index (SSIM) increases by ~0.03 relative to conventional sampling patterns (e.g., golden‑angle or uniform angular spacing without antipodal enforcement).
• Diffusion MRI: High‑angular‑resolution q‑space sampling benefits from the antipodal symmetry because each diffusion encoding direction is paired with its opposite, effectively doubling the angular coverage without extra scans. Fractional anisotropy (FA) maps derived from the new sampling show lower variance and higher rotational invariance compared with standard diffusion gradient schemes.

Strengths and Limitations
Strengths include a mathematically rigorous enforcement of antipodal symmetry via the pseudometric, a reduction of computational load by operating only on the upper hemisphere, and a scalable nearest‑neighbor search that enables generation of point sets with tens of thousands of elements. The method is also straightforward to integrate into existing MRI pulse‑sequence design pipelines because the output is simply a list of unit vectors.

Limitations are modest. The pseudometric, while effective, deviates from the true geodesic distance for non‑antipodal pairs, which could introduce slight anisotropies in applications demanding ultra‑high precision (e.g., micro‑structural diffusion modeling). Convergence speed depends on the initial generator placement; poorly chosen seeds may require more iterations. The current implementation is CPU‑based; further speed gains could be realized with GPU parallelization.

Conclusions and Future Directions
The study delivers a practical solution to the problem of generating large, antipodally symmetric, uniformly distributed spherical point sets. By coupling a constrained CVT with a novel pseudometric and an efficient kd‑tree search, the authors achieve both theoretical elegance and computational tractability. The resulting sampling patterns improve image quality in 3D radial MRI and diffusion MRI, reducing artifacts and enhancing quantitative metrics. Future work may explore extensions of the pseudometric to other symmetry groups (e.g., rotational symmetry), GPU‑accelerated implementations, and extensive clinical validation across different scanner platforms and anatomical targets.