SHARE
TWEET

B-mode linear transducer

a guest Oct 21st, 2019 124 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. % Simulating B-mode Ultrasound Images Example
  2. %
  3. % This example illustrates how k-Wave can be used for the simulation of
  4. % B-mode ultrasound images (including tissue harmonic imaging) analogous to
  5. % those produced by a modern diagnostic ultrasound scanner. It builds on
  6. % the Defining An Ultrasound Transducer, Simulating Ultrasound Beam
  7. % Patterns, and Using An Ultrasound Transducer As A Sensor examples.
  8. %
  9. % Note, this example generates a B-mode ultrasound image from a 3D
  10. % scattering phantom using kspaceFirstOrder3D. Compared to ray-tracing or
  11. % Field II, this approach is very general. In particular, it accounts for
  12. % nonlinearity, multiple scattering, power law acoustic absorption, and a
  13. % finite beam width in the elevation direction. However, it is also
  14. % computationally intensive. Using a modern GPU and the Parallel Computing
  15. % Toolbox (with 'DataCast' set to 'gpuArray-single'), each scan line takes
  16. % around 3 minutes to compute. Using a modern desktop CPU (with 'DataCast'
  17. % set to 'single'), this increases to around 30 minutes. In this example,
  18. % the final image is constructed using 96 scan lines. This makes the total
  19. % computational time around 4.5 hours using a single GPU, or 2 days using a
  20. % CPU.
  21. %
  22. % To allow the simulated scan line data to be processed multiple times with
  23. % different settings, the simulated RF data is saved to disk. This can be
  24. % reloaded by setting RUN_SIMULATION = false within the example m-file. The
  25. % data can also be downloaded from:
  26. %
  27. % http://www.k-wave.org/datasets/example_us_bmode_scan_lines.mat
  28. %
  29. % author: Bradley Treeby
  30. % date: 3rd August 2011
  31. % last update: 9th June 2017
  32. %  
  33. % This function is part of the k-Wave Toolbox (http://www.k-wave.org)
  34. % Copyright (C) 2011-2017 Bradley Treeby
  35.  
  36. % This file is part of k-Wave. k-Wave is free software: you can
  37. % redistribute it and/or modify it under the terms of the GNU Lesser
  38. % General Public License as published by the Free Software Foundation,
  39. % either version 3 of the License, or (at your option) any later version.
  40. %
  41. % k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
  42. % WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  43. % FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
  44. % more details.
  45. %
  46. % You should have received a copy of the GNU Lesser General Public License
  47. % along with k-Wave. If not, see <http://www.gnu.org/licenses/>.
  48.  
  49. %#ok<*UNRCH>
  50.  
  51. clearvars;
  52.  
  53. % simulation settings
  54. DATA_CAST       = 'single';     % set to 'single' or 'gpuArray-single' to speed up computations
  55. RUN_SIMULATION  = true;         % set to false to reload previous results instead of running simulation
  56.  
  57. % =========================================================================
  58. % DEFINE THE K-WAVE GRID
  59. % =========================================================================
  60.  
  61. % set the size of the perfectly matched layer (PML)
  62. pml_x_size = 20;                % [grid points]
  63. pml_y_size = 10;                % [grid points]
  64. pml_z_size = 10;                % [grid points]
  65.  
  66. % set total number of grid points not including the PML
  67. Nx = 256 - 2 * pml_x_size;      % [grid points]
  68. Ny = 128 - 2 * pml_y_size;      % [grid points]
  69. Nz = 128 - 2 * pml_z_size;      % [grid points]
  70.  
  71. % set desired grid size in the x-direction not including the PML
  72. x = 40e-3;                      % [m]
  73.  
  74. % calculate the spacing between the grid points
  75. dx = x / Nx;                    % [m]
  76. dy = dx;                        % [m]
  77. dz = dx;                        % [m]
  78.  
  79. % create the k-space grid
  80. kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
  81.  
  82. % =========================================================================
  83. % DEFINE THE MEDIUM PARAMETERS
  84. % =========================================================================
  85.  
  86. % define the properties of the propagation medium
  87. c0 = 1540;                      % [m/s]
  88. rho0 = 1000;                    % [kg/m^3]
  89. medium.alpha_coeff = 0.75;      % [dB/(MHz^y cm)]
  90. medium.alpha_power = 1.5;
  91. medium.BonA = 6;
  92.  
  93. % create the time array
  94. t_end = (Nx * dx) * 2.2 / c0;   % [s]
  95. kgrid.makeTime(c0, [], t_end);
  96.  
  97. % =========================================================================
  98. % DEFINE THE INPUT SIGNAL
  99. % =========================================================================
  100.  
  101. % define properties of the input signal
  102. source_strength = 1e6;          % [Pa]
  103. tone_burst_freq = 1.5e6;        % [Hz]
  104. tone_burst_cycles = 4;
  105.  
  106. % create the input signal using toneBurst
  107. input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
  108.  
  109. % scale the source magnitude by the source_strength divided by the
  110. % impedance (the source is assigned to the particle velocity)
  111. input_signal = (source_strength ./ (c0 * rho0)) .* input_signal;
  112.  
  113. % =========================================================================
  114. % DEFINE THE ULTRASOUND TRANSDUCER
  115. % =========================================================================
  116.  
  117. % physical properties of the transducer
  118. transducer.number_elements = 32;    % total number of transducer elements
  119. transducer.element_width = 2;       % width of each element [grid points]
  120. transducer.element_length = 24;     % length of each element [grid points]
  121. transducer.element_spacing = 0;     % spacing (kerf  width) between the elements [grid points]
  122. transducer.radius = inf;            % radius of curvature of the transducer [m]
  123.  
  124. % calculate the width of the transducer in grid points
  125. transducer_width = transducer.number_elements * transducer.element_width ...
  126.     + (transducer.number_elements - 1) * transducer.element_spacing;
  127.  
  128. % use this to position the transducer in the middle of the computational grid
  129. transducer.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]);
  130.  
  131. % properties used to derive the beamforming delays
  132. transducer.sound_speed = c0;                    % sound speed [m/s]
  133. transducer.focus_distance = 20e-3;              % focus distance [m]
  134. transducer.elevation_focus_distance = 19e-3;    % focus distance in the elevation plane [m]
  135. transducer.steering_angle = 0;                  % steering angle [degrees]
  136.  
  137. % apodization
  138. transducer.transmit_apodization = 'Hanning';    
  139. transducer.receive_apodization = 'Rectangular';
  140.  
  141. % define the transducer elements that are currently active
  142. number_active_elements = 32;
  143. transducer.active_elements = ones(transducer.number_elements, 1);
  144.  
  145. % append input signal used to drive the transducer
  146. transducer.input_signal = input_signal;
  147.  
  148. % create the transducer using the defined settings
  149. transducer = kWaveTransducer(kgrid, transducer);
  150.  
  151. % print out transducer properties
  152. transducer.properties;
  153.  
  154. % =========================================================================
  155. % DEFINE THE MEDIUM PROPERTIES
  156. % =========================================================================
  157.  
  158. % define a large image size to move across
  159. number_scan_lines = 96;
  160. Nx_tot = Nx;
  161. Ny_tot = Ny + number_scan_lines * transducer.element_width;
  162. Nz_tot = Nz;
  163.  
  164. % define a random distribution of scatterers for the medium
  165. background_map_mean = 1;
  166. background_map_std = 0.008;
  167. background_map = background_map_mean + background_map_std * randn([Nx_tot, Ny_tot, Nz_tot]);
  168.  
  169. % define a random distribution of scatterers for the highly scattering
  170. % region
  171. scattering_map = randn([Nx_tot, Ny_tot, Nz_tot]);
  172. scattering_c0 = c0 + 25 + 75 * scattering_map;
  173. scattering_c0(scattering_c0 > 1600) = 1600;
  174. scattering_c0(scattering_c0 < 1400) = 1400;
  175. scattering_rho0 = scattering_c0 / 1.5;
  176.  
  177. % define properties
  178. sound_speed_map = c0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;
  179. density_map = rho0 * ones(Nx_tot, Ny_tot, Nz_tot) .* background_map;
  180.  
  181. % define a sphere for a highly scattering region
  182. radius = 6e-3;      % [m]
  183. x_pos = 27.5e-3;    % [m]
  184. y_pos = 20.5e-3;    % [m]
  185. scattering_region1 = makeBall(Nx_tot, Ny_tot, Nz_tot, round(x_pos/dx), round(y_pos/dx), Nz_tot/2, round(radius/dx));
  186.  
  187. % assign region
  188. sound_speed_map(scattering_region1 == 1) = scattering_c0(scattering_region1 == 1);
  189. density_map(scattering_region1 == 1) = scattering_rho0(scattering_region1 == 1);
  190.  
  191. % define a sphere for a highly scattering region
  192. radius = 5e-3;      % [m]
  193. x_pos = 30.5e-3;    % [m]
  194. y_pos = 37e-3;      % [m]
  195. scattering_region2 = makeBall(Nx_tot, Ny_tot, Nz_tot, round(x_pos/dx), round(y_pos/dx), Nz_tot/2, round(radius/dx));
  196.  
  197. % assign region
  198. sound_speed_map(scattering_region2 == 1) = scattering_c0(scattering_region2 == 1);
  199. density_map(scattering_region2 == 1) = scattering_rho0(scattering_region2 == 1);
  200.  
  201. % define a sphere for a highly scattering region
  202. radius = 4.5e-3;    % [m]
  203. x_pos = 15.5e-3;    % [m]
  204. y_pos = 30.5e-3;    % [m]
  205. scattering_region3 = makeBall(Nx_tot, Ny_tot, Nz_tot, round(x_pos/dx), round(y_pos/dx), Nz_tot/2, round(radius/dx));
  206.  
  207. % assign region
  208. sound_speed_map(scattering_region3 == 1) = scattering_c0(scattering_region3 == 1);
  209. density_map(scattering_region3 == 1) = scattering_rho0(scattering_region3 == 1);
  210.  
  211. % =========================================================================
  212. % RUN THE SIMULATION
  213. % =========================================================================
  214.  
  215. % preallocate the storage
  216. scan_lines = zeros(number_scan_lines, kgrid.Nt);
  217.  
  218. % set the input settings
  219. input_args = {...
  220.     'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size], ...
  221.     'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', false};
  222.  
  223. % run the simulation if set to true, otherwise, load previous results from
  224. % disk
  225. if RUN_SIMULATION
  226.  
  227.     % set medium position
  228.     medium_position = 1;
  229.    
  230.     % loop through the scan lines
  231.     for scan_line_index = 1:number_scan_lines
  232.        
  233.         % update the command line status
  234.         disp('');
  235.         disp(['Computing scan line ' num2str(scan_line_index) ' of ' num2str(number_scan_lines)]);
  236.  
  237.         % load the current section of the medium
  238.         medium.sound_speed = sound_speed_map(:, medium_position:medium_position + Ny - 1, :);
  239.         medium.density = density_map(:, medium_position:medium_position + Ny - 1, :);
  240.        
  241.         % run the simulation
  242.         sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer, input_args{:});
  243.  
  244.         % extract the scan line from the sensor data
  245.         scan_lines(scan_line_index, :) = transducer.scan_line(sensor_data);
  246.        
  247.         % update medium position
  248.         medium_position = medium_position + transducer.element_width;
  249.  
  250.     end
  251.  
  252.     % save the scan lines to disk
  253.     save example_us_bmode_scan_lines scan_lines;
  254.    
  255. else
  256.    
  257.     % load the scan lines from disk
  258.     load example_us_bmode_scan_lines;
  259.    
  260. end
  261.  
  262. % =========================================================================
  263. % PROCESS THE RESULTS
  264. % =========================================================================
  265.  
  266. % -----------------------------
  267. % Remove Input Signal
  268. % -----------------------------
  269.  
  270. % create a window to set the first part of each scan line to zero to remove
  271. % interference from the input signal
  272. scan_line_win = getWin(kgrid.Nt * 2, 'Tukey', 'Param', 0.05).';
  273. scan_line_win = [zeros(1, length(input_signal) * 2), scan_line_win(1:end/2 - length(input_signal) * 2)];
  274.  
  275. % apply the window to each of the scan lines
  276. scan_lines = bsxfun(@times, scan_line_win, scan_lines);
  277.  
  278. % store a copy of the middle scan line to illustrate the effects of each
  279. % processing step
  280. scan_line_example(1, :) = scan_lines(end/2, :);
  281.  
  282. % -----------------------------
  283. % Time Gain Compensation
  284. % -----------------------------
  285.  
  286. % create radius variable assuming that t0 corresponds to the middle of the
  287. % input signal
  288. t0 = length(input_signal) * kgrid.dt / 2;
  289. r = c0 * ( (1:length(kgrid.t_array)) * kgrid.dt / 2 - t0);    % [m]
  290.  
  291. % create time gain compensation function based on attenuation value,
  292. % transmit frequency, and round trip distance
  293. tgc_alpha = 0.4;       % [dB/(MHz cm)]
  294. tgc = exp(2 * tgc_alpha * tone_burst_freq * 1e-6 * r * 100);
  295.  
  296. % apply the time gain compensation to each of the scan lines
  297. scan_lines = bsxfun(@times, tgc, scan_lines);
  298.  
  299. % store a copy of the middle scan line to illustrate the effects of each
  300. % processing step
  301. scan_line_example(2, :) = scan_lines(end/2, :);
  302.  
  303. % -----------------------------
  304. % Frequency Filtering
  305. % -----------------------------
  306.  
  307. % filter the scan lines using both the transmit frequency and the second
  308. % harmonic
  309. scan_lines_fund = gaussianFilter(scan_lines, 1/kgrid.dt, tone_burst_freq, 100, true);
  310. set(gca, 'XLim', [0, 6]);
  311. scan_lines_harm = gaussianFilter(scan_lines, 1/kgrid.dt, 2 * tone_burst_freq, 30, true);
  312. set(gca, 'XLim', [0, 6]);
  313.  
  314. % store a copy of the middle scan line to illustrate the effects of each
  315. % processing step
  316. scan_line_example(3, :) = scan_lines_fund(end/2, :);
  317.  
  318. % -----------------------------
  319. % Envelope Detection
  320. % -----------------------------
  321.  
  322. % envelope detection
  323. scan_lines_fund = envelopeDetection(scan_lines_fund);
  324. scan_lines_harm = envelopeDetection(scan_lines_harm);
  325.  
  326. % store a copy of the middle scan line to illustrate the effects of each
  327. % processing step
  328. scan_line_example(4, :) = scan_lines_fund(end/2, :);
  329.  
  330. % -----------------------------
  331. % Log Compression
  332. % -----------------------------
  333.  
  334. % normalised log compression
  335. compression_ratio = 3;
  336. scan_lines_fund = logCompression(scan_lines_fund, compression_ratio, true);
  337. scan_lines_harm = logCompression(scan_lines_harm, compression_ratio, true);
  338.  
  339. % store a copy of the middle scan line to illustrate the effects of each
  340. % processing step
  341. scan_line_example(5, :) = scan_lines_fund(end/2, :);
  342.  
  343. % -----------------------------
  344. % Scan Conversion
  345. % -----------------------------
  346.  
  347. % upsample the image using linear interpolation
  348. scale_factor = 2;
  349. scan_lines_fund = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_fund, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).');
  350. scan_lines_harm = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_harm, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).');
  351.  
  352. % =========================================================================
  353. % VISUALISATION
  354. % =========================================================================
  355.  
  356. % plot the medium, truncated to the field of view
  357. figure;
  358. imagesc((0:number_scan_lines * transducer.element_width - 1) * dy * 1e3, (0:Nx_tot-1) * dx * 1e3, sound_speed_map(:, 1 + Ny/2:end - Ny/2, Nz/2));
  359. axis image;
  360. colormap(gray);
  361. set(gca, 'YLim', [5, 40]);
  362. title('Scattering Phantom');
  363. xlabel('Horizontal Position [mm]');
  364. ylabel('Depth [mm]');
  365.  
  366. % plot the processing steps
  367. figure;
  368. stackedPlot(kgrid.t_array * 1e6, {'1. Beamformed Signal', '2. Time Gain Compensation', '3. Frequency Filtering', '4. Envelope Detection', '5. Log Compression'}, scan_line_example);
  369. xlabel('Time [\mus]');
  370. set(gca, 'XLim', [5, t_end * 1e6]);
  371.  
  372. % plot the processed b-mode ultrasound image
  373. figure;
  374. horz_axis = (0:length(scan_lines_fund(:, 1)) - 1) * transducer.element_width * dy / scale_factor * 1e3;
  375. imagesc(horz_axis, r * 1e3, scan_lines_fund.');
  376. axis image;
  377. colormap(gray);
  378. set(gca, 'YLim', [5, 40]);
  379. title('B-mode Image');
  380. xlabel('Horizontal Position [mm]');
  381. ylabel('Depth [mm]');
  382.  
  383. % plot the processed harmonic ultrasound image
  384. figure;
  385. imagesc(horz_axis, r * 1e3, scan_lines_harm.');
  386. axis image;
  387. colormap(gray);
  388. set(gca, 'YLim', [5, 40]);
  389. title('Harmonic Image');
  390. xlabel('Horizontal Position [mm]');
  391. ylabel('Depth [mm]');
  392.  
  393. % =========================================================================
  394. % VISUALISATION OF SIMULATION LAYOUT
  395. % =========================================================================
  396.  
  397. % % uncomment to generate a voxel plot of the simulation layout
  398. %
  399. % % physical properties of the transducer
  400. % transducer_plot.number_elements = 32 + number_scan_lines - 1;
  401. % transducer_plot.element_width = 2;
  402. % transducer_plot.element_length = 24;
  403. % transducer_plot.element_spacing = 0;
  404. % transducer_plot.radius = inf;
  405. %
  406. % % transducer position
  407. % transducer_plot.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]);
  408. %
  409. % % create expanded grid
  410. % kgrid_plot = kWaveGrid(Nx_tot, dx, Ny_tot, dy, Nz, dz);
  411. % kgrid_plot.setTime(1, 1);
  412. %
  413. % % create the transducer using the defined settings
  414. % transducer_plot = kWaveTransducer(kgrid_plot, transducer_plot);
  415. %
  416. % % create voxel plot of transducer mask and
  417. % voxelPlot(single(transducer_plot.active_elements_mask | scattering_region1 | scattering_region2 | scattering_region3));
  418. % view(26, 48);
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top