Coverage for kwave/utils/mapgen.py: 4%

1import logging

2import math

3import warnings

4from math import floor

6import matplotlib.pyplot as plt

7import numpy as np

8import scipy

9from beartype import beartype as typechecker

10from beartype.typing import List, Optional, Tuple, Union, cast

11from jaxtyping import Complex, Float, Int, Integer, Real

12from scipy import optimize

13from scipy.spatial.transform import Rotation

15import kwave.utils.typing as kt

16from kwave.utils.math import compute_linear_transform, compute_rotation_between_vectors

18from ..data import Vector

19from .conversion import db2neper, neper2db

20from .data import scale_SI

21from .math import Rx, Ry, Rz, compute_linear_transform, cosd, sind

22from .matlab import ind2sub, matlab_assign, matlab_find, sub2ind

23from .matrix import max_nd

24from .tictoc import TicToc

26# GLOBALS

27# define literals (ref: http://www.wolframalpha.com/input/?i=golden+angle)

28GOLDEN_ANGLE = 2.39996322972865332223155550663361385312499901105811504

29PACKING_NUMBER = 7 # 2*pi

32def make_cart_disc(

33 disc_pos: np.ndarray, radius: float, focus_pos: np.ndarray, num_points: int, plot_disc: bool = False, use_spiral: bool = False

34) -> np.ndarray:

35 """

36 Create evenly distributed Cartesian points covering a disc.

38 Args:

39 disc_pos: Cartesian position of the center of the disc given as a

40 two (2D) or three (3D) element tuple [m].

41 radius: Radius of the disc [m].

42 focus_pos: Any point on the beam axis of the disc given as a three

43 element tuple [fx, fy, fz] [m]. Can be set to None to define a disc in the x-y plane.

44 num_points: Number of points on the disc.

45 plot_disc: Boolean controlling whether the Cartesian points are plotted (default = False).

46 use_spiral: Boolean controlling whether the Cartesian points are chosen using a spiral sampling pattern.

47 Concentric sampling is used by default.

49 Returns:

50 disc: 2 x num_points or 3 x num_points array of Cartesian coordinates.

51 """

53 # check input values

54 if radius <= 0:

55 raise ValueError("The radius must be positive.")

57 def make_spiral_circle_points(num_points: int, radius: float) -> np.ndarray:

58 # compute spiral parameters

59 theta = lambda t: GOLDEN_ANGLE * t

60 C = np.pi * radius**2 / (num_points - 1)

61 r = lambda t: np.sqrt(C * t / np.pi)

63 # compute canonical spiral points

64 t = np.linspace(0, num_points - 1, num_points)

65 p0 = np.multiply(np.vstack((np.cos(theta(t)), np.sin(theta(t)))), r(t))

66 return p0

68 def make_concentric_circle_points(num_points: int, radius: float) -> Tuple[np.ndarray, int]:

69 assert num_points >= 1, "The number of points must be greater or equal to 1."

71 num_radial = int(np.ceil(np.sqrt(num_points / np.pi)))

73 try:

74 d_radial = radius / (num_radial - 1)

75 except ZeroDivisionError:

76 d_radial = float("inf")

78 r = np.arange(num_radial) * (radius - d_radial / 2) / (num_radial - 1)

80 # recompute the number of points that will be created below

81 num_points = 1

82 for k in range(2, num_radial + 1):

83 num_theta = round((k - 1) * PACKING_NUMBER)

84 num_points += num_theta

86 # compute canonical concentric circle points

87 points = np.full((2, num_points), np.nan)

88 points[:, 0] = [0, 0]

89 i_left = 1

90 for k in range(2, num_radial + 1):

91 num_theta = round((k - 1) * PACKING_NUMBER)

92 thetas = np.arange(num_theta) * 2 * np.pi / num_theta

93 p = r[k - 1] * np.vstack((np.cos(thetas), np.sin(thetas)))

94 i_right = i_left + num_theta

95 points[:, i_left:i_right] = p

96 i_left = i_right

97 return points, num_points

99 if use_spiral:

100 p0 = make_spiral_circle_points(num_points, radius)

101

102 else:

103 # otherwise use concentric circles (note that the num_points is increased

104 # to ensure a full set of concentric rings)

105

106 p0, num_points = make_concentric_circle_points(num_points, radius)

107

108 # add z-dimension points if in 3D

109 if len(disc_pos) == 3:

110 p0 = np.vstack((p0, np.zeros((1, num_points))))

111

112 # if 3D and focus_pos is defined, rotate the canonical points to give the

113 # specified disc

114 if len(disc_pos) == 3:

115 # check the focus position isn't coincident with the disc position

116 if all(disc_pos == focus_pos):

117 raise ValueError("The focus_pos must be different from the disc_pos.")

118

119 # compute rotation matrix and apply

120 R, _ = compute_linear_transform(disc_pos, focus_pos)

121 p0 = np.dot(R, p0)

122

123 # shift the disc to the appropriate center

124 disc = p0 + np.array(disc_pos).reshape(-1, 1)

125

126 # plot results

127 if plot_disc:

128 # select suitable axis scaling factor

129 _, scale, prefix, unit = scale_SI(np.max(disc))

130

131 # create the figure

132 if len(disc_pos) == 2:

133 plt.figure()

134 plt.plot(disc[1, :] * scale, disc[0, :] * scale, ".")

135 plt.gca().invert_yaxis()

136 plt.xlabel(f"y-position [{prefix}m]")

137 plt.ylabel(f"x-position [{prefix}m]")

138 plt.axis("equal")

139 else:

140 fig = plt.figure()

141 ax = fig.add_subplot(111, projection="3d")

142 ax.plot3D(disc[0, :] * scale, disc[1, :] * scale, disc[2, :] * scale, ".")

143 ax.set_xlabel(f"[{prefix}m]")

144 ax.set_ylabel(f"[{prefix}m]")

145 ax.set_zlabel(f"[{prefix}m]")

146 ax.axis("equal")

147 ax.grid(True)

148 ax.box(True)

149

150 return np.squeeze(disc)

151

152

153@typechecker

154def make_cart_bowl(

155 bowl_pos: np.ndarray, radius: float, diameter: float, focus_pos: np.ndarray, num_points: int, plot_bowl: Optional[bool] = False

156) -> Float[np.ndarray, "3 NumPoints"]:

157 """

158 Create evenly distributed Cartesian points covering a bowl.

159

160 Args:

161 bowl_pos: Cartesian position of the centre of the rear surface of

162 the bowl given as a three element vector [bx, by, bz] [m].

163 radius: Radius of curvature of the bowl [m].

164 diameter: Diameter of the opening of the bowl [m].

165 focus_pos: Any point on the beam axis of the bowl given as a three

166 element vector [fx, fy, fz] [m].

167 num_points: Number of points on the bowl.

168 plot_bowl: Boolean controlling whether the Cartesian points are

169 plotted.

170

171 Returns:

172 3 x num_points array of Cartesian coordinates.

173

174 Examples:

175 bowl = makeCartBowl([0, 0, 0], 1, 2, [0, 0, 1], 100)

176 bowl = makeCartBowl([0, 0, 0], 1, 2, [0, 0, 1], 100, True)

177 """

178

179 # check input values

180 if radius <= 0:

181 raise ValueError("The radius must be positive.")

182 if diameter <= 0:

183 raise ValueError("The diameter must be positive.")

184 if diameter > 2 * radius:

185 raise ValueError("The diameter of the bowl must be equal or less than twice the radius of curvature.")

186 if np.all(bowl_pos == focus_pos):

187 raise ValueError("The focus_pos must be different from the bowl_pos.")

188

189 # check for infinite radius of curvature, and call makeCartDisc instead

190 if np.isinf(radius):

191 # bowl = make_cart_disc(bowl_pos, diameter / 2, focus_pos, num_points, plot_bowl)

192 # return bowl

193 raise NotImplementedError("make_cart_disc")

194

195 # compute arc angle from chord (ref: https://en.wikipedia.org/wiki/Chord_(geometry))

196 varphi_max = np.arcsin(diameter / (2 * radius))

197

198 # compute spiral parameters

199 theta = lambda t: GOLDEN_ANGLE * t

200 C = 2 * np.pi * (1 - np.cos(varphi_max)) / (num_points - 1)

201 varphi = lambda t: np.arccos(1 - C * t / (2 * np.pi))

202

203 # compute canonical spiral points

204 t = np.linspace(0, num_points - 1, num_points)

205 p0 = np.array([np.cos(theta(t)) * np.sin(varphi(t)), np.sin(theta(t)) * np.sin(varphi(t)), np.cos(varphi(t))])

206 p0 = radius * p0

207

208 # linearly transform the canonical spiral points to give bowl in correct orientation

209 R, b = compute_linear_transform(bowl_pos, focus_pos, radius)

210

211 if b.ndim == 1:

212 b = np.expand_dims(b, axis=-1) # expand dims for broadcasting

213

214 bowl = R @ p0 + b

215

216 # plot results

217 if plot_bowl is True:

218 # select suitable axis scaling factor

219 _, scale, prefix, unit = scale_SI(np.max(bowl))

220

221 # create the figure

222 fig = plt.figure()

223 ax = fig.add_subplot(111, projection="3d")

224 ax.scatter(bowl[0, :] * scale, bowl[1, :] * scale, bowl[2, :] * scale)

225 ax.set_xlabel("[" + prefix + unit + "]")

226 ax.set_ylabel("[" + prefix + unit + "]")

227 ax.set_zlabel("[" + prefix + unit + "]")

228 ax.set_box_aspect([1, 1, 1])

229 plt.grid(True)

230 plt.show()

231

232 return bowl

233

234

235def get_spaced_points(start: float, stop: float, n: int = 100, spacing: str = "linear") -> np.ndarray:

236 """

237 Generate a row vector of either logarithmically or linearly spaced points between `start` and `stop`.

238

239 When `spacing` is set to 'linear', the function is identical to the inbuilt `np.linspace` function.

240 When `spacing` is set to 'log', the function is similar to the inbuilt `np.logspace` function, except

241 that `start` and `stop` define the start and end numbers, not decades. For logarithmically spaced

242 points, `start` must be > 0. If `n` < 2, `stop` is returned.

243

244 Args:

245 start: start value for the spaced points

246 stop: end value for the spaced points

247 n: number of points to generate

248 spacing: type of spacing to use, either 'linear' or 'log'

249

250 Returns:

251 points: row vector of spaced points

252

253 Raises:

254 ValueError: if `stop` <= `start` or `spacing` is not 'linear' or 'log'

255

256 """

257

258 if stop <= start:

259 raise ValueError("`stop` must be larger than `start`.")

260

261 if spacing == "linear":

262 return np.linspace(start, stop, num=n)

263 elif spacing == "log":

264 return np.geomspace(start, stop, num=n)

265 else:

266 raise ValueError(f"`spacing` {spacing} is not a valid argument. Choose from 'linear' or 'log'.")

267

268

269def fit_power_law_params(a0: float, y: float, c0: float, f_min: float, f_max: float, plot_fit: bool = False) -> Tuple[float, float]:

270 """

271 Calculate absorption parameters that fit a power law over a given frequency range.

272

273 This function calculates the absorption parameters that should be defined in the simulation functions

274 to achieve the desired power law absorption behavior defined by `a0` and `y`. This takes into account

275 the actual absorption behavior exhibited by the fractional Laplacian wave equation.

276

277 This fitting is required when using large absorption values or high frequencies, as the fractional

278 Laplacian wave equation solved in `kspaceFirstOrderND` and `kspaceSecondOrder` no longer encapsulates

279 absorption of the form `a = a0*f^y`.

280

281 The returned values should be used to define `medium.alpha_coeff` and `medium.alpha_power` within the

282 simulation functions. The absorption behavior over the frequency range `f_min`:`f_max` will then

283 follow the power law defined by `a0` and `y`.

284

285 Args:

286 a0: coefficient in the power law absorption equation

287 y: exponent in the power law absorption equation

288 c0: speed of sound in the medium

289 f_min: minimum frequency in the range to fit the power law

290 f_max: maximum frequency in the range to fit the power law

291 plot_fit: whether to plot the fit

292

293 Returns:

294 A tuple of the absorption coefficient and fitted exponent of the power law absorption equation.

295

296 """

297

298 # define frequency axis

299 f = get_spaced_points(f_min, f_max, 200)

300 w = 2 * np.pi * f

301 # convert user defined a0 to Nepers/((rad/s)^y m)

302 a0_np = db2neper(a0, y)

303

304 desired_absorption = a0_np * w**y

305

306 def abs_func(trial_vals):

307 """Second-order absorption error"""

308 a0_np_trial, y_trial = trial_vals

309

310 actual_absorption = (

311 a0_np_trial * w**y_trial / (1 - (y_trial + 1) * a0_np_trial * c0 * np.tan(np.pi * y_trial / 2) * w ** (y_trial - 1))

312 )

313

314 absorption_error = np.sqrt(np.sum((desired_absorption - actual_absorption) ** 2))

315

316 return absorption_error

317

318 a0_np_fit, y_fit = optimize.fmin(abs_func, [a0_np, y])

319

320 a0_fit = neper2db(a0_np_fit, y_fit)

321

322 if plot_fit:

323 raise NotImplementedError

324

325 return a0_fit, y_fit

326

327

328def power_law_kramers_kronig(w: np.ndarray, w0: float, c0: float, a0: float, y: float) -> np.ndarray:

329 """

330 Compute the variation in sound speed for an attenuating medium using the Kramers-Kronig for power law attenuation.

331

332 This function computes the variation in sound speed for an attenuating medium using the Kramers-Kronig

333 formula for power law attenuation, where `att = a0 * w^y`. The power law parameters must be in Nepers/m,

334 with the frequency in rad/s. The variation is given about the sound speed `c0` at a reference frequency `w0`.

335

336 Args:

337 w: input frequency array [rad/s]

338 w0: reference frequency [rad/s]

339 c0: sound speed at w0 [m/s]

340 a0: power law coefficient [Nepers/((rad/s)^y m)]

341 y: power law exponent, where 0 < y < 3

342

343 Returns:

344 Variation of sound speed with w [m/s]

345

346 """

347

348 if 0 >= y or y >= 3:

349 logging.log(logging.WARN, f"{UserWarning.__name__}: y must be within the interval (0,3)")

350 warnings.warn("y must be within the interval (0,3)", UserWarning)

351 c_kk = c0 * np.ones_like(w)

352 elif y == 1:

353 # Kramers-Kronig for y = 1

354 c_kk = 1 / (1 / c0 - 2 * a0 * np.log(w / w0) / np.pi)

355 else:

356 # Kramers-Kronig for 0 < y < 1 and 1 < y < 3

357 c_kk = 1 / (1 / c0 + a0 * np.tan(y * np.pi / 2) * (w ** (y - 1) - w0 ** (y - 1)))

358

359 return c_kk

360

361

362def water_absorption(f: float, temp: Union[float, kt.NP_DOMAIN]) -> Union[float, kt.NP_DOMAIN]:

363 """

364 Calculates the ultrasonic absorption in distilled

365 water at a given temperature and frequency using a 7 th order

366 polynomial fitted to the data given by Pinkerton (1949).

367

368

369 Args:

370 f: f frequency value [MHz]

371 T: water temperature value [degC]

372

373 Returns:

374 abs: absorption [dB / cm]

375

376 Examples:

377 >>> abs = waterAbsorption(f, T)

378

379 References:

380 [1] J. M. M. Pinkerton (1949) "The Absorption of Ultrasonic Waves in Liquids and its

381 Relation to Molecular Constitution,"

382 Proceedings of the Physical Society. Section B, 2, 129-141

383

384 """

385

386 NEPER2DB = 8.686

387 # check temperature is within range

388 if not np.all([np.all(temp >= 0.0), np.all(temp <= 60.0)]):

389 raise Warning("Temperature outside range of experimental data")

390

391 # conversion factor between Nepers and dB NEPER2DB = 8.686;

