Conversation
Implements the geometric "Toblerone" voxelisation method (Kirk et al., IEEE TMI 2020) as a selectable alternative to the existing brute-force partial volume estimator, exposed through a new -algorithm option on mesh2voxel and made the default. The shared first stage (surface-voxel intersection) is retained, after which Toblerone classifies whole voxels by a deterministic ray-casting parity fill and estimates edge-voxel fractions by adaptive sub-voxel subdivision, using exact convex-polytope clipping for clean planar cuts and finer reduced-ray resampling elsewhere. This replaces the brute-force method's dense fixed point lattice and its heuristic interior/exterior classification, improving both speed and robustness while leaving the original algorithm available and unchanged. Session prompts: 1. > File "toblerone.txt" provides the content of a scientific journal article that describes a method for computing, from a closed surface mesh in 3D, partial volume fractions for a 3D voxel grid. MRtrix3 already has functionality for this task, implemented in cpp/core/surface/algo/mesh2image.* and interfaced in command cpp/cmd/mesh2voxel.cpp. However this implementation is both slow (as it involves generating and testing a dense 3D point lattice for each intersected voxel) and bug-prone (relies on heuristic classification of non-intersected voxels as inside vs. outside surface). Attempt to draft a plan for implementation of the "Toblerone" algorithm alongside the existing brute-force algorithm. Comment on whether the failure to adequately propagate mathematical symbols into this text file, or the absence of demonstrative figure images, is prohibitive. Additional web searching is permitted. Generated-by: Claude Opus 4.8 <noreply@anthropic.com>
- Add manuscript describing method to mesh2voxel references. - Update test data and tests.
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
| + Argument ("template", "the template image").type_image_in() | ||
| + Argument ("output", "the output image").type_image_out(); | ||
|
|
||
| OPTIONS |
There was a problem hiding this comment.
warning: no header providing "MR::App::OPTIONS" is directly included [misc-include-cleaner]
OPTIONS
^|
|
||
| // Select the partial volume estimation algorithm | ||
| Surface::Algo::Mesh2ImageOptions options; | ||
| auto opt = get_options("algorithm"); |
There was a problem hiding this comment.
warning: no header providing "MR::App::get_options" is directly included [misc-include-cleaner]
auto opt = get_options("algorithm");
^| options.method = MR::Enum::from_name<Surface::Algo::Mesh2ImageMethod>(opt[0][0]); | ||
| opt = get_options("subvoxel"); | ||
| if (!opt.empty()) | ||
| options.subvoxel_mm = static_cast<default_type>(opt[0][0]); |
There was a problem hiding this comment.
warning: no header providing "MR::default_type" is directly included [misc-include-cleaner]
cpp/cmd/mesh2voxel.cpp:26:
+ #include "types.h"|
|
||
| void mesh2image(const Mesh &mesh_realspace, Image<float> &image) { | ||
| // For every edge voxel, stores those polygons that may intersect the voxel | ||
| using Vox2Poly = std::map<Vox, std::vector<size_t>>; |
There was a problem hiding this comment.
warning: no header providing "std::vector" is directly included [misc-include-cleaner]
cpp/core/surface/algo/mesh2image.cpp:24:
+ #include <vector>| namespace { | ||
|
|
||
| // Encoding of the per-voxel interior/exterior classification used by the Toblerone algorithm | ||
| enum vox_class_t : uint8_t { VC_OUTSIDE = 0, VC_INSIDE = 1, VC_ON_MESH = 2 }; |
There was a problem hiding this comment.
warning: no header providing "uint8_t" is directly included [misc-include-cleaner]
cpp/core/surface/algo/mesh2image.cpp:20:
- #include <map>
+ #include <cstdint>
+ #include <map>| } | ||
| for (size_t axis = 0; axis != 3; ++axis) { | ||
| if (std::min({v0[axis], v1[axis], v2[axis]}) > box_half[axis] || | ||
| std::max({v0[axis], v1[axis], v2[axis]}) < -box_half[axis]) |
There was a problem hiding this comment.
warning: narrowing conversion from 'size_t' (aka 'unsigned long') to signed type 'Index' (aka 'long') is implementation-defined [bugprone-narrowing-conversions]
std::max({v0[axis], v1[axis], v2[axis]}) < -box_half[axis])
^| } | ||
| for (size_t axis = 0; axis != 3; ++axis) { | ||
| if (std::min({v0[axis], v1[axis], v2[axis]}) > box_half[axis] || | ||
| std::max({v0[axis], v1[axis], v2[axis]}) < -box_half[axis]) |
There was a problem hiding this comment.
warning: narrowing conversion from 'size_t' (aka 'unsigned long') to signed type 'Index' (aka 'long') is implementation-defined [bugprone-narrowing-conversions]
std::max({v0[axis], v1[axis], v2[axis]}) < -box_half[axis])
^| return false; | ||
| } | ||
| if (separating(e0.cross(e1))) | ||
| return false; |
There was a problem hiding this comment.
warning: redundant boolean literal in conditional return statement [readability-simplify-boolean-expr]
cpp/core/surface/algo/mesh2image.cpp:187:
- if (separating(e0.cross(e1)))
- return false;
- return true;
+ return !separating(e0.cross(e1));|
|
||
| // Whether any polygon of the local patch overlaps the axis-aligned box (centre, half-extents) | ||
| bool patch_overlaps_box(const Vertex ¢re, const Vertex &half, const std::vector<size_t> &polys, const Mesh &mesh) { | ||
| for (const size_t poly_index : polys) { |
There was a problem hiding this comment.
warning: replace loop by 'std::any_of()' [readability-use-anyofallof]
for (const size_t poly_index : polys) {
^| // Volume of the part of tetrahedron `v` interior to the plane through q with outward normal n, | ||
| // i.e. the side where (x - q).n <= 0, by enumeration of the marching-tetrahedra cases. | ||
| default_type tetra_clip_interior_volume(const std::array<Vertex, 4> &v, const Vertex &q, const Vertex &n) { | ||
| std::array<default_type, 4> f; |
There was a problem hiding this comment.
warning: uninitialized record type: 'f' [cppcoreguidelines-pro-type-member-init]
| std::array<default_type, 4> f; | |
| std::array<default_type, 4> f{}; |
1 participant
Given I recently attacked the long-standing
mesh2voxelbug (#3391), #2732 came to mind. I had some recollection that it wasn't just the computation of partial volume for intersected voxels that was different, there was also some difference in how it dealt with labelling voxels inside vs. outside that surface, which potentially would have obviated thatmesh2voxelissue. So I thought I'd take the chance to throw an agent at it, givenmesh2voxelis often a performance bottleneck, particularly in5ttgen hsvs.~ 7x performance improvement, <= 0.025 PVE difference to my brute force approach (with no real guarantee that mine is in fact more accurate).
Reconsider
-subvoxeloption:In addition to defining the default value as a constexpr and printing it into the command help page rather than duplicating, it might be preferable to change how that default is determined, given the software may operate on non-human data. Perhaps it should instead be a voxel size fraction?
@Lestropie Full manual code review.