Pith. sign in

REVIEW 3 major objections 7 minor 22 references

Numerically robust tetrahedron-based tomographic forward and backward projectors on parallel architectures

T0 review · 3 major / 7 minor · reviewed 2026-08-14 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read A safety margin keeps tetrahedral CT rays from getting lost.

desk verdict Solid, useful tetrahedral-CT projector with open code; robustness guarantee is asserted, not proven, and validation is inverse-crime. read the letter →

arxiv 1908.06909 v1 pith:FTYPP2AL submitted 2019-08-19 cs.DC

Ander Biguri , Hossein Towsyfyan , Richard Boardman , Thomas Blumensath This is my paper · ORCID
classification cs.DC MSC 65R3268U1065D18 PACS 87.57.Q
keywords computedtomographytetrahedralmeshraytracingforwardprojectionbackprojectionGPUparallelizationMöller-TrumboreR*-tree
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper sets out to make tetrahedral-mesh images usable in X-ray computed tomography by building forward and backward projectors that survive floating-point error. It claims that three ingredients—a graph of face-adjacent tetrahedra, a ray-triangle intersection test with a retry-based safety parameter, and a tree search for the first hit—let a GPU trace rays through arbitrarily shaped convex tetrahedra without missed intersections or infinite loops. If the claim is right, CT reconstruction no longer needs voxel grids: object surfaces and variable resolution can live in the mesh itself, and iterative algorithms can run on far fewer elements. The accompanying reconstruction tests support this by showing that a known-surface tetrahedral mesh avoids the streak artefacts seen in comparable voxel reconstructions.

What carries the argument

The load-bearing mechanism is a propagation loop that walks a ray through a tetrahedral mesh represented as a graph of face-neighbour links. At each tetrahedron, a Möller–Trumbore ray-triangle test is run; when it returns fewer than two intersections, the algorithm retries with a slightly enlarged face—epsilon multiplied by ten each time—until two faces are found, then resets epsilon. The second intersection's face label selects the next tetrahedron, and an R*-tree of boundary tetrahedra finds the entry element so no exhaustive search is needed. This combination carries the claim that arbitrary convex tetrahedral meshes can be traversed reliably in parallel on GPUs.

What would settle it

Generate a convex tetrahedral mesh with one deliberately inverted element or one internal face shared by three tetrahedra; run the single-precision GPU projector over thousands of random view angles. The first ray that stops early or cycles forever marks where the stated robustness fails; a ray aimed exactly through a shared vertex is the sharpest test, since the method must continue without any brute-force search.

Watch

Extended reading notes

Core claim

The central discovery is that the classic Möller–Trumbore ray-triangle intersection test, normally fast but leaky under single-precision floating point, can be made tomography-safe by retrying with a slightly enlarged triangle whenever fewer than two intersection points are found on a tetrahedron. Because the ray is already known to cross the element before the retry, growing the face by a tiny epsilon cannot admit false positives, and the intersection parameter t stays unchanged. The paper demonstrates that this mechanism, combined with a graph of face-neighbour links to propagate the ray element by element and an R*-tree to find the entry tetrahedron, eliminates premature ray termination and infinite neighbourhood loops on high-aspect-ratio meshes. It further shows that double precision in the intersection kernel is required for meshes whose element sizes vary widely, while the rest of the pipeline can stay in single precision.

Load-bearing premise

The propagation walk assumes the mesh is a valid convex tetrahedralization with consistent face-neighbour adjacency—every interior face shared by exactly two tetrahedra, a closed boundary, no inverted or zero-volume elements—and the paper offers no proof that the epsilon retry loop always terminates on the correct pair of faces when a mesh is malformed.

Editorial extensions

