spm_write_plane

 Write a transverse plane of image data.
 FORMAT V = spm_write_plane(V,A,p)
 V   - data structure containing image information.
       - see spm_vol for a description.
 A   - the two dimensional image to write.
 p   - the plane number (beginning from 1).

 VO  - (possibly) modified data structure containing image information.
       It is possible that future versions of spm_write_plane may
       modify scalefactors (for example).

_______________________________________________________________________
 @(#)spm_write_plane.m    2.19 John Ashburner 03/07/16

0001 function V = spm_write_plane(V,A,p)
0002 % Write a transverse plane of image data.
0003 % FORMAT V = spm_write_plane(V,A,p)
0004 % V   - data structure containing image information.
0005 %       - see spm_vol for a description.
0006 % A   - the two dimensional image to write.
0007 % p   - the plane number (beginning from 1).
0008 %
0009 % VO  - (possibly) modified data structure containing image information.
0010 %       It is possible that future versions of spm_write_plane may
0011 %       modify scalefactors (for example).
0012 %
0013 %_______________________________________________________________________
0014 % @(#)spm_write_plane.m    2.19 John Ashburner 03/07/16
0015 
0016 if any(V.dim(1:2) ~= size(A)), error('Incompatible image dimensions');      end;
0017 if p>V.dim(3),                 error('Plane number too high');              end;
0018 
0019 % Write Analyze image by default
0020 V = write_analyze_plane(V,A,p);
0021 return;
0022 %_______________________________________________________________________
0023 
0024 %_______________________________________________________________________
0025 function V = write_analyze_plane(V,A,p)
0026 
0027 types   = [    2      4      8   16   64   130    132    136,   512   1024   2048 4096 16384 33280  33792  34816];
0028 maxval  = [2^8-1 2^15-1 2^31-1  Inf  Inf 2^7-1 2^16-1 2^32-1, 2^8-1 2^15-1 2^31-1  Inf   Inf 2^8-1 2^16-1 2^32-1];
0029 minval  = [    0  -2^15  -2^31 -Inf -Inf  -2^7      0      0,     0  -2^15  -2^31 -Inf  -Inf  -2^7      0      0];
0030 intt    = [    1      1      1    0    0     1      1      1,     1      1      1    0     0     1      1      1];
0031 prec = str2mat('uint8','int16','int32','float','double','int8','uint16','uint32','uint8','int16','int32','float','double','int8','uint16','uint32');
0032 swapped = [    0      0      0    0    0     0      0      0,     1      1      1    1     1     1      1      1];
0033 bits    = [    8     16     32   32   64     8     16     32,     8     16     32   32    64     8     16     32];
0034 
0035 dt      = find(types==V.dim(4));
0036 if isempty(dt), error('Unknown datatype'); end;
0037 
0038 A = double(A);
0039 
0040 % Rescale to fit appropriate range
0041 if intt(dt),
0042     A(isnan(A)) = 0;
0043     mxv         = maxval(dt);
0044     mnv         = minval(dt);
0045     A           = round(A*(1/V.pinfo(1)) - V.pinfo(2));
0046     A(A > mxv)  = mxv;
0047     A(A < mnv)  = mnv;
0048 end;
0049 
0050 if ~isfield(V,'private') | ~isfield(V.private,'fid') | isempty(V.private.fid),
0051     mach = 'native';
0052     if swapped(dt),
0053         if spm_platform('bigend'), mach = 'ieee-le'; else, mach = 'ieee-be'; end;
0054     end; 
0055     [pth,nam,ext] = fileparts(V.fname);
0056     fname         = fullfile(pth,[nam, '.img']);
0057     fid           = fopen(fname,'r+',mach);
0058     if fid == -1,
0059         fid   = fopen(fname,'w',mach);
0060         if fid == -1,
0061             error(['Error opening ' fname '. Check that you have write permission.']);
0062         end;
0063     end;
0064 else,
0065     if isempty(fopen(V.private.fid)),
0066         mach = 'native';
0067         if swapped(dt),
0068             if spm_platform('bigend'), mach = 'ieee-le'; else, mach = 'ieee-be'; end;
0069         end;
0070         V.private.fid = fopen(fname,'r+',mach);
0071         if V.private.fid == -1,
0072             error(['Error opening ' fname '. Check that you have write permission.']);
0073         end;
0074     end;
0075     fid = V.private.fid;
0076 end;
0077 
0078 % Seek to the appropriate offset
0079 datasize = bits(dt)/8;
0080 off   = (p-1)*datasize*prod(V.dim(1:2)) + V.pinfo(3,1);
0081 fseek(fid,0,'bof'); % A bug in Matlab 6.5 means that a rewind is needed
0082 if fseek(fid,off,'bof')==-1,
0083     % Need this because fseek in Matlab does not seek past the EOF
0084     fseek(fid,0,'bof'); % A bug in Matlab 6.5 means that a rewind is needed
0085     fseek(fid,0,'eof');
0086     curr_off = ftell(fid);
0087     blanks   = zeros(off-curr_off,1);
0088     if fwrite(fid,blanks,'uchar') ~= prod(size(blanks)),
0089         write_error_message(V.fname);
0090         error(['Error writing ' V.fname '.']);
0091     end;
0092     fseek(fid,0,'bof'); % A bug in Matlab 6.5 means that a rewind is needed
0093     if fseek(fid,off,'bof') == -1,
0094         write_error_message(V.fname);
0095         error(['Error writing ' V.fname '.']);
0096         return;
0097     end;
0098 end;
0099 
0100 if fwrite(fid,A,deblank(prec(dt,:))) ~= prod(size(A)),
0101     write_error_message(V.fname);
0102     error(['Error writing ' V.fname '.']);
0103 end;
0104 
0105 if ~isfield(V,'private') | ~isfield(V.private,'fid') | isempty(V.private.fid), fclose(fid); end;
0106 
0107 return;
0108 %_______________________________________________________________________
0109 
0110 %_______________________________________________________________________
0111 function write_error_message(q)
0112 str = {...
0113     'Error writing:',...
0114     ' ',...
0115     ['        ',spm_str_manip(q,'k40d')],...
0116     ' ',...
0117     'Check disk space / disk quota.'};
0118 spm('alert*',str,mfilename,sqrt(-1));
0119 
0120 return;
0121 %_______________________________________________________________________