k-Wave User Forum » Topic: image reconstruction for defects

k-Wave User Forum » Topic: image reconstruction for defects http://www.k-wave.org/forum/topic/image-reconstruction-for-defects Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:15:40 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php bencox on "image reconstruction for defects" http://www.k-wave.org/forum/topic/image-reconstruction-for-defects#post-1651 Fri, 22 Mar 2013 10:11:57 +0000 bencox 1651@http://www.k-wave.org/forum/ Hi oozunlu, The problem is with the forward modelling of the scattering, and specifically the very large contrast between the acoustic properties in the steel background and the air inclusion. k-Wave was written with only small heterogeneities in mind (say 10%), such as might be found between different soft tissue when doing ultrasound imaging. Large and sharp contrasts don't work well because the numerical method is a spectral method - it uses a Fourier transform to calculate the gradient. (Technically, the coefficients in the spatial Fourier transform of the acoustic properties do not decay fast enough, and so the solution quickly becomes inaccurate.) You could try smoothing the acoustic properties using the smooth function, or trying an inclusion closer in properties to the steel background, but you would probably be better trying a model more specifically designed for modelling scattering from what is essentially a pressure-release (sound soft) interface. ie. one in which you can actually specify that as a boundary condition. Another thing to note is that in practice you will have shear waves in steel, which are also not modelled by k-Wave which assumes the medium is a fluid. Ben oozunlu on "image reconstruction for defects" http://www.k-wave.org/forum/topic/image-reconstruction-for-defects#post-1615 Thu, 21 Mar 2013 15:38:38 +0000 oozunlu 1615@http://www.k-wave.org/forum/ I modified "2D Time Reversal Reconstruction For A Line Sensor Example" and made medium heterogeneous with the following addition: However, simulation results in an expected situation. Did I do something wrong? Is it possible to apply image reconstruction for defects? %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% (I added an air bubble into a steel block) disc_x_pos=round(Nx*0.7); disc_y_pos=round(Ny/4); disc_radius=5; air_buble = makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius); air_buble_reverse = (air_buble-1)*-1; medium.sound_speed = medium.sound_speed .* air_buble_reverse; medium.sound_speed = medium.sound_speed + air_buble * 340 ; medium.density = 7800*ones(Nx, Ny); % [kg/m^3] medium.density = medium.density .* air_buble_reverse; medium.density = medium.density + air_buble * 12.25 ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% therefore, the complete solution is: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; % ========================================================================= % SIMULATION % ========================================================================= % create the computational grid PML_size = 20; % size of the PML in grid points Nx = 128 - 2*PML_size; % number of grid points in the x (row) direction Ny = 256 - 2*PML_size; % 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] kgrid = makeGrid(Nx, dx, Ny, dy); % define the properties of the propagation medium medium.sound_speed = 1500; % [m/s] % create initial pressure distribution using makeDisc disc_magnitude = 5; % [au] disc_x_pos = 60; % [grid points] disc_y_pos = 140; % [grid points] disc_radius = 5; % [grid points] disc_2 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius); disc_x_pos = 30; % [grid points] disc_y_pos = 110; % [grid points] disc_radius = 8; % [grid points] disc_1 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius); medium.sound_speed = 1500*ones(Nx, Ny); % [m/s] disc_x_pos=round(Nx*0.7); disc_y_pos=round(Ny/4); disc_radius=5; air_buble = makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius); air_buble_reverse = (air_buble-1)*-1; medium.sound_speed = medium.sound_speed .* air_buble_reverse; medium.sound_speed = medium.sound_speed + air_buble * 340 ; medium.density = 7800*ones(Nx, Ny); % [kg/m^3] medium.density = medium.density .* air_buble_reverse; medium.density = medium.density + air_buble * 12.25 ; % smooth the initial pressure distribution and restore the magnitude %p0 = smooth(kgrid, disc_1 + disc_2, true); p0 = smooth(kgrid, disc_1 + disc_2, true); % assign to the source structure source.p0 = p0; % define a binary line sensor sensor.mask = zeros(Nx, Ny); sensor.mask(1, :) = 1; % create the time array [kgrid.t_array dt] = makeTime(kgrid, medium.sound_speed); % set the input arguements: force the PML to be outside the computational % grid; switch off p0 smoothing within kspaceFirstOrder2D input_args = {'PMLInside', false, 'PMLSize', PML_size, 'Smooth', false, 'PlotPML', false}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % reset the initial pressure source.p0 = 0; % assign the time reversal data sensor.time_reversal_boundary_data = sensor_data; % run the time reversal reconstruction p0_recon = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % add first order compensation for only recording over a half plane p0_recon = p0_recon*2; % repeat the FFT reconstruction for comparison p_xy = kspaceLineRecon(sensor_data.', dy, dt, medium.sound_speed, 'PosCond', true, 'Interp', 'linear'); % define a second k-space grid using the dimensions of p_xy [Nx_recon, Ny_recon] = size(p_xy); kgrid_recon = makeGrid(Nx_recon, dt*medium.sound_speed, Ny_recon, dy); % resample p_xy to be the same size as source.p0 p_xy_rs = interp2(kgrid_recon.y, kgrid_recon.x - min(kgrid_recon.x(:)), p_xy, kgrid.y, kgrid.x - min(kgrid.x(:))); % ========================================================================= % VISUALISATION % ========================================================================= % plot the initial pressure and sensor distribution figure; imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, p0 + sensor.mask*disc_magnitude, [-disc_magnitude disc_magnitude]); colormap(getColorMap); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; colorbar; % plot the reconstructed initial pressure figure; imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, p0_recon, [-disc_magnitude disc_magnitude]); colormap(getColorMap); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; colorbar; % apply a positivity condition p0_recon(p0_recon < 0) = 0; % plot the reconstructed initial pressure with positivity condition figure; imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, p0_recon, [-disc_magnitude disc_magnitude]); colormap(getColorMap); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; colorbar; % plot a profile for comparison figure; plot(kgrid.y_vec*1e3, p0(disc_x_pos, :), 'k-', ... kgrid.y_vec*1e3, p_xy_rs(disc_x_pos, :), 'r--', ... kgrid.y_vec*1e3, p0_recon(disc_x_pos, :), 'b:'); xlabel('y-position [mm]'); ylabel('Pressure'); legend('Initial Pressure', 'FFT Reconstruction', 'Time Reversal'); axis tight; set(gca, 'YLim', [0 5.1]);