From a35d06ca798ef3cfc64b5bdd357b964026cd1dd0 Mon Sep 17 00:00:00 2001 From: Walter Simson Date: Mon, 18 Dec 2023 14:10:52 -0800 Subject: [PATCH] update b-mode reconstruction example --- examples/bmode_reconstruction_example.py | 120 ++++++++++++++++++++--- kwave/ktransducer.py | 71 ++++++++------ kwave/reconstruction/beamform.py | 1 + kwave/reconstruction/tools.py | 2 +- kwave/utils/filters.py | 6 +- 5 files changed, 154 insertions(+), 46 deletions(-) diff --git a/examples/bmode_reconstruction_example.py b/examples/bmode_reconstruction_example.py index b52318b1..2d30edcd 100644 --- a/examples/bmode_reconstruction_example.py +++ b/examples/bmode_reconstruction_example.py @@ -13,10 +13,12 @@ from kwave.ktransducer import NotATransducer, kWaveTransducerSimple from kwave.options.simulation_execution_options import SimulationExecutionOptions from kwave.options.simulation_options import SimulationOptions -from kwave.reconstruction.beamform import beamform -from kwave.reconstruction.converter import build_channel_data from kwave.utils.dotdictionary import dotdict -from kwave.utils.signals import tone_burst +from kwave.utils.signals import tone_burst, get_win +from kwave.utils.filters import gaussian_filter +from kwave.reconstruction.tools import log_compression +from kwave.reconstruction.beamform import scan_conversion, envelope_detection + def main(): @@ -99,7 +101,7 @@ def main(): logging.log(logging.INFO, f"RUN_SIMULATION set to {RUN_SIMULATION}") # preallocate the storage set medium position - scan_lines = np.zeros((number_scan_lines, not_transducer.number_active_elements, kgrid.Nt)) + scan_lines = np.zeros((number_scan_lines, kgrid.Nt)) medium_position = 0 for scan_line_index in range(0, number_scan_lines): @@ -133,14 +135,14 @@ def main(): execution_options=SimulationExecutionOptions(is_gpu_simulation=True) ) - scan_lines[scan_line_index, :] = not_transducer.combine_sensor_data(sensor_data['p'].T) + scan_lines[scan_line_index, :] = not_transducer.scan_line(not_transducer.combine_sensor_data(sensor_data['p'].T)) # update medium position medium_position = medium_position + transducer.element_width if RUN_SIMULATION: simulation_data = scan_lines - # scipy.io.savemat('sensor_data.mat', {'sensor_data_all_lines': simulation_data}) + scipy.io.savemat('sensor_data.mat', {'sensor_data_all_lines': simulation_data}) else: logging.log(logging.INFO, "Downloading data from remote server...") @@ -150,17 +152,103 @@ def main(): simulation_data = scipy.io.loadmat(sensor_data_path)['sensor_data_all_lines'] - # temporary fix for dimensionality - simulation_data = simulation_data[None, :] - - sampling_frequency = 2.772e+07 - prf = 1e4 - focal_depth = 20e-3 # [m] - channel_data = build_channel_data(simulation_data, kgrid, not_transducer, - sampling_frequency, prf, focal_depth) + # temporary fix for dimensionality + import matplotlib.pyplot as plt + + scan_lines = simulation_data + + + # ========================================================================= + # PROCESS THE RESULTS + # ========================================================================= + + # ----------------------------- + # Remove Input Signal + # ----------------------------- + # Trim the delay offset from the scan line data + t0_offset = 85 # round(len(input_signal) / 2) + (not_transducer.stored_appended_zeros - not_transducer.beamforming_delays_offset) + scan_lines = scan_lines[:, t0_offset:] + + # Get the new length of the scan lines + Nt = len(scan_lines[0, :]) + + # TODO: fix this + scan_line_win, cg = get_win(N= Nt* 2,type_='Tukey', plot_win=False,param=0.05) + # scan_line_win = np.concatenate((np.zeros([1, t0_offset * 2]), scan_line_win.T[:, :Nt - t0_offset * 2]), axis=1) + + # ----------------------------- + # Time Gain Compensation + # ----------------------------- + + # Create radius variable + r = c0 * np.arange(1, Nt + 1) * kgrid.dt / 2 + + # Define absorption value and convert to correct units + tgc_alpha_db_cm = medium.alpha_coeff * (tone_burst_freq * 1e-6)**medium.alpha_power + tgc_alpha_np_m = tgc_alpha_db_cm / 8.686 * 100 + + # Create time gain compensation function + tgc = np.exp(tgc_alpha_np_m * 2 * r) + + # Apply the time gain compensation to each of the scan lines + scan_lines *= tgc + + # ----------------------------- + # Frequency Filtering + # ----------------------------- + + scan_lines_fund = gaussian_filter(scan_lines.T, 1/kgrid.dt, tone_burst_freq, 100) + scan_lines_harm = gaussian_filter(scan_lines.T, 1/kgrid.dt, 2 * tone_burst_freq, 30) # plotting was not impl. + + # ----------------------------- + # Envelope Detection + # ----------------------------- + + scan_lines_fund = envelope_detection(scan_lines_fund) + scan_lines_harm = envelope_detection(scan_lines_harm) + + # ----------------------------- + # Log Compression + # ----------------------------- + + compression_ratio = 3 + scan_lines_fund = log_compression(scan_lines_fund, compression_ratio, True) + scan_lines_harm = log_compression(scan_lines_harm, compression_ratio, True) + + # ========================================================================= + # VISUALISATION + # ========================================================================= + + # Set the desired size of the image + image_size = kgrid.size + + # Create the axis variables + x_axis = [0, image_size[0] * 1e3] # [mm] + y_axis = [-0.5 * image_size[0] * 1e3, 0.5 * image_size[1] * 1e3] # [mm] + # Plot the data before and after scan conversion + plt.figure() + # plot the sound speed map + plt.subplot(1, 3, 1) + plt.imshow(sound_speed_map[int(t0_offset/4):, 64:-64, int(grid_size_points.z / 2)], aspect='auto', + extent=[y_axis[0], y_axis[1], x_axis[1], x_axis[0]]) + plt.title('Sound Speed Map') + plt.xlabel('Image width [mm]') + plt.ylabel('Depth [mm]') + plt.subplot(1, 3, 2) + plt.imshow(scan_lines_fund, cmap='bone', aspect='auto', extent=[y_axis[0], y_axis[1], x_axis[1], x_axis[0]]) + plt.xlabel('Image width [mm]') + plt.title('Fundamental') + plt.yticks([]) + plt.subplot(1, 3, 3) + plt.imshow(scan_lines_harm, cmap='bone', aspect='auto', + extent=[y_axis[0], y_axis[1], x_axis[1], x_axis[0]]) + plt.yticks([]) + plt.xlabel('Image width [mm]') + plt.title('Harmonic') + + plt.show() + # plt.title('Raw Scan - logging.log(logging.INFO, "Beamforming channel data and reconstructing the image...") - beamform(channel_data) if __name__ == '__main__': diff --git a/kwave/ktransducer.py b/kwave/ktransducer.py index 2375b8b7..afca7089 100644 --- a/kwave/ktransducer.py +++ b/kwave/ktransducer.py @@ -670,34 +670,49 @@ def elevation_beamforming_delays(self): delay_times = np.zeros((1, self.transducer.element_length)) return delay_times - # # function to create a scan line based on the input sensor data and the current apodization and beamforming setting - # def scan_line(obj, sensor_data): - # - # # get the current apodization setting - # apodization = obj.get_receive_apodization - # - # # get the current beamforming weights and reverse - # delays = -obj.beamforming_delays - # - # # offset the received sensor_data by the beamforming delays and apply receive apodization - # for element_index in np.arange(0, obj.number_active_elements): - # if delays(element_index) > 0: - # # shift element data forwards - # sensor_data[element_index, :] = apodization[element_index] * [ - # sensor_data(element_index, 1 + delays(element_index):end), - # zeros(1, delays(element_index)) - # ]; - # - # elif delays(element_index) < 0 - # # shift element data backwards - # sensor_data[element_index, :] = apodization[element_index] *[ - # zeros(1, -delays(element_index)), - # sensor_data(element_index, 1:end + delays(element_index)) - # ]; - # - # # form the a-line summing across the elements - # line = sum(sensor_data) - # return lin + def get_receive_apodization(self): + """ + Get the current receive apodization setting. + """ + # Example implementation, adjust based on actual logic + if is_number(self.receive_apodization): + assert self.receive_apodization.size == self.number_active_elements, \ + 'The length of the receive apodization input must match the number of active elements' + return self.receive_apodization + else: + if self.number_active_elements > 1: + apodization, _ = get_win(int(self.number_active_elements), type_=self.receive_apodization) + else: + apodization = 1 + return np.array(apodization) + + def scan_line(self, sensor_data): + """ + Apply beamforming and apodization to the sensor data. + """ + # Get the current apodization setting + apodization = self.get_receive_apodization() + + # Get the current beamforming weights and reverse + delays = -self.beamforming_delays + + # Offset the received sensor_data by the beamforming delays and apply receive apodization + for element_index in range(self.number_active_elements): + if delays[element_index] > 0: + # Shift element data forwards + sensor_data[element_index, :] = np.pad(sensor_data[element_index, delays[element_index]:], + (0, delays[element_index]), + 'constant') * apodization[element_index] + elif delays[element_index] < 0: + # Shift element data backwards + sensor_data[element_index, :] = np.pad(sensor_data[element_index, :sensor_data.shape[1] + delays[element_index]], + (-delays[element_index], 0), + 'constant') * apodization[element_index] + + # Form the line summing across the elements + line = np.sum(sensor_data, axis=0) + return line + def combine_sensor_data(self, sensor_data): # check the data is the correct size diff --git a/kwave/reconstruction/beamform.py b/kwave/reconstruction/beamform.py index a99b5561..1cce9e68 100644 --- a/kwave/reconstruction/beamform.py +++ b/kwave/reconstruction/beamform.py @@ -128,6 +128,7 @@ def beamform(channel_data: ChannelData) -> None: filename = "example_bmode.png" plt.savefig(os.path.join(os.getcwd(), filename)) + plt.plot() logging.log(logging.INFO, f"Plot saved to {os.path.join(os.getcwd(), filename)}") pass diff --git a/kwave/reconstruction/tools.py b/kwave/reconstruction/tools.py index 2736c66b..3101190d 100644 --- a/kwave/reconstruction/tools.py +++ b/kwave/reconstruction/tools.py @@ -14,7 +14,7 @@ def log_compression(signal, cf, normalize=False): Returns: signal: log-compressed signal """ if normalize: - ms = max(signal) + ms = np.max(signal, axis=0) signal = ms * (np.log10(1 + cf * signal / ms) / np.log10(1 + cf)) else: signal = np.log10(1 + cf * signal) / np.log10(1 + cf) diff --git a/kwave/utils/filters.py b/kwave/utils/filters.py index 5f7f900b..660a3dc4 100644 --- a/kwave/utils/filters.py +++ b/kwave/utils/filters.py @@ -397,8 +397,12 @@ def gaussian_filter(signal: Union[np.ndarray, List[float]], # create double-sided Gaussain filter gfilter = np.fmax(gaussian(f, magnitude, mean, variance), gaussian(f, magnitude, -mean, variance)) + # add dimensions to filter to be broadcastable to signal shape + if len(signal.shape) == 2: + gfilter = gfilter[:, np.newaxis] + # apply filter - signal = np.real(ifft(ifftshift(gfilter * fftshift(fft(signal))))) + signal = np.real(ifft(ifftshift(gfilter.T * fftshift(fft(signal.T))))).T return signal