Building blocks

The primitives the pipelines are made of. Each has the same name and argument order on every backend, so a chain written against one runs on another.

See use the building blocks for worked examples.

On the processor

vidmag.cpu.ops

The same building blocks as :mod:vidmag.cuda.ops, computed with NumPy.

This module adds no arithmetic of its own. It is a naming layer: the functions here forward to the implementations in :mod:vidmag.cpu.pyramids and :mod:vidmag.cpu.filters, under the names and argument order the GPU operations use. Two reasons that is worth having.

A caller can write one chain of operations and choose the backend separately, because vidmag.cpu.ops.blur_dn and vidmag.cuda.ops.blur_dn take the same arguments in the same order.

And the conformance tests can run the same list of operations against every backend, comparing each against this one. That makes this module the reference the GPU is checked against, which is the role the whole project already gives the NumPy code.

level_sizes(height, width, levels)

The (height, width) of each pyramid level, finest first.

Source code in src/vidmag/cpu/ops.py

def level_sizes(height: int, width: int, levels: int) -> list[tuple[int, int]]:
    """The (height, width) of each pyramid level, finest first."""
    out, h, w = [], height, width
    for _ in range(levels):
        out.append((h, w))
        h, w = (h + 1) // 2, (w + 1) // 2
    return out

bgr_u8_to_ntsc(frames)

Convert (T, H, W, 3) 8-bit blue-green-red frames to NTSC, float32.

The per-frame helper takes 8-bit input and does the channel reversal and the divide by 255 itself, so frames are handed to it unchanged.

Source code in src/vidmag/cpu/ops.py

def bgr_u8_to_ntsc(frames: np.ndarray) -> np.ndarray:
    """Convert (T, H, W, 3) 8-bit blue-green-red frames to NTSC, float32.

    The per-frame helper takes 8-bit input and does the channel reversal and
    the divide by 255 itself, so frames are handed to it unchanged.
    """
    if frames.dtype != np.uint8:
        raise TypeError(f"bgr_u8_to_ntsc: expected uint8, got {frames.dtype}")
    return np.stack([_rgb_frame_to_ntsc(f) for f in frames], axis=0).astype(np.float32)

ntsc_to_bgr_u8(frames)

Convert NTSC frames back to 8-bit blue-green-red, with rounding.

Source code in src/vidmag/cpu/ops.py

def ntsc_to_bgr_u8(frames: np.ndarray) -> np.ndarray:
    """Convert NTSC frames back to 8-bit blue-green-red, with rounding."""
    return np.stack([_ntsc_to_bgr_uint8(f) for f in frames], axis=0)

blur_dn(frames, levels)

Blur and halve the resolution, levels times.

Source code in src/vidmag/cpu/ops.py

def blur_dn(frames: np.ndarray, levels: int) -> np.ndarray:
    """Blur and halve the resolution, ``levels`` times."""
    return np.stack([_pyramids.blur_dn_clr(f, levels) for f in frames], axis=0)

build_lpyr(frames, levels)

Build a Laplacian pyramid: one band per scale, finest first.

Returned bands are shaped (3, T, height, width) to match the GPU layout, which keeps each colour channel's planes together.

Source code in src/vidmag/cpu/ops.py

def build_lpyr(frames: np.ndarray, levels: int) -> list[np.ndarray]:
    """Build a Laplacian pyramid: one band per scale, finest first.

    Returned bands are shaped (3, T, height, width) to match the GPU layout,
    which keeps each colour channel's planes together.
    """
    # The per-frame helper returns (levels, pind); each entry of `levels` is
    # one pyramid band shaped (height, width, channel), so the channel axis has
    # to be moved to the front to match the GPU's plane-per-channel layout.
    per_frame = [_pyramids.laplacian_pyramid_channels(f, levels)[0] for f in frames]
    bands = []
    for level in range(levels):
        stacked = np.stack([per_frame[t][level] for t in range(len(frames))])
        bands.append(np.ascontiguousarray(np.moveaxis(stacked, 3, 0), dtype=np.float32))
    return bands