If this is right

  • Any iterative CT solver (SART, CGLS, OS-SART) can be run with a tetrahedral image basis, so object boundaries encoded in the mesh survive reconstruction.
  • X-ray projection can be simulated directly from CAD-derived volume meshes, enabling realistic multi-material transmission studies.
  • Known-surface meshes let tomography reconstruct from far fewer elements than voxel grids, and the paper's experiments show this reduces streak artefacts when the surface prior is accurate.
  • The R*-tree initialization grows logarithmically with boundary element count, keeping the first-hit search from dominating on large industrial meshes.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The paper does not prove that the epsilon retry loop always terminates; a natural extension would cap the number of retries and flag the ray for a conservative fallback, turning the heuristic into a bounded algorithm.
  • Because the method treats element size and aspect ratio as free, the same tetrahedral mesh could serve both CT reconstruction and downstream finite-element simulation, removing one segmentation step.
  • A testable prediction follows from the safety-parameter design: meshes with poor element quality should show elevated retry rates per million rays even when no ray actually fails, making the retry counter a cheap mesh-quality diagnostic.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 7 minor

Summary. This manuscript describes a GPU-oriented forward and backward projector for X-ray CT in which the attenuation volume is represented as an unstructured tetrahedral mesh. The algorithm stores tetrahedra with an explicit neighbor graph, initializes each ray through an R*-tree search over boundary elements, and then traverses the mesh using Möller-Trumbore triangle intersections with an adaptive safety parameter epsilon that is increased whenever fewer than two face intersections are found. The authors present reconstructions from simulated projections using tetrahedral and voxel bases, a CAD-model reconstruction, a floating-point precision comparison, and scaling measurements. The central claim is numerical robustness for convex tetrahedral meshes with high-aspect-ratio tetrahedra, with an open-source implementation.

Significance. If the robustness claim holds, the paper fills a practical gap: no currently available openly released tetrahedron-based CT projector handles high-aspect-ratio meshes on GPUs, and the authors do release their code. The experiments on challenging meshes and the observed single-precision failures support the motivation for the work. However, the safety-parameter mechanism is asserted rather than proven, and the main reconstruction validation is an inverse-crime setup, so the central robustness claim is not yet established. The performance characterization and the open-source implementation are useful contributions in themselves.

major comments (3)
  1. [Section II-C and Algorithm 2] The load-bearing claim that the safety parameter cannot introduce false positives is asserted, not proven. In Algorithm 2, epsilon is multiplied by 10 until TetraRayIntersection returns two intersections. For a ray through a vertex or along an edge, or for sliver tetrahedra, the widened barycentric acceptance region (-epsilon, 1+epsilon) can eventually accept a face whose infinite plane is crossed but whose geometric triangle is not crossed by the ray. The stopping condition is 'some face accepts', not 'the true entry/exit face pair was found', so the second t can be a false positive and the neighbor propagation in the following line can step to the wrong element. No bound on epsilon, no argument that the accepted t lies on a true boundary, and no termination proof are provided. Without these, the central claim of robustness for arbitrary tetrahedral meshes is unsupported. I request either a correctness proof with explicit bounds on epsilon, or a redesigned fallback with a bounded verification step, plus stress-test experiments on vertex- and edge-grazing rays and deliberately sliver-containing meshes.
  2. [Section III-A and Figure 6] The reconstruction experiment uses the same tetrahedral mesh to simulate the projections and to represent the reconstruction. As the authors acknowledge, this is an inverse-crime validation: it shows that the projector/backprojector pair is self-consistent, but it cannot validate the geometric accuracy of the ray-tracing against ground truth or against an independent projector. The comparison is also purely qualitative; no error metrics such as projection-domain RMSE, reconstruction RMSE, or boundary error are reported. Please add quantitative validation where data are simulated on one mesh or by an independent ray-tracer and reconstructed on a different mesh, and report per-element or per-region errors.
  3. [Section II-D] The handling of zero-length intersections ('if t2 = t1 check if they need to be swapped') is not specified. Figure 3 describes a configuration where a ray passes through a node shared by several tetrahedra, and the text states that an extra check is needed to avoid backtracking, but the check itself is not defined. This is load-bearing because a wrong choice of the neighbor at t2 propagates the ray along an incorrect path. Please provide the exact decision rule, including how the face index at t2 is chosen when t1 == t2, and why that rule prevents both missed propagation and infinite loops.
