Hypergeometric Grid Visibility

High-throughput full-field visibility for 2D grids and 3D voxel volumes.

The project computes one scalar visibility value from a source to every cell of a transmission field. It provides a native C++20 2D solver, parallel OpenMP and optional CUDA 3D solvers, and NumPy reference implementations. It also includes exact canonical-path enumeration, comparisons with continuous visibility under grid refinement, and the three 3D comparison methods used in the paper.

These implementations accompany Hypergeometric Bridges for Multidimensional Grid Visibility.

The public interface covers field computation, rendering, bounded path enumeration, accuracy checks, and timing. See the quick start below and docs/reproduction.md for the paper experiments.

3D example

The yellow sphere marks the source, and the pink shapes are opaque obstacles. The left panel shows cells below a display threshold throughout the volume, coloured by their computed invisibility (1 - V). The right panel shows the corresponding field on the floor, rear boundary, and right boundary.

Install

Python 3.10 or newer is required.

git clone https://github.com/IbrahimSquared/hypergeometric-grid-visibility.git
cd hypergeometric-grid-visibility
python -m pip install -e '.[visualization]'

Build the native C++ 2D and OpenMP 3D solvers with CMake 3.24+, a C++20 compiler, and OpenMP:

cmake --preset cpu-release
cmake --build --preset cpu-release

For the optional CUDA solver, install the CUDA toolkit and run:

cmake --preset cuda-release
cmake --build --preset cuda-release

Python API

The transmission field is represented by an array indexed in (x, y) or (x, y, z) order. A transmission value of zero blocks visibility, one preserves it, and an intermediate value attenuates it proportionally.

import numpy as np
from hypergrid_visibility import compute_visibility, render_visibility, save_result

transmission = np.ones((384, 288), dtype=np.float32)
transmission[168:188, 32:232] = 0.0
transmission[168:188, 124:160] = 1.0

result = compute_visibility(
    transmission,
    source=(56, 140),
    method="two-parent",
    implementation="cpp",
)
save_result("visibility.npz", result)
render_visibility(result, "visibility.pdf")

The bundled quick start computes two distinct 2D fields: one from a generated 384-by-288 map and one by loading the viewable binary map examples/maps/visibility_polygon.png. Each run uses the native C++ solver by default, saves the numeric field, renders the result, and prints its native solve time:

python examples/quickstart.py --output-dir outputs --format png

Use --implementation numpy to run the same examples with the NumPy recurrence. elapsed_seconds reports in-process computation time for NumPy and native solve time for C++, OpenMP, and CUDA. Native timings exclude executable launch, binary-file exchange, buffer allocation, and rendering. CUDA solve time also excludes transfers between CPU and GPU memory; the paper reproduction runner reports transfer-inclusive CUDA time separately. The Python API defaults to NumPy when no implementation is named; the quick start explicitly selects cpp.

result.visibility is the numeric field. load_map reads NPY, NPZ, or a 2D image. render_visibility renders a 2D field directly and renders the three source-centred coordinate slices of a 3D field. The same compute call selects a native 3D implementation explicitly:

result = compute_visibility(
    volume,
    source=(20, 30, 40),
    method="four-parent",
    implementation="openmp",
    threads=8,
)

The public names match the paper:

Scientific method	Dimension	Implementations	Accepted transmission values
two-parent	2D	cpp, numpy	binary or nonbinary
four-parent	3D	numpy, openmp, cuda	binary or nonbinary
edge-gated-six-neighbour	3D	openmp	binary
edge-gated-tetrahedral	3D	openmp	binary
controlled-tetrahedral	3D	openmp	binary

Use cpp for fast 2D computation and NumPy for direct numerical comparison. In 3D, use OpenMP for rectangular CPU volumes, CUDA for cubic GPU volumes, and NumPy for direct numerical comparison. The three comparison methods reproduce the alternative rules studied in the paper; they are not alternative implementations of the four-parent method.

The OpenMP default is schedule="persistent-shell". The octant schedule is also available. It recorded the highest CPU throughput for centred N=1000 volumes but fell below the persistent-shell schedule for the near-corner cases. The default initial_visibility=1 is the paper's unit-source convention. To use the transmitting-source convention, pass initial_visibility=float(transmission[source]). The optional alpha parameter is the paper's constant per-layer decay factor; its default value of one applies no added decay.

