ivp_binary_sensor_mask.py

"""
Using A Binary Sensor Mask Example

Ported from: k-Wave/examples/example_ivp_binary_sensor_mask.m

Demonstrates how to use a binary sensor mask for detection of the pressure
field generated by an initial pressure distribution within a two-dimensional
homogeneous propagation medium. It builds on the Homogeneous Propagation
Medium Example.

Instead of the Cartesian sensor used in the homogeneous example, this example
uses a binary arc created with make_circle. You should observe the same two
expanding wavefronts, now sampled at grid-aligned sensor positions along a
three-quarter arc.
"""
# %%
import numpy as np

from kwave.data import Vector
from kwave.kgrid import kWaveGrid
from kwave.kmedium import kWaveMedium
from kwave.ksensor import kSensor
from kwave.ksource import kSource
from kwave.kspaceFirstOrder import kspaceFirstOrder
from kwave.utils.mapgen import make_circle, make_disc


# %%
def setup():
    """Set up the simulation physics (grid, medium, source).

    Returns:
        tuple: (kgrid, medium, source)
    """

    # create the computational grid
    Nx = 128  # number of grid points in the x direction
    Ny = 128  # number of grid points in the y direction
    dx = 0.1e-3  # grid point spacing in the x direction [m]
    dy = 0.1e-3  # grid point spacing in the y direction [m]
    kgrid = kWaveGrid(Vector([Nx, Ny]), Vector([dx, dy]))

    # define the properties of the propagation medium
    medium = kWaveMedium(
        sound_speed=1500,  # [m/s]
        alpha_coeff=0.75,  # [dB/(MHz^y cm)]
        alpha_power=1.5,
    )

    # create initial pressure distribution using make_disc
    # -- disc 1 --
    disc_1_magnitude = 5  # [Pa]
    disc_1_x_pos = 50  # [grid points, 1-based]
    disc_1_y_pos = 50  # [grid points, 1-based]
    disc_1_radius = 8  # [grid points]
    disc_1 = disc_1_magnitude * make_disc(Vector([Nx, Ny]), Vector([disc_1_x_pos, disc_1_y_pos]), disc_1_radius)

    # -- disc 2 --
    disc_2_magnitude = 3  # [Pa]
    disc_2_x_pos = 80  # [grid points, 1-based]
    disc_2_y_pos = 60  # [grid points, 1-based]
    disc_2_radius = 5  # [grid points]
    disc_2 = disc_2_magnitude * make_disc(Vector([Nx, Ny]), Vector([disc_2_x_pos, disc_2_y_pos]), disc_2_radius)

    source = kSource()
    source.p0 = (disc_1 + disc_2).astype(float)

    # set time array
    kgrid.makeTime(1500)

    return kgrid, medium, source


# %%
def run(backend="python", device="cpu", quiet=True):
    """Run with the original binary arc sensor.

    Returns:
        dict: Simulation results with key 'p' (sensor_points x time_steps).
    """
    kgrid, medium, source = setup()

    Nx = 128  # number of grid points in the x direction
    Ny = 128  # number of grid points in the y direction

    # define a binary sensor mask (arc)
    sensor_x_pos = Nx // 2  # [grid points, 1-based]
    sensor_y_pos = Ny // 2  # [grid points, 1-based]
    sensor_radius = Nx // 2 - 22  # [grid points]
    sensor_arc_angle = 3 * np.pi / 2  # [radians]
    sensor_mask = make_circle(
        Vector([Nx, Ny]),
        Vector([sensor_x_pos, sensor_y_pos]),
        sensor_radius,
        sensor_arc_angle,
    )

    sensor = kSensor(mask=sensor_mask)
    sensor.record = ["p"]

    return kspaceFirstOrder(
        kgrid,
        medium,
        source,
        sensor,
        backend=backend,
        device=device,
        quiet=quiet,
        pml_inside=True,
    )


# %%
if __name__ == "__main__":
    import matplotlib.pyplot as plt

    result = run(quiet=False)
    p = np.asarray(result["p"])

    Nx, Ny = 128, 128
    kgrid, _, source = setup()

    # build sensor mask for overlay
    sensor_mask = make_circle(
        Vector([Nx, Ny]),
        Vector([Nx // 2, Ny // 2]),
        Nx // 2 - 22,
        3 * np.pi / 2,
    )

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))

    # plot the initial pressure and sensor distribution
    ax = axes[0]
    overlay = source.p0 + 5.0 * sensor_mask
    ax.imshow(
        overlay.T,
        extent=[
            kgrid.x_vec[0] * 1e3,
            kgrid.x_vec[-1] * 1e3,
            kgrid.y_vec[-1] * 1e3,
            kgrid.y_vec[0] * 1e3,
        ],
        vmin=-1,
        vmax=1,
        cmap="RdBu_r",
    )
    ax.set_ylabel("x-position [mm]")
    ax.set_xlabel("y-position [mm]")
    ax.set_title("Initial Pressure + Sensor Mask")
    ax.set_aspect("equal")

    # plot the simulated sensor data
    ax = axes[1]
    im = ax.imshow(p, aspect="auto", vmin=-1, vmax=1)
    ax.set_ylabel("Sensor Position")
    ax.set_xlabel("Time Step")
    ax.set_title("Sensor Data")
    fig.colorbar(im, ax=ax)

    fig.suptitle("Binary Sensor Mask Example")
    fig.tight_layout()
    plt.show()