diff --git a/k-Wave/examples/example_pr_2D_TR_iterative.m b/k-Wave/examples/example_pr_2D_TR_iterative.m index b0a8fac..a7ae729 100644 --- a/k-Wave/examples/example_pr_2D_TR_iterative.m +++ b/k-Wave/examples/example_pr_2D_TR_iterative.m @@ -5,9 +5,9 @@ % example builds on the 2D Time Reversal Reconstruction For A Line Sensor % Example. % -% author: Ben Cox and Bradley Treeby +% author: Ben Cox, Bradley Treeby and Felix Lucka % date: 22nd January 2012 -% last update: 29th July 2019 +% last update: 22nd October 2021 % % This function is part of the k-Wave Toolbox (http://www.k-wave.org) % Copyright (C) 2012-2019 Ben Cox and Bradley Treeby @@ -32,7 +32,7 @@ % ========================================================================= % define literals -NUMBER_OF_ITERATIONS = 3; % number of iterations +NUMBER_OF_ITERATIONS = 5; % number of iterations PML_SIZE = 20; % size of the perfectly matched layer in grid points % create the computational grid @@ -68,59 +68,33 @@ % set the input arguments: 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, 'PlotSim', true}; + 'PlotPML', false, 'PlotSim', false}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % ========================================================================= -% RECONSTRUCT AN IMAGE USING TIME REVERSAL +% RECONSTRUCT AN IMAGE USING ITERATIVE TIME REVERSAL % ========================================================================= -% remove the initial pressure field used in the simulation -source = rmfield(source, 'p0'); +% set the initial reconstructed image to be zeros +p0_estimate = zeros(Nx, Ny); -% use the sensor points as sources in time reversal -source.p_mask = sensor.mask; +% set the initial model times series to be zero +modelled_time_series = zeros(size(sensor_data)); -% time reverse and assign the data -source.p = fliplr(sensor_data); +% calculate the difference between the measured and modelled data +data_difference = sensor_data - modelled_time_series; -% enforce, rather than add, the time-reversed pressure values -source.p_mode = 'dirichlet'; - -% set the simulation to record the final image (at t = 0) -sensor.record = {'p_final'}; - -% run the time reversal reconstruction -p0_estimate = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); - -% apply a positivity condition -p0_estimate.p_final = p0_estimate.p_final .* (p0_estimate.p_final > 0); - -% store the latest image estimate -p0_1 = p0_estimate.p_final; +% track goodness of fit and relative error during the iteration +gof = ones(NUMBER_OF_ITERATIONS+1, 1); +rel_err = ones(NUMBER_OF_ITERATIONS+1, 1); % ========================================================================= % ITERATE TO IMPROVE THE IMAGE % ========================================================================= -for loop = 2:NUMBER_OF_ITERATIONS - - % remove the source used in the previous time reversal - source = rmfield(source, 'p'); - - % set the initial pressure to be the latest estimate of p0 - source.p0 = p0_estimate.p_final; - - % set the simulation to record the time series - sensor = rmfield(sensor, 'record'); - - % calculate the time series using the latest estimate of p0 - sensor_data2 = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); - - % calculate the error in the estimated time series - data_difference = sensor_data - sensor_data2; +for loop = 1:NUMBER_OF_ITERATIONS % assign the data_difference as a time-reversal source source.p_mask = sensor.mask; @@ -133,16 +107,36 @@ % run the time reversal reconstruction p0_update = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); - - % add the update to the latest image - p0_estimate.p_final = p0_estimate.p_final + p0_update.p_final; - + + % update the image + p0_estimate = p0_estimate + p0_update.p_final; + % apply a positivity condition - p0_estimate.p_final = p0_estimate.p_final .* (p0_estimate.p_final > 0); + p0_estimate = p0_estimate .* (p0_estimate >= 0); % store the latest image estimate - eval(['p0_' num2str(loop) ' = p0_estimate.p_final;']); + p0_iterates{loop} = p0_estimate; + + + % remove the source used in the previous time reversal + source = rmfield(source, 'p'); + % set the initial pressure to be the latest estimate of p0 + source.p0 = p0_estimate; + + % set the simulation to record the time series + sensor = rmfield(sensor, 'record'); + + % calculate the time series using the latest estimate of p0 + modelled_time_series = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); + + % calculate the error in the estimated time series + data_difference = sensor_data - modelled_time_series; + + % measure goodness of fit and relative error + gof(loop+1) = norm(data_difference(:))^2/norm(sensor_data(:))^2; + rel_err(loop+1) = norm(p0_estimate(:) - p0(:))^2/norm(p0(:))^2; + end % ========================================================================= @@ -168,7 +162,7 @@ % plot the first iteration figure; -imagesc(p0_1, c_axis); +imagesc(p0_iterates{1}, c_axis); axis image; set(gca, 'XTick', [], 'YTick', []); ylabel('x-position [mm]'); @@ -184,7 +178,7 @@ % plot the 2nd iteration figure; -imagesc(p0_2, c_axis); +imagesc(p0_iterates{2}, c_axis); axis image; set(gca, 'XTick', [], 'YTick', []); title('Time Reversal Reconstruction, 2 Iterations'); @@ -196,12 +190,12 @@ plot([Ny, 1], [1, 1], 'k-'); plot([1, 1], [1, Nx], 'k-'); -% plot the 3rd iteration +% plot the last iteration figure; -imagesc(p0_3, c_axis); +imagesc(p0_iterates{end}, c_axis); axis image; set(gca, 'XTick', [], 'YTick', []); -title('Time Reversal Reconstruction, 3 Iterations'); +title(['Time Reversal Reconstruction, ' int2str(NUMBER_OF_ITERATIONS) ' Iterations']); colorbar; scaleFig(1, 0.65); @@ -209,4 +203,10 @@ hold on; plot([Ny, 1], [1, 1], 'k-'); plot([1, 1], [1, Nx], 'k-'); -plot([Ny, 1], [1, Nx], 'k--'); \ No newline at end of file +plot([Ny, 1], [1, Nx], 'k--'); + +% plot gof and relative error +figure(); +plot(0:NUMBER_OF_ITERATIONS, gof, 'DisplayName', 'goodness of fit'); hold on +plot(0:NUMBER_OF_ITERATIONS, rel_err, 'DisplayName', 'relative error'); hold on +legend(gca,'show'); diff --git a/k-Wave/examples/example_pr_2D_adjoint.m b/k-Wave/examples/example_pr_2D_adjoint.m index 98debb7..a9808c7 100644 --- a/k-Wave/examples/example_pr_2D_adjoint.m +++ b/k-Wave/examples/example_pr_2D_adjoint.m @@ -11,9 +11,9 @@ % For a more detailed discussion of this example and the underlying % techniques, see Arridge et al. Inverse Problems 32, 115012 (2016). % -% author: Ben Cox and Bradley Treeby +% author: Ben Cox, Bradley Treeby and Felix Lucka % date: 29th May 2017 -% last update: 29th July 2019 +% last update: 22nd October 2021 % % This function is part of the k-Wave Toolbox (http://www.k-wave.org) % Copyright (C) 2017-2019 Ben Cox and Bradley Treeby @@ -38,7 +38,7 @@ % ========================================================================= % define literals -NUMBER_OF_ITERATIONS = 5; % number of iterations +NUMBER_OF_ITERATIONS = 5; % number of iterations PML_SIZE = 20; % size of the perfectly matched layer in grid points % create the computational grid @@ -89,6 +89,13 @@ % set the initial model times series to be zero modelled_time_series = zeros(size(sensor_data)); +% calculate the difference between the measured and modelled data +data_difference = modelled_time_series - sensor_data; + +% track goodness of fit and relative error during the iteration +gof = ones(NUMBER_OF_ITERATIONS+1, 1); +rel_err = ones(NUMBER_OF_ITERATIONS+1, 1); + % reconstruct p0 image iteratively using the adjoint to estimate the % functional gradient for loop = 1:NUMBER_OF_ITERATIONS @@ -106,33 +113,35 @@ % set the source type to act as an adjoint source source.p_mode = 'additive'; - % calculate the difference between the measured and modelled data - difference = modelled_time_series - sensor_data; - % assign the difference time series as an adjoint source - % (see Appendix B in Arridge et al. Inverse Problems 32, 115012 (2016)) - time_reversed_data = fliplr(difference); - source.p = [time_reversed_data(:, 1), time_reversed_data(:, 1), time_reversed_data(:, 1:end-1)] + ... - [zeros(size(time_reversed_data(:, 1))), time_reversed_data(:, 1:end-1), 2 * time_reversed_data(:, end)]; + % (see Appendix B, Eqn B.2 in Arridge et al. Inverse Problems 32, 115012 (2016)) + p_adj = [data_difference(:, end:-1:1), zeros(size(data_difference, 1), 1)] ... + + [zeros(size(data_difference, 1), 1), data_difference(:, end:-1:1)]; + p_adj(:, end-1) = p_adj(:, end-1) + p_adj(:, end); + p_adj = p_adj(:, 1:end-1); + source.p = p_adj; % send difference through adjoint model image_update = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); + % add smoothing (see Appendix B in Arridge et al. Inverse Problems 32, 115012 (2016)) + image_update = smooth(image_update.p_final, false); + % set the step length (this could be chosen by doing a line search) - nu = 0.25; + nu = 0.5; % update the image - p0_estimate = p0_estimate - nu * image_update.p_final; + p0_estimate = p0_estimate - nu * image_update; % apply a positivity condition p0_estimate = p0_estimate .* (p0_estimate >= 0); % store the latest image estimate - eval(['p0_' num2str(loop) ' = p0_estimate;']); + p0_iterates{loop} = p0_estimate; % set the latest image to be the initial pressure in the forward model source = rmfield(source, 'p'); - source.p0 = p0_estimate; + source.p0 = smooth(p0_estimate, true); % set the sensor to record time series (by default) sensor = rmfield(sensor, 'record'); @@ -140,6 +149,13 @@ % calculate the time series at the sensors using the latest estimate of p0 modelled_time_series = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); + % calculate the difference between the measured and modelled data + data_difference = modelled_time_series - sensor_data; + + % measure goodness of fit and relative error + gof(loop+1) = norm(data_difference(:))^2/norm(sensor_data(:))^2; + rel_err(loop+1) = norm(p0_estimate(:) - p0(:))^2/norm(p0(:))^2; + end % ========================================================================= @@ -165,7 +181,7 @@ % plot the first iteration figure; -imagesc(p0_1, c_axis); +imagesc(p0_iterates{1}, c_axis); axis image; set(gca, 'XTick', [], 'YTick', []); title('Gradient Descent, 1 Iteration'); @@ -179,7 +195,7 @@ % plot the 2nd iteration figure; -imagesc(p0_2, c_axis); +imagesc(p0_iterates{2}, c_axis); axis image; set(gca, 'XTick', [], 'YTick', []); title('Gradient Descent, 2 Iterations'); @@ -191,12 +207,12 @@ plot([Ny, 1], [1, 1], 'k-'); plot([1, 1], [1, Nx], 'k-'); -% plot the 5th iteration +% plot the last iteration figure; -imagesc(p0_5, c_axis); +imagesc(p0_iterates{end}, c_axis); axis image; set(gca, 'XTick', [], 'YTick', []); -title('Gradient Descent, 5 Iterations'); +title(['Gradient Descent, ' int2str(NUMBER_OF_ITERATIONS) ' Iterations']); colorbar; scaleFig(1, 0.65); @@ -204,4 +220,10 @@ hold on; plot([Ny, 1], [1, 1], 'k-'); plot([1, 1], [1, Nx], 'k-'); -plot([Ny, 1], [1, Nx], 'k--'); \ No newline at end of file +plot([Ny, 1], [1, Nx], 'k--'); + +% plot gof and relative error +figure(); +plot(0:NUMBER_OF_ITERATIONS, gof, 'DisplayName', 'goodness of fit'); hold on +plot(0:NUMBER_OF_ITERATIONS, rel_err, 'DisplayName', 'relative error'); hold on +legend(gca,'show'); diff --git a/k-Wave/testing/regression/generateRegressionData.m b/k-Wave/testing/regression/generateRegressionData.m index fbbfeb6..5816548 100644 --- a/k-Wave/testing/regression/generateRegressionData.m +++ b/k-Wave/testing/regression/generateRegressionData.m @@ -243,7 +243,7 @@ index = index + 1; test_examples{index}.name = 'example_pr_2D_adjoint'; -test_examples{index}.outputs = {'sensor_data', 'p0', 'p0_1', 'p0_2', 'p0_3', 'p0_4', 'p0_5'}; +test_examples{index}.outputs = {'sensor_data', 'p0_estimate'}; test_examples{index}.precision = 'double'; % ========================================================================= diff --git a/k-Wave/testing/unit/kWaveDiffusion_compare_1D_2D_3D_plane_waves.m b/k-Wave/testing/unit/kWaveDiffusion_compare_1D_2D_3D_plane_waves.m index d261f10..1140082 100644 --- a/k-Wave/testing/unit/kWaveDiffusion_compare_1D_2D_3D_plane_waves.m +++ b/k-Wave/testing/unit/kWaveDiffusion_compare_1D_2D_3D_plane_waves.m @@ -108,7 +108,7 @@ comparison_thresh = COMPARISON_THRESH_DOUBLE; case 2 comparison_thresh = COMPARISON_THRESH_SINGLE; - input_args = [input_args, {'DataCast', 'gpuArray-single'}]; %#ok + input_args = [input_args, {'DataCast', 'single'}]; %#ok end % loop through tests diff --git a/k-Wave/testing/unit/kWaveDiffusion_compare_3D_heterog_perfusion_with_exact.m b/k-Wave/testing/unit/kWaveDiffusion_compare_3D_heterog_perfusion_with_exact.m index b68f624..9b1af85 100644 --- a/k-Wave/testing/unit/kWaveDiffusion_compare_3D_heterog_perfusion_with_exact.m +++ b/k-Wave/testing/unit/kWaveDiffusion_compare_3D_heterog_perfusion_with_exact.m @@ -105,7 +105,7 @@ case 1 input_args = {'PlotSim', plot_simulations, 'PlotScale', [37, 38]}; case 2 - input_args = {'PlotSim', plot_simulations, 'PlotScale', [37, 38], 'DataCast', 'gpuArray-single'}; + input_args = {'PlotSim', plot_simulations, 'PlotScale', [37, 38], 'DataCast', 'single'}; end % create kWaveDiffusion object and take time steps diff --git a/k-Wave/testing/unit/pstdElastic2D_save_movie.m b/k-Wave/testing/unit/pstdElastic2D_save_movie.m index fdc1eb8..b35a910 100644 --- a/k-Wave/testing/unit/pstdElastic2D_save_movie.m +++ b/k-Wave/testing/unit/pstdElastic2D_save_movie.m @@ -27,11 +27,7 @@ test_pass = true; % folder to store temporary movies in -if nargin == 0 - folder = ''; -else - folder = tempdir; -end +folder = tempdir; % movie names movie_name_1 = 'test-movie-elastic-2D-avi';