% pinn2amira.m
%
% Ryan T Flynn 23 August 2007
%
% Reads in raw Pinnacle data and creates ROI mask files and an hx file that
% is readable by Amira.
%
% 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-defined parameters %%%%%

% locations of the Pinnacle files, getting rid of backslashes
imageheaderfile = regexprep('C:\data_temp\steveTPS\plan\HN\ImageSet_0.header','\\','/');
imagefile = regexprep('C:\Documents and Settings\Steve\My Documents\Research\treatmentPlanningSystem\patients\hn003\PinnacleData\ImageSet_0.img','\\','/');
roifile = regexprep('C:\Documents and Settings\Steve\My Documents\Research\treatmentPlanningSystem\patients\hn003\PinnacleData\Plan_0\plan.roi','\\','/');

amiraFolder = regexprep('C:\Documents and Settings\Steve\My Documents\Research\treatmentPlanningSystem\patients\hn003\AmiraData','\\','/');  % folder where Amira data will be saved

amiraImgDataName = 'CT';  % name to give the image data in Amira

amiraNetworkFileName = 'headAndNeck.hx';  % name of resulting Amira network file

%%%%% End of user-defined parameters

mkdir(amiraFolder);  % create Amira data folder if it doesn't exist already

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

% 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
    elseif length(findstr(tline,'bytes_pix'))
        eval(tline);  % number of bytes per pixel
    end
    tline = fgets(fid_imageheaderfile);
end

fclose(fid_imageheaderfile);

% calculate axes for the image
x = x_start + [0:x_dim-1]*x_pixdim;
y = y_start + [0:y_dim-1]*y_pixdim;
z = z_start + [0:z_dim-1]*z_pixdim;

if bytes_pix == 2
    dataType = 'short 1';
end

% create Amira script
fidAmira = fopen([amiraFolder '/' amiraNetworkFileName],'w');
fprintf(fidAmira,'# Amira Script\n');
fprintf(fidAmira,'remove -all\n');
fprintf(fidAmira,'remove ImageSet_0.img\n');
fprintf(fidAmira,'\n');
fprintf(fidAmira,'set hideNewModules 0\n');
fprintf(fidAmira,['[ load -raw ' imagefile ' little xfastest ' dataType ' ' num2str(x_dim) ' ' num2str(y_dim) ' ' num2str(z_dim) ...
     ' ' num2str(x(1)) ' ' num2str(x(end)) ' ' num2str(y(1)) ' ' num2str(y(end)) ' ' num2str(z(1)) ' ' num2str(z(end)) ...
     '] setLabel ' amiraImgDataName '\n']);
fprintf(fidAmira,[amiraImgDataName ' setIconPosition 20 10\n']);
fprintf(fidAmira,[amiraImgDataName ' flip 1\n']);  % flip the data along the y-axis
fprintf(fidAmira,[amiraImgDataName ' fire\n']);
fprintf(fidAmira,[amiraImgDataName ' setViewerMask 65535\n']);
fprintf(fidAmira,'set hideNewModules 0\n');

% done with image dataslice

% read in ROIS
fid_roifile = fopen(roifile,'r');

roinames = {};  % start a cell array of the roi names

% 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;

Nroi = length(ROIS);  % get number of ROIS

