-
Notifications
You must be signed in to change notification settings - Fork 0
/
read_ctf_mri.m
134 lines (115 loc) · 6.38 KB
/
read_ctf_mri.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
function [mri, hdr] = read_ctf_mri(filename);
% READ_CTF_MRI reads header and image data from CTF format MRI file
%
% [mri, hdr] = read_ctf_mri(filename)
%
% See also READ_CTF_MEG4, READ_CTF_RES4
% Copyright (C) 2003 Robert Oostenveld
%
fid = fopen(filename,'rb', 'ieee-be');
if fid<=0
error(sprintf('could not open MRI file: %s\n', filename));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READ THE IMAGE HEADER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
warning off
% general header information
hdr.identifierString = setstr(fread(fid,32,'char'))'; % CTF_MRI_FORMAT VER 2.2
hdr.imageSize = fread(fid,1,'int16'); % always = 256
hdr.dataSize = fread(fid,1,'int16'); % 1 or 2(bytes)
hdr.clippingRange = fread(fid,1,'int16'); % max.integer value of data
hdr.imageOrientation = fread(fid,1,'int16'); % eg., 0 = left on left, 1 = left on right
hdr.mmPerPixel_sagittal = fread(fid,1,'float'); % voxel dimensions in mm
hdr.mmPerPixel_coronal = fread(fid,1,'float'); % voxel dimensions in mm
hdr.mmPerPixel_axial = fread(fid,1,'float'); % voxel dimensions in mm
% HeadModel_Info specific header items
hdr.HeadModel.Nasion_Sag = fread(fid,1,'int16'); % fid.point coordinate(in voxels) for nasion - sagittal
hdr.HeadModel.Nasion_Cor = fread(fid,1,'int16'); % nasion - coronal
hdr.HeadModel.Nasion_Axi = fread(fid,1,'int16'); % nasion - axial
hdr.HeadModel.LeftEar_Sag = fread(fid,1,'int16'); % left ear - sagittal
hdr.HeadModel.LeftEar_Cor = fread(fid,1,'int16'); % left ear - coronal
hdr.HeadModel.LeftEar_Axi = fread(fid,1,'int16'); % left ear - axial
hdr.HeadModel.RightEar_Sag = fread(fid,1,'int16'); % right ear - sagittal
hdr.HeadModel.RightEar_Cor = fread(fid,1,'int16'); % right ear - coronal
hdr.HeadModel.RightEar_Axi = fread(fid,1,'int16'); % right ear - axial
fread(fid,2,'char'); % padding to 4 byte boundary
hdr.HeadModel.defaultSphereX = fread(fid,1,'float'); % sphere origin x coordinate(in mm)
hdr.HeadModel.defaultSphereY = fread(fid,1,'float'); % sphere origin y coordinate(in mm)
hdr.HeadModel.defaultSphereZ = fread(fid,1,'float'); % sphere origin z coordinate(in mm)
hdr.HeadModel.defaultSphereRadius = fread(fid,1,'float'); % default sphere radius(in mm)
% Image_Info specific header items
hdr.Image.modality = fread(fid,1,'int16'); % 0 = MRI, 1 = CT, 2 = PET, 3 = SPECT, 4 = OTHER
hdr.Image.manufacturerName = setstr(fread(fid,64,'char'))';
hdr.Image.instituteName = setstr(fread(fid,64,'char'))';
hdr.Image.patientID = setstr(fread(fid,32,'char'))';
hdr.Image.dateAndTime = setstr(fread(fid,32,'char'))';
hdr.Image.scanType = setstr(fread(fid,32,'char'))';
hdr.Image.contrastAgent = setstr(fread(fid,32,'char'))';
hdr.Image.imagedNucleus = setstr(fread(fid,32,'char'))';
fread(fid,2,'char'); % padding to 4 byte boundary
hdr.Image.Frequency = fread(fid,1,'float');
hdr.Image.FieldStrength = fread(fid,1,'float');
hdr.Image.EchoTime = fread(fid,1,'float');
hdr.Image.RepetitionTime = fread(fid,1,'float');
hdr.Image.InversionTime = fread(fid,1,'float');
hdr.Image.FlipAngle = fread(fid,1,'float');
hdr.Image.NoExcitations = fread(fid,1,'int16');
hdr.Image.NoAcquisitions = fread(fid,1,'int16');
hdr.Image.commentString = setstr(fread(fid,256,'char'))';
hdr.Image.forFutureUse = setstr(fread(fid,64,'char'))';
% continuation general header
hdr.headOrigin_sagittal = fread(fid,1,'float'); % voxel location of head origin
hdr.headOrigin_coronal = fread(fid,1,'float'); % voxel location of head origin
hdr.headOrigin_axial = fread(fid,1,'float'); % voxel location of head origin
% euler angles to align MR to head coordinate system(angles in degrees !)
hdr.rotate_coronal = fread(fid,1,'float'); % 1. rotate in coronal plane by this angle
hdr.rotate_sagittal = fread(fid,1,'float'); % 2. rotate in sagittal plane by this angle
hdr.rotate_axial = fread(fid,1,'float'); % 3. rotate in axial plane by this angle
hdr.orthogonalFlag = fread(fid,1,'int16'); % if set then image is orthogonal
hdr.interpolatedFlag = fread(fid,1,'int16'); % if set than image was interpolated
hdr.originalSliceThickness = fread(fid,1,'float'); % original spacing between slices before interpolation
transformMatrix = fread(fid,[4 4],'float')', % transformation matrix head->MRI[column][row]
fread(fid,204,'char'); % unused, padding to 1028 bytes
warning on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READ THE IMAGE DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if hdr.dataSize==1
mri = uint8(fread(fid, 256*256*256, 'uint8'));
elseif hdr.dataSize==2
mri = uint16(fread(fid, 256*256*256, 'uint16'));
else
error('unknown datasize in CTF mri file');
end
mri = reshape(mri, [256 256 256]);
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DO POST-PROCESSING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reorient the image data to obtain corresponding image data and transformation matrix
mri = permute(mri, [3 1 2]); % this was determined by trial and error
% reorient the image data and the transformation matrix along the left-right direction
% remember that the fiducials in voxel coordinates also have to be flipped (see down)
mri = flipdim(mri, 1);
flip = [-1 0 0 256
0 1 0 0
0 0 1 0
0 0 0 1 ];
transformMatrix = flip*transformMatrix;
% re-compute the homogeneous transformation matrices (apply voxel scaling)
scale = eye(4);
scale(1,1) = hdr.mmPerPixel_sagittal;
scale(2,2) = hdr.mmPerPixel_coronal;
scale(3,3) = hdr.mmPerPixel_axial;
hdr.transformHead2MRI = transformMatrix*inv(scale);
hdr.transformMRI2Head = scale*inv(transformMatrix);
% determine location of fiducials in MRI voxel coordinates
% flip the fiducials in voxel coordinates to correspond to the previous flip along left-right
hdr.fiducial.mri.nas = [256 - hdr.HeadModel.Nasion_Sag hdr.HeadModel.Nasion_Cor hdr.HeadModel.Nasion_Axi];
hdr.fiducial.mri.lpa = [256 - hdr.HeadModel.LeftEar_Sag hdr.HeadModel.LeftEar_Cor hdr.HeadModel.LeftEar_Axi];
hdr.fiducial.mri.rpa = [256 - hdr.HeadModel.RightEar_Sag hdr.HeadModel.RightEar_Cor hdr.HeadModel.RightEar_Axi];
% compute location of fiducials in MRI and HEAD coordinates
% hdr.fiducial.head.nas = warp3d(hdr.transformMRI2Head, hdr.fiducial.mri.nas, 'homogenous');
% hdr.fiducial.head.lpa = warp3d(hdr.transformMRI2Head, hdr.fiducial.mri.lpa, 'homogenous');
% hdr.fiducial.head.rpa = warp3d(hdr.transformMRI2Head, hdr.fiducial.mri.rpa, 'homogenous');