Skip to content

Commit

Permalink
Merge pull request e0404#825 from remocristoforetti/dev_FRED_ijVersio…
Browse files Browse the repository at this point in the history 
…nReaderUpdate
  • Loading branch information
@wahln
wahln authored Feb 14, 2025
2 parents fdf2906 + 39f749d commit 786ffb6
Show file tree
Hide file tree
Showing 6 changed files with 257 additions and 27 deletions.
15 changes: 15 additions & 0 deletions matRad/basedata/matRad_MCemittanceBaseData.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -541,6 +541,21 @@
mcDataOptics.Divergence2y = 0;
mcDataOptics.Correlation2y = 0;
mcDataOptics.FWHMatIso = 2.355 * sigmaSqIso;

% Parameters for parametrization of sigma squared model.
% sigma^2 = a + b*z + c*z^2;
mcDataOptics.sSQ_a = (sigmaSqIso/10)^2;
mcDataOptics.sSQ_b = (2*rho*sigmaT*sigmaSqIso)/10;
mcDataOptics.sSQ_c = sigmaT^2;

% Parameters for emittance model
mcDataOptics.twissEpsilonX = sqrt(mcDataOptics.sSQ_a*mcDataOptics.sSQ_c - (mcDataOptics.sSQ_b^2)/4);
mcDataOptics.twissAlphaX = - mcDataOptics.sSQ_b/(2*mcDataOptics.twissEpsilonX);
mcDataOptics.twissBetaX = mcDataOptics.sSQ_a/mcDataOptics.twissEpsilonX;

mcDataOptics.twissEpsilonY = mcDataOptics.twissEpsilonX;
mcDataOptics.twissAlphaY = mcDataOptics.twissAlphaX;
mcDataOptics.twissBetaY = mcDataOptics.twissBetaX;
end


Expand Down
17 changes: 5 additions & 12 deletions matRad/doseCalc/+DoseEngines/@matRad_ParticleFREDEngine/calcDose.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,12 +86,6 @@
% Use the emittance base data class to recover MC information
emittanceBaseData = matRad_MCemittanceBaseData(this.machine,stf);

% if this.calcDoseDirect
% cumulativeWeights = 0;
% for i=1:length(stf)
% cumulativeWeights = cumulativeWeights + sum([stf(i).ray.weight]);
% end
% end

% Loop over fields. FRED performs one single simulation for multiple
% fields
Expand Down Expand Up @@ -162,7 +156,7 @@
stfFred(i).emittanceRefPlaneDistance = [];

% Need to get the parameters for the model from MCemittance
for eIdx=emittanceBaseData.energyIndex
for eIdx=emittanceBaseData.energyIndex'
% Only using first focus index for now
tmpOpticsData = emittanceBaseData.fitBeamOpticsForEnergy(eIdx,1);
stfFred(i).emittanceX = [stfFred(i).emittanceX, tmpOpticsData.twissEpsilonX];
Expand All @@ -176,7 +170,7 @@
stfFred(i).sSQr_b = [];
stfFred(i).sSQr_c = [];

for eIdx=emittanceBaseData.energyIndex
for eIdx=emittanceBaseData.energyIndex'
tmpOpticsData = emittanceBaseData.fitBeamOpticsForEnergy(eIdx,1);
stfFred(i).sSQr_a = [stfFred(i).sSQr_a, tmpOpticsData.sSQ_a];
stfFred(i).sSQr_b = [stfFred(i).sSQr_b, tmpOpticsData.sSQ_b];
Expand Down Expand Up @@ -272,7 +266,6 @@
for j=1:numel(stfFred(i).nominalEnergies)
stfFred(i).energyLayer(j).rayPosX = stfFred(i).energyLayer(j).rayPosX/10;
stfFred(i).energyLayer(j).rayPosY = stfFred(i).energyLayer(j).rayPosY/10;
%stfFred(i).energyLayer(j).targetPoints = stfFred(i).energyLayer(j).targetPoints/10;
stfFred(i).energyLayer(j).nBixels = numel(stfFred(i).energyLayer(j).bixelNum);

if this.calcDoseDirect
Expand Down Expand Up @@ -330,7 +323,7 @@
systemCall = [this.cmdCall, ' -nogpu -f fred.inp'];
end

% printOutput to matLab console
% printOutput to matlab console
if this.printOutput
[status,~] = system(systemCall,'-echo');
else
Expand All @@ -346,14 +339,14 @@
end

% read simulation output
[doseCube, letdCube] = this.readSimulationOutput(this.MCrunFolder,this.calcDoseDirect, logical(this.calcLET));
[doseCube, letdCube] = this.readSimulationOutput(this.MCrunFolder,this.calcDoseDirect, 'calcLET', logical(this.calcLET), 'readFunctionHandle', this.dijReaderHandle);

otherwise % A path for loading has been provided