392 # coefficients for 7th order polynomial fit

393 a = [

394 56.723531840522710,

395 -2.899633796917384,

396 0.099253401567561,

397 -0.002067402501557,

398 2.189417428917596e-005,

399 -6.210860973978427e-008,

400 -6.402634551821596e-010,

401 3.869387679459408e-012,

402 ]

403

404 # compute absorption

405 a_on_fsqr = (

406 a[0] + a[1] * temp + a[2] * temp**2 + a[3] * temp**3 + a[4] * temp**4 + a[5] * temp**5 + a[6] * temp**6 + a[7] * temp**7

407 ) * 1e-17

408

409 abs = NEPER2DB * 1e12 * f**2 * a_on_fsqr

410 return abs

411

412

413def water_sound_speed(temp: Union[float, kt.NP_DOMAIN]) -> Union[float, kt.NP_DOMAIN]:

414 """

415 Calculate the sound speed in distilled water with temperature according to Marczak (1997)

416

417 Args:

418 temp: The temperature of the water in degrees Celsius.

419

420 Returns:

421 c: The sound speed in distilled water in m/s.

422

423 Raises:

424 ValueError: if `temp` is not between 0 and 95

425

426 References:

427 [1] R. Marczak (1997). "The sound velocity in water as a function of temperature".

428 Journal of Research of the National Institute of Standards and Technology, 102(6), 561-567.

429

430 """

431

432 # check limits

433 if not np.all([np.all(temp >= 0.0), np.all(temp <= 95.0)]):

434 raise ValueError("`temp` must be between 0 and 95.")

435

436 # find value

437 p = [2.787860e-9, -1.398845e-6, 3.287156e-4, -5.779136e-2, 5.038813, 1.402385e3]

438 c = np.polyval(p, temp)

439 return c

440

441

442def water_density(temp: Union[kt.NUMERIC, np.ndarray]) -> Union[kt.NUMERIC, np.ndarray]:

443 """

444 Calculate the density of air-saturated water with temperature.

445

446 This function calculates the density of air-saturated water at a given temperature using the 4th order polynomial

447 given by Jones [1].

448

449 Args:

450 temp: water temperature in the range 5 to 40 [degC]

451

452 Returns:

453 density: density of water [kg/m^3]

454

455 Raises:

456 ValueError: if `temp` is not between 5 and 40

457

458 References:

459 [1] F. E. Jones and G. L. Harris (1992) "ITS-90 Density of Water Formulation for Volumetric Standards Calibration,"

460 Journal of Research of the National Institute of Standards and Technology, 97(3), 335-340.

461

462 """

463

464 # check limits

465 if not np.all([np.all(np.asarray(temp) >= 5.0), np.all(np.asarray(temp) <= 40.0)]):

466 raise ValueError("`temp` must be between 5 and 40.")

467

468 # calculate density of air-saturated water

469 density = 999.84847 + 6.337563e-2 * temp - 8.523829e-3 * temp**2 + 6.943248e-5 * temp**3 - 3.821216e-7 * temp**4

470 return density

471

472

473def water_non_linearity(temp: Union[float, kt.NP_DOMAIN]) -> Union[float, kt.NP_DOMAIN]:

474 """

475 Calculates the parameter of nonlinearity B/A at a

476 given temperature using a fourth-order polynomial fitted to the data

477 given by Beyer (1960).

478

479 Args:

480 temp: water temperature in the range 0 to 100 [degC]

481

482 Returns:

483 BonA: parameter of nonlinearity

484

485 Examples:

486 >>> BonA = waterNonlinearity(T)

487

488 References:

489 [1] R. T. Beyer (1960) "Parameter of nonlinearity in fluids,"

490 J. Acoust. Soc. Am., 32(6), 719-721.

491

492 """

493

494 # check limits

495 if not np.all([np.all(temp >= 0.0), np.all(temp <= 100.0)]):

496 raise ValueError("`temp` must be between 0 and 100.")

497

498 # find value

499 p = [-4.587913769504693e-08, 1.047843302423604e-05, -9.355518377254833e-04, 5.380874771364909e-2, 4.186533937275504]

500 BonA = np.polyval(p, temp)

501 return BonA

502

503

504@typechecker

505def make_ball(

506 grid_size: Vector, ball_center: Vector, radius: int, plot_ball: bool = False, binary: bool = False

507) -> Union[kt.NP_ARRAY_INT_3D, kt.NP_ARRAY_BOOL_3D]:

508 """

509 Creates a binary map of a filled ball within a 3D grid.

510

511 Args:

512 grid_size: size of the 3D grid in [grid points].

513 ball_center: centre of the ball in [grid points]

514 radius: ball radius [grid points].

515 plot_ball: whether to plot the ball using voxelPlot (default = False).

516 binary: whether to return the ball map as a double precision matrix (False) or a logical matrix (True) (default = False).

517

518 Returns:

519 ball: 3D binary map of a filled ball.

520

521 """

522

523 # define literals

524 MAGNITUDE = 1

525 assert grid_size.shape == (3,), "grid_size must be a 3 element vector"

526 assert ball_center.shape == (3,), "ball_center must be a 3 element vector"

527

528 # force integer values

529 grid_size = cast(Vector, grid_size.astype(int))

530 ball_center = cast(Vector, ball_center.astype(int))

531

532 # check for zero values

533 for i in range(3):

534 if ball_center[i] == 0:

535 ball_center[i] = int(floor(grid_size[i] / 2)) + 1

536

537 # create empty matrix

538 ball = np.zeros(grid_size).astype(bool if binary else int)

539

540 # define np.pixel map

541 r = make_pixel_map(grid_size, shift=[0, 0, 0])

542

543 # create ball

544 ball[r <= radius] = MAGNITUDE

545

546 # shift centre

547 ball_center = ball_center - np.ceil(grid_size / 2).astype(int)

548 ball = np.roll(ball, ball_center, axis=(0, 1, 2))

549

550 # plot results

551 if plot_ball:

552 raise NotImplementedError

553 # voxelPlot(double(ball))

554 return ball

555

556

557@typechecker

558def make_cart_sphere(

559 radius: Union[float, int], num_points: int, center_pos: Vector = Vector([0, 0, 0]), plot_sphere: bool = False

560) -> Float[np.ndarray, "3 NumPoints"]:

561 """

562 Cart_sphere creates a set of points in Cartesian coordinates defining a sphere.

563

564 Args:

565 radius: the radius of the sphere.

566 num_points: the number of points to be generated.

567 center_pos: the coordinates of the center of the sphere. Defaults to (0, 0, 0).

568 plot_sphere: whether to plot the sphere. Defaults to False.

569

570 Returns:

571 The points on the sphere.

572

573 """

574 # generate angle functions using the Golden Section Spiral method

575 inc = np.pi * (3 - np.sqrt(5))

576 off = 2 / num_points

577 k = np.arange(0, num_points)

578 y = k * off - 1 + (off / 2)

579 r = np.sqrt(1 - (y**2))

580 phi = k * inc

581

582 if num_points <= 0:

583 raise ValueError("num_points must be greater than 0")

584

585 # create the sphere

586 sphere = radius * np.concatenate([np.cos(phi) * r[np.newaxis, :], y[np.newaxis, :], np.sin(phi) * r[np.newaxis, :]])

587

588 # offset if needed

589 sphere = sphere + center_pos[:, None]

590

591 # plot results

592 if plot_sphere:

593 # select suitable axis scaling factor

594 [x_sc, scale, prefix, _] = scale_SI(np.max(sphere))

595

596 # create the figure

597 plt.figure()

598 plt.style.use("seaborn-poster")

599 ax = plt.axes(projection="3d")

600 ax.plot3D(sphere[0, :] * scale, sphere[1, :] * scale, sphere[2, :] * scale, ".")

601 ax.set_xlabel(f"[{prefix} m]")

602 ax.set_ylabel(f"[{prefix} m]")

603 ax.set_zlabel(f"[{prefix} m]")

604 ax.axis("auto")

605 ax.grid()

606 plt.show()

607

608 return sphere.squeeze()

609

610

611def make_cart_circle(

612 radius: float, num_points: int, center_pos: Vector = Vector([0, 0]), arc_angle: float = 2 * np.pi, plot_circle: bool = False

613) -> Float[np.ndarray, "2 NumPoints"]:

614 """

615 Create a set of points in cartesian coordinates defining a circle or arc.

616

617 This function creates a set of points in cartesian coordinates defining a circle or arc.

618

619 Args:

620 radius: radius of the circle or arc

621 num_points: number of points to generate

622 center_pos: center position of the circle or arc

623 arc_angle: arc angle in radians

624 plot_circle: whether to plot the circle or arc

625

626 Returns:

627 2 x `num_points` array of cartesian coordinates

628

629 """

630

631 # check for arc_angle input

632 if arc_angle == 2 * np.pi:

633 full_circle = True

634 else:

635 full_circle = False

636

637 n_steps = num_points if full_circle else num_points - 1

638

639 # create angles

640 angles = np.arange(0, num_points) * arc_angle / n_steps + np.pi / 2

641

642 # create cartesian grid

643 circle = np.concatenate([radius * np.cos(angles[np.newaxis, :]), radius * np.sin(-angles[np.newaxis])])

644

645 # offset if needed

646 circle = circle + center_pos[:, None]

647

648 # plot results

649 if plot_circle:

650 # select suitable axis scaling factor

651 [_, scale, prefix, _] = scale_SI(np.max(abs(circle)))

652

653 # create the figure

654 plt.figure()

655 plt.plot(circle[1, :] * scale, circle[0, :] * scale, "b.")

656 plt.xlabel([f"y-position [{prefix} m]"])

657 plt.ylabel([f"x-position [{prefix} m]"])

658 plt.axis("equal")

659 plt.show()

660

661 return np.squeeze(circle)

662

663

664@typechecker

665def make_disc(grid_size: Vector, center: Vector, radius, plot_disc=False) -> kt.NP_ARRAY_BOOL_2D:

666 """

667 Create a binary map of a filled disc within a 2D grid.

668

669 This function creates a binary map of a filled disc within a two-dimensional grid. The disc position is denoted by 1's

670 in the matrix with 0's elsewhere. A single grid point is taken as the disc centre, so the total diameter of the disc

671 will always be an odd number of grid points. If used within a k-Wave grid where dx != dy, the disc will appear oval

672 shaped. If part of the disc overlaps the grid edge, the rest of the disc will wrap to the grid edge on the opposite

673 side.

674

675 Args:

676 grid_size: A 2D vector of the grid size in grid points.

677 center: A 2D vector of the disc centre in grid points.

678 radius: The radius of the disc.

679 plot_disc: If set to True, the disc will be plotted using Matplotlib.

680

681 Returns:

682 A binary map of the disc in the 2D grid.

683

684 """

685

686 assert len(grid_size) == 2, "Grid size must be 2D."

687 assert len(center) == 2, "Center must be 2D."

688

689 # define literals

690 MAGNITUDE = 1

691

692 # force integer values

693 grid_size = grid_size.round().astype(int)

694 center = center.round().astype(int)

695

696 # check for zero values

697 center.x = center.x if center.x != 0 else np.floor(grid_size.x / 2).astype(int) + 1

698 center.y = center.y if center.y != 0 else np.floor(grid_size.y / 2).astype(int) + 1

699

700 # check the inputs

701 assert np.all(0 < center) and np.all(center <= grid_size), "Disc center must be within grid."

702

703 # create empty matrix

704 disc = np.zeros(grid_size, dtype=bool)

705

706 # define np.pixel map

707 r = make_pixel_map(grid_size, shift=[0, 0])

708

709 # create disc

710 disc[r <= radius] = MAGNITUDE

711

712 # shift centre

713 center = center - np.ceil(grid_size / 2).astype(int)

714 disc = np.roll(disc, center, axis=(0, 1))

715

716 # create the figure

717 if plot_disc:

718 raise NotImplementedError

719 return disc

720

721

722@typechecker

723def make_circle(

724 grid_size: Vector, center: Vector, radius: Real[kt.ScalarLike, ""], arc_angle: Optional[float] = None, plot_circle: bool = False

725) -> kt.NP_ARRAY_INT_2D:

726 """

727 Create a binary map of a circle within a 2D grid.

728

729 This function creates a binary map of a circle (or arc) using the midpoint circle algorithm within a two-dimensional grid.

730 The circle position is denoted by 1's in the matrix with 0's elsewhere. A single grid point is taken as the circle

731 centre, so the total diameter will always be an odd number of grid points. The centre of the circle and the radius

732 are not constrained by the grid dimensions, so it is possible to create sections of circles or a blank image if none

733 of the circle intersects the grid.

734

735 Args:

736 grid_size: A 2D vector of the grid size in grid points.

737 center: A 2D vector of the circle centre in grid points.

738 radius: The radius of the circle.

739 arc_angle: The angle of the circular arc in degrees. If set to None, a full circle will be created.

740 plot_circle: If set to True, the circle will be plotted using Matplotlib.

741

742 Returns:

743 A binary map of the circle in the 2D grid.

744

745 """

746

747 assert len(grid_size) == 2, "Grid size must be 2D"

748 assert len(center) == 2, "Center must be 2D"

749

750 # define literals

751 MAGNITUDE = 1

752

753 if arc_angle is None:

754 arc_angle = 2 * np.pi

755 elif arc_angle > 2 * np.pi:

756 arc_angle = 2 * np.pi

757 elif arc_angle < 0:

758 arc_angle = 0

759

760 # force integer values

761 grid_size = grid_size.round().astype(int)

762 center = center.round().astype(int)

763 radius = int(round(radius))

764

765 # check for zero values

766 center.x = center.x if center.x != 0 else int(floor(grid_size.x / 2)) + 1

767 center.y = center.y if center.y != 0 else int(floor(grid_size.y / 2)) + 1

768 cx, cy = center

769

770 # create empty matrix

771 circle = np.zeros(grid_size, dtype=int)

772

773 # initialise loop variables

774 x = 0

775 y = radius

776 d = 1 - radius

777

778 if (cx >= 1) and (cx <= grid_size.x) and ((cy - y) >= 1) and ((cy - y) <= grid_size.y):

779 circle[cx - 1, cy - y - 1] = MAGNITUDE

780

781 # draw the remaining cardinal points

782 px = [cx, cx + y, cx - y]

783 py = [cy + y, cy, cy]

784 for point_index, (px_i, py_i) in enumerate(zip(px, py)):

785 # check whether the point is within the arc made by arc_angle, and lies

786 # within the grid

787 if (np.arctan2(px_i - cx, py_i - cy) + np.pi) <= arc_angle:

788 if (px_i >= 1) and (px_i <= grid_size.x) and (py_i >= 1) and (py_i <= grid_size.y):

789 circle[px_i - 1, py_i - 1] = MAGNITUDE

790

791 # loop through the remaining points using the midpoint circle algorithm

792 while x < (y - 1):

793 x = x + 1

794 if d < 0:

795 d = d + x + x + 1

796 else:

797 y = y - 1

798 a = x - y + 1

799 d = d + a + a

800

801 # setup point indices (break coding standard for readability)

802 px = [x + cx, y + cx, y + cx, x + cx, -x + cx, -y + cx, -y + cx, -x + cx]

803 py = [y + cy, x + cy, -x + cy, -y + cy, -y + cy, -x + cy, x + cy, y + cy]

804

805 # loop through each point

806 for point_index, (px_i, py_i) in enumerate(zip(px, py)):

807 # check whether the point is within the arc made by arc_angle, and

808 # lies within the grid

809 if (np.arctan2(px_i - cx, py_i - cy) + np.pi) <= arc_angle:

810 if (px_i >= 1) and (px_i <= grid_size.x) and (py_i >= 1) and (py_i <= grid_size.y):

811 circle[px_i - 1, py_i - 1] = MAGNITUDE

812

813 if plot_circle:

814 plt.imshow(circle, cmap="gray_r")