recon_lpyr(bands, height, width)

Sum a Laplacian pyramid back, returning one plane per channel.

Bands arrive shaped (channel, time, height, width), matching the GPU layout. The per-frame helper wants the opposite: a list of levels each shaped (height, width, channel), plus the table of level dimensions it calls pind. Both are rebuilt here.

Source code in src/vidmag/cpu/ops.py

def recon_lpyr(bands: list[np.ndarray], height: int, width: int) -> np.ndarray:
    """Sum a Laplacian pyramid back, returning one plane per channel.

    Bands arrive shaped (channel, time, height, width), matching the GPU
    layout. The per-frame helper wants the opposite: a list of levels each
    shaped (height, width, channel), plus the table of level dimensions it
    calls ``pind``. Both are rebuilt here.
    """
    T = bands[0].shape[1]
    pind = np.array([[b.shape[2], b.shape[3]] for b in bands], dtype=np.float64)
    out = np.empty((T * 3, height, width), dtype=np.float32)
    for t in range(T):
        levels = [np.moveaxis(b[:, t], 0, 2) for b in bands]
        frame = _pyramids.reconstruct_from_channels(levels, pind)
        for c in range(3):
            out[c * T + t] = frame[:, :, c]
    return out

ideal_bandpass(frames, fl, fh, sampling_rate)

Keep only frequencies strictly between fl and fh, along time.

Source code in src/vidmag/cpu/ops.py

def ideal_bandpass(
    frames: np.ndarray, fl: float, fh: float, sampling_rate: float
) -> np.ndarray:
    """Keep only frequencies strictly between ``fl`` and ``fh``, along time."""
    return _filters.ideal_bandpass(frames, fl, fh, sampling_rate, axis=0)

butter_bandpass(frames, fl, fh, sampling_rate, order=1)

First-order Butterworth bandpass, along time.

Source code in src/vidmag/cpu/ops.py

def butter_bandpass(
    frames: np.ndarray, fl: float, fh: float, sampling_rate: float, order: int = 1
) -> np.ndarray:
    """First-order Butterworth bandpass, along time."""
    return _filters.butter_bandpass(frames, fl, fh, sampling_rate, order=order, axis=0)

iir_bandpass(frames, r1, r2)

The difference of two exponential moving averages, along time.

Source code in src/vidmag/cpu/ops.py

def iir_bandpass(frames: np.ndarray, r1: float, r2: float) -> np.ndarray:
    """The difference of two exponential moving averages, along time."""
    return _filters.iir_bandpass(frames, r1, r2, axis=0)

apply_gain(frames, gain_y, gain_i, gain_q)

Scale the three NTSC channels independently.

Source code in src/vidmag/cpu/ops.py

def apply_gain(
    frames: np.ndarray, gain_y: float, gain_i: float, gain_q: float
) -> np.ndarray:
    """Scale the three NTSC channels independently."""
    return frames * np.array([gain_y, gain_i, gain_q], dtype=frames.dtype)

The pipelines, derived from the primitives

vidmag.backend.generic

The four magnification pipelines, written once against :class:Ops.

A backend only has to supply the primitive operations in :mod:vidmag.backend.ops — colour conversion, blur and downsample, pyramid build and reconstruct, three temporal filters, gain, upsample, quantize. The four pipelines then come from here for free, which is what makes adding hardware support a bounded job rather than a rewrite.

A backend may still replace any of these with its own version, and the two that matter for speed do: the hand-written CUDA code fuses stages and collapses kernel launches, which nothing expressed as a sequence of separate operations can match. The functions here are the correct-by-construction fallback, and the thing every new backend is measured against before it earns an override.

The arithmetic is the same as :mod:vidmag.cpu.magnify; only the spelling differs, because here every step goes through the protocol rather than calling NumPy directly.

