k-Wave User Forum » Topic: 2D image reconstruction using 'universal back-projection algorithm'

k-Wave User Forum » Topic: 2D image reconstruction using 'universal back-projection algorithm' http://www.k-wave.org/forum/topic/2d-image-reconstruction-using-universal-back-projection-algorithm Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 01:02:16 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php james luo on "2D image reconstruction using 'universal back-projection algorithm'" http://www.k-wave.org/forum/topic/2d-image-reconstruction-using-universal-back-projection-algorithm#post-8610 Sun, 11 Sep 2022 09:37:11 +0000 james luo 8610@http://www.k-wave.org/forum/ Dear Prof. Ben Cox and Prof. Bradley Treeby, I meet the same problem liked Tho and failed to use 2D universal back-projection image reconstruction successfully, could you please kindly share me how you did it？ My email: <a href="mailto:1336459011@qq.com">1336459011@qq.com</a> Thank you ! Best wishes, James. bencox on "2D image reconstruction using 'universal back-projection algorithm'" http://www.k-wave.org/forum/topic/2d-image-reconstruction-using-universal-back-projection-algorithm#post-7195 Wed, 22 Jan 2020 16:01:42 +0000 bencox 7195@http://www.k-wave.org/forum/ Hi Tho, Sorry for the delay in replying. There is a 2D version of the universal backprojection algorithm derived by Burgholzer et al in this paper: Burgholzer et al. "Temporal back-projection algorithms for photoacoustic tomography with integrating line detectors." Inverse Problems 23.6 (2007): S65. Hope that helps. Ben Tho on "2D image reconstruction using 'universal back-projection algorithm'" http://www.k-wave.org/forum/topic/2d-image-reconstruction-using-universal-back-projection-algorithm#post-7131 Mon, 18 Nov 2019 23:04:19 +0000 Tho 7131@http://www.k-wave.org/forum/ Dear Prof. Ben Cox and Prof. Bradley Treeby, Dear colleagues, I am doing 2D image reconstruction with point sensor ring array using 'universal back-projection algorithm (UBP). As you know, the UBP in Prof. Lihong Wang's article 'Universal back-projection algorithm for photoacoustic computed tomography' was obtained in 3D, so I did a simplification of his reconstruction equation from 3D to 2D (To be honest, I am not 100% sure if I did it correctly. But here on the forum, I cannot post equations and pics). I used k-wave tool box to generate the forward date, from which I did the image reconstruction. However, my UBP reconstruction has some wired fluctuations at the image transition parts compared to time reversal image reconstruction, which I really don't know the reason. Increasing the number of sensors and decreasing the grid size do not alleviate the above mention reconstruction fluctuations much. I am running out thoughts now! Below, I post my matlab code. Could anyone please have a look at my code to see where I may not do correctly? If anyone has done 2D universal back-projection image reconstruction successfully, could you please kindly share me how you did it？ In case you need further input, my email: <a href="mailto:isabelle.thomas19950321@gmail.com">isabelle.thomas19950321@gmail.com</a>. Thank you all! Best wishes, Isabelle %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MATLAB CODE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clearvars close all % load the initial pressure distribution as shepp-Logan phantom p0_magnitude = 1; grid_size_SL = 100; p = p0_magnitude.*phantom('Modified Shepp-Logan',grid_size_SL); tot_grid_size = grid_size_SL + 100; %% the extra 100 grid points are used for PML and placing transducers p0 = zeros(tot_grid_size); p0(51:tot_grid_size-50,51:tot_grid_size-50) = p; % assign the grid size and create the computational grid PML_size = 20; % size of the PML in grid points Nx = tot_grid_size; % number of grid points in the x direction Ny = tot_grid_size; % number of grid points in the y direction x= 10e-3; % total grid size [m] y = 10e-3; % total grid size [m] dx = x / Nx; % grid point spacing in the x direction [m] dy = y / Ny; % grid point spacing in the y direction [m] kgrid = kWaveGrid(Nx, dx, Ny, dy); % smooth the initial pressure distribution and restore the magnitude p0 = smooth(kgrid, p0, true); % assign to the source structure source.p0 = p0; % define the properties of the propagation medium sound_speed = 1500; medium.sound_speed = sound_speed; % [m/s] % define a centered Cartesian circular sensor sensor_radius = 4.5e-3; % [m] sensor_angle = 2*pi; % [rad] sensor_pos = [0, 0]; % [m] num_sensor_points = 100; cart_sensor_mask = makeCartCircle(sensor_radius, num_sensor_points, sensor_pos, sensor_angle); % assign to sensor structure sensor.mask = cart_sensor_mask; % create the time array kgrid.makeTime(medium.sound_speed); % set the input options input_args = {'Smooth', false, 'PMLInside', true, 'PlotPML', true}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); %% back projection t_array = kgrid.t_array; dt = kgrid.dt; fs = 1/dt; np = length(t_array); NFFT = 2^nextpow2(np); f_vec = fs*linspace(0,1,NFFT); w_vec = 2*pi*f_vec; k_vec = ifftshift(2*pi*f_vec/(sound_speed)); weighting_omega = 0; p0_recon = zeros(Nx, Ny); for ii = 1:num_sensor_points x_ii = cart_sensor_mask(1,ii); y_ii = cart_sensor_mask(2,ii); distance_x_square = (kgrid.x - x_ii).^2; distance_y_square = (kgrid.y - y_ii).^2; distance_xy = sqrt(distance_x_square + distance_y_square); distance_xy_index = round(distance_xy/(sound_speed*dt)); distance_xy_index(~distance_xy_index)=1; p_ii = sensor_data(ii,:); p_ii_derivative0 = ifft((1i*k_vec).*fft(p_ii,NFFT)); p_ii_derivative = real((p_ii_derivative0(1:np))); bp_ii = (p_ii - sound_speed.*t_array.*p_ii_derivative); p0_ii = bp_ii(distance_xy_index); %%% weighting angle weighting_arc_ii = (2*pi*sensor_radius)/num_sensor_points; weighting_phy_ii = 2*pi/num_sensor_points; weighting_omega_ii = weighting_arc_ii*(1./((distance_x_square + distance_y_square).*sensor_radius).*(-x_ii.*kgrid.x - y_ii.*kgrid.y + sensor_radius.^2)); weighting_omega_ii(isnan(weighting_omega_ii)) = weighting_phy_ii; weighting_omega_ii(isinf(weighting_omega_ii)) = weighting_phy_ii; p0_ii = weighting_omega_ii.*p0_ii; weighting_omega = weighting_omega + weighting_omega_ii; p0_recon = p0_recon + p0_ii; end p0_recon0 = p0_recon./(2*pi); figure(1) %% reconstructed image imagesc(kgrid.y_vec*1000,kgrid.x_vec*1000,p0_recon0,[-1,1]) %./max(p0_recon(:)) axis image colorbar ylabel('x(mm)'); xlabel('y(mm)'); slice_pos = 4.5e-3; %% comparison with original image amplitude along a horizontal line figure(2) plot(kgrid.y_vec,p0(round(slice_pos/kgrid.dx),:)) hold on plot(kgrid.y_vec,p0_recon0(round(slice_pos/kgrid.dx),:)) hold off axis tight xlabel('y(mm)'); ylabel('Reconstructed amplitude');