815 plt.ylabel("x-position [grid points]")

816 plt.xlabel("y-position [grid points]")

817 plt.show()

818

819 return circle

820

821

822def make_pixel_map(grid_size: Vector, shift=None, origin_size="single") -> np.ndarray:

823 """

824 Generates a matrix with values of the distance of each pixel from the center of a grid.

825

826 This function generates a matrix populated with values of how far each pixel in a grid is from the center. The center

827 can be a single pixel or a double pixel, and the optional input parameter 'OriginSize' controls this. For grids where

828 the dimension size and center pixel size are not both odd or even, the optional input parameter 'Shift' can be used to

829 control the location of the center point.

830

831 Args:

832 grid_size: A 2D or 3D vector of the grid size in grid points.

833 *args: additional optional arguments

834

835 Returns:

836 r: pixel-radius

837

838 Examples:

839

840 Single pixel origin size for odd and even (with 'Shift' = [1 1] and

841 [0 0], respectively) grid sizes:

842

843 x x x x x x x x x x x

844 x 0 x x x x x x 0 x x

845 x x x x x 0 x x x x x

846 x x x x x x x x

847

848 Double pixel origin size for even and odd (with 'Shift' = [1 1] and

849 [0 0], respectively) grid sizes:

850

851 x x x x x x x x x x x x x x

852 x 0 0 x x x x x x x 0 0 x x

853 x 0 0 x x x 0 0 x x 0 0 x x

854 x x x x x x 0 0 x x x x x x

855 x x x x x x x x x x

856

857 By default, a single pixel centre is used which is shifted towards

858 the final row and column.

859

860 """

861 assert len(grid_size) == 2 or len(grid_size) == 3, "Grid size must be a 2 or 3 element vector."

862

863 # define defaults

864 shift_def = 1

865

866 Nx = grid_size[0]

867 Ny = grid_size[1]

868 Nz = None

869 if len(grid_size) == 3:

870 Nz = grid_size[2]

871

872 # detect whether the inputs are for two or three dimensions

873 map_dimension = 2 if Nz is None else 3

874 if shift is None:

875 shift = [shift_def] * map_dimension

876

877 # catch input errors

878 assert origin_size in ["single", "double"], "Unknown setting for optional input Center."

879

880 assert (

881 len(shift) == map_dimension

882 ), f"Optional input Shift must have {map_dimension} elements for {map_dimension} dimensional input parameters."

883

884 if map_dimension == 2:

885 # create the maps for each dimension

886 nx = create_pixel_dim(Nx, origin_size, shift[0])

887 ny = create_pixel_dim(Ny, origin_size, shift[1])

888

889 # create plaid grids

890 r_x, r_y = np.meshgrid(nx, ny, indexing="ij")

891

892 # extract the pixel radius

893 r = np.sqrt(r_x**2 + r_y**2)

894 if map_dimension == 3:

895 # create the maps for each dimension

896 nx = create_pixel_dim(Nx, origin_size, shift[0])

897 ny = create_pixel_dim(Ny, origin_size, shift[1])

898 nz = create_pixel_dim(Nz, origin_size, shift[2])

899

900 # create plaid grids

901 r_x, r_y, r_z = np.meshgrid(nx, ny, nz, indexing="ij")

902

903 # extract the pixel radius

904 r = np.sqrt(r_x**2 + r_y**2 + r_z**2)

905 return r

906

907

908def create_pixel_dim(Nx: int, origin_size: float, shift: float) -> Tuple[np.ndarray, float]:

909 """

910 Create an array of pixel dimensions and a pixel size.

911

912 Args:

913 Nx: The number of pixels in the x-dimension.

914 origin_size: The size of the origin in the x-dimension.

915 shift: The shift of the pixels in the x-dimension.

916

917 Returns:

918 The pixel dimensions.

919

920 """

921

922 # Nested function to create the pixel radius variable

923

924 # grid dimension has an even number of points

925 if Nx % 2 == 0:

926 # pixel numbering has a single centre point

927 if origin_size == "single":

928 # centre point is shifted towards the final pixel

929 if shift == 1:

930 nx = np.arange(-Nx / 2, Nx / 2 - 1 + 1, 1)

931

932 # centre point is shifted towards the first pixel

933 else:

934 nx = np.arange(-Nx / 2 + 1, Nx / 2 + 1, 1)

935

936 # pixel numbering has a double centre point

937 else:

938 nx = np.hstack([np.arange(-Nx / 2 + 1, 0 + 1, 1), np.arange(0, Nx / 2 - 1 + 1, 1)])

939

940 # grid dimension has an odd number of points

941 else:

942 # pixel numbering has a single centre point

943 if origin_size == "single":

944 nx = np.arange(-(Nx - 1) / 2, (Nx - 1) / 2 + 1, 1)

945

946 # pixel numbering has a double centre point

947 else:

948 # centre point is shifted towards the final pixel

949 if shift == 1:

950 nx = np.hstack([np.arange(-(Nx - 1) / 2, 0 + 1, 1), np.arange(0, (Nx - 1) / 2 - 1 + 1, 1)])

951

952 # centre point is shifted towards the first pixel

953 else:

954 nx = np.hstack([np.arange(-(Nx - 1) / 2 + 1, 0 + 1, 1), np.arange(0, (Nx - 1) / 2 + 1, 1)])

955 return nx

956

957

958@typechecker

959def make_line(

960 grid_size: Vector,

961 startpoint: Union[Tuple[Int[kt.ScalarLike, ""], Int[kt.ScalarLike, ""]], Int[np.ndarray, "2"]],

962 endpoint: Optional[Union[Tuple[Int[kt.ScalarLike, ""], Int[kt.ScalarLike, ""]], Int[np.ndarray, "2"]]] = None,

963 angle: Optional[Float[kt.ScalarLike, ""]] = None,

964 length: Optional[Int[kt.ScalarLike, ""]] = None,

965) -> kt.NP_ARRAY_BOOL_2D:

966 """

967 Generate a line shape with a given start and end point, angle, or length.

968

969 Args:

970 grid_size: The size of the grid in pixels.

971 startpoint: The start point of the line, given as a tuple of x and y coordinates.

972 endpoint: The end point of the line, given as a tuple of x and y coordinates.

973 If not specified, the line is drawn from the start point at a given angle and length.

974 angle: The angle of the line in radians, measured counterclockwise from the x-axis.

975 If not specified, the line is drawn from the start point to the end point.

976 length: The length of the line in pixels.

977 If not specified, the line is drawn from the start point to the end point.

978

979 Returns:

980 line: A 2D array of the same size as the input parameters,

981 with a value of 1 for pixels that are part of the line and 0 for pixels that are not.

982 """

983 assert len(grid_size) == 2, "Grid size must be a 2-element vector."

984

985 startpoint = np.array(startpoint, dtype=int)

986 if endpoint is not None:

987 endpoint = np.array(endpoint, dtype=int)

988

989 if len(startpoint) != 2:

990 raise ValueError("startpoint should be a two-element vector.")

991

992 if np.any(startpoint < 1) or startpoint[0] > grid_size.x or startpoint[1] > grid_size.y:

993 ValueError("The starting point must lie within the grid, between [1 1] and [grid_size.x grid_size.y].")

994

995 # =========================================================================

996 # LINE BETWEEN TWO POINTS OR ANGLED LINE?

997 # =========================================================================

998

999 if endpoint is not None:

1000 linetype = "AtoB"

1001 a, b = startpoint, endpoint

1002

1003 # Addition => Fix Matlab2Python indexing

1004 a -= 1

1005 b -= 1

1006 else:

1007 linetype = "angled"

1008 angle, linelength = angle, length

1009

1010 # =========================================================================

1011 # MORE INPUT CHECKING

1012 # =========================================================================

1013

1014 if linetype == "AtoB":

1015 # a and b must be different points

1016 if np.all(a == b):

1017 raise ValueError("The first and last points cannot be the same.")

1018

1019 # end point must be a two-element row vector

1020 if len(b) != 2:

1021 raise ValueError("endpoint should be a two-element vector.")

1022

1023 # a and b must be within the grid

1024 xx = np.array([a[0], b[0]], dtype=int)

1025 yy = np.array([a[1], b[1]], dtype=int)

1026 if np.any(a < 0) or np.any(b < 0) or np.any(xx > grid_size.x - 1) or np.any(yy > grid_size.y - 1):

1027 raise ValueError("Both the start and end points must lie within the grid.")

1028

1029 if linetype == "angled":

1030 # angle must lie between -np.pi and np.pi

1031 angle = angle % (2 * np.pi)

1032 if angle > np.pi:

1033 angle = angle - (2 * np.pi)

1034 elif angle < -np.pi:

1035 angle = angle + (2 * np.pi)

1036

1037 # =========================================================================

1038 # CALCULATE A LINE FROM A TO B

1039 # =========================================================================

1040

1041 if linetype == "AtoB":

1042 # define an empty grid to hold the line

1043 line = np.zeros(grid_size, dtype=bool)

1044

1045 # find the equation of the line

1046 m = (b[1] - a[1]) / (b[0] - a[0]) # gradient of the line

1047 c = a[1] - m * a[0] # where the line crosses the y axis

1048

1049 if abs(m) < 1:

1050 # start at the end with the smallest value of x

1051 if a[0] < b[0]:

1052 x, y = a

1053 x_end = b[0]

1054 else:

1055 x, y = b

1056 x_end = a[0]

1057

1058 # fill in the first point

1059 line[x, y] = 1

1060

1061 while x < x_end:

1062 # next points to try are

1063 poss_x = [x, x, x + 1, x + 1, x + 1]

1064 poss_y = [y - 1, y + 1, y - 1, y, y + 1]

1065

1066 # find the point closest to the line

1067 true_y = m * poss_x + c

1068 diff = (poss_y - true_y) ** 2

1069 index = matlab_find(diff == min(diff))[0]

1070

1071 # the next point

1072 x = poss_x[index[0] - 1]

1073 y = poss_y[index[0] - 1]

1074

1075 # add the point to the line

1076 line[x - 1, y - 1] = 1

1077

1078 elif not np.isinf(abs(m)):

1079 # start at the end with the smallest value of y

1080 if a[1] < b[1]:

1081 x = a[0]

1082 y = a[1]

1083 y_end = b[1]

1084 else:

1085 x = b[0]

1086 y = b[1]

1087 y_end = a[1]

1088

1089 # fill in the first point

1090 line[x, y] = 1

1091

1092 while y < y_end:

1093 # next points to try are

1094 poss_y = [y, y, y + 1, y + 1, y + 1]

1095 poss_x = [x - 1, x + 1, x - 1, x, x + 1]

1096

1097 # find the point closest to the line

1098 true_x = (poss_y - c) / m

1099 diff = (poss_x - true_x) ** 2

1100 index = matlab_find(diff == min(diff))[0]

1101

1102 # the next point

1103 x = poss_x[index[0] - 1]

1104 y = poss_y[index[0] - 1]

1105

1106 # add the point to the line

1107 line[x, y] = 1

1108

1109 else: # m = +-Inf

1110 # start at the end with the smallest value of y

1111 if a[1] < b[1]:

1112 x = a[0]

1113 y = a[1]

1114 y_end = b[1]

1115 else:

1116 x = b[0]

1117 y = b[1]

1118 y_end = a[1]

1119

1120 # fill in the first point

1121 line[x, y] = 1

1122

1123 while y < y_end:

1124 # next point

1125 y = y + 1

1126

1127 # add the point to the line

1128 line[x, y] = 1

1129

1130 # =========================================================================

1131 # CALCULATE AN ANGLED LINE

1132 # =========================================================================

1133

1134 elif linetype == "angled":

1135 # define an empty grid to hold the line

1136 line = np.zeros(grid_size, dtype=bool)

1137

1138 # start at the atart

1139 x, y = startpoint

1140

1141 # fill in the first point

1142 line[x - 1, y - 1] = 1

1143

1144 # initialise the current length of the line

1145 line_length = 0

1146

1147 if abs(angle) == np.pi:

1148 while line_length < linelength:

1149 # next point

1150 y = y + 1

1151

1152 # stop the points incrementing at the edges

1153 if y > grid_size.y:

1154 break

1155

1156 # add the point to the line

1157 line[x - 1, y - 1] = 1

1158

1159 # calculate the current length of the line

1160 line_length = np.sqrt((x - startpoint[0]) ** 2 + (y - startpoint[1]) ** 2)

1161

1162 elif (angle < np.pi) and (angle > np.pi / 2):

1163 # define the equation of the line

1164 m = -np.tan(angle - np.pi / 2) # gradient of the line

1165 c = y - m * x # where the line crosses the y axis

1166

1167 while line_length < linelength:

1168 # next points to try are

1169 poss_x = np.array([x - 1, x - 1, x])

1170 poss_y = np.array([y, y + 1, y + 1])

1171

1172 # find the point closest to the line

1173 true_y = m * poss_x + c

1174 diff = (poss_y - true_y) ** 2

1175 index = matlab_find(diff == min(diff))[0]

1176

1177 # the next point

1178 x = poss_x[index[0] - 1]

1179 y = poss_y[index[0] - 1]

1180

1181 # stop the points incrementing at the edges

1182 if (x < 0) or (y > grid_size.y - 1):

1183 break

1184

1185 # add the point to the line

1186 line[x - 1, y - 1] = 1

1187

1188 # calculate the current length of the line

1189 line_length = np.sqrt((x - startpoint[0]) ** 2 + (y - startpoint[1]) ** 2)

1190

1191 elif angle == np.pi / 2:

1192 while line_length < linelength:

1193 # next point

1194 x = x - 1

1195

1196 # stop the points incrementing at the edges

1197 if x < 1:

1198 break

1199

1200 # add the point to the line

1201 line[x - 1, y - 1] = 1

1202

1203 # calculate the current length of the line

1204 line_length = np.sqrt((x - startpoint[0]) ** 2 + (y - startpoint[1]) ** 2)

1205

1206 elif (angle < np.pi / 2) and (angle > 0):

1207 # define the equation of the line

1208 m = np.tan(np.pi / 2 - angle) # gradient of the line

1209 c = y - m * x # where the line crosses the y axis

1210

1211 while line_length < linelength:

1212 # next points to try are

1213 poss_x = np.array([x - 1, x - 1, x])

1214 poss_y = np.array([y, y - 1, y - 1])

1215

1216 # find the point closest to the line

1217 true_y = m * poss_x + c

1218 diff = (poss_y - true_y) ** 2

1219 index = matlab_find(diff == min(diff))[0]

1220

1221 # the next point

1222 x = poss_x[index[0] - 1]

1223 y = poss_y[index[0] - 1]

1224

1225 # stop the points incrementing at the edges

1226 if (x < 1) or (y < 1):

1227 break

1228

1229 # add the point to the line

1230 line[x - 1, y - 1] = 1

1231

1232 # calculate the current length of the line

1233 line_length = np.sqrt((x - startpoint[0]) ** 2 + (y - startpoint[1]) ** 2)

1234

1235 elif angle == 0:

1236 while line_length < linelength:

1237 # next point

1238 y = y - 1

1239

1240 # stop the points incrementing at the edges

1241 if y < 1:

1242 break

1243

1244 # add the point to the line

1245 line[x - 1, y - 1] = 1

1246

1247 # calculate the current length of the line

1248 line_length = np.sqrt((x - startpoint[0]) ** 2 + (y - startpoint[1]) ** 2)

1249

1250 elif (angle < 0) and (angle > -np.pi / 2):

1251 # define the equation of the line

1252 m = -np.tan(np.pi / 2 + angle) # gradient of the line

1253 c = y - m * x # where the line crosses the y axis

1254

1255 while line_length < linelength:

1256 # next points to try are

1257 poss_x = np.array([x + 1, x + 1, x])

1258 poss_y = np.array([y, y - 1, y - 1])

1259

1260 # find the point closest to the line

1261 true_y = m * poss_x + c

1262 diff = (poss_y - true_y) ** 2

