You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: coreg_siemens.m
    path of file: ./repos/osprey/coreg
    
    the code of the file until where you have to start completion: % coreg_siemens.m
% Georg Oeltzschner, Johns Hopkins University 2019.
% 
% USAGE:
% [vol_mask] = coreg_siemens(in, vol_image, maskFile);
% 
% DESCRIPTION:
% Creates a SPM volume containing a voxel mask with the same dimensions 
% as the SPM volume containing a structural image. The voxel mask will be
% created with the geometry information stored in the FID-A structure "in".
% This routine works for the Siemens TWIX, DICOM, and RDA data types, where
% the voxel orientation is defined by a normal vector and a rotation around
% this normal vector.
% 
% CREDITS:
% The routine f
    
function [vol_mask, T1_max, voxel_ctr] = coreg_siemens(in, vol_image, maskFile)

% Deactivate MATLAB warnings and load geometry parameters
warning('off','MATLAB:nearlySingularMatrix');
warning('off','MATLAB:qhullmx:InternalWarning');
geom = in.geometry;

%%% 1. PREPARE THE STRUCTURAL IMAGE
% Create SPM volume and read in the NIfTI file with the structural image.
[T1,XYZ]    = spm_read_vols(vol_image);
T1_max      = max(T1(:));
%Shift imaging voxel coordinates by half an imaging voxel so that the XYZ matrix
%tells us the x,y,z coordinates of the MIDDLE of that imaging voxel.
[~,voxdim] = spm_get_bbox(vol_image,'fv');
voxdim = abs(voxdim)';
halfpixshift = -voxdim(1:3)/2;
halfpixshift(3) = -halfpixshift(3);
XYZ = XYZ + repmat(halfpixshift, [1 size(XYZ,2)]);


%%% 2. GENERATE THE COORDINATES OF THE VOXEL CORNERS
% Extract voxel position and rotation parameters
NormSag         = geom.rot.NormSag;
NormCor         = geom.rot.NormCor;
NormTra         = geom.rot.NormTra;
VoI_InPlaneRot  = geom.rot.VoI_InPlaneRot;
% Correct voxel offsets by table position (if field exists)
if isfield(geom.pos,'TablePosTra')
    VoxOffs = [geom.pos.PosSag+geom.pos.TablePosSag geom.pos.PosCor+geom.pos.TablePosCor geom.pos.PosTra+geom.pos.TablePosTra];
else
    VoxOffs = [geom.pos.PosSag geom.pos.PosCor geom.pos.PosTra];
end