ctf_write_mri.html

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

0086 <span class="comment">%    'ieee-be.l64' or 's' - IEEE floating point with big-endian byte</span>
0087 <span class="comment">%                            ordering and 64 bit long data type.</span>
0088 [fid,message] = fopen(file,<span class="string">'wb'</span>,<span class="string">'s'</span>);
0089 <span class="keyword">if</span> fid &lt; 0, error(<span class="string">'cannot open file for writing'</span>); <span class="keyword">end</span>
0090 
0091 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0092 <span class="comment">% write the file header</span>
0093 fprintf(<span class="string">'...writing header '</span>);
0094 <a href="#_sub1" class="code" title="subfunction Version_2_Header_write(fid,Version_2_Header),">Version_2_Header_write</a>(fid, mri.hdr);
0095 
0096 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0097 <span class="comment">% check header size, header should be 1028 bytes</span>
0098 header_bytes = ftell(fid);
0099 fprintf(<span class="string">'(wrote %d bytes)\n'</span>,header_bytes);
0100 <span class="keyword">if</span> header_bytes ~= 1028,
0101   msg = sprintf(<span class="string">'failed to write 1028 bytes into the header, wrote %d bytes'</span>,header_bytes);
0102   error(msg);
0103 <span class="keyword">end</span>
0104 
0105 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0106 <span class="comment">% seek beyond the header, to the beginning of the data matrix</span>
0107 <span class="comment">%fseek(fid,1028,'bof');</span>
0108 
0109 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0110 <span class="comment">% check if the data is 8 or 16 bits</span>
0111 <span class="keyword">switch</span> mri.hdr.dataSize,
0112   <span class="keyword">case</span> 1, <span class="comment">% we have 8 bit data</span>
0113     fprintf(<span class="string">'...writing 8 bit image data\n'</span>);
0114     precision = <span class="string">'uchar'</span>;
0115   <span class="keyword">case</span> 2, <span class="comment">% we have 16 bit data</span>
0116     fprintf(<span class="string">'...writing 16 bit image data\n'</span>);
0117     precision = <span class="string">'int16'</span>;
0118   <span class="keyword">otherwise</span>,
0119     msg = sprintf(<span class="string">'unknown mri.hdr.dataSize: %g'</span>,mri.hdr.dataSize);
0120     error(msg);
0121 <span class="keyword">end</span>
0122 
0123 
0124 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0125 <span class="comment">% now have the data array in a 3D matrix, we have to write it out in the</span>
0126 <span class="comment">% correct byte order</span>
0127 
0128 <span class="comment">% The CTF MRI File format used by MRIViewer consists of a binary file with a</span>
0129 <span class="comment">% 1,028 byte header. The MRI data can be in 8-bit (unsigned character) or 16-bit</span>
0130 <span class="comment">% (unsigned short integer) format and consists of 256 x 256 pixel slices, stored as</span>
0131 <span class="comment">% 256 contiguous sagittal slices from left to right (or right to left if head orientation</span>
0132 <span class="comment">% is &quot;left-on-right&quot;). Each slice is stored as individual pixels starting at the</span>
0133 <span class="comment">% top left corner and scanning downwards row by row. Therefore the coronal</span>
0134 <span class="comment">% position is fastest changing, axial position second fastest changing and sagittal</span>
0135 <span class="comment">% position slowest changing value in the file, always in the positive direction for</span>
0136 <span class="comment">% each axis (see section on Head Coordinate System for axis definitions). By</span>
0137 <span class="comment">% default CTF MRI files have the file extension &quot;.mri&quot;</span>
0138 
0139 <span class="comment">% MRIViewer uses these cardinal directions as axes in an internal coordinate system</span>
0140 <span class="comment">% where sagittal = X, coronal = Y and axial = Z forming an additional</span>
0141 <span class="comment">% right-handed coordinate system which is translated and rotated with respect to</span>
0142 <span class="comment">% the Head Coordinate System and has its origin at the upper left anterior corner</span>
0143 <span class="comment">% of the volume.</span>
0144 
0145 PixelDim = 256;
0146 RowDim   = 256;
0147 SliceDim = 256;
0148 
0149 <span class="comment">% imageOrientation, 0 = left on left, 1 = left on right</span>
0150 <span class="keyword">switch</span> mri.hdr.imageOrientation,
0151   <span class="keyword">case</span> 0,
0152     fprintf(<span class="string">'...sagittal slices are neurological orientation (left is on the left)\n'</span>);
0153     fprintf(<span class="string">'...+X left to right, +Y anterior to posterior, +Z superior to inferior\n'</span>);
0154   <span class="keyword">case</span> 1,
0155     fprintf(<span class="string">'...sagittal slices are radiological orientation (left is on the right)\n'</span>);
0156     fprintf(<span class="string">'...+X right to left, +Y anterior to posterior, +Z superior to inferior\n'</span>);
0157   <span class="keyword">otherwise</span>,
0158     msg = sprintf(<span class="string">'...unknown mri.hdr.imageOrientation: %d\n'</span>,mri.hdr.imageOrientation);
0159     error(msg);
0160 <span class="keyword">end</span>
0161 
0162 <span class="comment">% output into sagittal slices, with the fastest moving index being Y,</span>
0163 <span class="comment">% from anterior to posterior, then Z, from superior to inferior, then X,</span>
0164 <span class="comment">% from left to right (or vice versa; depending on input mri struct).</span>
0165 
0166 n = 1;
0167 y = 1:PixelDim; <span class="comment">% +Y is from anterior to posterior</span>
0168 
0169 <span class="keyword">for</span> x = 1:SliceDim, <span class="comment">% +X is from left to right (or vice versa)</span>
0170   <span class="keyword">for</span> z = 1:RowDim, <span class="comment">% +Z is from superior to inferior</span>
0171       
0172       count = fwrite(fid,mri.img(x,y,z),precision);
0173       
0174       <span class="keyword">if</span> count ~= PixelDim,
0175           error(<span class="string">'failed to output 256 data points'</span>);
0176       <span class="keyword">end</span>
0177       
0178   <span class="keyword">end</span>
0179 <span class="keyword">end</span>
0180 
0181 fclose(fid);
0182 
0183 t=toc; fprintf(<span class="string">'...done (%5.2f sec).\n\n'</span>,t);
0184 
0185 <span class="keyword">return</span>
0186 
0187 
0188 
0189 
0190 
0191 
0192 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0193 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0194 <a name="_sub1" href="#_subfunctions" class="code">function Version_2_Header_write(fid,Version_2_Header),</a>
0195 
0196 identifierString = sprintf(<span class="string">'%-32s'</span>, Version_2_Header.identifierString);
0197 <span class="keyword">if</span> length(identifierString) &lt; 32,
0198     paddingN = 32 - length(identifierString);
0199     padding = char(repmat(double(<span class="string">' '</span>),1,paddingN));
0200     identifierString = [identifierString,padding];
0201 <span class="keyword">end</span>
0202     
0203 fwrite(fid,identifierString(1:32),<span class="string">'char'</span>);
0204 
0205 fwrite(fid,Version_2_Header.imageSize           ,<span class="string">'short'</span>); <span class="comment">% always = 256</span>
0206 fwrite(fid,Version_2_Header.dataSize            ,<span class="string">'short'</span>); <span class="comment">% 1 or 2 (bytes), 8 or 16 bits</span>
0207 fwrite(fid,Version_2_Header.clippingRange       ,<span class="string">'short'</span>); <span class="comment">% max. integer value of data</span>
0208 fwrite(fid,Version_2_Header.imageOrientation    ,<span class="string">'short'</span>); <span class="comment">% eg., 0 = left on left, 1 = left on right</span>
0209 
0210 <span class="comment">% voxel dimensions in mm</span>
0211 fwrite(fid,Version_2_Header.mmPerPixel_sagittal ,<span class="string">'float'</span>);
0212 fwrite(fid,Version_2_Header.mmPerPixel_coronal  ,<span class="string">'float'</span>);
0213 fwrite(fid,Version_2_Header.mmPerPixel_axial    ,<span class="string">'float'</span>);
0214 
0215 <a href="#_sub2" class="code" title="subfunction headModel_write(fid,HeadModel_Info),">headModel_write</a>(fid,Version_2_Header.HeadModel_Info); <span class="comment">% defined below...</span>
0216 <a href="#_sub3" class="code" title="subfunction imageInfo_write(fid,Image_Info),">imageInfo_write</a>(fid,Version_2_Header.Image_Info);     <span class="comment">% defined below...</span>
0217 
0218 <span class="comment">% voxel location of head origin</span>
0219 fwrite(fid,Version_2_Header.headOrigin_sagittal ,<span class="string">'float'</span>);
0220 fwrite(fid,Version_2_Header.headOrigin_coronal  ,<span class="string">'float'</span>);
0221 fwrite(fid,Version_2_Header.headOrigin_axial    ,<span class="string">'float'</span>);
0222 
0223 <span class="comment">% euler angles to align MR to head coordinate system (angles in degrees!)</span>
0224 <span class="comment">% 1. rotate in coronal plane by this angle</span>
0225 <span class="comment">% 2. rotate in sagittal plane by this angle</span>
0226 <span class="comment">% 3. rotate in axial plane by this angle</span>
0227 fwrite(fid,Version_2_Header.rotate_coronal  ,<span class="string">'float'</span>);
0228 fwrite(fid,Version_2_Header.rotate_sagittal ,<span class="string">'float'</span>);
0229 fwrite(fid,Version_2_Header.rotate_axial    ,<span class="string">'float'</span>);
0230 
0231 fwrite(fid,Version_2_Header.orthogonalFlag   ,<span class="string">'short'</span>); <span class="comment">% if set then image is orthogonal</span>
0232 fwrite(fid,Version_2_Header.interpolatedFlag ,<span class="string">'short'</span>); <span class="comment">% if set than image was interpolated</span>
0233 
0234 <span class="comment">% original spacing between slices before interpolation to CTF format</span>
0235 fwrite(fid,Version_2_Header.originalSliceThickness ,<span class="string">'float'</span>);
0236 
0237 <span class="comment">% transformation matrix head-&gt;MRI [column][row]</span>
0238 fwrite(fid,Version_2_Header.transformMatrix' ,<span class="string">'float'</span>)';
0239 
0240 <span class="comment">% padding to 1028 bytes</span>
0241 <span class="comment">% according to CTF manual, this should</span>
0242 <span class="comment">%be 202, but it works out to the 1028 bytes with 204.</span>
0243 spaces = char(repmat(double(<span class="string">' '</span>),1,204));
0244 fwrite(fid,spaces,<span class="string">'uchar'</span>);