Hi everybody!
I would like to ask you some questions about the simualtion of ultrasound B-mode images starting from 'Simulating B-mode Ultrasound Images Example".
I have tried to reproduce it using only one ball . In my example the ball demarcates a region containing Cerebral Spinal Fluid so I define the typical density and speed of sound inside that region equal to --> density:1007 and speed of sound:1504.5;
Outside the ball I imagine having the properties of brain tissue so: speed of sound:1546.3 and Density :1046.
as regards the remaining part of the code I tried to follow the example provided.
My questions are :
1 ) I am not sure of having correctly selected all the parameters , in particular, I have a matrix that defines me the speed of sound in the medium and one that defines the density but when I create the input signal is required the use of a single value. In my example I used the maximum speed of sound value and the minimum density value inside the matrices but I'm not sure if it is correct .
this line : input_signal = (source_strength./(1546.3*1007)).*input_signal;
2 ) I do not know how to properly plot the resulting image . I tried to download the "scan_lines matrix " provided at the start of the example to figure out how to get the final image of the three balls but I do not get anything recognizable , i suppose something is wrong in my postprocessing.
3 ) the value of the attenuation coefficient alpha should be different between the two regions (inside and outside the ball ) or it should be a scalar , for example a mean value?
I put my code here :
Thank you very much for any help and suggestion.
clear all
clc
close all
% simulation settings
DATA_CAST = 'single';
% create the computational grid
Nx = 100;
Ny = 20;
Nz=40;
% number of grid points in the y (column) direction
dx = 0.1e-3; % grid point spacing in the x direction [m]
dy = 0.1e-3; % grid point spacing in the y direction [m]
dz=0.1e-3;
kgrid = makeGrid(Nx, dx, Ny,dy,Nz,dz);
object=makeBall(Nx, 100, Nz, Nx/2, Ny/2, Nz/2, 10, true);
medium.sound_speed = 1546.3*ones(Nx,100,Nz); % [m/s]
medium.density = 1046*ones(Nx,100,Nz); % [kg/m^3]
medium.sound_speed(object==1)=1504.5; %inside the ball
medium.density(object==1) =1007;
% define properties of the input signal
source_strength = 1e6; % [Pa]
tone_burst_freq = 0.5e6; % [Hz]
tone_burst_cycles = 5;
% create the time array
% kgrid.t_array = makeTime(kgrid, medium.sound_speed, [], t_end);
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
% create the input signal using toneBurst
input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles);
% scale the source magnitude by the source_strength divided by the
% impedance (the source is assigned to the particle velocity)
input_signal = (source_strength./(1546.3*1007)).*input_signal;
plot(input_signal);
% append input signal used to drive the transducer
transducer.input_signal = input_signal;
% physical properties of the transducer
transducer.number_elements =19; % total number of transducer elements
transducer.element_width = 1; % width of each element [grid points]
transducer.element_length = 8; % length of each element [grid points]
transducer.element_spacing = 0; % spacing (kerf width) between the elements [grid points]
transducer.radius = inf; % radius of curvature of the transducer [m]
% calculate the width of the transducer in grid points
transducer_width = transducer.number_elements*transducer.element_width ...
+ (transducer.number_elements - 1)*transducer.element_spacing;
% use this to position the transducer in the middle of the computational grid
transducer.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]);
% properties used to derive the beamforming delays
% properties used to derive the beamforming delays
transducer.sound_speed = 1504.5; % sound speed [m/s]%WHICH VALUE OF THE MATRIX SHOULD I PUT?
transducer.focus_distance = Nx/4; % focus distance [m]
transducer.elevation_focus_distance = 19e-3;% focus distance in the elevation plane [m]
transducer.steering_angle = 0; % steering angle [degrees]
transducer.transmit_apodization = 'Rectangular';
transducer.receive_apodization = 'Rectangular';
% define the transducer elements that are currently active
transducer.active_elements = zeros(transducer.number_elements, 1);
transducer.active_elements(1:19) = 1;
% % append input signal used to drive the transducer
transducer.input_signal = input_signal;
% create the transducer using the defined settings
transducer = makeTransducer(kgrid, transducer);
transducer.plot;
transducer.properties;
t=medium.sound_speed;
u= medium.density;
number_scan_lines=79;
medium_position=[1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2];
% run the simulation
for scan_line_index = 1:number_scan_lines
% update the command line status
disp('');
disp(['Computing scan line ' num2str(scan_line_index) ' of ' num2str(number_scan_lines)]);
% load the current section of the medium
medium.sound_speed = t(:, medium_position(2):medium_position(2)+ Ny-1 , :);
medium.density= u(:, medium_position(2):medium_position(2) + Ny-1 , :);
% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer)%,input_args{:});
% extract the scan line from the sensor data
scan_lines(scan_line_index, :) = transducer.scan_line(sensor_data);
% update medium position
medium_position = medium_position + [0 transducer.element_width 0];
end
c0=1540 ; %WHICH VALUE SHOULD I PUT?
% create radius variable assuming that t0 corresponds to the middle of the
% input signal
t0 = length(input_signal)*kgrid.dt/2;
r = c0*( (1:length(kgrid.t_array))*kgrid.dt/2 - t0); % [m]
% create time gain compensation function based on attenuation value,
% transmit frequency, and round trip distance
tgc_alpha = 0.6; % [dB/(MHz cm)] %IT's THE ATTENUATION OUTSIDE THE BALL, HOW CAN I TAKE IN ACCOUNT ALSO THE %ALPHA OF THE INSIDE MEDIUM?
tgc = exp(2*tgc_alpha*tone_burst_freq/1e6*r*100);
% apply the time gain compensation to each of the scan lines
scan_lines = bsxfun(@times, tgc, scan_lines);
% filter the scan lines using both the transmit frequency and the second harmonic
scan_lines_fund = gaussianFilter(scan_lines, 1/kgrid.dt, tone_burst_freq, 100, true);
scan_lines_harm = gaussianFilter(scan_lines, 1/kgrid.dt, 2*tone_burst_freq, 30, true);
% envelope detection
scan_lines_fund = envelopeDetection(scan_lines_fund);
scan_lines_harm = envelopeDetection(scan_lines_harm);
% normalised log compression
compression_ratio = 3;
scan_lines_fund = logCompression(scan_lines_fund, compression_ratio, true);
scan_lines_harm = logCompression(scan_lines_harm, compression_ratio, true);
% upsample the image using linear interpolation
scale_factor = 2;
scan_lines_fund = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_fund, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).');
scan_lines_harm = interp2(1:kgrid.Nt, (1:number_scan_lines).', scan_lines_harm, 1:kgrid.Nt, (1:1/scale_factor:number_scan_lines).');