pinn2matlab.asv 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492
  1. % pinn2matlab.m
  2. %
  3. % Ryan T Flynn 23 August 2007
  4. %
  5. % Reads in raw Pinnacle data to a MATLAB data structure. The patient CT,
  6. % ROIs, and Trial files can all be read in if available and desired. The
  7. % result is a structure called 'Geometry', which contains ROIS and Trial
  8. % fields.
  9. %
  10. % The CT image and all of the ROIs are flipped onto a coordinate system in
  11. % which the position of each voxel index (i,j,k) are given by:
  12. %
  13. % x(i) = x0 + i*delx; (i=0,...,M-1)
  14. % y(j) = y0 + j*delx; (j=0,...,N-1)
  15. % z(k) = z0 + k*delx; (k=0,...,Q-1)
  16. %
  17. % As opposed to the Pinnacle coordinate system, which
  18. %
  19. % x(i) = x0 + i*delx; (i=0,...,M-1)
  20. % y(j) = y0 + (N-1-j)*dely; (j=0,...,N-1)
  21. % z(k) = z0 + k*delz; (k=0,...,Q-1)
  22. %
  23. % Stephen R Bowen April 2009
  24. %
  25. % The MATLAB data structure is called Geometry with the following format:
  26. %
  27. % Geometry.voxel_size (dx dy dz)
  28. % .start (x y z)
  29. % .data (M x N x Q single float matrix)
  30. % .ROIS {1 x Nroi cell array} if readROIS flag on
  31. %
  32. % Within each index of the ROIS cell array, a data structure exists
  33. % containing the vector information from a Pinnacle mesh contour and indices
  34. % of the ROI mask used by WiscPlan:
  35. %
  36. % Geometry.ROIS{1}.name 'string'
  37. % .num_curves int
  38. % .curves {1 x num_curves cell array}
  39. % .ind (int32 data array)
  40. % .dim (M N Q)
  41. % .pix_dim (dx dy dz)
  42. % .start (x y z)
  43. % .start_ind 'string'
  44. %
  45. % This enables seamless translation between Pinnacle and MATLAB
  46. %%%%% Start of user supplied inputs %%%%%
  47. % Read-in flags
  48. readTrial = 'no'; % either 'yes' or 'no'
  49. readROIS = 'yes'; % either 'yes' or 'no'
  50. % locations of the Pinnacle files
  51. roifile = '../plan/HN003/plan.roi';
  52. imageheaderfile = '../plan/HN003/HN003_planningCT.header';
  53. imagefile = '../plan/HN003/HN003_planningCT.img'; % CT file
  54. save_file = '../plan/HN003/HN003.mat';
  55. %%%%% End of user supplied inputs %%%%%
  56. start_ind = '(1,1,1)'; % this tag will be added to the ROIS so it is known where the start voxel is
  57. clear roi ROIS Trial Geometry CTimage; % start with a clean slate of structures
  58. % extract geometric information from the header file
  59. fid_imageheaderfile = fopen(imageheaderfile,'r');
  60. tline = fgets(fid_imageheaderfile);
  61. while tline ~= -1
  62. % check the line for key words
  63. if length(findstr(tline,'byte_order'))
  64. eval(tline);
  65. if length(findstr(tline,'x_dim')) & ~length(findstr(tline,'fname_index_start'))
  66. eval(tline); % run the line to get x_dim
  67. elseif length(findstr(tline,'y_dim'))
  68. eval(tline); % run the line to get y_dim
  69. elseif length(findstr(tline,'z_dim'))
  70. eval(tline); % run the line to get z_dim
  71. elseif length(findstr(tline,'x_pixdim'))
  72. eval(tline); % run the line to get x_pixdim
  73. elseif length(findstr(tline,'y_pixdim'))
  74. eval(tline); % run the line to get y_pixdim
  75. elseif length(findstr(tline,'z_pixdim'))
  76. eval(tline); % run the line to get z_pixdim
  77. elseif length(findstr(tline,'x_start'))
  78. eval(tline); % run the line to get x_start
  79. elseif length(findstr(tline,'y_start'))
  80. eval(tline); % run the line to get x_start
  81. elseif length(findstr(tline,'z_start'))
  82. eval(tline); % run the line to get x_start
  83. end
  84. tline = fgets(fid_imageheaderfile);
  85. end
  86. fclose(fid_imageheaderfile);
  87. % Read in the CT data
  88. fid_image = fopen(imagefile,'rb');
  89. if byte_order == 0
  90. CTimage = fread(fid_image,'short=>int16');
  91. elseif byte_order == 1
  92. CTimage = fread(fid_image,'short=>int16','ieee-be');
  93. end
  94. CTimage = reshape(CTimage,x_dim,y_dim,z_dim);
  95. fclose(fid_image);
  96. % flip the CT image onto the non-Pinnacle coordinate system
  97. CTimage = flipdim(CTimage,2);
  98. % save the results to the Geometry structure
  99. Geometry.start = single([x_start y_start z_start]);
  100. Geometry.voxel_size = single([x_pixdim y_pixdim z_pixdim]);
  101. Geometry.data = single(CTimage)/1024; % convert the CT data to density data, Pinnacle style
  102. Geometry.data(Geometry.data < 0) = 0; % truncate anything less than zero
  103. clear CTimage;
  104. if strcmp(readROIS,'yes')
  105. % read in the roi file
  106. fid_roifile = fopen(roifile,'r');
  107. roinames = {}; % start a cell array of the roi names
  108. Nroi = 0; % number of rois read in
  109. % Flags to indicate which sets of angled brackets in the roi file tline is
  110. % inside.
  111. inroi = 0;
  112. incurve = 0;
  113. inpoints = 0;
  114. roi_num = 0; % current roi number
  115. tline = fgets(fid_roifile);
  116. while tline ~= -1
  117. % check the line for key words
  118. if length(findstr(tline,'roi={'))
  119. inroi = 1; % mark that we are now currently inside of an roi
  120. roi_num = roi_num + 1;
  121. % next line contains the roi name
  122. tline = fgets(fid_roifile);
  123. % pop off first token in line, the remainder of the line is the roi name
  124. [T,R] = strtok(tline);
  125. roi.name = strtrim(R);
  126. % pop off lines until we get to the number of curves in this roi
  127. while ~length(findstr(tline,'num_curve'))
  128. tline = fgets(fid_roifile);
  129. end
  130. % pop off the num_curve string
  131. [T,R] = strtok(tline);
  132. % pop off the equals sign
  133. [T,R] = strtok(R);
  134. % pop off the number of curves in this roi
  135. T = strtok(R,';');
  136. % save the number of curves to the roi stucture
  137. eval(['roi.num_curves = ' num2str(T) ';']);
  138. roi.curves = {}; % get the curves structure started
  139. % read in the next curve structure
  140. curve_num = 0; % number of the current curve
  141. while roi.num_curves > 0 & curve_num < roi.num_curves
  142. while ~length(findstr(tline,'curve={'));
  143. tline = fgets(fid_roifile);
  144. end
  145. curve_num = curve_num + 1;
  146. incurve = 1; % inside the curve structure now
  147. % find the number of points in this structure
  148. while ~length(findstr(tline,'num_points'))
  149. tline = fgets(fid_roifile);
  150. end
  151. % pop off the num_points string
  152. [T,R] = strtok(tline);
  153. % pop off the equals sign
  154. [T,R] = strtok(R);
  155. % pop off the number of points in this curve
  156. T = strtok(R,';');
  157. eval(['num_points = ' num2str(T) ';']);
  158. % find the points identifier
  159. while ~length(findstr(tline,'points={'))
  160. tline = fgets(fid_roifile);
  161. end
  162. inpoints = 1; % inside the points structure now
  163. % read in the block of points data
  164. block = fscanf(fid_roifile,'%g',[3 num_points]);
  165. % save the block of points to the roi stucture
  166. roi.curves{curve_num} = block';
  167. % read in the right parantheses for the points and curve
  168. % structures
  169. while ~length(findstr(tline,'};'))
  170. tline = fgets(fid_roifile);
  171. end
  172. inpoints = 0; % stepped outside of the points structure
  173. while ~length(findstr(tline,'};'))
  174. tline = fgets(fid_roifile);
  175. end
  176. incurve = 0; % stepped outside of the curves structure
  177. end
  178. % read until the roi closing bracket appears
  179. while ~length(findstr(tline,'};'))
  180. tline = fgets(fid_roifile);
  181. end
  182. inroi = 0; % we are now outside of the roi
  183. fprintf('read in %s roi\n',roi.name);
  184. ROIS{roi_num} = roi;
  185. end
  186. tline = fgets(fid_roifile);
  187. end
  188. fclose(fid_roifile);
  189. % ROIS have all been read-in, now just have to convert them to mask
  190. % matrices. We'll use roipoly for this.
  191. % Recall that we must use the Pinnacle coordinate system for this to work,
  192. % that is, (x,y,z) coordinates are given by:
  193. % x(i) = x_start + i*x_pixdim, i=0..x_dim-1
  194. % y(j) = y_start + (y_dim-1)*y_pixdim - j*y_pixdim, j=0..y_dim-1
  195. % z(k) = z_start + k*z_pixdim, k=0..z_dim-1
  196. %
  197. % Not all treatment planning systems use this type of coordinate system
  198. % definition, so it is very important to get them straight.
  199. %
  200. % To get around this we will manipulate the data such that we'll have:
  201. %
  202. % x(i) = x_start + i*x_pixdim, i=0..x_dim-1
  203. % y(j) = y_start + j*y_pixdim, j=0..y_dim-1
  204. % z(k) = z_start + k*z_pixdim, k=0..z_dim-1
  205. %
  206. % This can be done by a simple fliplr operation on all of the CT slices
  207. % define the coordinate system
  208. x = x_start:x_pixdim:x_start + (x_dim-1)*x_pixdim;
  209. y = y_start:y_pixdim:y_start + (y_dim-1)*y_pixdim;
  210. z = z_start:z_pixdim:z_start + (z_dim-1)*z_pixdim;
  211. % define the locations of the corners of the pixels in each slice
  212. xCorners = [x - x_pixdim/2 x(end) + x_pixdim/2];
  213. yCorners = [y - y_pixdim/2 y(end) + y_pixdim/2];
  214. % loop through all rois
  215. for roi_num=1:length(ROIS)
  216. % set up the roi mask
  217. roimask = zeros(x_dim,y_dim,z_dim,'int8');
  218. roimaskCorners = zeros(x_dim+1,y_dim+1,z_dim,'int8');
  219. % loop through all of the curves in the roi
  220. for curve_num=1:length(ROIS{roi_num}.curves)
  221. % make a copy of the curve for easy access
  222. current_curve = ROIS{roi_num}.curves{curve_num};
  223. % find the z-index (slice number) of this curve
  224. q = round((current_curve(1,3)-z_start)/z_pixdim) + 1;
  225. % put these index vectors into roipoly
  226. if q >= 1 & q <= z_dim
  227. roisliceCorners = double(roimaskCorners(:,:,q));
  228. % find which corners are inside the contour
  229. BWcorners = roipoly(yCorners,xCorners,roisliceCorners,current_curve(:,2),current_curve(:,1));
  230. % Mark all all pixels bordering corners that are inside the
  231. % contour
  232. roi_vox = sum(BWcorners(:)); % number of voxels in this ROI
  233. % find the voxel overlap between the current roi and BW:
  234. overlap_vox = sum(sum(BWcorners.*roisliceCorners));
  235. if overlap_vox == roi_vox
  236. roisliceCorners = roisliceCorners - BWcorners;
  237. else % if there is not perfect overlap, add the rois
  238. roisliceCorners = roisliceCorners + BWcorners;
  239. end
  240. roisliceCorners(roisliceCorners > 0) = 1; % make sure all mask values are unity
  241. roimaskCorners(:,:,q) = int8(roisliceCorners); % update the overall mask
  242. end
  243. end
  244. % save time be resampling only a subregion
  245. ind = find(roimaskCorners);
  246. [I,J,K] = ind2sub([x_dim+1 y_dim+1 z_dim],ind);
  247. indx = min(I)-3:max(I)+3;
  248. indy = min(J)-3:max(J)+3;
  249. indz = min(K)-3:max(K)+3;
  250. indx = indx(indx >= 1 & indx <= x_dim);
  251. indy = indy(indy >= 1 & indy <= y_dim);
  252. indz = indz(indz >= 1 & indz <= z_dim);
  253. % convert the corners to a 3-D roi mask
  254. roimask(indx,indy,indz) = ceil(gridResample3D(xCorners,yCorners,z,roimaskCorners,x(indx),y(indy),z(indz)));
  255. % save the indices of the roi mask
  256. ROIS{roi_num}.ind = int32(find(roimask ~= 0));
  257. ROIS{roi_num}.dim = [x_dim y_dim z_dim];
  258. ROIS{roi_num}.pixdim = [x_pixdim y_pixdim z_pixdim];
  259. ROIS{roi_num}.start = [x_start y_start z_start];
  260. ROIS{roi_num}.start_ind = start_ind;
  261. fprintf('Converted %s vectors to an roi mask.\n',ROIS{roi_num}.name);
  262. end
  263. % save ROIS to the Geometry structure
  264. Geometry.ROIS = ROIS;
  265. clear ROIS;
  266. end
  267. if strcmp(readTrial,'yes')
  268. % read in the Trial information
  269. fid_trialfile = fopen(trialfile,'r');
  270. tline = fgets(fid_trialfile);
  271. structstack = {}; % stack for determining which field we're inside
  272. structind = []; % structure index values for each field
  273. lines = 1;
  274. while tline ~= -1 % read through the trial file
  275. tline = regexprep(tline,'"',''''); % replace " with '
  276. tline = regexprep(tline,'[',''); % bad characters to remove
  277. tline = regexprep(tline,']',''); % bad characters to remove
  278. tline = regexprep(tline,'#','num');
  279. tline = regexprep(tline,'\\','''');
  280. % fprintf(tline);
  281. if length(findstr(tline,'{'))
  282. % mark the structure we just entered, removing any blank space
  283. structstack{length(structstack) + 1} = regexprep(strtok(tline,'='),' ','');
  284. if length(structstack) == 1
  285. structstack{1} = 'Trial'; % make sure structure name is always the same
  286. end
  287. if strcmp(structstack{end},'Beam') % found a beam structure
  288. % find the number of existing beams at this location
  289. eval_line = ['beams = getfield(' structstack{1}];
  290. for k=2:length(structstack)
  291. eval_line = [eval_line ',''' structstack{k} ''''];
  292. end
  293. eval_line = [eval_line ');'];
  294. try
  295. eval(eval_line);
  296. num_beams = length(beams) + 1; % this is the current beam number
  297. catch % there was an error, so there are no beams yet
  298. num_beams = 1; % we're on the first beam now
  299. end
  300. end
  301. elseif length(findstr(tline,'}'))
  302. structstack(end) = []; % step outside of the current field
  303. else % the line must be a subfield assignment
  304. if strcmp(structstack{end},'Points') % this case treated specially
  305. % read-in values until the next '}' is reached
  306. lines = lines + 1;
  307. mlc_vec = ['['];
  308. while ~length(findstr(tline,'}'))
  309. % build up the vector of MLC positions
  310. mlc_vec = [mlc_vec regexprep(tline,' ','')];
  311. tline = fgets(fid_trialfile);
  312. % fprintf(tline);
  313. lines = lines + 1;
  314. end
  315. mlc_vec = [mlc_vec ']'];
  316. % evaluate the points vector
  317. eval_statement = structstack{1};
  318. for k=2:length(structstack)
  319. if strcmp(structstack{k},'Beam')
  320. % include the beam number
  321. eval_statement = [eval_statement '.' structstack{k} '(' num2str(num_beams) ')'];
  322. else % non-beam structure, so don't worry about it
  323. eval_statement = [eval_statement '.' structstack{k}];
  324. end
  325. end
  326. eval_statement = [eval_statement ' = ' mlc_vec ';'];
  327. eval(eval_statement);
  328. structstack(end) = []; % step outside of the current Points[] field
  329. else
  330. % see if there are any subfields in the expression
  331. % evaluate the statement
  332. eval_statement = structstack{1};
  333. for k=2:length(structstack)
  334. if strcmp(structstack{k},'Beam')
  335. % include the beam number
  336. eval_statement = [eval_statement '.' structstack{k} '(' num2str(num_beams) ')'];
  337. else % non-beam structure, so don't worry about it
  338. eval_statement = [eval_statement '.' structstack{k}];
  339. end
  340. end
  341. eval_statement = [eval_statement '.' tline];
  342. eval(eval_statement);
  343. end
  344. end
  345. tline = fgets(fid_trialfile);
  346. lines = lines + 1;
  347. if isempty(structstack)
  348. break; % we're done if we're no longer inside any sets of curly brackets
  349. end
  350. end
  351. fclose(fid_trialfile);
  352. % dose grid dimensions
  353. xdim = Trial.DoseGrid.Dimension.X;
  354. ydim = Trial.DoseGrid.Dimension.Y;
  355. zdim = Trial.DoseGrid.Dimension.Z;
  356. siz = [xdim ydim zdim]; % dosegrid size vector
  357. % vector describing the start location of the dose grids
  358. xstart = Trial.DoseGrid.Origin.X;
  359. ystart = Trial.DoseGrid.Origin.Y;
  360. zstart = Trial.DoseGrid.Origin.Z;
  361. start = [xstart ystart zstart]; % dosegrid start vector
  362. % vector describing the size of each voxel
  363. dx = Trial.DoseGrid.VoxelSize.X;
  364. dy = Trial.DoseGrid.VoxelSize.Y;
  365. dz = Trial.DoseGrid.VoxelSize.Z;
  366. voxel_size = [dx dy dz];
  367. % Read in the beamlet dose distributions, DRRs, and fluences and store them in
  368. % Trial.BeamList.Beam(k), where k is the beam number
  369. for k=0:num_beams-1
  370. % read-in the dose distribution
  371. num_vec = num2str(k*5+4);
  372. % ensure that num_vec has a length of 3
  373. while length(num_vec) < 3
  374. num_vec = ['0' num_vec];
  375. end
  376. dose_file = [dose_file_base_name '.' num_vec];
  377. fid = fopen(dose_file,'rb','ieee-be');
  378. if fid ~= -1
  379. dose = fread(fid,'float=>float');
  380. fclose(fid);
  381. try
  382. dose = reshape(dose,siz);
  383. % flip the dose distribution in accordance with the CT image
  384. for j=1:siz(3)
  385. dose(:,:,j) = single(fliplr(double(dose(:,:,j))));
  386. end
  387. % save the dose
  388. Trial = setfield(Trial,'BeamList','Beam',{k+1},'dose',dose);
  389. catch
  390. fprintf('Dose for beam %g did not meet specifications so it was not read-in.\n',k+1);
  391. end
  392. end
  393. % read-in the DRR
  394. % find the size of this drr
  395. drr_x = Trial.BeamList.Beam(1).FilmImageList.FilmImage.Image.Dimension.X;
  396. drr_y = Trial.BeamList.Beam(1).FilmImageList.FilmImage.Image.Dimension.Y;
  397. siz_drr = [drr_x drr_y];
  398. num_vec = num2str(k*5+3);
  399. % ensure that num_vec has a length of 3
  400. while length(num_vec) < 3
  401. num_vec = ['0' num_vec];
  402. end
  403. drr_file = [dose_file_base_name '.' num_vec];
  404. fid = fopen(drr_file,'rb','ieee-be');
  405. if fid ~= -1
  406. drr = fread(fid,'float');
  407. fclose(fid);
  408. try % proceed only if the drr was read-in properly
  409. drr = reshape(drr,siz_drr);
  410. % save the DRR
  411. Trial = setfield(Trial,'BeamList','Beam',{k+1},'FilmImageList','FilmImage','Image','drr',drr);
  412. catch
  413. fprintf('DRR for beam %g did not meet specifications so it was not read-in.\n',k+1);
  414. end
  415. end
  416. % read in the fluence map
  417. num_vec = num2str(k*5+2);
  418. fluence_siz = [101 101 3]; % all fluence must have this size otherwise and error will be thrown
  419. % ensure that num_vec has a length of 3
  420. while length(num_vec) < 3
  421. num_vec = ['0' num_vec];
  422. end
  423. fluence_file = [dose_file_base_name '.' num_vec];
  424. fid = fopen(fluence_file,'rb','ieee-be');
  425. if fid ~= -1
  426. fluence = fread(fid,'float');
  427. fclose(fid);
  428. % check to see if fluence matches dose grid size
  429. if Trial.FluenceGridMatchesDoseGrid ~= 1
  430. error('Fluence grid does not match dose grid for beam %g, don''t know what to do here.',k+1);
  431. end
  432. try % proceed only if the fluence is read-in properly
  433. fluence = reshape(fluence,fluence_siz);
  434. % save the fluence
  435. Trial = setfield(Trial,'BeamList','Beam',{k+1},'FluenceMap',fluence);
  436. catch
  437. fprintf('Fluence map for beam %g did not meet specifications so it was not read-in.\n',k+1);
  438. end
  439. end
  440. end
  441. % add up all of the doses
  442. dose_tot = zeros(siz);
  443. for k=1:num_beams
  444. if isfield(Trial.BeamList.Beam(k),'dose') & ~isempty(Trial.BeamList.Beam(k).dose)
  445. dose_tot = dose_tot + double(Trial.BeamList.Beam(k).dose);
  446. end
  447. end
  448. Geometry.Trial = Trial;
  449. clear Trial;
  450. end
  451. save(save_file,'Geometry');