From 3a4f994541309cbdb8398907ed853d28e19afe01 Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Wed, 20 Apr 2016 15:32:44 +0200
Subject: [PATCH] gtMathsGetInterpolationIndices: extended to 1D and 2D and
simplified
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
zUtil_Maths/gtMathsGetInterpolationIndices.m | 112 +++++++++----------
1 file changed, 52 insertions(+), 60 deletions(-)
diff --git a/zUtil_Maths/gtMathsGetInterpolationIndices.m b/zUtil_Maths/gtMathsGetInterpolationIndices.m
index 11d01f6f..24398607 100644
--- a/zUtil_Maths/gtMathsGetInterpolationIndices.m
+++ b/zUtil_Maths/gtMathsGetInterpolationIndices.m
@@ -1,67 +1,59 @@
-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
--
GitLab