% pinn2matlab.m
%
%  Ryan T Flynn 23 August 2007
%
% Reads in raw Pinnacle data to a MATLAB data structure. The patient CT,
% ROIs, and Trial files can all be read in if available and desired. The
% result is a structure called 'Geometry', which contains ROIS and Trial
% fields.
%
% The CT image and all of the ROIs are flipped onto a coordinate system in
% which the position of each voxel index (i,j,k) are given by:
%
% x(i) = x0 + i*delx; (i=0,...,M-1)
% y(j) = y0 + j*delx; (j=0,...,N-1)
% z(k) = z0 + k*delx; (k=0,...,Q-1)
% 
% As opposed to the Pinnacle coordinate system, which
% 
% x(i) = x0 + i*delx;        (i=0,...,M-1)
% y(j) = y0 + (N-1-j)*dely;  (j=0,...,N-1)
% z(k) = z0 + k*delz;        (k=0,...,Q-1)

%%%%% Start of user supplied inputs %%%%%

% Read-in flags
readTrial = 'no';  % either 'yes' or 'no'
readROIS = 'yes';   % either 'yes' or 'no'

% locations of the Pinnacle files
roifile = 'C:\Documents and Settings\Steve\Desktop\steveTPS\plan\plan.roi';
imageheaderfile = 'C:\Documents and Settings\Steve\Desktop\steveTPS\plan\ImageSet_0.header';
imagefile = 'C:\Documents and Settings\Steve\Desktop\steveTPS\plan\ImageSet_0.img';  % CT file
save_file = 'C:\Documents and Settings\Steve\Desktop\steveTPS\plan\lung.mat';

%%%%% End of user supplied inputs %%%%%

start_ind = '(1,1,1)';  % this tag will be added to the ROIS so it is known where the start voxel is

clear roi ROIS Trial Geometry CTimage;  % start with a clean slate of structures

% extract geometric information from the header file
fid_imageheaderfile = fopen(imageheaderfile,'r');

tline = fgets(fid_imageheaderfile);
while tline ~= -1
    % check the line for key words
    if length(findstr(tline,'x_dim')) & ~length(findstr(tline,'fname_index_start'))
        eval(tline);  % run the line to get x_dim
    elseif length(findstr(tline,'y_dim'))
        eval(tline);  % run the line to get y_dim
    elseif length(findstr(tline,'z_dim'))
        eval(tline);  % run the line to get z_dim
    elseif length(findstr(tline,'x_pixdim'))
        eval(tline);  % run the line to get x_pixdim
    elseif length(findstr(tline,'y_pixdim'))
        eval(tline);  % run the line to get y_pixdim
    elseif length(findstr(tline,'z_pixdim'))
        eval(tline);  % run the line to get z_pixdim
    elseif length(findstr(tline,'x_start'))
        eval(tline);  % run the line to get x_start
    elseif length(findstr(tline,'y_start'))
        eval(tline);  % run the line to get x_start
    elseif length(findstr(tline,'z_start'))
        eval(tline);  % run the line to get x_start
    end
    tline = fgets(fid_imageheaderfile);
end

fclose(fid_imageheaderfile);

% Read in the CT data
fid_image = fopen(imagefile,'rb');
CTimage = fread(fid_image,'short=>int16','ieee-be');
CTimage = reshape(CTimage,x_dim,y_dim,z_dim);
fclose(fid_image);

% flip the CT image onto the non-Pinnacle coordinate system
CTimage = flipdim(CTimage,2);

% save the results to the Geometry structure
Geometry.start = single([x_start y_start z_start]);
Geometry.voxel_size = single([x_pixdim y_pixdim z_pixdim]);
Geometry.data = single(CTimage)/1024; % convert the CT data to density data, Pinnacle style
Geometry.data(Geometry.data < 0) = 0; % truncate anything less than zero

clear CTimage;

