Coverage for kwave/kWaveSimulation_helper/create_absorption_variables.py: 19%

1import logging 

2import math 

3 

4import numpy as np 

5 

6from kwave.kgrid import kWaveGrid 

7from kwave.kmedium import kWaveMedium 

8from kwave.utils.conversion import db2neper 

9 

10 

11def create_absorption_variables(kgrid: kWaveGrid, medium: kWaveMedium, equation_of_state): 

12 # define the lossy derivative operators and proportionality coefficients 

13 if equation_of_state == "absorbing": 

14 return create_absorbing_medium_variables(kgrid.k, medium) 

15 elif equation_of_state == "stokes": 

16 return create_stokes_medium_variables(medium) 

17 else: 

18 raise NotImplementedError 

19 

20 

21def create_absorbing_medium_variables(kgrid_k, medium: kWaveMedium): 

22 # convert the absorption coefficient to nepers.(rad/s)^-y.m^-1 

23 medium.alpha_coeff = db2neper(medium.alpha_coeff, medium.alpha_power) 

24 

25 absorb_nabla1, absorb_tau = get_absorbtion(kgrid_k, medium) 

26 absorb_nabla2, absorb_eta = get_dispersion(kgrid_k, medium) 

27 

28 absorb_nabla1, absorb_nabla2 = apply_alpha_filter(medium, absorb_nabla1, absorb_nabla2) 

29 

30 return absorb_nabla1, absorb_nabla2, absorb_tau, absorb_eta 

31 

32 

33def create_stokes_medium_variables(medium: kWaveMedium): 

34 # convert the absorption coefficient to nepers.(rad/s)^-2.m^-1 

35 medium.alpha_coeff = db2neper(medium.alpha_coeff, 2) 

36 

37 # compute the absorbing coefficient 

38 absorb_tau = compute_absorbing_coeff(medium) 

39 return None, None, absorb_tau, None 

40 

41 

42def get_absorbtion(kgrid_k, medium): 

43 # compute the absorbing fractional Laplacian operator and coefficient 

44 if medium.alpha_mode == "no_absorption": 

45 return 0, 0 

46 

47 nabla1 = kgrid_k ** (medium.alpha_power - 2) 

48 nabla1[np.isinf(nabla1)] = 0 

49 nabla1 = np.fft.ifftshift(nabla1) 

50 tau = compute_absorbing_coeff(medium) 

51 return nabla1, tau 

52 

53 

54def get_dispersion(kgrid_k, medium): 

55 # compute the dispersive fractional Laplacian operator and coefficient 

56 if medium.alpha_mode == "no_dispersion": 

57 return 0, 0 

58 nabla2 = kgrid_k ** (medium.alpha_power - 1) 

59 nabla2[np.isinf(nabla2)] = 0 

60 nabla2 = np.fft.ifftshift(nabla2) 

61 eta = compute_dispersive_coeff(medium) 

62 return nabla2, eta 

63 

64 

65def compute_absorbing_coeff(medium): 

66 tau = -2 * medium.alpha_coeff * (medium.sound_speed ** (medium.alpha_power - 1)) 

67 

68 # modify the sign of the absorption operator if alpha_sign is defined 

69 # (this is used for time-reversal photoacoustic image reconstruction 

70 # with absorption compensation) 

71 if medium.alpha_sign is not None: 

72 tau = np.sign(medium.alpha_sign[0]) * tau 

73 return tau 

74 

75 

76def compute_dispersive_coeff(medium): 

77 eta = 2 * medium.alpha_coeff * medium.sound_speed ** (medium.alpha_power) * math.tan(math.pi * medium.alpha_power / 2) 

78 

79 # modify the sign of the absorption operator if alpha_sign is defined 

80 # (this is used for time-reversal photoacoustic image reconstruction 

81 # with absorption compensation) 

82 if medium.alpha_sign is not None: 

83 eta = np.sign(medium.alpha_sign[1]) * eta 

84 return eta 

85 

86 

87def apply_alpha_filter(medium, nabla1, nabla2): 

88 # pre-filter the absorption parameters if alpha_filter is defined (this 

89 # is used for time-reversal photoacoustic image reconstruction 

90 # with absorption compensation) 

91 if medium.alpha_filter is None: 

92 return nabla1, nabla2 

93 

94 # update command line status 

95 logging.log(logging.INFO, " filtering absorption variables...") 

96 

97 # frequency shift the absorption parameters 

98 nabla1 = np.fft.fftshift(nabla1) 

99 nabla2 = np.fft.fftshift(nabla2) 

100 

101 # apply the filter 

102 nabla1 = nabla1 * medium.alpha_filter 

103 nabla2 = nabla2 * medium.alpha_filter 

104 

105 # shift the parameters back 

106 nabla1 = np.fft.ifftshift(nabla1) 

107 nabla2 = np.fft.ifftshift(nabla2) 

108 return nabla1, nabla2