Coverage for kwave/utils/conversion.py: 12%
214 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-03-24 12:06 -0700
1import logging
2import math
3from typing import Any
5import numpy as np
6from beartype import beartype as typechecker
7from beartype.typing import Tuple, Union
8from jaxtyping import Float, Num, Real
9from numpy import ndarray
11import kwave.utils.typing as kt
12from kwave.kgrid import kWaveGrid
13from kwave.utils.matlab import matlab_mask
14from kwave.utils.matrix import sort_rows
17@typechecker
18def db2neper(alpha: Real[kt.ArrayLike, "..."], y: Real[kt.ScalarLike, ""] = 1) -> Real[kt.ArrayLike, "..."]:
19 """
20 Convert decibels to nepers.
22 Args:
23 alpha: Attenuation in dB / (MHz ^ y cm).
24 y: Power law exponent (default=1).
26 Returns:
27 Attenuation in Nepers / ((rad / s) ^ y m).
29 """
31 # calculate conversion
32 alpha = 100.0 * alpha * (1e-6 / (2.0 * math.pi)) ** y / (20.0 * np.log10(np.exp(1)))
33 return alpha
36@typechecker
37def neper2db(alpha: Real[kt.ArrayLike, "..."], y: Real[kt.ScalarLike, ""] = 1) -> Real[kt.ArrayLike, "..."]:
38 """
39 Converts an attenuation coefficient in units of Nepers / ((rad / s) ^ y m) to units of dB / (MHz ^ y cm).
41 Args:
42 alpha: Attenuation in Nepers / ((rad / s) ^ y m)
43 y: Power law exponent (default=1)
45 Returns:
46 Attenuation in dB / (MHz ^ y cm)
48 """
50 # calculate conversion
51 alpha = 20 * math.log10(math.exp(1)) * alpha * (2 * math.pi * 1e6) ** y / 100
52 return alpha
55@typechecker
56def cast_to_type(data: Real[kt.ArrayLike, "..."], matlab_type: str) -> Any:
57 """
59 Args:
60 data: The data to cast.
61 matlab_type: The type to cast to.
63 Returns:
64 The cast data.
66 """
67 if not isinstance(data, np.ndarray):
68 data = np.array(data)
69 type_map = {
70 "single": np.float32,
71 "double": np.float64,
72 "uint64": np.uint64,
73 "uint32": np.uint32,
74 "uint16": np.uint16,
75 }
76 return data.astype(type_map[matlab_type])
79@typechecker
80def cart2pol(x: Real[kt.ArrayLike, "..."], y: Real[kt.ArrayLike, "..."]) -> Tuple[Real[kt.ArrayLike, "..."], Real[kt.ArrayLike, "..."]]:
81 """
82 Convert from cartesian to polar coordinates.
84 Args:
85 x: The x-coordinate of the point.
86 y: The y-coordinate of the point.
88 Returns:
89 A tuple containing the polar coordinates of the point.
91 """
93 rho = np.sqrt(x**2 + y**2)
94 phi = np.arctan2(y, x)
95 return phi, rho
98@typechecker
99def grid2cart(input_kgrid: kWaveGrid, grid_selection: ndarray) -> Tuple[ndarray, ndarray]:
100 """
101 Returns the Cartesian coordinates of the non-zero points of a binary grid.
103 Args:
104 input_kgrid: k-Wave grid object returned by kWaveGrid
105 grid_selection: binary grid with the same dimensions as the k-Wave grid kgrid
107 Returns:
108 cart_data: 1 x N, 2 x N, or 3 x N (for 1, 2, and 3 dimensions) array of Cartesian sensor points
109 order_index: returns a list of indices of the returned cart_data coordinates.
111 Raises:
112 ValueError: when input_kgrid.dim is not in [1, 2, 3]
114 """
116 grid_data = np.array((grid_selection != 0), dtype=bool)
117 cart_data = np.zeros((input_kgrid.dim, np.sum(grid_data)))
119 if input_kgrid.dim > 0:
120 cart_data[0, :] = input_kgrid.x[grid_data]
121 if input_kgrid.dim > 1:
122 cart_data[1, :] = input_kgrid.y[grid_data]
123 if input_kgrid.dim > 2:
124 cart_data[2, :] = input_kgrid.z[grid_data]
125 if 0 <= input_kgrid.dim > 3:
126 raise ValueError("kGrid with unsupported size passed.")
128 order_index = np.argwhere(grid_data.squeeze() != 0)
129 return cart_data.squeeze(), order_index
132@typechecker
133def freq2wavenumber(n: int, k_max: float, filter_cutoff: float, c: float, k_dim: Union[int, Tuple[int]]) -> Tuple[int, float]:
134 """
135 Convert the given frequency and maximum wavenumber to a wavenumber cutoff and filter size.
137 Args:
138 n: The size of the grid.
139 k_max: The maximum wavenumber.
140 filter_cutoff: The frequency to convert to a wavenumber cutoff.
141 c: The speed of sound.
142 k_dim: The dimensions of the wavenumber grid.
144 Returns:
145 A tuple containing the calculated filter size and wavenumber cutoff.
147 """
149 k_cutoff = 2 * np.pi * filter_cutoff / c
151 # set the alpha_filter size
152 filter_size = round(n * k_cutoff / k_dim[-1])
154 # check the alpha_filter size
155 if filter_size > n:
156 # set the alpha_filter size to be the same as the grid size
157 filter_size = n
158 filter_cutoff = k_max * c / (2 * np.pi)
159 return filter_size, filter_cutoff
162@typechecker
163def cart2grid(
164 kgrid: kWaveGrid,
165 cart_data: Union[Float[ndarray, "1 NumPoints"], Float[ndarray, "2 NumPoints"], Float[ndarray, "3 NumPoints"]],
166 axisymmetric: bool = False,
167) -> Tuple:
168 """
169 Interpolates the set of Cartesian points defined by
170 cart_data onto a binary matrix defined by the kWaveGrid object
171 kgrid using nearest neighbour interpolation. An error is returned if
172 the Cartesian points are outside the computational domain defined by
173 kgrid.
175 Args:
176 kgrid: simulation grid
177 cart_data: Cartesian sensor points
178 axisymmetric: set to True to use axisymmetric interpolation
180 Returns:
181 A binary grid
183 """
185 # check for axisymmetric input
186 if axisymmetric and kgrid.dim != 2:
187 raise AssertionError("Axisymmetric flag only supported in 2D.")
189 # detect whether the inputs are for one, two, or three dimensions
190 if kgrid.dim == 1:
191 # one-dimensional
192 data_x = cart_data[0, :]
194 # scale position values to grid centered pixel coordinates using
195 # nearest neighbour interpolation
196 data_x = np.round(data_x / kgrid.dx).astype(int)
198 # shift pixel coordinates to coincide with matrix indexing
199 data_x = data_x + np.floor(kgrid.Nx // 2).astype(int)
201 # check if the points all lie within the grid
202 if data_x.max() > kgrid.Nx or data_x.min() < 1:
203 raise AssertionError("Cartesian points must lie within the grid defined by kgrid.")
205 # create empty grid
206 grid_data = np.zeros((kgrid.Nx, 1))
208 # create index variable
209 point_index = np.arange(1, data_x.size + 1)
211 # map values
212 for data_index in range(data_x.size):
213 grid_data[data_x[data_index]] = point_index[data_index]
215 # extract reordering index
216 reorder_index = np.reshape(grid_data[grid_data != 0], (-1, 1))
218 elif kgrid.dim == 2:
219 # two-dimensional
220 data_x = cart_data[0, :]
221 data_y = cart_data[1, :]
223 # scale position values to grid centered pixel coordinates using
224 # nearest neighbour interpolation
225 data_x = np.round(data_x / kgrid.dx).astype(int)
226 data_y = np.round(data_y / kgrid.dy).astype(int)
228 # shift pixel coordinates to coincide with matrix indexing (leave
229 # y-direction = radial-direction if axisymmetric)
230 data_x = data_x + np.floor(kgrid.Nx // 2).astype(int)
231 if not axisymmetric:
232 data_y = data_y + np.floor(kgrid.Ny // 2).astype(int)
234 # check if the points all lie within the grid
235 if data_x.max() >= kgrid.Nx or data_y.max() >= kgrid.Ny or data_x.min() < 0 or data_y.min() < 0:
236 raise AssertionError("Cartesian points must lie within the grid defined by kgrid.")
238 # create empty grid
239 grid_data = -1 * np.ones((kgrid.Nx, kgrid.Ny))
241 # map values
242 for data_index in range(data_x.size):
243 grid_data[data_x[data_index], data_y[data_index]] = int(data_index)
245 # extract reordering index
246 reorder_index = grid_data.flatten(order="F")[grid_data.flatten(order="F") != -1]
247 reorder_index = reorder_index[:, None] + 1 # [N] => [N, 1]
249 elif kgrid.dim == 3:
250 # three dimensional
251 data_x = cart_data[0, :]
252 data_y = cart_data[1, :]
253 data_z = cart_data[2, :]
255 # scale position values to grid centered pixel coordinates using
256 # nearest neighbour interpolation
257 data_x = np.round(data_x / kgrid.dx).astype(int)
258 data_y = np.round(data_y / kgrid.dy).astype(int)
259 data_z = np.round(data_z / kgrid.dz).astype(int)
261 # shift pixel coordinates to coincide with matrix indexing
262 data_x = data_x + np.floor(kgrid.Nx // 2).astype(int)
263 data_y = data_y + np.floor(kgrid.Ny // 2).astype(int)
264 data_z = data_z + np.floor(kgrid.Nz // 2).astype(int)
266 # check if the points all lie within the grid
267 assert (
268 0 <= data_x.min()
269 and 0 <= data_y.min()
270 and 0 <= data_z.min()
271 and data_x.max() < kgrid.Nx
272 and data_y.max() < kgrid.Ny
273 and data_z.max() < kgrid.Nz
274 ), "Cartesian points must lie within the grid defined by kgrid."
276 # create empty grid
277 grid_data = -1 * np.ones((kgrid.Nx, kgrid.Ny, kgrid.Nz), dtype=int)
279 # create index variable
280 point_index = np.arange(1, data_x.size + 1)
282 # map values
283 for data_index in range(data_x.size):
284 grid_data[data_x[data_index], data_y[data_index], data_z[data_index]] = point_index[data_index]
286 # extract reordering index
287 reorder_index = grid_data.flatten(order="F")[grid_data.flatten(order="F") != -1]
288 reorder_index = reorder_index[:, None, None] # [N] => [N, 1, 1]
289 else:
290 raise ValueError("Input cart_data must be a 1, 2, or 3 dimensional matrix.")
292 # compute the reverse ordering index (i.e., what is the index of each point
293 # in the reordering vector)
294 order_index = np.ones((reorder_index.size, 2), dtype=int)
295 order_index[:, 0] = np.squeeze(reorder_index)
296 order_index[:, 1] = np.arange(1, reorder_index.size + 1)
297 order_index = sort_rows(order_index, 0)
298 order_index = order_index[:, 1]
299 order_index = order_index[:, None] # [N] => [N, 1]
301 # reset binary grid values
302 if kgrid.dim == 1:
303 grid_data[grid_data != 0] = 1
304 else:
305 grid_data[grid_data != -1] = 1
306 grid_data[grid_data == -1] = 0
308 # check if any Cartesian points have been mapped to the same grid point,
309 # thereby reducing the total number of points
310 num_discarded_points = cart_data.shape[1] - np.sum(grid_data)
311 if num_discarded_points != 0:
312 logging.log(logging.INFO, f" cart2grid: {num_discarded_points} Cartesian points mapped to overlapping grid points")
313 return grid_data.astype(int), order_index, reorder_index
316@typechecker
317def hounsfield2soundspeed(
318 ct_data: Union[Float[ndarray, "Dim1 Dim2"], Float[ndarray, "Dim1 Dim2 Dim3"]],
319) -> Union[Float[ndarray, "Dim1 Dim2"], Float[ndarray, "Dim1 Dim2 Dim3"]]:
320 """
321 Calculates the sound speed of a medium given a CT (computed tomography) of the medium.
322 For soft tissue, the approximate sound speed can also be returned using the empirical relationship
323 given by Mast [1].
325 Args:
326 ct_data: matrix of Hounsfield values.
328 Returns:
329 A matrix of sound speed values of size of ct_data.
331 References:
332 [1] Mast, T. D., "Empirical relationships between acoustic parameters in human soft tissues,"
333 Acoust. Res. Lett. Online, 1(2), pp. 37-42 (2000).
334 """
335 # calculate corresponding sound speed values if required using soft tissue relationship
336 # TODO confirm that this linear relationship is correct
337 sound_speed = (hounsfield2density(ct_data) + 349) / 0.893
339 return sound_speed
342@typechecker
343def hounsfield2density(
344 ct_data: Union[Float[ndarray, "Dim1 Dim2"], Float[ndarray, "Dim1 Dim2 Dim3"]], plot_fitting: bool = False
345) -> Union[Float[ndarray, "Dim1 Dim2"], Float[ndarray, "Dim1 Dim2 Dim3"]]:
346 """
347 Convert Hounsfield units in CT data to density values [kg / m ^ 3] based on experimental data.
349 Args:
350 ct_data: The CT data in Hounsfield units.
351 plot_fitting (bool, optional): Whether to plot the fitting curve (default: False).
353 Returns:
354 The density values in [kg / m ^ 3].
356 """
358 # create empty density matrix
359 density = np.zeros(ct_data.shape, like=ct_data)
361 # apply conversion in several parts using linear fits to the data
362 # Part 1: Less than 930 Hounsfield Units
363 density[ct_data < 930] = np.polyval([1.025793065681423, -5.680404011488714], ct_data[ct_data < 930])
365 # Part 2: Between 930 and 1098(soft tissue region)
366 index_selection = np.logical_and(930 <= ct_data, ct_data <= 1098)
367 density[index_selection] = np.polyval([0.9082709691264, 103.6151457847139], ct_data[index_selection])
369 # Part 3: Between 1098 and 1260(between soft tissue and bone)
370 index_selection = np.logical_and(1098 < ct_data, ct_data < 1260)
371 density[index_selection] = np.polyval([0.5108369316599, 539.9977189228704], ct_data[index_selection])
373 # Part 4: Greater than 1260(bone region)
374 density[ct_data >= 1260] = np.polyval([0.6625370912451, 348.8555178455294], ct_data[ct_data >= 1260])
376 if plot_fitting:
377 raise NotImplementedError("Plotting function not implemented in Python")
379 return density
382tol = None
383subs0 = None
386@typechecker
387def tol_star(
388 tolerance: kt.NUMERIC,
389 kgrid: kWaveGrid,
390 point: Union[Float[ndarray, "1"], Float[ndarray, "2"], Float[ndarray, "3"]],
391 debug,
392) -> Tuple[ndarray, ndarray, ndarray, ndarray]:
393 global tol, subs0
395 ongrid_threshold = kgrid.dx * 1e-3
397 kgrid_dim = kgrid.dim
398 if tol is None or tolerance != tol or len(subs0) != kgrid_dim:
399 tol = tolerance
401 decay_subs = int(np.ceil(1 / (np.pi * tol)))
403 lin_ind = np.arange(-decay_subs, decay_subs + 1)
405 if kgrid_dim == 1:
406 is0 = lin_ind
407 elif kgrid_dim == 2:
408 is0, js0 = np.meshgrid(lin_ind, lin_ind, indexing="ij")
409 elif kgrid_dim == 3:
410 is0, js0, ks0 = np.meshgrid(lin_ind, lin_ind, lin_ind, indexing="ij")
412 if kgrid_dim == 1:
413 subs0 = [is0]
414 elif kgrid_dim == 2:
415 instar = np.abs(is0 * js0) <= decay_subs
416 is0 = matlab_mask(is0, instar)
417 js0 = matlab_mask(js0, instar)
418 subs0 = [is0, js0]
419 elif kgrid_dim == 3:
420 instar = np.abs(is0 * js0 * ks0) <= decay_subs
421 is0 = matlab_mask(is0, instar).T
422 js0 = matlab_mask(js0, instar).T
423 ks0 = matlab_mask(ks0, instar).T
424 subs0 = [is0, js0, ks0]
426 is_ = subs0[0].copy()
427 js = subs0[1].copy() if kgrid_dim > 1 else []
428 ks = subs0[2].copy() if kgrid_dim > 2 else []
430 # returns python index value (0 start) not MATLAB index (1 start)
431 x_closest, x_closest_ind = find_closest(kgrid.x_vec, point[0])
433 if np.abs(x_closest - point[0]) < ongrid_threshold:
434 if kgrid_dim > 1:
435 js = js[is_ == 0]
436 if kgrid_dim > 2:
437 ks = ks[is_ == 0]
438 is_ = is_[is_ == 0]
440 if kgrid_dim > 1:
441 y_closest, y_closest_ind = find_closest(kgrid.y_vec, point[1])
442 if np.abs(y_closest - point[1]) < ongrid_threshold:
443 is_ = is_[js == 0]
444 if kgrid_dim > 2:
445 ks = ks[js == 0]
446 js = js[js == 0]
448 if kgrid_dim > 2:
449 z_closest, z_closest_ind = find_closest(kgrid.z_vec, point[2])
450 if np.abs(z_closest - point[2]) < ongrid_threshold:
451 is_ = is_[ks == 0]
452 js = js[ks == 0]
453 ks = ks[ks == 0]
455 # TODO: closest index is added to is_, +1 might not be correct
456 is_ += x_closest_ind + 1
457 if kgrid_dim > 1:
458 js += y_closest_ind + 1
459 if kgrid_dim > 2:
460 ks += z_closest_ind + 1
462 if kgrid_dim == 1:
463 inbounds = (1 <= is_) & (is_ <= kgrid.Nx)
464 # TODO: this should likely be matlabmask and indexes should be checked.
465 is_ = is_[inbounds]
466 elif kgrid_dim == 2:
467 inbounds = (1 <= is_) & (is_ <= kgrid.Nx) & (1 <= js) & (js <= kgrid.Ny)
468 is_ = is_[inbounds]
469 js = js[inbounds]
470 if kgrid_dim == 3:
471 inbounds = (1 <= is_) & (is_ <= kgrid.Nx) & (1 <= js) & (js <= kgrid.Ny) & (1 <= ks) & (ks <= kgrid.Nz)
472 is_ = is_[inbounds]
473 js = js[inbounds]
474 ks = ks[inbounds]
476 if kgrid_dim == 1:
477 lin_ind = is_
478 elif kgrid_dim == 2:
479 lin_ind = kgrid.Nx * (js - 1) + is_
480 elif kgrid_dim == 3:
481 lin_ind = kgrid.Nx * kgrid.Ny * (ks - 1) + kgrid.Nx * (js - 1) + is_
483 # TODO: in current form linear indexing matches MATLAB, but might break in 0 indexed python
484 return lin_ind, np.array(is_) - 1, np.array(js) - 1, np.array(ks) - 1 # -1 for mapping from Matlab indexing to Python indexing
487@typechecker
488def find_closest(array: ndarray, value: Num[kt.ScalarLike, ""]):
489 array = np.asarray(array)
490 idx = (np.abs(array - value)).argmin()
491 return array[idx], idx
494def create_index_at_dim(ndim: int, dim: int, index_value: Any) -> tuple:
495 """
496 Create a tuple of slice objects with a specific index value at the specified dimension.
498 Args:
499 ndim: Number of dimensions in the array
500 dim: The dimension where the specific index should be placed
501 index_value: The index value to place at the specified dimension
503 Returns:
504 A tuple of slice objects with the index value at the specified dimension
505 """
506 return tuple(index_value if i == dim else slice(None) for i in range(ndim))