minor comments (7)
  1. [Algorithm 1] The hard-coded threshold |a| < 10^-8 is not related to the precision of the mesh or to epsilon, and the comment 'Check if its zero' is unclear; please explain how this threshold is chosen and whether it is independent of coordinate scale.
  2. [Section II-D] The statement 'Our algorithm has only one constraint: the volumetric mesh must be convex' conflicts with the earlier claim of 'no constrains on the triangulation'; please restate the exact geometric and topological assumptions.
  3. [Section II-C1] The word 'imporatnat' should be 'important'.
  4. [Section III-D] The word 'compassion' should be 'comparison'.
  5. [Section II-E] The words 'initialzie' and 'spacial' should be 'initialize' and 'spatial'.
  6. [Algorithm 2] The pseudocode condition 'while not Intersection' is ambiguous because TetraRayIntersection is called in the body; please use a boolean variable and return the two intersection parameters explicitly.
  7. [Section III-C and Figure 9] The floating-point failure is shown only qualitatively as black dots; please report the number or percentage of rays that fail to terminate in single precision versus double precision.

Circularity Check

1 steps flagged · score 2.0 of 10

Projector derivation is self-contained; only an acknowledged inverse-crime validation introduces minor circularity.

  1. other [Section III-A, Image Reconstruction (Figure 6 discussion)]
    "On a tetrahedra basis reconstructions, the same mesh as the one used to generate the data is used... Note that we here find the best reconstruction using tetrahedra basis using the same mesh to generate the projection data and compute the tetrahedral reconstruction. Whilst this is not realistic in real applications, it is done here to 1) show that our method does accurately compute forward and backward projections and thus provides low error reconstructions in the ideal case"

    The synthetic projections are produced with the same tetrahedral forward projector and the same mesh that is then used as the reconstruction basis inside OS-SART. The forward operator used to generate the data and the forward model used for inversion are therefore identical by construction, which is the classic inverse crime. The resulting low-error reconstruction on the known mesh is a self-consistency check of the code rather than an independent validation of the projector's physical accuracy or numerical robustness. The favorable outcome is partly forced by the experimental design, since the data lie exactly in the range of the discrete forward operator being inverted.

full rationale

The core derivation of the tetrahedral projector is not circular. The Moller-Trumbore intersection with the safety parameter epsilon is introduced as a numerical remedy for missed intersections; it is not fitted to data and no parameter is renamed as a prediction. The graph-based neighborhood propagation and R*-tree initialization are standard algorithmic components, and the forward/backprojection pair is defined by the same intersection lengths without importing a uniqueness theorem or self-citation as load-bearing evidence. The paper does cite prior work by the authors, but only for context (e.g., a limited-angle tomography study) and for the TIGRE voxel toolbox used in comparisons; none of these citations carries the central claim. The repeated claim that marginally enlarging a triangle cannot introduce false positives is an unproven robustness assertion, but it is a correctness concern, not a circularity: the algorithm does not define its conclusion in terms of that assertion. The main circular element is the validation in Section III-A, where simulated projections are generated with the same mesh and forward projector used for reconstruction, an acknowledged inverse-crime setup. Because this is a validation limitation rather than a derivation step, and because the paper also tests against a different random mesh and voxel reconstructions, the overall circularity score is low.

Assumptions & free parameters 4 free parameters · 6 assumptions · 0 invented entities

No new physical entities are introduced. The free parameters are algorithmic hyperparameters, not data-fitted physical constants. The axioms are all standard geometry or tomography assumptions except the safety-parameter assertion, which is the load-bearing unproven numerical claim.