Command line

Maps can be .npy, .npz, or grayscale 2D images. White image pixels are free and black pixels are blocked. Image files use their usual top-left origin and are converted to the Cartesian lower-left origin used by the API.

# Compute, store, and render a field
hgv compute map.npy --source 20,30 --method two-parent \
  --implementation cpp --output field.npz --render field.png

# Use the OpenMP solver on a 3D volume
hgv compute volume.npz --source 20,30,40 \
  --method four-parent --implementation openmp --threads 8

# Use the CUDA solver on a cubic 3D volume
hgv compute cubic-volume.npz --source 64,64,64 \
  --method four-parent --implementation cuda

# Enumerate and draw a bounded canonical-path family
hgv paths --source 0,0 --target 5,2 --output paths.json --render paths.pdf

# Compare a built-in physical scene with continuous visibility under refinement
hgv accuracy --scene G3 --resolutions 16,32,64 \
  --method four-parent --implementation openmp

# Time selected methods on one unchanged map and source
hgv benchmark examples/maps/visibility_polygon.png --source 299,279 \
  --methods two-parent --implementation cpp --warmups 5 --repeats 21

hgv benchmark volume.npz --source 20,30,40 \
  --methods four-parent,edge-gated-six-neighbour,edge-gated-tetrahedral \
  --implementation openmp --repeats 7

Path enumeration uses integer schedules and rational arithmetic. The --max-paths limit prevents accidental enumeration of combinatorially large families. Rendering supports PDF, PNG, JPG, and JPEG; stored numeric results use compressed NPZ.

Verification

python -m pip install -e '.[visualization,test]'
cmake --preset cpu-release
cmake --build --preset cpu-release
ctest --test-dir build/cpu --output-on-failure
pytest -q

The tests compare recurrence values with the finite-path average computed by direct canonical-path enumeration in 2D and 3D, check the native C++ 2D solver against NumPy, require bitwise equality across the native 3D CPU schedules, and check numerical parity between NumPy, OpenMP, and CUDA when those binaries are available.

Reported performance results

The paper compares the four CPU methods on an Intel Core i9-13980HX using eight pinned OpenMP threads. The table reports throughput on the G4 oblique-box scene in million returned visibility values per second.

Method	N=64	N=128	N=256	N=512	Relative runtime at N=512
Four-parent	1036	1490	647	543	1.00×
Controlled tetrahedral	723	887	394	411	1.32×
Edge-gated six-neighbour	398	724	434	236	2.30×
Edge-gated tetrahedral	395	375	184	193	2.81×

The four-parent implementation recorded the highest throughput at every tested G0 and G4 resolution. A second study kept the four-parent method fixed and changed only the OpenMP or CUDA schedule. These results are geometric means over 32 combinations of G0 or G4, four source positions, and four resolutions. Throughput is in billion returned visibility values per second.

Execution scope	Earlier schedule	Current schedule	Earlier throughput	Current throughput	Gain
OpenMP solve	Thread team per layer	Persistent shell	0.745	0.797	6.9%
CUDA solve	Face kernels	Fused layer	2.40	4.28	78.6%
CUDA including transfers	Face kernels	Fused layer	0.829	0.996	20.2%

At N=1000, octant OpenMP recorded the highest CPU throughput for centred volumes, whereas persistent-shell OpenMP was faster near a corner. Fused CUDA had the highest device-resident throughput and remained faster than both CPU schedules near a corner when field transfers were included. The measurements are specific to the reported hardware and workloads. Checksummed summaries and reproduction commands are in docs/reproduction.md.

Scope

The solvers evaluate the declared recurrence in floating-point arithmetic. The recurrence represents the corresponding finite-path average. Continuous visibility is a separate geometric comparison, and cell transmission is a scalar attenuation model rather than complete optics. See docs/methods.md for method roles and qualifications, and docs/provenance.md for source lineage and checksums.

License and citation

Apache-2.0. See CITATION.cff for citation metadata.