Skip to content

Commit aecbc90

Browse files
mlapaevawahln
mlapaeva
and
wahln
authored
Synthetic CT example (#827)
* [synthetic CT example]: create dicoms for ct and sCT from liver phantom * [synthetic CT example]: fix bug with non-closing modal window, when dicom is loaded withour structure file * [synthetic CT example]: creating synthetic CT example for DVH diff calcs * Revert "[synthetic CT example]: create dicoms for ct and sCT from liver phantom" This reverts commit 9a92505. * [synthetic CT example]: DICOM import included, and changes from review are applied * [Synthetic CT example]: name is added to author list and citation * [synthetic CT example]: example is added to autotests --------- Co-authored-by: mlapaeva <[email protected]> Co-authored-by: Niklas Wahl <[email protected]>
1 parent ca05aed commit aecbc90

File tree

5 files changed

+329
-2
lines changed

5 files changed

+329
-2
lines changed

‎AUTHORS.txt‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
List of all matRad developers that contributed code (alphabetical)
1+
List of all matRad developers that contributed code (alphabetical)
22

33
* Nelly Abbani
44
* Nabe Al-Hasnawi
@@ -26,6 +26,7 @@
2626
* Navid Khaledi
2727
* Thomas Klinge
2828
* Jeremias Kunz
29+
* Mariia Lapaeva
2930
* Paul Anton Meder
3031
* Henning Mescher
3132
* Lucas-Raphael Müller

‎CITATION.cff‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,8 @@ authors:
6161
given-names: "Thomas"
6262
- family-names: "Kunz"
6363
given-names: "Jeremias"
64+
- family-names: "Lapaeva"
65+
given-names: "Mariia"
6466
- family-names: "Mairani"
6567
given-names: "Andrea"
6668
- family-names: "Meder"

‎examples/matRad_example18_CT_sCT_DVH_difference_photons.m‎

Lines changed: 321 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,321 @@
1+
%% Example: Photon Treatment Plan
2+
%
3+
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4+
%
5+
% Copyright 2017 the matRad development team.
6+
%
7+
% This file is part of the matRad project. It is subject to the license
8+
% terms in the LICENSE file found in the top-level directory of this
9+
% distribution and at https://github.com/e0404/matRad/LICENSE.md. No part
10+
% of the matRad project, including this file, may be copied, modified,
11+
% propagated, or distributed except according to the terms contained in the
12+
% LICENSE file.
13+
%
14+
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15+
16+
17+
%% In this example, we will show:
18+
% (0) how to export .mat data into DICOM format to illustrate a synthetic CT data generation,
19+
% as it is usually stored in DICOM format;
20+
% (i) how to load DICOM patient data into matRad and modify dosimetric
21+
% objectives and constraints;
22+
% (ii) how to set up a photon dose calculation and optimization based on the real CT;
23+
% (iii) how to translate the optimized plan to be applied to a synthetic CT image; and
24+
% (iv) how to calculate the DVH differences between the plans.
25+
26+
%% Preliminary step: DICOM data preparation and export from .mat format to .dcm for further use
27+
% Start with a clean MATLAB environment. We will export Liver phantom data
28+
% from .mat format to .dcm, consisting of 3D body volume slices as well as a structure file and use it as a real CT.
29+
% Additionally, we will slightly modify the intensities of the phantom to generate
30+
% synthetic CT images, illustrating dose difference.
31+
clear;
32+
matRad_cfg = matRad_rc; %If this throws an error, run it from the parent directory first to set the paths
33+
34+
% 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
48+
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);
56+
end
57+
58+
% Now, as the directories are created, let us create the DICOM data. Here, the DICOM
59+
% files will be created from the .mat format.
60+
load('Liver.mat');
61+
realCTct=ct;
62+
%review real CT volume
63+
matRadGUI;
64+
%saving as DICOMs
65+
dcmExpRealCT= matRad_DicomExporter; % create instance of matRad_DicomExporter
66+
dcmExpRealCT.dicomDir = patDirRealCT; % set the output path for the Dicom export
67+
dcmExpRealCT.cst = cst; % set the structure set for the Dicom export
68+
dcmExpRealCT.ct = realCTct; % set the image volume for the Dicom export
69+
dcmExpRealCT.matRad_exportDicom();
70+
71+
72+
% In order to create a fake or synthetic CT volume, we will use real CT data and re-use its HU values
73+
% for illustrative purposes. First, we extract the original 3D HU image data.
74+
fakeCTct=ct;
75+
fakeCTcubeHU = fakeCTct.cubeHU{1};
76+
77+
% Then, we define the range values to be coerced to create the fake CT volume. For this, we will
78+
% coerce the soft tissue range and scale values in the range [0,100] to
79+
% [100,300] and export them as .dcm files.
80+
oldMin = 0; oldMax = 100; % Original range
81+
newMin = 100; newMax = 300; % Target range
82+
%scaling
83+
mask = (fakeCTcubeHU >= oldMin) & (fakeCTcubeHU <= oldMax);
84+
fakeCTcubeHU(mask) = newMin + (fakeCTcubeHU(mask) - oldMin) * (newMax - newMin) / (oldMax - oldMin);
85+
fakeCTct.cubeHU{1} = fakeCTcubeHU;
86+
87+
%review fake CT volume
88+
ct= fakeCTct;
89+
matRadGUI;
90+
91+
%saving as DICOMs
92+
dcmExpFakeCT= matRad_DicomExporter; % create instance of matRad_DicomExporter
93+
dcmExpFakeCT.dicomDir = patDirFakeCT; % set the output path for the Dicom export
94+
dcmExpFakeCT.cst = []; % set the structure set for the Dicom export [we will keep it empty and use further the real CT cst]
95+
dcmExpFakeCT.ct = fakeCTct; % set the image volume for the Dicom export
96+
dcmExpFakeCT.matRad_exportDicom();
97+
98+
%clear all except of paths, close windows to start from clean space
99+
delete(matRadGUI);
100+
clearvars -except 'patDir' 'patDirFakeCT' 'patDirRealCT' 'matRad_cfg' 'name';
101+
102+
103+
%% Patient Data Import from DICOM
104+
% Here we are. Let us import prepared DICOM real CT data. If you
105+
% have your data you can try to replace it with your data in the folder,
106+
% mentioned above.
107+
% Functions for importing DICOM files will search for .dcm files in the directory
108+
% and automatically recognize 3D volumes and structure files.
109+
% The DICOM files will then be read into the workspace as 'ct' and 'cst',
110+
% representing the CT images and the structure set, respectively.
111+
% Ensure that the matRad root directory, along with all its subdirectories,
112+
% is added to the MATLAB search path.
113+
114+
impRealCT = matRad_DicomImporter(patDirRealCT);
115+
matRad_importDicom(impRealCT);
116+
% Review the exported file and automatically identified
117+
% optimization objectives and constraints.
118+
matRadGUI;
119+
120+
%% Modifying Plan Optimization Objectives
121+
% The constraints for all defined volumes of interest (VOIs) are stored in the cst cell.
122+
% When importing DICOM structure files, matRad uses matRad_createCst() to generate this cell.
123+
% Default regular expressions are used to define target objectives. However,
124+
% additional constraints for organs at risk (OARs) need to be added for the plan.
125+
% One of the methods to do this is by directly modifying the cst cell.
126+
% If different volumes with varying optimization objectives need to be loaded,
127+
% refer to the documentation for matRad_createCst() and choose suitable options
128+
% to create a custom implementation. In this case, the plan will be modified directly here.
129+
130+
for i = 1:size(cst, 1)
131+
% Accessing each optimization objective and constain
132+
% Read more about how to set it in matRad_DoseOptimizationFunction - Superclass for objectives and constraints.
133+
% This is the superclass for all objectives and constraints to enable easy
134+
% identification.
135+
136+
if cst{i,3} == "OAR"
137+
if ~isempty(regexpi(cst{i,2},'spinal', 'once'))
138+
objective = DoseObjectives.matRad_MaxDVH;
139+
objective.penalty = 1;
140+
objective.parameters = {12,1}; %dose, to volume
141+
cst{i,6} = {};
142+
cst{i,6}{1} = struct(objective);
143+
144+
elseif ~isempty(regexpi(cst{i,2},'stomach','once')) || ...
145+
~isempty(regexpi(cst{i,2},'duodenum', 'once'))
146+
objective = DoseObjectives.matRad_MaxDVH;
147+
objective.penalty = 1;
148+
objective.parameters = {30,1}; %dose, to volume
149+
cst{i,6} = {};
150+
cst{i,6}{1} = struct(objective);
151+
152+
end
153+
end
154+
end
155+
156+
% Verify that the new objectives have been added and are visible in the
157+
% user interface. We save it for further synthetic CT calculations
158+
matRadGUI;
159+
realCTcst=cst;
160+
161+
%% Treatment Plan
162+
% The next step is to define your treatment plan labeled as 'pln'. This
163+
% matlab structure requires input from the treatment planner and defines
164+
% the most important cornerstones of your treatment plan.
165+
166+
%%
167+
% First of all, we need to define what kind of radiation modality we would
168+
% like to use. Possible values are photons, protons or carbon. In this case
169+
% we want to use photons. Then, we need to define a treatment machine to
170+
% correctly load the corresponding base data. matRad includes base data for
171+
% generic photon linear accelerator called 'Generic'. By this means matRad
172+
% will look for 'photons_Generic.mat' in our root directory and will use
173+
% the data provided in there for dose calculation.
174+
% The number of fractions is set to 30. Internally, matRad considers the
175+
% fraction dose for optimization, however, objetives and constraints are
176+
% defined for the entire treatment.
177+
178+
pln.radiationMode = 'photons';
179+
pln.machine = 'Generic';
180+
pln.numOfFractions = 1;
181+
182+
%%
183+
% Define the biological model used for modeling biological dose (esp. for
184+
% particles).
185+
% Possible biological models are:
186+
% none: use no specific biological model
187+
% constRBE: use a constant RBE
188+
% MCN: use the variable RBE McNamara model for protons
189+
% WED: use the variable RBE Wedenberg model for protons
190+
% LEM: use the biophysical variable RBE Local Effect model for carbons
191+
% As we are using photons, we simply set the parameter to 'none'
192+
pln.bioModel = 'none';
193+
194+
%%
195+
% It is possible to request multiple error scenarios for robustness
196+
% analysis and optimization. Here, we just use the "nominal scenario"
197+
% (nomScen)
198+
pln.multScen = 'nomScen';
199+
200+
%%
201+
% Now we have to set some beam parameters. We can define multiple beam
202+
% angles for the treatment which might be used in generic fashion for quick dose calculation. All corresponding couch angles are set to 0 at this
203+
% point. Moreover, we set the bixelWidth to 4, which results in a beamlet
204+
% size of 4 x 4 mm in the isocenter plane.
205+
206+
pln.propStf.gantryAngles = [0 50 100 150 200 250 300];
207+
pln.propStf.couchAngles = zeros(1,numel(pln.propStf.gantryAngles));
208+
pln.propStf.bixelWidth = 4;
209+
210+
% Obtain the number of beams and voxels from the existing variables and
211+
% calculate the iso-center which is per default the center of gravity of
212+
% all target voxels.
213+
pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles);
214+
pln.propStf.isoCenter = matRad_getIsoCenter(cst,ct,0);
215+
216+
%% Dose calculation settings
217+
% 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]
221+
222+
%%
223+
% Enable sequencing and disable direct aperture optimization (DAO) for now.
224+
% A DAO optimization is shown in a seperate example.
225+
pln.propSeq.runSequencing = 1;
226+
pln.propOpt.runDAO = 0;
227+
228+
%% Generate Beam Geometry STF
229+
% The steering file struct comprises the complete beam geometry along with
230+
% ray position, pencil beam positions and energies, source to axis distance (SAD) etc.
231+
stf = matRad_generateStf(ct,cst,pln);
232+
233+
%% Dose Calculation for real CT
234+
% Let's generate dosimetric information by pre-computing dose influence
235+
% matrices for unit beamlet intensities for real CT of a patient. Having dose influences available
236+
% allows subsequent inverse optimization.
237+
dij = matRad_calcDoseInfluence(ct,cst,stf,pln);
238+
239+
%% Inverse Optimization for IMRT for real CT
240+
% The goal of the fluence optimization is to find a set of beamlet/pencil
241+
% beam weights which yield the best possible dose distribution according to
242+
% the clinical objectives and constraints underlying the radiation
243+
% treatment. Once the optimization has finished, trigger once the GUI to
244+
% visualize the optimized dose cubes.
245+
resultGUI = matRad_fluenceOptimization(dij,cst,pln);
246+
247+
% Get the weights of the optimized plan to apply them to the synthetic CT image of the patient later.
248+
% Call the matRad_planAnalysis function with the prepared arguments to extract DVH
249+
% parameters calculated for the real CT image.
250+
weights=resultGUI.w;
251+
resultGUI = matRad_planAnalysis(resultGUI,ct,cst,stf,pln);
252+
253+
%Get the plan parameters
254+
dvh = resultGUI.dvh;
255+
qi = resultGUI.qi;
256+
257+
% Save DVH parameters calculated on the real CT to the table
258+
% for further comparison with plans calculated on the synthetic CT.
259+
% 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);
268+
269+
%% Now clear the data from the real CT image, except for the plan parameters,
270+
% and load the synthetic (fake) CT image of the same patient.
271+
% It is important that your image sets are compatible (i.e., same number of CT slices,
272+
% same isocenter position, etc.). We will re-use the structure file from real CT with adjusted objectives
273+
% for dose calculations.
274+
delete(matRadGUI);
275+
clear resultGUI ct cst idx qi* dvh dij;
276+
277+
278+
impFakeCT = matRad_DicomImporter(patDirFakeCT);
279+
matRad_importDicom(impFakeCT);
280+
cst=realCTcst;
281+
% Review the exported file and previously identified
282+
% optimization objectives and constraints.
283+
matRadGUI;
284+
285+
%% Perform the dose calculation for synthetic CT by using the weights of plan, calculated on real CT
286+
% (i.e., obtain the dij variable) for the synthetic CT image.
287+
resultGUI = matRad_calcDoseDirect(ct,stf,pln,cst, weights);
288+
resultGUI = matRad_planAnalysis(resultGUI,ct,cst,stf,pln);
289+
matRadGUI;
290+
%Get the plan parameters calculated on synthetic CT
291+
dvh = resultGUI.dvh;
292+
qi = resultGUI.qi;
293+
294+
% 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);
300+
301+
%% 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";
309+
end
310+
end
311+
end
312+
313+
%% Check the difference
314+
disp(dvhTableDiff)
315+
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;

‎matRad/dicom/@matRad_DicomImporter/matRad_importDicom.m‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -189,7 +189,9 @@ function matRad_importDicom(obj)
189189
obj.resultGUI.w = [obj.resultGUI.w; [obj.stf(i).ray.weight]'];
190190
end
191191
end
192-
192+
if any(ishandle(h))
193+
close(h)
194+
end
193195
%% put ct, cst, pln, stf, resultGUI to the workspace
194196
ct = obj.ct;
195197
cst = obj.cst;

‎test/autoExampleTest/test_examples.m‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@
2323
'examples/matRad_example12_simpleParticleMonteCarlo.m',...
2424
'examples/matRad_example15_brachy.m',...
2525
'examples/matRad_example17_biologicalModels.m',...
26+
'examples/matRad_example18_CT_sCT_DVH_difference_photons.m',...
2627
'matRad.m',...
2728
};
2829

0 commit comments

Comments
 (0)