diff --git a/kwave/reconstruction/beamform.py b/kwave/reconstruction/beamform.py index 6b2487ea..712a4147 100644 --- a/kwave/reconstruction/beamform.py +++ b/kwave/reconstruction/beamform.py @@ -1,139 +1,11 @@ import logging -import os from typing import Tuple, Optional -from warnings import warn - import numpy as np import scipy -from matplotlib import pyplot as plt -from scipy.interpolate import interp1d -from scipy.signal import hilbert -from uff import ChannelData -from uff.position import Position from kwave.utils.conversion import cart2pol from kwave.utils.data import scale_time from kwave.utils.tictoc import TicToc -from .shifted_transform import ShiftedTransform -from .tools import make_time_vector, get_t0, get_origin_array, apodize - - -def beamform(channel_data: ChannelData) -> None: - """ - - Args: - channel_data: shape => (1, 96, 32, 1585) - - Returns: - - """ - warn("This method is no longer supported and will be depricated in v0.4.0 of k-wave-python.", DeprecationWarning, stacklevel=2) - f_number = 1.2 - num_px_z = 256 - imaging_depth = 40e-3 - # apodization_window = 'boxcar' - apodization_window = "none" - number_samples = np.size(channel_data.data, axis=-1) - - # create depth vector - z = np.linspace(0, imaging_depth, num_px_z) - - # allocate memory for beamformed image - beamformed_data = np.zeros((len(z), len(channel_data.sequence)), dtype=complex) - - # hilbert transform rf data to get envelope - channel_data.data = hilbert(channel_data.data, axis=3) - - # allocate memory for - wave_origin_x = np.empty(len(channel_data.sequence)) - - for e_id, event in enumerate(channel_data.sequence): - # todo event.event should be event.event_id or event.key - # todo make itteratable getter for events - event = channel_data.unique_events[e_id] - probe = event.receive_setup.probe - sampling_freq = event.receive_setup.sampling_frequency - # We assume one transmit wave per transmit event... hence 0 index - transmit_wave = event.transmit_setup.transmit_waves[0] - - # make time vector - time_vector = make_time_vector(num_samples=number_samples, sampling_freq=sampling_freq, - time_offset=event.receive_setup.time_offset) - - # todo: make indexing 0 min and not 1 min - wave_origin_x[e_id] = channel_data.unique_waves[transmit_wave.wave - 1].origin.position.x - - # todo: make position objects - pixel_positions = np.stack([wave_origin_x[e_id] * np.ones(len(z)), np.zeros(len(z)), z]).T - expanding_aperture = pixel_positions[:, 2] / f_number - - # time zero delays for spherical waves - origin = get_origin_array(channel_data, transmit_wave) - t0_point = get_t0(transmit_wave) - - # logging.log(logging.INFO, origin, t0_point) - - transmit_distance = np.sign(pixel_positions[:, 2] - origin[2]) * \ - np.sqrt(np.sum((pixel_positions - origin) ** 2, axis=1)) + \ - np.abs(1.2 * t0_point[0]) - # np.sqrt(np.sum((origin - t0_point) ** 2)) - - probe = channel_data.probes[probe - 1] - # todo: why are element positions saved as transforms and not positions? - transform = ShiftedTransform.deserialize(probe.transform.serialize()) - # todo: remove list from channel mapping. currently [[,]...] - - # dataset.channel_data.unique_waves[transmit_wave.wave - 1].origin.position.x - - # event.transmit_setup.channel_mapping = np.arange(1, 33) # Added by Farid - plt.plot(transmit_distance) - - for element_number in event.transmit_setup.channel_mapping: - element_number = element_number[0] # Changed by Farid - - # todo: why are element positions saved as transformations? - element_position = Position.deserialize( - probe.element[element_number - 1].transform.translation.serialize()) - element_location = Position.deserialize(transform(element_position).serialize()) - - pixel_element_lateral_distance = abs(pixel_positions[:, 0] - element_location[0]) - # logging.log(logging.INFO, pixel_element_lateral_distance) - receive_apodization = apodize(pixel_element_lateral_distance, expanding_aperture, apodization_window) - - # receive distance - receive_distance = np.sqrt(np.sum((pixel_positions - np.array(element_location)) ** 2, axis=1)) - - t0 = transmit_wave.time_offset - - # round trip delay - delay = (transmit_distance + receive_distance) / channel_data.sound_speed + t0 - - # beamformed data - chan_id = element_number - 1 - event.transmit_setup.channel_mapping[0][0] # tricky part - signal = np.squeeze(channel_data.data[:, e_id, chan_id, :]) - interp = interp1d(x=time_vector, y=signal, kind='cubic', bounds_error=False, fill_value=0) - beamformed_data[:, e_id] += np.squeeze(receive_apodization * interp(delay).T) - - # Envelope and plot - envelope_beamformed_data = np.absolute(beamformed_data) - compressed_beamformed_data = 20 * np.log10(envelope_beamformed_data / np.amax(envelope_beamformed_data) + 1e-12) - - plt.figure - x_dis = 1e3 * wave_origin_x - z_dis = 1e3 * z - plt.imshow(compressed_beamformed_data, vmin=-60, vmax=0, cmap='Greys_r', - extent=[min(x_dis), max(x_dis), max(z_dis), min(z_dis)]) - plt.xlabel('x[mm]', fontsize=12) - plt.ylabel('z[mm]', fontsize=12) - plt.title(channel_data.description) - plt.colorbar() - - 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 def focus(kgrid, input_signal, source_mask, focus_position, sound_speed): diff --git a/kwave/reconstruction/converter.py b/kwave/reconstruction/converter.py deleted file mode 100644 index 8a7f948d..00000000 --- a/kwave/reconstruction/converter.py +++ /dev/null @@ -1,175 +0,0 @@ -import time - -import numpy as np -import uff -from uff.linear_array import LinearArray - -from kwave.ktransducer import NotATransducer - - -def build_channel_data(sensor_data: np.ndarray, - kgrid, - not_transducer: NotATransducer, - sampling_frequency, - prf, - focal_depth): - number_scan_lines = sensor_data.shape[1] - input_signal = not_transducer.input_signal - transducer = not_transducer.transducer - - x0 = np.arange(0, 96) * transducer.element_width * kgrid.dy - x0 = x0 - (x0.max() / 2) - - unique_excitations = uff.excitation.Excitation( - waveform=input_signal, - sampling_frequency=sampling_frequency, - pulse_shape='Gaussian' - ) - - probe = LinearArray( - number_elements=transducer.number_elements * 4, - pitch=transducer.element_width * kgrid.dy + transducer.element_spacing * kgrid.dx, - element_width=transducer.element_width * kgrid.dy, - element_height=transducer.element_length * kgrid.dx, - transform=uff.transform.Transform(uff.transform.Translation(0, 0, 0), uff.transform.Rotation(0, 0, 0)) - ) - - ####### Move to LinearArray.update - - if probe.pitch and probe.number_elements: - - # for all lines in this block, probe.element_geometry[0] is hard coded as 0. Should be dynamic! - - # update perimeter - if probe.element_width and probe.element_height: - - if not probe.element_geometry: - probe.element_geometry = [uff.element_geometry.ElementGeometry( - uff.perimeter.Perimeter([ - uff.Position(0, 0, 0), - uff.Position(0, 0, 0), - uff.Position(0, 0, 0), - uff.Position(0, 0, 0), - ]) - )] - - if not probe.element_geometry[0].perimeter: - probe.element_geometry.perimeter = uff.perimeter.Perimeter([ - uff.Position(0, 0, 0), - uff.Position(0, 0, 0), - uff.Position(0, 0, 0), - uff.Position(0, 0, 0), - ]) - - probe.element_geometry[0].perimeter.position[0] = uff.Position(x=-probe.element_width / 2, - y=-probe.element_height / 2, z=0) - probe.element_geometry[0].perimeter.position[1] = uff.Position(x=probe.element_width / 2, - y=-probe.element_height / 2, z=0) - probe.element_geometry[0].perimeter.position[2] = uff.Position(x=probe.element_width / 2, - y=probe.element_height / 2, z=0) - probe.element_geometry[0].perimeter.position[3] = uff.Position(x=-probe.element_width / 2, - y=probe.element_height / 2, z=0) - - # compute element position in the x-axis - x0_local = np.arange(1, probe.number_elements + 1) * probe.pitch - x0_local = x0_local - x0_local.mean() - x0_local = x0_local.astype(np.float32) - - # create array of elements - probe.element = [] - for n in range(probe.number_elements): - element = uff.Element( - transform=uff.Transform( - translation=uff.Translation(x=x0_local[n], y=0, z=0), - rotation=uff.Rotation(x=0, y=0, z=0) - ), - # linear array can only have a single perimenter impulse_response - # impulse_response=[], # probe.element_impulse_response, - element_geometry=probe.element_geometry[0] - ) - - probe.element.append(element) - - ####### Move end - - unique_waves = [] - for n in range(number_scan_lines): - wave = uff.Wave( - wave_type=uff.WaveType.CONVERGING, - origin=uff.SphericalWaveOrigin( - position=uff.Position(x=x0[n], z=focal_depth) - ), - aperture=uff.Aperture( - origin=uff.Position(x=x0[n]), - fixed_size=transducer.element_length*kgrid.dx, # wrong, correct one is in below - # fixed_size=[active_aperture*el_pitch, el_height], - window='rectangular' - ) - ) - unique_waves.append(wave) - - unique_events = [] - lag = 1.370851370851371e-06 # (length(input_signal) / 2 + 1) * dt - downsample_factor = 1 - first_sample_time = 650e-9 * np.ones((number_scan_lines, 1)) # CHANGE - - for n in range(number_scan_lines): - # we select a time zero reference point: in this case the location of the first element fired - arg_min = np.argmax(np.sqrt((x0[n] - x0) ** 2 + focal_depth ** 2)) - time_zero_reference_point = x0[arg_min] - - transmit_wave = uff.TransmitWave( - wave=n, - time_zero_reference_point=uff.TimeZeroReferencePoint(x=time_zero_reference_point, y=0, z=0), - time_offset=lag, - weight=1 - ) - - # transmit and receive channel - # channel_mapping = np.arange(n, number_elements + n).tolist() - channel_mapping = np.arange(n, transducer.number_elements + n)[None, :].T # shape: [number_elements x 1] - transmit_setup = uff.TransmitSetup( - probe=1, # different from Matlab - channel_mapping=channel_mapping, - transmit_waves=[transmit_wave] # should be list! - ) - - # active_elment_positions = channel_data.probes{1}.get_elements_centre(); - - receive_setup = uff.ReceiveSetup( - probe=1, # different from Matlab - time_offset=first_sample_time[n], - channel_mapping=channel_mapping, - sampling_frequency=sampling_frequency / downsample_factor - ) - - unique_event = uff.Event( - transmit_setup=transmit_setup, - receive_setup=receive_setup - ) - - unique_events.append(unique_event) - - sequence = [] - for n in range(len(unique_events)): - sequence.append(uff.TimedEvent( - event=n, - time_offset=(n - 1) / prf - )) - - channel_data = uff.channel_data.ChannelData( - data=sensor_data, - probes=[probe], - unique_excitations=[unique_excitations], - unique_waves=unique_waves, - unique_events=unique_events, - description='US B-Mode Linear Transducer Scan, kWave.', - authors='kWave.py', - sound_speed=not_transducer.sound_speed, - system='kWave.py', - country_code='DE', - local_time=time.strftime('%Y-%m-%dTHH:MM:SS+00'), # fix - repetition_rate=prf / len(unique_events), - sequence=sequence - ) - return channel_data diff --git a/kwave/reconstruction/shifted_transform.py b/kwave/reconstruction/shifted_transform.py deleted file mode 100644 index 40f4a320..00000000 --- a/kwave/reconstruction/shifted_transform.py +++ /dev/null @@ -1,11 +0,0 @@ -from uff import Position, Transform - - -class ShiftedTransform(Transform): - def __call__(self, point): - if isinstance(point, Position): - # todo: add rotation - return Position(point.x + self.translation.x, point.y + self.translation.y, point.z + self.translation.z) - else: - raise TypeError(f"Type {type(point)} not recognized.") - pass diff --git a/kwave/reconstruction/tools.py b/kwave/reconstruction/tools.py index b05b92ca..2a06980f 100644 --- a/kwave/reconstruction/tools.py +++ b/kwave/reconstruction/tools.py @@ -1,5 +1,4 @@ import numpy as np -from uff import Position def log_compression(signal, cf, normalize=False): @@ -58,17 +57,6 @@ def apodize(distance, aperture, window): return apod -def get_t0(transmit_wave): - serialized_tx_wave = transmit_wave.time_zero_reference_point.serialize() - return np.array(Position.deserialize(serialized_tx_wave)) - - -def get_origin_array(channel_data, transmit_wave): - serialized_origin = channel_data.unique_waves[transmit_wave.wave - 1].origin.position.serialize() - return np.array( - Position.deserialize(serialized_origin)) - - def make_time_vector(num_samples, sampling_freq, time_offset): return np.linspace(start=0, num=num_samples, stop=num_samples / sampling_freq) + time_offset diff --git a/pyproject.toml b/pyproject.toml index 28531052..cef4fe23 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,7 +30,6 @@ dependencies = [ "matplotlib==3.7.2", "numpy>=1.22.2,<1.25.0", "scikit-image==0.21.0", - "uff-reader==0.0.2", ] [project.urls]