color_gdown_ideal_core(ops, frames_bgr_u8, fps, *, alpha, level, fl, fh, chrom_attenuation=1.0, sampling_rate=None)

Amplify colour change: blur down, bandpass over time, scale, add back.

This is the pipeline that makes a pulse visible. Reducing the resolution first is what suppresses noise: a heartbeat changes a whole region of skin together, so averaging over a region keeps the signal and discards most of what is random.

Source code in src/vidmag/backend/generic.py

def color_gdown_ideal_core(
    ops: Any,
    frames_bgr_u8: np.ndarray,
    fps: float,
    *,
    alpha: float,
    level: int,
    fl: float,
    fh: float,
    chrom_attenuation: float = 1.0,
    sampling_rate: float | None = None,
) -> np.ndarray:
    """Amplify colour change: blur down, bandpass over time, scale, add back.

    This is the pipeline that makes a pulse visible. Reducing the resolution
    first is what suppresses noise: a heartbeat changes a whole region of skin
    together, so averaging over a region keeps the signal and discards most of
    what is random.
    """
    rate = _resolve_rate(fps, sampling_rate)
    frames = ops.from_numpy(np.ascontiguousarray(frames_bgr_u8))

    ntsc = ops.bgr_u8_to_ntsc(frames)
    small = ops.blur_dn(ntsc, level)
    filtered = ops.ideal_bandpass(small, fl, fh, rate)
    amplified = ops.apply_gain(
        filtered, alpha, alpha * chrom_attenuation, alpha * chrom_attenuation
    )

    _, height, width, _ = ntsc.shape
    delta = ops.upsample_bilinear(amplified, height, width)
    out: np.ndarray = ops.to_numpy(ops.add_and_quantize(ntsc, delta))
    return out

motion_lpyr_ideal_core(ops, frames_bgr_u8, fps, *, alpha, lambda_c, fl, fh, chrom_attenuation=0.0, sampling_rate=None, exaggeration_factor=EXAGGERATION_FACTOR)

Amplify motion, selecting the frequency band with a Fourier transform.

Needs the whole clip, because the transform runs over all of time at once.

Source code in src/vidmag/backend/generic.py

def motion_lpyr_ideal_core(
    ops: Any,
    frames_bgr_u8: np.ndarray,
    fps: float,
    *,
    alpha: float,
    lambda_c: float,
    fl: float,
    fh: float,
    chrom_attenuation: float = 0.0,
    sampling_rate: float | None = None,
    exaggeration_factor: float = EXAGGERATION_FACTOR,
) -> np.ndarray:
    """Amplify motion, selecting the frequency band with a Fourier transform.

    Needs the whole clip, because the transform runs over all of time at once.
    """
    rate = _resolve_rate(fps, sampling_rate)
    return _motion_core(
        ops,
        frames_bgr_u8,
        fps,
        alpha=alpha,
        lambda_c=lambda_c,
        chrom_attenuation=chrom_attenuation,
        exaggeration_factor=exaggeration_factor,
        filter_band=lambda band: ops.ideal_bandpass(band, fl, fh, rate),
    )

motion_lpyr_butter_core(ops, frames_bgr_u8, fps, *, alpha, lambda_c, fl, fh, chrom_attenuation=0.0, sampling_rate=None, order=1, exaggeration_factor=EXAGGERATION_FACTOR)

Amplify motion, selecting the band with a Butterworth filter.

Runs forward in time only, so it also works on frames as they arrive.

Source code in src/vidmag/backend/generic.py

