diff --git a/5_reconstruction/GtAssembleVol3D.m b/5_reconstruction/GtAssembleVol3D.m
index 339c42aba91b145a1093fab37782307b03130ba9..c39c09d5576c82230941535aa2b530b8bce731e6 100644
--- a/5_reconstruction/GtAssembleVol3D.m
+++ b/5_reconstruction/GtAssembleVol3D.m
@@ -547,9 +547,10 @@ classdef GtAssembleVol3D < handle
                         single_dim_dmvol = gtCrop(single_dim_dmvol, segbb);
                         single_dim_dmvol(~seg_vol) = 0;
 
-                        dmvol(:, :, :, ii_dim) = gtPlaceSubVolume( ...
-                            dmvol(:, :, :, ii_dim), single_dim_dmvol, ...
-                            volBBox(1:3), [], obj.localPar.overlaps, false);
+                        % Turning the C functions off is way slower, but
+                        % for some reason they don't work in this case
+                        dmvol = gtPlaceSubVolume( dmvol, single_dim_dmvol, ...
+                            [volBBox(1:3), ii_dim-1], [], obj.localPar.overlaps, false);
                     end
                 catch mexc
                     gtPrintException(mexc, sprintf('\nGrain %d failed!!\n', ii));
diff --git a/zUtil_Imaging/gtPlaceSubVolume.m b/zUtil_Imaging/gtPlaceSubVolume.m
index 014821f408ddd0ef232dfc5b48bbd8810b858bf0..9008728f6b16f242570f57830a6062cc91a5dde3 100644
--- a/zUtil_Imaging/gtPlaceSubVolume.m
+++ b/zUtil_Imaging/gtPlaceSubVolume.m
@@ -32,12 +32,35 @@ function output = gtPlaceSubVolume(output, input, shift, index, assign_op, use_c
         assign_op = 'assign';
     end
 
+    inputSize = size(input);
+    outputSize = size(output);
+    num_input_dims = numel(inputSize);
+    num_output_dims = numel(outputSize);
+    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
-    [inputSize(1), inputSize(2), inputSize(3)] = size(input);
-    [outputSize(1), outputSize(2), outputSize(3)] = size(output);
+    inputSize((num_input_dims+1):3) = 1;
+    outputSize((num_output_dims+1):3) = 1;
+
+    shift((numel(shift)+1):max(num_output_dims, 3)) = 0;
+
+    % We add one because the shifting is done in matlab code!
+    extra_dims_shift = arrayfun(@(x){x+1}, shift(4:num_output_dims));
+
+    full_output = output;
+    if (~isempty(extra_dims_shift))
+        output = full_output(:, :, :, extra_dims_shift{:});
+    end
 
     % output and input volume limits
-    [outLims, inLims] = gtGetVolsIntersectLimits(outputSize, inputSize, shift);
+    [outLims, inLims] = gtGetVolsIntersectLimits(outputSize(1:3), inputSize, shift(1:3));
 
     input = input(  inLims(1, 1):inLims(2, 1), ...
                     inLims(1, 2):inLims(2, 2), ...
@@ -56,6 +79,15 @@ function output = gtPlaceSubVolume(output, input, shift, index, assign_op, use_c
             input = uint8(input);
         end
 
+        type_input = class(input);
+        type_output = class(output);
+        if (type_output ~= type_input)
+            warning('gtPlaceSubVolume:heterogeneous_types', ...
+                'Converting input (%s) to output type (%s), it may reduce performance', ...
+                type_input, type_output)
+            input = cast(input, type_output);
+        end
+
         switch (assign_op)
             case {'zero', 'conflict', 'adaptive'} % Mark overlapping...
                 output = internal_gtAssignGrainToVol_interf( output, input, index, lims );
@@ -97,4 +129,9 @@ function output = gtPlaceSubVolume(output, input, shift, index, assign_op, use_c
                 error('PLACE:wrong_argument', 'No option for "%s"', assign_op);
         end
     end
+
+    if (~isempty(extra_dims_shift))
+        full_output(:, :, :, extra_dims_shift{:}) = output;
+        output = full_output;
+    end
 end