1- %% Example: Photon Treatment Plan
1+ %% Example: Comparison of a CT and (fake) synthetic CT dose calculation
22%
33% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44%
2828% from .mat format to .dcm, consisting of 3D body volume slices as well as a structure file and use it as a real CT.
2929% Additionally, we will slightly modify the intensities of the phantom to generate
3030% synthetic CT images, illustrating dose difference.
31- clear ;
3231matRad_cfg = matRad_rc ; % If this throws an error, run it from the parent directory first to set the paths
3332
3433% Create directories to store DICOM real CT and synthetic (fake) CT data
35- patDir= " userdata\syntheticCT\" ; % If you want to export your data, use "userdata" folder as it is ignored by git
36- name= " pat001" ;
37- patDirRealCT= fullfile(patDir ,name ," realCT" );
38- patDirFakeCT= fullfile(patDir ,name ," fakeCT" );
39-
40- try
41- if ~exist(patDirRealCT , ' dir' )
42- mkdir(patDirRealCT );
43- end
44- catch ME
45- fprintf(' Error creating realCT directory: %s\n ' , patDirRealCT );
46- fprintf(' Error message: %s\n ' , ME .message );
47- end
34+ patDir = [matRad_cfg .primaryUserFolder filesep ' syntheticCT' filesep ]; % If you want to export your data, use "userdata" folder as it is ignored by git
35+ name = ' LIVER' ;
36+ patDirRealCT = [name ' _realCT' ];
37+ patDirFakeCT = [name ' _fakeCT' ];
4838
49- try
50- if ~exist(patDirFakeCT , ' dir' )
51- mkdir(patDirFakeCT );
52- end
53- catch ME
54- fprintf(' Error creating fakeCT directory: %s\n ' , patDirFakeCT );
55- fprintf(' Error message: %s\n ' , ME .message );
39+ if ~mkdir(patDir ,patDirRealCT ) || ~mkdir(patDir ,patDirFakeCT )
40+ matRad_cfg .dispError(' Error creating directories %s and %s in %s ' ,patDirRealCT ,patDirFakeCT ,patDir );
5641end
5742
43+ patDirRealCT = fullfile(patDir ,patDirRealCT );
44+ patDirFakeCT = fullfile(patDir ,patDirFakeCT );
45+
5846% Now, as the directories are created, let us create the DICOM data. Here, the DICOM
5947% files will be created from the .mat format.
6048load(' LIVER.mat' );
61- realCTct= ct ;
49+ indicators = {' mean' ,' std' ,' max' , ' min' , ' D_2' ,' D_5' , ' D_95' , ' D_98' };
50+ realCTct = ct ;
51+ realCTcst = cst ;
52+
6253% review real CT volume
6354matRadGUI ;
6455% saving as DICOMs
65- dcmExpRealCT= matRad_DicomExporter ; % create instance of matRad_DicomExporter
56+ dcmExpRealCT = matRad_DicomExporter ; % create instance of matRad_DicomExporter
6657dcmExpRealCT.dicomDir = patDirRealCT ; % set the output path for the Dicom export
6758dcmExpRealCT.cst = cst ; % set the structure set for the Dicom export
6859dcmExpRealCT.ct = realCTct ; % set the image volume for the Dicom export
7162
7263% In order to create a fake or synthetic CT volume, we will use real CT data and re-use its HU values
7364% for illustrative purposes. First, we extract the original 3D HU image data.
74- fakeCTct= ct ;
65+ fakeCTct = ct ;
7566fakeCTcubeHU = fakeCTct.cubeHU{1 };
7667
7768% Then, we define the range values to be coerced to create the fake CT volume. For this, we will
8576fakeCTct.cubeHU{1 } = fakeCTcubeHU ;
8677
8778% review fake CT volume
88- ct= fakeCTct ;
89- matRadGUI ;
79+ ct = fakeCTct ;
80+ hGUI = matRadGUI ;
9081
9182% saving as DICOMs
9283dcmExpFakeCT= matRad_DicomExporter ; % create instance of matRad_DicomExporter
9687dcmExpFakeCT .matRad_exportDicom();
9788
9889% clear all except of paths, close windows to start from clean space
99- delete(matRadGUI );
100- clearvars - except ' patDir' ' patDirFakeCT' ' patDirRealCT' ' matRad_cfg' ' name' ;
90+ delete(hGUI );
10191
10292
10393%% Patient Data Import from DICOM
113103
114104impRealCT = matRad_DicomImporter(patDirRealCT );
115105matRad_importDicom(impRealCT );
106+
107+ % StructureSet Import does not work in octave
108+ if matRad_cfg .isOctave
109+ cst = realCTcst ;
110+ end
111+
116112% Review the exported file and automatically identified
117113% optimization objectives and constraints.
118114matRadGUI ;
133129 % This is the superclass for all objectives and constraints to enable easy
134130 % identification.
135131
136- if cst{i ,3 } == " OAR"
132+ if strcmp( cst{i ,3 }, ' OAR' )
137133 if ~isempty(regexpi(cst{i ,2 },' spinal' , ' once' ))
138134 objective = DoseObjectives .matRad_MaxDVH ;
139135 objective.penalty = 1 ;
205201
206202pln.propStf.gantryAngles = [0 50 100 150 200 250 300 ];
207203pln.propStf.couchAngles = zeros(1 ,numel(pln .propStf .gantryAngles ));
208- pln.propStf.bixelWidth = 4 ;
204+ pln.propStf.bixelWidth = 5 ;
209205
210206% Obtain the number of beams and voxels from the existing variables and
211207% calculate the iso-center which is per default the center of gravity of
215211
216212%% Dose calculation settings
217213% set resolution of dose calculation and optimization
218- pln.propDoseCalc.doseGrid.resolution.x = 4 ; % [mm]
219- pln.propDoseCalc.doseGrid.resolution.y = 4 ; % [mm]
220- pln.propDoseCalc.doseGrid.resolution.z = 4 ; % [mm]
214+ pln.propDoseCalc.doseGrid.resolution.x = 7 ; % [mm]
215+ pln.propDoseCalc.doseGrid.resolution.y = 7 ; % [mm]
216+ pln.propDoseCalc.doseGrid.resolution.z = 7 ; % [mm]
221217
222218%%
223219% Enable sequencing and disable direct aperture optimization (DAO) for now.
257253% Save DVH parameters calculated on the real CT to the table
258254% for further comparison with plans calculated on the synthetic CT.
259255% Include the patient number and indicate the CT type as "real"
260- % to facilitate delta calculations between the plans later.
261- dvhTableReal= struct2table(qi );
262- % Select only DVH parameters from QI table you are interested in comparison
263- dvhTableReal= dvhTableReal(: ,1 : 9 );
264- dvhTableReal.patient= repmat(char(name ),length(qi ),1 );
265- dvhTableReal.ct_type = repmat(char(" real" ),length(qi ),1 );
266- % Check DVH table for real CT
267- disp(dvhTableReal );
256+ % to facilitate delta calculations between the plans later.
257+ if matRad_cfg .isMatlab % tables not supported by Octave
258+ dvhTableReal= struct2table(qi );
259+ % Select only DVH parameters from QI table you are interested in comparison
260+ dvhTableReal= dvhTableReal(: ,horzcat({' name' },indicators ));
261+ dvhTableReal.patient= repmat(char(name ),length(qi ),1 );
262+ dvhTableReal.ct_type = repmat(' real' ,length(qi ),1 );
263+ % Check DVH table for real CT
264+ disp(dvhTableReal );
265+ end
268266
269267%% Now clear the data from the real CT image, except for the plan parameters,
270268% and load the synthetic (fake) CT image of the same patient.
292290qi = resultGUI .qi ;
293291
294292% In the similar fashion, save DVH parameters calculated on the synthetic CT to the table
295- dvhTableFake= struct2table(qi );
296- % Select the same DVH parameters for further comparison
297- dvhTableFake= dvhTableFake(: ,1 : 9 );
298- dvhTableFake.patient= repmat(char(name ),length(qi ),1 );
299- dvhTableFake.ct_type = repmat(char(" fake" ),length(qi ),1 );
293+ if matRad_cfg .isMatlab % tables not supported by Octave
294+ dvhTableFake= struct2table(qi );
295+ % Select the same DVH parameters for further comparison
296+ dvhTableFake= dvhTableFake(: ,horzcat({' name' },indicators ));
297+ dvhTableFake.patient= repmat(char(name ),length(qi ),1 );
298+ dvhTableFake.ct_type = repmat(' fake' ,length(qi ),1 );
299+ end
300300
301301%% Calculate the difference in between of DVH calculated on real CT (dvh_table_real) and synthetic CT (dvh_table_fake)
302- dvhTableDiff= dvhTableReal ;
303- for i = 1 : height(dvhTableReal )
304- for j= [" mean" ," std" ," max" , " min" , " D_2" ," D_5" , " D_95" , " D_98" ]
305- if cell2mat(dvhTableReal{i ," name" })==cell2mat(dvhTableFake{i ," name" })
306- dvhTableDiff{i ,j }=(dvhTableReal{i ,j }-dvhTableFake{i ,j })/dvhTableReal{i ,j }*100 ;
307- else
308- dvhTableDiff{i ,j }=" error" ;
302+ if matRad_cfg .isMatlab
303+ dvhTableDiff = dvhTableReal ;
304+ for i = 1 : height(dvhTableReal )
305+ for j = indicators
306+ if cell2mat(dvhTableReal{i ,' name' }) == cell2mat(dvhTableFake{i ,' name' })
307+ dvhTableDiff{i ,j } = (dvhTableReal{i ,j } - dvhTableFake{i ,j }) / dvhTableReal{i ,j }*100 ;
308+ else
309+ dvhTableDiff{i ,j } = ' error' ;
310+ end
309311 end
310312 end
311- end
312313
313- %% Check the difference
314- disp(dvhTableDiff )
314+ disp(dvhTableDiff )
315315
316- % %% Save the results to CSV file
317- % file_path_dvh_diff = "YOUR PATH"
318- % writetable(dvh_table_diff, file_path_dvh_diff);
319- delete(matRadGUI );
320- close all ;
321- clear ;
316+ % %% Save the results to CSV file
317+ % file_path_dvh_diff = "YOUR PATH"
318+ % writetable(dvh_table_diff, file_path_dvh_diff);
319+ end
0 commit comments