1263 index = matlab_find(diff == min(diff))[0]

1264

1265 # the next point

1266 x = poss_x[index[0] - 1]

1267 y = poss_y[index[0] - 1]

1268

1269 # stop the points incrementing at the edges

1270 if (x > grid_size.x) or (y < 1):

1271 break

1272

1273 # add the point to the line

1274 line[x - 1, y - 1] = 1

1275

1276 # calculate the current length of the line

1277 line_length = np.sqrt((x - startpoint[0]) ** 2 + (y - startpoint[1]) ** 2)

1278

1279 elif angle == -np.pi / 2:

1280 while line_length < linelength:

1281 # next point

1282 x = x + 1

1283

1284 # stop the points incrementing at the edges

1285 if x > grid_size.x:

1286 break

1287

1288 # add the point to the line

1289 line[x - 1, y - 1] = 1

1290

1291 # calculate the current length of the line

1292 line_length = np.sqrt((x - startpoint[0]) ** 2 + (y - startpoint[1]) ** 2)

1293

1294 elif (angle < -np.pi / 2) and (angle > -np.pi):

1295 # define the equation of the line

1296 m = np.tan(-angle - np.pi / 2) # gradient of the line

1297 c = y - m * x # where the line crosses the y axis

1298

1299 while line_length < linelength:

1300 # next points to try are

1301 poss_x = np.array([x + 1, x + 1, x])

1302 poss_y = np.array([y, y + 1, y + 1])

1303

1304 # find the point closest to the line

1305 true_y = m * poss_x + c

1306 diff = (poss_y - true_y) ** 2

1307 index = matlab_find(diff == min(diff))[0]

1308

1309 # the next point

1310 x = poss_x[index[0] - 1]

1311 y = poss_y[index[0] - 1]

1312

1313 # stop the points incrementing at the edges

1314 if (x > grid_size.x) or (y > grid_size.y):

1315 break

1316

1317 # add the point to the line

1318 line[x - 1, y - 1] = 1

1319

1320 # calculate the current length of the line

1321 line_length = np.sqrt((x - startpoint[0]) ** 2 + (y - startpoint[1]) ** 2)

1322

1323 return line

1324

1325

1326@typechecker

1327def make_arc(

1328 grid_size: Vector, arc_pos: np.ndarray, radius: Real[kt.ScalarLike, ""], diameter: Int[kt.ScalarLike, ""], focus_pos: Vector

1329) -> Union[kt.NP_ARRAY_INT_2D, kt.NP_ARRAY_BOOL_2D]:

1330 """

1331 Generates an arc shape with a given radius, diameter, and focus position.

1332

1333 Args:

1334 grid_size: The size of the grid, given as a 1D array with the number of pixels in each dimension.

1335 arc_pos: The position of the arc, given as a 1D array with the coordinates in each dimension.

1336 radius: The radius of the arc.

1337 diameter: The diameter of the arc.

1338 focus_pos: The position of the focus, given as a 1D array with the coordinates in each dimension.

1339

1340 Returns:

1341 np.ndarray: A 2D array with the arc shape.

1342 """

1343 assert len(grid_size) == 2, "The grid size must be a 2D vector."

1344 assert len(arc_pos) == 2, "The arc position must be a 2D vector."

1345 assert len(focus_pos) == 2, "The focus position must be a 2D vector."

1346

1347 # force integer input values

1348 grid_size = grid_size.round().astype(int)

1349 arc_pos = arc_pos.round().astype(int)

1350 diameter = int(round(diameter))

1351 focus_pos = focus_pos.round().astype(int)

1352

1353 try:

1354 radius = int(radius)

1355 except OverflowError:

1356 radius = float(radius)

1357

1358 # check the input ranges

1359 if np.any(grid_size < 1):

1360 raise ValueError("The grid size must be positive.")

1361 if radius <= 0:

1362 raise ValueError("The radius must be positive.")

1363

1364 if diameter <= 0:

1365 raise ValueError("The diameter must be positive.")

1366

1367 if np.any(arc_pos < 1) or np.any(arc_pos > grid_size):

1368 raise ValueError("The centre of the arc must be within the grid.")

1369

1370 if diameter > 2 * radius:

1371 raise ValueError("The diameter of the arc must be less than twice the radius of curvature.")

1372

1373 if diameter % 2 != 1:

1374 raise ValueError("The diameter must be an odd number of grid points.")

1375

1376 if np.all(arc_pos == focus_pos):

1377 raise ValueError("The focus_pos must be different to the arc_pos.")

1378

1379 # assign variable names to vector components

1380 Nx, Ny = grid_size

1381 ax, ay = arc_pos

1382 fx, fy = focus_pos

1383

1384 # =========================================================================

1385 # CREATE ARC

1386 # =========================================================================

1387

1388 if not np.isinf(radius):

1389 # find half the arc angle

1390 half_arc_angle = np.arcsin(diameter / 2 / radius)

1391

1392 # find centre of circle on which the arc lies

1393 distance_cf = np.sqrt((ax - fx) ** 2 + (ay - fy) ** 2)

1394 cx = round(radius / distance_cf * (fx - ax) + ax)

1395 cy = round(radius / distance_cf * (fy - ay) + ay)

1396 c = np.array([cx, cy])

1397

1398 # create circle

1399 arc = make_circle(grid_size, Vector([cx, cy]), radius)

1400

1401 # form vector from the geometric arc centre to the arc midpoint

1402 v1 = arc_pos - c

1403

1404 # calculate length of vector

1405 l1 = np.sqrt(sum((arc_pos - c) ** 2))

1406

1407 # extract all points that form part of the arc

1408 arc_ind = matlab_find(arc, mode="eq", val=1)

1409

1410 # loop through the arc points

1411 for arc_ind_i in arc_ind:

1412 # extract the indices of the current point

1413 x_ind, y_ind = ind2sub([Nx, Ny], arc_ind_i)

1414 p = np.array([x_ind, y_ind])

1415

1416 # form vector from the geometric arc centre to the current point

1417 v2 = p - c

1418

1419 # calculate length of vector

1420 l2 = np.sqrt(sum((p - c) ** 2))

1421

1422 # find the angle between the two vectors using the dot product,

1423 # normalised using the vector lengths

1424 theta = np.arccos(sum(v1 * v2 / (l1 * l2)))

1425

1426 # if the angle is greater than the half angle of the arc, remove

1427 # it from the arc

1428 if theta > half_arc_angle:

1429 arc = matlab_assign(arc, arc_ind_i - 1, 0)

1430 else:

1431 # calculate arc direction angle, then rotate by 90 degrees

1432 ang = np.arctan((fx - ax) / (fy - ay)) + np.pi / 2

1433

1434 # draw lines to create arc with infinite radius

