diff --git a/examples/rca/bmode/display.py b/examples/rca/bmode/display.py index b84a0fe..9960b87 100644 --- a/examples/rca/bmode/display.py +++ b/examples/rca/bmode/display.py @@ -2,8 +2,8 @@ # Display configuration file. displays = { - "OXZ": Display2D( - title=f"OXZ B-mode", + "XY": Display2D( + title=f"y = 0", layers=( Layer2D( value_range=(45, 80), @@ -13,8 +13,8 @@ ), ax_labels=("OX (m)", "OZ (m)") ), - "OYZ": Display2D( - title=f"B-mode", + "YX": Display2D( + title=f"x = 0", layers=( Layer2D( value_range=(45, 80), diff --git a/examples/rca/bmode/env.py b/examples/rca/bmode/env.py index 1f04ea4..2409929 100644 --- a/examples/rca/bmode/env.py +++ b/examples/rca/bmode/env.py @@ -1,323 +1,101 @@ -import numpy as np import arrus.medium -import arrus_rca_utils.probe_params as probe_params -from arrus.devices.probe import ProbeModel +import numpy as np from arrus.ops.imaging import * from arrus.ops.us4r import * -from arrus.ops.us4r import TxRx, Tx, Rx, Aperture, Pulse -from arrus.ops.us4r import ( - TxRxSequence -) from arrus.utils.imaging import * from arrus_rca_utils.reconstruction import ( - ReconstructHriRca, - GetFramesForRange, - Concatenate, - get_frame_ranges, - get_rx_aperture_size, - PipelineSequence, - SelectBatch, - Slice + get_reconstruction_graph +) +from arrus_rca_utils.sequence import ( + get_pw_sequence ) -from arrus_rca_utils.sequence import convert_to_system_sequence, RcaSequence from gui4us.model.envs.arrus import ( UltrasoundEnv, ArrusEnvConfiguration, Curve ) -# DEMO parameters: Vermon 128+128 elements -ARRAY_X = probe_params.ProbeArray( - n_elements=128, - pitch=0.2e-3, - start=0, - arrangement="ox" -) -ARRAY_Y = probe_params.ProbeArray( - n_elements=128, - pitch=0.2e-3, - start=128, - arrangement="oy" -) - - -# ------------------------------------------ SEQUENCE -def create_sequence( - medium: arrus.medium.Medium, - center_frequency: float, - n_periods: float, - angles: np.ndarray, - pri: float, - sample_range: Tuple[int, int], -) -> Tuple[TxRxSequence, TxRxSequence]: - """ - Returns TX/RX sequences for RCA probe. - - This function returns two sequences: - - transmit with OX elements, receive with OY elements len(angles) times - - transmit with OY elements, receive with OX elements len(angles) times - - First TX=OX, RX=OY len(angles) times, then TX=OY, RX=OX len(angles) times. - - NOTE: this function assumes the parameters as defined in the probe_params. - You have to adjust the values there in order to get the proper sequence - for your RCA probe. - """ - # Transmit with OX elements, receive with OY elements - sequence_xy = TxRxSequence( - ops=[ - TxRx( - tx=Tx( - # Transmit with all elements. - # Center == 0.0 means the center of the probe. - aperture=Aperture(center=0.0, size=ARRAY_X.n_elements), - excitation=Pulse( - center_frequency=center_frequency, - n_periods=n_periods, - inverse=False - ), - # Plane wave (focus depth = inf). - focus=np.inf, - angle=angle, - speed_of_sound=medium.speed_of_sound - ), - rx=Rx( - # Receive with all elements. - aperture=Aperture(center=0.0, size=ARRAY_Y.n_elements), - sample_range=sample_range, - ), - pri=pri - ) - for angle in angles - ], - ) - sequence_yx = TxRxSequence( - ops=[ - TxRx( - tx=Tx( - # Transmit with all elements. - # Center == 0.0 means the center of the probe. - aperture=Aperture(center=0.0, size=ARRAY_Y.n_elements), - excitation=Pulse( - center_frequency=center_frequency, - n_periods=n_periods, - inverse=False - ), - # Plane wave (focus depth = inf). - focus=np.inf, - angle=angle, - speed_of_sound=medium.speed_of_sound - ), - rx=Rx( - # Receive with all elements. - aperture=Aperture(center=0.0, size=ARRAY_X.n_elements), - sample_range=sample_range, - ), - pri=pri - ) - for angle in angles - ], - ) - return sequence_xy, sequence_yx - - -def get_system_sequence( - sequence_xy: TxRxSequence, sequence_yx: TxRxSequence, - probe_model: ProbeModel, - device_sampling_frequency: float -): - return convert_to_system_sequence( - sequences=[ - RcaSequence( - sequence=sequence_xy, - tx=ARRAY_X, rx=ARRAY_Y, - ), - RcaSequence( - sequence=sequence_yx, - tx=ARRAY_Y, rx=ARRAY_X, - ), - ], - probe_model=probe_model, - device_sampling_frequency=device_sampling_frequency - ) - - -# ------------------------------------------ RECONSTRUCTION -def reorder_rf(frames, aperture_size): - return ( - GetFramesForRange(frames=frames, aperture_size=aperture_size), - RemapToLogicalOrder() - ) - - -def to_hri( - y_grid, x_grid, z_grid, - array_tx, array_rx, - sequence, - fir_taps, name=None, - aperture_arrangement="alternate" -): - return ( - Transpose(axes=(0, 1, 3, 2)), - FirFilter(taps=fir_taps), - DigitalDownConversion(decimation_factor=4, fir_order=64), - ReconstructHriRca( - y_grid=y_grid, x_grid=x_grid, z_grid=z_grid, - probe_tx=array_tx, - probe_rx=array_rx, - sequence=sequence, - min_tang=-0.5, max_tang=0.5, - name=name, - arrangement=aperture_arrangement, - ), - ) - - -def to_bmode(dr_min, dr_max): - return ( - # Concatenate along TX axis - Concatenate(axis=0), - Mean(axis=0), - Squeeze(), - EnvelopeDetection(), - LogCompression(), - DynamicRangeAdjustment(min=dr_min, max=dr_max), - Squeeze(), - ) - -def branch(steps, name=None): - return Pipeline(steps=steps, name=name, placement="/GPU:0") - - -def get_pwi_reconstruction( - y_grid, x_grid, z_grid, - array_x: probe_params.ProbeArray, array_y: probe_params.ProbeArray, - sequence_xy: TxRxSequence, - sequence_yx: TxRxSequence, - fir_taps, - dr_min=20, dr_max=80, -): - seqs = sequence_xy, sequence_yx - range_xy, range_yx = get_frame_ranges(*seqs) - ap_size_xy, ap_size_yx = get_rx_aperture_size(*seqs) +def configure(session: arrus.Session): + us4r = session.get_device("/Us4R:0") - reconstruction = Pipeline( - steps=( - branch( - name="b", - steps=( - # TX=OY, RX=OX - *reorder_rf(frames=range_yx, aperture_size=ap_size_yx), - *to_hri( - fir_taps=fir_taps, - y_grid=y_grid, - x_grid=x_grid, - z_grid=z_grid, - array_tx=array_y, - array_rx=array_x, - sequence=sequence_yx, - ), - )), - # TX=OX, RX=OY - *reorder_rf(frames=range_xy, aperture_size=ap_size_xy), - *to_hri( - fir_taps=fir_taps, - y_grid=y_grid, - x_grid=x_grid, - z_grid=z_grid, - array_tx=array_x, array_rx=array_y, - sequence=sequence_xy, - ) - ), - placement="/GPU:0" - ) - # Expected input: - bmodes = Pipeline( - name="bmode", - steps=( - SelectBatch([0, 1]), - *to_bmode( - dr_min=dr_min, - dr_max=dr_max - ), - # (y, x, z) - branch( - steps=( - Slice(axis=0), - Squeeze(), - Transpose(), - )), - branch( - steps=( - Slice(axis=1), - Squeeze(), - Transpose(), - )) - ), - placement="/GPU:0" - ) - pipelines = [reconstruction, bmodes] - return PipelineSequence(pipelines) + # The ordinal number of the OX and OY sub-arrays. See us4r.prototxt + ox_ordinal = 0 + oy_ordinal = 1 + array_ox = us4r.get_probe_model(ordinal=ox_ordinal) + array_oy = us4r.get_probe_model(ordinal=oy_ordinal) + array_ox_id = f"Probe:{ox_ordinal}" + array_oy_id = f"Probe:{oy_ordinal}" + n_samples = 4*1024 -def configure(session: arrus.Session): - us4r = session.get_device("/Us4R:0") - medium = arrus.medium.Medium(name="ats549", speed_of_sound=1450) + medium = arrus.medium.Medium(name="ats560h", speed_of_sound=1480) angles = np.linspace(-10, 10, 32) * np.pi / 180 # [rad] center_frequency = 6e6 # [Hz] n_periods = 2 - sample_range = (0, 4*1024) + sample_range = (0, n_samples) pri = 400e-6 - # Imaging grid. - y_grid = np.arange(-8e-3, 8e-3, 0.2e-3) - x_grid = np.arange(-8e-3, 8e-3, 0.2e-3) - z_grid = np.arange(10e-3, 40e-3, 0.2e-3) - # Initial TGC curve. - tgc_sampling_points = np.linspace(np.min(z_grid), np.max(z_grid), 10) - tgc_values = np.linspace(34, 54, 10) + # Imaging grid. + y_grid = np.arange(-10e-3, 10e-3, 0.2e-3) + x_grid = np.arange(-10e-3, 10e-3, 0.2e-3) + z_grid = np.arange(5e-3, 45e-3, 0.2e-3) - # TX/RX PW sequence - sequence_xy, sequence_yx = create_sequence( + common_parameters = dict( medium=medium, angles=angles, n_periods=n_periods, center_frequency=center_frequency, sample_range=sample_range, - pri=pri + pri=pri, + ) + # Initial TGC curve. + tgc_sampling_points = np.linspace(np.min(z_grid), np.max(z_grid), 10) + tgc_values = np.linspace(34, 54, 10) + + # TX with OX elements, RX with OY elements + sequence_xy = get_pw_sequence( + array_tx_id=array_ox_id, + array_tx=array_ox, + array_rx_id=array_oy_id, + array_rx=array_oy, + name="XY", + **common_parameters + ) + # TX with OY elements, RX with OX elements + sequence_yx = get_pw_sequence( + array_tx_id=array_oy_id, + array_tx=array_oy, + array_rx_id=array_ox_id, + array_rx=array_ox, + name="YX", + **common_parameters ) - # Image reconstruction. fir_taps = scipy.signal.firwin( numtaps=64, cutoff=np.array([0.5, 1.5]) * center_frequency, pass_zero="bandpass", fs=us4r.current_sampling_frequency ) - pipeline = get_pwi_reconstruction( - array_x=ARRAY_X, - array_y=ARRAY_Y, - y_grid=y_grid, - x_grid=x_grid, - z_grid=z_grid, + + reconstruction = get_reconstruction_graph( fir_taps=fir_taps, + x_grid=x_grid, y_grid=y_grid, z_grid=z_grid, sequence_xy=sequence_xy, sequence_yx=sequence_yx, - dr_min=20, dr_max=80, + slice=True ) scheme = Scheme( - tx_rx_sequence=get_system_sequence( - sequence_xy=sequence_xy, - sequence_yx=sequence_yx, - probe_model=us4r.get_probe_model(), - device_sampling_frequency=us4r.current_sampling_frequency - ), - processing=pipeline + tx_rx_sequence=[ + sequence_xy, + sequence_yx + ], + processing=reconstruction ) return ArrusEnvConfiguration( medium=medium, scheme=scheme, tgc=Curve(points=tgc_sampling_points, values=tgc_values), - voltage=20 + voltage=40 ) diff --git a/examples/rca/bmode/us4r.prototxt b/examples/rca/bmode/us4r.prototxt index e7844a7..9bf2b18 100644 --- a/examples/rca/bmode/us4r.prototxt +++ b/examples/rca/bmode/us4r.prototxt @@ -1,35 +1,78 @@ us4r: { - probe: - { - id: { - manufacturer: "Vermon" - name: "RCA128" - } - n_elements: 256, - pitch: 0.2e-3, - tx_frequency_range: { - begin: 1e6, - end: 15e6 + probe: [ + { + id: { + manufacturer: "Vermon" + name: "RCA128-OX" + } + n_elements: 128, + pitch: 0.2e-3, + tx_frequency_range: { + begin: 1e6, + end: 15e6 + }, + voltage_range { + begin: 0, + end: 50 + } }, - voltage_range { - begin: 0, - end: 50 + { + id: { + manufacturer: "Vermon" + name: "RCA128-OY" + } + n_elements: 128, + pitch: 0.2e-3, + tx_frequency_range: { + begin: 1e6, + end: 15e6 + }, + voltage_range { + begin: 0, + end: 50 + } } - } + ] adapter_id: { manufacturer: "us4us" name: "dlp408r" } - probe_to_adapter_connection: { - channel_mapping_ranges: [ + probe_to_adapter_connection: [ + { + probe_model_id: { + manufacturer: "Vermon" + name: "RCA128-OX" + } + probe_adapter_model_id: { + manufacturer: "us4us", + name: "dlp408r" + } + channel_mapping_ranges: [ + { + begin: 0 + end: 127 + } + ] + }, { - begin: 0 - end: 255 + probe_model_id: { + manufacturer: "Vermon" + name: "RCA128-OY" + } + probe_adapter_model_id: { + manufacturer: "us4us", + name: "dlp408r" + } + channel_mapping_ranges: [ + { + begin: 128 + end: 255 + } + ] } ] - } # Default initial values. rx_settings: { @@ -42,19 +85,9 @@ us4r: { hv: { model_id { manufacturer: "us4us" - name: "hv256" + name: "us4oemhvps" } } - - channels_mask: { - channels: [] - } - - # us4oem channels mask are redundant here to minimize the risk of changing masking by mistake - us4oem_channels_mask: [ - {channels: []}, - {channels: []} - ] } diff --git a/examples/rca/holo_display/dataset.py b/examples/rca/holo_display/dataset.py deleted file mode 100644 index c327d30..0000000 --- a/examples/rca/holo_display/dataset.py +++ /dev/null @@ -1,17 +0,0 @@ -import os.path -import urllib.request -import pickle - - -def load_rca_cyst_dataset(): - filename = "rca_cyst.pickle" - if not os.path.exists(filename): - url = "https://www.dropbox.com/scl/fi/6vpt9uxne7pilhbhbtw24/vermon_rca_dataset.pickle?rlkey=pik6wku00su0b9um3p5xdb6hf&dl=1" - print("Downloading dataset...") - urllib.request.urlretrieve(url, filename) - print("... dataset downloaded.") - data = pickle.load(open(filename, "rb")) - rf = data["rf"] - metadata = data["metadata"][0] - return rf, metadata - diff --git a/examples/rca/holo_display/probe_params.py b/examples/rca/holo_display/probe_params.py deleted file mode 100644 index b0a114d..0000000 --- a/examples/rca/holo_display/probe_params.py +++ /dev/null @@ -1,16 +0,0 @@ -from arrus_rca_utils.probe_params import * - - -# DEMO parameters: Vermon 128+128 elements -APERTURE_X = ProbeArray( - n_elements=128, - pitch=0.2e-3, - start=0, - arrangement="ox" -) -APERTURE_Y = ProbeArray( - n_elements=128, - pitch=0.2e-3, - start=128, - arrangement="oy" -) diff --git a/examples/rca/holo_display/reconstruction.py b/examples/rca/holo_display/reconstruction.py deleted file mode 100644 index 00afba7..0000000 --- a/examples/rca/holo_display/reconstruction.py +++ /dev/null @@ -1,123 +0,0 @@ -from arrus_rca_utils.reconstruction import ( - ReconstructHriRca, - GetFramesForRange, - Concatenate, - get_frame_ranges, - get_rx_aperture_size, - PipelineSequence, - SelectBatch -) -from arrus.ops.us4r import ( - TxRxSequence -) -from arrus.utils.imaging import * -import probe_params - - -def reorder_rf(frames, aperture_size): - return ( - GetFramesForRange(frames=frames, aperture_size=aperture_size), - RemapToLogicalOrder() - ) - - -def to_hri( - y_grid, x_grid, z_grid, - array_tx, array_rx, - sequence, - fir_taps, name=None, - aperture_arrangement="alternate" -): - return ( - Transpose(axes=(0, 1, 3, 2)), - FirFilter(taps=fir_taps), - DigitalDownConversion(decimation_factor=8, fir_order=64), - ReconstructHriRca( - y_grid=y_grid, x_grid=x_grid, z_grid=z_grid, - probe_tx=array_tx, - probe_rx=array_rx, - sequence=sequence, - min_tang=-0.5, max_tang=0.5, - name=name, - arrangement=aperture_arrangement, - ), - ) - - -def to_bmode(dr_min, dr_max): - return ( - # Concatenate along TX axis - Concatenate(axis=0), - Mean(axis=0), - Squeeze(), - EnvelopeDetection(), - LogCompression(), - DynamicRangeAdjustment(min=dr_min, max=dr_max), - Squeeze(), - ) - - -def branch(steps, name=None): - return Pipeline(steps=steps, name=name, placement="/GPU:0") - - -def get_pwi_reconstruction( - y_grid, x_grid, z_grid, - array_x: probe_params.ProbeArray, array_y: probe_params.ProbeArray, - sequence_xy: TxRxSequence, - sequence_yx: TxRxSequence, - fir_taps, - dr_min=20, dr_max=80, -): - seqs = sequence_xy, sequence_yx - range_xy, range_yx = get_frame_ranges(*seqs) - ap_size_xy, ap_size_yx = get_rx_aperture_size(*seqs) - - reconstruction = Pipeline( - steps=( - branch( - name="b", - steps=( - # TX=OY, RX=OX - *reorder_rf(frames=range_yx, aperture_size=ap_size_yx), - *to_hri( - fir_taps=fir_taps, - y_grid=y_grid, - x_grid=x_grid, - z_grid=z_grid, - array_tx=array_y, - array_rx=array_x, - sequence=sequence_yx, - ), - )), - # TX=OX, RX=OY - *reorder_rf(frames=range_xy, aperture_size=ap_size_xy), - *to_hri( - fir_taps=fir_taps, - y_grid=y_grid, - x_grid=x_grid, - z_grid=z_grid, - array_tx=array_x, array_rx=array_y, - sequence=sequence_xy, - ) - ), - placement="/GPU:0" - ) - # Expected input: - bmodes = Pipeline( - name="bmode", - steps=( - SelectBatch([0, 1]), - *to_bmode( - dr_min=dr_min, - dr_max=dr_max - ), - - ), - placement="/GPU:0" - ) - - pipelines = [reconstruction, bmodes] - - pipeline = PipelineSequence(pipelines) - return pipeline diff --git a/examples/rca/holo_display/requirements.txt b/examples/rca/holo_display/requirements.txt deleted file mode 100644 index db85083..0000000 --- a/examples/rca/holo_display/requirements.txt +++ /dev/null @@ -1,3 +0,0 @@ -git+ssh://git@github.com/us4useu/arrus-toolkit.git@master#subdirectory=examples/rca/utils -vtk==9.2.2 -vtk-lookingglass==1.0.0 diff --git a/examples/rca/holo_display/run_3d_us4r.py b/examples/rca/holo_display/run_3d_us4r.py deleted file mode 100644 index f806441..0000000 --- a/examples/rca/holo_display/run_3d_us4r.py +++ /dev/null @@ -1,90 +0,0 @@ -import arrus -import arrus.medium -import arrus.ops.tgc -import arrus.session -import arrus.utils.imaging -import arrus.utils.us4r -import scipy.signal -from arrus.ops.us4r import Scheme -from arrus.utils.imaging import * - -import probe_params -from reconstruction import get_pwi_reconstruction -from sequence import ( - create_sequence, - get_system_sequence -) -from visualizer import VTKVisualizer - - -def main(): - with arrus.Session("us4r.prototxt") as sess: - us4r = sess.get_device("/Us4R:0") - us4r.set_hv_voltage(30) - MEDIUM = arrus.medium.Medium(name="tissue", speed_of_sound=1540) - angles = np.linspace(-10, 10, 32) * np.pi / 180 # [rad] - center_frequency = 6e6 # [Hz] - n_periods = 2 - sample_range = (0, 5 * 1024) - pri = 400e-6 - - # TX/RX PW sequence - sequence_xy, sequence_yx = create_sequence( - medium=MEDIUM, - angles=angles, - n_periods=n_periods, - center_frequency=center_frequency, - sample_range=sample_range, - pri=pri - ) - - - # Image reconstruction. - fir_taps = scipy.signal.firwin( - numtaps=64, cutoff=np.array([0.5, 1.5]) * center_frequency, - pass_zero="bandpass", fs=us4r.current_sampling_frequency - ) - - pipeline = get_pwi_reconstruction( - array_x=probe_params.APERTURE_X, - array_y=probe_params.APERTURE_Y, - y_grid=np.arange(-6e-3, 6e-3, 0.3e-3), - x_grid=np.arange(-6e-3, 6e-3, 0.3e-3), - z_grid=np.arange(20e-3, 43e-3, 0.3e-3), - fir_taps=fir_taps, - sequence_xy=sequence_xy, - sequence_yx=sequence_yx, - dr_min=5, dr_max=120, - ) - # TODO(pjarosik) avoid using the below wrapper - pipeline_wrapper = Pipeline( - steps=( - Lambda( - lambda data: pipeline.process(data)[0], - prepare_func=lambda metadata: pipeline.prepare(metadata)[0] - ), - ), - placement="/GPU:0" - ) - scheme = Scheme( - tx_rx_sequence=get_system_sequence( - sequence_xy=sequence_xy, - sequence_yx=sequence_yx, - probe_model=us4r.get_probe_model(), - device_sampling_frequency=us4r.current_sampling_frequency - ), - processing=pipeline_wrapper - ) - buffer, output_metadata = sess.upload(scheme) - visualizer = VTKVisualizer(output_metadata.input_shape, use_lgf=False) - print("Press CTRL+C multiple times to stop the example") - - us4r.set_tgc([50] * 26) - sess.start_scheme() - while True: - volume = buffer.get()[0] - visualizer.update(volume) - - -if __name__ == "__main__": - main() diff --git a/examples/rca/holo_display/sequence.py b/examples/rca/holo_display/sequence.py deleted file mode 100644 index 7690f35..0000000 --- a/examples/rca/holo_display/sequence.py +++ /dev/null @@ -1,116 +0,0 @@ -import numpy as np -import arrus.medium -from arrus.ops.us4r import TxRxSequence, TxRx, Tx, Rx, Aperture, Pulse -from arrus.devices.probe import ProbeModel -from typing import Tuple - -from arrus_rca_utils.sequence import convert_to_system_sequence, RcaSequence -import probe_params - - -def create_sequence( - medium: arrus.medium.Medium, - center_frequency: float, - n_periods: float, - angles: np.ndarray, - pri: float, - sample_range: Tuple[int, int], -) -> Tuple[TxRxSequence, TxRxSequence]: - """ - Returns TX/RX sequences for RCA probe. - - This function returns two sequences: - - transmit with OX elements, receive with OY elements len(angles) times - - transmit with OY elements, receive with OX elements len(angles) times - - First TX=OX, RX=OY len(angles) times, then TX=OY, RX=OX len(angles) times. - - NOTE: this function assumes the parameters as defined in the probe_params. - You have to adjust the values there in order to get the proper sequence - for your RCA probe. - """ - # Transmit with OX elements, receive with OY elements - sequence_xy = TxRxSequence( - ops=[ - TxRx( - tx=Tx( - # Transmit with all elements. - # Center == 0.0 means the center of the probe. - aperture=Aperture(center=0.0, - size=probe_params.APERTURE_X.n_elements), - excitation=Pulse( - center_frequency=center_frequency, - n_periods=n_periods, - inverse=False - ), - # Plane wave (focus depth = inf). - focus=np.inf, - angle=angle, - speed_of_sound=medium.speed_of_sound - ), - rx=Rx( - # Receive with all elements. - aperture=Aperture(center=0.0, - size=probe_params.APERTURE_Y.n_elements), - sample_range=sample_range, - ), - pri=pri - ) - for angle in angles - ], - ) - - # Transmit with OY elements, receive with OX elements. - sequence_yx = TxRxSequence( - ops=[ - TxRx( - tx=Tx( - # Transmit with all elements. - # Center == 0.0 means the center of the probe. - aperture=Aperture( - center=0.0, - size=probe_params.APERTURE_Y.n_elements), - excitation=Pulse( - center_frequency=center_frequency, - n_periods=n_periods, - inverse=False - ), - # Plane wave (focus depth = inf). - focus=np.inf, - angle=angle, - speed_of_sound=medium.speed_of_sound - ), - rx=Rx( - # Receive with all elements. - aperture=Aperture( - center=0.0, - size=probe_params.APERTURE_X.n_elements), - sample_range=sample_range, - ), - pri=pri - ) - for angle in angles - ], - ) - return sequence_xy, sequence_yx - - -def get_system_sequence( - sequence_xy: TxRxSequence, sequence_yx: TxRxSequence, - probe_model: ProbeModel, - device_sampling_frequency: float -): - return convert_to_system_sequence( - sequences=[ - RcaSequence( - sequence=sequence_xy, - tx=probe_params.APERTURE_X, rx=probe_params.APERTURE_Y, - ), - RcaSequence( - sequence=sequence_yx, - tx=probe_params.APERTURE_Y, rx=probe_params.APERTURE_X, - ), - ], - probe_model=probe_model, - device_sampling_frequency=device_sampling_frequency - ) diff --git a/examples/rca/holo_display/us4r.prototxt b/examples/rca/holo_display/us4r.prototxt deleted file mode 100644 index e7844a7..0000000 --- a/examples/rca/holo_display/us4r.prototxt +++ /dev/null @@ -1,60 +0,0 @@ -us4r: { - probe: - { - id: { - manufacturer: "Vermon" - name: "RCA128" - } - n_elements: 256, - pitch: 0.2e-3, - tx_frequency_range: { - begin: 1e6, - end: 15e6 - }, - voltage_range { - begin: 0, - end: 50 - } - } - - adapter_id: { - manufacturer: "us4us" - name: "dlp408r" - } - - probe_to_adapter_connection: { - channel_mapping_ranges: [ - { - begin: 0 - end: 255 - } - ] - } - - # Default initial values. - rx_settings: { - lna_gain: 24 - pga_gain: 30 - lpf_cutoff: 15000000 - active_termination: 200 - } - - hv: { - model_id { - manufacturer: "us4us" - name: "hv256" - } - } - - channels_mask: { - channels: [] - } - - # us4oem channels mask are redundant here to minimize the risk of changing masking by mistake - us4oem_channels_mask: [ - {channels: []}, - {channels: []} - ] -} - - diff --git a/examples/rca/holo_display/visualizer.py b/examples/rca/holo_display/visualizer.py deleted file mode 100644 index de211df..0000000 --- a/examples/rca/holo_display/visualizer.py +++ /dev/null @@ -1,120 +0,0 @@ -import vtkmodules.vtkInteractionStyle -from vtkmodules.vtkCommonColor import vtkNamedColors -from vtkmodules.vtkCommonDataModel import vtkPiecewiseFunction -from vtkmodules.vtkIOLegacy import vtkStructuredPointsReader -from vtkmodules.vtkRenderingCore import ( - vtkColorTransferFunction, - vtkRenderWindow, - vtkRenderWindowInteractor, - vtkRenderer, - vtkVolume, - vtkVolumeProperty -) -from vtkmodules.vtkRenderingCore import ( - vtkActor, - vtkActor2D, - vtkPolyDataMapper, - vtkImageMapper, - vtkRenderWindow, - vtkRenderWindowInteractor, - vtkRenderer, - vtkCamera, - vtkWindowToImageFilter -) -from vtkmodules.vtkIOImage import vtkPNGWriter -from vtkmodules.vtkRenderingVolume import vtkFixedPointVolumeRayCastMapper -from vtkmodules.vtkCommonDataModel import vtkImageData -import vtk.util.numpy_support -import pickle -import numpy as np -import sys -import matplotlib.pyplot as plt -import time -import scipy.ndimage -import vtk -from vtkmodules.vtkRenderingVolumeOpenGL2 import vtkOpenGLRayCastImageDisplayHelper -from vtk import vtkRenderingLookingGlass - - -class VTKVisualizer: - - def __init__(self, dimensions, use_lgf=True): - colors = vtkNamedColors() - # Create a looking glass render window - ren1 = vtk.vtkRenderer() - iren = vtkRenderWindowInteractor() - - if use_lgf: - self.renWin = vtkRenderingLookingGlass.vtkLookingGlassInterface.CreateLookingGlassRenderWindow() - if self.renWin.GetDeviceType() == "standard": - # This looks better on large settings - self.renWin.SetDeviceType("large") - num_tiles = np.prod(self.renWin.GetInterface().GetQuiltTiles()) - iren.SetRenderWindow(self.renWin) - iren.SetDesiredUpdateRate(iren.GetDesiredUpdateRate()*num_tiles) - else: - self.renWin = vtkRenderWindow() - iren.SetRenderWindow(self.renWin) - - self.renWin.AddRenderer(ren1) - - self.data = vtkImageData() - dimensions = tuple(reversed(dimensions)) - self.data.SetDimensions(*dimensions) - self.data.SetSpacing([1, 1, 1]) - self.data.SetOrigin([0, 0, 0]) - - opacityTransferFunction = vtkPiecewiseFunction() - opacityTransferFunction.AddPoint(0, 0.0) - opacityTransferFunction.AddPoint(40, 0.0) - opacityTransferFunction.AddPoint(120, 1.0) - - # Create transfer mapping scalar value to color. - colorTransferFunction = vtkColorTransferFunction() - colorTransferFunction.AddRGBPoint(0.0, 1.0, 1.0, 1.0) - - # The property describes how the data will look. - volumeProperty = vtkVolumeProperty() - volumeProperty.SetColor(colorTransferFunction) - volumeProperty.SetScalarOpacity(opacityTransferFunction) - volumeProperty.ShadeOn() - volumeProperty.SetInterpolationTypeToLinear() - - # The mapper / ray cast function know how to render the data. - volumeMapper = vtk.vtkGPUVolumeRayCastMapper() - # or vtk.vtkOpenGLGPUVolumeRayCastMapper() - # or vtkFixedPointVolumeRayCastMapper() - volumeMapper.SetInputData(self.data) - - # The volume holds the mapper and the property and - # can be used to position/orient the volume. - volume = vtkVolume() - volume.SetMapper(volumeMapper) - volume.SetProperty(volumeProperty) - - ren1.AddVolume(volume) - ren1.SetBackground(colors.GetColor3d('Wheat')) - ren1.GetActiveCamera().Azimuth(-10) - ren1.GetActiveCamera().Elevation(-10) - ren1.GetActiveCamera().Roll(-90) - ren1.ResetCameraClippingRange() - ren1.ResetCamera() - iren_style = vtk.vtkInteractorStyleTrackballCamera() - iren.SetInteractorStyle(iren_style) - - def update(self, volume): - frames = volume - gain = np.linspace(1, 3, frames.shape[-1]) - gain = gain.reshape(1, 1, -1) - # frames = frames*gain - # frames[frames == -np.inf] = np.max(frames) - frames = -frames - frames -= np.min(frames) - frames = -frames - frames = frames - np.min(frames) - frames = frames[:, :, :] - frames = scipy.ndimage.median_filter(frames, size=5) - volume_data = vtk.util.numpy_support.numpy_to_vtk(frames.ravel(), deep=True) - self.data.GetPointData().SetScalars(volume_data) - self.renWin.Render() - diff --git a/examples/rca/holo_display_dataset/dataset.py b/examples/rca/holo_display_dataset/dataset.py deleted file mode 100644 index c327d30..0000000 --- a/examples/rca/holo_display_dataset/dataset.py +++ /dev/null @@ -1,17 +0,0 @@ -import os.path -import urllib.request -import pickle - - -def load_rca_cyst_dataset(): - filename = "rca_cyst.pickle" - if not os.path.exists(filename): - url = "https://www.dropbox.com/scl/fi/6vpt9uxne7pilhbhbtw24/vermon_rca_dataset.pickle?rlkey=pik6wku00su0b9um3p5xdb6hf&dl=1" - print("Downloading dataset...") - urllib.request.urlretrieve(url, filename) - print("... dataset downloaded.") - data = pickle.load(open(filename, "rb")) - rf = data["rf"] - metadata = data["metadata"][0] - return rf, metadata - diff --git a/examples/rca/holo_display_dataset/probe_params.py b/examples/rca/holo_display_dataset/probe_params.py deleted file mode 100644 index b0a114d..0000000 --- a/examples/rca/holo_display_dataset/probe_params.py +++ /dev/null @@ -1,16 +0,0 @@ -from arrus_rca_utils.probe_params import * - - -# DEMO parameters: Vermon 128+128 elements -APERTURE_X = ProbeArray( - n_elements=128, - pitch=0.2e-3, - start=0, - arrangement="ox" -) -APERTURE_Y = ProbeArray( - n_elements=128, - pitch=0.2e-3, - start=128, - arrangement="oy" -) diff --git a/examples/rca/holo_display_dataset/reconstruction.py b/examples/rca/holo_display_dataset/reconstruction.py deleted file mode 100644 index 00afba7..0000000 --- a/examples/rca/holo_display_dataset/reconstruction.py +++ /dev/null @@ -1,123 +0,0 @@ -from arrus_rca_utils.reconstruction import ( - ReconstructHriRca, - GetFramesForRange, - Concatenate, - get_frame_ranges, - get_rx_aperture_size, - PipelineSequence, - SelectBatch -) -from arrus.ops.us4r import ( - TxRxSequence -) -from arrus.utils.imaging import * -import probe_params - - -def reorder_rf(frames, aperture_size): - return ( - GetFramesForRange(frames=frames, aperture_size=aperture_size), - RemapToLogicalOrder() - ) - - -def to_hri( - y_grid, x_grid, z_grid, - array_tx, array_rx, - sequence, - fir_taps, name=None, - aperture_arrangement="alternate" -): - return ( - Transpose(axes=(0, 1, 3, 2)), - FirFilter(taps=fir_taps), - DigitalDownConversion(decimation_factor=8, fir_order=64), - ReconstructHriRca( - y_grid=y_grid, x_grid=x_grid, z_grid=z_grid, - probe_tx=array_tx, - probe_rx=array_rx, - sequence=sequence, - min_tang=-0.5, max_tang=0.5, - name=name, - arrangement=aperture_arrangement, - ), - ) - - -def to_bmode(dr_min, dr_max): - return ( - # Concatenate along TX axis - Concatenate(axis=0), - Mean(axis=0), - Squeeze(), - EnvelopeDetection(), - LogCompression(), - DynamicRangeAdjustment(min=dr_min, max=dr_max), - Squeeze(), - ) - - -def branch(steps, name=None): - return Pipeline(steps=steps, name=name, placement="/GPU:0") - - -def get_pwi_reconstruction( - y_grid, x_grid, z_grid, - array_x: probe_params.ProbeArray, array_y: probe_params.ProbeArray, - sequence_xy: TxRxSequence, - sequence_yx: TxRxSequence, - fir_taps, - dr_min=20, dr_max=80, -): - seqs = sequence_xy, sequence_yx - range_xy, range_yx = get_frame_ranges(*seqs) - ap_size_xy, ap_size_yx = get_rx_aperture_size(*seqs) - - reconstruction = Pipeline( - steps=( - branch( - name="b", - steps=( - # TX=OY, RX=OX - *reorder_rf(frames=range_yx, aperture_size=ap_size_yx), - *to_hri( - fir_taps=fir_taps, - y_grid=y_grid, - x_grid=x_grid, - z_grid=z_grid, - array_tx=array_y, - array_rx=array_x, - sequence=sequence_yx, - ), - )), - # TX=OX, RX=OY - *reorder_rf(frames=range_xy, aperture_size=ap_size_xy), - *to_hri( - fir_taps=fir_taps, - y_grid=y_grid, - x_grid=x_grid, - z_grid=z_grid, - array_tx=array_x, array_rx=array_y, - sequence=sequence_xy, - ) - ), - placement="/GPU:0" - ) - # Expected input: - bmodes = Pipeline( - name="bmode", - steps=( - SelectBatch([0, 1]), - *to_bmode( - dr_min=dr_min, - dr_max=dr_max - ), - - ), - placement="/GPU:0" - ) - - pipelines = [reconstruction, bmodes] - - pipeline = PipelineSequence(pipelines) - return pipeline diff --git a/examples/rca/holo_display_dataset/requirements.txt b/examples/rca/holo_display_dataset/requirements.txt deleted file mode 100644 index 8147aa3..0000000 --- a/examples/rca/holo_display_dataset/requirements.txt +++ /dev/null @@ -1,3 +0,0 @@ -git+ssh://git@github.com/us4useu/arrus-toolkit.git@ref-ARRUS-224#subdirectory=examples/rca/utils -vtk==9.2.2 -vtk-lookingglass==1.0.0 diff --git a/examples/rca/holo_display_dataset/run_3d_dataset.py b/examples/rca/holo_display_dataset/run_3d_dataset.py deleted file mode 100644 index f382f6d..0000000 --- a/examples/rca/holo_display_dataset/run_3d_dataset.py +++ /dev/null @@ -1,67 +0,0 @@ -import cupy as cp -import scipy.signal -import numpy as np - -import arrus.medium - -import probe_params -from dataset import load_rca_cyst_dataset -from sequence import create_sequence -from reconstruction import get_pwi_reconstruction - -from visualizer import VTKVisualizer - - -def main(): - rf, metadata = load_rca_cyst_dataset() - MEDIUM = arrus.medium.Medium(name="tissue", speed_of_sound=1540) - angles = np.linspace(-10, 10, 64) * np.pi / 180 # [rad] - center_frequency = 6e6 # [Hz] - n_periods = 2 - sample_range = (0, 5 * 1024) - pri = 400e-6 - - # TX/RX PW sequence - sequence_xy, sequence_yx = create_sequence( - medium=MEDIUM, - angles=angles, - n_periods=n_periods, - center_frequency=center_frequency, - sample_range=sample_range, - pri=pri - ) - - # Image reconstruction. - fir_taps = scipy.signal.firwin( - numtaps=64, cutoff=np.array([0.5, 1.5]) * center_frequency, - pass_zero="bandpass", fs=metadata.context.device.sampling_frequency - ) - - pipeline = get_pwi_reconstruction( - array_x=probe_params.APERTURE_X, - array_y=probe_params.APERTURE_Y, - y_grid=np.arange(-6e-3, 6e-3, 0.4e-3), - x_grid=np.arange(-6e-3, 6e-3, 0.4e-3), - z_grid=np.arange(25e-3, 43e-3, 0.4e-3), - fir_taps=fir_taps, - sequence_xy=sequence_xy, - sequence_yx=sequence_yx, - dr_min=-5, dr_max=120, - ) - output_metadata = pipeline.prepare(metadata) - visualizer = VTKVisualizer(output_metadata[0].input_shape, use_lgf=False) - - print("Press CTRL+C multiple times to stop the example") - - n_frames = len(rf) - - i = 0 - while True: - output = pipeline.process(cp.asarray(rf[i % n_frames])) - vol = output[0].get() - visualizer.update(vol) - i += 1 - - -if __name__ == "__main__": - main() diff --git a/examples/rca/holo_display_dataset/sequence.py b/examples/rca/holo_display_dataset/sequence.py deleted file mode 100644 index af90102..0000000 --- a/examples/rca/holo_display_dataset/sequence.py +++ /dev/null @@ -1,95 +0,0 @@ -import numpy as np -import arrus.medium -from arrus.ops.us4r import TxRxSequence, TxRx, Tx, Rx, Aperture, Pulse -from typing import Tuple - -import probe_params - - -def create_sequence( - medium: arrus.medium.Medium, - center_frequency: float, - n_periods: float, - angles: np.ndarray, - pri: float, - sample_range: Tuple[int, int], -) -> Tuple[TxRxSequence, TxRxSequence]: - """ - Returns TX/RX sequences for RCA probe. - - This function returns two sequences: - - transmit with OX elements, receive with OY elements len(angles) times - - transmit with OY elements, receive with OX elements len(angles) times - - First TX=OX, RX=OY len(angles) times, then TX=OY, RX=OX len(angles) times. - - NOTE: this function assumes the parameters as defined in the probe_params. - You have to adjust the values there in order to get the proper sequence - for your RCA probe. - """ - # Transmit with OX elements, receive with OY elements - sequence_xy = TxRxSequence( - ops=[ - TxRx( - tx=Tx( - # Transmit with all elements. - # Center == 0.0 means the center of the probe. - aperture=Aperture(center=0.0, - size=probe_params.APERTURE_X.n_elements), - excitation=Pulse( - center_frequency=center_frequency, - n_periods=n_periods, - inverse=False - ), - # Plane wave (focus depth = inf). - focus=np.inf, - angle=angle, - speed_of_sound=medium.speed_of_sound - ), - rx=Rx( - # Receive with all elements. - aperture=Aperture(center=0.0, - size=probe_params.APERTURE_Y.n_elements), - sample_range=sample_range, - ), - pri=pri - ) - for angle in angles - ], - tgc_curve=[] # Will be set later. - ) - - # Transmit with OY elements, receive with OX elements. - sequence_yx = TxRxSequence( - ops=[ - TxRx( - tx=Tx( - # Transmit with all elements. - # Center == 0.0 means the center of the probe. - aperture=Aperture( - center=0.0, - size=probe_params.APERTURE_Y.n_elements), - excitation=Pulse( - center_frequency=center_frequency, - n_periods=n_periods, - inverse=False - ), - # Plane wave (focus depth = inf). - focus=np.inf, - angle=angle, - speed_of_sound=medium.speed_of_sound - ), - rx=Rx( - # Receive with all elements. - aperture=Aperture( - center=0.0, - size=probe_params.APERTURE_X.n_elements), - sample_range=sample_range, - ), - pri=pri - ) - for angle in angles - ], - tgc_curve=[] # Will be set later. - ) - return sequence_xy, sequence_yx diff --git a/examples/rca/holo_display_dataset/visualizer.py b/examples/rca/holo_display_dataset/visualizer.py deleted file mode 100644 index de211df..0000000 --- a/examples/rca/holo_display_dataset/visualizer.py +++ /dev/null @@ -1,120 +0,0 @@ -import vtkmodules.vtkInteractionStyle -from vtkmodules.vtkCommonColor import vtkNamedColors -from vtkmodules.vtkCommonDataModel import vtkPiecewiseFunction -from vtkmodules.vtkIOLegacy import vtkStructuredPointsReader -from vtkmodules.vtkRenderingCore import ( - vtkColorTransferFunction, - vtkRenderWindow, - vtkRenderWindowInteractor, - vtkRenderer, - vtkVolume, - vtkVolumeProperty -) -from vtkmodules.vtkRenderingCore import ( - vtkActor, - vtkActor2D, - vtkPolyDataMapper, - vtkImageMapper, - vtkRenderWindow, - vtkRenderWindowInteractor, - vtkRenderer, - vtkCamera, - vtkWindowToImageFilter -) -from vtkmodules.vtkIOImage import vtkPNGWriter -from vtkmodules.vtkRenderingVolume import vtkFixedPointVolumeRayCastMapper -from vtkmodules.vtkCommonDataModel import vtkImageData -import vtk.util.numpy_support -import pickle -import numpy as np -import sys -import matplotlib.pyplot as plt -import time -import scipy.ndimage -import vtk -from vtkmodules.vtkRenderingVolumeOpenGL2 import vtkOpenGLRayCastImageDisplayHelper -from vtk import vtkRenderingLookingGlass - - -class VTKVisualizer: - - def __init__(self, dimensions, use_lgf=True): - colors = vtkNamedColors() - # Create a looking glass render window - ren1 = vtk.vtkRenderer() - iren = vtkRenderWindowInteractor() - - if use_lgf: - self.renWin = vtkRenderingLookingGlass.vtkLookingGlassInterface.CreateLookingGlassRenderWindow() - if self.renWin.GetDeviceType() == "standard": - # This looks better on large settings - self.renWin.SetDeviceType("large") - num_tiles = np.prod(self.renWin.GetInterface().GetQuiltTiles()) - iren.SetRenderWindow(self.renWin) - iren.SetDesiredUpdateRate(iren.GetDesiredUpdateRate()*num_tiles) - else: - self.renWin = vtkRenderWindow() - iren.SetRenderWindow(self.renWin) - - self.renWin.AddRenderer(ren1) - - self.data = vtkImageData() - dimensions = tuple(reversed(dimensions)) - self.data.SetDimensions(*dimensions) - self.data.SetSpacing([1, 1, 1]) - self.data.SetOrigin([0, 0, 0]) - - opacityTransferFunction = vtkPiecewiseFunction() - opacityTransferFunction.AddPoint(0, 0.0) - opacityTransferFunction.AddPoint(40, 0.0) - opacityTransferFunction.AddPoint(120, 1.0) - - # Create transfer mapping scalar value to color. - colorTransferFunction = vtkColorTransferFunction() - colorTransferFunction.AddRGBPoint(0.0, 1.0, 1.0, 1.0) - - # The property describes how the data will look. - volumeProperty = vtkVolumeProperty() - volumeProperty.SetColor(colorTransferFunction) - volumeProperty.SetScalarOpacity(opacityTransferFunction) - volumeProperty.ShadeOn() - volumeProperty.SetInterpolationTypeToLinear() - - # The mapper / ray cast function know how to render the data. - volumeMapper = vtk.vtkGPUVolumeRayCastMapper() - # or vtk.vtkOpenGLGPUVolumeRayCastMapper() - # or vtkFixedPointVolumeRayCastMapper() - volumeMapper.SetInputData(self.data) - - # The volume holds the mapper and the property and - # can be used to position/orient the volume. - volume = vtkVolume() - volume.SetMapper(volumeMapper) - volume.SetProperty(volumeProperty) - - ren1.AddVolume(volume) - ren1.SetBackground(colors.GetColor3d('Wheat')) - ren1.GetActiveCamera().Azimuth(-10) - ren1.GetActiveCamera().Elevation(-10) - ren1.GetActiveCamera().Roll(-90) - ren1.ResetCameraClippingRange() - ren1.ResetCamera() - iren_style = vtk.vtkInteractorStyleTrackballCamera() - iren.SetInteractorStyle(iren_style) - - def update(self, volume): - frames = volume - gain = np.linspace(1, 3, frames.shape[-1]) - gain = gain.reshape(1, 1, -1) - # frames = frames*gain - # frames[frames == -np.inf] = np.max(frames) - frames = -frames - frames -= np.min(frames) - frames = -frames - frames = frames - np.min(frames) - frames = frames[:, :, :] - frames = scipy.ndimage.median_filter(frames, size=5) - volume_data = vtk.util.numpy_support.numpy_to_vtk(frames.ravel(), deep=True) - self.data.GetPointData().SetScalars(volume_data) - self.renWin.Render() - diff --git a/examples/rca/offline/6464.prototxt b/examples/rca/offline/6464.prototxt deleted file mode 100644 index a324a43..0000000 --- a/examples/rca/offline/6464.prototxt +++ /dev/null @@ -1,63 +0,0 @@ -us4r: { - probe: - { - id: { - manufacturer: "vermon" - name: "RCA3/64+64-1792" - } - n_elements: 128, - pitch: 0.28e-3, - tx_frequency_range: { - begin: 1e6, - end: 15e6 - }, - voltage_range { - begin: 0, - end: 50 - } - } - - adapter_id: { - manufacturer: "us4us" - name: "dlp408r" - } - - probe_to_adapter_connection: { - channel_mapping: [4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 34, 32, 30, 28, 26, 24, - 54, 52, 50, 48, 46, 44, 56, 58, 60, 62, 64, 66, 68, 70, 72, 74, - 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 35, 33, 31, 29, 27, 25, - 55, 53, 51, 49, 47, 45, 57, 59, 61, 63, 65, 67, 69, 71, 73, 75, - 84, 86, 88, 90, 92, 94, 96, 98, 100, 102, 114, 112, 110, 108, 106, 104, - 134, 132, 130, 128, 126, 124, 136, 138, 140, 142, 144, 146, 148, 150, 152, 154, - 85, 87, 89, 91, 93, 95, 97, 99, 101, 103, 115, 113, 111, 109, 107, 105, - 135, 133, 131, 129, 127, 125, 137, 139, 141, 143, 145, 147, 149, 151, 153, 155] - } - - # Default initial values. - rx_settings: { - lna_gain: 24 - pga_gain: 30 - lpf_cutoff: 15000000 - active_termination: 200 - # dtgc_attenuation: 18 - } - - hv: { - model_id { - manufacturer: "us4us" - name: "us4oemhvps" - } - } - - channels_mask: { - channels: [] - } - - # us4oem channels mask are redundant here to minimize the risk of changing masking by mistake - us4oem_channels_mask: [ - {channels: []}, - {channels: []} - ] -} - - diff --git a/examples/rca/offline/README.md b/examples/rca/offline/README.md deleted file mode 100644 index b80e4ba..0000000 --- a/examples/rca/offline/README.md +++ /dev/null @@ -1,82 +0,0 @@ -# RF Data Acquisition Scripts for us4R-lite with RCA Probe - -This repository contains a set of scripts for collecting large RF channel data sets from the **us4R-lite ultrasound system** with a **Vermon RCA 64+64 probe**. - -The scripts are designed to acquire long RF frame sequences **without on-the-fly processing**, so that the **maximum frame rate is not limited by computation**. -The frame rate is constrained only by: -- ultrasound wave propagation time, -- data transfer time, -- host-to-host memory copy time. - ---- - -## Requirements - -- **ARRUS 0.10.x** (recommended: `0.10.6`) -- **us4R-lite+** with **Vermon RCA 64+64 probe** - ---- - -## Installation - -1. Clone this repository or download this folder. -2. Install additional Python tools (required for the `display.ipynb` notebook): - ```bash - pip install git+https://github.com/pjarosik/pjtools - ``` - ---- - -## Usage - -### 1. Acquiring RF Data - -To perform a measurement, use **`acquire_offline.py`**: - -1. Adjust acquisition parameters defined in the `__main__` section of the script (e.g., pulse length, voltage, etc.). In particular: - - `batch_size`, i.e. the number of frames that will be acquired into the ultrasound system's internal DDR memory, before sending it to the host PC RAM memory. Increasing this value increases the amount of data transferred from the system to the computer, thereby reducing the impact of transfer overhead on the frame rate. - - `nframes`: total number of frames to acquire. NOTE: this must be a multiple of the `batch_size`. - - `pri`: pulse repetition interval, i.e. the time interval between two physical TX/RX events. NOTE: for the us4R-lite+ and the Vermon RCA 64+64 probe, each full RX aperture acquisition requires two TX/RXs (due to the quite unusual pin-mapping between the probe and the ultrasound system). - - In order to determine the maximum frame rate, try reducing the `pri` until you encounter buffer overflow errors (see the description below). -2. Run the script: - ```bash - python acquire_offline.py dataset_name - ``` - - -NOTE: if you get buffer overflow errors, e.g.: -``` -WARNING: Detected RX data overflow. -``` - -it is likely that the data transfer from the ultrasound system to the host PC cannot keep up with the ultrasound frame acquisition rate (according to the given PRI). To mitigate this issue, consider increasing the PRI or enlarging the batch size. - -If successful, the script will generate `dataset_name_raw_dataset.pkl`: - - the contains *physically* raw RF data (i.e. exactly in the format as produced by the ultrasound system), - - metadata, including frame rate (`time_intervals` – the time between consecutive batches). NOTE: this is the time between consecutive batch acquistions, to get the actual frame rate, divide it by the `batch_size` (see the example in the `acquire_ofline.py` script), - - see the `__main__` section in `acquire_offline.py` for the full list of stored information. - ---- - -### 2. Reconstructing High-Resolution Images - -To process the acquired dataset into a **High-Resolution Image (HRI)**, use **`reconstruct.py`**: - -1. Adjust reconstruction parameters in the `__main__` section of the script (e.g., grid resolution). - > Note: here you are not limited by real-time constraints – feel free to use optimal parameters. -2. Run the script: - ```bash - python reconstruct.py dataset_name - ``` - -If successful, the script will generate `dataset_name_beamformed_dataset.h5`, which contains: - - high-resolution images (i.e. I/Q beamformed data), - - RF data in the logical order, i.e. with shape: (frame number, transmit number, sample count, receive channels), - ---- - -### 3. Visualizing RF and HRI Data - -You can visualize both raw RF data and the reconstructed HRI using the `display.ipynb` Jupyter notebook. - ---- diff --git a/examples/rca/offline/acquire_offline.py b/examples/rca/offline/acquire_offline.py deleted file mode 100644 index 59e4c6a..0000000 --- a/examples/rca/offline/acquire_offline.py +++ /dev/null @@ -1,195 +0,0 @@ -# import libs.... -# -## we use the acquire function to acquire data bypassing the GUI with a given -## parameters configuration and to save the IQ beamformed data as .mat files - -from collections import deque -import time -import pickle -import threading -import sys -from utils import * - - -class CaptureBuffer: - - def __init__(self, size): - self.buffer = [] - self.size = size - self.buffer_is_ready = threading.Event() - self.last_time = time.time() - - def callback(self, element): - now = time.time() - dt = now-self.last_time - print(f"Data size: {element.data.nbytes/(2**20)} MiB", end="\r") - - if len(self.buffer) < self.size: - # Copy the numpy array to the buffer - self.buffer.append(element.data.copy()) - elif len(self.buffer) == self.size: - # If this is the last element in the buffer -- unblock the main thread. - self.buffer_is_ready.set() - element.release() - - def wait_for_data(self): - self.buffer_is_ready.wait() - - def get_data(self): - return np.stack(self.buffer) - - -def acquire(medium, nframes, batch_size, sequence_xy, sequence_yx, tgc_sampling_points, tgc_values, voltage, - system_buffer_size=4): - """ - Acquires raw RF data (before IQ demodulation) for the given parameters (sequences, medium, etc.). - - :return: pair: capture buffer (with the RF data) and ARRUS metadata - """ - - if nframes % batch_size != 0: - raise ValueError("n frames should be a multiple of batch size") - - - with arrus.Session("6464.prototxt", medium=medium) as sess: - us4r = sess.get_device("/Us4R:0") - us4r.set_hv_voltage(voltage) - - system_sequence = get_system_sequence( - sequence_xy=sequence_xy, - sequence_yx=sequence_yx, - probe_model=us4r.get_probe_model(), - device_sampling_frequency=us4r.current_sampling_frequency, - ) - - system_sequence = dataclasses.replace(system_sequence, n_repeats=batch_size) - - scheme = Scheme( - tx_rx_sequence=system_sequence, - rx_buffer_size=system_buffer_size, - output_buffer=DataBufferSpec(type="FIFO", n_elements=system_buffer_size), - work_mode="SYNC" - ) - - # Upload the scheme on the us4r-lite device. - device_buffer, metadata = sess.upload(scheme) - capture_buffer = CaptureBuffer(size=nframes / batch_size) - device_buffer.append_on_new_data_callback(capture_buffer.callback) - us4r.set_tgc((tgc_sampling_points, tgc_values)) - - sess.start_scheme() - capture_buffer.wait_for_data() - sess.stop_scheme() - return capture_buffer, metadata - - -def get_time_intervals(data, metadata): - """ - This function determines the actual time intervals (frame rate) between BATCHES, i.e. scales it properly - according to the given number of batches. - - :return: returns pair: beamformed HRIs and the time intervals between consecutive frames (calculated per batch) - """ - batch_size = metadata.context.sequence.n_repeats - start_sample, end_sample = metadata.context.sequence.ops[0].rx.sample_range - n_samples = end_sample-start_sample - fs = metadata.context.device.sampling_frequency - - n_frames, total_n_samples, n_channels = data.shape - - # Determine timestamps -> time intervals. - frame_metadata = data[..., 0, 4:8] - timestamps = frame_metadata.view(np.uint64).copy() - timestamps = timestamps.flatten() - timestamps = timestamps / fs - timestamps = timestamps.squeeze() - time_intervals = np.diff(timestamps) - return time_intervals / batch_size - - -if __name__ == "__main__": - if len(sys.argv) < 2: - raise ValueError("The name of the output dataset should be provided") - - name = sys.argv[1].strip() - medium = arrus.medium.Medium(name="tissue", speed_of_sound=1450) - angles = np.linspace(-20, 20, 9) * np.pi / 180 # [rad] NUMBER OF PLANE WAVES FOR COMPOUNDING - center_frequency = 7e6 # [Hz] - n_periods = 12 - sample_range = (0, 8*1024) # @65 MHz #5*1024 # DEPTH DIMENSION - pri = 240e-6 # 119e-6 - batch_size = 14 - nframes = batch_size*10 # Should be a multiple of batch_size. - voltage = 5 - - - # Image reconstruction parameters. - y_grid = np.arange(-10e-3, 10e-3, 0.2e-3) # np.arange(-10e-3, 10e-3, 0.2e-3) ###### VOLUME of VIEW - x_grid = np.arange(-10e-3, 10e-3, 0.2e-3) # np.arange(-10e-3, 10e-3, 0.2e-3) - z_grid = np.arange(0e-3, 70e-3, 0.2e-3) # np.arange(0e-3, 30e-3, 0.1e-3) - - tgc_sampling_points = np.linspace(np.min(z_grid), np.max(z_grid), 10) - tgc_values = np.linspace(24, 34, 10) - - sequence_xy, sequence_yx = create_sequence( - medium=medium, - angles=angles, - n_periods=n_periods, - center_frequency=center_frequency, - sample_range=sample_range, - pri=pri - ) - - # times 2, because 2 TX/RXs must be performed per angle. - expected_acquisition_time = 2 * (np.sum([op.pri for op in sequence_xy.ops]) + np.sum([op.pri for op in sequence_yx.ops])) - expected_frame_rate = 1 / expected_acquisition_time - - - print(f"Acquiring data, expected frame: {expected_frame_rate:.2f} Hz") - capture_buffer, metadata = acquire( - medium=medium, - nframes=nframes, - sequence_xy=sequence_xy, - sequence_yx=sequence_yx, - voltage=voltage, - tgc_sampling_points=tgc_sampling_points, - tgc_values=tgc_values, - batch_size=batch_size, - ) - - rf = capture_buffer.get_data() - - print("Acquisition done.") - time_intervals = get_time_intervals( - data=rf, - metadata=metadata - ) - - actual_frame_rate_per_batch = 1 / time_intervals - - print(f"Actual frame rate per batch (minimum, maxium): {np.min(actual_frame_rate_per_batch):.2f}, {np.max(actual_frame_rate_per_batch):.2f} Hz") - - - filename = f"{name}_raw_dataset.pkl" - print(f"Saving data to file: {filename}...") - - with open(filename, 'wb') as f: - pickle.dump( - { - "expected_frame_rate": expected_frame_rate, - "time_intervals": time_intervals, - "rf": rf, - "metadata": metadata, - "sequence_xy": sequence_xy, - "sequence_yx": sequence_yx, - "tgc_sampling_points": tgc_sampling_points, - "tgc_values": tgc_values, - "voltage": voltage - }, - f - ) - print("Done") - - - - diff --git a/examples/rca/offline/display.ipynb b/examples/rca/offline/display.ipynb deleted file mode 100644 index efefe0b..0000000 --- a/examples/rca/offline/display.ipynb +++ /dev/null @@ -1,132 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "id": "d7106f7d-e67d-4be6-a1f6-e6a76963c1ce", - "metadata": {}, - "outputs": [], - "source": [ - "# only for the `create_animation` function\n", - "!pip install git+https://github.com/pjarosik/pjtools" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "933d0411-773a-4fa7-9ea5-cdb61c6e3d81", - "metadata": {}, - "outputs": [], - "source": [ - "import h5py\n", - "import numpy as np\n", - "from pjtools.visualization import create_animation\n", - "from IPython.display import HTML" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c4f88645-39d9-453a-9adc-7522c05f2229", - "metadata": {}, - "outputs": [], - "source": [ - "filename = \"test_phantom_beamformed_dataset.h5\"" - ] - }, - { - "cell_type": "markdown", - "id": "939ccf28-5b2c-49ef-8e25-933605b9599d", - "metadata": {}, - "source": [ - "# RF data" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e4457fc8-8baf-413f-a49b-6d1eb5981901", - "metadata": {}, - "outputs": [], - "source": [ - "f = h5py.File(filename, \"r\")\n", - "rfs = f[\"rf\"]\n", - "rfs.shape # (number of frames, number of transmits, number of samples, number of RX elements)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1a03f0cd-1212-4073-91d3-9e48705bdbde", - "metadata": {}, - "outputs": [], - "source": [ - "anim = create_animation(rfs[0, :, :600, :], value_range=(-100, 100), cmap=\"viridis\")\n", - "HTML(anim.to_jshtml())" - ] - }, - { - "cell_type": "markdown", - "id": "0d86e873-ccf1-4552-9b6c-518205910278", - "metadata": {}, - "source": [ - "# HRI/B-mode" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c5fbb22b-9f61-4411-b368-de70071739b2", - "metadata": {}, - "outputs": [], - "source": [ - "hri = f[\"hri\"] \n", - "n_frames, ny, nx, nz = hri.shape # (number of frames, n pixels oy, ox, oz)\n", - "hri.shape" - ] - }, - { - "cell_type": "markdown", - "id": "063c4141-5916-458f-9ac6-26619b01aece", - "metadata": {}, - "source": [ - "Slice y = 0." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1ea6e1a8-4501-4f3a-bb68-e1b39cc24f68", - "metadata": {}, - "outputs": [], - "source": [ - "s = hri[:, ny//2, ...] # (n frames, nx, nz)\n", - "s = s.transpose(0, 2, 1) # (n frames, nz, nx)\n", - "bmodes = 20*np.log10(np.abs(s)+1e-9)\n", - "anim = create_animation(bmodes, value_range=(20, 80), cmap=\"gray\")\n", - "HTML(anim.to_jshtml())" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.13" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/examples/rca/offline/reconstruct.py b/examples/rca/offline/reconstruct.py deleted file mode 100644 index 5af8370..0000000 --- a/examples/rca/offline/reconstruct.py +++ /dev/null @@ -1,182 +0,0 @@ -from utils import * -import cupy as cp -import sys -import h5py - - -class BatchedReconstructHriRca(Operation): - def __init__(self, **kwargs): - self.op = ReconstructHriRca(**kwargs) - - def prepare(self, metadata): - single_frame_input_shape = list(metadata.input_shape) - single_frame_input_shape[0] = 1 # batch size - single_frame_input_shape = tuple(single_frame_input_shape) - - single_batch_metadata = metadata.copy(input_shape=single_frame_input_shape) - output_metadata = self.op.prepare(single_batch_metadata) - return output_metadata.copy(input_shape=(metadata.input_shape[0], ) + output_metadata.input_shape[1:]) - - def process(self, data): - result = [] - for frame_rf in data: - frame = self.op.process(frame_rf)[0] - result.append(frame) - return cp.stack(result) - - -def to_hri( - y_grid, x_grid, z_grid, - array_tx, array_rx, - sequence, - fir_taps, name=None, - aperture_arrangement="alternate" -): - return ( - Transpose(axes=(0, 1, 3, 2)), - FirFilter(taps=fir_taps), - DigitalDownConversion(decimation_factor=4, fir_order=128), - BatchedReconstructHriRca( - y_grid=y_grid, x_grid=x_grid, z_grid=z_grid, - probe_tx=array_tx, - probe_rx=array_rx, - sequence=sequence, - min_tang=-0.5, max_tang=0.5, - name=name, - arrangement=aperture_arrangement, - ), - ) - - -def get_hri_reconstruction( - y_grid, x_grid, z_grid, - array_x: probe_params.ProbeArray, array_y: probe_params.ProbeArray, - sequence_xy: TxRxSequence, - sequence_yx: TxRxSequence, - fir_taps -): - """ - Creates a pipeline for HRI reconstruction. - """ - seqs = sequence_xy, sequence_yx - range_xy, range_yx = get_frame_ranges(*seqs) - ap_size_xy, ap_size_yx = get_rx_aperture_size(*seqs) - - - reconstruction = Pipeline( - steps=( - branch( - name="b", - steps=( - # TX=OY, RX=OX - *reorder_rf(frames=range_yx, aperture_size=ap_size_yx), - Output(), - *to_hri( - fir_taps=fir_taps, - y_grid=y_grid, - x_grid=x_grid, - z_grid=z_grid, - array_tx=array_y, - array_rx=array_x, - sequence=sequence_yx, - ), - )), - # TX=OX, RX=OY - *reorder_rf(frames=range_xy, aperture_size=ap_size_xy), - Output(), - *to_hri( - fir_taps=fir_taps, - y_grid=y_grid, - x_grid=x_grid, - z_grid=z_grid, - array_tx=array_x, array_rx=array_y, - sequence=sequence_xy, - ) - ), - placement="/GPU:0" - ) - return reconstruction - - -def reconstruct_and_save(output_hdf5, dataset, x_grid, y_grid, z_grid, fir_taps): - rf = dataset["rf"] # (n batches, ...) - # NOTE: this is a list of metadata, each subsequence metadata corresponds to separate probe - # (in that case -- we have only a single probe connected, so just take element [0]). - metadata = dataset["metadata"] - sequence_xy = dataset["sequence_xy"] - sequence_yx = dataset["sequence_yx"] - pipeline = get_hri_reconstruction( - y_grid=y_grid, x_grid=x_grid, z_grid=z_grid, - array_x=ARRAY_X, array_y=ARRAY_Y, # source: utils.py - sequence_xy=sequence_xy, sequence_yx=sequence_yx, - fir_taps=fir_taps - ) - - # As the ReconstructHriRca does not support batched processing, we need to modify - # metadata a little bit in order to - metadata = pipeline.prepare(metadata) - - hri_shape = metadata[0].input_shape # (batch size, nz, ny, nx) - hri_dtype = metadata[0].dtype - - rf_shape = list(metadata[1].input_shape) # (batch size, nTX, nsamples, n elements) - - rf_shape[1] *= 2 - rf_shape = tuple(rf_shape) # (batch size, 2*nTX, nsamples, n elements), 2*nTx, because we have XY and YX transmits - rf_dtype = metadata[1].dtype - - batch_size = metadata[0].input_shape[0] - n_batches = rf.shape[0] - n_frames = batch_size*n_batches - - rfs = output_hdf5.create_dataset("rf", shape=(n_frames, ) + rf_shape[1:], dtype=rf_dtype) - hris = output_hdf5.create_dataset("hri", shape=(n_frames, ) + hri_shape[1:], dtype=hri_dtype) - - current_size = 0 - - for i, batch in enumerate(rf): - print(f"Batch: {i}/{len(rf)}", end="\r") - output_batch = pipeline.process(cp.asarray(batch)) - hri_xy, hri_yx = output_batch[0].get(), output_batch[2].get() - rf_xy, rf_yx = output_batch[1].get(), output_batch[3].get() - # Average HRI for XY and YX transmissions. - hri = np.mean((hri_xy, hri_yx), axis=0) # (batch size, nz, ny, nx) - # Concatenate RFs for XY and YX transmissions. - rf_logical = np.stack((rf_xy, rf_yx)) # (2, batch size, nTx, nsamples, n elements) - rf_logical = rf_logical.transpose((1, 0, 2, 3, 4)) # (batch size, 2, nTX, n samples, n elements) - batch_size, _, n_tx, n_samples, n_elements = rf_logical.shape - rf_logical = rf_logical.reshape(batch_size, -1, n_samples, n_elements) - # Save to file. - print("Saving batch...", end="\r") - rfs[current_size:(current_size+batch_size)] = rf_logical - hris[current_size:(current_size+batch_size)] = hri - current_size += batch_size - - -if __name__ == "__main__": - y_grid = np.arange(-10e-3, 10e-3, 0.2e-3) # np.arange(-10e-3, 10e-3, 0.2e-3) ###### VOLUME of VIEW - x_grid = np.arange(-10e-3, 10e-3, 0.2e-3) # np.arange(-10e-3, 10e-3, 0.2e-3) - z_grid = np.arange(0e-3, 70e-3, 0.2e-3) # np.arange(0e-3, 30e-3, 0.1e-3) - - if len(sys.argv) < 2: - raise ValueError("The name of the raw dataset should be provided.") - - print("Reading data") - name = sys.argv[1] - filename = f"{name}_raw_dataset.pkl" - with open(filename, "rb") as f: - dataset = pickle.load(f) - - center_frequency = dataset["metadata"].context.sequence.ops[0].tx.excitation.center_frequency - sampling_frequency = dataset["metadata"].context.device.sampling_frequency - fir_taps = scipy.signal.firwin( - numtaps=64, cutoff=np.array([0.5, 1.5]) * center_frequency, - pass_zero="bandpass", fs=sampling_frequency - ) - - output_filename = f"{name}_beamformed_dataset.h5" - print(f"Reconstructing and saving to file: {output_filename}") - with h5py.File(output_filename, "w") as f: - reconstruct_and_save(output_hdf5=f, dataset=dataset, x_grid=x_grid, y_grid=y_grid, z_grid=z_grid, fir_taps=fir_taps) - print() - print("Done") diff --git a/examples/rca/offline/utils.py b/examples/rca/offline/utils.py deleted file mode 100644 index c645481..0000000 --- a/examples/rca/offline/utils.py +++ /dev/null @@ -1,273 +0,0 @@ -import time -import numpy as np -import arrus.medium -import arrus_rca_utils.probe_params as probe_params -from arrus.devices.probe import ProbeModel -from arrus.ops.imaging import * -from arrus.ops.us4r import * -from arrus.ops.us4r import TxRx, Tx, Rx, Aperture, Pulse -from arrus.ops.us4r import ( - TxRxSequence -) -import pickle -from arrus.utils.imaging import * -from arrus_rca_utils.reconstruction import ( - ReconstructHriRca, - GetFramesForRange, - Concatenate, - get_frame_ranges, - get_rx_aperture_size, - PipelineSequence, - SelectBatch, - Slice -) -from collections import deque -from arrus_rca_utils.sequence import convert_to_system_sequence, RcaSequence -from gui4us.model.envs.arrus import ( - UltrasoundEnv, - ArrusEnvConfiguration, - Curve -) - - -ARRAY_X = probe_params.ProbeArray( - n_elements=64, - pitch=0.28e-3, - start=0, - arrangement="ox" -) -ARRAY_Y = probe_params.ProbeArray( - n_elements=64, - pitch=0.28e-3, - start=64, - arrangement="oy" -) - - -# ------------------------------------------ SEQUENCE -def create_sequence( - medium: arrus.medium.Medium, - center_frequency: float, - n_periods: float, - angles: np.ndarray, - pri: float, - sample_range: Tuple[int, int], -) -> Tuple[TxRxSequence, TxRxSequence]: - """ - Returns TX/RX sequences for RCA probe. - - This function returns two sequences: - - transmit with OX elements, receive with OY elements len(angles) times - - transmit with OY elements, receive with OX elements len(angles) times - - First TX=OX, RX=OY len(angles) times, then TX=OY, RX=OX len(angles) times. - - NOTE: this function assumes the parameters as defined in the probe_params. - You have to adjust the values there in order to get the proper sequence - for your RCA probe. - """ - # Transmit with OX elements, receive with OY elements - sequence_xy = TxRxSequence( - ops=[ - TxRx( - tx=Tx( - # Transmit with all elements. - # Center == 0.0 means the center of the probe. - aperture=Aperture(center=0.0, size=ARRAY_X.n_elements), - excitation=Pulse( - center_frequency=center_frequency, - n_periods=n_periods, - inverse=False - ), - # Plane wave (focus depth = inf). - focus=np.inf, - angle=angle, - speed_of_sound=medium.speed_of_sound - ), - rx=Rx( - # Receive with all elements. - aperture=Aperture(center=0.0, size=ARRAY_Y.n_elements), - sample_range=sample_range, - ), - pri=pri - ) - for angle in angles - ], - ) - sequence_yx = TxRxSequence( - ops=[ - TxRx( - tx=Tx( - # Transmit with all elements. - # Center == 0.0 means the center of the probe. - aperture=Aperture(center=0.0, size=ARRAY_Y.n_elements), - excitation=Pulse( - center_frequency=center_frequency, - n_periods=n_periods, - inverse=False - ), - # Plane wave (focus depth = inf). - focus=np.inf, - angle=angle, - speed_of_sound=medium.speed_of_sound - ), - rx=Rx( - # Receive with all elements. - aperture=Aperture(center=0.0, size=ARRAY_X.n_elements), - sample_range=sample_range, - ), - pri=pri - ) - for angle in angles - ], - ) - return sequence_xy, sequence_yx - - -def get_system_sequence( - sequence_xy: TxRxSequence, sequence_yx: TxRxSequence, - probe_model: ProbeModel, - device_sampling_frequency: float -): - return convert_to_system_sequence( - sequences=[ - RcaSequence( - sequence=sequence_xy, - tx=ARRAY_X, rx=ARRAY_Y, - ), - RcaSequence( - sequence=sequence_yx, - tx=ARRAY_Y, rx=ARRAY_X, - ), - ], - probe_model=probe_model, - device_sampling_frequency=device_sampling_frequency - ) - - -# ------------------------------------------ RECONSTRUCTION -def reorder_rf(frames, aperture_size): - return ( - GetFramesForRange(frames=frames, aperture_size=aperture_size), - RemapToLogicalOrder() - ) - -def to_hri( - y_grid, x_grid, z_grid, - array_tx, array_rx, - sequence, - fir_taps, name=None, - aperture_arrangement="alternate" -): - return ( - Transpose(axes=(0, 1, 3, 2)), - FirFilter(taps=fir_taps), - DigitalDownConversion(decimation_factor=4, fir_order=128), - ReconstructHriRca( - y_grid=y_grid, x_grid=x_grid, z_grid=z_grid, - probe_tx=array_tx, - probe_rx=array_rx, - sequence=sequence, - min_tang=-0.5, max_tang=0.5, - name=name, - arrangement=aperture_arrangement, - ), - ) - - -def to_bmode(dr_min, dr_max): - return ( - # Concatenate along TX axis - Concatenate(axis=0), - Mean(axis=0), - Squeeze(), - EnvelopeDetection(), - LogCompression(), - # DynamicRangeAdjustment(min=dr_min, max=dr_max), - Squeeze(), - ) - - -def branch(steps, name=None): - return Pipeline(steps=steps, name=name, placement="/GPU:0") - - -def get_pwi_reconstruction( - y_grid, x_grid, z_grid, - array_x: probe_params.ProbeArray, array_y: probe_params.ProbeArray, - sequence_xy: TxRxSequence, - sequence_yx: TxRxSequence, - fir_taps, - dr_min=20, dr_max=80, -): - seqs = sequence_xy, sequence_yx - range_xy, range_yx = get_frame_ranges(*seqs) - ap_size_xy, ap_size_yx = get_rx_aperture_size(*seqs) - - - reconstruction = Pipeline( - steps=( - branch( - name="b", - steps=( - # TX=OY, RX=OX - *reorder_rf(frames=range_yx, aperture_size=ap_size_yx), - *to_hri( - fir_taps=fir_taps, - y_grid=y_grid, - x_grid=x_grid, - z_grid=z_grid, - array_tx=array_y, - array_rx=array_x, - sequence=sequence_yx, - ), - )), - # TX=OX, RX=OY - *reorder_rf(frames=range_xy, aperture_size=ap_size_xy), - *to_hri( - fir_taps=fir_taps, - y_grid=y_grid, - x_grid=x_grid, - z_grid=z_grid, - array_tx=array_x, array_rx=array_y, - sequence=sequence_xy, - ) - ), - placement="/GPU:0" - ) - - last_time = time.time() - def log_frame_rate(data): - nonlocal last_time - print(f"Current frame rate: {1.0/(time.time() - last_time)}") - last_time = time.time() - return data - - # Expected input: - bmodes = Pipeline( - name="bmode", - steps=( - SelectBatch([0, 1]), - *to_bmode( - dr_min=dr_min, - dr_max=dr_max - ), - # (y, x, z) - Output(), - branch( - steps=( - Slice(axis=0), - Squeeze(), - Transpose(), - )), - branch( - steps=( - Slice(axis=1), - Squeeze(), - Transpose(), - )) - ), - placement="/GPU:0" - ) - pipelines = [reconstruction, bmodes] - return PipelineSequence(pipelines) diff --git a/examples/rca/rf/app.py b/examples/rca/rf/app.py new file mode 100644 index 0000000..c255b83 --- /dev/null +++ b/examples/rca/rf/app.py @@ -0,0 +1,3 @@ +# General application settings. + +CAPTURE_BUFFER_SIZE = 3 diff --git a/examples/rca/rf/display.py b/examples/rca/rf/display.py new file mode 100644 index 0000000..61d5369 --- /dev/null +++ b/examples/rca/rf/display.py @@ -0,0 +1,29 @@ +from gui4us.cfg.display import * + +# Display configuration file. +displays = { + "XY": Display2D( + title=f"XY sequence (one frame)", + layers=( + Layer2D( + value_range=(-100, 100), + cmap="viridis", + input=StreamDataId("default", 0), + ), + ), + ax_labels=("channel #", "sample #") + ), + "YX": Display2D( + title=f"YX sequence (one frame)", + layers=( + Layer2D( + value_range=(-100, 100), + cmap="viridis", + input=StreamDataId("default", 1), + ), + ), + ax_labels=("channel #", "sample #") + ), +} + +VIEW_CFG = ViewCfg(displays) diff --git a/examples/rca/rf/env.py b/examples/rca/rf/env.py new file mode 100644 index 0000000..cec5be9 --- /dev/null +++ b/examples/rca/rf/env.py @@ -0,0 +1,163 @@ +from typing import Tuple + +import numpy as np +import arrus.medium +from arrus.devices.probe import ProbeModel +from arrus.ops.imaging import * +from arrus.ops.us4r import * +from arrus.ops.us4r import TxRx, Tx, Rx, Aperture, Pulse +from arrus.ops.us4r import ( + TxRxSequence +) +from arrus.utils.imaging import * +from gui4us.model.envs.arrus import ( + UltrasoundEnv, + ArrusEnvConfiguration, + Curve +) + +# ------------------------------------------ SEQUENCE +def create_pw_sequence( + medium: arrus.medium.Medium, + center_frequency: float, + n_periods: float, + angles: np.ndarray, + pri: float, + sample_range: Tuple[int, int], + array_tx, array_rx, + array_tx_id: str, array_rx_id: str, + name: str +) -> Tuple[TxRxSequence, TxRxSequence]: + """ + Returns PW sequence for RCA probe. + """ + # Transmit with OX elements, receive with OY elements + return TxRxSequence( + ops=[ + TxRx( + tx=Tx( + # Transmit with all elements. + # Center == 0.0 means the center of the probe. + aperture=Aperture(center=0.0, size=array_tx.n_elements), + excitation=Pulse( + center_frequency=center_frequency, + n_periods=n_periods, + inverse=False + ), + # Plane wave (focus depth = inf). + focus=np.inf, + angle=angle, + speed_of_sound=medium.speed_of_sound, + placement=array_tx_id + ), + rx=Rx( + # Receive with all elements. + aperture=Aperture(center=0.0, size=array_tx.n_elements), + sample_range=sample_range, + placement=array_tx_id + ), + pri=pri + ) + for angle in angles + ], + name=name + ) + + +def create_processing(name: str, angles): + return Pipeline( + steps=( + RemapToLogicalOrder(), + Squeeze(), + SelectFrames([len(angles)//2]), + Squeeze(), + ), + placement="/GPU:0", + name=name + ) + + + +def configure(session: arrus.Session): + us4r = session.get_device("/Us4R:0") + + # The ordinal number of the OX and OY sub-arrays. See us4r.prototxt + ox_ordinal = 0 + oy_ordinal = 1 + array_ox = us4r.get_probe_model(ordinal=ox_ordinal) + array_oy = us4r.get_probe_model(ordinal=oy_ordinal) + array_ox_id = f"Probe:{ox_ordinal}" + array_oy_id = f"Probe:{oy_ordinal}" + + sampling_frequency = us4r.sampling_frequency + n_samples = 4*1024 + + medium = arrus.medium.Medium(name="ats549", speed_of_sound=1450) + angles = np.linspace(-10, 10, 32) * np.pi / 180 # [rad] + center_frequency = 6e6 # [Hz] + n_periods = 2 + sample_range = (0, n_samples) + pri = 400e-6 + + common_parameters = dict( + medium=medium, + angles=angles, + n_periods=n_periods, + center_frequency=center_frequency, + sample_range=sample_range, + pri=pri, + ) + + # Initial TGC curve. + tgc_sampling_points = np.linspace(0, n_samples / sampling_frequency, 10) + tgc_values = np.linspace(14, 44, 10) + + # TX with OX elements, RX with OY elements + sequence_xy = create_pw_sequence( + array_tx_id=array_ox_id, + array_tx=array_ox, + array_rx_id=array_oy_id, + array_rx=array_oy, + name="XY", + **common_parameters + ) + sequence_yx = create_pw_sequence( + array_tx_id=array_oy_id, + array_tx=array_oy, + array_rx_id=array_ox_id, + array_rx=array_ox, + name="YX", + **common_parameters + ) + + preprocess_rf_xy = create_processing("preprocess_XY", angles) + preprocess_rf_yx = create_processing("preprocess_YX", angles) + + processing = Graph( + operations={preprocess_rf_xy, preprocess_rf_yx}, + dependencies={ + "Output:0": preprocess_rf_xy.name, + "Output:1": preprocess_rf_yx.name, + preprocess_rf_xy.name: sequence_xy.name, + preprocess_rf_yx.name: sequence_yx.name + } + ) + scheme = Scheme( + tx_rx_sequence=[ + sequence_xy, + sequence_yx + ], + processing=processing + ) + return ArrusEnvConfiguration( + medium=medium, + scheme=scheme, + tgc=Curve(points=tgc_sampling_points, values=tgc_values), + voltage=5 + ) + + +ENV = UltrasoundEnv( + session_cfg=str(Path(__file__).parent / "us4r.prototxt"), + configure=configure, +) diff --git a/examples/rca/rf/us4r.prototxt b/examples/rca/rf/us4r.prototxt new file mode 100644 index 0000000..9bf2b18 --- /dev/null +++ b/examples/rca/rf/us4r.prototxt @@ -0,0 +1,93 @@ +us4r: { + probe: [ + { + id: { + manufacturer: "Vermon" + name: "RCA128-OX" + } + n_elements: 128, + pitch: 0.2e-3, + tx_frequency_range: { + begin: 1e6, + end: 15e6 + }, + voltage_range { + begin: 0, + end: 50 + } + }, + { + id: { + manufacturer: "Vermon" + name: "RCA128-OY" + } + n_elements: 128, + pitch: 0.2e-3, + tx_frequency_range: { + begin: 1e6, + end: 15e6 + }, + voltage_range { + begin: 0, + end: 50 + } + } + ] + + adapter_id: { + manufacturer: "us4us" + name: "dlp408r" + } + + probe_to_adapter_connection: [ + { + probe_model_id: { + manufacturer: "Vermon" + name: "RCA128-OX" + } + probe_adapter_model_id: { + manufacturer: "us4us", + name: "dlp408r" + } + channel_mapping_ranges: [ + { + begin: 0 + end: 127 + } + ] + }, + { + probe_model_id: { + manufacturer: "Vermon" + name: "RCA128-OY" + } + probe_adapter_model_id: { + manufacturer: "us4us", + name: "dlp408r" + } + channel_mapping_ranges: [ + { + begin: 128 + end: 255 + } + ] + } + ] + + # Default initial values. + rx_settings: { + lna_gain: 24 + pga_gain: 30 + lpf_cutoff: 15000000 + active_termination: 200 + } + + hv: { + model_id { + manufacturer: "us4us" + name: "us4oemhvps" + } + } +} + + diff --git a/examples/rca/utils/arrus_rca_utils/probe_params.py b/examples/rca/utils/arrus_rca_utils/probe_params.py deleted file mode 100644 index d4a4794..0000000 --- a/examples/rca/utils/arrus_rca_utils/probe_params.py +++ /dev/null @@ -1,25 +0,0 @@ -from dataclasses import dataclass -from arrus.devices.probe import ProbeModel, ProbeModelId - - -@dataclass(frozen=True) -class ProbeArray: - """ - RCA probe array description. - - :param n_elements: number of elements in the given array - :param pitch: the distance between consecutive elements - :param start: the number of system channel connected to the first element - of the array - """ - n_elements: int - pitch: float - start: int - arrangement: str - - def to_arrus_probe(self): - return ProbeModel( - model_id=ProbeModelId("us4us", f"RCA_{self.start}"), - n_elements=self.n_elements, - pitch=self.pitch, curvature_radius=0 - ) \ No newline at end of file diff --git a/examples/rca/utils/arrus_rca_utils/reconstruction.py b/examples/rca/utils/arrus_rca_utils/reconstruction.py index 982cd59..7ece947 100644 --- a/examples/rca/utils/arrus_rca_utils/reconstruction.py +++ b/examples/rca/utils/arrus_rca_utils/reconstruction.py @@ -1,101 +1,175 @@ -from arrus_rca_utils.sequence import ( - _get_system_sequence_alternate_arrangement, - _get_system_sequence_same_arrangement -) -import copy +from typing import Tuple, List + from arrus.utils.imaging import * -import itertools import cupy as cp +import numpy as np from arrus.ops.us4r import TxRxSequence -import arrus_rca_utils.probe_params as probe_params import dataclasses import arrus.metadata - - -def get_frame_ranges(*seqs): - lengths = [len(s.ops) for s in seqs] - r = list(itertools.accumulate(lengths)) - l = [0] + r[:-1] - return zip(l, r) - - -def get_rx_aperture_size(*seqs): - return tuple(s.ops[0].rx.aperture.size for s in seqs) - - -def _get_param_name(op_name: str, param_name: str): - param_name = param_name.strip() - if not param_name.startswith("/"): - param_name = f"/{param_name}" - return f"/{op_name}{param_name}" - - -class GetFramesForRange(Operation): - - def __init__(self, frames, aperture_size): - super().__init__() - self.frame_numbers = frames - self.aperture_size = aperture_size - - def prepare(self, const_metadata): - # new FCM - old_fcm = const_metadata.data_description.custom["frame_channel_mapping"] - frame_slice = slice(*self.frame_numbers) - frames = old_fcm.frames[frame_slice, :self.aperture_size] - channels = old_fcm.channels[frame_slice, :self.aperture_size] - us4oems = old_fcm.us4oems[frame_slice, :self.aperture_size] - new_fcm = dataclasses.replace(old_fcm, frames=frames, - channels=channels, us4oems=us4oems) - new_custom_fields = copy.deepcopy(const_metadata.data_description.custom) - new_custom_fields["frame_channel_mapping"] = new_fcm - new_data_desc = dataclasses.replace( - const_metadata.data_description, custom=new_custom_fields) - return const_metadata.copy(data_desc=new_data_desc) - - def process(self, data): - return data +import arrus +import scipy.signal +import arrus.devices.us4r + + + +def get_reconstruction_graph(fir_taps: np.ndarray, + x_grid: np.ndarray, y_grid: np.ndarray, + z_grid: np.ndarray, + sequence_xy: arrus.ops.us4r.TxRxSequence, + sequence_yx: arrus.ops.us4r.TxRxSequence, + slice: bool = False): + """ + Returns RCA volume reconstruction graph. + + The assumption is that the TX=OX, RX=OY and TX=OY, RX=OX sequences are + used. The output of the processing will be a 3D volume with + (y_grid, x_grid, z_grid) coordinates. + + :param us4r: a handle to the Us4R device + :param fir_taps: FIR filter taps to be used on the pre-processing step + :param x_grid: OX coordinates + :param y_grid: OY coordinates + :param z_grid: OZ coordinates + :param sequence_xy: sequence where OX elements transmits, OY elements receive + :param sequence_yx: sequence where OY elements transmits, OX elements receive + :param slice: True if two 2D B-mode images (x = 0 and y = 0) slices should be produced, False otherwise + """ + + preprocess_xy = create_pre_processing(name="preprocess_xy") + preprocess_yx = create_pre_processing(name="preprocess_yx") + + to_hri_xy = to_hri( + name="to_hri_xy", + x_grid=x_grid, y_grid=y_grid, z_grid=z_grid, + fir_taps=fir_taps, + tx_orientation="ox", + rx_orientation="oy", + ) + to_hri_yx = to_hri( + name="to_hri_yx", + x_grid=x_grid, y_grid=y_grid, z_grid=z_grid, + fir_taps=fir_taps, + tx_orientation="oy", + rx_orientation="ox", + ) + to_bmode_all = to_bmode(name="to_bmode") + + if slice: + to_slice_xz = to_slice(axis=0, name="to_slice_xz") + to_slice_yz = to_slice(axis=1, name="to_slice_yz") + + operations = { + preprocess_xy, preprocess_yx, + to_bmode_all, + to_hri_xy, to_hri_yx, + } + deps = { + "Output:0": to_bmode_all.name, + to_bmode_all.name: [to_hri_xy.name, to_hri_yx.name], + to_hri_yx.name: preprocess_yx.name, + to_hri_xy.name: preprocess_xy.name, + preprocess_xy.name: sequence_xy.name, + preprocess_yx.name: sequence_yx.name + } + if slice: + # Append volume slicing for 2D display. + operations.add(to_slice_xz) + operations.add(to_slice_yz) + + deps["Output:0"] = to_slice_xz.name + deps["Output:1"] = to_slice_yz.name + deps[to_slice_yz.name] = to_bmode_all.name + deps[to_slice_xz.name] = to_bmode_all.name + + return Graph( + operations=operations, + dependencies=deps + ) + + +def create_pre_processing(name: str): + return Pipeline( + steps=( + RemapToLogicalOrder(), + ), + placement="/GPU:0", + name=name + ) + +def to_hri( + y_grid, x_grid, z_grid, + fir_taps, + tx_orientation: str, + rx_orientation: str, + name: str, + zeros=False +): + return Pipeline( + steps=( + Transpose(axes=(0, 1, 3, 2)), + FirFilter(taps=fir_taps), + DigitalDownConversion(decimation_factor=4, fir_order=64), + ReconstructHriRca( + y_grid=y_grid, x_grid=x_grid, z_grid=z_grid, + min_tang=-0.5, max_tang=0.5, + name=name, + tx_orientation=tx_orientation, + rx_orientation=rx_orientation + ), + Lambda(lambda data: cp.zeros_like(data) if zeros else data) + ), + placement="/GPU:0", + name=name + ) + + +def to_bmode(name: str): + return Pipeline( + steps=( + Concatenate(axis=0, name="concat"), # Concatenate XY and YX sequences + Mean(axis=0), # Average HRIs (coherently) + Squeeze(), + EnvelopeDetection(), + LogCompression(), + Squeeze(), + ), + placement="/GPU:0", + name=name + ) + + +def to_slice(axis, name: str): + return Pipeline( + steps=( + Slice(axis=axis), + Transpose() + ), + placement="/GPU:0", + name=name + ) class Concatenate(Operation): - def __init__(self, axis): - super().__init__() + def __init__(self, axis, name: str = None): + super().__init__(name=name) self.axis = axis - def prepare(self, const_metadata): + def prepare(self, const_metadata: List[arrus.metadata.ConstMetadata]): # TODO verify that all the metadata objects are compatible, i.e. have the same data type, etc. axis_total_size = 0 for cm in const_metadata: axis_total_size += cm.input_shape[self.axis] output_shape = list(const_metadata[0].input_shape) output_shape[self.axis] = axis_total_size - return const_metadata[0].copy(input_shape=tuple(output_shape)) + res = const_metadata[0].copy(input_shape=tuple(output_shape)) + print("!!!!!!!!!!!!!!!!!!!! METADATA") + print(res) + return res - def process(self, data): + def process(self, data: Tuple[cp.ndarray, cp.ndarray]): return cp.concatenate(data, axis=self.axis) -class SelectBatch(Operation): - - def __init__(self, indices): - super().__init__() - if not isinstance(indices, Iterable): - indices = [indices] - self.indices = indices - - def prepare(self, const_metadata): - result = [const_metadata[i] for i in self.indices] - if len(result) == 1: - return result[0] - else: - return result - - def process(self, data): - result = [data[i] for i in self.indices] - if len(result) == 1: - return result[0] - return result - - class Slice(Operation): def __init__(self, axis, position=None): @@ -136,78 +210,36 @@ def process(self, data): return data[self.slicing] -class PipelineSequence(Pipeline): - - def __init__(self, pipelines): - # Intentionally not calling super constructor. - self.pipelines = pipelines - self._set_names() - self._determine_params() - - def prepare(self, metadata): - m = metadata - for p in self.pipelines: - m = p.prepare(m) - return m - - def process(self, data): - d = data - for p in self.pipelines: - d = p.process(d) - return d - - def __call__(self, data): - return self.process(data) - - def set_parameter(self, key: str, value: Sequence[Number]): - pipeline, pipeline_param_name = self._pipeline_params[key] - pipeline.set_parameter(pipeline_param_name, value) - - def get_parameter(self, key: str) -> Sequence[Number]: - pipeline, pipeline_param_name = self._pipeline_params[key] - return pipeline.get_parameter(pipeline_param_name) - - def get_parameters(self) -> Dict[str, ParameterDef]: - return self._param_defs - - def _set_names(self): - for i, pipeline in enumerate(self.pipelines): - if not hasattr(pipeline, "name") or pipeline.name is None: - pipeline.name = f"Pipeline:{i}" - - def _determine_params(self): - self._pipeline_params = {} - self._param_defs = {} - for pipeline in self.pipelines: - name = pipeline.name - params = pipeline.get_parameters() - for k, param_def in params.items(): - prefixed_k = _get_param_name(name, k) - self._pipeline_params[prefixed_k] = pipeline, k - self._param_defs[prefixed_k] = param_def - - class ReconstructHriRca(Operation): def __init__( - self, y_grid, x_grid, z_grid, - probe_tx: probe_params.ProbeArray, - probe_rx: probe_params.ProbeArray, - sequence, arrangement="alternate", + self, + y_grid: np.ndarray, + x_grid: np.ndarray, + z_grid: np.ndarray, + tx_orientation: str, + rx_orientation: str, min_tang=-0.5, max_tang=0.5, name=None ): + """ + RCA beamformer. + + The beamformer reconstructs High-resolution Image (i.e. + compounds LRIs). + + :param tx_orientation: the orientation of TX probe, one of: "ox", "oy" + :param rx_orientation: the orientation of RX probe, one of: "ox", "oy" + """ super().__init__(name) self.y_grid = y_grid self.x_grid = x_grid self.z_grid = z_grid - self.sequence = sequence - self.array_rx = probe_rx - self.array_tx = probe_tx + self.array_rx_orientation = rx_orientation + self.array_tx_orientation = tx_orientation if min_tang > max_tang: - raise ValueError( - f"min tang {min_tang} should be <= max tang {max_tang}") + raise ValueError(f"min tang {min_tang} " + f"should be <= max tang {max_tang}") self.min_tang = min_tang self.max_tang = max_tang - self.arrangement = arrangement self.n_sigma = 3.0 self.num_pkg = cp @@ -219,8 +251,8 @@ def set_pkgs(self, num_pkg, **kwargs): def set_parameter(self, key: str, value: Sequence[Number]): if not hasattr(self, key): raise ValueError(f"{type(self).__name__} has no {key} parameter.") - # TODO assumming scalar parameters only - setattr(self, key, self.num_pkg.float32(value)) + # TODO assuming scalar parameters only + setattr(self, key, cp.float32(value)) def get_parameter(self, key: str) -> Sequence[Number]: if not hasattr(self, key): @@ -251,32 +283,23 @@ def prepare(self, const_metadata): current_dir = os.path.dirname(os.path.join(os.path.abspath(__file__))) kernel_path = Path(current_dir) / "reconstruct.cu" kernel_source = kernel_path.read_text() - self.kernel_module = self.num_pkg.RawModule(code=kernel_source) + self.kernel_module = cp.RawModule(code=kernel_source) self.kernel = self.kernel_module.get_function("iqRaw2Hri") - # INPUT PARAMETERS. # Input data shape. self.n_seq, self.n_tx, self.n_rx, self.n_samples = const_metadata.input_shape if self.n_seq > 1: - raise ValueError("At most 1 sequence is supported.") - - if self.arrangement == "alternate": - raw_seq, tx_cent_delay = _get_system_sequence_alternate_arrangement( - sequence=self.sequence, - array_tx=self.array_tx, array_rx=self.array_rx, - device_sampling_frequency=const_metadata.context.device.sampling_frequency - ) - elif self.arrangement in {"same_x", "same_y"}: - raw_seq, tx_cent_delay = _get_system_sequence_same_arrangement( - sequence=self.sequence, - array=self.array_tx, - arrangement=self.arrangement, - metadata=const_metadata - ) - else: - raise ValueError(f"Unknown aperture orientation: {self.arrangement}") + raise ValueError("At most 1 TX/RX sequence is supported by this beamformer") + self.sequence = const_metadata.context.sequence + self.raw_sequence = const_metadata.context.raw_sequence reference_op = self.sequence.ops[0] + array_tx_id = self.sequence.get_tx_probe_id_unique() + array_rx_id = self.sequence.get_rx_probe_id_unique() + + self.array_tx = const_metadata.context.device.get_probe_by_id(array_tx_id).model + self.array_rx = const_metadata.context.device.get_probe_by_id(array_rx_id).model + start_sample = reference_op.rx.sample_range[0] speed_of_sound = reference_op.tx.speed_of_sound pulse = reference_op.tx.excitation @@ -287,8 +310,8 @@ def prepare(self, const_metadata): output_shape = (1, self.y_size, self.x_size, self.z_size) - self.output_buffer = self.num_pkg.zeros( - output_shape, dtype=self.num_pkg.complex64 + self.output = cp.zeros( + output_shape, dtype=cp.complex64 ) y_block_size = min(self.y_size, 8) x_block_size = min(self.x_size, 8) @@ -299,75 +322,60 @@ def prepare(self, const_metadata): int((self.x_size - 1) // x_block_size + 1), int((self.y_size - 1) // y_block_size + 1)) - self.y_pix = self.num_pkg.asarray( - self.y_grid, dtype=self.num_pkg.float32) - self.x_pix = self.num_pkg.asarray( - self.x_grid, dtype=self.num_pkg.float32) - self.z_pix = self.num_pkg.asarray( - self.z_grid, dtype=self.num_pkg.float32) + self.y_pix = cp.asarray(self.y_grid, dtype=cp.float32) + self.x_pix = cp.asarray(self.x_grid, dtype=cp.float32) + self.z_pix = cp.asarray(self.z_grid, dtype=cp.float32) # System and transmit properties. - self.sos = self.num_pkg.float32(speed_of_sound) - self.fs = self.num_pkg.float32( - const_metadata.data_description.sampling_frequency) - self.fn = self.num_pkg.float32(pulse.center_frequency) + self.sos = cp.float32(speed_of_sound) + self.fs = cp.float32(const_metadata.data_description.sampling_frequency) + self.fn = cp.float32(pulse.center_frequency) - self.probe_tx_pitch = self.num_pkg.float32(self.array_tx.pitch) - self.probe_tx_n_elements = self.num_pkg.int32(self.array_tx.n_elements) - tx_arrangement = self.get_arrangement_code(self.array_tx) - self.array_tx_arrangement = self.num_pkg.uint8(tx_arrangement) + self.probe_tx_pitch = cp.float32(self.array_tx.pitch) + self.probe_tx_n_elements = cp.int32(self.array_tx.n_elements) + tx_orientation = self.get_arrangement_code(self.array_tx_orientation) + self.array_tx_orientation = cp.uint8(tx_orientation) - self.probe_rx_pitch = self.num_pkg.float32(self.array_rx.pitch) - self.probe_rx_n_elements = self.num_pkg.int32(self.array_rx.n_elements) - rx_arrangement = self.get_arrangement_code(self.array_rx) - self.array_rx_arrangement = self.num_pkg.uint8(rx_arrangement) + self.probe_rx_pitch = cp.float32(self.array_rx.pitch) + self.probe_rx_n_elements = cp.int32(self.array_rx.n_elements) + rx_arrangement = self.get_arrangement_code(self.array_rx_orientation) + self.array_rx_orientation = cp.uint8(rx_arrangement) # TX focus and angle tx_focus = (op.tx.focus for op in self.sequence.ops) - tx_focus = [foc if foc is not None else np.inf for foc in tx_focus] - self.tx_focus = self.num_pkg.asarray( - tx_focus, - dtype=self.num_pkg.float32 - ) - + # TODO: this assumption, i.e. foc None means inf should be moved to ARRUS imaging + tx_focus = [foc if foc is not None else cp.inf for foc in tx_focus] + self.tx_focus = cp.asarray(tx_focus, dtype=cp.float32) + # TODO: this assumption, i.e. foc None means inf should be moved to ARRUS imaging tx_angle = (op.tx.angle for op in self.sequence.ops) tx_angle = [angle if angle is not None else 0.0 for angle in tx_angle] - self.tx_angle = self.num_pkg.asarray( - tx_angle, - dtype=self.num_pkg.float32 - ) + self.tx_angle = cp.asarray(tx_angle, dtype=cp.float32) # TX aperture centers. - tx_ap_cent = [ - self.get_ap_center(op.tx.aperture, self.array_tx) - for op in self.sequence.ops - ] - self.tx_ap_cent = self.num_pkg.asarray( - tx_ap_cent, - dtype=self.num_pkg.float32 - ) + tx_ap_cent = [self.get_ap_center(op.tx.aperture, self.array_tx) + for op in self.sequence.ops] + self.tx_ap_cent = cp.asarray(tx_ap_cent, dtype=cp.float32) + # first/last probe element in TX aperture - tx_ap_begin, tx_ap_end = self.get_tx_ap_bounds(raw_seq) - self.tx_ap_begin = self.num_pkg.asarray( - tx_ap_begin, - dtype=self.num_pkg.float32 - ) - self.tx_ap_end = self.num_pkg.asarray( - tx_ap_end, - dtype=self.num_pkg.float32 - ) - rx_ap_origin = self.get_rx_ap_origin(raw_seq) - self.rx_ap_origin = self.num_pkg.asarray( - rx_ap_origin, - dtype=self.num_pkg.int32) + tx_ap_begin, tx_ap_end = self.get_tx_ap_bounds(self.raw_sequence) + self.tx_ap_begin = cp.asarray(tx_ap_begin, dtype=cp.float32) + self.tx_ap_end = cp.asarray(tx_ap_end, dtype=cp.float32) + + rx_ap_origin = self.get_rx_ap_origin(self.raw_sequence) + self.rx_ap_origin = cp.asarray(rx_ap_origin, dtype=cp.int32) # Min/max tang - self.min_tang = self.num_pkg.float32(self.min_tang) - self.max_tang = self.num_pkg.float32(self.max_tang) + self.min_tang = cp.float32(self.min_tang) + self.max_tang = cp.float32(self.max_tang) burst_factor = pulse.n_periods / (2 * self.fn) - self.initial_delay = -start_sample/const_metadata.context.device.sampling_frequency - self.initial_delay = burst_factor + tx_cent_delay - self.initial_delay = self.num_pkg.float32(self.initial_delay) + tx_cent_delay = arrus.kernels.tx_rx_sequence.get_center_delay( + sequence=self.sequence, + probe_tx=self.array_tx, + probe_rx=self.array_rx, + ) + self.init_delay = -start_sample / const_metadata.context.device.sampling_frequency + self.init_delay = burst_factor + tx_cent_delay + self.init_delay = cp.float32(self.init_delay) # Output metadata new_signal_description = dataclasses.replace( const_metadata.data_description, @@ -381,9 +389,9 @@ def prepare(self, const_metadata): ) def process(self, data): - data = self.num_pkg.ascontiguousarray(data) + data = cp.ascontiguousarray(data) params = ( - self.output_buffer, + self.output, data, self.n_tx, self.n_samples, self.n_rx, self.z_pix, self.z_size, @@ -393,32 +401,31 @@ def process(self, data): self.tx_focus, self.tx_angle, self.tx_ap_cent, self.probe_tx_pitch, self.probe_tx_n_elements, - self.array_tx_arrangement, + self.array_tx_orientation, self.probe_rx_pitch, self.probe_rx_n_elements, - self.array_rx_arrangement, + self.array_rx_orientation, self.tx_ap_begin, self.tx_ap_end, self.rx_ap_origin, self.min_tang, self.max_tang, - self.initial_delay, + self.init_delay, self.n_sigma ) self.kernel(self.grid_size, self.block_size, params) - return self.output_buffer + return self.output - def get_arrangement_code(self, array: probe_params.ProbeArray): - return 0 if array.arrangement == "ox" else 1 + def get_arrangement_code(self, orientation: str): + return 0 if orientation == "ox" else 1 def get_tx_ap_bounds(self, raw_sequence: TxRxSequence): first_els = [] last_els = [] for op in raw_sequence.ops: ap_elems = np.argwhere(op.tx.aperture) - first = np.min(ap_elems)-self.array_tx.start - last = np.max(ap_elems)-self.array_tx.start + first, last = np.min(ap_elems), np.max(ap_elems) first_els.append(first) last_els.append(last) first_els, last_els = np.asarray(first_els), np.asarray(last_els) - probe_model = self.array_tx.to_arrus_probe() + probe_model = self.array_tx pos_x = probe_model.element_pos_x return pos_x[first_els], pos_x[last_els] @@ -426,15 +433,17 @@ def get_rx_ap_origin(self, raw_sequence: TxRxSequence): origins = [] for op in raw_sequence.ops: ap_elems = np.argwhere(op.rx.aperture) - first = np.min(ap_elems) - self.array_rx.start + first = np.min(ap_elems) origins.append(first) return np.asarray(origins) - def get_ap_center( - self, - aperture: arrus.ops.us4r.Aperture, - probe: probe_params.ProbeArray - ): + def get_ap_center(self, aperture: arrus.ops.us4r.Aperture, probe): + """ + Returns the physical location [m] of the aperture center for the + given probe. + + :param probe: probe model description (arrus.devices.probe.ProbeModel) + """ if aperture.center is None and aperture.center_element is None: return 0.0 # default - center of the probe if aperture.center is not None: diff --git a/examples/rca/utils/arrus_rca_utils/sequence.py b/examples/rca/utils/arrus_rca_utils/sequence.py index eb77de4..24bc87f 100644 --- a/examples/rca/utils/arrus_rca_utils/sequence.py +++ b/examples/rca/utils/arrus_rca_utils/sequence.py @@ -1,221 +1,56 @@ -from typing import Tuple, Iterable -import dataclasses -from dataclasses import dataclass - -import numpy as np +from typing import Tuple import arrus.medium -from arrus.ops.us4r import ( - TxRxSequence, TxRx, Tx, Rx, Aperture, Pulse -) -from arrus.devices.probe import ProbeModel, ProbeModelId -from arrus.kernels.tx_rx_sequence import ( - convert_to_us4r_sequence -) -import arrus_rca_utils.probe_params as probe_params - +from arrus.ops.us4r import * +import numpy as np -@dataclass(frozen=True) -class RcaSequence: - """ - :param arrangement: available values: alternate (xy or yx), same_x, same_y +def get_pw_sequence( + medium: arrus.medium.Medium, + center_frequency: float, + n_periods: float, + angles: np.ndarray, + pri: float, + sample_range: Tuple[int, int], + array_tx, array_rx, + array_tx_id: str, array_rx_id: str, + name: str +) -> Tuple[TxRxSequence, TxRxSequence]: """ - sequence: TxRxSequence - tx: probe_params.ProbeArray - rx: probe_params.ProbeArray - arrangement: str = "alternate" + Returns PW sequence for RCA probe. - -def convert_to_system_sequence_for_metadata( - sequences: Iterable[RcaSequence], - metadata -): - return convert_to_system_sequence( - sequences=sequences, - probe_model=metadata.context.device.probe.model, - device_sampling_frequency=metadata.context.device.sampling_frequency - ) - - -def convert_to_system_sequence( - sequences: Iterable[RcaSequence], - probe_model: ProbeModel, - device_sampling_frequency, -): - seqs = [] - for s in sequences: - if s.arrangement in {"same_x", "same_y"}: - raw_sequence, _ = _get_system_sequence_same_arrangement( - sequence=s.sequence, array=s.tx, - arrangement=s.arrangement, - device_sampling_frequency=device_sampling_frequency, - probe_model=probe_model - ) - elif s.arrangement == "alternate": - raw_sequence, _ = _get_system_sequence_alternate_arrangement( - sequence=s.sequence, - array_tx=s.tx, array_rx=s.rx, - device_sampling_frequency=device_sampling_frequency + :param array_tx: the probe to be used on TX (e.g. us4r.get_probe_model(i)) + :param array_rx: the probe to be used on RX (e.g. us4r.get_probe_model(i)) + :param array_tx_id: the id of the probe to be used on TX (ordinal number, e.g. `i`) + :param array_rx_id: the id of the probe to be used on RX (ordinal number e.g. `i`) + """ + return TxRxSequence( + ops=[ + TxRx( + tx=Tx( + # Transmit with all elements. + # Center == 0.0 means the center of the probe. + aperture=Aperture(center=0.0, size=array_tx.n_elements), + excitation=Pulse( + center_frequency=center_frequency, + n_periods=n_periods, + inverse=False + ), + # Plane wave (focus depth = inf). + focus=np.inf, + angle=angle, + speed_of_sound=medium.speed_of_sound, + placement=array_tx_id + ), + rx=Rx( + # Receive with all elements. + aperture=Aperture(center=0.0, size=array_rx.n_elements), + sample_range=sample_range, + placement=array_rx_id + ), + pri=pri ) - else: - raise ValueError(f"Unknown orientation: {s.arrangement}") - seqs.append(raw_sequence) - system_sequence = __concatenate_sequences_all(seqs) - return __equalize_rx_aperture(system_sequence) - - -def _get_system_sequence_same_arrangement( - sequence: TxRxSequence, - array: probe_params.ProbeArray, - arrangement: str, - probe_model: ProbeModel, - device_sampling_frequency: float, -): - total_n_elements = probe_model.n_elements - arrus_pm = array.to_arrus_probe() - raw_sequence, tx_delay = convert_to_us4r_sequence( - sequence=sequence, - probe_model=arrus_pm, - fs=device_sampling_frequency - ) - if arrangement == "same_x": - should_append_left = True - elif arrangement == "same_y": - should_append_left = False - else: - raise ValueError(f"Unknown arrangement: {arrangement}") - raw_sequence = __extend_aperture( - sequence=raw_sequence, - n_missing_elements=total_n_elements - array.n_elements, - is_append_left=should_append_left - ) - return raw_sequence, tx_delay - - -def _get_system_sequence_alternate_arrangement( - sequence, - array_tx: probe_params.ProbeArray, - array_rx: probe_params.ProbeArray, - device_sampling_frequency: float -): - seq_array_tx, seq_array_rx = __split_sequence_between_arrays( - sequence=sequence, array_tx=array_tx, array_rx=array_rx - ) - arrus_array_tx = array_tx.to_arrus_probe() - arrus_array_rx = array_rx.to_arrus_probe() - raw_seq_tx, init_delay = convert_to_us4r_sequence( - sequence=seq_array_tx, - probe_model=arrus_array_tx, - fs=device_sampling_frequency - ) - raw_seq_rx, _ = convert_to_us4r_sequence( - sequence=seq_array_rx, - probe_model=arrus_array_rx, - fs=device_sampling_frequency - ) - if array_rx.start < array_tx.start: - seqs = (raw_seq_rx, raw_seq_tx) - else: - seqs = (raw_seq_tx, raw_seq_rx) - return __merge_apertures(*seqs), init_delay - - -def __split_sequence_between_arrays( - sequence: TxRxSequence, - array_tx: probe_params.ProbeArray, - array_rx: probe_params.ProbeArray -): - array_tx_ops = [] - array_rx_ops = [] - for op in sequence.ops: - tx: Tx = op.tx - rx: Rx = op.rx - - empty_rx_aperture = np.zeros(array_tx.n_elements).astype(bool) - empty_rx = dataclasses.replace(rx, aperture=empty_rx_aperture) - array_tx_op = dataclasses.replace(op, tx=tx, rx=empty_rx) - array_tx_ops.append(array_tx_op) - - empty_tx_aperture = np.zeros(array_rx.n_elements).astype(bool) - empty_tx = dataclasses.replace(tx, aperture=empty_tx_aperture) - probe_rx_op = dataclasses.replace(op, tx=empty_tx, rx=rx) - array_rx_ops.append(probe_rx_op) - - seq_array_tx = dataclasses.replace( - sequence, - ops=array_tx_ops, + for angle in angles + ], + name=name ) - seq_array_rx = dataclasses.replace( - sequence, - ops=array_rx_ops - ) - return seq_array_tx, seq_array_rx - - -def __merge_apertures(a: TxRxSequence, b: TxRxSequence): - new_ops = [] - assert len(a.ops) == len(b.ops) - - for op1, op2 in zip(a.ops, b.ops): - new_rx = dataclasses.replace( - op1.rx, - aperture=np.concatenate((op1.rx.aperture, op2.rx.aperture)) - ) - if np.sum(op1.tx.aperture) > 0: - base_tx = op1.tx - else: - base_tx = op2.tx - new_tx = dataclasses.replace( - base_tx, - aperture=np.concatenate((op1.tx.aperture, op2.tx.aperture)) - ) - new_op = dataclasses.replace(op1, tx=new_tx, rx=new_rx) - new_ops.append(new_op) - - return dataclasses.replace(a, ops=new_ops) - - -def __extend_aperture( - sequence: TxRxSequence, n_missing_elements: int, is_append_left: bool -): - result = [] - empty_aperture = np.zeros(n_missing_elements).astype(bool) - - for op in sequence.ops: - if is_append_left: - rx_ap = np.concatenate((empty_aperture, op.rx.aperture)) - tx_ap = np.concatenate((empty_aperture, op.tx.aperture)) - else: - rx_ap = np.concatenate((op.rx.aperture, empty_aperture)) - tx_ap = np.concatenate((op.tx.aperture, empty_aperture)) - - new_rx = dataclasses.replace(op.rx, aperture=rx_ap) - new_tx = dataclasses.replace(op.tx, aperture=tx_ap) - new_op = dataclasses.replace(op, tx=new_tx, rx=new_rx) - result.append(new_op) - return dataclasses.replace(sequence, ops=result) - - -def __concatenate_sequences_all(seqs): - ops = [] - for s in seqs: - ops.extend(s.ops) - return dataclasses.replace(seqs[0], ops=ops) - - -def __equalize_rx_aperture(seq: TxRxSequence): - new_ops = [] - max_ap_size = np.max([np.sum(op.rx.aperture) for op in seq.ops]) - for op in seq.ops: - rx = op.rx - rx_ap_n_elements = np.sum(rx.aperture) - padding_right = max_ap_size - rx_ap_n_elements - assert rx.padding == (0, 0), "Currently padded RX is not supported" - new_padding = (0, padding_right) - new_rx = dataclasses.replace(rx, padding=new_padding) - new_op = dataclasses.replace(op, rx=new_rx) - new_ops.append(new_op) - return dataclasses.replace(seq, ops=new_ops) - -