CASE STUDY · PARALLEL COMPUTING
Parallel 3D Voxel Pattern Matching v2
High-performance 3D CT voxel pattern matching with optimized sequential, POSIX Threads and MPI implementations.
A high-performance 3D CT voxel pattern-matching system designed to identify all ordered pairs of identical 4 × 4 × 4 voxel blocks within the first 64 slices of a 1024 × 1024 × 314 CT dataset. The solution replaces an impractical quadratic comparison strategy with compact 64-bit bitmask encoding, radix sorting, and linear run aggregation. It includes optimized sequential, POSIX Threads, and MPI implementations, all verified against the same result of 970,777,711,592 matching ordered pairs.
- ROLE
- Lead Developer · Algorithm Design, Optimization, Parallelization, Benchmarking, and Technical Documentation
- YEAR
- 2026
TECH_STACK
Technologies
MEDIA
Project media
CONTEXT
Technical context
The project was developed around a raw CT volume containing 1024 × 1024 × 314 voxels, with each voxel stored as an unsigned 8-bit value. Only the first 64 slices are processed, producing a 1024 × 1024 × 64 working volume. Each voxel is thresholded at a value of 25: values less than or equal to the threshold are classified as inactive, while higher values are classified as active. The volume is then partitioned into non-overlapping 4 × 4 × 4 blocks, resulting in 1,048,576 individual voxel patterns. Because each block contains exactly 64 binary voxel states, it can be represented as a single `uint64_t` bitmask. The task is to count every ordered pair of identical patterns, meaning that A → B and B → A are treated as separate matches. This requirement creates a large computational workload and makes the project well suited for algorithmic optimization, shared-memory parallelism, and distributed-memory processing.
PROBLEM
Problem
The original implementation compared every generated bitmask with every other bitmask, producing approximately 1.1 trillion comparisons for 1,048,576 patterns. Although functionally correct, this quadratic O(n²) approach caused excessive execution time, high memory usage, and poor scalability. The main challenge was to preserve the exact matching semantics while replacing brute-force pairwise comparison with a substantially more efficient algorithm. The optimized solution also needed to support three execution models—sequential, POSIX Threads, and MPI—while ensuring that every version returned the same count of 970,777,711,592 ordered matching pairs and measured computation time only after the input data had been loaded from disk.
OVERVIEW
Project overview
The project processes the first 64 slices of a 1024 × 1024 × 314 raw CT dataset and converts the selected volume into 1,048,576 non-overlapping 4 × 4 × 4 voxel patterns. Each block is thresholded at a value of 25 and encoded as a single 64-bit bitmask.
The original solution relied on quadratic pairwise comparison, requiring approximately 1.1 trillion comparisons. The redesigned architecture uses an 8-pass LSD radix sort to group identical masks and calculates the number of ordered matching pairs through linear run aggregation.
Three implementations were developed using the same core algorithm: an optimized sequential baseline, a shared-memory POSIX Threads version, and a distributed-memory MPI version. Each implementation preserves identical matching semantics and produces the verified result of 970,777,711,592 ordered pairs.
The POSIX Threads implementation parallelizes bitmask generation and stable radix sorting, while the MPI implementation distributes generation and sorting across processes before compressing and redistributing local mask counts for final aggregation.
Every benchmark configuration was executed 20 times after the CT data had been loaded from disk. The best result was achieved by the 16-thread POSIX Threads implementation with a median computation time of 0.012105 seconds and a 3.95× speedup over the optimized sequential baseline.
APPROACH
Technical approach
- Loaded the first 64 slices of the 1024 × 1024 × 314 CT dataset before starting computation timing
- Applied a threshold of 25 to classify each voxel as inactive or active
- Partitioned the working volume into 1,048,576 non-overlapping 4 × 4 × 4 voxel blocks
- Encoded every 64-voxel block directly as a compact `uint64_t` bitmask
- Replaced the original quadratic pairwise comparison with an 8-pass LSD radix sort
- Counted contiguous runs of identical masks and calculated ordered pairs using `k × (k − 1)`
- Established an optimized sequential implementation as the reference baseline
- Parallelized bitmask generation and stable radix sorting with POSIX Threads
- Distributed mask generation, local sorting, run compression, ownership redistribution, and final aggregation with MPI
- Validated every implementation and benchmark configuration against the same result of 970,777,711,592 ordered matching pairs
- Repeated each production benchmark 20 times and evaluated performance using median, mean, minimum, and maximum execution times
FOCUS
Technical focus
- Quadratic-to-linear algorithm redesign
- 3D CT voxel thresholding and block extraction
- Compact 64-bit voxel-pattern encoding
- Eight-pass LSD radix sorting
- Ordered matching-pair aggregation
- Sequential performance optimization
- POSIX Threads shared-memory parallelization
- MPI distributed-memory parallelization
- Compressed mask-count redistribution
- Synchronization and communication overhead analysis
- Memory-layout and allocation optimization
- Repeatable multi-configuration benchmarking
- Cross-implementation correctness validation
ARCHITECTURE
Architecture notes
- Contiguous 64 MiB buffer for the processed CT volume
- Direct 4 × 4 × 4 block-to-uint64_t encoding without intermediate pattern matrices
- Shared core matching logic across sequential, POSIX Threads, and MPI implementations
- Eight stable byte-wise radix passes over 64-bit masks
- Linear run aggregation using k × (k − 1) for ordered pair counting
- Thread-partitioned bitmask generation in the POSIX Threads implementation
- Barrier-synchronized shared-memory radix sorting without result mutexes
- Rank-local mask generation and radix sorting in the MPI implementation
- Local run compression into compact mask-and-count records
- Hash-based ownership assignment for identical masks across MPI ranks
- MPI_Alltoallv redistribution of compressed records to owner ranks
- Distributed count aggregation followed by MPI_Reduce for the final result
- Computation timing isolated from input file loading
- Identical correctness target across every execution model and benchmark configuration
HIGHLIGHTS
Engineering highlights
- Reduced approximately 1.1 trillion pairwise comparisons to a radix-based processing pipeline
- Processed 1,048,576 non-overlapping 4 × 4 × 4 voxel patterns
- Encoded every voxel pattern into a single 64-bit value
- Verified 970,777,711,592 ordered matching pairs across all implementations
- Achieved a 0.047756-second median with the optimized sequential implementation
- Reached the best overall performance at 16 POSIX threads with a 0.012105-second median
- Delivered a 3.95× speedup with shared-memory parallelization
- Reached the best MPI performance at 8 processes with a 0.029341-second median
- Delivered a 1.63× speedup with distributed-memory processing
- Validated every thread and process configuration through 20 benchmark runs
- Reduced memory usage through contiguous storage and compact bitmask representation
- Documented the complete algorithm, architecture, validation process, and benchmark results in English and Slovak
OUTCOME
Outcome
The optimized solution replaced approximately 1.1 trillion pairwise comparisons with a radix-based processing pipeline and reduced execution time from tens of seconds or minutes to milliseconds. All sequential, POSIX Threads, and MPI implementations produced the same verified result of 970,777,711,592 ordered matching pairs across every production benchmark run. The sequential implementation achieved a median computation time of 0.047756 seconds, while the POSIX Threads version reached the best overall performance with 16 threads and a median of 0.012105 seconds, delivering a 3.95× speedup. The MPI implementation performed best with 8 processes at a median of 0.029341 seconds and a 1.63× speedup. Higher thread and process counts provided no further improvement because synchronization, memory bandwidth, and communication overhead began to dominate the short computation.