gtMathsGetInterpolationIndices: extended to 1D and 2D and simplified

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

zUtil_Maths/gtMathsGetInterpolationIndices.m

-function [inds8, ints8] = gtMathsGetInterpolationIndices(pos, ints, data_type, mode)
+function [inds, ints] = gtMathsGetInterpolationIndices(pos, pows, data_type)
+% FUNCTION [inds, ints] = gtMathsGetInterpolationIndices(pos, pows, data_type)
+%
 
     if (~exist('data_type', 'var') || isempty(data_type))
         data_type = 'single';
     end
 
-    if (~exist('mode', 'var') || isempty(mode))
-        mode = 'Nicola';
+    pos = cast(pos, data_type);
+    pows = cast(pows, data_type);
+
+    pos_l = floor(pos); % lower position indices
+    pos_h = pos_l + 1;  % higher position indices
+
+    mu_l = pos_h - pos;
+    mu_h = 1 - mu_l;
+
+    n_dims = size(pos, 2);
+    switch (n_dims)
+        case 3
+            % uvw indices of 8 nearest neighbours ((nx8) x 3)
+            inds = [ pos_l; ...
+                [ pos_l(:, [1, 2]), pos_h(:, 3) ]; ...
+                [ pos_l(:, 1),      pos_h(:, 2),     pos_l(:, 3) ]; ...
+                [ pos_l(:, 1),      pos_h(:, [2, 3]) ]; ...
+                [ pos_h(:, 1),      pos_l(:, [2, 3]) ]; ...
+                [ pos_h(:, 1),      pos_l(:, 2),     pos_h(:, 3) ]; ...
+                [ pos_h(:, [1, 2]), pos_l(:, 3) ]; ...
+                pos_h ];
+
+            ints = [ ...
+                mu_l(:, 1) .* mu_l(:, 2) .* mu_l(:, 3) .* pows;
+                mu_l(:, 1) .* mu_l(:, 2) .* mu_h(:, 3) .* pows;
+                mu_l(:, 1) .* mu_h(:, 2) .* mu_l(:, 3) .* pows;
+                mu_l(:, 1) .* mu_h(:, 2) .* mu_h(:, 3) .* pows;
+                mu_h(:, 1) .* mu_l(:, 2) .* mu_l(:, 3) .* pows;
+                mu_h(:, 1) .* mu_l(:, 2) .* mu_h(:, 3) .* pows;
+                mu_h(:, 1) .* mu_h(:, 2) .* mu_l(:, 3) .* pows;
+                mu_h(:, 1) .* mu_h(:, 2) .* mu_h(:, 3) .* pows ];
+        case 2
+            inds = [ pos_l; ...
+                [ pos_l(:, 1), pos_h(:, 2) ]; ...
+                [ pos_h(:, 1), pos_l(:, 2) ]; ...
+                pos_h ];
+
+            ints = [ ...
+                mu_l(:, 1) .* mu_l(:, 2) .* pows;
+                mu_l(:, 1) .* mu_h(:, 2) .* pows;
+                mu_h(:, 1) .* mu_l(:, 2) .* pows;
+                mu_h(:, 1) .* mu_h(:, 2) .* pows ];
+        case 1
+            inds = [ pos_l; pos_h ];
+
+            ints = [ mu_l .* pows; mu_h .* pows ];
+        otherwise
+            error('gtMathsGetInterpolationIndices:wrong_argument', ...
+                'Interpolation can only be done for 1D, 2D and 3D positions')
     end
-
-    pos_l = floor(pos); % lower pixel indices
-    pos_h = pos_l + 1;      % higher pixel indices
-
-    % uvw indices of 8 nearest neighbours ((nx8) x 3)
-    inds8 = [ pos_l; ...
-        [ pos_l(:, [1, 2]), pos_h(:, 3) ]; ...
-        [ pos_l(:, 1),      pos_h(:, 2),     pos_l(:, 3) ]; ...
-        [ pos_l(:, 1),      pos_h(:, [2, 3]) ]; ...
-        [ pos_h(:, 1),      pos_l(:, [2, 3]) ]; ...
-        [ pos_h(:, 1),      pos_l(:, 2),     pos_h(:, 3) ]; ...
-        [ pos_h(:, [1, 2]), pos_l(:, 3) ]; ...
-        pos_h ];
-
-    ints8 = zeros(numel(ints), 8, data_type);
-
-    % Output is identical, but Nicola and Peter like to do the same thing
-    % in a different way :)
-    switch (mode)
-        case 'Peter'
-            % This way of rounding is correct also for round values when
-            % mu=mod... is used above:
-            % Modulus of uvw position
-            mu = mod(pos, 1);
-            mu1 = mu(:, 1);
-            mu2 = mu(:, 2);
-            mu3 = mu(:, 3);
-
-            mu12  = mu1 .* mu2;
-            mu13  = mu1 .* mu3;
-            mu23  = mu2 .* mu3;
-            mu123 = mu12 .* mu3;
-
-            ints8(:, 1) = mu12 + mu13 + mu23 - mu1 - mu2 - mu3 + - mu123 + 1;
-            ints8(:, 2) = mu3 - mu13 - mu23 + mu123;
-            ints8(:, 3) = mu2 - mu12 - mu23 + mu123;
-            ints8(:, 4) = mu23 - mu123;
-            ints8(:, 5) = mu1 - mu12 - mu13 + mu123;
-            ints8(:, 6) = mu13 - mu123;
-            ints8(:, 7) = mu12 - mu123;
-            ints8(:, 8) = mu123;
-
-            ints8 = ints8 .* ints(:, [1 1 1 1 1 1 1 1]);
-        case 'Nicola'
-            mu_l = pos_h - pos;
-            mu_h = 1 - mu_l;
-
-            ints8(:, 1) = mu_l(:, 1) .* mu_l(:, 2) .* mu_l(:, 3) .* ints;
-            ints8(:, 2) = mu_l(:, 1) .* mu_l(:, 2) .* mu_h(:, 3) .* ints;
-            ints8(:, 3) = mu_l(:, 1) .* mu_h(:, 2) .* mu_l(:, 3) .* ints;
-            ints8(:, 4) = mu_l(:, 1) .* mu_h(:, 2) .* mu_h(:, 3) .* ints;
-            ints8(:, 5) = mu_h(:, 1) .* mu_l(:, 2) .* mu_l(:, 3) .* ints;
-            ints8(:, 6) = mu_h(:, 1) .* mu_l(:, 2) .* mu_h(:, 3) .* ints;
-            ints8(:, 7) = mu_h(:, 1) .* mu_h(:, 2) .* mu_l(:, 3) .* ints;
-            ints8(:, 8) = mu_h(:, 1) .* mu_h(:, 2) .* mu_h(:, 3) .* ints;
-    end
-    ints8 = reshape(ints8, [], 1);
 end