Gradient/divergence: Reduced their performance impact

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

5_reconstruction/gtAstra3DTV.m

@@ -219,9 +219,9 @@ function output = projection(algo_id, output_id, input_id, input)
 end
 
 function d = gradient(vol)
-    dx = cat(1, vol(2:end, :, :) - vol(1:end-1, :, :), -vol(end, :, :));
-    dy = cat(2, vol(:, 2:end, :) - vol(:, 1:end-1, :), -vol(:, end, :));
-    dz = cat(3, vol(:, :, 2:end) - vol(:, :, 1:end-1), -vol(:, :, end));
+    dx = diff(padarray(vol, [1 0 0], 'post'), 1, 1);
+    dy = diff(padarray(vol, [0 1 0], 'post'), 1, 2);
+    dz = diff(padarray(vol, [0 0 1], 'post'), 1, 3);
     d = cat(4, dx, dy, dz);
 end
 
@@ -229,22 +229,8 @@ function vol = divergence(d)
     dx = d(:, :, :, 1);
     dy = d(:, :, :, 2);
     dz = d(:, :, :, 3);
-    ddx = cat(1, dx(1, :, :), dx(2:end, :, :) - dx(1:end-1, :, :));
-    ddy = cat(2, dy(:, 1, :), dy(:, 2:end, :) - dy(:, 1:end-1, :));
-    ddz = cat(3, dz(:, :, 1), dz(:, :, 2:end) - dz(:, :, 1:end-1));
+    ddx = diff(padarray(dx, [1 0 0], 'pre'), 1, 1);
+    ddy = diff(padarray(dy, [0 1 0], 'pre'), 1, 2);
+    ddz = diff(padarray(dz, [0 0 1], 'pre'), 1, 3);
     vol = ddx + ddy + ddz;
 end
-
-% function d = gradient2D(vol)
-%     dx = cat(1, vol(1:end-1, :, :) - vol(2:end, :, :), vol(end, :, :));
-%     dy = cat(2, vol(:, 1:end-1, :) - vol(:, 2:end, :), vol(:, end, :));
-%     d = cat(4, dx, dy, dz);
-% end
-% 
-% function vol = divergence2D(d)
-%     dx = d(:, :, :, 1);
-%     dy = d(:, :, :, 2);
-%     ddx = cat(1, -dx(1, :, :), dx(1:end-1, :, :) - dx(2:end, :, :));
-%     ddy = cat(2, -dy(:, 1, :), dy(:, 1:end-1, :) - dy(:, 2:end, :));
-%     vol = ddx + ddy;
-% end

zUtil_Maths/gtMathsDivergence.m

-function dd = gtMathsDivergence(d)
-% FUNCTION d = gtMathsDivergence(d)
-% Computes the divergence of a 1D-2D-3D object.
+function d = gtMathsDivergence(x)
+% FUNCTION d = gtMathsDivergence(x)
+% Computes the divergence of an n-dimensional vector object.
 
-    if (numel(d) >= 3)
-        dz = d{3};
-        dd{3} = cat(3, dz(:, :, 1), dz(:, :, 2:end) - dz(:, :, 1:end-1));
-    end
-    if (numel(d) >= 2)
-        dy = d{2};
-        dd{2} = cat(2, dy(:, 1, :), dy(:, 2:end, :) - dy(:, 1:end-1, :));
+    num_dims = numel(x);
+    d = cell(num_dims, 1);
+    for ii_d = 1:num_dims
+        padsize = zeros(1, num_dims);
+        padsize(ii_d) = 1;
+
+        d{ii_d} = diff(padarray(x{ii_d}, padsize, 'pre'), 1, ii_d);
     end
-    dx = d{1};
-    dd{1} = cat(1, dx(1, :, :), dx(2:end, :, :) - dx(1:end-1, :, :));
 
-    dd = gtMathsSumCellVolumes(dd);
+    d = gtMathsSumCellVolumes(d);
 end
\ No newline at end of file

zUtil_Maths/gtMathsGradient.m

 function d = gtMathsGradient(x)
-% FUNCTION d = gtMathsGradient(vol)
-% Computes the gradient of a 1D-2D-3D object and outputs a cell
+% FUNCTION d = gtMathsGradient(x)
+% Computes the gradient of an n-dimensional object and outputs a cell
 % concatenation of the directional gradients.
 
-    if (size(x, 3) > 1)
-        d{3} = cat(3, x(:, :, 2:end) - x(:, :, 1:end-1), -x(:, :, end));
-    end
-    if (size(x, 2) > 1)
-        d{2} = cat(2, x(:, 2:end, :) - x(:, 1:end-1, :), -x(:, end, :));
+    dims = size(x);
+    num_dims = numel(dims);
+    d = cell(num_dims, 1);
+    for ii_d = 1:num_dims
+        padsize = zeros(1, num_dims);
+        padsize(ii_d) = 1;
+
+        d{ii_d} = diff(padarray(x, padsize, 'post'), 1, ii_d);
     end
-    d{1} = cat(1, x(2:end, :, :) - x(1:end-1, :, :), -x(end, :, :));
 end
\ No newline at end of file