diff --git a/CHANGELOG.md b/CHANGELOG.md index 9850a1b2c..48b76863f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed - possible negative doses in finesampling engine due to extrapolation in kernel interpolation - correct parsing of all optional arguments of the `traceCube` function for `matRad_RayTracer` +- corrected an inconsistency where analytical dose calculations used an RSP cube wit `ignoreOutsideDensities` applied, while MC dose engines often converted materials directly from the HU cube without the same masking. Added the variables `ignoreOutsideDensities` and `useGivenEqDensityCube` to handle this also consistently between STF generation and dose engines. ### Changed - New version of photons_Generic.mat basedata file can now be provided, allowing a "version" field alongside "meta" and "data" files within the machine struct. Version 2 requires correct kernel normalization (without implying a spacing in the convolution integral). photons_Generic.mat has been updated to version 2 with correct kernel normalization. diff --git a/docs/datastructures/pln.rst b/docs/datastructures/pln.rst index cf05b041a..7853e5687 100644 --- a/docs/datastructures/pln.rst +++ b/docs/datastructures/pln.rst @@ -89,9 +89,12 @@ pln.propOpt Setting this value to ``true`` will enable sequencing algorithms run. +.. _pln_additional_properties: + Additional adjustable properties -------------------------------- + The following properties of the pln struct can additionally be adjusted. If they are not explicitly set, default values are used. The default values are handled by the `MatRad_Config class `_. **pln.propStf.longitudinalSpotSpacing** @@ -102,6 +105,18 @@ The following properties of the pln struct can additionally be adjusted. If they If this property is set to *true*, the target is expanded for beamlet finding. Default: *true*. +**pln.propStf.ignoreOutsideDensities** + If this property is set to *true*, the stf generation will ignore densities outside of the patient, this is important for the selection of energies for the particles. Default: *false*. + +**pln.propStf.useGivenEqDensityCube** + If this property is set to *true*, the stf generation calculation will use the given equivalent density cube. Default: *false*. + +**pln.propDoseCalc.ignoreOutsideDensities** + If this property is set to *true*, the dose calculation will ignore densities outside of the patient. This should be the same as **pln.propStf.ignoreOutsideDensities**. Default: *false*. + +**pln.propDoseCalc.useGivenEqDensityCube** + If this property is set to *true*, the dose calculation will use the given equivalent density cube. This should be the same as **pln.propStf.useGivenEqDensityCube**. Default: *false*. + **pln.propDoseCalc.doseGrid.resolution** Specifies the resolution for the dose calculation. Default: x direction: *3* mm, y direction: *3* mm, z direction: *3* mm. diff --git a/docs/datastructures/stf.rst b/docs/datastructures/stf.rst index 4181f1319..f87c4922d 100644 --- a/docs/datastructures/stf.rst +++ b/docs/datastructures/stf.rst @@ -34,7 +34,8 @@ Screenshot of the stf structure: - *Photons*: each ray is a single bixel. *Particles*: :ref:`number of bixels / dose spots along each ray `. * - **total number of bixels** - Total number of bixels for each beam. - + * - **props** + - Additional properties used during generation of the stf struct, used to check consistency to the dose calculation. .. _ray: stf.ray substructure @@ -92,3 +93,10 @@ Screenshot of the stf.numOfBixelsPerRay field: .. image:: /images/stfStructNumOfBixelsScreenshot.png :alt: stf.numOfBixelsPerRay field screenshot + + +stf.props field +=========================== + +The *stf.props* field contains additional properties used during generation of the stf struct, used to check consistency to the dose calculation. +The parameters ignoreOutsideDensities and useGivenEqDensityCube are saved (:ref:`pln_additional_properties`) \ No newline at end of file diff --git a/matRad/MatRad_Config.m b/matRad/MatRad_Config.m index e4b2e66db..0b2e2124e 100644 --- a/matRad/MatRad_Config.m +++ b/matRad/MatRad_Config.m @@ -19,65 +19,65 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - properties - %Logging - logLevel = 4; %1 = only Errors, 2 = with Warnings, 3 = Info output, 4 = deprecation warnings, 5 = debug information - keepLog = false; %Stores the full log in memory - writeLog = false; %Writes the log to a file on-the-fly + % Logging + logLevel = 4 % 1 = only Errors, 2 = with Warnings, 3 = Info output, 4 = deprecation warnings, 5 = debug information + keepLog = false % Stores the full log in memory + writeLog = false % Writes the log to a file on-the-fly - defaults; + defaults - %Disable GUI - disableGUI = false; + % Disable GUI + disableGUI = false - devMode = false; - eduMode = false; + devMode = false + eduMode = false - gui; + gui - %User folders - userfolders; %Cell array of user folders containing machines, patients, hluts. Default contains the userdata folder in the matRad root directory + % User folders + userfolders % Cell array of user folders containing machines, patients, hluts. Default contains the userdata folder in the matRad root directory end - %Deprecated properties referencing a newer one - properties (Dependent,SetAccess = private) - %Default Properties - propDoseCalc; - propOpt; - propStf; + % Deprecated properties referencing a newer one + properties (Dependent, SetAccess = private) + % Default Properties + propDoseCalc + propOpt + propStf end properties (SetAccess = private) - messageLog = {}; - logFileHandle; + messageLog = {} + logFileHandle - %For storing the Environment & its version - env; - envVersion; - isOctave; %Helper bool to check for Octave - isMatlab; %Helper bool to check for Matlab - matRad_version; %MatRad version string + % For storing the Environment & its version + env + envVersion + isOctave % Helper bool to check for Octave + isMatlab % Helper bool to check for Matlab + matRad_version % MatRad version string - matRadRoot; %Path to matRadRoot + matRadRoot % Path to matRadRoot end properties (SetAccess = private, Dependent) - matRadSrcRoot; %Path to matRadSrcRoot ("matRad" subfolder of matRadRoot) - primaryUserFolder; %Points to the first entry in userfolders - exampleFolder; %Contains examples - thirdPartyFolder; %Contains third party tools/libraries used in matRad + matRadSrcRoot % Path to matRadSrcRoot ("matRad" subfolder of matRadRoot) + primaryUserFolder % Points to the first entry in userfolders + exampleFolder % Contains examples + thirdPartyFolder % Contains third party tools/libraries used in matRad end methods (Access = private) + function obj = MatRad_Config() - %MatRad_Config Constructs an instance of this class. + % MatRad_Config Constructs an instance of this class. % The configuration is implemented as a singleton and used globally % Therefore its constructor is private % For instantiation, use the static MatRad_Config.instance(); - %Set + % Set if isdeployed obj.matRadRoot = [ctfroot filesep 'matRad']; else @@ -89,103 +89,105 @@ obj.updateUserfolders(); - %set version + % set version obj.getEnvironment(); obj.matRad_version = matRad_version(obj.matRadRoot); - %Configure Environment + % Configure Environment obj.configureEnvironment(); - %Just to catch people messing with the properties in the file + % Just to catch people messing with the properties in the file if ~isempty(obj.writeLog) && obj.writeLog logFile = [obj.matRadRoot filesep 'matRad.log']; - obj.logFileHandle = fopen(logFile,'a'); + obj.logFileHandle = fopen(logFile, 'a'); end - %Call the reset function for remaining inatialization + % Call the reset function for remaining inatialization obj.reset(); end function delete(~) - %might not be desired by users - %rmpath(genpath(matRad_cfg.matRadRoot)); + % might not be desired by users + % rmpath(genpath(matRad_cfg.matRadRoot)); end - function displayToConsole(obj,type,formatSpec,varargin) - %displayToConsole lowest-level logging function for matRad. + function displayToConsole(obj, type, formatSpec, varargin) + % displayToConsole lowest-level logging function for matRad. % Display to console will be called from the public wrapper % functions dispError, dispWarning, dispInfo, dispDebug % % input: - % type: type of the log information. + % type: type of the log information. % Needs to be one of 'error', 'warning', 'info' or 'debug'. - % formatSpec: string to print using format specifications similar to fprintf - % varargin: variables according to formatSpec + % formatSpec: string to print using format specifications similar to fprintf + % varargin: variables according to formatSpec if nargin < 4 forwardArgs = {formatSpec}; else - forwardArgs = [{formatSpec},varargin(:)']; + forwardArgs = [{formatSpec}, varargin(:)']; end if obj.keepLog - obj.messageLog{end+1,1} = upper(type); - obj.messageLog{end,2} = sprintf(forwardArgs{:}); + obj.messageLog{end + 1, 1} = upper(type); + obj.messageLog{end, 2} = sprintf(forwardArgs{:}); end switch type - case{'info'} + case {'info'} if obj.logLevel >= 3 fprintf(forwardArgs{:}); end - case{'debug'} + case {'debug'} if obj.logLevel >= 5 forwardArgs{1} = ['DEBUG: ' forwardArgs{1}]; fprintf(forwardArgs{:}); end - case{'dep'} + case {'dep'} if obj.logLevel >= 4 forwardArgs{1} = ['DEPRECATION WARNING: ' forwardArgs{1}]; - warning('matRad:Deprecated',forwardArgs{:}); + warning('matRad:Deprecated', forwardArgs{:}); end - case{'warning'} + case {'warning'} if obj.logLevel >= 2 - warning('matRad:Warning',forwardArgs{:}); + warning('matRad:Warning', forwardArgs{:}); end case {'error'} if obj.logLevel >= 1 - %We create an error structure to later clean the - %stack trace from the last two files/lines (i.e., - %this function / file) + % We create an error structure to later clean the + % stack trace from the last two files/lines (i.e., + % this function / file) err.message = sprintf(forwardArgs{:}); err.identifier = 'matRad:Error'; - %err.stack = dbstack(2); + % err.stack = dbstack(2); error(err); end otherwise - error('Log type %s not defined!',type); + error('Log type %s not defined!', type); end if obj.writeLog - fprintf(obj.logFileHandle,forwardArgs{:}); + fprintf(obj.logFileHandle, forwardArgs{:}); end end + end methods + function reset(obj) - %Set all default properties for matRad's computations + % Set all default properties for matRad's computations obj.setDefaultProperties(); obj.setDefaultGUIProperties(); end function setDefaultProperties(obj) - %setDefaultProperties set matRad's default computation + % setDefaultProperties set matRad's default computation % properties - %Default machines + % Default machines obj.defaults.machine.photons = 'Generic'; obj.defaults.machine.protons = 'Generic'; obj.defaults.machine.helium = 'Generic'; @@ -194,8 +196,7 @@ function setDefaultProperties(obj) obj.defaults.machine.fallback = 'Generic'; obj.defaults.machine.VHEE = 'Generic'; - - %Default Bio Model + % Default Bio Model obj.defaults.bioModel.photons = 'none'; obj.defaults.bioModel.protons = 'constRBE'; obj.defaults.bioModel.helium = 'HEL'; @@ -204,23 +205,25 @@ function setDefaultProperties(obj) obj.defaults.bioModel.fallback = 'none'; obj.defaults.bioModel.VHEE = 'none'; - %Default Steering/Geometry Properties - obj.defaults.propStf.generator = {'PhotonIMRT','ParticleIMPT','SimpleBrachy','VHEE'}; + % Default Steering/Geometry Properties + obj.defaults.propStf.generator = {'PhotonIMRT', 'ParticleIMPT', 'SimpleBrachy', 'VHEE'}; obj.defaults.propStf.longitudinalSpotSpacing = 2; - obj.defaults.propStf.addMargin = true; %expand target for beamlet finding + obj.defaults.propStf.addMargin = true; % expand target for beamlet finding obj.defaults.propStf.bixelWidth = 5; - - %Dose Calculation Options - obj.defaults.propDoseCalc.engine = {'SVDPB','HongPB'}; %Names for default engines used when no other is given - obj.defaults.propDoseCalc.doseGrid.resolution = struct('x',3,'y',3,'z',3); %[mm] - obj.defaults.propDoseCalc.dosimetricLateralCutOff = 0.995; %[rel.] - obj.defaults.propDoseCalc.geometricLateralCutOff = 50; %[mm] - obj.defaults.propDoseCalc.kernelCutOff = Inf; %[mm] - obj.defaults.propDoseCalc.ssdDensityThreshold = 0.05; %[rel.] - obj.defaults.propDoseCalc.useGivenEqDensityCube = false; %Use the given density cube ct.cube and omit conversion from cubeHU. - obj.defaults.propDoseCalc.ignoreOutsideDensities = true; %Ignore densities outside of cst contours - obj.defaults.propDoseCalc.useCustomPrimaryPhotonFluence = false; %Use a custom primary photon fluence - obj.defaults.propDoseCalc.calcLET = true; %calculate LETs for particles + obj.defaults.propStf.ignoreOutsideDensities = false; % Ignore densities outside of cst contours + obj.defaults.propStf.useGivenEqDensityCube = false; % Use a custom primary photon fluence + + % Dose Calculation Options + obj.defaults.propDoseCalc.engine = {'SVDPB', 'HongPB'}; % Names for default engines used when no other is given + obj.defaults.propDoseCalc.doseGrid.resolution = struct('x', 3, 'y', 3, 'z', 3); % [mm] + obj.defaults.propDoseCalc.dosimetricLateralCutOff = 0.995; % [rel.] + obj.defaults.propDoseCalc.geometricLateralCutOff = 50; % [mm] + obj.defaults.propDoseCalc.kernelCutOff = Inf; % [mm] + obj.defaults.propDoseCalc.ssdDensityThreshold = 0.05; % [rel.] + obj.defaults.propDoseCalc.useGivenEqDensityCube = false; % Use the given density cube ct.cube and omit conversion from cubeHU. + obj.defaults.propDoseCalc.ignoreOutsideDensities = false; % Ignore densities outside of cst contours + obj.defaults.propDoseCalc.useCustomPrimaryPhotonFluence = false; % Use a custom primary photon fluence + obj.defaults.propDoseCalc.calcLET = true; % calculate LETs for particles obj.defaults.propDoseCalc.selectVoxelsInScenarios = 'all'; obj.defaults.propDoseCalc.airOffsetCorrection = true; @@ -228,23 +231,21 @@ function setDefaultProperties(obj) obj.defaults.propDoseCalc.fineSampling.sigmaSub = 1; obj.defaults.propDoseCalc.fineSampling.N = 2; obj.defaults.propDoseCalc.fineSampling.method = 'fitCircle'; - %Monte Carlo options + % Monte Carlo options obj.defaults.propDoseCalc.numHistoriesPerBeamlet = 2e4; obj.defaults.propDoseCalc.numHistoriesDirect = 1e6; obj.defaults.propDoseCalc.outputMCvariance = true; - %Optimization Options + % Optimization Options obj.defaults.propOpt.optimizer = 'IPOPT'; obj.defaults.propOpt.maxIter = 500; obj.defaults.propOpt.runDAO = 0; obj.defaults.propOpt.clearUnusedVoxels = false; obj.defaults.propOpt.enableGPU = false; - %Sequencing Options + % Sequencing Options obj.defaults.propSeq.sequencer = 'siochi'; - - obj.disableGUI = false; obj.defaults.samplingScenarios = 25; @@ -254,30 +255,30 @@ function setDefaultProperties(obj) end - %%For testing + %% For testing function setDefaultPropertiesForTesting(obj) - %setDefaultPropertiesForTesting sets matRad's default - %properties during testing to reduce computational load + % setDefaultPropertiesForTesting sets matRad's default + % properties during testing to reduce computational load obj.setDefaultProperties(); - obj.logLevel = 1; %Omit output except errors + obj.logLevel = 1; % Omit output except errors - %Default Steering/Geometry Properties + % Default Steering/Geometry Properties obj.defaults.propStf.longitudinalSpotSpacing = 20; obj.defaults.propStf.bixelWidth = 20; - %Dose Calculation Options - obj.defaults.propDoseCalc.doseGrid.resolution = struct('x',5,'y',6,'z',7); %[mm] + % Dose Calculation Options + obj.defaults.propDoseCalc.doseGrid.resolution = struct('x', 5, 'y', 6, 'z', 7); % [mm] obj.defaults.propDoseCalc.geometricLateralCutOff = 20; - obj.defaults.propDoseCalc.dosimetricLateralCutOff = 0.8; - obj.defaults.propDoseCalc.kernelCutOff = 20; %[mm] + obj.defaults.propDoseCalc.dosimetricLateralCutOff = 0.9; + obj.defaults.propDoseCalc.kernelCutOff = 20; % [mm] - %Monte Carlo options + % Monte Carlo options obj.defaults.propDoseCalc.numHistoriesPerBeamlet = 100; obj.defaults.propDoseCalc.numHistoriesDirect = 100; - %Optimization Options + % Optimization Options obj.defaults.propOpt.maxIter = 10; obj.defaults.samplingScenarios = 2; @@ -288,20 +289,20 @@ function setDefaultPropertiesForTesting(obj) obj.eduMode = false; end - %%for edu mode + %% for edu mode function setDefaultPropertiesForEduMode(obj) obj.setDefaultProperties(); obj.logLevel = 2; - %Default Steering/Geometry Properties + % Default Steering/Geometry Properties obj.defaults.propStf.longitudinalSpotSpacing = 3; - %Dose calculation options - obj.defaults.propDoseCalc.resolution = struct('x',4,'y',4,'z',4); %[mm] - obj.defaults.propDoseCalc.lateralCutOff = 0.975; %[rel.] + % Dose calculation options + obj.defaults.propDoseCalc.resolution = struct('x', 4, 'y', 4, 'z', 4); % [mm] + obj.defaults.propDoseCalc.lateralCutOff = 0.975; % [rel.] - %Optimization Options + % Optimization Options obj.defaults.propOpt.maxIter = 500; obj.disableGUI = false; @@ -311,19 +312,19 @@ function setDefaultPropertiesForEduMode(obj) end function setDefaultGUIProperties(obj) - %Detect current theme + % Detect current theme light = false; try if ispc - light = logical(winqueryreg('HKEY_CURRENT_USER','Software\\Microsoft\\Windows\\CurrentVersion\\Themes\\Personalize','AppsUseLightTheme')); + light = logical(winqueryreg('HKEY_CURRENT_USER', 'Software\\Microsoft\\Windows\\CurrentVersion\\Themes\\Personalize', 'AppsUseLightTheme')); elseif ismac out = system('defaults read -g AppleInterfaceStyle'); - if ~strcmp(out,'Dark') + if ~strcmp(out, 'Dark') light = true; end else out = system('gsettings get org.gnome.desktop.interface color-scheme'); - if strcmp(out,'prefer-light') + if strcmp(out, 'prefer-light') light = true; end end @@ -340,35 +341,35 @@ function setDefaultGUIProperties(obj) obj.gui = struct(theme); end - function dispDebug(obj,formatSpec,varargin) - %dispDebug print debug messages (log level >= 4) + function dispDebug(obj, formatSpec, varargin) + % dispDebug print debug messages (log level >= 4) % % input: - % formatSpec: string to print using format specifications similar to fprintf - % varargin: variables according to formatSpec + % formatSpec: string to print using format specifications similar to fprintf + % varargin: variables according to formatSpec - obj.displayToConsole('debug',formatSpec,varargin{:}); + obj.displayToConsole('debug', formatSpec, varargin{:}); end - function dispInfo(obj,formatSpec,varargin) - %dispInfo print information console output (log level >= 3) + function dispInfo(obj, formatSpec, varargin) + % dispInfo print information console output (log level >= 3) % % input: - % formatSpec: string to print using format specifications similar to fprintf - % varargin: variables according to formatSpec - obj.displayToConsole('info',formatSpec,varargin{:}); + % formatSpec: string to print using format specifications similar to fprintf + % varargin: variables according to formatSpec + obj.displayToConsole('info', formatSpec, varargin{:}); end - function dispError(obj,formatSpec,varargin) - %dispError print errors (forwarded to "error" that will stop the program) (log level >= 1) + function dispError(obj, formatSpec, varargin) + % dispError print errors (forwarded to "error" that will stop the program) (log level >= 1) % % input: - % formatSpec: string to print using format specifications + % formatSpec: string to print using format specifications % similar to 'error' - % varargin: variables according to formatSpec + % varargin: variables according to formatSpec try - obj.displayToConsole('error',formatSpec,varargin{:}); + obj.displayToConsole('error', formatSpec, varargin{:}); catch ME if obj.isOctave ME.stack = ME.stack(3:end); % Removes the dispError and dispToConsole from the stack @@ -379,94 +380,94 @@ function dispError(obj,formatSpec,varargin) end end - function dispWarning(obj,formatSpec,varargin) - %dispError print warning (forwarded to 'warning') (log level >= 2) + function dispWarning(obj, formatSpec, varargin) + % dispError print warning (forwarded to 'warning') (log level >= 2) % % input: - % formatSpec: string to print using format specifications + % formatSpec: string to print using format specifications % similar to 'warning' - % varargin: variables according to formatSpec - obj.displayToConsole('warning',formatSpec,varargin{:}); + % varargin: variables according to formatSpec + obj.displayToConsole('warning', formatSpec, varargin{:}); end - function dispDeprecationWarning(obj,formatSpec,varargin) - %dispDeprecationWarning wrapper for deprecation warnings forwarded to displayToConsole - obj.displayToConsole('dep',formatSpec,varargin{:}); + function dispDeprecationWarning(obj, formatSpec, varargin) + % dispDeprecationWarning wrapper for deprecation warnings forwarded to displayToConsole + obj.displayToConsole('dep', formatSpec, varargin{:}); end - function obj = writeLogToFile(obj,filename) - %writeLogToFile writes the log kept in MatRad_Config to file. + function obj = writeLogToFile(obj, filename) + % writeLogToFile writes the log kept in MatRad_Config to file. % Note that the switch keepLog must be enabled for MatRad_Config to store all logging output. singleString = '%s: %s\n'; - fID = fopen(filename,'w'); - fprintf(fID,repmat(singleString,1,size(obj.messageLog,1)),obj.messageLog{:}); + fID = fopen(filename, 'w'); + fprintf(fID, repmat(singleString, 1, size(obj.messageLog, 1)), obj.messageLog{:}); fclose(fID); end - function set.logLevel(obj,newLogLevel) - %%Property set methods for logLevel + function set.logLevel(obj, newLogLevel) + %% Property set methods for logLevel minLevel = 1; maxLevel = 5; if newLogLevel >= minLevel && newLogLevel <= maxLevel obj.logLevel = newLogLevel; else - obj.dispError('Invalid log level. Value must be between %d and %d',minLevel,maxLevel); + obj.dispError('Invalid log level. Value must be between %d and %d', minLevel, maxLevel); end end - function set.userfolders(obj,userfolders) + function set.userfolders(obj, userfolders) oldFolders = obj.userfolders; - %Check if folders need to be created + % Check if folders need to be created for f = 1:numel(userfolders) if ~isfolder(userfolders{f}) [status, msg] = mkdir(userfolders{f}); if status == 0 - obj.dispWarning('Userfolder %s not added because it could not be created: %s',userfolders{f},msg); + obj.dispWarning('Userfolder %s not added because it could not be created: %s', userfolders{f}, msg); else - subfolders = {'hluts','machines','patients','scripts'}; - [status,msgs] = cellfun(@(sub) mkdir([userfolders{f} filesep sub]),subfolders,'UniformOutput',false); + subfolders = {'hluts', 'machines', 'patients', 'scripts'}; + [status, msgs] = cellfun(@(sub) mkdir([userfolders{f} filesep sub]), subfolders, 'UniformOutput', false); if any(cell2mat(status) ~= 1) - obj.dispWarning('Problem when creating subfolder in Userfolder %s!',userfolders{f}) + obj.dispWarning('Problem when creating subfolder in Userfolder %s!', userfolders{f}); end end end end - %We do this to verify folders - nonWorkingFolders = cellfun(@isempty,userfolders); + % We do this to verify folders + nonWorkingFolders = cellfun(@isempty, userfolders); userfolders(nonWorkingFolders) = []; - allNewFolders = cellfun(@dir, userfolders,'UniformOutput',false); + allNewFolders = cellfun(@dir, userfolders, 'UniformOutput', false); if isempty(allNewFolders) obj.dispWarning('No user folders specified. Defaulting to userdata folder in matRad root directory.'); if ~isdeployed - allNewFolders = {[fileparts(mfilename('fullpath')) filesep 'userdata' filesep]}; %We don't access obj.matRadRoot here because of Matlab's weird behavior with properties + allNewFolders = {[fileparts(mfilename('fullpath')) filesep 'userdata' filesep]}; % We don't access obj.matRadRoot here because of Matlab's weird behavior with properties else - allNewFolders = {[ctfroot filesep 'userdata' filesep]}; %We don't access obj.matRadRoot here because of Matlab's weird behavior with properties + allNewFolders = {[ctfroot filesep 'userdata' filesep]}; % We don't access obj.matRadRoot here because of Matlab's weird behavior with properties end end - cleanedNewFolders = cellfun(@(x) x(1).folder,allNewFolders,'UniformOutput',false); + cleanedNewFolders = cellfun(@(x) x(1).folder, allNewFolders, 'UniformOutput', false); % Identify newly added folder paths and add them to path if ~isdeployed - if ~isempty(oldFolders) %if statement for octave compatibility + if ~isempty(oldFolders) % if statement for octave compatibility addedFolders = setdiff(cleanedNewFolders, oldFolders); else addedFolders = cleanedNewFolders; end - addedFolders = cellfun(@genpath,addedFolders,'UniformOutput',false); - addedFolders = strjoin(addedFolders,pathsep); + addedFolders = cellfun(@genpath, addedFolders, 'UniformOutput', false); + addedFolders = strjoin(addedFolders, pathsep); addpath(addedFolders); end % Identify removed folder paths - if ~isempty(oldFolders) %if statement for octave compatibility + if ~isempty(oldFolders) % if statement for octave compatibility removedFolders = setdiff(oldFolders, cleanedNewFolders); - removedFolders = cellfun(@genpath,removedFolders,'UniformOutput',false); - removedFolders = strjoin(removedFolders,pathsep); + removedFolders = cellfun(@genpath, removedFolders, 'UniformOutput', false); + removedFolders = strjoin(removedFolders, pathsep); if ~isdeployed rmpath(removedFolders); end @@ -521,10 +522,10 @@ function dispDeprecationWarning(obj,formatSpec,varargin) end end - function set.writeLog(obj,writeLog) + function set.writeLog(obj, writeLog) if writeLog logFile = [obj.matRadRoot filesep 'matRad.log']; - obj.logFileHandle = fopen(logFile,'a'); + obj.logFileHandle = fopen(logFile, 'a'); obj.writeLog = true; else fclose(obj.logFileHandle); @@ -550,9 +551,9 @@ function getEnvironment(obj) end end - function pln = getDefaultProperties(obj,pln,fields) + function pln = getDefaultProperties(obj, pln, fields) % Function to load all non-set parameters into pln struct - standardFields = {'propDoseCalc','propOpt','propStf'}; + standardFields = {'propDoseCalc', 'propOpt', 'propStf'}; % Check if only one argument was given if ~iscell(fields) @@ -562,12 +563,12 @@ function getEnvironment(obj) for i = 1:length(fields) currField = fields{i}; - if ismember(currField,standardFields) + if ismember(currField, standardFields) % Get defaults for standard fields that can easily be read from set default values - if ~isfield(pln,currField) + if ~isfield(pln, currField) pln.(currField) = obj.defaults.(currField); else - pln.(currField) = matRad_recursiveFieldAssignment(pln.(currField),obj.defaults.(currField),false); + pln.(currField) = matRad_recursiveFieldAssignment(pln.(currField), obj.defaults.(currField), false); end end end @@ -575,13 +576,13 @@ function getEnvironment(obj) function configureEnvironment(obj) if obj.isOctave - struct_levels_to_print(0); %Disables full printing of struct array fields - warning("off","Octave:data-file-in-path"); %Disables warning of loading patients from the data folder + struct_levels_to_print(0); % Disables full printing of struct array fields + warning("off", "Octave:data-file-in-path"); % Disables warning of loading patients from the data folder end end function userfolders = updateUserfolders(obj) - %obtain userfolders + % obtain userfolders if ispc userdir = getenv('USERPROFILE'); else @@ -594,7 +595,7 @@ function configureEnvironment(obj) else userfolders = {[obj.matRadRoot filesep 'userdata' filesep]}; if isfolder(userfolderInHomeDir) - userfolders{end+1} = userfolderInHomeDir; + userfolders{end + 1} = userfolderInHomeDir; end end @@ -606,26 +607,27 @@ function configureEnvironment(obj) for i = 1:numel(envPaths) p = strtrim(envPaths{i}); if ~isempty(p) && isfolder(p) - userfolders{end+1} = p; + userfolders{end + 1} = p; end end end obj.userfolders = userfolders; end + end - %methods (Access = private) + % methods (Access = private) % function renameFields(obj) - %end + % end methods (Static) function obj = instance() - %instance creates a singleton instance of MatRad_Config + % instance creates a singleton instance of MatRad_Config % In MatRad_Config, the constructor is private to make sure only on global instance exists. % Call this static function to get or create an instance of the matRad configuration class - persistent uniqueInstance; + persistent uniqueInstance if isempty(uniqueInstance) obj = MatRad_Config(); @@ -647,7 +649,7 @@ function configureEnvironment(obj) fields = fieldnames(basic_struct); for k = 1:length(fields) disp(fields{k}); - if(isfield(changed_struct, fields{k})) + if isfield(changed_struct, fields{k}) if isstruct(changed_struct.(fields{k})) && isstruct(basic_struct.(fields{i})) basic_struct.(fields{k}) = mergeStructs(basic_struct.(fields{k}), changed_struct.(fields{i})); else @@ -671,14 +673,14 @@ function configureEnvironment(obj) % Throw warning if the version differs and remove the % matRad_version field from the loaded struct, in order to % not overwrite the version later - if (isfield(sobj, 'matRad_version') && ~(strcmp(obj.matRad_version, sobj.matRad_version))) + if isfield(sobj, 'matRad_version') && ~(strcmp(obj.matRad_version, sobj.matRad_version)) warning('MatRad version or git Branch of the loaded object differs from the current version!'); sobj = rmfield(sobj, 'matRad_version'); end % Iterate over the properties of the newly created MatRad_Config object for i = 1:length(props) % check if the field exists in the loaded object - if(isfield(sobj,props{i})) + if isfield(sobj, props{i}) objField = obj.(props{i}); sobjField = sobj.(props{i}); % If field from loaded object and from the newly @@ -686,8 +688,8 @@ function configureEnvironment(obj) % value of the loaded object and if it's a struct % check it's field recursively if ~(isequal(sobjField, objField)) - if (isstruct(sobjField) && isstruct(objField)) - retStruct = mergeStructs(objField,sobjField); + if isstruct(sobjField) && isstruct(objField) + retStruct = mergeStructs(objField, sobjField); obj.(props{i}) = retStruct; else obj.(props{i}) = sobjField; @@ -700,6 +702,5 @@ function configureEnvironment(obj) end end - end end diff --git a/matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/matRad_DoseEngineBase.m b/matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/matRad_DoseEngineBase.m index c1f9cda17..7ad8086f0 100644 --- a/matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/matRad_DoseEngineBase.m +++ b/matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/matRad_DoseEngineBase.m @@ -36,6 +36,8 @@ precision = 'double'; % floating point precision for the dij and computations. enableGPU = false; % whether to use GPU arrays (experimental) for dose calculation (if supported by subclass implementation). %bioModel; % name of the biological model + ignoreOutsideDensities % Ignore densities outside of cst contours + useGivenEqDensityCube; % Use the given density cube ct.cube and omit conversion from cubeHU. end % Protected properties with public get access @@ -215,6 +217,7 @@ function assignBioModelPropertiesFromPln(this, plnModel, warnWhenPropertyChanged %Set direct dose calculation and compute "dij" this.directWeights = w; this.calcDoseDirect = true; + ct = this.preprocessCT(ct,cst,stf); dij = this.calcDose(ct,cst,stf); dij = this.finalizeDose(dij); @@ -275,17 +278,51 @@ function assignBioModelPropertiesFromPln(this, plnModel, warnWhenPropertyChanged function dij = calcDoseInfluence(this,ct,cst,stf) this.calcDoseDirect = false; + ct = this.preprocessCT(ct,cst,stf); dij = this.calcDose(ct,cst,stf); dij = this.finalizeDose(dij); end function setDefaults(this) % future code for property validation on creation here matRad_cfg = MatRad_Config.instance(); - %Assign default parameters from MatRad_Config - this.doseGrid = matRad_cfg.defaults.propDoseCalc.doseGrid; this.multScen = 'nomScen'; + this.doseGrid = matRad_cfg.defaults.propDoseCalc.doseGrid; this.selectVoxelsInScenarios = matRad_cfg.defaults.propDoseCalc.selectVoxelsInScenarios; + this.ignoreOutsideDensities = matRad_cfg.defaults.propDoseCalc.ignoreOutsideDensities; + this.useGivenEqDensityCube = matRad_cfg.defaults.propDoseCalc.useGivenEqDensityCube; + end + + function ct = preprocessCT(this,ct,cst,stf) + matRad_cfg = MatRad_Config.instance(); + % process ct + V = [cst{:,4}]; + V = unique(vertcat(V{:})); + % check consistent with stf + if isfield(stf(1),'props') && this.ignoreOutsideDensities ~= stf(1).props.ignoreOutsideDensities + matRad_cfg.dispWarning('The parameter ignoreOutsideDensities is inconsistant between stf generation and dose calculation') + end + if isfield(stf(1),'props') && this.useGivenEqDensityCube ~= stf(1).props.useGivenEqDensityCube + matRad_cfg.dispWarning('The parameter useGivenEqDensityCube is inconsistant between stf generation and dose calculation') + end + if this.ignoreOutsideDensities + % ignore densities outside of contours + eraseCtDensMask = ones(prod(ct.cubeDim), 1); + eraseCtDensMask(V) = 0; + for i = 1:ct.numOfCtScen + ct.cubeHU{i}(eraseCtDensMask == 1) = -1000; + end + end + if this.useGivenEqDensityCube && ~isfield(ct,'cube') + matRad_cfg.dispWarning('HU Conversion requested to be omitted but no ct.cube exists! Will override and do the conversion anyway!'); + this.useGivenEqDensityCube = false; + end + + if this.useGivenEqDensityCube + matRad_cfg.dispInfo('Omitting HU to rED/rSP conversion and using existing ct.cube!\n'); + else + ct = matRad_calcWaterEqD(ct, stf(1).radiationMode); % Maybe we can avoid duplicating the CT here? + end end end diff --git a/matRad/doseCalc/+DoseEngines/@matRad_ParticleFREDEngine/matRad_ParticleFREDEngine.m b/matRad/doseCalc/+DoseEngines/@matRad_ParticleFREDEngine/matRad_ParticleFREDEngine.m index 1deb506e2..a25e779fe 100644 --- a/matRad/doseCalc/+DoseEngines/@matRad_ParticleFREDEngine/matRad_ParticleFREDEngine.m +++ b/matRad/doseCalc/+DoseEngines/@matRad_ParticleFREDEngine/matRad_ParticleFREDEngine.m @@ -74,7 +74,6 @@ printOutput primaryMass numOfNucleons - ignoreOutsideDensities workingDir forceDijFormatVersion end @@ -536,7 +535,6 @@ function setDefaults(this) this.numOfNucleons = 1; this.outputMCvariance = false; this.constantRBE = NaN; - this.ignoreOutsideDensities = false; end diff --git a/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m b/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m index c835449cd..a29af436f 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m +++ b/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m @@ -25,8 +25,6 @@ dosimetricLateralCutOff; %relative dosimetric cut-off (in fraction of values calculated) ssdDensityThreshold; % Threshold for SSD computation - useGivenEqDensityCube; % Use the given density cube ct.cube and omit conversion from cubeHU. - ignoreOutsideDensities; % Ignore densities outside of cst contours numOfDijFillSteps = 10; % Number of times during dose calculation the temporary containers are moved to a sparse matrix @@ -171,32 +169,12 @@ function setDefaults(this) dij = initDoseCalc@DoseEngines.matRad_DoseEngineBase(this,ct,cst,stf); - % calculate rED or rSP from HU or take provided wedCube - if this.useGivenEqDensityCube && ~isfield(ct,'cube') - matRad_cfg.dispWarning('HU Conversion requested to be omitted but no ct.cube exists! Will override and do the conversion anyway!'); - this.useGivenEqDensityCube = false; - end - - if this.useGivenEqDensityCube - matRad_cfg.dispInfo('Omitting HU to rED/rSP conversion and using existing ct.cube!\n'); - else - ct = matRad_calcWaterEqD(ct, stf); % Maybe we can avoid duplicating the CT here? - end this.cubeWED = cellfun(@(x) cast(x,this.precision),ct.cube, 'UniformOutput',false); if isfield(ct,'hlut') this.hlut = cast(ct.hlut,this.precision); end - % ignore densities outside of contours - if this.ignoreOutsideDensities - eraseCtDensMask = ones(prod(ct.cubeDim),1); - eraseCtDensMask(this.VctGrid) = 0; - for i = 1:ct.numOfCtScen - this.cubeWED{i}(eraseCtDensMask == 1) = 0; - end - end - % Allocate memory for quantity containers dij = this.allocateQuantityMatrixContainers(dij,{'physicalDose'}); diff --git a/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m b/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m index bf339955e..13dddc10b 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m @@ -29,7 +29,6 @@ properties hlut; - useGivenEqDensityCube; % Use the given density cube ct.cube and omit conversion from cubeHU. calcLET = false; calcBioDose = false; prescribedDose = []; @@ -189,9 +188,6 @@ function setDefaults(this) this.setDefaults@DoseEngines.matRad_MonteCarloEngineAbstract(); matRad_cfg = MatRad_Config.instance(); %Instance of matRad configuration class - - this.useGivenEqDensityCube = matRad_cfg.defaults.propDoseCalc.useGivenEqDensityCube; - % Default execution paths are set here this.topasFolder = [matRad_cfg.matRadSrcRoot filesep 'doseCalc' filesep 'topas' filesep]; this.workingDir = [matRad_cfg.primaryUserFolder filesep 'TOPAS' filesep]; @@ -709,18 +705,6 @@ function writeAllFiles(obj,ct,cst,stf,machine,w) dij = this.initDoseCalc@DoseEngines.matRad_MonteCarloEngineAbstract(ct,cst,stf); matRad_cfg = MatRad_Config.instance(); - % calculate rED or rSP from HU or take provided wedCube - if this.useGivenEqDensityCube && ~isfield(ct,'cube') - matRad_cfg.dispWarning('HU Conversion requested to be omitted but no ct.cube exists! Will override and do the conversion anyway!'); - this.useGivenEqDensityCube = false; - end - - if this.useGivenEqDensityCube - matRad_cfg.dispInfo('Omitting HU to rED/rSP conversion and using existing ct.cube!\n'); - else - ct = matRad_calcWaterEqD(ct, stf); % Maybe we can avoid duplicating the CT here? - end - if isfield(ct,'hlut') this.hlut = ct.hlut; else diff --git a/matRad/steering/matRad_StfGeneratorBase.m b/matRad/steering/matRad_StfGeneratorBase.m index e269d7671..39cf5de40 100644 --- a/matRad/steering/matRad_StfGeneratorBase.m +++ b/matRad/steering/matRad_StfGeneratorBase.m @@ -33,6 +33,8 @@ radiationMode; %Radiation Mode machine; %Machine enableGPU = false; %Enable computation on the GPU (experimenta, default false) + ignoreOutsideDensities % Ignore densities outside of cst contours + useGivenEqDensityCube; % Use the given density cube ct.cube and omit conversion from cubeHU. end properties (Access = protected) @@ -192,6 +194,10 @@ function assignPropertiesFromPln(this,pln,warnWhenPropertyChanged) this.initialize(); this.createPatientGeometry(); stf = this.generateSourceGeometry(); + propStruct = this.getProperties(); + for i = 1:size(stf,2) + stf(i).props = propStruct; + end end end @@ -282,13 +288,25 @@ function createPatientGeometry(this) % take only voxels inside patient V = [this.cst{:,4}]; V = unique(vertcat(V{:})); + if this.ignoreOutsideDensities + % ignore densities outside of contours + eraseCtDensMask = ones(prod(this.ct.cubeDim), 1); + eraseCtDensMask(V) = 0; + for i = 1:this.ct.numOfCtScen + this.ct.cubeHU{i}(eraseCtDensMask == 1) = -1000; + end + end + if this.useGivenEqDensityCube && ~isfield(this.ct,'cube') + matRad_cfg.dispWarning('HU Conversion requested to be omitted but no ct.cube exists! Will override and do the conversion anyway!'); + this.useGivenEqDensityCube = false; + end - % ignore densities outside of contours - eraseCtDensMask = ones(prod(this.ct.cubeDim), 1); - eraseCtDensMask(V) = 0; - for i = 1:this.ct.numOfCtScen - this.ct.cube{i}(eraseCtDensMask == 1) = 0; + if this.useGivenEqDensityCube + matRad_cfg.dispInfo('Omitting HU to rED/rSP conversion and using existing ct.cube!\n'); + else + this.ct = matRad_calcWaterEqD(this.ct, this.radiationMode); % Maybe we can avoid duplicating the CT here? end + end function pbMargin = getPbMargin(this) @@ -305,7 +323,12 @@ function createPatientGeometry(this) % make the whole calculation more modular) throw(MException('MATLAB:class:AbstractMember','Abstract function generateSourceGeometry of your StfGenerator needs to be implemented!')); end - end + + function s = getProperties(this) + s.ignoreOutsideDensities = this.ignoreOutsideDensities; + s.useGivenEqDensityCube = this.useGivenEqDensityCube; + end + end methods (Static) function generator = getGeneratorFromPln(pln, warnDefault) @@ -505,5 +528,6 @@ function createPatientGeometry(this) function machine = loadMachine(radiationMode,machineName) machine = matRad_loadMachine(struct('radiationMode',radiationMode,'machine',machineName)); end + end end diff --git a/matRad/util/matRad_interp3.m b/matRad/util/matRad_interp3.m index 57e7da202..f1b5db34a 100644 --- a/matRad/util/matRad_interp3.m +++ b/matRad/util/matRad_interp3.m @@ -1,5 +1,5 @@ function y = matRad_interp3(xi,yi,zi,x,xq,yq,zq,mode,extrapVal) -% interpolates 3-D data (table lookup) +% interpolates 3-D data (table lookup) % % call: % y = matRad_interp3(xi,yi,zi,x,xq,yq,zy) @@ -7,14 +7,14 @@ % y = matRad_interp3(xi,yi,zi,x,xq,yq,zy,mode,extrapVal) % % input: -% xi,yi,zi: grid vectors +% xi,yi,zi: grid vectors % x: data % xq,yq,zq: coordinates of quer points as a grid % mode: optional interpolation mode (default linear) % extrapVal: (optional) value for extrapolation -% +% % output: -% y: interpolated data +% y: interpolated data % % References % - @@ -22,12 +22,12 @@ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Copyright 2019-2026 the matRad development team. -% -% This file is part of the matRad project. It is subject to the license -% terms in the LICENSE file found in the top-level directory of this -% distribution and at https://github.com/e0404/matRad/LICENSE.md. No part -% of the matRad project, including this file, may be copied, modified, -% propagated, or distributed except according to the terms contained in the +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSE.md. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the % LICENSE file. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -44,19 +44,22 @@ extrapVal = NaN; end -switch env - case 'MATLAB' - y = interp3(xi,yi,zi,x,xq,yq,zq,mode,extrapVal); - case 'OCTAVE' - %If we do a vector query with similar sizes only don't create a meshgrid - if isequal(size(xq),size(yq),size(zq)) +if isequal(xi,xq) && isequal(yi,yq) && isequal(zi,zq) + y = x; +else + switch env + case 'MATLAB' y = interp3(xi,yi,zi,x,xq,yq,zq,mode,extrapVal); - else - %Here we require a meshgrid to force octave to return the correct size - %Maybe the same thing could be achieved with a reshape? - [xqMesh,yqMesh,zqMesh] = meshgrid(xq,yq,zq); - y = interp3(xi,yi,zi,x,xqMesh,yqMesh,zqMesh,mode,extrapVal); - end - %Old implementation + case 'OCTAVE' + %If we do a vector query with similar sizes only don't create a meshgrid + if isequal(size(xq),size(yq),size(zq)) + y = interp3(xi,yi,zi,x,xq,yq,zq,mode,extrapVal); + else + %Here we require a meshgrid to force octave to return the correct size + %Maybe the same thing could be achieved with a reshape? + [xqMesh,yqMesh,zqMesh] = meshgrid(xq,yq,zq); + y = interp3(xi,yi,zi,x,xqMesh,yqMesh,zqMesh,mode,extrapVal); + end + end end diff --git a/test/doseCalc/test_TopasMCEngine.m b/test/doseCalc/test_TopasMCEngine.m index a19ce6660..41067f21f 100644 --- a/test/doseCalc/test_TopasMCEngine.m +++ b/test/doseCalc/test_TopasMCEngine.m @@ -322,7 +322,7 @@ assertTrue(isfile([folderName filesep 'matRad_cube.dat'])); assertTrue(isfile([folderName filesep 'matRad_cube.txt'])); assertTrue(isfile([folderName filesep 'MCparam.mat'])); - for j = 1:pln.propStf.numOfBeams + for j = 1:numel(pln.propStf.gantryAngles) assertTrue(isfile([folderName filesep 'beamSetup_matRad_plan_field' num2str(j) '.txt'])); assertTrue(isfile([folderName filesep 'matRad_plan_field' num2str(j) '_run1.txt'])); diff --git a/test/doseCalc/test_ignoreOutsideDensities.m b/test/doseCalc/test_ignoreOutsideDensities.m new file mode 100644 index 000000000..9aedf8502 --- /dev/null +++ b/test/doseCalc/test_ignoreOutsideDensities.m @@ -0,0 +1,25 @@ +function test_suite = test_ignoreOutsideDensities + +test_functions = localfunctions(); + +initTestSuite; + +function test_higherEnergiesInStf +testData = load('protons_testData.mat'); + +testData.ct.cubeHU{1}(1,:,:) = -500; + +stf_noIgnore = matRad_generateStf(testData.ct,testData.cst,testData.pln); + +testData.pln.propStf.ignoreOutsideDensities = true; + +stf_Ignore = matRad_generateStf(testData.ct,testData.cst,testData.pln); + +% check if energies for beam 1 are higher +for i = 1:size(stf_noIgnore(1).ray,2) + assertTrue( stf_Ignore(1).ray(i).energy(1)