gtPlaceSubVolumes: added function to place many volumes in one call

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

zUtil_Imaging/gtGetVolsIntersectLimits.m

@@ -3,15 +3,18 @@ function [lims1, lims2] = gtGetVolsIntersectLimits(size1, size2, shift)
 % according to the shift.
 %   [lims1, lims2] = gtGetVolsIntersectLimits(size1, size2, shift)
 
-    size1((numel(size1)+1):3) = 1;
-    size2((numel(size2)+1):3) = 1;
-    if (numel(size1) > 3)
+    num_dims_1 = size(size1, 2);
+    num_dims_2 = size(size2, 2);
+
+    size1(:, (num_dims_1+1):3) = 1;
+    size2(:, (num_dims_2+1):3) = 1;
+    if (num_dims_1 > 3)
         warning('gtGetVolsIntersectLimits:wrong_argument', ...
-            'first volume size refers to a %d-dimensional volume', numel(size1))
+            'first volume size refers to a %d-dimensional volume', num_dims_1)
     end
-    if (numel(size2) > 3)
+    if (num_dims_2 > 3)
         warning('gtGetVolsIntersectLimits:wrong_argument', ...
-            'second volume size refers to a %d-dimensional volume', numel(size2))
+            'second volume size refers to a %d-dimensional volume', num_dims_2)
     end
 
     % volume1 limits
@@ -22,6 +25,10 @@ function [lims1, lims2] = gtGetVolsIntersectLimits(size1, size2, shift)
     lims2_min = max(1, -shift +1);
     lims2_max = min(size2, size1 - shift);
 