def motion_lpyr_butter_core(
    ops: Any,
    frames_bgr_u8: np.ndarray,
    fps: float,
    *,
    alpha: float,
    lambda_c: float,
    fl: float,
    fh: float,
    chrom_attenuation: float = 0.0,
    sampling_rate: float | None = None,
    order: int = 1,
    exaggeration_factor: float = EXAGGERATION_FACTOR,
) -> np.ndarray:
    """Amplify motion, selecting the band with a Butterworth filter.

    Runs forward in time only, so it also works on frames as they arrive.
    """
    rate = _resolve_rate(fps, sampling_rate)
    return _motion_core(
        ops,
        frames_bgr_u8,
        fps,
        alpha=alpha,
        lambda_c=lambda_c,
        chrom_attenuation=chrom_attenuation,
        exaggeration_factor=exaggeration_factor,
        filter_band=lambda band: ops.butter_bandpass(band, fl, fh, rate, order),
    )

motion_lpyr_iir_core(ops, frames_bgr_u8, fps, *, alpha, lambda_c, r1, r2, chrom_attenuation=0.1, exaggeration_factor=EXAGGERATION_FACTOR)

Amplify motion, selecting the band by subtracting two running averages.

The cheapest of the three and the only one that needs no history beyond the previous frame, which is what makes it the one a live stream can use.

Source code in src/vidmag/backend/generic.py

def motion_lpyr_iir_core(
    ops: Any,
    frames_bgr_u8: np.ndarray,
    fps: float,
    *,
    alpha: float,
    lambda_c: float,
    r1: float,
    r2: float,
    chrom_attenuation: float = 0.1,
    exaggeration_factor: float = EXAGGERATION_FACTOR,
) -> np.ndarray:
    """Amplify motion, selecting the band by subtracting two running averages.

    The cheapest of the three and the only one that needs no history beyond the
    previous frame, which is what makes it the one a live stream can use.
    """
    return _motion_core(
        ops,
        frames_bgr_u8,
        fps,
        alpha=alpha,
        lambda_c=lambda_c,
        chrom_attenuation=chrom_attenuation,
        exaggeration_factor=exaggeration_factor,
        filter_band=lambda band: ops.iir_bandpass(band, r1, r2),
    )

bind(ops)

Give a set of primitive operations all four pipelines.

This is what makes supporting new hardware a bounded job: implement the operations in :mod:vidmag.backend.ops, call this, and the backend is complete.

Source code in src/vidmag/backend/generic.py

def bind(ops: Any) -> _BoundPipelines:
    """Give a set of primitive operations all four pipelines.

    This is what makes supporting new hardware a bounded job: implement the
    operations in :mod:`vidmag.backend.ops`, call this, and the backend is
    complete.
    """
    return _BoundPipelines(ops)

Video reading and writing

vidmag.io.video

Video I/O helpers.

Loads a video into a single float32 array of shape (T, H, W, C) with values in [0, 1] and writes one back out. Keeping the whole clip in memory is fine for the baseline (the EVM temporal filters need random access to all frames anyway) and makes the algorithm easy to read; the CUDA port will stream frames through device memory instead.

VideoInfo dataclass

Metadata needed to reconstruct an output file from a float array.

Source code in src/vidmag/io/video.py

@dataclass
class VideoInfo:
    """Metadata needed to reconstruct an output file from a float array."""

    fps: float
    width: int
    height: int
    frame_count: int
    is_color: bool

load_video(path)

Load a video as a float32 array of shape (T, H, W, C) in [0, 1].