free parameters (4)
  • epsilon_initial = 1e-9
    Initial safety parameter in Algorithm 1 and 2; chosen by hand, increased by factor 10 until two intersections are found. No sensitivity analysis is provided.
  • epsilon_multiplier = 10
    Multiplicative increment for the safety parameter; chosen by hand. The paper notes the increment should be gradual but does not justify this specific factor.
  • R_tree_children = min 4, max 10
    R*-tree node capacity chosen by hand to balance tree depth and leaf size; the authors state this choice minimizes the number of tetrahedra checked but give no systematic study.
  • GPU_block_size = 8x8
    Thread block size found empirically to work best on GTX 10XX GPUs; the authors recommend re-tuning for other GPUs.
assumptions (6)
  • standard math Moller-Trumbore ray-triangle intersection correctly computes barycentric coordinates and the ray parameter t in exact arithmetic.
    The algorithm relies on this standard geometric test; in floating point it can leak, which the paper addresses.
  • domain assumption The volumetric mesh is convex and represented as a graph with complete, consistent neighbor information.
    Section II-B and II-D: propagation moves to the neighbor of the face at the largest t; this requires a closed, oriented surface and a convex domain.
  • ad hoc to paper Increasing the triangle face size by the safety parameter epsilon does not change the intersection parameter t in a way that affects the integral, nor does it introduce false positives for elements already known to be intersected.
    Section II-C, after Algorithm 1: this is asserted with a plausibility argument, not a proof or error bound.
  • domain assumption The R*-tree initialization returns the earliest boundary intersection for each ray.
    Section II-E: the depth-first search checks all tetrahedra in the leaf node and keeps the minimum t; this assumes the R*-tree bounding boxes are met correctly and the boundary is closed.
  • domain assumption Beer-Lambert line integral model in Eq. (2) and (3) is a valid linear model for X-ray CT.
    Section II-A: the system matrix entries are line intersection lengths; this is the standard CT forward model.
  • domain assumption OS-SART converges to a useful reconstruction for both voxel and tetrahedral bases with the chosen parameters.
    Section III-A: OS-SART is used for all reconstructions; convergence is not analyzed for the tetrahedral basis.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Numerically robust tetrahedron-based tomographic forward and backward projectors on parallel architectures." pith.science (2026). https://pith.science/paper/FTYPP2AL

@misc{pith2026190806909,
  author       = {Pith},
  title        = {Pith review of: Numerically robust tetrahedron-based tomographic forward and backward projectors on parallel architectures},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/FTYPP2AL}},
  note         = {Machine review of arXiv:1908.06909}
}
read the original abstract

X-ray tomographic reconstruction typically uses voxel basis functions to represent volumetric images. Due to the structure in voxel basis representations, efficient ray-tracing methods exist allowing fast, GPU accelerated implementations. Tetrahedral mesh basis functions are a valuable alternative to voxel based image representations as they provide flexible, inhomogeneous partitionings which can be used to provide reconstructions with reduced numbers of elements or with arbitrarily fine object surface representations. We thus present a robust parallelizable ray-tracing method for volumetric tetrahedral domains developed specifically for Computed Tomography image reconstruction. Tomographic image reconstruction requires algorithms that are robust to numerical errors in floating point arithmetic whilst typical data sizes encountered in tomography require the algorithm to be parallelisable in GPUs which leads to additional constraints on algorithm choices. Based on these considerations, this article presents numerical solutions to the design of efficient ray-tracing algorithms for the projection and backprojection operations. Initial reconstruction results using CAD data to define a triangulation of the domain demonstrate the advantages of our method and contrast tetrahedral mesh based reconstructions to voxel based methods.

Figures

Figures reproduced from arXiv: 1908.06909 by the authors.

