gtPlaceSubVolume.m

function output = gtPlaceSubVolume(output, input, shift, index, assign_op, use_c_functions)
% GTPLACESUBVOLUME  Analogous to gtPlaceSubImage
%     output = gtPlaceSubVolume(output, input, shift)
%     -----------------------------------------------
%     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.
%
%     Note - this used to add on volume to another, rather than place one image
%     in the other.  I have modified so it no longer adds.  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'))
        use_c_functions = true;
    end
    if (~exist('index', 'var'))
        index = 0;
    end
    if (~exist('assign_op', 'var'))
        assign_op = 'assign';
    end

    [inputSize(1), inputSize(2), inputSize(3)] = size(input);
    [outputSize(1), outputSize(2), outputSize(3)] = size(output);

    % output volume limits
    outMinLims = max(1, shift +1);
    outMaxLims = min(outputSize, shift + inputSize);

    % input volume limits
    inMinLims = max(1, -shift +1);
    inMaxLims = min(inputSize, outputSize - shift);

    input = input(  inMinLims(1):inMaxLims(1), ...
                    inMinLims(2):inMaxLims(2), ...
                    inMinLims(3):inMaxLims(3)  );
    lims = [  outMinLims(1) outMaxLims(1) ...
                outMinLims(2) outMaxLims(2) ...
                outMinLims(3) outMaxLims(3)  ];

    if (use_c_functions)
        switch (assign_op)
            case {'zero', 'conflict', 'adaptive'} % Mark overlapping...
                output = internal_gtAssignGrainToVol_interf( output, input, index, lims );
            case 'summed'
                output = internal_gtAssignGrainToVol_sum( output, input, index, lims );
            case 'assign'
                output = internal_gtAssignGrainToVol( output, input, index, lims );
            otherwise
                error('PLACE:wrong_argument', 'No option for "%s"', assign_op);
        end
    else
        switch (assign_op)
            case {'zero', 'conflict', 'adaptive'}
                % Matlab sugars to do the same that 'internal_gtAssignGrainToVol_interf'
                % does: (please keep them around, in case you break the C function)
                temp_out = output(lims(1):lims(2), lims(3):lims(4), lims(5):lims(6));

                temp_out(temp_out & input) = -1;
                indexes = ((temp_out ==0) & input);
                temp_out(indexes) = input(indexes);

                output(lims(1):lims(2), lims(3):lims(4), lims(5):lims(6)) = temp_out;
            case 'summed'
                % 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
    end
end