k-Wave User Forum » Topic: US probe frequency

k-Wave User Forum » Topic: US probe frequency http://www.k-wave.org/forum/topic/us-probe-frequency Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:38:08 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "US probe frequency" http://www.k-wave.org/forum/topic/us-probe-frequency#post-6233 Tue, 12 Dec 2017 22:39:57 +0000 Bradley Treeby 6233@http://www.k-wave.org/forum/ Hi Amanda, Unfortunately, the geometry functions in k-Wave assume the grid is uniform. You could write your own function to get around this without too much effort by using the Cartesian coordinates of the grid points, rather than the grid coordinates. However, I would suggest keeping dx = dy = dz, otherwise the grid will support a different maximum frequency in each direction. Hope that helps, Brad. Amanda on "US probe frequency" http://www.k-wave.org/forum/topic/us-probe-frequency#post-6218 Thu, 16 Nov 2017 02:01:07 +0000 Amanda 6218@http://www.k-wave.org/forum/ Hello :) I recently started studies related to reconstruction of ultrasound images and I am using kwave software (the example of the b-mode using linear transducer) to implement functions as spatial compound and pulse inversion harmonic image techniques. However, when I did the harmonic image simulations for a frequency of 3.5MHz, I got results that i think are wrong. Reading the comments in the forum and the user manual (which I'm still having difficulties to understand, because the content present there goes beyond my knowledge, and I'm still learning), I realized that there was a need to change the value of dx (mainly), but changing only dx the phantom was distorted (there were no more spheres but ellipses), even so, I could not increase the value of dy as much as I increased dx because otherwise the simulations would take too long, even using the GPU. I would like to confirm that what I understand so far is correct about the grid changes and if this distortion of the phantom really happens for different values of dx and dy? Here's the part of my code for grid changes: '% set total number of grid points not including the PML Nx = 512 - 2 * pml_x_size; % [grid points] Ny = 128 - 2 * pml_y_size; % [grid points] Nz = 128 - 2 * pml_z_size; % [grid points] % set desired grid size in the x-direction not including the PML x = 40e-3; % [m] dx_old = x/(256 - 2 * pml_x_size); % [grid points] % calculate the spacing between the grid points dx = x / Nx; % [m] dy = dx_old; % [m] dz = dx_old; % [m] ' and about the part of the code to make the scattering region: % define a sphere for a highly scattering region radius = 6e-3; % [m] x_pos = 27.5e-3; % [m] y_pos = 20.5e-3; % [m] scattering_region1 = makeBall(Nx_tot, Ny_tot, Nz_tot, round(x_pos/dx), round(y_pos/dx_old), Nz_tot/2, round(radius/dx_old)); % assign region sound_speed_map(scattering_region1 == 1) = scattering_c0(scattering_region1 == 1); density_map(scattering_region1 == 1) = scattering_rho0(scattering_region1 == 1); % define a sphere for a highly scattering region radius = 5e-3; % [m] x_pos = 30.5e-3; % [m] y_pos = 37e-3; % [m] scattering_region2 = makeBall(Nx_tot, Ny_tot, Nz_tot, round(x_pos/dx), round(y_pos/dx_old), Nz_tot/2, round(radius/dx_old)); % assign region sound_speed_map(scattering_region2 == 1) = scattering_c0(scattering_region2 == 1); density_map(scattering_region2 == 1) = scattering_rho0(scattering_region2 == 1); % define a sphere for a highly scattering region radius = 4.5e-3; % [m] x_pos = 15.5e-3; % [m] y_pos = 30.5e-3; % [m] scattering_region3 = makeBall(Nx_tot, Ny_tot, Nz_tot, round(x_pos/dx), round(y_pos/dx_old), Nz_tot/2, round(radius/dx_old)); ' thank you! Amanda marlenep on "US probe frequency" http://www.k-wave.org/forum/topic/us-probe-frequency#post-5576 Tue, 05 Jul 2016 15:20:28 +0000 marlenep 5576@http://www.k-wave.org/forum/ Hi Bradley, thank you very much for your advice I share my code and the resulting scan lines matrix with you at this link: <a href="https://imperialcollegelondon.box.com/s/vu0obiih6kn1watkojfcc8x2sheezy1o" rel="nofollow">https://imperialcollegelondon.box.com/s/vu0obiih6kn1watkojfcc8x2sheezy1o</a> I would like to obtain a simulated ultrasound image similar to real ultrasound image ("real_us.bmp"). I already tried playing with different parameters and the best result I obtained is the one you can find in the folder. Do you think there is a way to improve it or the limitation depends only on the mri segmentation which of course approximates the real brain tissues differentiation? thank you very much for your time Bradley Treeby on "US probe frequency" http://www.k-wave.org/forum/topic/us-probe-frequency#post-5557 Mon, 04 Jul 2016 12:50:06 +0000 Bradley Treeby 5557@http://www.k-wave.org/forum/ Hi marlenep, You should set the total simulation time to be long enough to capture the reflections of interest. The line <code>t_end = (Nx*dx)*2.2/c0;</code> will set the time long enough for a single round trip the length of the domain. Regarding the time step, I would refer to Section 2.7 in the k-Wave manual, where this is discussed in detail. In particular, the sound speed you pass to <code>makeTime</code> (which sets <code>dt</code>) and the value of <code>medium.sound_speed_ref</code> (which sets the value of sound speed in the k-space operator; this defaults to <code>max(medium.sound_speed(:))</code>) will both affect the simulation. >> I don't have such a good result from my simulation Can you say a bit more about what the problem is? Brad. marlenep on "US probe frequency" http://www.k-wave.org/forum/topic/us-probe-frequency#post-5529 Thu, 09 Jun 2016 10:53:08 +0000 marlenep 5529@http://www.k-wave.org/forum/ thank you Brad, I increased the PPW because and seems to be better. How can i set a proper time step and total simulation time? Is still ok if I use this command? t_end = (Nx*dx)*2.2/c0; [kgrid.t_array delta] = makeTime(kgrid, c0, [], t_end); I am trying to simulate an ultrasound image from a phantom with similar properties to those of the soft tissues and with two holes in the middle fill with water. I have acquired real ultrasound images using a linear transducer probe, in particular is the ultrasonic L14-5/38. I don't have such a good result from my simulation, can you give me some advice to improve it? thank you very much for any help This is my code ( in this code I am still using 2 PPW) clear all clc close all % simulation settings DATA_CAST = 'gpuArray-single'; %DATA_CAST='single'; RUN_SIMULATION =true; points_per_wavelength =2; freq=5e6; dx = 1482.3/(points_per_wavelength*freq); dy=dx; dz=dy; Nx = round((46e-3)/dx); % phantom dimensions Ny_tot=round((50.2e-3)/dy); Ny=66; Nz=round((50.2e-3)/dz); % set the size of the perfectly matched layer (PML) pml_x_size = 30; % [grid points] pml_y_size = 10; % [grid points] pml_z_size = 30; % [grid points] % % set total number of grid points not including the PML sc = 1; % % create the k-space grid kgrid = makeGrid(Nx, dx,Ny, dy, Nz, dz); o=round(36.5e-03/dx); p=round(36.9e-03/dx); q=round(12.9e-03/dx); r=round(7.5e-03/dx); circle =makeDisc(Nx,Ny_tot,o,q,r,false)+(makeDisc(Nx,Ny_tot,o,p,r,false)); %3.65,1.29,0.75//36.5e-03 36.9e-03 7.5e-03 rect=[]; for z=1:Nz for x=1:Nx for y=1:Ny_tot if (circle(x,y)==0) rect(x,y,z)=0; else rect(x,y,z)=1; end end end end figure,imshow(rect(:,:,20)); % define the properties of the propagation medium %first medium c0 = 1540; % [m/s]water rho0 = 1068; % [kg/m^3] medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)] medium.alpha_power = 1.5; medium.BonA = 6; % create the time array t_end = (Nx*dx)*2.2/c0; [kgrid.t_array delta] = makeTime(kgrid, c0, [], t_end); %define properties of the input signal source_strength = 0.25e3;%1e6; % [Pa] tone_burst_freq = freq/sc; % [Hz] tone_burst_cycles = 5; % 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./(c0*rho0)).*input_signal; transducer.input_signal = input_signal; % DEFINE THE ULTRASOUND TRANSDUCER % ========================================================================= uno=round((38e-03/128)/dx); % transducer size len=round(10e-03/dx); %transducer heigth % physical properties of the transducer transducer.number_elements = 32; % total number of transducer elements transducer.element_width = uno; % width of each element [grid points] transducer.element_length = len; % 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 transducer.sound_speed = c0; % sound speed [m/s] transducer.focus_distance = 6e-3; % focus distance [m] transducer.elevation_focus_distance = 19e-3;% focus distance in the elevation plane [m] transducer.steering_angle = 0; % transducer,element_pitch=0.3e-03;% steering angle [degrees] %transducer.steering_angle_max = 32; % apodization transducer.transmit_apodization = 'Hanning'; transducer.receive_apodization = 'Rectangular'; % define the transducer elements that are currently active number_active_elements = 32; transducer.active_elements = ones(transducer.number_elements, 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); % print out transducer properties transducer.properties; transducer.plot; % ========================================================================= % DEFINE THE MEDIUM PROPERTIES % ========================================================================= % define a large image size to move across number_scan_lines = 96; Nx_tot = Nx; %Ny_tot = Ny + number_scan_lines*transducer.element_width; Nz_tot = Nz; mu0_map = zeros(Nx, Ny_tot, Nz); mu1_map = zeros(Nx, Ny_tot, Nz); sigma0_map = zeros(Nx, Ny_tot, Nz); %background sound_speed_map = 1482.3*ones(Nx, Ny_tot, Nz);%.*background_map_s; density_map = 994*ones(Nx, Ny_tot, Nz);%.*background_map_s; alpha_coeff_map= ones(Nx,Ny_tot,Nz); %water sound_speed_map(rect==1)= 1482.3;%.*background_map(Vol==6); %giusto?? density_map(rect==1) = 994;%.*background_map(Vol==6); alpha_coeff_map(rect==1) =0.2; mu0_map(rect==1) = 0;%0.2; mu1_map(rect==1) = 0;%0.4; sigma0_map(rect==1) = 0; %soft tissue sound_speed_map(rect==0)= 1540;%.*background_map(Vol==6); %giusto?? density_map(rect==0) = 1068;%.*background_map(Vol==6); alpha_coeff_map(rect==0) =0.54; mu0_map(rect==0) = 0.53; mu1_map(rect==0) = 0.5; sigma0_map(rect==0) = 0; %scattering %generate noise maps T0 = randn(Nx, Ny_tot, Nz); T1 = randn(Nx, Ny_tot, Nz); % calc noise to add S = T0.*sigma0_map + mu0_map; rho_noise = zeros(Nx, Ny_tot, Nz); rho_noise(T1 <= mu1_map) = S(T1 <= mu1_map); density_map = density_map.*(1+rho_noise); % assign to the medium inputs medium.sound_speed = sound_speed_map; medium.density = density_map; medium.alpha_coeff=alpha_coeff_map; % ========================================================================= % RUN THE SIMULATION % ========================================================================= % preallocate the storage scan_lines = zeros(number_scan_lines, kgrid.Nt); % set the input settings input_args = {... 'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size], ... 'DataCast', DATA_CAST, 'DataRecast', true, 'PlotSim', false}; % run the simulation if set to true, otherwise, load previous results from % disk if RUN_SIMULATION % set medium position medium_position= round((Ny_tot-(128*transducer.element_width))/2); %medium_position = 146; %%%%%%%%%% 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 = sound_speed_map(:, medium_position:medium_position + Ny - 1, :); medium.density = density_map(:, medium_position:medium_position + Ny - 1, :); medium.alpha_coeff = alpha_coeff_map(:, medium_position:medium_position + Ny - 1, :); % run the simulation % sensor_data = kspaceFirstOrder3D(kgrid, medium, transducer, transducer,'DataCast','gpuArray-single', input_args{:}); sensor_data = kspaceFirstOrder3DG(kgrid, medium, transducer, transducer,'DataCast','gpuArray-single', 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 + transducer.element_width; end % save the scan lines to disk save example_us_bmode_scan_lines scan_lines; else % load the scan lines from disk load example_us_bmode_scan_lines end % ========================================================================= % PROCESS THE RESULTS % ========================================================================= % ----------------------------- % Remove Input Signal % ----------------------------- % create a window to set the first part of each scan line to zero to remove % interference from the input signal scan_line_win = getWin(kgrid.Nt*2, 'Tukey', 'Param', 0.05).'; scan_line_win = [zeros(1, length(input_signal)*2), scan_line_win(1:end/2 - length(input_signal)*2)]; % apply the window to each of the scan lines scan_lines = bsxfun(@times, scan_line_win, scan_lines); % store a copy of the middle scan line to illustrate the effects of each % processing step scan_line_example(1, :) = scan_lines(end/2, :); % ----------------------------- % Time Gain Compensation % ----------------------------- % 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.4; % [dB/(MHz cm)] 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); % store a copy of the middle scan line to illustrate the effects of each % processing step scan_line_example(2, :) = scan_lines(end/2, :); % ----------------------------- % Frequency Filtering % ----------------------------- % 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); % store a copy of the middle scan line to illustrate the effects of each % processing step scan_line_example(3, :) = scan_lines_fund(end/2, :); % ----------------------------- % Envelope Detection % ----------------------------- % envelope detection scan_lines_fund = envelopeDetection(scan_lines_fund); scan_lines_harm = envelopeDetection(scan_lines_harm); % store a copy of the middle scan line to illustrate the effects of each % processing step scan_line_example(4, :) = scan_lines_fund(end/2, :); % ----------------------------- % Log Compression % ----------------------------- % 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); % store a copy of the middle scan line to illustrate the effects of each % processing step scan_line_example(5, :) = scan_lines_fund(end/2, :); % ----------------------------- % Scan Conversion % ----------------------------- % 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).'); % ========================================================================= % VISUALISATION % ========================================================================= % plot the medium, truncated to the field of view figure; 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)); axis image; colormap(gray); set(gca, 'YLim', [0, 40]); title('Scattering Phantom'); xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); % plot the processing steps figure; stackedPlot(kgrid.t_array*1e6, {'1. Beamformed Signal', '2. Time Gain Compensation', '3. Frequency Filtering', '4. Envelope Detection', '5. Log Compression'}, scan_line_example); xlabel('Time [\mus]'); set(gca, 'XLim', [5 t_end*1e6]); % plot the processed b-mode ultrasound image figure; horz_axis = (0:length(scan_lines_fund(:, 1))-1)*transducer.element_width*dy/scale_factor*1e3; imagesc(horz_axis, r*1e3, scan_lines_fund.'); axis image; colormap(gray); set(gca, 'YLim', [5, 40]); title('B-mode Image'); xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); % plot the processed harmonic ultrasound image figure; imagesc(horz_axis, r*1e3, scan_lines_harm.'); axis image; colormap(gray); set(gca, 'YLim', [5, 40]); title('Harmonic Image'); xlabel('Horizontal Position [mm]'); ylabel('Depth [mm]'); Bradley Treeby on "US probe frequency" http://www.k-wave.org/forum/topic/us-probe-frequency#post-5526 Wed, 08 Jun 2016 20:57:09 +0000 Bradley Treeby 5526@http://www.k-wave.org/forum/ Hi marlenep, Possible things to check are whether the source frequency is still supported by the grid (I'd recommend using at least 3 PPW at the highest frequency contained in the tone burst), that you have also scaled the time step appropriately so the simulation is still stable, and that the total simulation time is still appropriate to capture the behaviour you are trying to model (e.g., scattering). The source strength or number of cycles don't need to be changed. Brad. marlenep on "US probe frequency" http://www.k-wave.org/forum/topic/us-probe-frequency#post-5523 Tue, 07 Jun 2016 10:33:23 +0000 marlenep 5523@http://www.k-wave.org/forum/ Hi! I would want to change the frequency of the US source from 1.5 to 3MHz. When I change it, if i plot the signal it seems that I get a signal that is not correct. In my example I define the grid with respect to the frequency so I don't think the problem could be related to the grid size. Do you know if maybe some other parameters have to be changed when you change the frequency as for example the source strenght or the number of cycles ? thank you for any help