Figure 1
Figure 1. Graph representation of a triangular mesh in 2D. Each element [PITH_FULL_IMAGE:figures/full_fig_p002_1.png] view at source ↗
Figure 2
Figure 2. The intersection algorithm computes t1 and t2 and the face labels denoting the faces that are intersected. Ray-triangle intersection algorithms have been widely stud￾ied, as the computer graphics field requires fast intersections for rendering with ray-tracing. The most common method is the Moller Trumbore[13] algorithm as it is one of the fastest ¨ method available. A faster algorithm is available by Baldwin and We… view at source ↗
Figure 4
Figure 4. Surface representation of the tetrahedra mesh used for validation of [PITH_FULL_IMAGE:figures/full_fig_p005_4.png] view at source ↗
Figures from the paper (5 more)
Figure 5
Figure 5. Figure 5: Simulated X-ray projections for 0 ◦ and 90◦ rotation positions. mesh, and a 56×56×56 mesh (the same amount of elements as the original tetrahedra mesh) spanning the same volumetric area is used. OS-SART[19] is here used for reconstruction. The algorithm is run with the…
Figure 6
Figure 6. Figure 6: Iterative reconstruction using the OS-SART algorithm with 50 itera [PITH_FULL_IMAGE:figures/full_fig_p006_6.png]
Figure 7
Figure 7. Figure 7: CAD model designed for CT image quality and feature detection [PITH_FULL_IMAGE:figures/full_fig_p007_7.png]
Figure 8
Figure 8. Figure 8: Iterative reconstruction using the OS-SART algorithm with 50 [PITH_FULL_IMAGE:figures/full_fig_p007_8.png]
Figure 10
Figure 10. Figure 10: Kernel initialization and ray-propagation times for a single [PITH_FULL_IMAGE:figures/full_fig_p008_10.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

22 extracted references · 22 canonical work pages

  1. [1]

    A fast forward projection using multithreads for multirays on GPUs in medical image reconstruction,

    Cheng-Ying Chou, Yi-Yen Chuo, Yukai Hung, and Weichung Wang, “A fast forward projection using multithreads for multirays on GPUs in medical image reconstruction,” Medical Physics, vol. 38, no. 7, pp. 4052–4065, 2011. 9

    work page 2011

  2. [2]

    Systematic performance optimiza- tion of cone-beam back-projection on the Kepler architecture,

    Timo Zinsser and Benjamin Keck, “Systematic performance optimiza- tion of cone-beam back-projection on the Kepler architecture,”

    work page

  3. [3]

    Tomographic image reconstruction based on a content-adaptive mesh model,

    Jovan G Brankov, Yongyi Yang, and Miles N Wernick, “Tomographic image reconstruction based on a content-adaptive mesh model,” IEEE Transactions on medical imaging , vol. 23, no. 2, pp. 202–212, 2004

    work page 2004

  4. [4]

    Tomo- graphic reconstruction using an adaptive tetrahedral mesh defined by a point cloud,

    Arkadiusz Sitek, Ronald H. Huesman, and Grant T. Gullberg, “Tomo- graphic reconstruction using an adaptive tetrahedral mesh defined by a point cloud,” IEEE Transactions on Medical Imaging , vol. 25, no. 9, pp. 1172–1179, sep 2006

    work page 2006

  5. [5]

    The sinogram polygonizer for reconstructing 3D shapes,

    Daiki Yamanaka, Yutaka Ohtake, and Hiromasa Suzuki, “The sinogram polygonizer for reconstructing 3D shapes,” IEEE Transactions on Visualization and Computer Graphics , 2013

    work page 2013

  6. [6]

    Michele Arcangelo Quinto, M´ethode de reconstruction adaptive en tomographie par rayons X: optimisation sur architectures parall `eles de type GPU, Grenoble, Apr 2013

    work page 2013

  7. [7]

    Tetrahedral volume reconstruction in X-ray tomography using GPU architecture,

    Michele Arcangelo Quinto, Dominique Houzet, and Fanny Buyens, “Tetrahedral volume reconstruction in X-ray tomography using GPU architecture,” 2013 Conference on Design and Architectures for Signal and Image Processing , pp. 334–339, 2013

    work page 2013

  8. [8]

    Practical cone-beam algorithm,

    L.A. Feldkamp, L.C. Davis, and J.W. Kress, “Practical cone-beam algorithm,” J. Opt. Soc. Am. A , vol. 1, no. 6, pp. 612–619, Jun 1984

    work page 1984

Show all 22 references
  1. [9]

    Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm,

    A.H. Andersen and A.C. Kak, “Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm,” Ultrasonic imaging, vol. 6, no. 1, pp. 81–94, 1984

    1984

  2. [10]

    Bjorck, Numerical Methods for Least Squares Problems , Society for Industrial and Applied Mathematics, 1996

    A. Bjorck, Numerical Methods for Least Squares Problems , Society for Industrial and Applied Mathematics, 1996

    1996

  3. [11]

    TIGRE: A MATLAB-GPU toolbox for CBCT image reconstruction,

    Ander Biguri, Steven Hancock, Manjit Dosanjh, and Manuchehr Soleimani, “TIGRE: A MATLAB-GPU toolbox for CBCT image reconstruction,” Biomedical Physics & Engineering Express , vol. 2, no. 5, pp. 055010, 2016

    2016

  4. [12]

    Fast and flexible x-ray tomography using the ASTRA toolbox,

    Wim van Aarle, Willem Jan Palenstijn, Jeroen Cant, Eline Janssens, Folkert Bleichrodt, Andrei Dabravolski, Jan De Beenhouwer, K Joost Batenburg, and Jan Sijbers, “Fast and flexible x-ray tomography using the ASTRA toolbox,” Optics express, vol. 24, no. 22, pp. 25129–25147, 2016

    2016

  5. [13]

    Fast, minimum storage ray/triangle intersection,

    Tomas M ¨oller and Ben Trumbore, “Fast, minimum storage ray/triangle intersection,” in ACM SIGGRAPH 2005 Courses . ACM, 2005, p. 7

    2005

  6. [14]

    Fast ray-triangle intersections by coordinate transformation,

    Doug Baldwin and Michael Weber, “Fast ray-triangle intersections by coordinate transformation,” Journal of Computer Graphics Techniques , vol. 5, no. 3, 2016

    2016

  7. [15]

    IEEE standard for binary floating-point arithmetic, Institute of Electrical and Electronics Engineers, New York, 1985, Note: Standard 754–1985

    1985

  8. [16]

    Watertight ray/triangle intersection,

    Sven Woop, Carsten Benthin, and Ingo Wald, “Watertight ray/triangle intersection,” Journal of Computer Graphics Techniques , vol. 2, no. 1, 2013

    2013

  9. [17]

    The r*-tree: an efficient and robust access method for points and rectangles,

    Norbert Beckmann, Hans-Peter Kriegel, Ralf Schneider, and Bernhard Seeger, “The r*-tree: an efficient and robust access method for points and rectangles,” in Acm Sigmod Record . Acm, 1990, vol. 19, pp. 322– 331

    1990

  10. [18]

    14, ACM, 1984

    Antonin Guttman, R-trees: A dynamic index structure for spatial searching, vol. 14, ACM, 1984

    1984

  11. [19]

    Block-iterative algorithms with diagonally scaled oblique projections for the linear feasibility problem,

    Y . Censor and T. Elfving, “Block-iterative algorithms with diagonally scaled oblique projections for the linear feasibility problem,” SIAM Journal on Matrix Analysis and Applications , vol. 24, no. 1, pp. 40–58, 2002

    2002

  12. [20]

    Thick- ness measurement of metal plate using CT projection images and nominal shape,

    Tasuku Ito, Yutaka Ohtake, Yukie Nagai, and Hiromasa Suzuki, “Thick- ness measurement of metal plate using CT projection images and nominal shape,” in 9th Conference on Industrial Computed Tomography. Feb. 2019, NDT

    2019

  13. [21]

    On packing r-trees,

    Ibrahim Kamel and Christos Faloutsos, “On packing r-trees,” University of Maryland

  14. [22]

    Initial result on the use of tetrahedra-based imaging for limited angle tomography,

    Ander Biguri, Hossein Towsyfyan, Richard Broadman, and Thomas Blumensath, “Initial result on the use of tetrahedra-based imaging for limited angle tomography,” in 9th Conference on Industrial Computed Tomography. Feb. 2019, NDT

    2019

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.