|
| 1 | +"""Phantom generation and loading commands.""" |
| 2 | + |
| 3 | +import click |
| 4 | +import numpy as np |
| 5 | + |
| 6 | +from kwave.cli.main import pass_session |
| 7 | +from kwave.cli.schema import CLIError, CLIResponse, ValidationError, json_command |
| 8 | + |
| 9 | + |
| 10 | +@click.group("phantom") |
| 11 | +def phantom(): |
| 12 | + """Define the simulation phantom (medium + initial pressure).""" |
| 13 | + pass |
| 14 | + |
| 15 | + |
| 16 | +@phantom.command() |
| 17 | +@click.option("--type", "phantom_type", required=True, type=click.Choice(["disc", "spherical", "layered"])) |
| 18 | +@click.option("--grid-size", required=True, help="Grid dimensions, e.g. 128,128") |
| 19 | +@click.option("--spacing", required=True, type=float, help="Grid spacing in meters, e.g. 0.1e-3") |
| 20 | +@click.option("--sound-speed", type=float, default=1500, help="Medium sound speed (m/s)") |
| 21 | +@click.option("--density", type=float, default=1000, help="Medium density (kg/m^3)") |
| 22 | +@click.option("--disc-center", default=None, help="Disc center, e.g. 64,64") |
| 23 | +@click.option("--disc-radius", type=int, default=5, help="Disc radius in grid points") |
| 24 | +@pass_session |
| 25 | +@json_command("phantom.generate") |
| 26 | +def generate(sess, phantom_type, grid_size, spacing, sound_speed, density, disc_center, disc_radius): |
| 27 | + """Generate an analytical phantom.""" |
| 28 | + sess.load() |
| 29 | + |
| 30 | + grid_n = tuple(int(x) for x in grid_size.split(",")) |
| 31 | + ndim = len(grid_n) |
| 32 | + grid_spacing = (spacing,) * ndim |
| 33 | + |
| 34 | + if phantom_type == "disc": |
| 35 | + if ndim != 2: |
| 36 | + raise ValidationError( |
| 37 | + CLIError( |
| 38 | + code="DISC_REQUIRES_2D", |
| 39 | + field="grid_size", |
| 40 | + value=grid_size, |
| 41 | + constraint="disc phantom requires 2D grid", |
| 42 | + suggestion="Use --grid-size Nx,Ny (two dimensions)", |
| 43 | + ) |
| 44 | + ) |
| 45 | + from kwave.data import Vector |
| 46 | + from kwave.utils.mapgen import make_disc |
| 47 | + |
| 48 | + if disc_center is None: |
| 49 | + center = Vector([n // 2 for n in grid_n]) |
| 50 | + else: |
| 51 | + center = Vector([int(x) for x in disc_center.split(",")]) |
| 52 | + |
| 53 | + p0 = make_disc(Vector(list(grid_n)), center, disc_radius).astype(float) |
| 54 | + |
| 55 | + elif phantom_type == "spherical": |
| 56 | + # Simple spherical inclusion centered in grid |
| 57 | + center = np.array([n // 2 for n in grid_n]) |
| 58 | + coords = np.mgrid[tuple(slice(0, n) for n in grid_n)] |
| 59 | + dist = np.sqrt(sum((c - cn) ** 2 for c, cn in zip(coords, center))) |
| 60 | + p0 = (dist <= disc_radius).astype(float) |
| 61 | + |
| 62 | + elif phantom_type == "layered": |
| 63 | + p0 = np.zeros(grid_n) |
| 64 | + layer_pos = grid_n[0] // 4 |
| 65 | + p0[layer_pos, ...] = 1.0 |
| 66 | + |
| 67 | + # Save arrays |
| 68 | + p0_path = sess.save_array("p0", p0) |
| 69 | + |
| 70 | + # Update session |
| 71 | + sess.update( |
| 72 | + "grid", |
| 73 | + { |
| 74 | + "N": list(grid_n), |
| 75 | + "spacing": list(grid_spacing), |
| 76 | + "sound_speed_for_time": sound_speed, |
| 77 | + }, |
| 78 | + ) |
| 79 | + sess.update( |
| 80 | + "medium", |
| 81 | + { |
| 82 | + "sound_speed": sound_speed, |
| 83 | + "density": density, |
| 84 | + }, |
| 85 | + ) |
| 86 | + sess.update( |
| 87 | + "source", |
| 88 | + { |
| 89 | + "type": "initial-pressure", |
| 90 | + "p0_path": p0_path, |
| 91 | + }, |
| 92 | + ) |
| 93 | + |
| 94 | + return CLIResponse( |
| 95 | + result={ |
| 96 | + "phantom_type": phantom_type, |
| 97 | + "grid_size": list(grid_n), |
| 98 | + "spacing": list(grid_spacing), |
| 99 | + "p0_shape": list(p0.shape), |
| 100 | + "p0_max": float(p0.max()), |
| 101 | + "sound_speed": sound_speed, |
| 102 | + "density": density, |
| 103 | + }, |
| 104 | + derived={ |
| 105 | + "ndim": ndim, |
| 106 | + "grid_points": int(np.prod(grid_n)), |
| 107 | + }, |
| 108 | + ) |
0 commit comments