Color videos come back as 3-channel BGR (OpenCV's native order, so we can hand frames straight back to save_video without permuting channels). Grayscale videos come back as (T, H, W, 1) so downstream code can assume a trailing channel axis unconditionally.

Source code in src/vidmag/io/video.py

def load_video(path: str | Path) -> tuple[np.ndarray, VideoInfo]:
    """Load a video as a float32 array of shape ``(T, H, W, C)`` in ``[0, 1]``.

    Color videos come back as 3-channel BGR (OpenCV's native order, so we can
    hand frames straight back to ``save_video`` without permuting channels).
    Grayscale videos come back as ``(T, H, W, 1)`` so downstream code can assume
    a trailing channel axis unconditionally.
    """
    path = str(path)
    cap = cv2.VideoCapture(path)
    if not cap.isOpened():
        raise FileNotFoundError(f"Could not open video: {path!r}")

    fps = cap.get(cv2.CAP_PROP_FPS) or 30.0
    width = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH))
    height = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT))

    frames: list[np.ndarray] = []
    while True:
        ok, frame = cap.read()
        if not ok:
            break
        if frame.ndim == 2:
            frame = frame[:, :, None]
        frames.append(frame)
    cap.release()

    if not frames:
        raise ValueError(f"No frames decoded from {path!r}")

    stack = np.stack(frames, axis=0).astype(np.float32) / 255.0
    info = VideoInfo(
        fps=float(fps),
        width=width or int(stack.shape[2]),
        height=height or int(stack.shape[1]),
        frame_count=int(stack.shape[0]),
        is_color=stack.shape[3] == 3,
    )
    return stack, info

save_video(frames, path, fps, *, codec='libx264')

Write a float32 array in [0, 1] back to an H.264 MP4.

frames is (T, H, W, C) with C in {1, 3}. Single-channel arrays are broadcast to 3 channels for the encoder. Values are clipped to the valid range and converted to uint8 (BGR), then handed to :func:encode_video.

Source code in src/vidmag/io/video.py

def save_video(
    frames: np.ndarray,
    path: str | Path,
    fps: float,
    *,
    codec: str = "libx264",
) -> None:
    """Write a float32 array in ``[0, 1]`` back to an H.264 MP4.

    ``frames`` is ``(T, H, W, C)`` with C in {1, 3}. Single-channel arrays are
    broadcast to 3 channels for the encoder. Values are clipped to the valid
    range and converted to ``uint8`` (BGR), then handed to :func:`encode_video`.
    """
    if frames.ndim != 4:
        raise ValueError(f"Expected (T,H,W,C); got shape {frames.shape!r}")

    c = frames.shape[3]
    clipped = np.clip(frames, 0.0, 1.0)
    scaled = np.round(clipped * 255.0).astype(np.uint8)
    if c == 1:
        # H.264/yuv420p needs a 3-channel source; replicate grayscale to BGR.
        scaled = np.repeat(scaled, 3, axis=3)

    encode_video(scaled, path, fps, codec=codec)

rgb_to_yiq(rgb)

Convert RGB float in [0, 1] to YIQ, matching MATLAB rgb2ntsc.

Uses the exact transform matrix documented for MATLAB's rgb2ntsc so the luminance/chrominance split is identical to the MIT reference. The row order is (Y, I, Q). Input shape (..., 3).

Source code in src/vidmag/io/video.py

def rgb_to_yiq(rgb: np.ndarray) -> np.ndarray:
    """Convert RGB float in ``[0, 1]`` to YIQ, matching MATLAB ``rgb2ntsc``.

    Uses the exact transform matrix documented for MATLAB's ``rgb2ntsc`` so the
    luminance/chrominance split is identical to the MIT reference. The row order
    is ``(Y, I, Q)``. Input shape ``(..., 3)``.
    """
    m = np.array(
        [
            [0.299, 0.587, 0.114],
            [-0.168736, -0.331264, 0.5],
            [0.5, -0.418688, -0.081312],
        ],
        dtype=np.float32,
    )
    return rgb @ m.T

yiq_to_rgb(yiq)

Inverse of :func:rgb_to_yiq, matching MATLAB ntsc2rgb.

Source code in src/vidmag/io/video.py

def yiq_to_rgb(yiq: np.ndarray) -> np.ndarray:
    """Inverse of :func:`rgb_to_yiq`, matching MATLAB ``ntsc2rgb``."""
    m = np.array(
        [
            [1.0, -1.21889419e-06, 1.40199959],
            [1.0, -3.44135678e-01, -7.14136156e-01],
            [1.0, 1.77200007, 4.06298063e-07],
        ],
        dtype=np.float32,
    )
    return yiq @ m.T