function Geometry = GeometryGenerator()

%% % % GeometryGenerator % % % % %
%  This function builds the geometry files needed for the WiscPlan dose 
%  optimization.  Requires that you have my amiraStructReader function. You 
%  need contours in any amira format. Contours should ALL BE CROPPED TO THE
%  SAME DIMENSION.
%
%  DLB - 6/5/09

s = pwd;

% Get the planning CT
[planningCT, CT_pathname, filterindex] = uigetfile({'*.am','Amira file (*.am)'},'Select the planning CT.','MultiSelect','off');
cd(CT_pathname);

% Generate a list of contours to include in Geometry.ROIS
[files, pathname, filterindex] = uigetfile({'*.am','Amira file (*.am)'; '*.*','All files (*.*)'},'Select all the contours for the Geometry file (use Ctrl).','MultiSelect','on');
cd(pathname);

% Set the initial Geometry structure members
data = amiraStructReader([CT_pathname planningCT]);
Geometry.voxel_size = [data.xpix data.ypix data.zpix];
Geometry.start = [data.boundingbox(1) data.boundingbox(3) data.boundingbox(5)];
Geometry.data = single(data.image);

%% Adding ROIS to Geometry
for i = 1:length(files)
    data = amiraStructReader([pathname files{i}]);
    % remove the .am file extension
    name = strrep(files{i},'.am','');
    % remove everything before the first underscore
    name = name(findstr(name,'_')+1:length(name));
    % anything that remotely has the letters GTV is renamed GTV
    if ~isempty(findstr(name,'GTV'))
            name = 'GTV';
    end
    
    % assign parameters
    Geometry.ROIS{i}.name = name;
    Geometry.ROIS{i}.ind = int32(find(data.image));
    Geometry.ROIS{i}.dim = [data.dims(1) data.dims(2) data.dims(3)];
    Geometry.ROIS{i}.pix_dim = [data.xpix data.ypix data.zpix];
    Geometry.ROIS{i}.start = [data.boundingbox(1) data.boundingbox(3) data.boundingbox(5)];
    Geometry.ROIS{i}.start_ind = '(1,1,1)';
end


cd(pathname);

% save the Geometry file in the same location as your contours.
save('Geometry.mat','Geometry');

% return to the original directory
cd(s);


%% This is the amiraStructReader in case you don't have it. %%
function data = amiraStructReader(varargin)

%% % % amiraStructReader % % % % %
%  Built for when you might need additional amira header info in matlab.
%  Same functionality as amiraReader, but adds stuff in structure format 
%
%  DLB - 5/5/09  Cinco de Mayo!!!  But another 12 hour day :*(
%

warning off all;  % you'll get truncation warnings because you read inf in line 24.

cd('C:\Data\');

switch nargin 
    case 1
        filename = varargin{1};
    case 0
        [file, pathname] = uigetfile({'*.am'},'Please locate Amira file.','Multiselect','off');
        filename = [pathname file];
end

fid = fopen(filename,'r','ieee-be');
datachar = char(fread(fid,inf,'char'))'; % this will be the char data in the header
switch isempty(findstr(datachar,'BINARY-LITTLE-ENDIAN'))
    case 1
        endianness = 'be';
    case 0
        endianness = 'le';
end
lattice = findstr(datachar,'define Lattice');
parameter = findstr(datachar,'Parameters');
start_data_position = findstr(datachar,'@1');  % note there will be two @1, so take the second one!

fseek(fid,lattice+14,'bof');
data.dims = str2num(fgetl(fid));

bounding = findstr(datachar,'Bounding');
fseek(fid,bounding + 11,'bof');
data.boundingbox = str2num(fgetl(fid));


data.x = (data.boundingbox(2)-data.boundingbox(1));
data.y = (data.boundingbox(4)-data.boundingbox(3));
data.z = (data.boundingbox(6)-data.boundingbox(5))*(1+1/(data.dims(3)-1));

data.xpix = data.x/(data.dims(1)-1);
data.ypix = data.y/(data.dims(2)-1);
data.zpix = data.z/data.dims(3);

%% Determine if we're dealing with a labelfield or regular data
switch isempty(findstr(datachar,'HxByteRLE'))
    case 1  % regular data case
        data_type = char(datachar(findstr(datachar,'Lattice { ')+10:findstr(datachar,'Data } ')-2));
        switch data_type
            case 'byte'
                out_type = 'int8';
            case 'short'
                out_type = 'int16';
            case 'ushort'
                out_type = 'uint16'; 
            case 'int32'
                out_type = 'int32';
            case 'float'
                out_type = 'float';
            case 'double' 
                out_type = 'double';
        end
        fid = fopen(filename,'r',['ieee-' endianness]);
        fseek(fid,start_data_position(2)+2,'bof');
        data.image = reshape(fread(fid,data.dims(1)*data.dims(2)*data.dims(3),out_type),data.dims(1),data.dims(2),data.dims(3));
 
    case 0 % labelfield case
        ByteRLE = findstr(datachar,'HxByteRLE,');
        DataSection = findstr(datachar,'# Data section follows');
        
        fseek(fid,ByteRLE+9,'bof');
        byte_length = str2num(char(fread(fid,DataSection - ByteRLE-13,'char'))');

        switch endianness
            case 'le'
                fclose(fid);
                fid = fopen(filename,'r',['ieee-' endianness]);
                fseek(fid,start_data_position(2)+2,'bof');
                data.image = fread(fid,byte_length,'int8');
            case 'be'
                fseek(fid,start_data_position(2)+2,'bof');
                data.image = fread(fid,byte_length,'int8');
        end

        out_line = zeros(data.dims(1)*data.dims(2)*data.dims(3),1);

        final_offset = sum(data.image(find(data.image<0))+127);

        net_counter = 0;  % this sifts through the out_line variable
        offset = 0;  % this adds an offset to data_counter when writing VERBATIM

        for i = 1:floor((byte_length-final_offset)/2) % 
            data_counter = 2*i+offset;  % this will sift through the data variable
            if ((data.image(data_counter) == 0) && (data.image(data_counter-1) > 0)) % check if the value is 0 AND that we're not writing VERBATIM
                net_counter = data.image(data_counter-1)+net_counter;
            else
                if data.image(data_counter-1) < 0 % check if we're writing VERBATIM
                    for j = 1:128+data.image(data_counter-1)
                        out_line(j+net_counter) = data.image(data_counter+j-1);
                    end
                    net_counter = 128+data.image(data_counter-1)+net_counter;
                    offset = offset + 127+ data.image(data_counter-1);
                else  % assume that we're copying value data(data_counter), data(data_counter-1) times.
                    for j = 1:data.image(data_counter-1)
                        out_line(j+net_counter) = data.image(data_counter);
                    end
                    net_counter = data.image(data_counter-1)+net_counter;
                end
            end
        end 

        data.image = reshape(out_line,data.dims(1),data.dims(2),data.dims(3));
end