avw_img_read_4d.html

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

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