%% 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) % % Stephen R Bowen April 2009 % % The MATLAB structure Geometry has the following format: % % Geometry.voxel_size (dx dy dz) % .start (x y z) % .data (M x N x Q single float matrix) % .ROIS {1 x Nroi cell array} if readROIS flag on % % Within each index of the ROIS cell array, a data structure exists % containing the vector information from each Pinnacle mesh contour and indices % of the corresponding ROI mask used by WiscPlan: % % Geometry.ROIS{1}.name 'string' % .num_curves int % .curves {1 x num_curves cell array} % .ind (int32 data array) % .dim (M N Q) % .pix_dim (dx dy dz) % .start (x y z) % .start_ind 'string' % % This enables seamless translation between Pinnacle and MATLAB. %% %%%%% 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 = '../plan/HN003/plan.roi'; imageheaderfile = '../plan/HN003/HN003_planningCT.header'; imagefile = '../plan/HN003/HN003_planningCT.img'; % CT file save_file = '../plan/HN003/HN003.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,'byte_order')) eval(tline); elseif 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'); if byte_order == 0 CTimage = fread(fid_image,'short=>int16'); elseif byte_order == 1 CTimage = fread(fid_image,'short=>int16','ieee-be'); end 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');