% 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
% 
% % loop through all rois
% for roi_num=1:Nroi
%     % set up the roi mask
%     roimask = zeros(x_dim,y_dim,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;
%         % convert the x and y vectors to indices
%         xind = (current_curve(:,1) - x_start)/x_pixdim + 1;
%         yind = (current_curve(:,2) - y_start)/y_pixdim + 1;
%         % put these index vectors into roipoly
%         if q >= 1 & q <= z_dim
%             roislice = double(roimask(:,:,q));
%             BW = roipoly(roislice,yind,xind);
%             roi_vox = sum(BW(:));  % number of voxels in this ROI
%             % find the voxel overlap between the current roi and BW:
%             overlap_vox = sum(sum(BW.*roislice));
%             if overlap_vox == roi_vox
%                 roislice = roislice - BW;
%             else  % if there is not perfect overlap, add the rois
%                 roislice = roislice + BW;
%             end
%             roislice(roislice > 0) = 1;  % make sure all mask values are unity
%             roimask(:,:,q) = int8(roislice);
%         end      
%     end
%     % 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

% now that we have the ROIS in mask form, we can save them into an
% Amira-readable format.

% icon position vectors
xIcon = [0:2]*200 + 50;
yIcon = [0:1]*50;
yShift = [0:Nroi]*60 + 400;

% get the names of the rois
roi_names = cell(1,Nroi);
for k=1:Nroi
    roi_names{k} = ROIS{k}.name;
end

mask.start = [x(1) y(1) z(1)];
mask.voxel_size = [x_pixdim y_pixdim z_pixdim];

% create files for each mask region
for k=1:Nroi
    mask_filename = [regexprep(roi_names{k},{'/',' '},'') '.am'];
    mask_filename = regexprep(mask_filename,' ','');  % remove spaces
    mask.data = zeros([x_dim y_dim z_dim],'int8');
    mask.data(ROIS{k}.ind) = 1;  % fill in the mask
    geom2am(mask,[amiraFolder '/' mask_filename]);  % save the structure to an Amira file
    
    fprintf(fidAmira,'set hideNewModules 0\n');
    fprintf(fidAmira,['[ load ${SCRIPTDIR}/' mask_filename ' ] setLabel ' mask_filename '\n']);
    fprintf(fidAmira,[mask_filename ' setIconPosition ' num2str(xIcon(1)) ' ' num2str(yIcon(1)+yShift(k)) '\n']);
    fprintf(fidAmira,[mask_filename ' fire\n']);
    fprintf(fidAmira,[mask_filename ' fire\n']);
    fprintf(fidAmira,[mask_filename ' setViewerMask 65535\n']);
    fprintf(fidAmira,'\n');
    
    % set up a CastField module for a LabelField for the mask
    fprintf(fidAmira,'set hideNewModules 0\n');
    castFieldName = ['CastField' num2str(k)];
    fprintf(fidAmira,['create HxCastField {' castFieldName '}\n']);
    fprintf(fidAmira,[castFieldName ' setIconPosition ' num2str(xIcon(2)) ' ' num2str(yIcon(1)+yShift(k)) '\n']);
    fprintf(fidAmira,[castFieldName ' data connect ' mask_filename '\n']);
    fprintf(fidAmira,[castFieldName ' colormap setDefaultColor 1 0.8 0.5\n']);
    fprintf(fidAmira,[castFieldName ' colormap setDefaultAlpha 0.500000\n']);
    fprintf(fidAmira,[castFieldName ' fire\n']);
    fprintf(fidAmira,[castFieldName ' outputType setValue 0 6\n']);
    fprintf(fidAmira,[castFieldName ' scaling setMinMax 0 -1.00000001384843e+024 1.00000001384843e+024\n']);
    fprintf(fidAmira,[castFieldName ' scaling setValue 0 1\n']);
    fprintf(fidAmira,[castFieldName ' scaling setMinMax 1 -1.00000001384843e+024 1.00000001384843e+024\n']);
    fprintf(fidAmira,[castFieldName ' scaling setValue 1 0\n']);
    fprintf(fidAmira,[castFieldName ' voxelGridOptions setValue 0 1\n']);
    fprintf(fidAmira,[castFieldName ' colorFieldOptions setValue 0 0\n']);
    fprintf(fidAmira,[castFieldName ' fire\n']);
    fprintf(fidAmira,[castFieldName ' setViewerMask 65535\n']);
    fprintf(fidAmira,[castFieldName ' deselect\n']);
    fprintf(fidAmira,'\n');
    
    % set up mask plotting routines
    mask_labelfield_name = regexprep(mask_filename,'am','Labelfield'); 
    fprintf(fidAmira,'set hideNewModules 0\n');
    fprintf(fidAmira,['[ {' castFieldName '} create\n']);
    fprintf(fidAmira,['    ] setLabel {' mask_labelfield_name '}\n']);
    fprintf(fidAmira,[mask_labelfield_name ' setIconPosition ' num2str(xIcon(1)) ' ' num2str(yIcon(2)+yShift(k)) '\n']);
    fprintf(fidAmira,[mask_labelfield_name ' master connect ' castFieldName '\n']);
    fprintf(fidAmira,[mask_labelfield_name ' fire\n']);
    fprintf(fidAmira,[mask_labelfield_name ' primary setValue 0 0\n']);
    fprintf(fidAmira,[mask_labelfield_name ' fire\n']);
    fprintf(fidAmira,[mask_labelfield_name ' setViewerMask 65535\n']);
    fprintf(fidAmira,'\n');
    
    surfaceGenName = ['SurfaceGen' num2str(k)];
    fprintf(fidAmira,'set hideNewModules 0\n');
    fprintf(fidAmira,['create HxGMC {' surfaceGenName '}\n']);
    fprintf(fidAmira,[surfaceGenName ' setIconPosition ' num2str(xIcon(3)) ' ' num2str(yIcon(1)+yShift(k)) '\n']);
    fprintf(fidAmira,[surfaceGenName ' data connect ' mask_labelfield_name '\n']);
    fprintf(fidAmira,[surfaceGenName ' fire\n']);
    fprintf(fidAmira,[surfaceGenName ' smoothing setValue 0 3\n']);
    fprintf(fidAmira,[surfaceGenName ' options setValue 0 1\n']);
    fprintf(fidAmira,[surfaceGenName ' options setValue 1 0\n']);
    fprintf(fidAmira,[surfaceGenName ' border setValue 0 1\n']);
    fprintf(fidAmira,[surfaceGenName ' border setValue 1 0\n']);
    fprintf(fidAmira,[surfaceGenName ' minEdgeLength setMinMax 0 0 0.800000011920929\n']);
    fprintf(fidAmira,[surfaceGenName ' minEdgeLength setValue 0 0\n']);
    fprintf(fidAmira,[surfaceGenName ' materialList setValue 0 0\n']);
    fprintf(fidAmira,[surfaceGenName ' fire\n']);
    fprintf(fidAmira,[surfaceGenName ' setViewerMask 65535\n']);
    fprintf(fidAmira,'\n');
    
    mask_surf_name = regexprep(mask_filename,'am','surf'); 
    fprintf(fidAmira,'set hideNewModules 0\n');
    fprintf(fidAmira,['[ {' surfaceGenName '} create {' mask_surf_name '}\n']);
    fprintf(fidAmira,['    ] setLabel {' mask_surf_name '}\n']);
    fprintf(fidAmira,[mask_surf_name ' setIconPosition ' num2str(xIcon(2)) ' ' num2str(yIcon(2)+yShift(k)) '\n']);
    fprintf(fidAmira,[mask_surf_name ' master connect ' surfaceGenName '\n']);
    fprintf(fidAmira,[mask_surf_name ' fire\n']);
    fprintf(fidAmira,[mask_surf_name ' LevelOfDetail setMinMax -1 -1\n']);
    fprintf(fidAmira,[mask_surf_name ' LevelOfDetail setButtons 1\n']);
    fprintf(fidAmira,[mask_surf_name ' LevelOfDetail setIncrement 1\n']);
    fprintf(fidAmira,[mask_surf_name ' LevelOfDetail setValue -1\n']);
    fprintf(fidAmira,[mask_surf_name ' LevelOfDetail setSubMinMax -1 -1\n']);
    fprintf(fidAmira,[mask_surf_name ' fire\n']);
    fprintf(fidAmira,[mask_surf_name ' setViewerMask 65535\n']);
    fprintf(fidAmira,[mask_surf_name ' deselect\n']);
    fprintf(fidAmira,'\n');
    fprintf(fidAmira,'set hideNewModules 0\n');
    fprintf(fidAmira,'\n');
    
    surfaceViewName = ['SurfaceView' num2str(k)];
    fprintf(fidAmira,'set hideNewModules 0\n');
    fprintf(fidAmira,['create HxDisplaySurface {' surfaceViewName '}\n']);
    fprintf(fidAmira,[surfaceViewName ' setIconPosition ' num2str(xIcon(3)) ' ' num2str(yIcon(2)+yShift(k)) '\n']);
    fprintf(fidAmira,[surfaceViewName ' data connect ' mask_surf_name '\n']);
    fprintf(fidAmira,[surfaceViewName ' colormap setDefaultColor 1 0.1 0.1\n']);
    fprintf(fidAmira,[surfaceViewName ' colormap setDefaultAlpha 0.500000\n']);
    fprintf(fidAmira,[surfaceViewName ' fire\n']);
    fprintf(fidAmira,[surfaceViewName '  drawStyle setValue 4\n']);
    fprintf(fidAmira,[surfaceViewName '  drawStyle setSpecularLighting 1\n']);
    fprintf(fidAmira,[surfaceViewName '  drawStyle setTexture 0\n']);
    fprintf(fidAmira,[surfaceViewName '  drawStyle setAlphaMode 3\n']);
    fprintf(fidAmira,[surfaceViewName '  drawStyle setNormalBinding 0\n']);
    fprintf(fidAmira,[surfaceViewName '  drawStyle setCullingMode 0\n']);
    fprintf(fidAmira,[surfaceViewName '  selectionMode setValue 0 0\n']);
    fprintf(fidAmira,[surfaceViewName '  Patch setMinMax 0 1\n']);
    fprintf(fidAmira,[surfaceViewName '  Patch setButtons 1\n']);
    fprintf(fidAmira,[surfaceViewName '  Patch setIncrement 1\n']);
    fprintf(fidAmira,[surfaceViewName '  Patch setValue 1\n']);
    fprintf(fidAmira,[surfaceViewName '  Patch setSubMinMax 0 1\n']);
    fprintf(fidAmira,[surfaceViewName '  BoundaryId setValue 0 -1\n']);
    fprintf(fidAmira,[surfaceViewName '  materials setValue 0 1\n']);
    fprintf(fidAmira,[surfaceViewName '  materials setValue 1 0\n']);
    fprintf(fidAmira,[surfaceViewName '  colorMode setValue 0\n']);
    fprintf(fidAmira,[surfaceViewName '  baseTrans setMinMax 0 1\n']);
    fprintf(fidAmira,[surfaceViewName '  baseTrans setButtons 0\n']);
    fprintf(fidAmira,[surfaceViewName '  baseTrans setIncrement 0.1\n']);
    fprintf(fidAmira,[surfaceViewName '  baseTrans setValue 0.53\n']);
    fprintf(fidAmira,[surfaceViewName '  baseTrans setSubMinMax 0 1\n']);
    fprintf(fidAmira,[surfaceViewName '  VRMode setValue 0 0\n']);
    fprintf(fidAmira,[surfaceViewName '  fire\n']);
    fprintf(fidAmira,[surfaceViewName '  hideBox 1\n']);
    fprintf(fidAmira,['{' surfaceViewName '} selectTriangles zab HIJMPLPPHPGAABDPBAADAACGHIICIN\n']);
    fprintf(fidAmira,[surfaceViewName '  fire\n']);
    fprintf(fidAmira,[surfaceViewName '  setViewerMask 65535\n']);
    fprintf(fidAmira,[surfaceViewName '  deselect\n']);
    fprintf(fidAmira,'\n');
    
    % remove unneeded junk
    fprintf(fidAmira,['remove ' mask_filename '\n']);
    fprintf(fidAmira,['remove ' surfaceGenName '\n']);
end

fclose(fidAmira);