img_phy_sim.ism

   1"""
   2**Image Source Method (ISM) Ray Propagation and Visibility Maps**
   3
   4This module implements a fast 2D Image Source Method (ISM) solver to simulate
   5specular reflections and line-of-sight propagation inside environments described
   6by images (raster maps) that contain wall/obstacle pixels.
   7
   8The scene is interpreted as a 2D grid where specific pixel values represent
   9walls. From this raster representation, geometric wall boundary segments are
  10extracted and used to construct image sources for reflection sequences up to a
  11given order. For each receiver position on a grid, candidate reflection paths
  12are constructed, validated with a raster-based visibility test, and accumulated
  13into a path-count map.
  14
  15Core pipeline:
  161. Build a wall mask from an input image
  172. Extract wall boundary segments (OpenCV contours -> polyline edges)
  183. Build a binary occlusion map for raster visibility tests
  194. Enumerate reflection sequences and precompute image source positions
  205. Evaluate receiver points on a grid:
  21   - Construct candidate reflection paths (segment intersection + reflection)
  22   - Validate segment visibility against the occlusion raster
  23   - Accumulate valid paths into an output map
  24
  25As ASCII model:
  26```text
  27                 ┌──────────────────────────────┐
  28                 │        Input: Scene          │
  29                 │    (image, walls, source)    │
  30                 └──────────────┬───────────────┘
  31                                │
  32                                v
  33                 ┌──────────────────────────────┐
  34                 │   Raster → Wall Geometry     │
  35                 │   Extract wall segments      │
  36                 └──────────────┬───────────────┘
  37                                │
  38                                v
  39                 ┌──────────────────────────────┐
  40                 │     Build Occlusion Map      │
  41                 │   (for visibility checks)    │
  42                 └──────────────┬───────────────┘
  43                                │
  44                                v
  45                 ┌──────────────────────────────┐
  46                 │  Precompute Image Sources    │
  47                 │  (all reflection sequences)  │
  48                 └──────────────┬───────────────┘
  49                                │
  50                                v
  51                 ┌──────────────────────────────┐
  52                 │     Define Receiver Grid     │
  53                 └──────────────┬───────────────┘
  54                                │
  55                                v
  56                 ┌──────────────────────────────────────────┐
  57                 │ For each Receiver position R on the grid │
  58                 └──────────────┬───────────────────────────┘
  59                                │
  60                                v
  61                 ┌──────────────────────────────────────────┐
  62                 │ For each Image Source (reflection seq)   │
  63                 └──────────────┬───────────────────────────┘
  64                                │
  65                                v
  66                 ┌──────────────────────────────┐
  67                 │  Construct reflection path   │
  68                 │   (geometry / intersections) │
  69                 └──────────────┬───────────────┘
  70                                │
  71                                v
  72                 ┌──────────────────────────────┐
  73                 │     Check path visibility    │
  74                 │    (raster occlusion test)   │
  75                 └──────────────┬───────────────┘
  76                                │
  77                                v
  78                 ┌──────────────────────────────┐
  79                 │  Accumulate contribution     │
  80                 └──────────────┬───────────────┘
  81                                │
  82                                v
  83                 ┌──────────────────────────────┐
  84                 │          Output Map          │
  85                 └──────────────────────────────┘
  86```
  87
  88<br><br>
  89
  90Numba/parallelization notes:<br>
  91The hot loop is executed inside Numba kernels.
  92
  93- `eval_receivers_to_map_kernel`:<br>
  94  Computes and writes directly into `output_map` inside Numba. This is safe when
  95  receiver coordinates are unique (e.g., a regular grid where each receiver maps to
  96  a distinct pixel).
  97
  98Reflection sequences are packed into flat arrays (`seq_data`, `seq_off`,
  99``seq_len``) to avoid Python containers and to keep Numba in nopython mode.
 100
 101Example:
 102```python
 103reflection_map = ips.ism.compute_reflection_map(
 104    source_rel=(0.5, 0.5),
 105    img=input_,
 106    wall_values=[0],   
 107    wall_thickness=1,
 108    max_order=1,
 109    step_px=1
 110)
 111```
 112
 113Dependencies:
 114- numpy
 115- OpenCV (cv2)
 116- numba
 117
 118
 119Public API:
 120- compute_reflection_map(...)
 121- build_wall_mask(...)
 122- get_wall_segments_from_mask(...)
 123- build_occlusion_from_wallmask(...)
 124- precompute_image_sources(...)
 125- enumerate_wall_sequences_indices(...)
 126
 127Main internal kernels/utilities:
 128- reflect_point_across_infinite_line(...)
 129- reflect_point_across_infinite_line_numba(...)
 130- seg_seg_intersection_xy(...)
 131- is_visible_raster(...)
 132- build_path_for_sequence_packed(...)
 133- check_path_visibility_raster_arr(...)
 134- pack_precomputed(...)
 135- build_receivers_grid(...)
 136- eval_receivers_to_map_kernel(...)
 137"""
 138
 139
 140
 141# ---------------
 142# >>> Imports <<<
 143# ---------------
 144from __future__ import annotations
 145
 146import math
 147from dataclasses import dataclass
 148from typing import List, Tuple, Optional
 149
 150import numpy as np
 151import cv2
 152
 153# optimization
 154import numba
 155
 156
 157
 158# --------------
 159# >>> Helper <<<
 160# --------------
 161
 162def reflect_point_across_infinite_line(P: Tuple[float, float], A: Tuple[float, float], B: Tuple[float, float]) -> Tuple[float, float]:
 163    """
 164    Reflect a 2D point across an infinite line.
 165
 166    Computes the mirror image of point P with respect to the infinite line
 167    passing through points A and B.
 168
 169                  P
 170                  *
 171                   \
 172                    \
 173    -----------------*----------------   Wand (A-B)
 174                    proj
 175                      \
 176                       \
 177                        *
 178                        P'
 179
 180    Parameters:
 181    - P (Tuple[float, float]): <br>
 182        The point to reflect (x, y).
 183    - A (Tuple[float, float]): <br>
 184        First point defining the line (x, y).
 185    - B (Tuple[float, float]): <br>
 186        Second point defining the line (x, y).
 187
 188    Returns:
 189    - Tuple[float, float]: <br>
 190        The reflected point (x, y) as floats.
 191    """
 192    # var setup
 193    #    -> A/B are both points on the wall
 194    #    -> P is the point we want to mirror (not receiver, source ,...)
 195    x1, y1 = A
 196    x2, y2 = B
 197    px, py = P
 198
 199    # get direction of the wall
 200    #    -> from point to the wall
 201    ABx = x2 - x1
 202    ABy = y2 - y1
 203    # norm this wall direction
 204    n = math.hypot(ABx, ABy) + 1e-12
 205    ux, uy = ABx / n, ABy / n
 206
 207    # direction/vector from the wall to the point
 208    APx = px - x1
 209    APy = py - y1
 210
 211    # projection of point to wall direction to the wall direction
 212    #   -> finding sweetspot where we are directly under P on wall line/direction
 213    t = APx * ux + APy * uy
 214    # point on wall direction/line which is directly under P
 215    projx = x1 + t * ux
 216    projy = y1 + t * uy
 217
 218    # calc reflection/mirror
 219    rx = projx + (projx - px)
 220    ry = projy + (projy - py)
 221
 222    return (float(rx), float(ry))
 223
 224
 225
 226@numba.njit(cache=True, fastmath=True)
 227def reflect_point_across_infinite_line_numba(px, py, ax, ay, bx, by):
 228    """
 229    Reflect a 2D point across an infinite line defined by two points.
 230
 231    Computes the mirror image of point P = (px, py) with respect to the
 232    infinite line passing through A = (ax, ay) and B = (bx, by).
 233    The function first projects P orthogonally onto the line AB to obtain
 234    the foot point (projection), and then reflects P across this point
 235    using the relation:
 236
 237        P' = 2 * proj - P
 238
 239    This implementation is written in a Numba-friendly style and avoids
 240    expensive operations such as `math.hypot`, using only basic arithmetic
 241    and a square root for normalization.
 242
 243    Parameters:
 244    - px (float): <br>
 245        x-coordinate of the point to reflect.
 246    - py (float): <br>
 247        y-coordinate of the point to reflect.
 248    - ax (float): <br>
 249        x-coordinate of the first point defining the line.
 250    - ay (float): <br>
 251        y-coordinate of the first point defining the line.
 252    - bx (float): <br>
 253        x-coordinate of the second point defining the line.
 254    - by (float): <br>
 255        y-coordinate of the second point defining the line.
 256
 257    Returns:
 258    - Tuple[float, float]: <br>
 259        The reflected point (rx, ry) across the infinite line AB.
 260    """
 261    # direction vector of the wall segment A -> B
 262    abx = bx - ax
 263    aby = by - ay
 264
 265    # normalize to unit direction vector
 266    # small epsilon avoids division by zero for degenerate segments
 267    n = math.sqrt(abx * abx + aby * aby) + 1e-12
 268    ux = abx / n
 269    uy = aby / n
 270
 271    # vector from A to the point P
 272    apx = px - ax
 273    apy = py - ay
 274
 275    # projection of AP onto the wall direction
 276    t = apx * ux + apy * uy
 277    projx = ax + t * ux
 278    projy = ay + t * uy
 279
 280    # reflect point across the infinite line defined by A-B
 281    # P' = 2 * projection - P
 282    rx = 2.0 * projx - px
 283    ry = 2.0 * projy - py
 284
 285    return rx, ry
 286
 287
 288
 289def reflection_map_to_img(reflection_map):
 290    """
 291    Convert an reflection map to a uint8 visualization image.
 292
 293    Normalizes the input reflection_map to [0, 255] by dividing by its maximum
 294    value (with an epsilon for numerical stability), and converts the result
 295    to uint8.
 296
 297    Parameters:
 298    - reflection_map (np.ndarray): <br>
 299        A numeric array representing reflection values.
 300
 301    Returns:
 302    - np.ndarray: <br>
 303        A uint8 image array with values in [0, 255].
 304    """
 305    vis = reflection_map.copy()
 306    vis = vis / (vis.max() + 1e-9)
 307    return (vis * 255).astype(np.float64)  # .astype(np.uint8)
 308
 309
 310
 311# -----------------------------------------
 312# >>> Segment representation & geometry <<<
 313# -----------------------------------------
 314
 315@numba.njit(cache=True, fastmath=True)
 316def seg_seg_intersection_xy(x1, y1, x2, y2, x3, y3, x4, y4, eps=1e-9):
 317    """
 318    Compute the intersection point of two 2D line segments.
 319
 320    Determines whether the segment P1→P2 intersects with the segment Q1→Q2.
 321    The computation is performed using a parametric form and 2D cross products.
 322    If a unique intersection point exists within both segment bounds
 323    (including small numerical tolerances), the intersection coordinates are
 324    returned together with a success flag.
 325
 326    Parallel or colinear segments are treated as non-intersecting for the
 327    purposes of the Image Source Model, since such cases are ambiguous for
 328    specular reflection path construction.
 329
 330    Parameters:
 331    - x1 (float): <br>
 332        x-coordinate of the first endpoint of segment P.
 333    - y1 (float): <br>
 334        y-coordinate of the first endpoint of segment P.
 335    - x2 (float): <br>
 336        x-coordinate of the second endpoint of segment P.
 337    - y2 (float): <br>
 338        y-coordinate of the second endpoint of segment P.
 339    - x3 (float): <br>
 340        x-coordinate of the first endpoint of segment Q.
 341    - y3 (float): <br>
 342        y-coordinate of the first endpoint of segment Q.
 343    - x4 (float): <br>
 344        x-coordinate of the second endpoint of segment Q.
 345    - y4 (float): <br>
 346        y-coordinate of the second endpoint of segment Q.
 347    - eps (float): <br>
 348        Numerical tolerance used to handle floating point inaccuracies
 349        when testing for parallelism and segment bounds.
 350
 351    Returns:
 352    - Tuple[float, float, bool]: <br>
 353        (ix, iy, ok) where (ix, iy) is the intersection point if one exists,
 354        and `ok` indicates whether a valid intersection was found.
 355    """
 356    # r = segment p1 -> p2
 357    rx = x2 - x1
 358    ry = y2 - y1
 359
 360    # s = segment q1 -> q2
 361    sx = x4 - x3
 362    sy = y4 - y3
 363
 364    # cross product r x s used to detect parallelism
 365    rxs = rx * sy - ry * sx
 366
 367    # vector from p1 to q1
 368    qpx = x3 - x1
 369    qpy = y3 - y1
 370
 371    # cross product (q1 - p1) x r
 372    qpxr = qpx * ry - qpy * rx
 373
 374    # if r x s is near zero the segments are parallel or colinear
 375    # for the reflection model this is treated as an invalid intersection
 376    if abs(rxs) <= eps:
 377        # parallel or colinear -> treat as invalid for your ISM
 378        return 0.0, 0.0, False
 379
 380    # solve intersection parameters for the parametric segment equations
 381    t = (qpx * sy - qpy * sx) / rxs
 382    u = (qpx * ry - qpy * rx) / rxs
 383
 384    # check whether intersection lies within both segments (with epsilon tolerance)
 385    if (-eps <= t <= 1.0 + eps) and (-eps <= u <= 1.0 + eps):
 386        ix = x1 + t * rx
 387        iy = y1 + t * ry
 388        return ix, iy, True
 389
 390    # no valid intersection on the finite segments
 391    return 0.0, 0.0, False
 392
 393
 394
 395# -------------------------------
 396# >>> Raster-based visibility <<<
 397# -------------------------------
 398
 399@numba.njit(cache=True, fastmath=True)
 400def is_visible_raster(p1x:float, p1y:float, p2x:float, p2y:float, occ:np.ndarray, ignore_ends:int=1) -> bool:
 401    """
 402    Check line-of-sight visibility between two points using a raster occlusion map.
 403
 404    This function traces a discrete line between two pixel positions using
 405    Bresenham's line algorithm and tests whether any raster cell along this
 406    line is marked as occluded in the provided occlusion map `occ`.
 407    The algorithm operates entirely in integer pixel space and is written in
 408    a Numba-friendly style without temporary allocations.
 409
 410    To allow valid reflection paths that touch wall endpoints, a configurable
 411    number of pixels at the start and end of the line can be ignored via
 412    `ignore_ends`. This prevents false occlusion detections when rays start
 413    or end directly on wall pixels.
 414
 415    The Bresenham traversal is executed twice:
 416    1. First pass: counts the number of raster points along the line to
 417    determine the valid range after removing the ignored endpoints.
 418    2. Second pass: checks only the relevant interior pixels for occlusion.
 419
 420    Taken from https://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm -> last code shown on the website.
 421
 422    Parameters:
 423    - p1x (float): <br>
 424        x-coordinate of the start point in pixel space.
 425    - p1y (float): <br>
 426        y-coordinate of the start point in pixel space.
 427    - p2x (float): <br>
 428        x-coordinate of the end point in pixel space.
 429    - p2y (float): <br>
 430        y-coordinate of the end point in pixel space.
 431    - occ (np.ndarray): <br>
 432        Occlusion map where nonzero values represent blocked pixels.
 433        Expected indexing is `occ[y, x]`.
 434    - ignore_ends (int): <br>
 435        Number of raster pixels to ignore at both ends of the line.
 436
 437    Returns:
 438    - bool: <br>
 439        True if the path between the two points is free of occlusions,
 440        otherwise False.
 441    """
 442    H = occ.shape[0]
 443    W = occ.shape[1]
 444
 445    # round always up -> so + 0.5
 446    x0 = int(p1x + 0.5)
 447    y0 = int(p1y + 0.5)
 448    x1 = int(p2x + 0.5)
 449    y1 = int(p2y + 0.5)
 450
 451    # clamp endpoints -> should be inside
 452    if x0 < 0: x0 = 0
 453    if x0 >= W: x0 = W - 1
 454    if y0 < 0: y0 = 0
 455    if y0 >= H: y0 = H - 1
 456    if x1 < 0: x1 = 0
 457    if x1 >= W: x1 = W - 1
 458    if y1 < 0: y1 = 0
 459    if y1 >= H: y1 = H - 1
 460
 461    # preperation for breseham algorithm
 462    # err -> accumulator to decide to go in x or y direction
 463    dx = abs(x1 - x0)
 464    dy = abs(y1 - y0)
 465
 466    sx = 1 if x0 < x1 else -1
 467    sy = 1 if y0 < y1 else -1
 468
 469    err = dx - dy
 470
 471    # 1. calc length (amount of grid points -> for ignore ends
 472    # => count grid points between the 2 points
 473    # => later we can then say ignore the first and last N pixels
 474    tx, ty = x0, y0
 475    terr = err
 476    npts = 1
 477    while not (tx == x1 and ty == y1):
 478        # bresenham algorithm to move in pixel-space
 479        e2 = terr * 2
 480        if e2 > -dy:
 481            terr -= dy
 482            tx += sx
 483        if e2 < dx:
 484            terr += dx
 485            ty += sy
 486        # increase grid points
 487        npts += 1
 488
 489    # not enough grid points -> visible!
 490    # only start and end point, there cant be something in between
 491    if npts <= 2:
 492        return True
 493
 494    # if the ignore parts is cut away
 495    # is there even something left to check?
 496    start_i = ignore_ends
 497    end_i = npts - ignore_ends
 498    if start_i >= end_i:
 499        return True
 500    
 501    # 2. run again -> but this time check occlusion
 502    tx, ty = x0, y0
 503    terr = err
 504    i = 0
 505    while True:
 506        if i >= start_i and i < end_i:
 507            # wall between detected?
 508            if occ[ty, tx] != 0:
 509                return False
 510
 511        # end of line reached?
 512        if tx == x1 and ty == y1:
 513            break
 514
 515        # breseham -> continue line
 516        e2 = terr * 2
 517        if e2 > -dy:
 518            terr -= dy
 519            tx += sx
 520        if e2 < dx:
 521            terr += dx
 522            ty += sy
 523
 524        i += 1
 525
 526    return True
 527
 528
 529# -------------------------------------
 530# >>> Wall extraction from an image <<<
 531# -------------------------------------
 532
 533@dataclass(frozen=True)
 534class Segment:
 535    """
 536    Represent a 2D line segment.
 537
 538    A lightweight immutable segment representation used for wall geometry
 539    and intersection tests.
 540
 541    Attributes:
 542    - ax (float): x-coordinate of the first endpoint.
 543    - ay (float): y-coordinate of the first endpoint.
 544    - bx (float): x-coordinate of the second endpoint.
 545    - by (float): y-coordinate of the second endpoint.
 546
 547    Properties:
 548    - A (Tuple[float, float]): <br>
 549        First endpoint (ax, ay).
 550    - B (Tuple[float, float]): <br>
 551        Second endpoint (bx, by).
 552    """
 553    ax: float
 554    ay: float
 555    bx: float
 556    by: float
 557
 558    @property
 559    def A(self): return (self.ax, self.ay)
 560
 561    @property
 562    def B(self): return (self.bx, self.by)
 563
 564
 565
 566@numba.njit(cache=True, fastmath=True)
 567def _max_u8(img_u8):
 568    """
 569    Compute the maximum uint8 value in a 2D image using explicit loops
 570
 571    This helper avoids NumPy reductions so it can be used inside
 572    Numba-compatible code paths
 573
 574    Parameters:
 575    - img_u8 (np.ndarray): <br>
 576        2D uint8 image
 577
 578    Returns:
 579    - int: <br>
 580        Maximum pixel value found in the image
 581    """
 582    m = 0
 583    H, W = img_u8.shape[0], img_u8.shape[1]
 584    for y in range(H):
 585        for x in range(W):
 586            v = int(img_u8[y, x])
 587            if v > m:
 588                m = v
 589    return m
 590
 591
 592
 593@numba.njit(cache=True, fastmath=True)
 594def build_wall_mask_numba_no_values(img_u8):
 595    """
 596    Build a 0/255 wall mask from a uint8 image when no explicit wall values are given.
 597
 598    This function assumes that the input image is already mask-like. If the
 599    maximum pixel value is small (e.g. < 64), the image is treated as a
 600    binary or low-range mask and scaled to the range {0, 255}. Otherwise,
 601    the values are copied directly, assuming the image already represents
 602    a proper wall mask.
 603
 604    The implementation is written in a Numba-friendly style and operates
 605    purely with explicit loops for maximum compatibility and performance.
 606
 607    Parameters:
 608    - img_u8 (np.ndarray): <br>
 609        2D uint8 image interpreted as a mask-like input.
 610
 611    Returns:
 612    - np.ndarray: <br>
 613        A uint8 wall mask with values 0 (free space) and 255 (wall).
 614    """
 615    # img_u8: uint8 2D
 616    H, W = img_u8.shape[0], img_u8.shape[1]
 617
 618    # determine maximum value in the image to infer its encoding
 619    m = _max_u8(img_u8)
 620
 621    out = np.empty((H, W), dtype=np.uint8)
 622
 623    if m < 64:
 624        # image likely represents a 0/1 style mask
 625        # scale any non-zero value up to 255 for consistency
 626        for y in range(H):
 627            for x in range(W):
 628                out[y, x] = 255 if img_u8[y, x] != 0 else 0
 629    else:
 630        # image already looks like a proper uint8 mask or label image
 631        # simply copy values over
 632        for y in range(H):
 633            for x in range(W):
 634                out[y, x] = img_u8[y, x]
 635    
 636    return out
 637
 638
 639
 640@numba.njit(cache=True, fastmath=True)
 641def build_wall_mask_numba_values(img_u8, wall_values_i32):
 642    """
 643    Build a 0/255 wall mask from a uint8 image using explicit wall label values.
 644
 645    Each pixel in the input image is compared against a list of wall label
 646    values. If a pixel matches any of the provided values, it is marked as
 647    a wall (255), otherwise it is marked as free space (0).
 648
 649    This replaces the typical `np.isin` operation with explicit loops to
 650    remain fully compatible with Numba's nopython mode.
 651
 652    Parameters:
 653    - img_u8 (np.ndarray): <br>
 654        2D uint8 image containing label values.
 655    - wall_values_i32 (np.ndarray): <br>
 656        1D array of int32 values that should be interpreted as walls.
 657
 658    Returns:
 659    - np.ndarray: <br>
 660        A uint8 wall mask with values 0 (free space) and 255 (wall).
 661    """
 662    # img_u8: uint8 2D
 663    # wall_values_i32: int32 1D, e.g. np.array([3,7,9], np.int32)
 664    H, W = img_u8.shape[0], img_u8.shape[1]
 665    out = np.empty((H, W), dtype=np.uint8)
 666    # number of wall values to compare against
 667    nvals = wall_values_i32.shape[0]
 668
 669    for y in range(H):
 670        for x in range(W):
 671            # current pixel value
 672            v = int(img_u8[y, x])
 673
 674            # check if this pixel matches any of the wall values
 675            is_wall = False
 676            for k in range(nvals):
 677                if v == int(wall_values_i32[k]):
 678                    is_wall = True
 679                    break
 680
 681            # write wall mask pixel
 682            out[y, x] = 255 if is_wall else 0
 683    return out
 684
 685
 686
 687def build_wall_mask(img: np.ndarray, wall_values=None) -> np.ndarray:
 688    """
 689    Build a 0/255 wall mask from an input image.
 690
 691    This function acts as a Python wrapper that prepares the input data and
 692    delegates the actual computation to Numba-optimized kernels.
 693
 694    If `wall_values` is None, the image is assumed to be mask-like and
 695    processed accordingly. Otherwise, the provided label values are used
 696    to determine which pixels represent walls.
 697
 698    The input image must be two-dimensional. If it is not of type uint8,
 699    it is converted before processing.
 700
 701    Parameters:
 702    - img (np.ndarray): <br>
 703        2D input image representing either a mask-like image or a label map.
 704    - wall_values (optional): <br>
 705        Iterable of label values that should be interpreted as walls.
 706
 707    Returns:
 708    - np.ndarray: <br>
 709        A uint8 wall mask with values 0 (free space) and 255 (wall).
 710    """
 711    if img.ndim != 2:
 712        raise ValueError("build_wall_mask: erwartet 2D img (H,W). RGB vorher konvertieren.")
 713
 714    if img.dtype != np.uint8:
 715        img_u8 = img.astype(np.uint8)
 716    else:
 717        img_u8 = img
 718
 719    # if no explicit wall values are given, fall back to a default wall detection routine
 720    if wall_values is None:
 721        return build_wall_mask_numba_no_values(img_u8)
 722
 723    # convert wall values to a compact int array for numba-friendly lookup
 724    wall_vals = np.asarray(wall_values, dtype=np.int32)
 725
 726    # build wall mask by checking pixels against the provided wall values
 727    return build_wall_mask_numba_values(img_u8, wall_vals)
 728
 729
 730
 731def get_wall_segments_from_mask(mask_255: np.ndarray, thickness: int = 1, approx_epsilon: float = 1.5) -> List[Segment]:
 732    """
 733    Extract wall boundary segments from a binary wall mask.
 734
 735    Runs edge detection (Canny) and contour extraction to find wall boundaries,
 736    then approximates contours to polylines and converts edges into Segment
 737    instances.
 738
 739    Parameters:
 740    - mask_255 (np.ndarray): <br>
 741        Binary wall mask with values 0 and 255.
 742    - thickness (int): <br>
 743        Optional dilation thickness applied to edges before contour extraction.
 744        Values > 1 thicken edges to improve contour continuity.
 745    - approx_epsilon (float): <br>
 746        Epsilon value for contour polygon approximation (cv2.approxPolyDP).
 747        Higher values simplify contours more aggressively.
 748
 749    Returns:
 750    - List[Segment]: <br>
 751        List of wall boundary segments in pixel coordinates.
 752    """
 753    # detect wall boundaries from the binary mask using Canny edge detection
 754    edges = cv2.Canny(mask_255, 100, 200)
 755
 756    # optionally thicken edges to merge small gaps and stabilize contour extraction
 757    if thickness and thickness > 1:
 758        k = np.ones((thickness, thickness), np.uint8)
 759        edges = cv2.dilate(edges, k, iterations=1)
 760
 761    # extract external contours of the edge map
 762    contours, _ = cv2.findContours(edges, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
 763
 764    walls: List[Segment] = []
 765    for c in contours:
 766        # simplify contour to a polygon with fewer points
 767        c2 = cv2.approxPolyDP(c, epsilon=approx_epsilon, closed=True)
 768
 769        # convert polygon points to float tuples (x, y)
 770        pts = [tuple(p[0].astype(float)) for p in c2]
 771        if len(pts) < 2:
 772            continue
 773
 774        # create segments along the polygon polyline
 775        for i in range(len(pts) - 1):
 776            (x1, y1), (x2, y2) = pts[i], pts[i + 1]
 777
 778            # skip degenerate segments
 779            if x1 == x2 and y1 == y2:
 780                continue
 781            walls.append(Segment(x1, y1, x2, y2))
 782
 783        # close the loop by connecting the last point back to the first
 784        # only do this if we have a valid polygon-like shape
 785        if len(pts) >= 3:
 786            (x1, y1), (x2, y2) = pts[-1], pts[0]
 787            if x1 != x2 or y1 != y2:
 788                walls.append(Segment(x1, y1, x2, y2))
 789
 790    return walls
 791
 792
 793
 794def segments_to_walls4(walls):
 795    """
 796    Convert a list of wall segments into a compact float32 array representation.
 797
 798    Transforms a Python list of `Segment` objects (with endpoints ax/ay and bx/by)
 799    into a contiguous NumPy array of shape (N, 4), where each row encodes one wall
 800    segment as:
 801
 802        [ax, ay, bx, by]
 803
 804    This representation is convenient for passing wall geometry into Numba-compiled
 805    kernels, which cannot efficiently work with Python objects or dataclasses.
 806
 807    Parameters:
 808    - walls (List[Segment]): <br>
 809        A list of wall segments represented as `Segment` objects.
 810
 811    Returns:
 812    - np.ndarray: <br>
 813        A float32 array of shape (N, 4) containing wall endpoints.
 814    """
 815    walls4 = np.empty((len(walls), 4), dtype=np.float32)
 816
 817    for i, w in enumerate(walls):
 818        walls4[i,0] = w.ax 
 819        walls4[i,1] = w.ay
 820        walls4[i,2] = w.bx
 821        walls4[i,3] = w.by
 822
 823    return walls4
 824
 825
 826
 827@numba.njit(cache=True, fastmath=True)
 828def dilate_binary_square(src, k):
 829    """
 830    Perform binary dilation using a square structuring element.
 831
 832    Dilates a binary image `src` (values 0/1) using a kxk square kernel for a
 833    single iteration. A pixel in the output is set to 1 if any pixel in its
 834    kxk neighborhood in the input is nonzero.
 835
 836    This implementation is written in a Numba-friendly style (explicit loops,
 837    no OpenCV dependency) and is useful for thickening walls in an occlusion map.
 838
 839    Parameters:
 840    - src (np.ndarray): <br>
 841        2D uint8 array containing binary values {0, 1}.
 842    - k (int): <br>
 843        Kernel size (wall thickness). Values <= 1 return a copy of `src`.
 844
 845    Returns:
 846    - np.ndarray: <br>
 847        A uint8 binary array (0/1) of the same shape as `src` after dilation.
 848    """
 849    # src: uint8 0/1
 850    # k: wall_thickness (>=1)
 851
 852    # for thickness 1 nothing to do, just return a copy
 853    if k <= 1:
 854        return src.copy()
 855
 856    H, W = src.shape[0], src.shape[1]
 857
 858    # radius of the square neighborhood around each pixel
 859    r = k // 2
 860
 861    # output mask after dilation
 862    out = np.zeros((H, W), dtype=np.uint8)
 863
 864    for y in range(H):
 865        # clamp vertical neighborhood to image bounds
 866        y0 = y - r
 867        y1 = y + r
 868        if y0 < 0: y0 = 0
 869        if y1 >= H: y1 = H - 1
 870
 871        for x in range(W):
 872            # clamp horizontal neighborhood to image bounds
 873            x0 = x - r
 874            x1 = x + r
 875            if x0 < 0: x0 = 0
 876            if x1 >= W: x1 = W - 1
 877
 878            # check if any pixel in the k x k neighborhood is set
 879            # this effectively performs a square dilation of the wall mask
 880            found = 0
 881            for yy in range(y0, y1 + 1):
 882                for xx in range(x0, x1 + 1):
 883                    if src[yy, xx] != 0:
 884                        found = 1
 885                        break
 886                if found == 1:
 887                    break
 888
 889            # mark output pixel as wall if any neighbor was a wall
 890            out[y, x] = found
 891    
 892    return out
 893
 894
 895
 896@numba.njit(cache=True, fastmath=True)
 897def apply_mask_to_binary(mask_255):
 898    """
 899    Convert a 0/255 wall mask into a binary occlusion map.
 900
 901    Creates a binary occlusion map from a wall mask where nonzero pixels indicate
 902    walls. The output contains 0 for free space and 1 for occluded pixels.
 903
 904    This function is implemented using explicit loops for Numba compatibility.
 905
 906    Parameters:
 907    - mask_255 (np.ndarray): <br>
 908        2D uint8 wall mask, typically with values {0, 255}.
 909
 910    Returns:
 911    - np.ndarray: <br>
 912        A uint8 binary array of shape (H, W) with values:
 913        - 0 for free space
 914        - 1 for occluded (wall) pixels
 915    """
 916    H, W = mask_255.shape[0], mask_255.shape[1]
 917    occ = np.empty((H, W), dtype=np.uint8)
 918    for y in range(H):
 919        for x in range(W):
 920            occ[y, x] = 1 if mask_255[y, x] != 0 else 0
 921    return occ
 922
 923
 924
 925@numba.njit(cache=True, fastmath=True)
 926def build_occlusion_from_wallmask(mask_255, wall_thickness=1):
 927    """
 928    Build a binary occlusion map from a 0/255 wall mask, with optional wall dilation.
 929
 930    First converts the input wall mask (0/255) into a binary occlusion map (0/1).
 931    If `wall_thickness` is greater than 1, the binary map is dilated using a
 932    square kernel to increase the effective wall thickness for subsequent
 933    visibility checks.
 934
 935    All computation is performed in Numba-compatible kernels without relying
 936    on OpenCV operations.
 937
 938    Parameters:
 939    - mask_255 (np.ndarray): <br>
 940        2D uint8 wall mask, typically with values {0, 255}.
 941    - wall_thickness (int): <br>
 942        If > 1, applies a single binary dilation step with a square kernel of size
 943        (wall_thickness x wall_thickness).
 944
 945    Returns:
 946    - np.ndarray: <br>
 947        A uint8 binary occlusion map with values:
 948        - 0 for free space
 949        - 1 for occluded (wall) pixels
 950    """
 951    occ = apply_mask_to_binary(mask_255)
 952    if wall_thickness > 1:
 953        occ = dilate_binary_square(occ, wall_thickness)
 954    return occ
 955
 956
 957
 958# -----------------------------------------
 959# >>> ISM sequence enumeration & precompute
 960# -----------------------------------------
 961
 962def enumerate_wall_sequences_indices(n_walls: int, max_order: int, forbid_immediate_repeat: bool = True) -> List[Tuple[int, ...]]:
 963    """
 964    Enumerate reflection sequences over wall indices up to a given order.
 965
 966    Generates all sequences of wall indices representing reflection orders
 967    from 0 up to `max_order`. The empty sequence () represents the direct path.
 968
 969    Parameters:
 970    - n_walls (int): <br>
 971        Number of available wall segments.
 972    - max_order (int): <br>
 973        Maximum reflection order (length of sequence).
 974    - forbid_immediate_repeat (bool): <br>
 975        If True, prevents sequences with the same wall repeated consecutively,
 976        e.g., (..., 3, 3) is disallowed.
 977
 978    Returns:
 979    - List[Tuple[int, ...]]: <br>
 980        List of sequences, each a tuple of wall indices.
 981    """
 982    # start with the empty sequence which represents the direct path (no reflections)
 983    seqs = [()]
 984
 985    # iteratively build sequences up to the desired reflection order
 986    for _ in range(max_order):
 987        new = []
 988
 989        # extend all previously known sequences by one additional wall index
 990        for seq in seqs:
 991            for wi in range(n_walls):
 992
 993                # optionally forbid using the same wall twice in a row
 994                # this avoids degenerate mirror-back reflections
 995                if forbid_immediate_repeat and len(seq) > 0 and seq[-1] == wi:
 996                    continue
 997
 998                # create a new sequence with one more reflection
 999                new.append(seq + (wi,))
1000
1001        # append the newly created sequences to the master list
1002        seqs += new
1003
1004    # contains sequences of length 0..max_order
1005    return seqs
1006
1007
1008def precompute_image_sources(
1009    source_xy: Tuple[float, float],
1010    walls: List[Segment],
1011    max_order: int,
1012    forbid_immediate_repeat: bool = True,
1013    max_candidates: Optional[int] = None,
1014) -> List[Tuple[Tuple[int, ...], Tuple[float, float]]]:
1015    """
1016    Precompute image source positions for reflection sequences.
1017
1018    For every reflection sequence up to `max_order`, the source point is
1019    reflected across the corresponding wall lines to produce an image source
1020    position S_img.
1021
1022    Parameters:
1023    - source_xy (Tuple[float, float]): <br>
1024        Source position in pixel coordinates (x, y).
1025    - walls (List[Segment]): <br>
1026        Wall segments used for reflection.
1027    - max_order (int): <br>
1028        Maximum reflection order to consider.
1029    - forbid_immediate_repeat (bool): <br>
1030        If True, prevents immediate repetition of the same wall in sequences.
1031    - max_candidates (Optional[int]): <br>
1032        If provided, truncates the generated sequence list to at most this many
1033        candidates (useful as a speed cap).
1034
1035    Returns:
1036    - List[Tuple[Tuple[int, ...], Tuple[float, float]]]: <br>
1037        A list of tuples (seq, S_img) where: 
1038        - seq is the wall-index sequence
1039        - S_img is the resulting image source position
1040    """
1041    # enumerate all candidate wall index sequences up to the desired reflection order
1042    seqs = enumerate_wall_sequences_indices(len(walls), max_order, forbid_immediate_repeat)
1043
1044    # optionally cap how many sequences we keep to limit runtime and memory
1045    if max_candidates is not None and len(seqs) > max_candidates:
1046        seqs = seqs[:max_candidates]
1047
1048    # build a list of (sequence, image_source_xy) pairs
1049    # the image source is obtained by reflecting the real source across each wall
1050    # in the sequence in forward order
1051    pre: List[Tuple[Tuple[int, ...], Tuple[float, float]]] = []
1052    for seq in seqs:
1053        # start from the real source position
1054        S_img = source_xy
1055
1056        # apply successive reflections across each wall line
1057        for wi in seq:
1058            w = walls[wi]
1059            S_img = reflect_point_across_infinite_line(S_img, w.A, w.B)
1060
1061        # store the sequence and its final image source position
1062        pre.append((seq, S_img))
1063
1064    return pre
1065
1066
1067
1068def pack_precomputed(pre_py):
1069    """
1070    Pack Python reflection sequences into flat arrays for Numba kernels.
1071
1072    The input is a list of tuples where each entry contains:<br>
1073        (sequence_of_wall_indices, (image_source_x, image_source_y))
1074
1075    Since Numba cannot efficiently work with Python lists of variable-length
1076    tuples, this function converts the structure into contiguous flat arrays
1077    plus offset and length metadata so each sequence can be accessed quickly.
1078
1079    The result allows any sequence `i` to be reconstructed as:<br>
1080        seq_data[ seq_off[i] : seq_off[i] + seq_len[i] ]
1081
1082    Parameters:
1083    - pre_py (List[Tuple[Tuple[int], Tuple[float, float]]]): <br>
1084        Python list where each element contains a wall index sequence and
1085        the corresponding precomputed image source position
1086
1087    Returns:
1088    - seq_data (np.ndarray): <br>
1089        Flat int32 array containing all wall indices of all sequences
1090        concatenated back-to-back
1091    - seq_off (np.ndarray): <br>
1092        int32 array giving the start offset of each sequence in `seq_data`
1093    - seq_len (np.ndarray): <br>
1094        int32 array giving the length of each sequence
1095    - simg_x (np.ndarray): <br>
1096        float32 array of image source x-coordinates per sequence
1097    - simg_y (np.ndarray): <br>
1098        float32 array of image source y-coordinates per sequence
1099    """
1100    n = len(pre_py)
1101
1102    # allocate metadata arrays per sequence
1103    seq_len = np.empty(n, dtype=np.int32)
1104    seq_off = np.empty(n, dtype=np.int32)
1105    simg_x  = np.empty(n, dtype=np.float32)
1106    simg_y  = np.empty(n, dtype=np.float32)
1107
1108    # first pass computes lengths, offsets and image source positions
1109    total = 0
1110    for i, (seq_tuple, S_img) in enumerate(pre_py):
1111        L = len(seq_tuple)
1112        seq_len[i] = L
1113        seq_off[i] = total
1114        total += L
1115        simg_x[i] = float(S_img[0])
1116        simg_y[i] = float(S_img[1])
1117
1118    # allocate the flat wall-index array
1119    seq_data = np.empty(total, dtype=np.int32)
1120
1121    # second pass writes all sequences consecutively into seq_data
1122    k = 0
1123    for (seq_tuple, _S_img) in pre_py:
1124        for wi in seq_tuple:
1125            seq_data[k] = int(wi)
1126            k += 1
1127
1128    return seq_data, seq_off, seq_len, simg_x, simg_y
1129
1130
1131
1132# ----------------------------------
1133# >>> ISM path building (no shapely)
1134# ----------------------------------
1135
1136@numba.njit(cache=True, fastmath=True)
1137def build_path_for_sequence_packed(
1138        sx, sy, rx, ry,
1139        seq_data, off0, seq_len,
1140        s_imgx, s_imgy,
1141        walls4,
1142        path_out
1143    ):
1144    """
1145    Construct a specular reflection path for a packed wall-index sequence.
1146
1147    Builds the geometric reflection path between a real source position and
1148    a receiver position using a precomputed image source and a sequence of
1149    wall indices describing the reflection order.
1150
1151    The algorithm works backwards from the receiver by repeatedly
1152    1. Intersecting the line from the image source to the current virtual
1153      receiver with the corresponding wall segment
1154    2. Reflecting the virtual receiver across that wall line
1155
1156    This produces reflection points in reverse order which are then reversed
1157    in-place to form the final path
1158
1159        [source, r1, r2, ..., receiver]
1160
1161    The function is designed for Numba and writes into a preallocated array
1162    `path_out` instead of creating Python lists.
1163
1164    Parameters:
1165    - sx (float): <br>
1166        x-coordinate of the real source position
1167    - sy (float): <br>
1168        y-coordinate of the real source position
1169    - rx (float): <br>
1170        x-coordinate of the receiver position
1171    - ry (float): <br>
1172        y-coordinate of the receiver position
1173    - seq_data (np.ndarray): <br>
1174        Flat int array containing all sequences concatenated
1175    - off0 (int): <br>
1176        Start offset of this sequence inside `seq_data`
1177    - seq_len (int): <br>
1178        Number of wall indices in this sequence
1179    - s_imgx (float): <br>
1180        x-coordinate of the precomputed image source for this sequence
1181    - s_imgy (float): <br>
1182        y-coordinate of the precomputed image source for this sequence
1183    - walls4 (np.ndarray): <br>
1184        Array of shape (N,4) containing wall segments as [ax, ay, bx, by]
1185    - path_out (np.ndarray): <br>
1186        Preallocated float array of shape (>= seq_len+2, 2) used to store
1187        the resulting path points
1188
1189    Returns:
1190    - Tuple[int, bool]: <br>
1191        (path_length, ok) where `path_length` is the number of valid points
1192        written into `path_out`, and `ok` indicates whether a valid path
1193        could be constructed
1194    """
1195    # direct path with no reflections
1196    if seq_len == 0:
1197        path_out[0, 0] = sx; path_out[0, 1] = sy
1198        path_out[1, 0] = rx; path_out[1, 1] = ry
1199        return 2, True
1200
1201    # start from the real receiver as the current virtual receiver
1202    rvirt_x = rx
1203    rvirt_y = ry
1204
1205    for kk in range(seq_len - 1, -1, -1):
1206        # wall index for this bounce
1207        wi = seq_data[off0 + kk]
1208
1209        # wall segment endpoints
1210        ax = walls4[wi, 0]; ay = walls4[wi, 1]
1211        bx = walls4[wi, 2]; by = walls4[wi, 3]
1212
1213        # intersect ray from image source to virtual receiver with the wall segment
1214        ix, iy, ok = seg_seg_intersection_xy(
1215            s_imgx, s_imgy, rvirt_x, rvirt_y,
1216            ax, ay, bx, by
1217        )
1218        if not ok:
1219            return 0, False
1220
1221        # write intersection as a reflection point in reverse order
1222        path_out[1 + (seq_len - 1 - kk), 0] = ix
1223        path_out[1 + (seq_len - 1 - kk), 1] = iy
1224
1225        # update the virtual receiver by reflecting it across the infinite wall line
1226        rvirt_x, rvirt_y = reflect_point_across_infinite_line_numba(
1227            rvirt_x, rvirt_y, ax, ay, bx, by
1228        )
1229
1230    # reverse reflection points so they go from first bounce to last bounce
1231    i = 1
1232    j = seq_len
1233    while i < j:
1234        tmpx = path_out[i, 0]; tmpy = path_out[i, 1]
1235        path_out[i, 0] = path_out[j, 0]; path_out[i, 1] = path_out[j, 1]
1236        path_out[j, 0] = tmpx;          path_out[j, 1] = tmpy
1237        i += 1
1238        j -= 1
1239
1240    # insert real endpoints to finalize the path
1241    path_out[0, 0] = sx; path_out[0, 1] = sy
1242    path_out[seq_len + 1, 0] = rx
1243    path_out[seq_len + 1, 1] = ry
1244    return seq_len + 2, True
1245
1246
1247
1248@numba.njit(cache=True, fastmath=True)
1249def check_path_visibility_raster_arr(path_xy: np.ndarray, path_len: int, occ: np.ndarray, ignore_ends: int) -> bool:
1250    """
1251    Check whether a reconstructed reflection path is fully visible
1252    in a rasterized occlusion map.
1253
1254    The path consists of consecutive line segments:
1255        [source -> r1 -> r2 -> ... -> receiver]
1256
1257    Each segment is tested against the binary occlusion grid `occ`
1258    using a raster visibility test. If any segment intersects an
1259    occupied pixel (wall), the whole path is considered invalid.
1260
1261    This function is designed for Numba and operates on a preallocated
1262    path buffer without creating Python objects.
1263
1264    Parameters:
1265    - path_xy (np.ndarray): <br>
1266        Array of shape (max_order+2, 2) containing path points as (x, y).
1267        Only the first `path_len` entries are valid.
1268    - path_len (int): <br>
1269        Number of valid points in `path_xy` describing the path.
1270    - occ (np.ndarray): <br>
1271        2D uint8 occlusion map where non-zero values represent walls.
1272    - ignore_ends (int): <br>
1273        Number of pixels near segment endpoints to ignore during
1274        raster testing. Helps avoid false occlusion from grazing
1275        or endpoint-touching cases.
1276
1277    Returns:
1278    - bool: <br>
1279        True if all path segments are visible, False otherwise.
1280    """
1281    # paths with 0 or 1 point cannot be occluded
1282    if path_len <= 1:
1283        return True
1284    
1285    # check each consecutive segment of the path
1286    for i in range(path_len - 1):
1287        p1x = path_xy[i, 0]
1288        p1y = path_xy[i, 1]
1289        p2x = path_xy[i + 1, 0]
1290        p2y = path_xy[i + 1, 1]
1291
1292        # raster visibility test for this segment
1293        if not is_visible_raster(p1x, p1y, p2x, p2y, occ, ignore_ends):
1294            return False
1295        
1296    # all segments are visible
1297    return True
1298
1299
1300
1301def build_receivers_grid(H: int, W: int, step_px: int) -> np.ndarray:
1302    """
1303    Build a sparse grid of receiver positions over the image domain.
1304
1305    The receivers are placed on a regular grid with spacing `step_px`
1306    in both x and y directions. Each receiver is positioned at the
1307    pixel center using the convention (x + 0.5, y + 0.5), which matches
1308    how raster pixels are typically interpreted in the reflection kernel.
1309
1310    The result is a flat list of 2D receiver coordinates that can be
1311    iterated efficiently inside Numba kernels without additional
1312    meshgrid logic.
1313
1314    Parameters:
1315    - H (int): <br>
1316        Image height in pixels.
1317    - W (int): <br>
1318        Image width in pixels.
1319    - step_px (int): <br>
1320        Spacing between receivers in pixels. Larger values reduce the
1321        number of receivers and speed up computation at the cost of
1322        spatial resolution.
1323
1324    Returns:
1325    - np.ndarray: <br>
1326        Array of shape (N, 2) with dtype float32 containing receiver
1327        coordinates as (x+0.5, y+0.5).
1328    """
1329    # generate x and y sample positions at pixel centers
1330    xs = (np.arange(0, W, step_px, dtype=np.float32) + 0.5)
1331    ys = (np.arange(0, H, step_px, dtype=np.float32) + 0.5)
1332    
1333    # create a full grid of receiver coordinates
1334    xx, yy = np.meshgrid(xs, ys)
1335
1336    # pack into a flat (N,2) array for efficient iteration
1337    receivers = np.empty((xx.size, 2), dtype=np.float32)
1338    receivers[:, 0] = xx.reshape(-1)
1339    receivers[:, 1] = yy.reshape(-1)
1340    return receivers
1341
1342
1343
1344# -------------------------------
1345# >>> Main: noise / hit maps
1346# -------------------------------
1347
1348@numba.njit(cache=True, fastmath=True, parallel=True)
1349def eval_receivers_to_map_kernel(
1350    receivers: np.ndarray,      # (N,2) float32; contains x+0.5, y+0.5
1351    sx: float, sy: float,
1352    walls4: np.ndarray,         # (nWalls,4) float32
1353    occ: np.ndarray,            # (H,W) uint8
1354    max_order: int,
1355    ignore_zero_order: int,     # 0/1
1356    ignore_ends: int,
1357    seq_data: np.ndarray,       # (total,) int32
1358    seq_off: np.ndarray,        # (nSeq,)  int32  (start offset per seq)
1359    seq_len: np.ndarray,        # (nSeq,)  int32
1360    simg_x: np.ndarray,         # (nSeq,)  float32
1361    simg_y: np.ndarray,         # (nSeq,)  float32
1362    start_seq: int,
1363    end_seq: int,
1364    output_map: np.ndarray      # (H,W) float32
1365):
1366    H = output_map.shape[0]
1367    W = output_map.shape[1]
1368    nR = receivers.shape[0]
1369
1370    # iterate over all receiver positions in parallel
1371    for ri in numba.prange(nR):
1372        rx = float(receivers[ri, 0])
1373        ry = float(receivers[ri, 1])
1374
1375        # convert receiver center coordinate back to integer pixel index
1376        ix = int(rx)
1377        iy = int(ry)
1378
1379        # skip if receiver lies outside the image bounds
1380        if ix < 0 or ix >= W or iy < 0 or iy >= H:
1381            continue
1382
1383        # accumulate contribution for this receiver pixel
1384        val = 0.0
1385
1386        # local buffer for reconstructed reflection path
1387        # size is max_order reflections + source + receiver
1388        path_out = np.empty((max_order + 2, 2), dtype=np.float32)
1389
1390        # iterate over the subset of precomputed reflection sequences
1391        for si in range(start_seq, end_seq):
1392            L = int(seq_len[si])
1393            # optionally skip direct line-of-sight paths
1394            if ignore_zero_order == 1 and L == 0:
1395                continue
1396
1397            # offset into the packed wall-index array for this sequence
1398            off0 = int(seq_off[si])
1399
1400            # precomputed image source position for this wall sequence
1401            s_imgx = float(simg_x[si])
1402            s_imgy = float(simg_y[si])
1403
1404            # reconstruct the geometric reflection path for this receiver
1405            path_len, ok = build_path_for_sequence_packed(
1406                sx, sy,
1407                rx, ry,
1408                seq_data, off0, L,
1409                s_imgx, s_imgy,
1410                walls4,
1411                path_out
1412            )
1413            if not ok:
1414                continue
1415
1416            # check whether the path is occluded by any wall pixels
1417            if not check_path_visibility_raster_arr(path_out, path_len, occ, ignore_ends):
1418                continue
1419
1420            # valid visible reflection path contributes to this pixel
1421            val += 1.0
1422
1423        # safe to write because each receiver maps to a unique pixel
1424        output_map[iy, ix] += val
1425
1426
1427
1428def compute_reflection_map(
1429    source_rel,
1430    img,
1431    wall_values,
1432    wall_thickness:int=1,
1433    approx_epsilon:float=1.5,
1434    max_order:int=1,
1435    ignore_zero_order=False,
1436    step_px:int=1,
1437    forbid_immediate_repeat: bool = True,
1438    max_candidates=None,
1439    ignore_ends:int=1,
1440    iterative_tracking=False,
1441    iterative_steps=None
1442):
1443    """
1444    Compute a 2D "reflection contribution" map over an image grid using
1445    precomputed specular image sources and wall reflections.
1446
1447    The function:
1448    1. Converts the input image to grayscale (if needed)
1449    2. Builds a wall mask from the image based on provided wall pixel values
1450    3. Extracts wall segments from that mask and builds an occlusion structure
1451    4. Precomputes valid wall-reflection sequences (up to `max_order`) and their
1452       corresponding image-source positions
1453    5. Evaluates reflection paths from a real source to a grid of receiver points,
1454       accumulating per-receiver contributions into `output_map`
1455
1456    Optionally, the computation can be performed in chunks and intermediate
1457    snapshots of the partially accumulated `output_map` can be returned
1458    (`iterative_tracking=True`).
1459
1460    Parameters:
1461    - source_rel (Tuple[float, float] | np.ndarray): <br>
1462        Source position in *relative* image coordinates (x_rel, y_rel) in [0, 1].
1463        Converted internally to pixel coordinates via (x_rel*W, y_rel*H).
1464    - img (np.ndarray): <br>
1465        Input image. If 3D (H, W, C), only the first channel is used.
1466        The wall mask is derived from this image.
1467    - wall_values (Iterable[int] | np.ndarray): <br>
1468        Pixel values in `img` that should be treated as "wall" when building the
1469        wall mask (exact matching logic depends on `build_wall_mask`).
1470    - wall_thickness (int): <br>
1471        Thickness (in pixels) used when extracting wall segments and building
1472        the occlusion representation.
1473    - approx_epsilon (float): <br>
1474        Polygon/segment simplification tolerance used during wall extraction
1475        (`get_wall_segments_from_mask`).
1476    - max_order (int): <br>
1477        Maximum reflection order (number of wall bounces) to consider when
1478        precomputing image sources.
1479    - ignore_zero_order (bool): <br>
1480        If True, exclude the direct (0-bounce / line-of-sight) contribution.
1481    - step_px (int): <br>
1482        Receiver grid spacing in pixels. Larger values reduce runtime by
1483        evaluating fewer receiver positions.
1484    - forbid_immediate_repeat (bool): <br>
1485        If True, disallow reflection sequences that repeat the same wall in
1486        consecutive steps (helps avoid degenerate/near-degenerate paths).
1487    - max_candidates (int | None): <br>
1488        Optional cap on the number of candidate sequences/image sources kept
1489        during precomputation.
1490    - ignore_ends (int): <br>
1491        Number of wall endpoints to ignore during intersection tests (used to
1492        stabilize grazing / endpoint-hitting cases).
1493    - iterative_tracking (bool): <br>
1494        If True, compute in chunks and return intermediate `output_map` snapshots
1495        as a list of arrays. Useful for debugging/profiling progress.
1496    - iterative_steps (int | None): <br>
1497        Controls how many snapshots to collect when `iterative_tracking=True`:
1498        - None  -> treated as -1 (snapshot every chunk of size 1 sequence)
1499        - -1    -> snapshot after every sequence (finest granularity)
1500        - 1     -> snapshot only once at the end (same as full run)
1501        - k>1   -> snapshot roughly `k` times over the sequence list
1502
1503    Returns:
1504    - np.ndarray | List[np.ndarray]: <br>
1505        If `iterative_tracking` is False: returns `output_map` (H, W) float32.<br>
1506        If `iterative_tracking` is True: returns a list of snapshot maps, each
1507        shaped (H, W) float32, showing accumulation progress.
1508    """
1509    # If we received a color image, use only the first channel (expected: walls encoded per-pixel)
1510    if img.ndim == 3:
1511        img = img[..., 0]
1512
1513    # Image dimensions and source location in absolute pixel coordinates
1514    H, W = img.shape[:2]
1515    sx = float(source_rel[0] * W)
1516    sy = float(source_rel[1] * H)
1517
1518    # Build a binary wall mask (0/255), then extract simplified wall segments
1519    wall_mask_255 = build_wall_mask(img, wall_values=wall_values)
1520    walls = get_wall_segments_from_mask(
1521        wall_mask_255,
1522        thickness=wall_thickness,
1523        approx_epsilon=approx_epsilon
1524    )
1525
1526    # Precompute an occlusion structure derived from the wall mask
1527    occ = build_occlusion_from_wallmask(wall_mask_255, wall_thickness=wall_thickness)
1528
1529    # Convert wall segments into a packed (N,4) representation [ax, ay, bx, by]
1530    #    -> for numba usage
1531    walls4 = segments_to_walls4(walls)
1532
1533    # Precompute image sources + wall reflection sequences up to `max_order`
1534    # Each sequence describes the ordered wall indices to reflect against
1535    pre_py = precompute_image_sources(
1536        source_xy=(sx, sy),
1537        walls=walls,
1538        max_order=max_order,
1539        forbid_immediate_repeat=forbid_immediate_repeat,
1540        max_candidates=max_candidates,
1541    )
1542
1543    # Pack variable-length sequences into flat arrays for fast/Numba-friendly access
1544    seq_data, seq_off, seq_len, simg_x, simg_y = pack_precomputed(pre_py)
1545
1546    # Create a grid of receiver points
1547    receivers = build_receivers_grid(H, W, step_px)
1548
1549    # Accumulated output map (same resolution as the input image)
1550    output_map = np.zeros((H, W), dtype=np.float32)
1551
1552    # Optional progressive snapshots of accumulation -> iterative output
1553    if iterative_tracking:
1554        snapshots = []
1555
1556    # Default: snapshot every sequence (finest granularity) if tracking is enabled.
1557    if iterative_steps is None:
1558        iterative_steps = -1
1559
1560    nseq = len(pre_py)
1561
1562    # Decide chunk size (how many sequences per kernel call) based on snapshot settings
1563    if (not iterative_tracking) or (iterative_tracking and iterative_steps == 1):
1564        # Single call covering all sequences (or one snapshot at the end)
1565        snapshot_every_x_steps = nseq
1566    elif iterative_tracking and iterative_steps == -1:
1567        # Snapshot after every sequence
1568        snapshot_every_x_steps = 1
1569    else:
1570        # Snapshot about `iterative_steps` times, compute chunk size accordingly
1571        snapshot_every_x_steps = max(1, int(nseq // iterative_steps))
1572
1573    # Convert boolean to int for Numba kernel
1574    ignore_zero_order_i = 1 if ignore_zero_order else 0
1575
1576    # Evaluate reflection contributions in chunks of sequences for better control over snapshots
1577    for start_idx in range(0, nseq, snapshot_every_x_steps):
1578        end_idx = min(nseq, start_idx + snapshot_every_x_steps)
1579
1580        # Core kernel: for each receiver and each sequence in [start_idx, end_idx),
1581        # validate/specular-trace the path (including occlusion) and accumulate into output_map
1582        eval_receivers_to_map_kernel(
1583            receivers,
1584            sx, sy,
1585            walls4,
1586            occ,
1587            max_order,
1588            ignore_zero_order_i,
1589            ignore_ends,
1590            seq_data,
1591            seq_off,
1592            seq_len,
1593            simg_x,
1594            simg_y,
1595            start_idx,
1596            end_idx,
1597            output_map
1598        )
1599
1600        # Store a copy so later updates don't mutate previous snapshots
1601        if iterative_tracking:
1602            snapshots.append(output_map.copy())
1603
1604    # Return either the final map or the progressive snapshots
1605    if not iterative_tracking:
1606        return output_map
1607    else:
1608        return snapshots