matRad_cfg.dispInfo(['Reading simulation data from: ', strrep(this.MCrunFolder,'\','\\'), '\n']);

% read simulation output
[doseCube, letdCube, loadFileName] = this.readSimulationOutput(this.MCrunFolder,this.calcDoseDirect, logical(this.calcLET));
[doseCube, letdCube, loadFileName] = this.readSimulationOutput(this.MCrunFolder,this.calcDoseDirect, 'calcLET',logical(this.calcLET),'readFunctionHandle', this.dijReaderHandle);

dij.externalCalculationLodPath = loadFileName;

Expand Down
207 changes: 200 additions & 7 deletions matRad/doseCalc/+DoseEngines/@matRad_ParticleFREDEngine/matRad_ParticleFREDEngine.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,15 +54,17 @@

defaultHUtable = 'matRad_default_FredMaterialConverter';
AvailableSourceModels = {'gaussian', 'emittance', 'sigmaSqrModel'};
defaultDijFormatVersion = '20';

calcBioDose;
currentVersion;
availableVersions = {'3.70.0'}; % Interface requires latest FRED version
availableVersions = {'3.70.0'}; % Or higher.
radiationMode;

end

properties
dijFormatVersion;
externalCalculation;
useGPU;
calcLET;
Expand Down Expand Up @@ -100,6 +102,7 @@
inputFolder;
regionsFolder;
planFolder;
dijReaderHandle;
end

methods
Expand Down Expand Up @@ -316,6 +319,7 @@ function setDefaults(this)
this.outputMCvariance = false;
this.constantRBE = NaN;
this.ignoreOutsideDensities = false;
this.dijFormatVersion = this.defaultDijFormatVersion;

end

Expand Down Expand Up @@ -473,6 +477,7 @@ function writeHlut(this,hLutFile)
end
catch
matRad_cfg.dispWarning('Something wrong occured in checking FRED installation. Please check correct FRED installation');
version = [];
end
end

Expand Down Expand Up @@ -517,14 +522,70 @@ function writeHlut(this,hLutFile)
numVox = fread(f,1,"int32");

bixelCounter = bixelCounter +1;
% FRED adds 1000000 when new field is added
%bixNum = bixNum - (10^6*(i-1));

colIndices(end+1:end+numVox) = bixelCounter;
currVoxelIndices = fread(f,numVox,"uint32") + 1;
tmpValues = fread(f,numVox*nComponents,"float32");
valuesNom = tmpValues(1:nComponents:end);

if nComponents == 2
valuesDen = tmpValues(nComponents:nComponents:end);
values(end+1:end+numVox) = valuesNom./valuesDen;
else
values(end+1:end+numVox) = valuesNom;
end

% x and y components have been permuted in CT
[indY, indX, indZ] = ind2sub(dims, currVoxelIndices);

voxelIndices(end+1:end+numVox) = sub2ind(dims([2,1,3]), indX, indY, indZ);
matRad_cfg.dispInfo("\tRead beamlet %d, %d voxels...\n",bixNum,numVox);
end
dijMatrix = sparse(voxelIndices,colIndices,values,prod(dims),numberOfBixels);

fclose(f);
catch
fclose(f);
matRad_cfg.dispError('unable to load file: %s',fName);
end
end

function dijMatrix = readSparseDijBin_v21(fName)

matRad_cfg = MatRad_Config.instance();

f = fopen(fName,'r','l');

try
%Header
fileFormatVerison = fread(f,1,"int32");
dims = fread(f,3,"int32");
res = fread(f,3,"float32");
offset = fread(f,3,"float32");
orientation = fread(f,9,"float32");
nComponents = fread(f,1,"int32");
numberOfBixels = fread(f,1,"int32");

values = [];
valuesDen = [];
voxelIndices = [];
colIndices = [];
valuesNom = [];

matRad_cfg.dispInfo("Reading %d number of beamlets in %d voxels (%dx%dx%d)\n",numberOfBixels,prod(dims),dims(1),dims(2),dims(3));

bixelCounter = 0;
for i = 1:numberOfBixels
%Read Beamlet
bixNum = fread(f,1,"uint32");
numVox = fread(f,1,"int32");

bixelCounter = bixelCounter +1;

colIndices(end+1:end+numVox) = bixelCounter; %bixNum + 1;
colIndices(end+1:end+numVox) = bixelCounter;
currVoxelIndices = fread(f,numVox,"uint32") + 1;
tmpValues = fread(f,numVox*nComponents,"float32");
valuesNom = tmpValues(1:nComponents:end);%tmpValues(nComponents:nComponents:end);
% values(end+1:end+numVox) = tmpValuess(1:nComponents:end);%tmpValues(nComponents:nComponents:end);
valuesNom = tmpValues(1:nComponents:end);

if nComponents == 2
valuesDen = tmpValues(nComponents:nComponents:end);
Expand All @@ -546,12 +607,78 @@ function writeHlut(this,hLutFile)
fclose(f);
matRad_cfg.dispError('unable to load file: %s',fName);
end
end

function dijMatrix = readSparseDijBin_v31(fName)

matRad_cfg = MatRad_Config.instance();

f = fopen(fName,'r','l');

try
fileFormatVerison = fread(f,1,"int32");
dims = fread(f,3,"int32");
res = fread(f,3,"float32");
offset = fread(f,3,"float32");
orientation = fread(f,9,"float32");
nComponents = fread(f,1,"uint32");
numberOfBixels = fread(f,1,"uint32");

values = [];
valuesDen = [];
voxelIndices = [];
colIndices = [];
valuesNom = [];

matRad_cfg.dispInfo("Reading %d number of beamlets in %d voxels (%dx%dx%d)\n",numberOfBixels,prod(dims),dims(1),dims(2),dims(3));

allBixelMeta = fread(f, 3*numberOfBixels, "uint32");
allBixelMeta = reshape(allBixelMeta, 3,numberOfBixels)'; % (#bixels, PBidx, FID, PBID)

% if nComponents>1
% matRad_cfg.dispWarning('!! Only last component will be read!!')
% end
% Components header
for i = 1:nComponents
componentDataSize(i) = fread(f,1,"uint32");
end

% Data
for compIdx=1:nComponents
% PBidx and voxelIndices should be always the same for
% each component?
PBidxs = fread(f, componentDataSize(compIdx), "uint32")+1;
voxelIndices = fread(f, componentDataSize(compIdx), "uint32")+1;
tmpValues{compIdx} = fread(f, componentDataSize(compIdx), "float32");

end

% For now we only have Dose and LET scorers, if nComponents
% == 2, then it's LET

if nComponents>1
values = tmpValues{1}./tmpValues{2};
else
values = tmpValues{1};
end

% x and y components have been permuted in CT
[indY, indX, indZ] = ind2sub(dims, voxelIndices);

voxelIndices = sub2ind(dims([2,1,3]), indX, indY, indZ); % + (nComponents-1)*prod(dims);
% This will probably not work for multiple components
dijMatrix = sparse(voxelIndices,PBidxs,values,prod(dims), numberOfBixels);

fclose(f);

catch
fclose(f);
matRad_cfg.dispError('unable to load file: %s',fName);

end
end

[doseCubeV, letdCubeV, fileName] = readSimulationOutput(runFolder,calcDoseDirect, calcLET);
[doseCubeV, letdCubeV, fileName] = readSimulationOutput(runFolder,calcDoseDirect, varargin);

end

Expand Down Expand Up @@ -595,6 +722,25 @@ function updatePaths(obj, rootFolder)

matRad_cfg.dispWarning('Selected radiation modality: %s with primary mass: %2.3f', radiationMode, this.primaryMass);
end

function isHigher = isVersionHigher(this,version)
isHigher = false;

% This function directly looks at FRED installation, not at
% the current FRED version stored in the class property.
fredVersion = this.getVersion();

if ~isempty(fredVersion)
% Decompose the current version for comparison
v1 = sscanf(fredVersion, '%d.%d.%d')';
v2 = sscanf(version, '%d.%d.%d')';

if (v1(1) >= v2(1)) && (v1(2) >= v2(2)) && (v1(3) > v2(3))
isHigher = true;
end
end

end

end

Expand Down Expand Up @@ -626,6 +772,53 @@ function updatePaths(obj, rootFolder)

end

function set.dijFormatVersion(this,value)

matRad_cfg = MatRad_Config.instance();

if ~this.isVersionHigher('3.70.0')
% FRED version < 3.70.0 does not allow dij version
% selection and only works with ifFormatVersion < 21

this.dijFormatVersion = this.defaultDijFormatVersion;

if ~strcmp(value, this.defaultDijFormatVersion)
matRad_cfg.dispWarning(sprintf('ijFormat: %s not available for FRED version < 3.70.0', value));
end

else
%if this.useGPU
if strcmp(value, '20')
this.dijFormatVersion = '21';
else
this.dijFormatVersion = value;
end
%else
% this.dijFormatVersion = '20';
%end
end
end

function readerHandle = get.dijReaderHandle(this)

matRad_cfg = MatRad_Config.instance();

switch this.dijFormatVersion

case '20'
readerHandle = @(lFile) DoseEngines.matRad_ParticleFREDEngine.readSparseDijBin(lFile);
case '21'
readerHandle = @(lFile) DoseEngines.matRad_ParticleFREDEngine.readSparseDijBin_v21(lFile);
case '31'
readerHandle = @(lFile) DoseEngines.matRad_ParticleFREDEngine.readSparseDijBin_v31(lFile);
otherwise
matRad_cfg.dispWarning(sprintf('Unable to read dij format version: %s, using default: %s', this.dijFormatVersion, this.defaultDijFormatVersion));
readerHandle = @(lFile) DoseEngines.matRad_ParticleFREDEngine.readSparseDijBin(lFile);
end


end

function set.radiationMode(this,value)
if ischar(value)
if ~isempty(this.radiationMode) && ~strcmp(this.radiationMode, value)
Expand Down
Loading

0 comments on commit 786ffb6

Please sign in to comment.