gtPlaceSubVolume.m

function output = gtPlaceSubVolume(output, input, shift, index, 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 = false;
    end
    if (~exist('index', 'var') || isempty(index))
        index = 0;
    end
    if (~exist('assign_op', 'var') || isempty(assign_op))
        assign_op = 'assign';
    end

    if (index ~= 0)
        input = index .* input;
    end

    if (all(shift == round(shift)))
        output = place_sub_volume(output, input, shift, assign_op, use_c_functions);
    else
        inds = shift ~= round(shift);
        pos_inds = find(inds);
        output = place_sub_volume_at_inds(pos_inds, 1, output, input, shift, assign_op, use_c_functions);
    end
end

function output = place_sub_volume_at_inds(pos_inds, coeff, output, input, shift, assign_op, use_c_functions)
    if (isempty(pos_inds))
        output = place_sub_volume(output, coeff * input, shift, assign_op, use_c_functions);
    else
        ind = pos_inds(1);
        l_s_ind = floor(shift(ind));
        u_s_ind = l_s_ind + 1;

        l_c_ind = u_s_ind - shift(ind);
        u_c_ind = 1 - l_c_ind;

        l_shift = shift;
        l_shift(ind) = l_s_ind;
        u_shift = shift;
        u_shift(ind) = u_s_ind;

        if (l_c_ind > 0)
            output = place_sub_volume_at_inds(pos_inds(2:end), l_c_ind * coeff, output, input, l_shift, assign_op, use_c_functions);
        end
        if (u_c_ind > 0)
            output = place_sub_volume_at_inds(pos_inds(2:end), u_c_ind * coeff, output, input, u_shift, assign_op, use_c_functions);
        end
    end
end

function output = place_sub_volume(output, input, shift, assign_op, use_c_functions)
    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((num_input_dims+1):3) = 1;
    outputSize((num_output_dims+1):3) = 1;

    shift((numel(shift)+1):num_output_dims) = 0;

    % output and input volume limits
    [outLims, inLims] = gtGetVolsIntersectLimits(outputSize(1:3), inputSize, shift(1:3));

    if (use_c_functions)
        shifts_op = [outLims(1, :) - 1, shift(4:end)];
        shifts_ip = inLims(1, :) - 1;
        dims = outLims(2, :) - outLims(1, :) + 1;

        % Logicals are not handled in the C++ function
        is_out_logical = islogical(output);
        if (is_out_logical)
            output = uint8(output);
        end
        if (islogical(input))
            input = uint8(input);
        end

        type_input = class(input);
        type_output = class(output);
        if (~strcmpi(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 = gtCxxPlaceSubVolumeInterf(output, input, shifts_op, shifts_ip, dims);
            case 'summed'
                output = gtCxxPlaceSubVolumeSum(output, input, shifts_op, shifts_ip, dims);
            case {'assign', 'parent'}
                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
        % please keep them around, in case you break the C function
        input = input( ...
            inLims(1, 1):inLims(2, 1), ...
            inLims(1, 2):inLims(2, 2), ...
            inLims(1, 3):inLims(2, 3) );
        lims = [ ...
            outLims(1, 1), outLims(2, 1), ...
            outLims(1, 2), outLims(2, 2), ...
            outLims(1, 3), outLims(2, 3) ];

        % 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

        temp_out = output(lims(1):lims(2), lims(3):lims(4), lims(5):lims(6));

        switch (assign_op)
            case {'zero', 'conflict', 'adaptive'}
                temp_out(temp_out & input) = -1;
                indexes = ((temp_out == 0) & input);
                temp_out(indexes) = input(indexes);
            case 'summed'
                temp_out = temp_out + input;
            case {'assign', 'parent'}
                indexes = input ~= 0;
                temp_out(indexes) = input(indexes);
            otherwise
                error('PLACE:wrong_argument', 'No option for "%s"', assign_op);
        end

        output(lims(1):lims(2), lims(3):lims(4), lims(5):lims(6)) = temp_out;

        if (~isempty(extra_dims_shift))
            full_output(:, :, :, extra_dims_shift{:}) = output;
            output = full_output;
        end
    end
end