if strcmp(readROIS,'yes')
    % read in the roi file
    fid_roifile = fopen(roifile,'r');

    roinames = {};  % start a cell array of the roi names
    Nroi = 0;  % number of rois read in

    % Flags to indicate which sets of angled brackets in the roi file tline is
    % inside.
    inroi = 0;
    incurve = 0;
    inpoints = 0;

    roi_num = 0;  % current roi number

    tline = fgets(fid_roifile);
    while tline ~= -1
        % check the line for key words
        if length(findstr(tline,'roi={'))
            inroi = 1;   % mark that we are now currently inside of an roi
            roi_num = roi_num + 1;
            % next line contains the roi name
            tline = fgets(fid_roifile);
            % pop off first token in line, the remainder of the line is the roi name
            [T,R] = strtok(tline);
            roi.name = strtrim(R);
            % pop off lines until we get to the number of curves in this roi
            while ~length(findstr(tline,'num_curve'))
                tline = fgets(fid_roifile);
            end
            % pop off the num_curve string
            [T,R] = strtok(tline);
            % pop off the equals sign
            [T,R] = strtok(R);
            % pop off the number of curves in this roi
            T = strtok(R,';');
            % save the number of curves to the roi stucture
            eval(['roi.num_curves = ' num2str(T) ';']);
            roi.curves = {};  % get the curves structure started
            % read in the next curve structure
            curve_num = 0;  % number of the current curve
            while roi.num_curves > 0 & curve_num < roi.num_curves
                while ~length(findstr(tline,'curve={'));
                    tline = fgets(fid_roifile);
                end
                curve_num = curve_num + 1;
                incurve = 1;  % inside the curve structure now
                % find the number of points in this structure
                while ~length(findstr(tline,'num_points'))
                    tline = fgets(fid_roifile);
                end
                % pop off the num_points string
                [T,R] = strtok(tline);
                % pop off the equals sign
                [T,R] = strtok(R);
                % pop off the number of points in this curve
                T = strtok(R,';');
                eval(['num_points = ' num2str(T) ';']);
                % find the points identifier
                while ~length(findstr(tline,'points={'))
                    tline = fgets(fid_roifile);
                end
                inpoints = 1;  % inside the points structure now
                % read in the block of points data
                block = fscanf(fid_roifile,'%g',[3 num_points]);
                % save the block of points to the roi stucture
                roi.curves{curve_num} = block';
                % read in the right parantheses for the points and curve
                % structures
                while ~length(findstr(tline,'};'))
                    tline = fgets(fid_roifile);
                end
                inpoints = 0;   % stepped outside of the points structure
                while ~length(findstr(tline,'};'))
                    tline = fgets(fid_roifile);
                end
                incurve = 0;  % stepped outside of the curves structure
            end
            % read until the roi closing bracket appears
            while ~length(findstr(tline,'};'))
                tline = fgets(fid_roifile);
            end
            inroi = 0;  % we are now outside of the roi
            fprintf('read in %s roi\n',roi.name);
            ROIS{roi_num} = roi;
        end
        tline = fgets(fid_roifile);
    end

    fclose(fid_roifile);

    % ROIS have all been read-in, now just have to convert them to mask
    % matrices.  We'll use roipoly for this.
    % Recall that we must use the Pinnacle coordinate system for this to work,
    % that is, (x,y,z) coordinates are given by:
    % x(i) = x_start + i*x_pixdim, i=0..x_dim-1
    % y(j) = y_start + (y_dim-1)*y_pixdim - j*y_pixdim, j=0..y_dim-1
    % z(k) = z_start + k*z_pixdim, k=0..z_dim-1
    %
    % Not all treatment planning systems use this type of coordinate system
    % definition, so it is very important to get them straight.
    %
    % To get around this we will manipulate the data such that we'll have:
    %
    % x(i) = x_start + i*x_pixdim, i=0..x_dim-1
    % y(j) = y_start + j*y_pixdim, j=0..y_dim-1
    % z(k) = z_start + k*z_pixdim, k=0..z_dim-1
    %
    % This can be done by a simple fliplr operation on all of the CT slices

    % define the coordinate system
    x = x_start:x_pixdim:x_start + (x_dim-1)*x_pixdim;
    y = y_start:y_pixdim:y_start + (y_dim-1)*y_pixdim;
    z = z_start:z_pixdim:z_start + (z_dim-1)*z_pixdim;

    % define the locations of the corners of the pixels in each slice
    xCorners = [x - x_pixdim/2 x(end) + x_pixdim/2];
    yCorners = [y - y_pixdim/2 y(end) + y_pixdim/2];
    
    % loop through all rois
    for roi_num=1:length(ROIS)
        % set up the roi mask
        roimask = zeros(x_dim,y_dim,z_dim,'int8');
        roimaskCorners = zeros(x_dim+1,y_dim+1,z_dim,'int8');
        % loop through all of the curves in the roi
        for curve_num=1:length(ROIS{roi_num}.curves)
            % make a copy of the curve for easy access
            current_curve = ROIS{roi_num}.curves{curve_num};
            % find the z-index (slice number) of this curve
            q = round((current_curve(1,3)-z_start)/z_pixdim) + 1;
            
            % put these index vectors into roipoly
            if q >= 1 & q <= z_dim
                roisliceCorners = double(roimaskCorners(:,:,q));
                % find which corners are inside the contour
                BWcorners = roipoly(yCorners,xCorners,roisliceCorners,current_curve(:,2),current_curve(:,1)); 
                % Mark all all pixels bordering corners that are inside the
                % contour
                roi_vox = sum(BWcorners(:));  % number of voxels in this ROI
                % find the voxel overlap between the current roi and BW:
                overlap_vox = sum(sum(BWcorners.*roisliceCorners));
                if overlap_vox == roi_vox
                    roisliceCorners = roisliceCorners - BWcorners;
                else  % if there is not perfect overlap, add the rois
                    roisliceCorners = roisliceCorners + BWcorners;
                end
                roisliceCorners(roisliceCorners > 0) = 1;  % make sure all mask values are unity
                roimaskCorners(:,:,q) = int8(roisliceCorners);  % update the overall mask
            end
        end
        
        % save time be resampling only a subregion
        ind = find(roimaskCorners);
        [I,J,K] = ind2sub([x_dim+1 y_dim+1 z_dim],ind);
        indx = min(I)-3:max(I)+3;
        indy = min(J)-3:max(J)+3;
        indz = min(K)-3:max(K)+3;
        
        indx = indx(indx >= 1 & indx <= x_dim);
        indy = indy(indy >= 1 & indy <= y_dim);
        indz = indz(indz >= 1 & indz <= z_dim);
        
        % convert the corners to a 3-D roi mask
        roimask(indx,indy,indz) = ceil(gridResample3D(xCorners,yCorners,z,roimaskCorners,x(indx),y(indy),z(indz)));
        
        % save the indices of the roi mask
        ROIS{roi_num}.ind = int32(find(roimask ~= 0));
        ROIS{roi_num}.dim = [x_dim y_dim z_dim];
        ROIS{roi_num}.pixdim = [x_pixdim y_pixdim z_pixdim];
        ROIS{roi_num}.start = [x_start y_start z_start];
        ROIS{roi_num}.start_ind = start_ind;
        fprintf('Converted %s vectors to an roi mask.\n',ROIS{roi_num}.name);
    end
    
    % save ROIS to the Geometry structure
    Geometry.ROIS = ROIS;
    clear ROIS;