-    lims1 = [reshape(lims1_min, 1, []); reshape(lims1_max, 1, [])];
-    lims2 = [reshape(lims2_min, 1, []); reshape(lims2_max, 1, [])];
+    lims1 = [ ...
+        reshape(lims1_min', 1, num_dims_1, []); ...
+        reshape(lims1_max', 1, num_dims_1, []) ];
+    lims2 = [ ...
+        reshape(lims2_min', 1, num_dims_2, []); ...
+        reshape(lims2_max', 1, num_dims_2, []) ];
 end

zUtil_Imaging/gtPlaceSubVolumes.m

+function output = gtPlaceSubVolumes(output, input, shift, assign_op, use_c_functions)
+% GTPLACESUBVOLUME  Analogous to gtPlaceSubImage
+%     output = gtPlaceSubVolume(output, input, shift, index, assign_op, use_c_functions)
+%     ----------------------------------------------------------------------------------
+%     places the input volume in the output vol, with the origin of the
+%     input volume at point defined by shift = [x y z] in the output
+%     volume.
+%     Any part of the input volume falling outside of the output volume is
+%     cropped.
+%     Origin can contain negative coordinates, again anything outside the
+%     output volume is cropped.
+%     So... shift = [0 0 0] means that the voxel (a, b, c) in input will be
+%     at point (a, b, c) in the output.
+%     shift=[2 -2 0] means that voxel (a, b, c) will be at (a+2, b-2, c) in
+%     output.
+%
+%     if the output volume has more than 3 dimensions, the shifts have to
+%     be extended in order to select the position in the next dimensions.
+%     Otherwhise ones will be placed
+%     However, the placing will still happen in the first 3 dimensions
+%
+%     Note - this is used to add on volume to another one, rather than
+%     placing one image in the other.
+%     I have modified so it no longer adds it. I hope this doesn't break
+%     anything.
+%
+%     Note 2 - If index is 0, the volume will be copied as it is, otherwise
+%     all the voxels will have value of the index
+
+    if (~exist('use_c_functions', 'var') || isempty(use_c_functions))
+        use_c_functions = true;
+    end
+    if (~exist('assign_op', 'var') || isempty(assign_op))
+        assign_op = 'sum';
+    end
+
+    num_input = numel(input);
+
+    for ii = num_input:-1:1
+        try
+            input_size(ii, :) = size(input{ii});
+            input_type{ii} = class(input{ii});
+        catch mexc
+            warning('gtPlaceSubVolumes:wrong_argument', ...
+                'All inputs should have the same size')
+            rethrow(mexc)
+        end
+    end
+
+    outputSize = size(output);
+    num_input_dims = size(input_size, 2);
+    num_output_dims = size(outputSize, 2);
+    if (num_input_dims > 3)
+        error('gtPlaceSubVolume:wrong_argument', ...
+            'the subvolume to place should have dimensionality only up to 3D')
+    end
+    if (num_input_dims > num_output_dims)
+        warning('gtPlaceSubVolume:wrong_argument', ...
+            'the subvolume has a bigger dimensionality than the output volume')
+    end
+
+    % Force third dimension to be explicitly given
+    input_size(:, (num_input_dims+1):3) = 1;
+    outputSize((num_output_dims+1):3) = 1;
+    outputSize = outputSize(ones(num_input, 1), :);
+
+    num_shift_dims = size(shift, 2);
+    shift(:, (num_shift_dims+1):num_output_dims) = 0;
+
+    % output and input volume limits
+    [outLims, inLims] = gtGetVolsIntersectLimits(outputSize(:, 1:3), input_size, shift(:, 1:3));
+
+    if (use_c_functions)
+        shifts_op = [permute(outLims(1, :, :) - 1, [3 2 1]), shift(:, 4:end)];
+        shifts_ip = permute(inLims(1, :, :) - 1, [3 2 1]);
+        dims = permute(outLims(2, :, :) - outLims(1, :, :) + 1, [3 2 1]);
+
+        % Logicals are not handled in the C++ function
+        is_out_logical = islogical(output);
+        if (is_out_logical)
+            output = uint8(output);
+        end
+        % This will take care of 'logical' input as well
+        output_type = class(output);
+
+
+        wrong_type = ~strcmpi(output_type, input_type);
+        for ii = reshape(find(wrong_type), 1, [])
+            warning('gtPlaceSubVolume:heterogeneous_types', ...
+                'Converting input (%s) to output type (%s), it may reduce performance', ...
+                input_type{ii}, output_type)
+            input{ii} = cast(input{ii}, output_type);
+        end
+
+        switch (assign_op)
+            case 'sum'
+                for ii = 1:num_input
+                    output = gtCxxPlaceSubVolumeSum(output, input{ii}, shifts_op(ii, :), shifts_ip(ii, :), dims(ii, :));
+                end
+%                 output = gtCxxPlaceSubVolumeSum(output, input, shifts_op, shifts_ip, dims);
+            case 'assign'
+                for ii = 1:num_input
+                    output = gtCxxPlaceSubVolumeAssign(output, input{ii}, shifts_op(ii, :), shifts_ip(ii, :), dims(ii, :));
+                end
+%                 output = gtCxxPlaceSubVolumeAssign(output, input, shifts_op, shifts_ip, dims);
+            otherwise
+                    error('PLACE:wrong_argument', 'No option for "%s"', assign_op);
+        end
+
+        if (is_out_logical)
+            output = logical(output);
+        end
+    else
+        for ii = 1:num_input
+            input{ii} = input{ii}( ...
+                inLims(1, 1, ii):inLims(2, 1, ii), ...
+                inLims(1, 2, ii):inLims(2, 2, ii), ...
+                inLims(1, 3, ii):inLims(2, 3, ii) );
+            lims = [ ...
+                outLims(1, 1, ii), outLims(2, 1, ii), ...
+                outLims(1, 2, ii), outLims(2, 2, ii), ...
+                outLims(1, 3, ii), outLims(2, 3, ii) ];
+
+            % We add one because the shifting is done in matlab code!
+            extra_dims_shift = arrayfun(@(x){x+1}, shift(ii, 4:num_output_dims));
+
+            full_output = output;
+            if (~isempty(extra_dims_shift))
+                output = full_output(:, :, :, extra_dims_shift{:});
+            end
+
+            switch (assign_op)
+                case 'sum'
+                    % Matlab sugar to do the same that 'internal_gtAssignGrainToVol_sum'
+                    % does: (please keep them around, in case you break the C function)
+                    output(lims(1):lims(2), lims(3):lims(4), lims(5):lims(6)) ...
+                            = output(lims(1):lims(2), lims(3):lims(4), lims(5):lims(6)) ...
+                                + input;
+                case 'assign'
+                    % Matlab sugar to do the same that 'internal_gtAssignGrainToVol'
+                    % does: (please keep them around, in case you break the C function)
+                    output(lims(1):lims(2), lims(3):lims(4), lims(5):lims(6)) = input;
+
+                otherwise
+                    error('PLACE:wrong_argument', 'No option for "%s"', assign_op);
+            end
+
+            if (~isempty(extra_dims_shift))
+                full_output(:, :, :, extra_dims_shift{:}) = output;
+                output = full_output;
+            end
+        end
+    end
+end
+