% write_Space_volume( vol, filestem )
% --------------------------------------------
%
% Writes a 3-D array into a Space Volume Type 1 file, readable
% by Space software as an 8 bit grayscale volume. Values are scaled
% to the range [ 0, 255 ] based on the minimum and maximum array 
% values. The array will be displayed by Space in an x/y/z coordinate
% system corresponding to the first/second/third array dimensions.
%
% Input:  
%  
%    vol           - a 3-D array of type double
%
%    filestem   - a character string used to form the stem of the output
%                       file name, "[filestem].vol".
%
% Output:
%
%     "[filestem].vol"                
%                   - A file containing an 8 bit grayscale volume that conforms
%                     to the Space Volume Type 1 file specification.
%
%
%    Optional inputs:
%
%        global VOX (double)
%                    - If this variable is globally defined, and is a 3 element
%                      vector of positive values, the values are used as the x/y/z 
%                      physical dimensions of each voxel. Otherwise, the physical 
%                      dimensions are set to 0.
%
%       global TITLE (char)
%                    - If this variable is globally defined, up to 151 characters will
%                      be used as the title in the Space volume.
%
%       global DESCRIPTION (char)
%                    - If this variable is globally defined, up to 4900 characters will
%                      be used as the title in the Space volume.
%
% For details regarding Space software and the Space native file type, see:
%
% http://lcni.uoregon.edu/~mark/Space_program.html
%
% last modified: June 2004
% Mark Dow
% University of Oregon
% Brain Development Lab / Lewis Center for Neuroimaging
% dow@uoregon.edu
%

function write_Space_volume( vol, filestem )

global VOX;
global TITLE;
global DESCRIPTION;

outfilename = [ filestem '.vol'];
fprintf('\nWriting 3-D array to Space volume file...\n\n'); 

% Scale and shift array values to between 0 and 255.
vol = double(vol);
min_value = min(min(min(vol)));
max_value = max(max(max(vol)));
vol = vol - min_value;
vol = vol./max(max(max(vol)));
vol = vol.*255;
fprintf([ '    minimum input array value  ' num2str(min_value) ' scaled and shifted to 0\n' ]); 
fprintf([ '    maximum input array value  ' num2str(max_value) ' scaled and shifted to 255\n' ]); 

% Open file for write with little-endian byte ordering.
fid = fopen( outfilename, 'w', 'l' );  

% Write identifier for any Space volume type (5 char's - 5 x 8 bits)
file_id = 'mdvol';
count = fwrite( fid, file_id, 'char' );

% Write Space Volume Type file version indicator (1 char - 8 bits)
file_version = '1';
count = fwrite( fid, file_version, 'char' );

% Write header length (1 int - 1 x 32 bits)
header_length = 10000;
count = fwrite( fid, header_length, 'int' );

% Write volume (3-D array) dimensions (3 int's - 3 x 32 bits)
dim = size(vol);
fprintf( [ '\n    Volume dimensions (voxels):  ' num2str(dim(1)) ' x ' num2str(dim(2)) ' x ' num2str(dim(3)) '\n' ] ); 
dim_reordered(1) = dim(1);
dim_reordered(2) = dim(3);
dim_reordered(3) = dim(2);
count = fwrite( fid, dim_reordered, 'int' );
  
% Write pixel sizes (3 float's - 3 x 32 bits)
voxel_size = [0,0,0];
if length( VOX ) == 3  &  VOX(1) > 0  &  VOX(2) > 0  &  VOX(3) > 0
    voxel_size = VOX;
    fprintf( [ '\n    Voxel size (mm):  ' num2str(voxel_size(1)) ' x ' num2str(voxel_size(2)) ' x ' num2str(voxel_size(3)) '\n' ] ); 
else
    fprintf( [ '\n    No voxel size specified.\n' ] );
end

count = fwrite( fid, voxel_size, 'float32' );

% Write 8 bit grayscale data format code (3 char's - 3 x 8 bits)
data_format = 'g08';
count = fwrite( fid, data_format, 'char' );

% Write initial black and white point for display (2 float - 2 x 32 bits)
bw_point = [ 0.0, 1.0 ];
count = fwrite( fid, bw_point, 'float32' );
  
% Write initial gamma for display (1 float - 1 x 32 bits)
gamma = 1.0;
count = fwrite( fid, gamma, 'float32' );
  
% Write file format description text (4900 char's - 1500 x 8 bits)
for i = 1:4900
    file_description(i) = ' ';
end
text = '  This is a simple file, of Space Volume Type 1, consisting of a header followed by a set of image value bytes, for storing volume image data and basic image display parameters. The header (which should be in little-endian byte order) is a total of 10000 bytes long. The first five characters (mdvol, 5 x 1 byte) identify the file as volume with this file type. One character, 1 byte, serves as a file type version indicator -- this file is of type 1. One integer, 4 bytes, is the total length of the header. Three integers are dimensions, xzy in pixels, of the volume. Three floats, 3x4 bytes, are physical dimensions, xzy in mm, of each voxel. Two floats indicate a reasonable black point and white point for the initial display color transform, with a range [0,1]. One float indicates a reasonable gamma (>0) for the initial display. Three characters, the color code, specifies the color format that the image values represent: gray 8 bit (g08), indexed 8 bit (i08) or color 24 bit (c24). 4900 characters, 4900 x 1 byte, of text is this description of the file format. 151 characters are a text title. 4900 characters are a text description of the volume image content. The rest are image value bytes (8 bit integers), describing the volume image. If the color code is g08, there are x*y*z bytes of image values, in -x/-z/-y order. If the color code is i08, there are x*y*z bytes of image values, in -x/-z/-y order. If the color code is c24, there are 3*x*y*z bytes of image color values, in rgb/-x/-z/-y order.  Specified by Mark Dow, May 2004. ';

file_description(1:length(text) ) = text;
count = fwrite( fid, file_description, 'char' );

% Write volume title text (151 char's - 151 x 8 bits)
for i = 1:151
    volume_description(i) = ' ';
end
if length(TITLE) > 0
    volume_description( 1:min( length(TITLE), 151 ) ) = TITLE( 1:min( length(TITLE), 151 ) );
    fprintf( [ '\n    Volume title:  ' volume_description '\n' ] ); 
end
count = fwrite( fid, volume_description, 'char' );

% Write volume description text (4900 char's - 3 x 8 bits)
for i = 1:4900
    volume_description(i) = ' ';
end
if length(DESCRIPTION) > 0
    volume_description( 1:min( length(DESCRIPTION), 4900 ) ) = DESCRIPTION( 1:min( length(DESCRIPTION), 4900 ) );
    fprintf( [ '\n    Volume description:  ' volume_description '\n' ] ); 
end
count = fwrite( fid, volume_description, 'char' );

fprintf( [ '\n    Ordering and writing bytes...\n' ] );

% Rearrange bytes in -x/-z/-y order.
vol = permute( vol, [1 3 2] );
vol = flipdim( vol, 1 );
vol = flipdim( vol, 2 );
vol = flipdim( vol, 3 );

% Write grayscale bytes (uxvxw uchar's - uxvxw x 8 bits)
count = fwrite( fid, vol, 'char' );
 
fprintf( [ '\nDone writing file named: ' outfilename '\n\n' ] );

fclose(fid);