end

if strcmp(readTrial,'yes')
    % read in the Trial information
    fid_trialfile = fopen(trialfile,'r');
    tline = fgets(fid_trialfile);

    structstack = {};   % stack for determining which field we're inside
    structind = [];    % structure index values for each field
    lines = 1;

    while tline ~= -1   % read through the trial file
        tline = regexprep(tline,'"','''');  % replace " with '
        tline = regexprep(tline,'[','');  % bad characters to remove
        tline = regexprep(tline,']','');  % bad characters to remove
        tline = regexprep(tline,'#','num');
        tline = regexprep(tline,'\\','''');
        % fprintf(tline);
        if length(findstr(tline,'{'))
            % mark the structure we just entered, removing any blank space
            structstack{length(structstack) + 1} = regexprep(strtok(tline,'='),' ','');
            if length(structstack) == 1
                structstack{1} = 'Trial';  % make sure structure name is always the same
            end
            if strcmp(structstack{end},'Beam')  % found a beam structure
                % find the number of existing beams at this location
                eval_line = ['beams = getfield(' structstack{1}];
                for k=2:length(structstack)
                    eval_line = [eval_line ',''' structstack{k} ''''];
                end
                eval_line = [eval_line ');'];
                try
                    eval(eval_line);
                    num_beams = length(beams) + 1;  % this is the current beam number
                catch  % there was an error, so there are no beams yet
                    num_beams = 1;  % we're on the first beam now
                end
            end
        elseif length(findstr(tline,'}'))
            structstack(end) = [];  % step outside of the current field
        else  % the line must be a subfield assignment
            if strcmp(structstack{end},'Points')  % this case treated specially
                % read-in values until the next '}' is reached
                lines = lines + 1;
                mlc_vec = ['['];
                while ~length(findstr(tline,'}'))
                    % build up the vector of MLC positions
                    mlc_vec = [mlc_vec regexprep(tline,' ','')];
                    tline = fgets(fid_trialfile);
                    % fprintf(tline);
                    lines = lines + 1;
                end
                mlc_vec = [mlc_vec ']'];
                % evaluate the points vector
                eval_statement = structstack{1};
                for k=2:length(structstack)
                    if strcmp(structstack{k},'Beam')
                        % include the beam number
                        eval_statement = [eval_statement '.' structstack{k} '(' num2str(num_beams) ')'];
                    else  % non-beam structure, so don't worry about it
                        eval_statement = [eval_statement '.' structstack{k}];
                    end
                end
                eval_statement = [eval_statement ' = ' mlc_vec ';'];
                eval(eval_statement);
                structstack(end) = [];  % step outside of the current Points[] field
            else
                % see if there are any subfields in the expression
                % evaluate the statement
                eval_statement = structstack{1};
                for k=2:length(structstack)
                    if strcmp(structstack{k},'Beam')
                        % include the beam number
                        eval_statement = [eval_statement '.' structstack{k} '(' num2str(num_beams) ')'];
                    else  % non-beam structure, so don't worry about it
                        eval_statement = [eval_statement '.' structstack{k}];
                    end
                end
                eval_statement = [eval_statement '.' tline];
                eval(eval_statement);
            end
        end
        tline = fgets(fid_trialfile);
        lines = lines + 1;
        if isempty(structstack)
            break;   % we're done if we're no longer inside any sets of curly brackets
        end
    end
    fclose(fid_trialfile);

    % dose grid dimensions
    xdim = Trial.DoseGrid.Dimension.X;
    ydim = Trial.DoseGrid.Dimension.Y;
    zdim = Trial.DoseGrid.Dimension.Z;
    siz = [xdim ydim zdim];  % dosegrid size vector

    % vector describing the start location of the dose grids
    xstart = Trial.DoseGrid.Origin.X;
    ystart = Trial.DoseGrid.Origin.Y;
    zstart = Trial.DoseGrid.Origin.Z;
    start = [xstart ystart zstart];  % dosegrid start vector

    % vector describing the size of each voxel
    dx = Trial.DoseGrid.VoxelSize.X;
    dy = Trial.DoseGrid.VoxelSize.Y;
    dz = Trial.DoseGrid.VoxelSize.Z;
    voxel_size = [dx dy dz];

    % Read in the beamlet dose distributions, DRRs, and fluences and store them in
    % Trial.BeamList.Beam(k), where k is the beam number
    for k=0:num_beams-1
        % read-in the dose distribution
        num_vec = num2str(k*5+4);
        % ensure that num_vec has a length of 3
        while length(num_vec) < 3
            num_vec = ['0' num_vec];
        end
        dose_file = [dose_file_base_name '.' num_vec];
        fid = fopen(dose_file,'rb','ieee-be');
        if fid ~= -1
            dose = fread(fid,'float=>float');
            fclose(fid);
            try
                dose = reshape(dose,siz);
                % flip the dose distribution in accordance with the CT image
                for j=1:siz(3)
                    dose(:,:,j) = single(fliplr(double(dose(:,:,j))));
                end
                % save the dose
                Trial = setfield(Trial,'BeamList','Beam',{k+1},'dose',dose);
            catch
                fprintf('Dose for beam %g did not meet specifications so it was not read-in.\n',k+1);
            end
        end

        % read-in the DRR
        % find the size of this drr
        drr_x = Trial.BeamList.Beam(1).FilmImageList.FilmImage.Image.Dimension.X;
        drr_y = Trial.BeamList.Beam(1).FilmImageList.FilmImage.Image.Dimension.Y;
        siz_drr = [drr_x drr_y];

        num_vec = num2str(k*5+3);
        % ensure that num_vec has a length of 3
        while length(num_vec) < 3
            num_vec = ['0' num_vec];
        end
        drr_file = [dose_file_base_name '.' num_vec];
        fid = fopen(drr_file,'rb','ieee-be');
        if fid ~= -1
            drr = fread(fid,'float');
            fclose(fid);
            try  % proceed only if the drr was read-in properly
                drr = reshape(drr,siz_drr);
                % save the DRR
                Trial = setfield(Trial,'BeamList','Beam',{k+1},'FilmImageList','FilmImage','Image','drr',drr);
            catch
                fprintf('DRR for beam %g did not meet specifications so it was not read-in.\n',k+1);
            end
        end

        % read in the fluence map
        num_vec = num2str(k*5+2);
        fluence_siz = [101 101 3];  % all fluence must have this size otherwise and error will be thrown
        % ensure that num_vec has a length of 3
        while length(num_vec) < 3
            num_vec = ['0' num_vec];
        end
        fluence_file = [dose_file_base_name '.' num_vec];
        fid = fopen(fluence_file,'rb','ieee-be');
        if fid ~= -1
            fluence = fread(fid,'float');
            fclose(fid);
            % check to see if fluence matches dose grid size
            if Trial.FluenceGridMatchesDoseGrid ~= 1
                error('Fluence grid does not match dose grid for beam %g, don''t know what to do here.',k+1);
            end
            try  % proceed only if the fluence is read-in properly
                fluence = reshape(fluence,fluence_siz);
                % save the fluence
                Trial = setfield(Trial,'BeamList','Beam',{k+1},'FluenceMap',fluence);
            catch
                fprintf('Fluence map for beam %g did not meet specifications so it was not read-in.\n',k+1);
            end
        end
    end

    % add up all of the doses
    dose_tot = zeros(siz);
    for k=1:num_beams
        if isfield(Trial.BeamList.Beam(k),'dose') & ~isempty(Trial.BeamList.Beam(k).dose)
            dose_tot = dose_tot + double(Trial.BeamList.Beam(k).dose);
        end
    end
    Geometry.Trial = Trial;
    clear Trial;
end

save(save_file,'Geometry');