WiscPlan_olderMV/photonDoseCalc/scripts_m/am2geom.m

function Geometry = am2geom(AmiraMeshFile)
% Geometry = am2geom(AmiraMeshFile) converts an AmiraMesh (.am) file to a 
% geometry structure, as long as the .am file is written in the recognized 
% format.  The Geometry structure has the form:
%
% Geometry = 
%          start: [3x1 double]
%     voxel_size: [3x1 double]
%           data: [MxNxQ double]
%
% Example acceptable .am file format:
%
% # AmiraMesh 3D BINARY 2.0
% 
% # CreationDate: Thu Jun 01 15:42:59 2006
% 
% 
% define Lattice 509 509 3
% 
% Parameters {
% Content "509x509x36 short, uniform coordinates",
% BoundingBox -25.48 25.32 -25.48 25.32 9 26.5,
% CoordType "uniform"
% }
% 
% Lattice { short Data } @1
% 
% # Data section follows
% @1
%
% Ryan T Flynn 23 August 2007
% 

% ensure that the user input a string
if ~ischar(AmiraMeshFile)
    error('AmiraMeshFile must be a character string.');
end

% open the file
fid = fopen(AmiraMeshFile,'r');
if fid == -1
    error('Could not open %s.',AmiraMeshFile);
end

% Things to look for in the .am file:
% 1) The 'BINARY' or 'ASCII' words
% 2) The lattice size after the 'define Lattice' words
% 3) The BoundingBox line, which is used to calculate start and voxel_size
% 4) CoordType "uniform" line

tline = fgetl(fid);

% continue to read the file until the line containing only '@1' is found
while ~strncmp(tline,'@1',2)
    if ~isempty(regexp(tline,'AmiraMesh 3D'))
        % found the BINARY or ASCII specification line
        p = regexp(tline,'AmiraMesh 3D');  % find the position of the key words
        r = tline(p+length('AmiraMesh 3D'):end);  % get the remainder of the string
        binORascii = strtok(r);  % pop off the data type (either BINARY or ASCII)
    end
    if ~isempty(regexp(tline,'define Lattice'))
        % found the lattice definition line
        p = regexp(tline,'define Lattice');  % find the position of the key words
        r = tline(p+length('define Lattice'):end);  % get the remainder of the string
        siz = str2num(r);  % convert the remainder to the lattice dimensions
        if length(siz) ~= 3
            error('Data lattice must be 3D.');
        end
    end
    
    if ~isempty(regexp(tline,'BoundingBox'))
        % found the bounding box line
        p = regexp(tline,'BoundingBox');  % position of BoundingBox word
        r = tline(p+length('BoundingBox'):end); % remainder of the string
        box = str2num(r);  % convert the remainder to the box dimensions
        if length(box) ~= 6
            error('Number of parameters in BoundingBox line must be 6.');
        end
    end
    % Lattice { short Data } @1
    if ~isempty(regexp(tline,'Lattice {'))
        % found line that contains the data type for the grid
        p = regexp(tline,'Lattice {');  % position of BoundingBox word
        r = tline(p+length('Lattice {'):end); % remainder of the string
        AmiraDataType = strtok(r);  % pop off the data type
        % Figure out what the data type is
        for k=1:6
            if strcmp(AmiraDataType,'byte')
                MatlabDataType = 'int8';
            elseif strcmp(AmiraDataType,'short')
                MatlabDataType = 'int16';
            elseif strcmp(AmiraDataType,'ushort')
                MatlabDataType = 'uint16';
            elseif strcmp(AmiraDataType,'int32')
                MatlabDataType = 'int32';
            elseif strcmp(AmiraDataType,'float')
                MatlabDataType = 'single';
            elseif strcmp(AmiraDataType,'double')
                MatlabDataType = 'double';
            else
                error('Amira data type in %s is not supported by Matlab.',AmiraMeshFile);
            end
        end
    end
    tline = fgetl(fid);
end

% ensure that all of the required lines were found
if ~exist('binORascii')
    error('Could not determine of data are stored in binary or ascii format.');
end

if ~exist('siz')
    error('Could not find the data dimensions during the reading process.');
end

if ~exist('box')
    error('Could not find the bounding box in the reading process.');
end

if ~exist(MatlabDataType)
    error('Could not find the data type for file %s.',AmiraMeshFile);
end

% ensure that the data are stored in binary format
if ~strcmp(binORascii,'BINARY')
    error('Data in %s must be stored in binary format.',AmiraMeshFile);
end

% calculate the voxel_size vector
if siz(1) == 1
    Geometry.voxel_size(1) = box(2) - box(1);
elseif siz(1) > 1
    Geometry.voxel_size(1) = (box(2) - box(1))/(siz(1)-1);
else
    error('First dimension of data lattice must be greater than zero.');
end

if siz(2) == 1
    Geometry.voxel_size(2) = box(4) - box(3);
elseif siz(2) > 1
    Geometry.voxel_size(2) = (box(4) - box(3))/(siz(2)-1);
else
    error('Second dimension of data lattice must be greater than zero.');
end

if siz(3) == 1
    Geometry.voxel_size(3) = box(6) - box(5);
elseif siz(3) > 1
    Geometry.voxel_size(3) = (box(6) - box(5))/(siz(3)-1);
else
    error('Third dimension of data lattice must be greater than zero.');
end

% calculate the start vector
Geometry.start(1) = box(1);
Geometry.start(2) = box(3);
Geometry.start(3) = box(5);

% allocate space for the data
Geometry.data = zeros(siz,MatlabDataType);

% read in the data
Geometry.data = fread(fid,[MatlabDataType '=>' MatlabDataType],'ieee-be');

% If reading byte-by-byte, throw out last byte -- it's just a carriage return
if strcmp(MatlabDataType,'int8')
    Geometry.data(end,:) = [];
end

Geometry.data = reshape(Geometry.data,siz(1),siz(2),siz(3));

% close the AmiraMesh file
fclose(fid);