API reference¶
sonosk.Grid¶
- class sonosk.Grid(shape, spacing, origin=0.0)[source]¶
Regular grid.
- Parameters:
shape (int or tuple of int) – Number of points along each axis.
spacing (float or sequence of float) – Point spacing along each axis [m].
origin (float or sequence of float, optional) – Coordinate of the first point on each axis [m]. Default 0.
- classmethod centered(shape, spacing, depth_start=0.0)[source]¶
Grid centred on the depth axis; depth starts at depth_start.
- property extent¶
Length covered along each axis [m].
sonosk.transducer¶
Transducer geometries.
A transducer is a set of elements. Each element is split into small patches; a patch has a centre, an inward normal pointing into the medium and an area. Field and imaging code sum over the patches, so the patch size sets the accuracy. Patches no larger than a quarter wavelength are a safe choice for fields close to the source.
Coordinates are in metres. The array face is near z = 0 and the beam goes towards +z. For arrays, x is the array (azimuth) direction and y the elevation direction.
- class sonosk.transducer.Transducer(pos, normal, area, elem, width=0.0, height=0.0, name='')[source]¶
Collection of elements made of patches.
- Parameters:
pos ((N, 3) array) – Patch centres [m].
normal ((N, 3) array) – Unit normals pointing into the medium.
area ((N,) array) – Patch areas [m^2].
elem ((N,) int array) – Element index of each patch.
width (float) – Element size along x and y [m], used for element directivity.
height (float) – Element size along x and y [m], used for element directivity.
name (str)
- property centers¶
Area-weighted element centres, (E, 3).
- property element_normals¶
Mean unit normal of each element, (E, 3).
- sonosk.transducer.apodization(n, window='hann')[source]¶
Element weights for a window over n elements.
- sonosk.transducer.bowl(radius, roc, dx=0.0001)[source]¶
Spherically focused bowl, apex at the origin, focus at (0, 0, roc).
- sonosk.transducer.annular(radii, roc=None, dx=0.0001)[source]¶
Annular array. radii is a list of (inner, outer) radii, one per ring.
- sonosk.transducer.linear_array(n, pitch, width=None, height=0.005, elev_focus=None, dx=0.0001)[source]¶
Linear array of n elements along x, centred on the origin.
width is the element width (default 0.9 pitch, so kerf 0.1 pitch). elev_focus is the radius of the elevation lens [m]; None for a flat element.
- sonosk.transducer.phased_array(n, pitch, width=None, height=0.01, elev_focus=None, dx=0.0001)[source]¶
Phased array: a linear array with a small pitch, usually half a wavelength.
- sonosk.transducer.convex_array(n, pitch, radius, width=None, height=0.005, elev_focus=None, dx=0.0001)[source]¶
Convex (curved) array on an arc of given radius.
The centre of curvature is at (0, 0, -radius), so the middle of the array touches z = 0. Small radii give endocavity probes (transrectal, transvaginal).
sonosk.field¶
Continuous-wave pressure fields.
rayleigh() sums the Rayleigh integral over the patches of a transducer,
p(r) / p0 = (j k / 2 pi) sum_n w_n A_n exp(-(j k + a) R_n) / R_n,
where p0 = rho c u0, w_n is the complex patch weight (apodization and delay), A_n the patch area, R_n the distance and a the amplitude attenuation of a homogeneous medium. For curved sources this is the usual approximation that treats each patch as part of a baffled plane.
angular_spectrum() propagates a field given on a plane by the angular spectrum method.
Time convention: exp(j omega t).
References
O’Neil H.T. Theory of focusing radiators. J. Acoust. Soc. Am. 21(5):516-526, 1949.
Zeng X., McGough R.J. Evaluation of the angular spectrum approach for simulations of near-field pressures. J. Acoust. Soc. Am. 123(1):68-76, 2008.
- sonosk.field.rayleigh(transducer, points, f, c=1540.0, alpha0=0.0, y=1.0, max_pairs=4000000)[source]¶
Complex pressure p/p0 at points.
- Parameters:
transducer (Transducer)
points ((M, 3) array) – Field points [m].
f (float) – Frequency [Hz].
c (float) – Sound speed [m/s].
alpha0 (float) – Power-law attenuation of the medium, alpha0 in dB/(cm MHz^y).
y (float) – Power-law attenuation of the medium, alpha0 in dB/(cm MHz^y).
max_pairs (int) – Number of point-patch pairs handled at once; lower it if memory is short.
- sonosk.field.grid_points(grid, y=0.0)[source]¶
Field points for a grid: 1-D on the z axis, 2-D in the x-z plane at y, 3-D everywhere. Returns an (N, 3) array.
- sonosk.field.rayleigh_grid(transducer, grid, f, c=1540.0, alpha0=0.0, y=1.0, plane_y=0.0, **kw)[source]¶
rayleigh() evaluated on a Grid; result has the grid shape.
- sonosk.field.angular_spectrum(p0, dx, f, z, c=1540.0, alpha0=0.0, y=1.0, pad=2)[source]¶
Propagate a complex field from the plane z = 0 to the depths z.
- Parameters:
p0 (1-D or 2-D complex array) – Field on a line (x) or a plane (x, y). A 1-D input is a field that does not change along y, so the result is a 2-D (x, z) problem.
dx (float) – Sample spacing of p0 [m], the same along both axes.
f (float) – Frequency [Hz].
z (array) – Depths [m], measured from the input plane.
pad (int) – Zero padding factor against wrap-around.
- Return type:
array of shape (len(z),) + p0.shape
sonosk.medium¶
Materials and property maps.
Attenuation follows a power law, alpha(f) = alpha0 f^y, with alpha0 in the clinical unit dB/(cm MHz^y). Inside the solvers it is converted to an amplitude coefficient in Np/m.
- sonosk.medium.attenuation_np(alpha0, f, y=1.0)[source]¶
Amplitude attenuation [Np/m] at frequency f [Hz].
alpha0 is in dB/(cm MHz^y).
- class sonosk.medium.Material(name: str, c: float, rho: float, alpha0: float = 0.0, y: float = 1.0, ba: float = 0.0, k: float = 0.0, cp: float = 0.0, perfusion: float = 0.0, echo: float = 0.0, source: str = '')[source]¶
Acoustic and thermal properties of one material.
c [m/s], rho [kg/m^3], alpha0 [dB/(cm MHz^y)], y power-law exponent, ba nonlinearity B/A, k thermal conductivity [W/(m K)], cp specific heat [J/(kg K)], perfusion blood mass flow per unit volume [kg/(m^3 s)], echo backscatter level relative to a reference tissue [dB], source a reference for the values.
sonosk.phantom¶
Building blocks for numerical phantoms.
Shapes return boolean masks on a Grid. Label maps are integer arrays where each value selects a Material (see sonosk.medium). Point scatterers for imaging are drawn from a label map, with an amplitude set by the echo level of each material.
The depth axis is the last grid axis.
- sonosk.phantom.layers(grid, thickness, start=0.0)[source]¶
Label map of planar layers stacked in depth.
thickness is a list of layer thicknesses [m]. Points above start get label 0, the n-th layer gets label n, and everything below the last layer gets label len(thickness) + 1. Use paint() to change labels.
- sonosk.phantom.sphere(grid, center, radius)[source]¶
Sphere in 3-D, disc in the x-z plane in 2-D, segment in 1-D.
center is (x, y, z) [m]; y is ignored in 2-D and x, y in 1-D.
- sonosk.phantom.ellipsoid(grid, center, semi_axes, angle=0.0)[source]¶
Ellipsoid with semi-axes (ax, ay, az), rotated by angle [rad] about y.
- sonosk.phantom.cylinder(grid, point, direction, radius)[source]¶
Infinite cylinder (a vessel) through point along direction.
- sonosk.phantom.box(grid, low, high)[source]¶
Axis-aligned box between corners low and high, given as (x, y, z).
- sonosk.phantom.paint(labels, mask, label)[source]¶
Copy of labels with label written where mask is true.
- sonosk.phantom.scatterers(medium, density, rng=None, region=None)[source]¶
Random point scatterers for pulse-echo simulation.
- Parameters:
medium (Medium) – Its grid sets the region; the echo level of each material sets the mean scattering amplitude, 10**(echo/20).
density (float) – Scatterers per mm^2 for a 2-D grid, per mm^3 for a 3-D grid. Ten or more per resolution cell give fully developed speckle.
rng (numpy Generator or int, optional)
region (bool array, optional) – Only place scatterers where this mask is true.
- Returns:
pos ((N, 3) array [m])
amp ((N,) array)
sonosk.imaging¶
Pulse-echo simulation and image formation.
The forward model is linear (first-order Born scattering) in a homogeneous medium with sound speed c:
the transmit wave is described by a virtual source: a plane wave with a steering angle, a diverging wave from a point behind the array, or a wave focused on a point in front of it;
each receive element is a point receiver with the directivity of a rectangular element in a soft baffle, sinc(w sin(theta) / lambda) cos(theta);
spreading is 1/R on the receive path and 1/R from a diverging source; attenuation is applied at the centre frequency along both paths.
The echo of every scatterer is a copy of the pulse delayed by its travel time. Signals are kept in analytic (complex) form, so the envelope is the modulus of the beamformed image.
Beamforming is delay-and-sum with a receive f-number. Images from several transmits can be summed (coherent compounding).
References
Montaldo G., Tanter M., Bercoff J., Benech N., Fink M. Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 56(3):489-506, 2009.
Garcia D. SIMUS: an open-source simulator for medical ultrasound imaging. Part I: theory & examples. Comput. Methods Programs Biomed. 218:106726, 2022.
- class sonosk.imaging.Pulse(f0: float, bandwidth: float = 0.6)[source]¶
Gaussian-modulated pulse.
f0 centre frequency [Hz], bandwidth fractional -6 dB bandwidth of the spectrum (0.6 means 60 %).
- property duration¶
Support used for the simulation, 8 sigma [s].
- class sonosk.imaging.Transmit(kind, angle=0.0, source=None)[source]¶
Transmit wave described by a virtual source.
Use the constructors plane(), diverging() and focused().
- sonosk.imaging.simulate(array, transmit, pos, amp, pulse, c=1540.0, fs=None, alpha0=0.0, depth=None, oversample=4, element_chunk=16)[source]¶
Simulate the received signals of one transmit.
Every echo is placed as an impulse on a time axis oversample times finer than fs, split between the two nearest samples, then convolved with the pulse and decimated to fs.
- Parameters:
array (Transducer) – Linear, phased, convex or matrix array.
transmit (Transmit)
pos (arrays) – Scatterer positions (N, 3) [m] and amplitudes (N,).
amp (arrays) – Scatterer positions (N, 3) [m] and amplitudes (N,).
pulse (Pulse) – Two-way pulse (the electrical excitation convolved with both transducer responses).
c (float) – Sound speed [m/s].
fs (float, optional) – Sampling frequency of the output [Hz], default 8 f0.
alpha0 (float) – Attenuation in dB/(cm MHz), applied at f0.
depth (float, optional) – Record at least down to this depth [m].
oversample (int) – Ratio of the internal sampling rate to fs.
- Returns:
rf (complex array (n_elements, n_samples)) – Analytic received signals; rf.real is the RF signal.
fs (float)
- sonosk.imaging.beamform(rf, fs, array, transmit, x, z, c=1540.0, f_number=1.5, window='hann', y=0.0)[source]¶
Delay-and-sum image on the grid x by z (2-D, at elevation y).
Returns the complex image of shape (len(z), len(x)); its modulus is the envelope.
- sonosk.imaging.compound(rfs, fs, array, transmits, x, z, c=1540.0, **kw)[source]¶
Sum of beamformed images over several transmits.
sonosk.thermal¶
Heating by ultrasound, bioheat transfer and thermal dose.
The temperature T follows the Pennes bioheat equation,
rho cp dT/dt = div(k grad T) + w_b c_b (T_a - T) + Q,
with w_b the blood perfusion [kg/(m^3 s)], c_b the specific heat of blood and T_a the arterial temperature. The heat source of a CW field with pressure amplitude p is Q = 2 alpha I = alpha p^2 / (rho c), where alpha is the amplitude absorption coefficient [Np/m].
The equation is solved with an explicit finite-difference scheme on the grid of the medium, in 1, 2 or 3 dimensions, with the outer boundary held at the body temperature. The time step is chosen inside the stability limit.
Thermal dose is the cumulative equivalent minutes at 43 C,
CEM43 = sum dt R^(43 - T), R = 0.5 for T >= 43 C, 0.25 below,
with dt in minutes. A dose of 240 min is the usual threshold for thermal necrosis in soft tissue.
References
Pennes H.H. Analysis of tissue and arterial blood temperatures in the resting human forearm. J. Appl. Physiol. 1(2):93-122, 1948.
Sapareto S.A., Dewey W.C. Thermal dose determination in cancer therapy. Int. J. Radiat. Oncol. Biol. Phys. 10(6):787-800, 1984.
- sonosk.thermal.heat_source(p, alpha, rho, c)[source]¶
Absorbed power density [W/m^3] for pressure amplitude p [Pa].
alpha is the absorption coefficient [Np/m]; with no separate absorption data the attenuation is often used, which overestimates heating where scattering is large.
- sonosk.thermal.cem43(T, dt)[source]¶
Thermal dose [min] of a temperature history.
T has time on its first axis, sampled every dt seconds.
- class sonosk.thermal.BioheatResult(t: numpy.ndarray, T: numpy.ndarray, T_max: numpy.ndarray, dose: numpy.ndarray, dt: float, probe: numpy.ndarray = None)[source]¶
- sonosk.thermal.bioheat(medium, q, duration, schedule=None, T0=37.0, T_blood=37.0, c_blood=3617.0, dt=None, max_dt=0.1, snapshots=None, probe=None)[source]¶
Solve the Pennes equation on the grid of medium.
- Parameters:
medium (Medium) – Uses rho, cp, k and perfusion.
q (array) – Heat source [W/m^3] with the grid shape, at full power.
duration (float) – Simulated time [s].
schedule (callable, optional) – schedule(t) gives the power fraction between 0 and 1 at time t, for pulsed or interrupted exposures. Default: on all the time.
T0 (float) – Initial and arterial temperature [C]. The grid boundary stays at T0.
T_blood (float) – Initial and arterial temperature [C]. The grid boundary stays at T0.
dt (float, optional) – Time step [s]; by default 0.9 of the stability limit, at most max_dt. The schedule is sampled once per step, so dt must be short against the on and off times of a pulsed exposure.
snapshots (array, optional) – Times [s] at which to store the temperature field.
probe (list of index tuples, optional) – Grid points whose temperature is stored at every step.
- sonosk.thermal.treatment_report(lesion, medium, target)[source]¶
Volumes treated and spared.
- Parameters:
lesion (bool array) – Points above the dose threshold.
medium (Medium)
target (str or list of str) – Names of the target materials (for example the tumour).
- Returns:
dict with the target volume, ablated and residual target volume,
fraction of the target ablated, damaged volume outside the target and
the damaged volume of each material. Volumes in m^3 (3-D), m^2 (2-D)
or m (1-D).
sonosk.stack¶
One-dimensional model of a thickness-mode piezoelectric transducer.
The piezoelectric plate is described by the Mason equivalent circuit. Backing and load are semi-infinite media given by their acoustic impedance. Matching layers and a lens are transmission lines. The model gives the electrical input impedance, the transmit response and the shape of the pulse-echo response of a stack, which is the first step in designing a probe element.
Impedances of media are specific acoustic impedances in Rayl (kg m^-2 s^-1). Forces and velocities on the plate are per unit area.
References
Mason W.P. Electromechanical Transducers and Wave Filters, 2nd ed. Van Nostrand, 1948.
Kino G.S. Acoustic Waves: Devices, Imaging, and Analog Signal Processing. Prentice-Hall, 1987.
Desilets C.S., Fraser J.D., Kino G.S. The design of efficient broad-band piezoelectric transducers. IEEE Trans. Sonics Ultrason. 25(3):115-125, 1978.
- class sonosk.stack.Layer(thickness: float, c: float, rho: float, alpha0: float = 0.0, y: float = 1.0)[source]¶
Passive layer: thickness [m], sound speed [m/s], density [kg/m^3], attenuation alpha0 [dB/(cm MHz^y)] and exponent y.
- class sonosk.stack.Piezo(thickness: float, c: float, rho: float, kt: float, eps_r: float, area: float)[source]¶
Piezoelectric plate in thickness mode.
thickness [m], c stiffened longitudinal sound speed [m/s], rho [kg/m^3], kt thickness coupling factor, eps_r clamped relative permittivity eps33^S / eps0, area electrode area [m^2].
- property c0¶
Clamped capacitance [F].
- property f_half_wave¶
Frequency at which the plate is half a wavelength thick [Hz].
- class sonosk.stack.Stack(piezo: Piezo, backing: float = 400.0, layers: list = <factory>, load: float = 1500000.0)[source]¶
Piezo plate with backing, front layers (from the plate outwards) and load.
- transmit(f, r_source=50.0)[source]¶
Pressure in the load per volt of the source [Pa/V].
The source has internal resistance r_source [ohm].
sonosk.safety¶
Acoustic output quantities and regulatory limits for diagnostic ultrasound.
Derated (“.3”) values use the 0.3 dB/(cm MHz) model of IEC 62359 and the FDA guidance. Inputs are water values on the beam axis.
References
IEC 62359:2010+A1:2017. Ultrasonics - Field characterization - Test methods for the determination of thermal and mechanical indices related to medical diagnostic ultrasonic fields.
FDA. Marketing Clearance of Diagnostic Ultrasound Systems and Transducers. Guidance for Industry and FDA Staff, February 2023 (FDA-2017-D-5372).
- sonosk.safety.derate(p, f_mhz, z_cm, alpha=0.3)[source]¶
Derated pressure for frequency f_mhz [MHz] at depth z_cm [cm].
- sonosk.safety.derate_intensity(i, f_mhz, z_cm, alpha=0.3)[source]¶
Derated intensity (or PII) at depth z_cm [cm].
- sonosk.safety.pulse_intensity_integral(p, dt, rho=1000.0, c=1500.0)[source]¶
Pulse intensity integral [J/m^2] of pressure waveform p [Pa].
The time axis is the last axis of p; dt is the sample interval [s].
- sonosk.safety.pulse_duration(p, dt)[source]¶
Pulse duration [s]: 1.25 times the 10 % to 90 % interval of the running pulse intensity integral (IEC 62359).
- sonosk.safety.isppa(p, dt, rho=1000.0, c=1500.0)[source]¶
Pulse-average intensity [W/m^2] of a single waveform.
- sonosk.safety.mechanical_index(z_cm, pr_mpa, f_mhz, pii=None)[source]¶
Mechanical index.
MI = pr.3(z_sp) / sqrt(fc), where z_sp is the point on the beam axis where the derated pulse intensity integral is largest.
- Parameters:
z_cm (array) – Depths on the beam axis [cm].
pr_mpa (array) – Peak rarefactional pressure in water at z_cm [MPa], positive.
f_mhz (float) – Centre (acoustic working) frequency [MHz].
pii (array, optional) – Pulse intensity integral in water at z_cm. If omitted, pr_mpa**2 is used, which assumes the pulse shape does not change with depth.
- Returns:
mi (float)
z_sp (float) – Depth used [cm].
- sonosk.safety.check_limits(mi=None, ispta3=None, isppa3=None, ti=None, application='general')[source]¶
Compare values with the FDA Track 3 limits.
Intensities are given in clinical units: ispta3 in mW/cm^2, isppa3 in W/cm^2. For general use either MI or ISPPA.3 must be within its limit.
Returns a dict name -> (value, limit, within_limit) and an overall flag under the key “ok”.
sonosk.quality¶
Image quality measures.
Inputs are envelope (amplitude) values, not log-compressed images.
References
Burckhardt C.B. Speckle in ultrasound B-mode scans. IEEE Trans. Sonics Ultrason. 25(1):1-6, 1978.
Rodriguez-Molares A. et al. The generalized contrast-to-noise ratio: a formal definition for lesion detectability. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 67(4):745-759, 2020.
- sonosk.quality.width(profile, x, level_db=-6.0)[source]¶
Width of the main lobe of profile at level_db below its peak.
profile is an amplitude profile sampled at positions x. Edges are found by linear interpolation. Use level_db=-6 for the usual -6 dB resolution and -20 dB for the beam width at that level.
sonosk.analytic¶
Closed-form on-axis pressure of circular sources in a lossless fluid.
These solutions are used as references for the numerical field code. Results are pressure amplitudes normalised by p0 = rho c u0, where u0 is the normal velocity amplitude of the source surface.
References
Kinsler L.E., Frey A.R., Coppens A.B., Sanders J.V. Fundamentals of Acoustics, 4th ed. Wiley, 2000.
O’Neil H.T. Theory of focusing radiators. J. Acoust. Soc. Am. 21(5):516-526, 1949.
- sonosk.analytic.piston_axial(z, a, f, c=1500.0)[source]¶
Pressure amplitude p/p0 on the axis of a baffled flat piston.
z distance from the piston face [m], a piston radius [m], f frequency [Hz], c sound speed [m/s].
- sonosk.analytic.piston_last_maximum(a, f, c=1500.0)[source]¶
Position of the last on-axis maximum of a flat piston [m].
- sonosk.analytic.bowl_depth(a, roc)[source]¶
Depth h of a spherical bowl of aperture radius a and curvature roc.
sonosk.water¶
Properties of pure water.
References
Marczak W. Water as a standard in the measurements of speed of sound in liquids. J. Acoust. Soc. Am. 102(5):2776-2779, 1997.
Kell G.S. Density, thermal expansivity, and compressibility of liquid water from 0 to 150 C. J. Chem. Eng. Data 20(1):97-105, 1975.
sonosk.plot¶
Figures for fields, images and thermal results. Needs matplotlib.
- sonosk.plot.bmode(image, x, z, dynamic_range=60.0, ax=None, title=None)[source]¶
Log-compressed B-mode image. image is complex or envelope, (len(z), len(x)).
- sonosk.plot.field_map(p, x, z, db=True, dynamic_range=40.0, ax=None, title=None, cmap='inferno')[source]¶
Pressure amplitude map in the x-z plane. p has shape (len(x), len(z)).
- sonosk.plot.temperature_map(T, x, z, ax=None, title=None, lesion=None, vmin=37.0, vmax=None)[source]¶
Temperature in the x-z plane, optionally with the lesion outline. T and lesion have shape (len(x), len(z)).