1435 arc = np.logical_or(

1436 make_line(grid_size, arc_pos, endpoint=None, angle=ang, length=(diameter - 1) // 2),

1437 make_line(grid_size, arc_pos, endpoint=None, angle=(ang + np.pi), length=(diameter - 1) // 2),

1438 )

1439 return arc

1440

1441

1442def make_pixel_map_point(grid_size: Vector, centre_pos: np.ndarray) -> np.ndarray:

1443 """

1444 Generates a map of the distance of each pixel from a given centre position.

1445

1446 Args:

1447 grid_size: The size of the grid, given as a 1D array with the number of pixels in each dimension.

1448 centre_pos: The position of the centre, given as a 1D array with the coordinates in each dimension.

1449

1450 Returns:

1451 np.ndarray: A 2D array with the distance of each pixel from the given centre position.

1452

1453 Raises:

1454 ValueError: If `grid_size` and `centre_pos` do not have the same number of dimensions.

1455 """

1456 # check for number of dimensions

1457 num_dim = len(grid_size)

1458

1459 # check that centre_pos has the same dimensions

1460 if len(grid_size) != len(centre_pos):

1461 raise ValueError("The inputs centre_pos and grid_size must have the same number of dimensions.")

1462

1463 if num_dim == 2:

1464 # assign inputs and force to be integers

1465 Nx, Ny = grid_size.astype(int)

1466 cx, cy = centre_pos.astype(int)

1467

1468 # generate index vectors in each dimension

1469 nx = np.arange(0, Nx) - cx + 1

1470 ny = np.arange(0, Ny) - cy + 1

1471

1472 # combine index matrices

1473 pixel_map = np.zeros((Nx, Ny))

1474 pixel_map += (nx**2)[:, None]

1475 pixel_map += (ny**2)[None, :]

1476 pixel_map = np.sqrt(pixel_map)

1477

1478 elif num_dim == 3:

1479 # assign inputs and force to be integers

1480 Nx, Ny, Nz = grid_size.astype(int)

1481 cx, cy, cz = centre_pos.astype(int)

1482

1483 # generate index vectors in each dimension

1484 nx = np.arange(0, Nx) - cx + 1

1485 ny = np.arange(0, Ny) - cy + 1

1486 nz = np.arange(0, Nz) - cz + 1

1487

1488 # combine index matrices

1489 pixel_map = np.zeros((Nx, Ny, Nz))

1490 pixel_map += (nx**2)[:, None, None]

1491 pixel_map += (ny**2)[None, :, None]

1492 pixel_map += (nz**2)[None, None, :]

1493 pixel_map = np.sqrt(pixel_map)

1494

1495 else:

1496 # throw error

1497 raise ValueError("Grid size must be 2 or 3D.")

1498

1499 return pixel_map

1500

1501

1502def make_pixel_map_plane(grid_size: Vector, normal: np.ndarray, point: np.ndarray) -> np.ndarray:

1503 """

1504 Generates a pixel map of a plane with given normal vector and point.

1505

1506 Args:

1507 grid_size: The size of the grid as a NumPy array [Nx, Ny, Nz].

1508 normal: The normal vector of the plane as a NumPy array [nx, ny, nz].

1509 point: A point on the plane as a NumPy array [px, py, pz].

1510

1511 Returns:

1512 pixel_map: A 2D array with the distance of each pixel from the given plane.

1513

1514 Raises:

1515 ValueError: If the normal vector is zero.

1516 """ # error checking

1517 if np.all(normal == 0):

1518 raise ValueError("Normal vector should not be zero.")

1519

1520 # check for number of dimensions

1521 num_dim = len(grid_size)

1522

1523 if num_dim == 2:

1524 # assign inputs and force to be integers

1525 Nx = round(grid_size[0])

1526 Ny = round(grid_size[1])

1527

1528 # create coordinate meshes

1529 [px, py] = np.meshgrid(Nx, Ny)

1530 [pointx, pointy] = np.meshgrid(np.ones((1, Nx)) * point[0], np.ones(1, Ny) * point[1])

1531 [nx, ny] = np.meshgrid(np.ones((1, Nx)) * normal[0], np.ones(1, Ny) * normal[2])

1532

1533 # calculate distance according to Eq. (6) at

1534 # http://mathworld.wolfram.com/Point-PlaneDistance.html

1535 pixel_map = np.abs((px - pointx) * nx + (py - pointy) * ny) / np.sqrt(sum(normal**2))

1536

1537 elif num_dim == 3:

1538 # assign inputs and force to be integers

1539 Nx = np.round(grid_size[0])

1540 Ny = np.round(grid_size[1])

1541 Nz = np.round(grid_size[2])

1542

1543 # create coordinate meshes

1544 px, py, pz = np.meshgrid(np.arange(1, Nx + 1), np.arange(1, Ny + 1), np.arange(1, Nz + 1), indexing="ij")

1545 pointx, pointy, pointz = np.meshgrid(np.ones(Nx) * point[0], np.ones(Ny) * point[1], np.ones(Nz) * point[2], indexing="ij")

1546 nx, ny, nz = np.meshgrid(np.ones(Nx) * normal[0], np.ones(Ny) * normal[1], np.ones(Nz) * normal[2], indexing="ij")

1547

1548 # calculate distance according to Eq. (6) at

1549 # http://mathworld.wolfram.com/Point-PlaneDistance.html

1550 pixel_map = np.abs((px - pointx) * nx + (py - pointy) * ny + (pz - pointz) * nz) / np.sqrt(sum(normal**2))

1551

1552 else:

1553 # throw error

1554 raise ValueError("Grid size must be 2 or 3D.")

1555

1556 return pixel_map

1557

1558

1559@typechecker

1560def make_bowl(

1561 grid_size: Vector,

1562 bowl_pos: Vector,

1563 radius: Union[int, float],

1564 diameter: Real[kt.ScalarLike, ""],

1565 focus_pos: Vector,

1566 binary: bool = False,

1567 remove_overlap: bool = False,

1568) -> Union[kt.NP_ARRAY_BOOL_3D, kt.NP_ARRAY_INT_3D]:

1569 """

1570 Generate a matrix representing a bowl-shaped object in 3D space.

1571

1572 This function generates a 3D matrix representing a bowl-shaped object with the specified radius and diameter. The position

1573 of the bowl and the focus point can be specified, as well as whether to return a binary matrix or a matrix with

1574 continuous values. The optional parameter 'remove_overlap' can be used to remove any overlap with the surrounding grid.

1575

1576 Args:

1577 grid_size: size of the grid in each dimension

1578 bowl_pos: position of the bowl in the grid

1579 radius: radius of the bowl

1580 diameter: diameter of the bowl

1581 focus_pos: position of the focus point in the grid

1582 binary: whether to return a binary matrix

1583 remove_overlap: whether to remove overlap with the surrounding grid

1584

1585 Returns:

1586 matrix: 3D matrix representing the bowl-shaped object

1587

1588 Raises:

1589 ValueError: if any of the input arguments are outside the valid range

1590 """

1591 assert len(grid_size) == 3, "Grid size must be 3D."

1592 assert len(bowl_pos) == 3, "Bowl position must be 3D."

1593 assert len(focus_pos) == 3, "Focus position must be 3D."

1594

1595 # =========================================================================

1596 # DEFINE LITERALS

1597 # =========================================================================

1598

1599 # threshold used to find the closest point to the radius

1600 THRESHOLD = 0.5

1601

1602 # number of grid points to expand the bounding box compared to

1603 # sqrt(2)*diameter

1604 BOUNDING_BOX_EXP = 2

1605

1606 # =========================================================================

1607 # INPUT CHECKING

1608 # =========================================================================

1609

1610 # force integer input values

1611 grid_size = np.round(grid_size).astype(int)

1612 bowl_pos = np.round(bowl_pos).astype(int)

1613 focus_pos = np.round(focus_pos).astype(int)

1614 diameter = np.round(diameter)

1615 radius = np.round(radius)

1616

1617 # check the input ranges

1618 if np.any(grid_size < 1):

1619 raise ValueError("The grid size must be positive.")

1620 if np.any(bowl_pos < 1) or np.any(bowl_pos > grid_size):

1621 raise ValueError("The centre of the bowl must be within the grid.")

1622 if radius <= 0:

1623 raise ValueError("The radius must be positive.")

1624 if diameter <= 0:

1625 raise ValueError("The diameter must be positive.")

1626 if diameter > (2 * radius):

1627 raise ValueError("The diameter of the bowl must be less than twice the radius of curvature.")

1628 if diameter % 2 == 0:

1629 raise ValueError("The diameter must be an odd number of grid points.")

1630 if np.all(bowl_pos == focus_pos):

1631 raise ValueError("The focus_pos must be different to the bowl_pos.")

1632

1633 # =========================================================================

1634 # BOUND THE GRID TO SPEED UP CALCULATION

1635 # =========================================================================

1636

1637 # create bounding box slightly larger than bowl diameter * sqrt(2)

1638 Nx = np.round(np.sqrt(2) * diameter).astype(int) + BOUNDING_BOX_EXP

1639 Ny = Nx

1640 Nz = Nx

1641 grid_size_sm = Vector([Nx, Ny, Nz])

1642

1643 # set the bowl position to be the centre of the bounding box

1644 bx = np.ceil(Nx / 2).astype(int)

1645 by = np.ceil(Ny / 2).astype(int)

1646 bz = np.ceil(Nz / 2).astype(int)

1647 bowl_pos_sm = np.array([bx, by, bz])

1648

1649 # set the focus position to be in the direction specified by the user

1650 fx = bx + (focus_pos[0] - bowl_pos[0])

1651 fy = by + (focus_pos[1] - bowl_pos[1])

1652 fz = bz + (focus_pos[2] - bowl_pos[2])

1653 focus_pos_sm = [fx, fy, fz]

1654

1655 # preallocate storage variable

1656 if binary:

1657 bowl_sm = np.zeros((Nx, Ny, Nz), dtype=bool)

1658 else:

1659 bowl_sm = np.zeros((Nx, Ny, Nz))

1660

1661 # =========================================================================

1662 # CREATE DISTANCE MATRIX

1663 # =========================================================================

1664

1665 if not np.isinf(radius):

1666 # find half the arc angle

1667 half_arc_angle = np.arcsin(diameter / (2 * radius))

1668

1669 # find centre of sphere on which the bowl lies

1670 distance_cf = np.sqrt((bx - fx) ** 2 + (by - fy) ** 2 + (bz - fz) ** 2)

1671 cx = round(radius / distance_cf * (fx - bx) + bx)

1672 cy = round(radius / distance_cf * (fy - by) + by)

1673 cz = round(radius / distance_cf * (fz - bz) + bz)

1674 c = np.array([cx, cy, cz])

1675

1676 # generate matrix with distance from the centre

1677 pixel_map = make_pixel_map_point(grid_size_sm, c)

1678

1679 # set search radius to bowl radius

1680 search_radius = radius

1681

1682 else:

1683 # generate matrix with distance from the centre

1684 pixel_map = make_pixel_map_plane(grid_size_sm, bowl_pos_sm - focus_pos_sm, bowl_pos_sm)

1685

1686 # set search radius to 0 (the disc is flat)

1687 search_radius = 0

1688

1689 # calculate distance from search radius

1690 pixel_map = np.abs(pixel_map - search_radius)

1691

1692 # =========================================================================

1693 # DIMENSION 1

1694 # =========================================================================

1695

1696 # find the grid point that corresponds to the outside of the bowl in the

1697 # first dimension in both directions (the index gives the distance along

1698 # this dimension)

1699 value_forw, index_forw = pixel_map.min(axis=0), pixel_map.argmin(axis=0)

1700 value_back, index_back = np.flip(pixel_map, axis=0).min(axis=0), np.flip(pixel_map, axis=0).argmin(axis=0)

1701

1702 # extract the linear index in the y-z plane of the values that lie on the

1703 # bowl surface

1704 yz_ind_forw = matlab_find(value_forw < THRESHOLD)

1705 yz_ind_back = matlab_find(value_back < THRESHOLD)

1706

1707 # use these subscripts to extract the x-index of the grid points that lie

1708 # on the bowl surface

1709 x_ind_forw = index_forw.flatten(order="F")[yz_ind_forw - 1] + 1

1710 x_ind_back = index_back.flatten(order="F")[yz_ind_back - 1] + 1

1711

1712 # convert the linear index to equivalent subscript values

1713 y_ind_forw, z_ind_forw = ind2sub([Ny, Nz], yz_ind_forw)

1714 y_ind_back, z_ind_back = ind2sub([Ny, Nz], yz_ind_back)

1715

1716 # combine x-y-z indices into a linear index

1717 linear_index_forw = sub2ind([Nx, Ny, Nz], x_ind_forw - 1, y_ind_forw - 1, z_ind_forw - 1) + 1

1718 linear_index_back = sub2ind([Nx, Ny, Nz], Nx - x_ind_back, y_ind_back - 1, z_ind_back - 1) + 1

1719

1720 # assign these values to the bowl

1721 bowl_sm = matlab_assign(bowl_sm, linear_index_forw - 1, 1)

1722 bowl_sm = matlab_assign(bowl_sm, linear_index_back - 1, 1)

1723

1724 # set existing bowl values to a distance of zero in the pixel map (this

1725 # avoids problems with overlapping pixels)

1726 pixel_map[bowl_sm == 1] = 0

1727

1728 # =========================================================================

1729 # DIMENSION 2

1730 # =========================================================================

1731

1732 # find the grid point that corresponds to the outside of the bowl in the

1733 # second dimension in both directions (the pixel map is first re-ordered to

1734 # [X, Y, Z] -> [Y, Z, X])

1735 pixel_map_temp = np.transpose(pixel_map, (1, 2, 0))

1736 value_forw, index_forw = pixel_map_temp.min(axis=0), pixel_map_temp.argmin(axis=0)

1737 value_back, index_back = np.flip(pixel_map_temp, axis=0).min(axis=0), np.flip(pixel_map_temp, axis=0).argmin(axis=0)

1738 del pixel_map_temp

1739

1740 # extract the linear index in the y-z plane of the values that lie on the

1741 # bowl surface

1742 zx_ind_forw = matlab_find(value_forw < THRESHOLD)

1743 zx_ind_back = matlab_find(value_back < THRESHOLD)

1744

1745 # use these subscripts to extract the y-index of the grid points that lie

1746 # on the bowl surface

1747 y_ind_forw = index_forw.flatten(order="F")[zx_ind_forw - 1] + 1

1748 y_ind_back = index_back.flatten(order="F")[zx_ind_back - 1] + 1

1749

1750 # convert the linear index to equivalent subscript values

1751 z_ind_forw, x_ind_forw = ind2sub([Nz, Nx], zx_ind_forw)

1752 z_ind_back, x_ind_back = ind2sub([Nz, Nx], zx_ind_back)

1753

1754 # combine x-y-z indices into a linear index

1755 linear_index_forw = sub2ind([Nx, Ny, Nz], x_ind_forw - 1, y_ind_forw - 1, z_ind_forw - 1) + 1

1756 linear_index_back = sub2ind([Nx, Ny, Nz], x_ind_back - 1, Ny - y_ind_back, z_ind_back - 1) + 1

1757

1758 # assign these values to the bowl

1759 bowl_sm = matlab_assign(bowl_sm, linear_index_forw - 1, 1)

1760 bowl_sm = matlab_assign(bowl_sm, linear_index_back - 1, 1)

1761

1762 # set existing bowl values to a distance of zero in the pixel map (this

1763 # avoids problems with overlapping pixels)

1764 pixel_map[bowl_sm == 1] = 0

1765

1766 # =========================================================================

1767 # DIMENSION 3

1768 # =========================================================================

1769

1770 # find the grid point that corresponds to the outside of the bowl in the

1771 # third dimension in both directions (the pixel map is first re-ordered to

1772 # [X, Y, Z] -> [Z, X, Y])

1773 pixel_map_temp = np.transpose(pixel_map, (2, 0, 1))

1774 value_forw, index_forw = pixel_map_temp.min(axis=0), pixel_map_temp.argmin(axis=0)

1775 value_back, index_back = np.flip(pixel_map_temp, axis=0).min(axis=0), np.flip(pixel_map_temp, axis=0).argmin(axis=0)

1776 del pixel_map_temp

1777

1778 # extract the linear index in the y-z plane of the values that lie on the

1779 # bowl surface

1780 xy_ind_forw = matlab_find(value_forw < THRESHOLD)

1781 xy_ind_back = matlab_find(value_back < THRESHOLD)

1782

1783 # use these subscripts to extract the z-index of the grid points that lie

1784 # on the bowl surface

1785 z_ind_forw = index_forw.flatten(order="F")[xy_ind_forw - 1] + 1

1786 z_ind_back = index_back.flatten(order="F")[xy_ind_back - 1] + 1

1787

1788 # convert the linear index to equivalent subscript values

1789 x_ind_forw, y_ind_forw = ind2sub([Nx, Ny], xy_ind_forw)

1790 x_ind_back, y_ind_back = ind2sub([Nx, Ny], xy_ind_back)

1791

1792 # combine x-y-z indices into a linear index

1793 linear_index_forw = sub2ind([Nx, Ny, Nz], x_ind_forw - 1, y_ind_forw - 1, z_ind_forw - 1) + 1

1794 linear_index_back = sub2ind([Nx, Ny, Nz], x_ind_back - 1, y_ind_back - 1, Nz - z_ind_back) + 1

1795

1796 # assign these values to the bowl

1797 bowl_sm = matlab_assign(bowl_sm, linear_index_forw - 1, 1)

1798 bowl_sm = matlab_assign(bowl_sm, linear_index_back - 1, 1)

1799

1800 # =========================================================================

1801 # RESTRICT SPHERE TO BOWL

1802 # =========================================================================

1803

1804 # remove grid points within the sphere that do not form part of the bowl

1805 if not np.isinf(radius):

1806 # form vector from the geometric bowl centre to the back of the bowl

1807 v1 = bowl_pos_sm - c

1808

1809 # calculate length of vector

1810 l1 = np.sqrt(sum((bowl_pos_sm - c) ** 2))

1811

1812 # loop through the non-zero elements in the bowl matrix

1813 bowl_ind = matlab_find(bowl_sm == 1)[:, 0]

1814 for bowl_ind_i in bowl_ind:

1815 # extract the indices of the current point

1816 x_ind, y_ind, z_ind = ind2sub([Nx, Ny, Nz], bowl_ind_i)

1817 p = np.array([x_ind, y_ind, z_ind])

1818

1819 # form vector from the geometric bowl centre to the current point

1820 # on the bowl

1821 v2 = p - c

1822

1823 # calculate length of vector

1824 l2 = np.sqrt(sum((p - c) ** 2))

1825

1826 # find the angle between the two vectors using the dot product,

1827 # normalised using the vector lengths

1828 theta = np.arccos(sum(v1 * v2 / (l1 * l2)))

1829

1830 # # alternative calculation normalised using radius of curvature

1831 # theta2 = acos(sum( v1 .* v2 ./ radius**2 ))

1832

1833 # if the angle is greater than the half angle of the bowl, remove

1834 # it from the bowl

1835 if theta > half_arc_angle:

1836 bowl_sm = matlab_assign(bowl_sm, bowl_ind_i - 1, 0)

1837

1838 else:

1839 # form a distance map from the centre of the disc

1840 pixelMapPoint = make_pixel_map_point(grid_size_sm, bowl_pos_sm)

1841

1842 # set all points in the disc greater than the diameter to zero

1843 bowl_sm[pixelMapPoint > (diameter / 2)] = 0

1844

1845 # =========================================================================

1846 # REMOVE OVERLAPPED POINTS

1847 # =========================================================================

1848

1849 if remove_overlap:

1850 # define the shapes that capture the overlapped points, along with the

1851 # corresponding mask of which point to delete

1852 overlap_shapes = []

1853 overlap_delete = []

1854

1855 shape = np.zeros((3, 3, 3))

1856 shape[0, 0, :] = 1

1857 shape[1, 1, :] = 1

1858 shape[2, 2, :] = 1

1859 shape[0, 1, 1] = 1

1860 shape[1, 2, 1] = 1

1861 overlap_shapes.append(shape)

1862

1863 delete = np.zeros((3, 3, 3))

1864 delete[0, 1, 1] = 1

1865 delete[0, 2, 1] = 1

1866 overlap_delete.append(delete)

1867

1868 shape = np.zeros((3, 3, 3))

1869 shape[0, 0, :] = 1

1870 shape[1, 1, :] = 1

1871 shape[2, 2, :] = 1

1872 shape[0, 1, 1] = 1

1873 overlap_shapes.append(shape)

1874

1875 delete = np.zeros((3, 3, 3))

1876 delete[0, 1, 1] = 1

1877 overlap_delete.append(delete)

1878

1879 shape = np.zeros((3, 3, 3))

1880 shape[0:2, 0, :] = 1

1881 shape[2, 1, :] = 1

1882 shape[1, 1, 1] = 1

1883 overlap_shapes.append(shape)

1884

1885 delete = np.zeros((3, 3, 3))

1886 delete[1, 1, 1] = 1

1887 overlap_delete.append(delete)

1888

1889 shape = np.zeros((3, 3, 3))

1890 shape[0, 0, :] = 1

1891 shape[1, 1, :] = 1

1892 shape[2, 2, :] = 1

1893 shape[0, 1, 0] = 1

1894 overlap_shapes.append(shape)

1895

1896 delete = np.zeros((3, 3, 3))

1897 delete[0, 1, 0] = 1

1898 overlap_delete.append(delete)

1899

1900 shape = np.zeros((3, 3, 3))

1901 shape[0:2, 1, :] = 1

1902 shape[2, 2, :] = 1

1903 shape[2, 1, 0] = 1

1904 overlap_shapes.append(shape)

1905

1906 delete = np.zeros((3, 3, 3))

1907 delete[2, 1, 0] = 1

1908 overlap_delete.append(delete)

1909

1910 shape = np.zeros((3, 3, 3))

1911 shape[0, :, 2] = 1

1912 shape[1, :, 1] = 1

1913 shape[1, :, 0] = 1

1914 shape[2, 1, 0] = 1

1915 overlap_shapes.append(shape)

1916

1917 delete = np.zeros((3, 3, 3))

1918 delete[2, 1, 0] = 1

1919 overlap_delete.append(delete)

1920

1921 shape = np.zeros((3, 3, 3))

1922 shape[0, 2, :] = 1

1923 shape[1, 0:2, :] = 1

1924 shape[2, 0, 0] = 1

1925 overlap_shapes.append(shape)

1926

1927 delete = np.zeros((3, 3, 3))

1928 delete[2, 0, 0] = 1

1929 overlap_delete.append(delete)

1930

1931 shape = np.zeros((3, 3, 3))

1932 shape[:, :, 1] = 1

1933 shape[0, 0, 0] = 1

1934 overlap_shapes.append(shape)

1935

1936 delete = np.zeros((3, 3, 3))

1937 delete[0, 0, 0] = 1

1938 overlap_delete.append(delete)

1939

1940 shape = np.zeros((3, 3, 3))

1941 shape[0, :, 0] = 1

1942 shape[1, :, 1] = 1

1943 shape[1, :, 2] = 1

1944 shape[1, 1, 0] = 1

1945 overlap_shapes.append(shape)

1946

1947 delete = np.zeros((3, 3, 3))

1948 delete[1, 1, 0] = 1

1949 overlap_delete.append(delete)

1950

1951 shape = np.zeros((3, 3, 3))

1952 shape[1:3, 2, 0] = 1

1953 shape[0, 2, 1:3] = 1

1954 shape[0, 1, 2] = 1

1955 shape[1, 1, 1] = 1

1956 shape[2, 1, 0] = 1

1957 shape[1:3, 0, 1] = 1

1958 shape[1, 0, 2] = 1

1959 overlap_shapes.append(shape)

1960

1961 delete = np.zeros((3, 3, 3))

1962 delete[1, 0, 1] = 1

1963 overlap_delete.append(delete)

1964

1965 # set loop flag

1966 points_remaining = True

1967

1968 # initialise deleted point counter

1969 deleted_points = 0

1970

1971 # set list of possible permutations

1972 perm_list = [[0, 1, 2], [0, 2, 1], [1, 2, 0], [1, 0, 2], [2, 0, 1], [2, 1, 0]]

1973

1974 while points_remaining:

1975 # get linear index of non-zero bowl elements

1976 index_mat = matlab_find(bowl_sm > 0)[:, 0]

1977

1978 # set Boolean delete variable

1979 delete_point = False

1980

1981 # loop through all points on the bowl, and find the all the points with

1982 # more than 8 neighbours

1983 index = 0

1984 for index, index_mat_i in enumerate(index_mat):

1985 # extract subscripts for current point

1986 cx, cy, cz = ind2sub([Nx, Ny, Nz], index_mat_i)

1987

1988 # ignore edge points

1989 if (cx > 1) and (cx < Nx) and (cy > 1) and (cy < Ny) and (cz > 1) and (cz < Nz):

1990 # extract local region around current point

1991 region = bowl_sm[cx - 1 : cx + 1, cy - 1 : cy + 1, cz - 1 : cz + 1] # FARID might not work

1992

1993 # if there's more than 8 neighbours, check the point for

1994 # deletion

1995 if (region.sum() - 1) > 8:

1996 # loop through the different shapes

1997 for shape_index in range(len(overlap_shapes)):

1998 # check every permutation of the shape, and apply the

1999 # deletion mask if the pattern matches

2000

2001 # loop through possible shape permutations

2002 for ind1 in range(len(perm_list)):

2003 # get shape and delete mask

2004 overlap_s = overlap_shapes[shape_index]

2005 overlap_d = overlap_delete[shape_index]

2006

2007 # permute

2008 overlap_s = np.transpose(overlap_s, perm_list[ind1])

2009 overlap_d = np.transpose(overlap_d, perm_list[ind1])

2010

2011 # loop through possible shape reflections

2012 for ind2 in range(1, 8):

2013 # flipfunc the shape

2014 if ind2 == 2:

2015 overlap_s = np.flip(overlap_s, axis=0)

2016 overlap_d = np.flip(overlap_d, axis=0)

2017 elif ind2 == 3:

2018 overlap_s = np.flip(overlap_s, axis=1)

2019 overlap_d = np.flip(overlap_d, axis=1)

2020 elif ind2 == 4:

2021 overlap_s = np.flip(overlap_s, axis=2)

2022 overlap_d = np.flip(overlap_d, axis=2)

2023 elif ind2 == 5:

2024 overlap_s = np.flip(np.flip(overlap_s, axis=0), axis=1)

2025 overlap_d = np.flip(np.flip(overlap_d, axis=0), axis=1)

2026 elif ind2 == 6:

2027 overlap_s = np.flip(np.flip(overlap_s, axis=0), axis=2)

2028 overlap_d = np.flip(np.flip(overlap_d, axis=0), axis=2)

2029 elif ind2 == 7:

2030 overlap_s = np.flip(np.flip(overlap_s, axis=1), axis=2)

2031 overlap_d = np.flip(np.flip(overlap_d, axis=1), axis=2)

2032

2033 # rotate the shape 4 x 90 degrees

2034 for ind3 in range(4):

2035 # check if the shape matches

2036 if np.all(overlap_s == region):

2037 delete_point = True

2038

2039 # break from loop if a match is found

2040 if delete_point:

2041 break

2042

2043 # rotate shape

2044 overlap_s = np.rot90(overlap_s)

2045 overlap_d = np.rot90(overlap_d)

2046

2047 # break from loop if a match is found

2048 if delete_point:

2049 break

2050

2051 # break from loop if a match is found

2052 if delete_point:

2053 break

2054

2055 # remove point from bowl if required, and update

2056 # counter

2057 if delete_point:

2058 bowl_sm[cx - 1 : cx + 1, cy - 1 : cy + 1, cz - 1 : cz + 1] = bowl_sm[

2059 cx - 1 : cx + 1, cy - 1 : cy + 1, cz - 1 : cz + 1

2060 ] * np.bitwise_not(overlap_d).astype(float) # Farid won't work probably

2061 deleted_points = deleted_points + 1

2062 break

2063

2064 # break from loop if a match is found

2065 if delete_point:

2066 break

2067

2068 # break from while loop if the outer for loop has completed

2069 # without deleting a point

2070 if index == (len(index_mat) - 1):

2071 points_remaining = False

2072

2073 # display status

2074 if deleted_points:

2075 logging.log(logging.INFO, "{deleted_points} overlapped points removed from bowl")

2076

2077 # =========================================================================

2078 # PLACE BOWL WITHIN LARGER GRID

2079 # =========================================================================

2080

2081 # preallocate storage variable

2082 if binary:

2083 bowl = np.zeros(grid_size, dtype=bool)

2084 else:

2085 bowl = np.zeros(grid_size, dtype=int)

2086

2087 # calculate position of bounding box within larger grid

2088 x1 = bowl_pos[0] - bx

2089 x2 = x1 + Nx

2090 y1 = bowl_pos[1] - by

2091 y2 = y1 + Ny

2092 z1 = bowl_pos[2] - bz

2093 z2 = z1 + Nz

2094

2095 # truncate bounding box if it falls outside the grid

2096 if x1 < 0:

2097 bowl_sm = bowl_sm[abs(x1) :, :, :]

2098 x1 = 0

2099 if y1 < 0:

2100 bowl_sm = bowl_sm[:, abs(y1) :, :]

2101 y1 = 0

2102 if z1 < 0:

2103 bowl_sm = bowl_sm[:, :, abs(z1) :]

2104 z1 = 0

2105 if x2 >= grid_size[0]:

2106 to_delete = x2 - grid_size[0]

2107 bowl_sm = bowl_sm[:-to_delete, :, :]

2108 x2 = grid_size[0]

2109 if y2 >= grid_size[1]:

2110 to_delete = y2 - grid_size[1]

2111 bowl_sm = bowl_sm[:, :-to_delete, :]

2112 y2 = grid_size[1]

2113 if z2 >= grid_size[2]:

2114 to_delete = z2 - grid_size[2]

2115 bowl_sm = bowl_sm[:, :, :-to_delete]

2116 z2 = grid_size[2]

2117

2118 # place bowl into grid

2119 bowl[x1:x2, y1:y2, z1:z2] = bowl_sm

2120

2121 return bowl

2122

2123

2124def make_multi_bowl(

2125 grid_size: Vector,

2126 bowl_pos: List[Tuple[int, int]],

2127 radius: int,

2128 diameter: int,

2129 focus_pos: Tuple[int, int],

2130 binary: bool = False,

2131 remove_overlap: bool = False,

2132) -> Tuple[np.ndarray, List[np.ndarray]]:

2133 """

2134 Generates a multi-bowl mask for an image given the size of the grid, the positions of the bowls,

2135 the radius of each bowl, the diameter of the bowls, and the position of the focus.

2136

2137 Args:

2138 grid_size: The size of the grid (assumed to be square).

2139 bowl_pos: A list of tuples containing the (x, y) coordinates of the center of each bowl.

2140 radius: The radius of each bowl.

2141 diameter: The diameter of the bowls.

2142 focus_pos: The (x, y) coordinates of the focus.

2143 binary: Whether to return a binary mask (default: False).

2144 remove_overlap: Whether to remove overlap between the bowls (default: False).

2145

2146 Returns:

2147 bowls:

2148 bowls_labeled:

2149 """

2150

2151 # check inputs

2152 if bowl_pos.shape[-1] != 3:

2153 raise ValueError("bowl_pos should contain 3 columns, with [bx, by, bz] in each row.")

2154

2155 if len(radius) != 1 and len(radius) != bowl_pos.shape[0]:

2156 raise ValueError("The number of rows in bowl_pos and radius does not match.")

2157

2158 if len(diameter) != 1 and len(diameter) != bowl_pos.shape[0]:

2159 raise ValueError("The number of rows in bowl_pos and diameter does not match.")

2160

2161 # force integer grid size values

2162 grid_size = np.round(grid_size).astype(int)

2163 bowl_pos = np.round(bowl_pos).astype(int)

2164 focus_pos = np.round(focus_pos).astype(int)

2165 diameter = np.round(diameter)

2166 radius = np.round(radius)

2167

2168 # =========================================================================

2169 # CREATE BOWLS

2170 # =========================================================================

2171

2172 # preallocate output matrices

2173 if binary:

2174 bowls = np.zeros(grid_size, dtype=bool)

2175 else:

2176 bowls = np.zeros(grid_size)

2177

2178 bowls_labelled = np.zeros(grid_size)

2179

2180 # loop for calling make_bowl

2181 for bowl_index in range(bowl_pos.shape[0]):

2182 # update command line status

2183 if bowl_index == 1:

2184 TicToc.tic()

2185 else:

2186 TicToc.toc(reset=True)

2187 logging.log(logging.INFO, f"Creating bowl {bowl_index} of {bowl_pos.shape[0]} ... ")

2188

2189 # get parameters for current bowl

2190 if bowl_pos.shape[0] > 1:

2191 bowl_pos_k = bowl_pos[bowl_index]

2192 else:

2193 bowl_pos_k = bowl_pos

2194 bowl_pos_k = Vector(bowl_pos_k)

2195

2196 if len(radius) > 1:

2197 radius_k = radius[bowl_index]

2198 else:

2199 radius_k = radius

2200

2201 if len(diameter) > 1:

2202 diameter_k = diameter[bowl_index]

2203 else:

2204 diameter_k = diameter

2205

2206 if focus_pos.shape[0] > 1:

2207 focus_pos_k = focus_pos[bowl_index]

2208 else:

2209 focus_pos_k = focus_pos

2210 focus_pos_k = Vector(focus_pos_k)

2211

2212 # create new bowl

2213 new_bowl = make_bowl(grid_size, bowl_pos_k, radius_k, diameter_k, focus_pos_k, remove_overlap=remove_overlap, binary=binary)

2214

2215 # add bowl to bowl matrix

2216 bowls = bowls + new_bowl

2217

2218 # add new bowl to labelling matrix

2219 bowls_labelled[new_bowl == 1] = bowl_index

2220

2221 TicToc.toc()

2222

2223 # check if any of the bowls are overlapping

2224 max_nd_val, _ = max_nd(bowls)

2225 if max_nd_val > 1:

2226 # display warning

2227 logging.log(logging.WARN, f"{max_nd_val - 1} bowls are overlapping")

2228

2229 # force the output to be binary

2230 bowls[bowls != 0] = 1

2231

2232 return bowls, bowls_labelled

2233

2234

2235@typechecker

2236def make_multi_arc(

2237 grid_size: Vector, arc_pos: np.ndarray, radius: Union[int, np.ndarray], diameter: Union[int, np.ndarray], focus_pos: np.ndarray

2238) -> Tuple[kt.NP_ARRAY_FLOAT_2D, kt.NP_ARRAY_FLOAT_2D]:

2239 """

2240 Generates a multi-arc mask for an image given the size of the grid,

2241 the positions and properties of the arcs, and the position of the focus.

2242

2243 Args:

2244 grid_size: The size of the grid (assumed to be square).

2245 arc_pos: An array containing the (x, y) coordinates of the center of each arc.

2246 radius: The radius of each arc. Can be a single value or an array with one value for each arc.

2247 diameter: The diameter of the arcs. Can be a single value or an array with one value for each arc.

2248 focus_pos: The (x, y) coordinates of the focus.

2249

2250 Returns:

2251 arcs: A binary mask of the arcs.

2252 arcs_labelled: A labelled mask of the arcs.

2253

2254 Raises:

2255 ValueError: If the shape of arc_pos is not (N, 2), if the number of rows in arc_pos and radius do not match,

2256 or if the number of rows in arc_pos and diameter do not match.

2257 """

2258 # check inputs

2259 if arc_pos.shape[-1] != 2:

2260 raise ValueError("arc_pos should contain 2 columns, with [ax, ay] in each row.")

2261

2262 if len(radius) != 1 and len(radius) != arc_pos.shape[0]:

2263 raise ValueError("The number of rows in arc_pos and radius does not match.")

2264

2265 if len(diameter) != 1 and len(diameter) != arc_pos.shape[0]:

2266 raise ValueError("The number of rows in arc_pos and diameter does not match.")

2267

2268 # force integer grid size values

2269 grid_size = grid_size.round().astype(int)

2270 arc_pos = arc_pos.round().astype(int)

2271 diameter = diameter.round()

2272 focus_pos = focus_pos.round().astype(int)

2273 radius = radius.round()

2274

2275 # =========================================================================

2276 # CREATE ARCS

2277 # =========================================================================

2278

2279 # create empty matrix

2280 arcs = np.zeros(grid_size)

2281 arcs_labelled = np.zeros(grid_size)

2282

2283 # loop for calling make_arc

2284 for k in range(arc_pos.shape[0]):

2285 # get parameters for current arc

2286 if arc_pos.shape[0] > 1:

2287 arc_pos_k = arc_pos[k]

2288 else:

2289 arc_pos_k = arc_pos

2290

2291 if len(radius) > 1:

2292 radius_k = radius[k]

2293 else:

2294 radius_k = radius

2295

2296 if len(diameter) > 1:

2297 diameter_k = diameter[k]

2298 else:

2299 diameter_k = diameter

2300

2301 if focus_pos.shape[0] > 1:

2302 focus_pos_k = focus_pos[k]

2303 else:

2304 focus_pos_k = focus_pos

2305

2306 # create new arc

2307 new_arc = make_arc(grid_size, arc_pos_k, radius_k, diameter_k, Vector(focus_pos_k))

2308

2309 # add arc to arc matrix

2310 arcs = arcs + new_arc

2311

2312 # add new arc to labelling matrix

2313 arcs_labelled[new_arc == 1] = k

2314

2315 # check if any of the arcs are overlapping

2316 max_nd_val, _ = max_nd(arcs)

2317 if max_nd_val > 1:

2318 # display warning

2319 logging.log(logging.WARN, f"{max_nd_val - 1} arcs are overlapping")

2320

2321 # force the output to be binary

2322 arcs[arcs != 0] = 1

2323

2324 return arcs, arcs_labelled

2325

2326

2327@typechecker

2328def make_sphere(

2329 grid_size: Vector, radius: Union[float, int], plot_sphere: bool = False, binary: bool = False

2330) -> Union[kt.NP_ARRAY_INT_3D, kt.NP_ARRAY_BOOL_3D]:

2331 """

2332 Generates a sphere mask for a 3D grid given the dimensions of the grid, the radius of the sphere,

2333 and optional flags to plot the sphere and/or return a binary mask.

2334

2335 Args:

2336 grid_size: The size of the grid (assumed to be cubic).

2337 radius: The radius of the sphere.

2338 plot_sphere: Whether to plot the sphere (default: False).

2339 binary: Whether to return a binary mask (default: False).

2340

2341 Returns:

2342 sphere: The sphere mask as a NumPy array.

2343 """

2344 assert len(grid_size) == 3, "grid_size must be a 3D vector"

2345

2346 # enforce a centered sphere

2347 center = np.floor(grid_size / 2).astype(int) + 1

2348

2349 # preallocate the storage variable

2350 if binary:

2351 sphere = np.zeros(grid_size, dtype=bool)

2352 else:

2353 sphere = np.zeros(grid_size, dtype=int)

2354

2355 # create a guide circle from which the individual radii can be extracted

2356 guide_circle = make_circle(np.flip(grid_size[:2]), np.flip(center[:2]), radius)

2357

2358 # step through the guide circle points and create partially filled discs

2359 centerpoints = np.arange(center.x - radius, center.x + 1)

2360 reflection_offset = np.arange(len(centerpoints), 1, -1)

2361 for centerpoint_index in range(len(centerpoints)):

2362 # extract the current row from the guide circle

2363 row_data = guide_circle[:, centerpoints[centerpoint_index] - 1]

2364

2365 # add an index to the grid points in the current row

2366 row_index = row_data * np.arange(1, len(row_data) + 1)

2367

2368 # calculate the radius

2369 swept_radius = (row_index.max() - row_index[row_index != 0].min()) / 2

2370

2371 # create a circle to add to the sphere

2372 circle = make_circle(grid_size[1:], center[1:], swept_radius)

2373

2374 # make an empty fill matrix

2375 if binary:

2376 circle_fill = np.zeros(grid_size[1:], dtype=bool)

2377 else:

2378 circle_fill = np.zeros(grid_size[1:])

2379

2380 # fill in the circle line by line

2381 fill_centerpoints = np.arange(center.z - swept_radius, center.z + swept_radius + 1).astype(int)

2382 for fill_centerpoints_i in fill_centerpoints:

2383 # extract the first row

2384 row_data = circle[:, fill_centerpoints_i - 1]

2385

2386 # add an index to the grid points in the current row

2387 row_index = row_data * np.arange(1, len(row_data) + 1)

2388

2389 # calculate the diameter

2390 start_index = row_index[row_index != 0].min()

2391 stop_index = row_index.max()

2392

2393 # count how many points on the line

2394 num_points = sum(row_data)

2395

2396 # fill in the line

2397 if start_index != stop_index and (stop_index - start_index) >= num_points:

2398 circle_fill[(start_index + num_points // 2) - 1 : stop_index - (num_points // 2), fill_centerpoints_i - 1] = 1

2399

2400 # remove points from the filled circle that existed in the previous

2401 # layer

2402 if centerpoint_index == 0:

2403 sphere[centerpoints[centerpoint_index] - 1, :, :] = circle + circle_fill

2404 prev_circle = circle + circle_fill

2405 else:

2406 prev_circle_alt = circle + circle_fill

2407 circle_fill = circle_fill - prev_circle

2408 circle_fill[circle_fill < 0] = 0

2409 sphere[centerpoints[centerpoint_index] - 1, :, :] = circle + circle_fill

2410 prev_circle = prev_circle_alt

2411

2412 # create the other half of the sphere at the same time

2413 if centerpoint_index != len(centerpoints) - 1:

2414 sphere[center.x + reflection_offset[centerpoint_index] - 2, :, :] = sphere[centerpoints[centerpoint_index] - 1, :, :]

2415

2416 # plot results

2417 if plot_sphere:

2418 raise NotImplementedError

2419 return sphere

2420

2421

2422@typechecker

2423def make_spherical_section(

2424 radius: Real[kt.ScalarLike, ""],

2425 height: Real[kt.ScalarLike, ""],

2426 width: Optional[Real[kt.ScalarLike, ""]] = None,

2427 plot_section: bool = False,

2428 binary: bool = False,

2429) -> Tuple:

2430 """

2431 Generates a spherical section mask given the radius and height of the section and

2432 optional parameters to specify the width and/or plot and return a binary mask.

2433

2434 Args:

2435 radius: The radius of the spherical section.

2436 height: The height of the spherical section.

2437 width: The width of the spherical section (default: height).

2438 plot_section: Whether to plot the spherical section (default: False).

2439 binary: Whether to return a binary mask (default: False).

2440

2441 Returns:

2442 ss: The spherical section mask as a NumPy array.

2443 dist_map: The distance map of the spherical section mask as a NumPy array.

2444

2445 Raises:

2446 ValueError: If the width is not an odd number.

2447 NotImplementedError: Plotting not currently supported.

2448 """

2449 use_spherical_sections = True

2450

2451 # force inputs to be integers

2452 radius = int(radius)

2453 height = int(height)

2454

2455 use_width = width is not None

2456 if use_width:

2457 width = int(width)

2458 if width % 2 == 0:

2459 raise ValueError("Input width must be an odd number.")

2460

2461 # calculate minimum grid dimensions to fit entire sphere

2462 Nx = 2 * radius + 1

2463

2464 # create sphere

2465 ss = make_sphere(Vector([Nx] * 3), radius, False, binary)

2466

2467 # truncate to given height

2468 if use_spherical_sections:

2469 ss = ss[:height, :, :]

2470 else:

2471 ss = np.transpose(ss[:, :height, :], [1, 2, 0])

2472

2473 # flatten transducer and store the maximum and indices

2474 mx = np.squeeze(np.max(ss, axis=0))

2475

2476 # calculate the total length/width of the transducer

2477 length = mx[(len(mx) + 1) // 2].sum()

2478

2479 # truncate transducer grid based on length (removes empty rows and columns)

2480 offset = int((Nx - length) / 2)

2481 ss = ss[:, offset:-offset, offset:-offset]

2482

2483 # also truncate to given width if defined by user

2484 if use_width:

2485 # check the value is appropriate

2486 if width > length:

2487 raise ValueError("Input for width must be less than or equal to transducer length.")

2488

2489 # calculate offset

2490 offset = int((length - width) / 2)

2491

2492 # truncate transducer grid

2493 ss = ss[:, offset:-offset, :]

2494

2495 # compute average distance between each grid point and its contiguous

2496

2497 # calculate x-index of each grid point in the spherical section, create

2498 # mask and remove singleton dimensions

2499 mx, mx_ind = np.max(ss, axis=0), ss.argmax(axis=0) + 1

2500 mask = np.squeeze(mx != 0)

2501 mx_ind = np.squeeze(mx_ind) * mask

2502

2503 # double check there there is only one value of spherical section in

2504 # each matrix column

2505 if mx.sum() != ss.sum():

2506 raise ValueError("mean neighbour distance cannot be calculated uniquely due to overlapping points in the x-direction")

2507

2508 # calculate average distance to grid point neighbours in the flat case

2509 x_dist = np.tile([1, 0, 1], [3, 1])

2510 y_dist = x_dist.T

2511 flat_dist = np.sqrt(x_dist**2 + y_dist**2)

2512 flat_dist = np.mean(flat_dist)

2513

2514 # compute distance map

2515 dist_map = np.zeros(mx_ind.shape)

2516 sz = mx_ind.shape

2517 for m in range(sz[0]):

2518 for n in range(sz[1]):

2519 # clear map

2520 local_heights = np.zeros((3, 3))

2521

2522 # extract the height (x-distance) of the 8 neighbouring grid

2523 # points

2524 if m == 0 and n == 0:

2525 local_heights[1:3, 1:3] = mx_ind[m : m + 2, n : n + 2]

2526 elif m == (sz[0] - 1) and n == (sz[1] - 1):

2527 local_heights[0:2, 0:2] = mx_ind[m - 1 : m + 1, n - 1 : n + 1]

2528 elif m == 0 and n == (sz[1] - 1):

2529 local_heights[1:3, 0:2] = mx_ind[m : m + 2, n - 1 : n + 1]

2530 elif m == (sz[0] - 1) and n == 0:

2531 local_heights[0:2, 1:3] = mx_ind[m - 1 : m + 1, n : n + 2]

2532 elif m == 0:

2533 local_heights[1:3, :] = mx_ind[m : m + 2, n - 1 : n + 2]

2534 elif m == (sz[0] - 1):

2535 local_heights[0:2, :] = mx_ind[m - 1 : m + 1, n - 1 : n + 2]

2536 elif n == 0:

2537 local_heights[:, 1:3] = mx_ind[m - 1 : m + 2, n : n + 2]

2538 elif n == (sz[1] - 1):

2539 local_heights[:, 0:2] = mx_ind[m - 1 : m + 2, n - 1 : n + 1]

2540 else:

2541 local_heights = mx_ind[m - 1 : m + 2, n - 1 : n + 2]

2542

2543 # compute average variation from center

2544 local_heights_var = abs(local_heights - local_heights[1, 1])

2545

2546 # threshold no neighbours

2547 local_heights_var[local_heights == 0] = 0

2548

2549 # calculate total distance from centre

2550 dist = np.sqrt(x_dist**2 + y_dist**2 + local_heights_var**2)

2551

2552 # average and store as a ratio

2553 dist_map[m, n] = 1 + (np.mean(dist) - flat_dist) / flat_dist

2554

2555 # threshold out the non-transducer grid points

2556 dist_map[mask != 1] = 0

2557

2558 # plot if required

2559 if plot_section:

2560 raise NotImplementedError

2561

2562 return ss, dist_map

2563

2564

2565@typechecker

2566def make_cart_rect(

2567 rect_pos,

2568 Lx: Real[kt.ScalarLike, ""],

2569 Ly: Real[kt.ScalarLike, ""],

2570 theta: Optional[Union[Real[kt.ScalarLike, ""], List, Integer[np.ndarray, "..."], Float[np.ndarray, "..."]]] = None,

2571 num_points: int = 0,

2572 plot_rect: bool = False,

2573) -> Union[kt.NP_ARRAY_FLOAT_2D, kt.NP_ARRAY_FLOAT_3D]:

2574 """

2575 Create evenly distributed Cartesian points covering a rectangle.

2576

2577 Args:

2578 rect_pos : List. Cartesian position of the centre of the rectangle.

2579 Lx : Float. Height of the rectangle (along the x-coordinate before rotation).

2580 Ly : Float. Width of the rectangle (along the y-coordinate before rotation).

2581 theta : Float or List (Optional). Specifies the orientation of the rectangle.

2582 num_points : int (Optional). Approximate number of points on the rectangle.

2583 plot_rect : bool (Optional). Boolean controlling whether the Cartesian points are plotted.

2584

2585 Returns:

2586 np.array. 2 x num_points* or 3 x num_points* array of Cartesian coordinates.

2587 """

2588

2589 # Find number of points in along each axis

2590 npts_x = math.ceil(np.sqrt(num_points * Lx / Ly))

2591 npts_y = math.ceil(num_points / npts_x)

2592

2593 # Recalculate the true number of points

2594 num_points = npts_x * npts_y

2595

2596 # Distance between points in each dimension

2597 d_x = 2 / npts_x

2598 d_y = 2 / npts_y

2599

2600 # Compute canonical rectangle points ([-1, 1] x [-1, 1], z=0 plane)

2601 p_x = np.linspace(-1 + d_x / 2, 1 - d_x / 2, npts_x)

2602 p_y = np.linspace(-1 + d_y / 2, 1 - d_y / 2, npts_y)

2603 P_x, P_y = np.meshgrid(p_x, p_y, indexing="ij")

2604 p0 = np.stack((P_x.flatten(), P_y.flatten()), axis=0)

2605

2606 # Add z-dimension points if in 3D

2607 if len(rect_pos) == 3:

2608 p0 = np.vstack((p0, np.zeros(num_points)))

2609

2610 # Transform the canonical rectangle points to give the specified rectangle

2611 if len(rect_pos) == 2:

2612 # Scaling transformation

2613 S = np.array([[Lx, 0], [0, Ly]]) / 2

2614

2615 # Rotation

2616 if theta is None:

2617 R = np.eye(2)

2618 else:

2619 R = np.array([[cosd(theta), -sind(theta)], [sind(theta), cosd(theta)]])

2620

2621 else:

2622 # Scaling transformation

2623 S = np.array([[Lx, 0, 0], [0, Ly, 0], [0, 0, 2]]) / 2

2624

2625 # Rotation

2626 if theta is None:

2627 # No rotation

2628 R = np.eye(3)

2629 else:

2630 # Using intrinsic rotations chain from right to left (xyz rotations)

2631 R = Rotation.from_euler("xyz", theta, degrees=True).as_matrix()

2632

2633 # Combine scaling and rotation matrices

2634 A = np.dot(R, S)

2635

2636 # Apply this transformation to the canonical points

2637 p0 = np.dot(A, p0)

2638

2639 # Shift the rectangle to the appropriate centre

2640 rect = p0 + np.expand_dims(np.array(rect_pos), axis=1)

2641

2642 return rect

2643

2644

2645@typechecker

2646def focused_bowl_oneil(

2647 radius: kt.NUMERIC,

2648 diameter: kt.NUMERIC,

2649 velocity: kt.NUMERIC,

2650 frequency: kt.NUMERIC,

2651 sound_speed: kt.NUMERIC,

2652 density: kt.NUMERIC,

2653 axial_positions: Optional[Union[kt.NP_ARRAY_FLOAT_1D, float, List]] = None,

2654 lateral_positions: Optional[Union[kt.NP_ARRAY_FLOAT_1D, float, List]] = None,

2655) -> Tuple[Optional[kt.NP_ARRAY_FLOAT_1D], Optional[kt.NP_ARRAY_FLOAT_1D], Optional[kt.NP_ARRAY_COMPLEX_1D]]:

2656 """

2657 Calculates O'Neil's solution for the axial and lateral pressure amplitude generated by a focused bowl transducer.

2658

2659 Args:

2660 radius: The radius of the transducer.

2661 diameter: The diameter of the transducer.

2662 velocity: The normal surface velocity of the transducer.

2663 frequency: The driving frequency of the sinusoid.

2664 sound_speed: The sound speed in the medium.

2665 density: The density of the medium.

2666 axial_positions: The positions along the beam axis where the pressure is evaluated

2667 (0 corresponds to the transducer surface). Set to [] to return only lateral pressure.

2668 lateral_positions: The lateral positions through the geometric focus where the pressure is evaluated

2669 (0 corresponds to the beam axis). Set to [] to return only axial pressure.

2670

2671 Returns:

2672 p_axial: pressure amplitude at the positions specified by axial_position [Pa]

2673 p_lateral: pressure amplitude at the positions specified by lateral_position [Pa]

2674 p_axial_complex: complex pressure amplitude at the positions specified by axial_position [Pa]

2675

2676 Example:

2677 # define transducer parameters

2678 radius = 140e-3 # [m]

2679 diameter = 120e-3 # [m]

2680 velocity = 100e-3 # [m / s]

2681 frequency = 1e6 # [Hz]

2682 sound_speed = 1500 # [m / s]

2683 density = 1000 # [kg / m ^ 3]

2684

2685 # define position vectors

2686 axial_position = np.arange(0, 250e-3 + 1e-4, 1e-4) # [m]

2687 lateral_position = np.arange(-15e-3, 15e-3 + 1e-4, 1e-4) # [m]

2688

2689 # evaluate pressure

2690 p_axial, p_lateral, p_axial_complex = focused_bowl_oneil(radius, diameter,

2691 velocity, frequency, sound_speed, density,

2692 axial_position, lateral_position)

2693

2694 References:

2695 O'Neil, H. (1949). Theory of focusing radiators. J. Acoust. Soc. Am., 21(5), 516-526.

2696 """

2697

2698 float_eps = np.finfo(float).eps

2699

2700 # @typechecker => could not figure out what's wrong with type annotation here, revisit in the future

2701 def calculate_axial_pressure() -> Tuple[Float[np.ndarray, "N"], Complex[np.ndarray, "N"]]:

2702 # calculate distances

2703 B = np.sqrt((axial_positions - h) ** 2 + (diameter / 2) ** 2)

2704 d = B - axial_positions

2705 E = 2 / (1 - axial_positions / radius)

2706 M = (B + axial_positions) / 2

2707

2708 # compute pressure

2709 P = E * np.sin(k * d / 2)

2710

2711 # replace values where axial_position is equal to the radius with limit

2712 P[np.abs(axial_positions - radius) < float_eps] = k * h

2713

2714 # calculate magnitude of the on - axis pressure (Eq. 3.0)

2715 axial_pressure = density * sound_speed * velocity * np.abs(P)

2716 # calculate complex magnitude of the on - axis pressure assuming t = 0 (Eq. 3.0)

2717 complex_axial_pressure = density * sound_speed * velocity * P * 1j * np.exp(-1j * k * M)

2718

2719 return axial_pressure, complex_axial_pressure

2720

2721 @typechecker

2722 def calculate_lateral_pressure() -> Float[np.ndarray, "N"]:

2723 # calculate magnitude of the lateral pressure at the geometric focus

2724 Z = k * lateral_positions * diameter / (2 * radius)

2725 # TODO: this should work

2726 # assert np.all(Z) > 0, 'Z must be greater than 0'

2727 lateral_pressure = 2.0 * density * sound_speed * velocity * k * h * scipy.special.jv(1, Z) / Z

2728

2729 # replace origin with limit

2730 lateral_pressure[lateral_positions == 0] = density * sound_speed * velocity * k * h

2731 return lateral_pressure

2732

2733 # wave number

2734 k = 2 * np.pi * frequency / sound_speed

2735

2736 # height of rim

2737 h = radius - np.sqrt(radius**2 - (diameter / 2) ** 2)

2738

2739 p_axial = None

2740 p_lateral = None

2741 p_axial_complex = None

2742

2743 if lateral_positions is not None:

2744 p_lateral = calculate_lateral_pressure()

2745 if axial_positions is not None:

2746 p_axial, p_axial_complex = calculate_axial_pressure()

2747 return p_axial, p_lateral, p_axial_complex

2748

2749

2750@typechecker

2751def focused_annulus_oneil(

2752 radius: float,

2753 diameter: Union[Float[np.ndarray, "NumElements 2"], Float[np.ndarray, "2 NumElements"]],

2754 amplitude: Float[np.ndarray, "NumElements"],

2755 phase: Float[np.ndarray, "NumElements"],

2756 frequency: kt.NUMERIC,

2757 sound_speed: kt.NUMERIC,

2758 density: kt.NUMERIC,

2759 axial_positions: Union[kt.NP_ARRAY_FLOAT_1D, float, list],

2760) -> Union[kt.NP_ARRAY_FLOAT_1D, float]:

2761 """Compute axial pressure for focused annulus transducer using O'Neil's solution

2762

2763 focused_annulus_oneil calculates the axial pressure for a focused

2764 annular transducer using O'Neil's solution (O'Neil, H. Theory of

2765 focusing radiators. J. Acoust. Soc. Am., 21(5), 516-526, 1949). The

2766 annuluar elements are uniformly driven by a continuous wave sinusoid

2767 at a given frequency and normal surface velocity.

2768

2769 The solution is evaluated at the positions (along the beam axis) given

2770 by axial_position. Where 0 corresponds to the transducer surface.

2771

2772 Note, O'Neil's formulae are derived under the assumptions of the

2773 Rayleigh integral, which are valid when the transducer diameter is

2774 large compared to both the transducer height and the acoustic

2775 wavelength.

2776

2777 Example:

2778 # define transducer parameters

2779 radius = 140e-3 # [m]

2780 diameter = 120e-3 # [m]

2781 velocity = 100e-3 # [m / s]

2782 frequency = 1e6 # [Hz]

2783 sound_speed = 1500 # [m / s]

2784 density = 1000 # [kg / m^3]

2785

2786 # define position vectors

2787 axial_position = np.arange(0, 250e-3 + 1e-4, 1e-4) # [m]

2788 p_axial = focused_annulus_oneil(radius, diameter, amplitude, phase, frequency, sound_speed, density, axial_position)

2789

2790 Args:

2791 radius: transducer radius of curvature [m]

2792 diameter: 2 x num_elements array containing pairs of inner and outer aperture diameter

2793 (diameter of opening) [m].

2794 amplitude: array containing the normal surface velocities for each element [m/s]

2795 phase: array containing the phase for each element [rad]

2796 frequency: driving frequency [Hz]

2797 sound_speed: speed of sound in the propagating medium [m/s]

2798 density: density in the propagating medium [kg/m^3]

2799 axial_positions: vector of positions along the beam axis where the

2800 pressure amplitude is calculated [m]

2801

2802 Returns:

2803 p_axial: pressure amplitude at the positions specified by axial_position [Pa]

2804

2805 References:

2806 O'Neil, H. (1949). Theory of focusing radiators. J. Acoust. Soc. Am., 21(5), 516-526.

2807

2808 See also focused_bowl_oneil.

2809 """

2810

2811 if not np.greater_equal(diameter, np.zeros_like(diameter)).all() or not np.isreal(diameter).all() or not np.isfinite(diameter).all():

2812 raise ValueError("wrong values in diameter object")

2813

2814 if not np.all(np.isfinite(amplitude)):

2815 raise ValueError("amplitude contains an np.inf")

2816 if not np.all(np.isfinite(frequency)):

2817 raise ValueError("frequency contains an np.inf")

2818

2819 # set the number of elements in annular array

2820 num_elements: int = np.size(amplitude)

2821

2822 if (radius <= 0.0) or not np.isreal(radius) or not np.isfinite(radius):

2823 raise ValueError("radius is incorrect")

2824

2825 if ((phase < -np.pi).any() or (phase > np.pi).any()) or not np.isreal(phase).any() or not np.isfinite(phase).all():

2826 raise ValueError("phase is incorrect")

2827

2828 if np.shape(diameter)[0] != 2:

2829 diameter = np.transpose(diameter)

2830

2831 # pre-allocate output

2832 p_axial = np.zeros(np.shape(axial_positions))

2833

2834 # loop over elements and sum fields

2835 for ind in range(num_elements):

2836 # get complex pressure for bowls with inner and outer aperture diameter

2837 if diameter[0, ind] == 0:

2838 p_el_inner = 0.0 + 0.0j

2839 else:

2840 _, _, p_el_inner = focused_bowl_oneil(

2841 radius, diameter[0, ind], amplitude[ind], frequency, sound_speed, density, axial_positions=axial_positions

2842 )

2843

2844 _, _, p_el_outer = focused_bowl_oneil(

2845 radius, diameter[1, ind], amplitude[ind], frequency, sound_speed, density, axial_positions=axial_positions

2846 )

2847

2848 # pressure for annular element

2849 p_el = p_el_outer - p_el_inner

2850

2851 # account for phase

2852 p_el = np.abs(p_el) * np.exp(1.0j * (np.angle(p_el) + phase[ind]))

2853

2854 # add to complete response

2855 p_axial = p_axial + p_el

2856

2857 # take magnitude of complete response

2858 return np.abs(p_axial)

2859

2860

2861def ndgrid(*args):

2862 return np.array(np.meshgrid(*args, indexing="ij"))

2863

2864

2865def trim_cart_points(kgrid, points: np.ndarray):

2866 """

2867 trim_cart_points filters a dim x num_points array of Cartesian points

2868 so that only those within the bounds of a given kgrid remain.

2869 :param kgrid: Object of the kWaveGrid class defining the Cartesian

2870 and k-space grid fields.

2871 :param points: dim x num_points array of Cartesian coordinates to trim [m].

2872 :return: dim x num_points array of Cartesian coordinates that lie within the grid defined by kgrid [m].

2873 """

2874

2875 # find indices for points within the simulation domain

2876 ind_x = (points[0, :] >= kgrid.x_vec[0]) & (points[0, :] <= kgrid.x_vec[-1])

2877

2878 if kgrid.dim > 1:

2879 ind_y = (points[1, :] >= kgrid.y_vec[0]) & (points[1, :] <= kgrid.y_vec[-1])

2880

2881 if kgrid.dim > 2:

2882 ind_z = (points[2, :] >= kgrid.z_vec[0]) & (points[2, :] <= kgrid.z_vec[-1])

2883

2884 # combine indices

2885 if kgrid.dim == 1:

2886 ind = ind_x

2887 elif kgrid.dim == 2:

2888 ind = ind_x & ind_y

2889 elif kgrid.dim == 3:

2890 ind = ind_x & ind_y & ind_z

2891

2892 # output only valid points

2893 points = points[:, ind]

2894

2895 return points

2896

2897

2898@typechecker

2899def make_cart_arc(

2900 arc_pos: Float[np.ndarray, "2"],

2901 radius: Union[float, int],

2902 diameter: Union[float, int],

2903 focus_pos: Float[np.ndarray, "2"],

2904 num_points: int,

2905 plot_arc: bool = False,

2906) -> Float[np.ndarray, "2 NumPoints"]:

2907 """

2908 make_cart_arc creates a 2 x num_points array of the Cartesian

2909 coordinates of points evenly distributed over an arc. The midpoint of

2910 the arc is set by arc_pos. The orientation of the arc is set by

2911 focus_pos, which corresponds to any point on the axis of the arc

2912 (note, this must not be equal to arc_pos). It is assumed that the arc

2913 angle is equal to or less than pi radians. If the radius is set to

2914 inf, a line is generated.

2915

2916 Note, the first and last points do not lie on the end-points of the

2917 underlying canonical arc, but are offset by half the angular spacing

2918 between the points.

2919

2920 Args:

2921 arc_pos: center of the rear surface (midpoint) of the arc given as a two element vector [bx, by] [m].

2922 radius: radius of curvature of the arc [m].

2923 diameter: diameter of the opening of the arc (length of line connecting arc endpoints) [m].

2924 focus_pos: Any point on the beam axis of the arc given as a two element vector [fx, fy] [m].

2925 num_points: number of points to generate along the arc.

2926 plot_arc: boolean flag to plot the arc.

2927

2928 Returns:

2929 2 x num_points array of Cartesian coordinates of points along the arc [m].

2930

2931 """

2932

2933 # Check input values

2934 if radius <= 0:

2935 raise ValueError("The radius must be positive.")

2936 if diameter <= 0:

2937 raise ValueError("The diameter must be positive.")

2938 if diameter > 2 * radius:

2939 raise ValueError("The diameter of the arc must be less than twice the radius of curvature.")

2940 if np.all(arc_pos == focus_pos):

2941 raise ValueError("The focus_pos must be different from the arc_pos.")

2942

2943 # Check for infinite radius of curvature, and make a large finite number

2944 if np.isinf(radius):

2945 radius = 1e10 * diameter

2946

2947 # Compute arc angle from chord

2948 varphi_max = np.arcsin(diameter / (2 * radius))

2949

2950 # Angle between points

2951 dvarphi = 2 * varphi_max / num_points

2952

2953 # Compute canonical arc points where arc is centered on the origin, and its

2954 # back is placed at a distance of radius along the positive y-axis

2955 t = np.linspace(-varphi_max + dvarphi / 2, varphi_max - dvarphi / 2, num_points)

2956 p0 = radius * np.array([np.sin(t), np.cos(t)])

2957

2958 # Linearly transform canonical points to give arc in correct orientation

2959 R, b = compute_linear_transform2D(arc_pos, radius, focus_pos)

2960 arc = np.asarray(np.dot(R, p0) + b)

2961

2962 # Plot results

2963 if plot_arc:

2964 # Select suitable axis scaling factor

2965 _, scale, prefix, _ = scale_SI(np.max(np.abs(arc)))

2966

2967 # Create the figure

2968 plt.figure()

2969 plt.plot(arc[1, :] * scale, arc[0, :] * scale, "b.")

2970 plt.gca().invert_yaxis()

2971 plt.xlabel(f"y-position [{prefix}m]")

2972 plt.ylabel(f"x-position [{prefix}m]")

2973 plt.axis("equal")

2974 plt.show()

2975

2976 return arc

2977

2978

2979def compute_linear_transform2D(arc_pos: Vector, radius: float, focus_pos: Vector) -> Tuple[np.ndarray, np.ndarray]:

2980 """

2981

2982 Compute a rotation matrix to transform the computed arc points to the orientation

2983 specified by the arc and focus positions

2984

2985 Args:

2986 arc_pos:

2987 radius:

2988 focus_pos:

2989

2990 Returns:

2991 The rotation matrix R and an offset b

2992

2993 """

2994

2995 # Vector pointing from arc_pos to focus_pos

2996 beam_vec = focus_pos - arc_pos

2997

2998 # Normalize to give unit beam vector

2999 beam_vec = beam_vec / np.linalg.norm(beam_vec)

3000

3001 # Canonical normalized beam_vec (canonical arc_pos is [0, 1])

3002 beam_vec0 = np.array([0, -1])

3003

3004 # Find the angle between canonical and specified beam_vec

3005 theta = np.arctan2(beam_vec[1], beam_vec[0]) - np.arctan2(beam_vec0[1], beam_vec0[0])

3006

3007 # Convert to a rotation matrix

3008 R = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])

3009

3010 # Compute an offset for the arc, where arc_centre = move from arc_pos

3011 # towards focus by radius

3012 b = arc_pos.reshape(-1, 1) + radius * beam_vec.reshape(-1, 1)

3013

3014 return R, b

3015

3016

3017@typechecker

3018def make_cart_spherical_segment(

3019 bowl_pos: Float[np.ndarray, "3"],

3020 radius: Union[float, int],

3021 inner_diameter: Union[float, int],

3022 outer_diameter: Union[float, int],

3023 focus_pos: Float[np.ndarray, "3"],

3024 num_points: int,

3025 plot_bowl: Optional[bool] = False,

3026 num_points_inner: int = 0,

3027) -> Float[np.ndarray, "3 NumPoints"]:

3028 """

3029 Create evenly distributed Cartesian points covering a spherical segment.

3030

3031 Args:

3032 bowl_pos: Cartesian position of the centre of the rear surface

3033 of the underlying bowl on which the spherical segment lies given as a

3034 three element vector [bx, by, bz] [m].

3035 radius: Radius of curvature of the underlying bowl [m].

3036 inner_diameter: Inner aperture diameter of the spherical segment [m].

3037 outer_diameter: Outer aperture diameter of the spherical segment [m].

3038 focus_pos: Any point on the beam axis of the underlying bowl

3039 given as a three element vector [fx, fy, fz] [m].

3040 num_points: Number of points on the spherical segment.

3041 plot_bowl: Boolean controlling whether the Cartesian points are plotted.

3042 num_points_inner: If constructing an annular array with contiguous

3043 elements (no kerf), the positions of the points will not exactly match

3044 makeCartBowl, as each element has no knowledge of the number of points on

3045 the internal elements. To force the points to match, specify the total number

3046 of points used on all annular segments internal to the current one.

3047

3048 Returns:

3049 numpy.ndarray: 3 x num_points array of Cartesian coordinates.

3050 """

3051

3052 # check input values

3053 if radius <= 0:

3054 raise ValueError("The radius must be positive.")

3055 if inner_diameter < 0:

3056 raise ValueError("The inner diameter must be positive.")

3057 if inner_diameter >= outer_diameter:

3058 raise ValueError("The inner diameter must be less than the outer diameter.")

3059 if outer_diameter <= 0:

3060 raise ValueError("The outer diameter must be positive.")

3061 if outer_diameter > 2 * radius:

3062 raise ValueError("The outer diameter of the bowl must be equal or less than twice the radius of curvature.")

3063 if np.all(bowl_pos == focus_pos):

3064 raise ValueError("The focus_pos must be different from the bowl_pos.")

3065

3066 # check for infinite radius of curvature

3067 if np.isinf(radius):

3068 raise ValueError("Annular disc (infinite radius) not yet supported.")

3069

3070 # compute arc angle from chord

3071 varphi_min = np.arcsin(inner_diameter / (2 * radius))

3072 varphi_max = np.arcsin(outer_diameter / (2 * radius))

3073

3074 # compute spiral parameters over annulus

3075 theta = lambda t: GOLDEN_ANGLE * t

3076 if num_points_inner > 0:

3077 C = (1 - np.cos(varphi_max)) / (num_points + num_points_inner - 1)

3078 varphi = lambda t: np.arccos(1 - C * t)

3079 t_start = int(np.ceil((1 - np.cos(varphi_min)) / C))

3080 t = np.linspace(t_start, num_points_inner + num_points - 1, num_points)

3081 else:

3082 C = (1 - np.cos(varphi_max)) / (num_points - 1)

3083 varphi = lambda t: np.arccos(1 - C * t)

3084 t_start = int(np.ceil((1 - np.cos(varphi_min)) / C))

3085 t = np.linspace(t_start, num_points - 1, num_points)

3086

3087 # compute canonical spiral points

3088 p0 = np.array([np.cos(theta(t)) * np.sin(varphi(t)), np.sin(theta(t)) * np.sin(varphi(t)), np.cos(varphi(t))])

3089 p0 = radius * p0

3090

3091 # linearly transform the canonical spiral points to give bowl in correct orientation

3092 R, b = compute_linear_transform(bowl_pos, focus_pos, radius)

3093 if np.shape(b) == (3,):

3094 b = np.expand_dims(b, axis=1) # expand dims for broadcasting

3095

3096 segment = R @ p0 + b

3097

3098 # plot results

3099 if plot_bowl is True:

3100 _, scale, prefix, unit = scale_SI(np.max(segment))

3101

3102 # create the figure

3103 fig = plt.figure()

3104 ax = fig.add_subplot(111, projection="3d")

3105 ax.scatter(segment[0] * scale, segment[1] * scale, segment[2] * scale)

3106 ax.set_xlabel("[" + prefix + "m]")

3107 ax.set_ylabel("[" + prefix + "m]")

3108 ax.set_zlabel("[" + prefix + "m]")

3109 ax.set_box_aspect([1, 1, 1])

3110 plt.grid(True)

3111 plt.show()

3112

3113 return segment