ctf_read_mri.html

mri_toolbox是一个工具用来MRI. 来自于SourceForge, 我上传这个软件,希望能结识对医疗软件感兴趣的兄弟.

0127 RowDim   = 256;
0128 SliceDim = 256;
0129 
0130 <span class="comment">% imageOrientation, 0 = left on left, 1 = left on right</span>
0131 <span class="keyword">switch</span> mri.hdr.imageOrientation,
0132   <span class="keyword">case</span> 0,
0133     fprintf(<span class="string">'...sagittal slices are neurological orientation (left is on the left)\n'</span>);
0134     fprintf(<span class="string">'...+X left to right, +Y anterior to posterior, +Z superior to inferior\n'</span>);
0135   <span class="keyword">case</span> 1,
0136     fprintf(<span class="string">'...sagittal slices are radiological orientation (left is on the right)\n'</span>);
0137     fprintf(<span class="string">'...+X right to left, +Y anterior to posterior, +Z superior to inferior\n'</span>);
0138   <span class="keyword">otherwise</span>,
0139     msg = sprintf(<span class="string">'...unknown mri.hdr.imageOrientation: %d\n'</span>,mri.hdr.imageOrientation);
0140     error(msg);
0141 <span class="keyword">end</span>
0142 
0143 mri.img = zeros(SliceDim,PixelDim,RowDim);
0144 
0145 n = 1;
0146 y = 1:PixelDim; <span class="comment">% +Y is from anterior to posterior</span>
0147 
0148 <span class="keyword">for</span> x = 1:SliceDim, <span class="comment">% +X is from left to right</span>
0149   <span class="keyword">for</span> z = 1:RowDim, <span class="comment">% +Z is from superior to inferior</span>
0150     mri.img(x,y,z) = data(n:n+(PixelDim-1));
0151     n = n + PixelDim;
0152   <span class="keyword">end</span>
0153 <span class="keyword">end</span>
0154 
0155 
0156 
0157 t=toc; fprintf(<span class="string">'...done (%5.2f sec).\n\n'</span>,t);
0158 
0159 <span class="keyword">return</span>
0160 
0161 
0162 
0163 
0164 
0165 
0166 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0167 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0168 <a name="_sub1" href="#_subfunctions" class="code">function Version_2_Header = Version_2_Header_read(fid),</a>
0169 
0170 tmp = fread(fid,[1,32],<span class="string">'char'</span>);
0171 tmp(find(tmp&lt;0)) = 0;
0172 Version_2_Header.identifierString = char( tmp );
0173 
0174 <span class="comment">% check the header format version is 2.2</span>
0175 <span class="keyword">if</span> strmatch(<span class="string">'CTF_MRI_FORMAT VER 2.2'</span>,Version_2_Header.identifierString),
0176   <span class="comment">% OK this function should read this</span>
0177 <span class="keyword">else</span>
0178   msg = sprintf(<span class="string">'this function is not designed to read this format:\n%s\n'</span>,Version_2_Header.identifierString);
0179   warning(msg);
0180 <span class="keyword">end</span>
0181 
0182 Version_2_Header.imageSize           = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% always = 256</span>
0183 Version_2_Header.dataSize            = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% 1 or 2 (bytes), 8 or 16 bits</span>
0184 Version_2_Header.clippingRange       = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% max. integer value of data</span>
0185 Version_2_Header.imageOrientation    = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% eg., 0 = left on left, 1 = left on right</span>
0186 
0187 <span class="comment">% voxel dimensions in mm</span>
0188 Version_2_Header.mmPerPixel_sagittal = fread(fid,1,<span class="string">'float'</span>);
0189 Version_2_Header.mmPerPixel_coronal  = fread(fid,1,<span class="string">'float'</span>);
0190 Version_2_Header.mmPerPixel_axial    = fread(fid,1,<span class="string">'float'</span>);
0191 
0192 Version_2_Header.HeadModel_Info = <a href="#_sub2" class="code" title="subfunction HeadModel_Info = headModel(fid),">headModel</a>(fid); <span class="comment">% defined below...</span>
0193 Version_2_Header.Image_Info = <a href="#_sub3" class="code" title="subfunction Image_Info = imageInfo(fid),">imageInfo</a>(fid);     <span class="comment">% defined below...</span>
0194 
0195 <span class="comment">% voxel location of head origin</span>
0196 Version_2_Header.headOrigin_sagittal = fread(fid,1,<span class="string">'float'</span>);
0197 Version_2_Header.headOrigin_coronal  = fread(fid,1,<span class="string">'float'</span>);
0198 Version_2_Header.headOrigin_axial    = fread(fid,1,<span class="string">'float'</span>);
0199 
0200 <span class="comment">% euler angles to align MR to head coordinate system (angles in degrees!)</span>
0201 <span class="comment">% 1. rotate in coronal plane by this angle</span>
0202 <span class="comment">% 2. rotate in sagittal plane by this angle</span>
0203 <span class="comment">% 3. rotate in axial plane by this angle</span>
0204 Version_2_Header.rotate_coronal  = fread(fid,1,<span class="string">'float'</span>);
0205 Version_2_Header.rotate_sagittal = fread(fid,1,<span class="string">'float'</span>);
0206 Version_2_Header.rotate_axial    = fread(fid,1,<span class="string">'float'</span>);
0207 
0208 Version_2_Header.orthogonalFlag   = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% if set then image is orthogonal</span>
0209 Version_2_Header.interpolatedFlag = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% if set than image was interpolated</span>
0210 
0211 <span class="comment">% original spacing between slices before interpolation to CTF format</span>
0212 Version_2_Header.originalSliceThickness = fread(fid,1,<span class="string">'float'</span>);
0213 
0214 <span class="comment">% transformation matrix head-&gt;MRI [column][row]</span>
0215 Version_2_Header.transformMatrix = fread(fid,[4 4],<span class="string">'float'</span>)';
0216 
0217 
0218 Version_2_Header.transformMatrixHead2MRI = Version_2_Header.transformMatrix;
0219 Version_2_Header.transformMatrixMRI2Head = inv(Version_2_Header.transformMatrix);
0220 
0221 <span class="comment">% padding to 1028 bytes</span>
0222 <span class="comment">%tmp = fread(fid,[1,202],'uchar'); % according to CTF manual, this should</span>
0223 <span class="comment">%be 202, but it reads out to the 1028 bytes with 204.  So maybe there are a</span>
0224 <span class="comment">%few bytes in the file read operations above that are missed?</span>
0225 tmp = fread(fid,[1,204],<span class="string">'uchar'</span>);
0226 tmp = zeros(size(tmp));
0227 Version_2_Header.unused = char( tmp ); <span class="comment">% uchar</span>
0228 
0229 <span class="keyword">return</span>
0230 
0231 
0232 
0233 
0234 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0235 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0236 <a name="_sub2" href="#_subfunctions" class="code">function HeadModel_Info = headModel(fid),</a>
0237 
0238 <span class="comment">% this function is called from Version_2_Header_read</span>
0239 
0240 <span class="comment">% fid. point coordinate (in voxels)</span>
0241 HeadModel_Info.Nasion_Sag = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% nasion - sagittal</span>
0242 HeadModel_Info.Nasion_Cor = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% nasion - coronal</span>
0243 HeadModel_Info.Nasion_Axi = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% nasion - axial</span>
0244 HeadModel_Info.LeftEar_Sag = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% left ear - sagittal</span>
0245 HeadModel_Info.LeftEar_Cor = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% left ear - coronal</span>
0246 HeadModel_Info.LeftEar_Axi = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% left ear - axial</span>
0247 HeadModel_Info.RightEar_Sag = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% right ear - sagittal</span>
0248 HeadModel_Info.RightEar_Cor = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% right ear - coronal</span>
0249 HeadModel_Info.RightEar_Axi = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% right ear - axial</span>
0250 
0251 fread(fid,2,<span class="string">'char'</span>); <span class="comment">% padding to 4 byte boundary - from Robert Oostenveld</span>
0252 
0253 <span class="comment">% default sphere origin</span>
0254 HeadModel_Info.defaultSphereX = fread(fid,1,<span class="string">'float'</span>); <span class="comment">% sphere origin x coordinate ( in mm )</span>
0255 HeadModel_Info.defaultSphereY = fread(fid,1,<span class="string">'float'</span>); <span class="comment">% sphere origin y coordinate ( in mm )</span>
0256 HeadModel_Info.defaultSphereZ = fread(fid,1,<span class="string">'float'</span>); <span class="comment">% sphere origin z coordinate ( in mm )</span>
0257 HeadModel_Info.defaultSphereRadius = fread(fid,1,<span class="string">'float'</span>); <span class="comment">% default sphere radius ( in mm )</span>
0258 
0259 <span class="keyword">return</span>
0260 
0261 
0262 
0263 
0264 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0265 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0266 <a name="_sub3" href="#_subfunctions" class="code">function Image_Info = imageInfo(fid),</a>
0267 
0268 <span class="comment">% this function is called from Version_2_Header_read</span>
0269 
0270 Image_Info.modality         = fread(fid,1,<span class="string">'short'</span>); <span class="comment">% 0 = MRI, 1 = CT, 2 = PET, 3 = SPECT, 4 = OTHER</span>
0271 Image_Info.manufacturerName = char( fread(fid,[1,64],<span class="string">'char'</span>) );
0272 Image_Info.instituteName    = char( fread(fid,[1,64],<span class="string">'char'</span>) );
0273 Image_Info.patientID        = char( fread(fid,[1,32],<span class="string">'char'</span>) );
0274 Image_Info.dateAndTime      = char( fread(fid,[1,32],<span class="string">'char'</span>) );
0275 Image_Info.scanType         = char( fread(fid,[1,32],<span class="string">'char'</span>) );
0276 Image_Info.contrastAgent    = char( fread(fid,[1,32],<span class="string">'char'</span>) );
0277 Image_Info.imagedNucleus    = char( fread(fid,[1,32],<span class="string">'char'</span>) );
0278 
0279 fread(fid,2,<span class="string">'char'</span>); <span class="comment">% padding to 4 byte boundary - from Robert Oostenveld</span>
0280 
0281 Image_Info.Frequency        = fread(fid,1,<span class="string">'float'</span>);
0282 Image_Info.FieldStrength    = fread(fid,1,<span class="string">'float'</span>);
0283 Image_Info.EchoTime         = fread(fid,1,<span class="string">'float'</span>);
0284 Image_Info.RepetitionTime   = fread(fid,1,<span class="string">'float'</span>);
0285 Image_Info.InversionTime    = fread(fid,1,<span class="string">'float'</span>);
0286 Image_Info.FlipAngle        = fread(fid,1,<span class="string">'float'</span>);
0287 Image_Info.NoExcitations    = fread(fid,1,<span class="string">'short'</span>);
0288 Image_Info.NoAcquisitions   = fread(fid,1,<span class="string">'short'</span>);
0289 
0290 tmp = fread(fid,[1,256],<span class="string">'char'</span>);
0291 tmp = zeros(size(tmp));
0292 Image_Info.commentString = char( tmp );
0293 
0294 tmp = fread(fid,[1,64],<span class="string">'char'</span>);
0295 tmp = zeros(size(tmp));
0296 Image_Info.forFutureUse  = char( tmp );
0297 
0298 <span class="keyword">return</span>
0299 
0300 
0301 
0302 
0303 
0304 
0305 
0306 
0307 <span class="comment">% The CTF MRI File format used by MRIViewer consists of a binary file with a</span>
0308 <span class="comment">% 1,028 byte header. The MRI data can be in 8-bit (unsigned character) or 16-bit</span>
0309 <span class="comment">% (unsigned short integer) format and consists of 256 x 256 pixel slices, stored as</span>
0310 <span class="comment">% 256 contiguous sagittal slices from left to right (or right to left if head orientation</span>
0311 <span class="comment">% is 搇eft-on-right