#33 updated adjoint and TR examples (#34)

* updated adjoint and TR examples (thank you Felix Lucka)
* Fixed bug in generateRegressionData.m due to change in adjoint example
* Fixed cast to single bug in unit test for diffusion model
* Fixed tempdir folder bug in pstdElastic movie saving unit test

Author: bencox

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,17 +190,23 @@
 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);
 
 % add lines indicating the sensor mask and 'visible region'
 hold on;
 plot([Ny, 1], [1, 1], 'k-');
 plot([1, 1], [1, Nx], 'k-');
-plot([Ny, 1], [1, Nx], 'k--');
+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');

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,40 +113,49 @@
     % 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');
 
     % 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,17 +207,23 @@
 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);
 
 % add lines indicating the sensor mask and 'visible region'
 hold on;
 plot([Ny, 1], [1, 1], 'k-');
 plot([1, 1], [1, Nx], 'k-');
-plot([Ny, 1], [1, Nx], 'k--');
+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');

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';
 
 % =========================================================================

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<AGROW> 
+            input_args = [input_args, {'DataCast', 'single'}]; %#ok<AGROW> 
     end
 
     % loop through tests

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

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';