avw_img_read.html

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

0092 
0093 <span class="keyword">return</span>
0094 
0095 
0096 
0097 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0098 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0099 <a name="_sub1" href="#_subfunctions" class="code">function [ avw ] = read_image(avw,IMGorient,machine)</a>
0100 
0101 fid = fopen(sprintf(<span class="string">'%s.img'</span>,avw.fileprefix),<span class="string">'r'</span>,machine);
0102 <span class="keyword">if</span> fid &lt; 0,
0103   msg = sprintf(<span class="string">'...cannot open file %s.img\n\n'</span>,avw.fileprefix);
0104   error(msg);
0105 <span class="keyword">end</span>
0106 
0107 ver = <span class="string">'[$Revision: 1.9 $]'</span>;
0108 fprintf(<span class="string">'\nAVW_IMG_READ [v%s]\n'</span>,ver(12:16));  tic;
0109 
0110 <span class="comment">% short int bitpix;    /* Number of bits per pixel; 1, 8, 16, 32, or 64. */</span>
0111 <span class="comment">% short int datatype      /* Datatype for this image set */</span>
0112 <span class="comment">% /*Acceptable values for datatype are*/</span>
0113 <span class="comment">% #define DT_NONE             0</span>
0114 <span class="comment">% #define DT_UNKNOWN          0    /*Unknown data type*/</span>
0115 <span class="comment">% #define DT_BINARY           1    /*Binary             ( 1 bit per voxel)*/</span>
0116 <span class="comment">% #define DT_UNSIGNED_CHAR    2    /*Unsigned character ( 8 bits per voxel)*/</span>
0117 <span class="comment">% #define DT_SIGNED_SHORT     4    /*Signed short       (16 bits per voxel)*/</span>
0118 <span class="comment">% #define DT_SIGNED_INT       8    /*Signed integer     (32 bits per voxel)*/</span>
0119 <span class="comment">% #define DT_FLOAT           16    /*Floating point     (32 bits per voxel)*/</span>
0120 <span class="comment">% #define DT_COMPLEX         32    /*Complex (64 bits per voxel; 2 floating point numbers)/*</span>
0121 <span class="comment">% #define DT_DOUBLE          64    /*Double precision   (64 bits per voxel)*/</span>
0122 <span class="comment">% #define DT_RGB            128    /*A Red-Green-Blue datatype*/</span>
0123 <span class="comment">% #define DT_ALL            255    /*Undocumented*/</span>
0124 
0125 <span class="keyword">switch</span> double(avw.hdr.dime.bitpix),
0126   <span class="keyword">case</span>  1,   precision = <span class="string">'bit1'</span>;
0127   <span class="keyword">case</span>  8,   precision = <span class="string">'uchar'</span>;
0128   <span class="keyword">case</span> 16,   precision = <span class="string">'int16'</span>;
0129   <span class="keyword">case</span> 32,
0130     <span class="keyword">if</span>     isequal(avw.hdr.dime.datatype, 8), precision = <span class="string">'int32'</span>;
0131     <span class="keyword">else</span>                                      precision = <span class="string">'single'</span>;
0132     <span class="keyword">end</span>
0133   <span class="keyword">case</span> 64,   precision = <span class="string">'double'</span>;
0134   <span class="keyword">otherwise</span>,
0135     precision = <span class="string">'uchar'</span>;
0136     fprintf(<span class="string">'...precision undefined in header, using ''uchar''\n'</span>);
0137 <span class="keyword">end</span>
0138 
0139 <span class="comment">% read the whole .img file into matlab (faster)</span>
0140 fprintf(<span class="string">'...reading %s Analyze %s image format.\n'</span>,machine,precision);
0141 fseek(fid,0,<span class="string">'bof'</span>);
0142 <span class="comment">% adjust for matlab version</span>
0143 ver = version;
0144 ver = str2num(ver(1));
0145 <span class="keyword">if</span> ver &lt; 6,
0146   tmp = fread(fid,inf,sprintf(<span class="string">'%s'</span>,precision));
0147 <span class="keyword">else</span>,
0148   tmp = fread(fid,inf,sprintf(<span class="string">'%s=&gt;double'</span>,precision));
0149 <span class="keyword">end</span>
0150 fclose(fid);
0151 
0152 <span class="comment">% Update the global min and max values</span>
0153 avw.hdr.dime.glmax = max(double(tmp));
0154 avw.hdr.dime.glmin = min(double(tmp));
0155 
0156 
0157 <span class="comment">%---------------------------------------------------------------</span>
0158 <span class="comment">% Now partition the img data into xyz</span>
0159 
0160 <span class="comment">% --- first figure out the size of the image</span>
0161 
0162 <span class="comment">% short int dim[ ];      /* Array of the image dimensions */</span>
0163 <span class="comment">%</span>
0164 <span class="comment">% dim[0]      Number of dimensions in database; usually 4.</span>
0165 <span class="comment">% dim[1]      Image X dimension;  number of pixels in an image row.</span>
0166 <span class="comment">% dim[2]      Image Y dimension;  number of pixel rows in slice.</span>
0167 <span class="comment">% dim[3]      Volume Z dimension; number of slices in a volume.</span>
0168 <span class="comment">% dim[4]      Time points; number of volumes in database.</span>
0169 
0170 PixelDim = double(avw.hdr.dime.dim(2));
0171 RowDim   = double(avw.hdr.dime.dim(3));
0172 SliceDim = double(avw.hdr.dime.dim(4));
0173 
0174 PixelSz  = double(avw.hdr.dime.pixdim(2));
0175 RowSz    = double(avw.hdr.dime.pixdim(3));
0176 SliceSz  = double(avw.hdr.dime.pixdim(4));
0177 
0178 
0179 
0180 
0181 
0182 <span class="comment">% ---- NON STANDARD ANALYZE...</span>
0183 
0184 <span class="comment">% Some Analyze files have been found to set -ve pixdim values, eg</span>
0185 <span class="comment">% the MNI template avg152T1_brain in the FSL etc/standard folder,</span>
0186 <span class="comment">% perhaps to indicate flipped orientation?  If so, this code below</span>
0187 <span class="comment">% will NOT handle the flip correctly!</span>
0188 <span class="keyword">if</span> PixelSz &lt; 0,
0189   warning(<span class="string">'X pixdim &lt; 0 !!! resetting to abs(avw.hdr.dime.pixdim(2))'</span>);
0190   PixelSz = abs(PixelSz);
0191   avw.hdr.dime.pixdim(2) = single(PixelSz);
0192 <span class="keyword">end</span>
0193 <span class="keyword">if</span> RowSz &lt; 0,
0194   warning(<span class="string">'Y pixdim &lt; 0 !!! resetting to abs(avw.hdr.dime.pixdim(3))'</span>);
0195   RowSz = abs(RowSz);
0196   avw.hdr.dime.pixdim(3) = single(RowSz);
0197 <span class="keyword">end</span>
0198 <span class="keyword">if</span> SliceSz &lt; 0,
0199   warning(<span class="string">'Z pixdim &lt; 0 !!! resetting to abs(avw.hdr.dime.pixdim(4))'</span>);
0200   SliceSz = abs(SliceSz);
0201   avw.hdr.dime.pixdim(4) = single(SliceSz);
0202 <span class="keyword">end</span>
0203 
0204 <span class="comment">% ---- END OF NON STANDARD ANALYZE</span>
0205 
0206 
0207 
0208 
0209 
0210 <span class="comment">% --- check the orientation specification and arrange img accordingly</span>
0211 <span class="keyword">if</span> ~isempty(IMGorient),
0212   <span class="keyword">if</span> ischar(IMGorient),
0213     avw.hdr.hist.orient = uint8(str2num(IMGorient));
0214   <span class="keyword">else</span>
0215     avw.hdr.hist.orient = uint8(IMGorient);
0216   <span class="keyword">end</span>
0217 <span class="keyword">end</span>,
0218 
0219 <span class="keyword">if</span> isempty(avw.hdr.hist.orient),
0220   msg = [ <span class="string">'...unspecified avw.hdr.hist.orient, using default 0\n'</span>,<span class="keyword">...</span>
0221       <span class="string">'   (check image and try explicit IMGorient option).\n'</span>];
0222   fprintf(msg);
0223   avw.hdr.hist.orient = uint8(0);
0224 <span class="keyword">end</span>
0225 
0226 <span class="keyword">switch</span> double(avw.hdr.hist.orient),
0227   
0228   <span class="keyword">case</span> 0, <span class="comment">% transverse unflipped</span>
0229     
0230     <span class="comment">% orient = 0:  The primary orientation of the data on disk is in the</span>
0231     <span class="comment">% transverse plane relative to the object scanned.  Most commonly, the fastest</span>
0232     <span class="comment">% moving index through the voxels that are part of this transverse image would</span>
0233     <span class="comment">% span the right-left extent of the structure imaged, with the next fastest</span>
0234     <span class="comment">% moving index spanning the posterior-anterior extent of the structure.  This</span>
0235     <span class="comment">% 'orient' flag would indicate to Analyze that this data should be placed in</span>
0236     <span class="comment">% the X-Y plane of the 3D Analyze Coordinate System, with the Z dimension</span>
0237     <span class="comment">% being the slice direction.</span>
0238     
0239     <span class="comment">% For the 'transverse unflipped' type, the voxels are stored with</span>
0240     <span class="comment">% Pixels in 'x' axis (varies fastest) - from patient right to left</span>
0241     <span class="comment">% Rows in   'y' axis                  - from patient posterior to anterior</span>
0242     <span class="comment">% Slices in 'z' axis                  - from patient inferior to superior</span>
0243     
0244     fprintf(<span class="string">'...reading axial unflipped orientation\n'</span>);
0245     
0246     avw.img = zeros(PixelDim,RowDim,SliceDim);
0247     
0248     n = 1;
0249     x = 1:PixelDim;
0250     <span class="keyword">for</span> z = 1:SliceDim,
0251       <span class="keyword">for</span> y = 1:RowDim,
0252         <span class="comment">% load Y row of X values into Z slice avw.img</span>
0253         avw.img(x,y,z) = tmp(n:n+(PixelDim-1));
0254         n = n + PixelDim;
0255       <span class="keyword">end</span>
0256     <span class="keyword">end</span>
0257     
0258     <span class="comment">% no need to rearrange avw.hdr.dime.dim or avw.hdr.dime.pixdim</span>
0259     
0260     
0261   <span class="keyword">case</span> 1, <span class="comment">% coronal unflipped</span>
0262     
0263     <span class="comment">% orient = 1:  The primary orientation of the data on disk is in the coronal</span>
0264     <span class="comment">% plane relative to the object scanned.  Most commonly, the fastest moving</span>
0265     <span class="comment">% index through the voxels that are part of this coronal image would span the</span>
0266     <span class="comment">% right-left extent of the structure imaged, with the next fastest moving</span>
0267     <span class="comment">% index spanning the inferior-superior extent of the structure.  This 'orient'</span>
0268     <span class="comment">% flag would indicate to Analyze that this data should be placed in the X-Z</span>
0269     <span class="comment">% plane of the 3D Analyze Coordinate System, with the Y dimension being the</span>
0270     <span class="comment">% slice direction.</span>
0271     
0272     <span class="comment">% For the 'coronal unflipped' type, the voxels are stored with</span>
0273     <span class="comment">% Pixels in 'x' axis (varies fastest) - from patient right to left</span>
0274     <span class="comment">% Rows in   'z' axis                  - from patient inferior to superior</span>
0275     <span class="comment">% Slices in 'y' axis                  - from patient posterior to anterior</span>
0276     
0277     fprintf(<span class="string">'...reading coronal unflipped orientation\n'</span>);