From 3da02de78630af8f9fb925c58bde878e5cbd569a Mon Sep 17 00:00:00 2001 From: JenHardt <84377305+JenHardt@users.noreply.github.com> Date: Mon, 2 Feb 2026 16:33:04 +0100 Subject: [PATCH 1/2] Overworked 3d conformal and sequencing 3d conformal now no longer errors out when used together with sequencing, the issue here was the inconsonant number of bixels in the stf, which is why the collapseStf function was added. Also changes to the sequencer were made to accomondate for example a beam with 0 weight vector. --- matRad/gui/widgets/matRad_WorkflowWidget.m | 2 + matRad/optimization/matRad_collapseDij.m | 21 ++++--- matRad/optimization/matRad_collapseStf.m | 36 ++++++++++++ .../sequencing/matRad_engelLeafSequencing.m | 37 ++++++++---- .../sequencing/matRad_siochiLeafSequencing.m | 56 ++++++++++++------- matRad/sequencing/matRad_xiaLeafSequencing.m | 38 +++++++++---- 6 files changed, 137 insertions(+), 53 deletions(-) create mode 100644 matRad/optimization/matRad_collapseStf.m diff --git a/matRad/gui/widgets/matRad_WorkflowWidget.m b/matRad/gui/widgets/matRad_WorkflowWidget.m index e66759b76..25f915940 100644 --- a/matRad/gui/widgets/matRad_WorkflowWidget.m +++ b/matRad/gui/widgets/matRad_WorkflowWidget.m @@ -424,9 +424,11 @@ function btnCalcDose_Callback(this, hObject, eventdata) % prepare dij for 3d conformal if isfield(pln,'propOpt') && isfield(pln.propOpt,'conf3D') && pln.propOpt.conf3D dij = matRad_collapseDij(dij); + stf = matRad_collapseStf(stf); end % assign results to base worksapce assignin('base','dij',dij); + assignin('base','stf',stf); catch ME diff --git a/matRad/optimization/matRad_collapseDij.m b/matRad/optimization/matRad_collapseDij.m index 03c66895e..945a8b639 100644 --- a/matRad/optimization/matRad_collapseDij.m +++ b/matRad/optimization/matRad_collapseDij.m @@ -28,14 +28,14 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -dijNew.totalNumOfBixels = 1; -dijNew.totalNumOfRays = 1; -dijNew.numOfBeams = 1; -dijNew.numOfRaysPerBeam = 1; +dijNew.totalNumOfBixels = dij.numOfBeams; +dijNew.totalNumOfRays = dij.numOfBeams; +dijNew.numOfBeams = dij.numOfBeams; +dijNew.numOfRaysPerBeam = ones(dij.numOfBeams,1); -dijNew.beamNum = 1; -dijNew.bixelNum = 1; -dijNew.rayNum = 1; +dijNew.beamNum = [1:dij.numOfBeams]'; +dijNew.bixelNum = [1:dij.numOfBeams]'; +dijNew.rayNum = [1:dij.numOfBeams]'; dijNew.doseGrid = dij.doseGrid; dijNew.ctGrid = dij.ctGrid; @@ -43,6 +43,11 @@ dijNew.numOfScenarios = dij.numOfScenarios; for i = 1:dij.numOfScenarios - dijNew.physicalDose{i} = sum(dij.physicalDose{i},2); + tmp = sparse(dij.doseGrid.numOfVoxels,dij.numOfBeams); % initialize sparse matrix + for j = 1:dij.numOfBeams + % Sum only the columns corresponding to beam j + tmp(:, j) = sum(dij.physicalDose{i}(:, dij.beamNum == j), 2); + end + dijNew.physicalDose{i} = tmp; end diff --git a/matRad/optimization/matRad_collapseStf.m b/matRad/optimization/matRad_collapseStf.m new file mode 100644 index 000000000..7213d6bdf --- /dev/null +++ b/matRad/optimization/matRad_collapseStf.m @@ -0,0 +1,36 @@ +function stf = matRad_collapseStf(stf) +% matRad collapse dij function for simulation of 3D conformal treatments. +% Function to supress intensity-modulation for photons in order to simulate +% 3D conformal treatments. +% +% call +% dijNew = matRad_collapseDij(dij) +% +% input +% dij: dose influence matrix +% +% output +% dijNew: collapsed dose influence matrix +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2018 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 +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% dummy collapse so that it works with sequencing + +for i = 1:size(stf,2) + stf(i).numOfRays = 1; + stf(i).totalNumOfBixels = 1; +end \ No newline at end of file diff --git a/matRad/sequencing/matRad_engelLeafSequencing.m b/matRad/sequencing/matRad_engelLeafSequencing.m index a1d944374..de6462632 100644 --- a/matRad/sequencing/matRad_engelLeafSequencing.m +++ b/matRad/sequencing/matRad_engelLeafSequencing.m @@ -57,12 +57,12 @@ numOfRaysPerBeam = stf(i).numOfRays; % get relevant weights for current beam - wOfCurrBeams = resultGUI.w(1+offset:numOfRaysPerBeam+offset); + wOfCurrBeams = resultGUI.w(1+offset:numOfRaysPerBeam+offset).* ones(size(stf(i).ray,2),1); - X = ones(numOfRaysPerBeam,1)*NaN; - Z = ones(numOfRaysPerBeam,1)*NaN; + X = ones(size(stf(i).ray,2),1)*NaN; %this way it also works with3dconformal + Z = ones(size(stf(i).ray,2),1)*NaN; - for j=1:stf(i).numOfRays + for j=1:size(stf(i).ray,2) X(j) = stf(i).ray(j).rayPos_bev(:,1); Z(j) = stf(i).ray(j).rayPos_bev(:,3); end @@ -352,15 +352,28 @@ clear rightIntLimit; end - - sequencing.beam(i).numOfShapes = k; - sequencing.beam(i).shapes = shapes(:,:,1:k); - sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac; - sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; - sequencing.beam(i).fluence = D_0; - - sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = D_0(indInFluenceMx)/numOfLevels*calFac; + if sum(wOfCurrBeams)>0 + + sequencing.beam(i).numOfShapes = k; + sequencing.beam(i).shapes = shapes(:,:,1:k); + sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac; + sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; + sequencing.beam(i).fluence = D_0; + else + sequencing.beam(i).numOfShapes = 1; + sequencing.beam(i).shapes = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).shapesWeight = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; + sequencing.beam(i).fluence = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).sum = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + end + + if numOfRaysPerBeam >1 + sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = sequencing.beam(i).sum(indInFluenceMx); + else + sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = wOfCurrBeams(1); + end offset = offset + numOfRaysPerBeam; end diff --git a/matRad/sequencing/matRad_siochiLeafSequencing.m b/matRad/sequencing/matRad_siochiLeafSequencing.m index 4240a6c9a..02841d5b0 100644 --- a/matRad/sequencing/matRad_siochiLeafSequencing.m +++ b/matRad/sequencing/matRad_siochiLeafSequencing.m @@ -70,12 +70,12 @@ numOfRaysPerBeam = stf(i).numOfRays; % get relevant weights for current beam - wOfCurrBeams = wUnsequenced(1+offset:numOfRaysPerBeam+offset);%REVIEW OFFSET + wOfCurrBeams = wUnsequenced(1+offset:numOfRaysPerBeam+offset).* ones(size(stf(i).ray,2),1);%REVIEW OFFSET - X = ones(numOfRaysPerBeam,1)*NaN; - Z = ones(numOfRaysPerBeam,1)*NaN; + X = ones(size(stf(i).ray,2),1)*NaN; %this way it also works with3dconformal + Z = ones(size(stf(i).ray,2),1)*NaN; - for j = 1:stf(i).numOfRays + for j = 1:size(stf(i).ray,2) X(j) = stf(i).ray(j).rayPos_bev(:,1); Z(j) = stf(i).ray(j).rayPos_bev(:,3); end @@ -137,27 +137,41 @@ D_k_MinX = min(D_k_X); D_k_MaxX = max(D_k_X); - %Decompose the port, do rod pushing - [tops, bases] = matRad_siochiDecomposePort(D_k,dimOfFluenceMxZ,dimOfFluenceMxX,D_k_MinZ,D_k_MaxZ,D_k_MinX,D_k_MaxX); - %Form segments with and without visualization - if visBool - [shapes,shapesWeight,k,D_k]=matRad_siochiConvertToSegments(shapes,shapesWeight,k,tops,bases,visBool,i,D_k,numOfLevels,seqFig,seqSubPlots); + if sum(wOfCurrBeams)>0 + %Decompose the port, do rod pushing + [tops, bases] = matRad_siochiDecomposePort(D_k,dimOfFluenceMxZ,dimOfFluenceMxX,D_k_MinZ,D_k_MaxZ,D_k_MinX,D_k_MaxX); + %Form segments with and without visualization + if visBool + [shapes,shapesWeight,k,D_k]=matRad_siochiConvertToSegments(shapes,shapesWeight,k,tops,bases,visBool,i,D_k,numOfLevels,seqFig,seqSubPlots); + else + [shapes,shapesWeight,k]=matRad_siochiConvertToSegments(shapes,shapesWeight,k,tops,bases); + end + + sequencing.beam(i).numOfShapes = k; + sequencing.beam(i).shapes = shapes(:,:,1:k); + sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac; + sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; + sequencing.beam(i).fluence = D_0; + sequencing.beam(i).sum = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + + for j = 1:k + sequencing.beam(i).sum = sequencing.beam(i).sum+sequencing.beam(i).shapes(:,:,j)*sequencing.beam(i).shapesWeight(j); + end + else - [shapes,shapesWeight,k]=matRad_siochiConvertToSegments(shapes,shapesWeight,k,tops,bases); + sequencing.beam(i).numOfShapes = 1; + sequencing.beam(i).shapes = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).shapesWeight = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; + sequencing.beam(i).fluence = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).sum = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); end - sequencing.beam(i).numOfShapes = k; - sequencing.beam(i).shapes = shapes(:,:,1:k); - sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac; - sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; - sequencing.beam(i).fluence = D_0; - sequencing.beam(i).sum = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); - - for j = 1:k - sequencing.beam(i).sum = sequencing.beam(i).sum+sequencing.beam(i).shapes(:,:,j)*sequencing.beam(i).shapesWeight(j); + if numOfRaysPerBeam >1 + sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = sequencing.beam(i).sum(indInFluenceMx); + else + sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = wOfCurrBeams(1); end - sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = sequencing.beam(i).sum(indInFluenceMx); - offset = offset + numOfRaysPerBeam; end diff --git a/matRad/sequencing/matRad_xiaLeafSequencing.m b/matRad/sequencing/matRad_xiaLeafSequencing.m index aaa893324..c26d618ea 100644 --- a/matRad/sequencing/matRad_xiaLeafSequencing.m +++ b/matRad/sequencing/matRad_xiaLeafSequencing.m @@ -60,12 +60,12 @@ numOfRaysPerBeam = stf(i).numOfRays; % get relevant weights for current beam - wOfCurrBeams = resultGUI.w(1+offset:numOfRaysPerBeam+offset);%REVIEW OFFSET + wOfCurrBeams = resultGUI.w(1+offset:numOfRaysPerBeam+offset).* ones(size(stf(i).ray,2),1);%;%REVIEW OFFSET - X = ones(numOfRaysPerBeam,1)*NaN; - Z = ones(numOfRaysPerBeam,1)*NaN; - - for j=1:stf(i).numOfRays + X = ones(size(stf(i).ray,2),1)*NaN; %this way it also works with3dconformal + Z = ones(size(stf(i).ray,2),1)*NaN; + + for j=1:size(stf(i).ray,2) X(j) = stf(i).ray(j).rayPos_bev(:,1); Z(j) = stf(i).ray(j).rayPos_bev(:,3); end @@ -239,15 +239,29 @@ L_k = max(D_k(:)); % eq 5 end + + if sum(wOfCurrBeams)>0 - sequencing.beam(i).numOfShapes = k; - sequencing.beam(i).shapes = shapes(:,:,1:k); - sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac; - sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; - sequencing.beam(i).fluence = D_0; + sequencing.beam(i).numOfShapes = k; + sequencing.beam(i).shapes = shapes(:,:,1:k); + sequencing.beam(i).shapesWeight = shapesWeight(1:k)/numOfLevels*calFac; + sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; + sequencing.beam(i).fluence = D_0; - sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = D_0(indInFluenceMx)/numOfLevels*calFac; - + else + sequencing.beam(i).numOfShapes = 1; + sequencing.beam(i).shapes = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).shapesWeight = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; + sequencing.beam(i).fluence = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + sequencing.beam(i).sum = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); + end + + if numOfRaysPerBeam >1 + sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = sequencing.beam(i).sum(indInFluenceMx); + else + sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = wOfCurrBeams(1); + end offset = offset + numOfRaysPerBeam; end From e9299e74056195310f19ecf267725f53e85ca479 Mon Sep 17 00:00:00 2001 From: JenHardt <84377305+JenHardt@users.noreply.github.com> Date: Tue, 3 Feb 2026 09:23:14 +0100 Subject: [PATCH 2/2] removing of sum --- matRad/sequencing/matRad_engelLeafSequencing.m | 8 ++------ matRad/sequencing/matRad_xiaLeafSequencing.m | 9 +++------ 2 files changed, 5 insertions(+), 12 deletions(-) diff --git a/matRad/sequencing/matRad_engelLeafSequencing.m b/matRad/sequencing/matRad_engelLeafSequencing.m index de6462632..bc1436574 100644 --- a/matRad/sequencing/matRad_engelLeafSequencing.m +++ b/matRad/sequencing/matRad_engelLeafSequencing.m @@ -366,14 +366,10 @@ sequencing.beam(i).shapesWeight = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; sequencing.beam(i).fluence = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); - sequencing.beam(i).sum = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); end - if numOfRaysPerBeam >1 - sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = sequencing.beam(i).sum(indInFluenceMx); - else - sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = wOfCurrBeams(1); - end + sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = D_0(indInFluenceMx)/numOfLevels*calFac; + offset = offset + numOfRaysPerBeam; end diff --git a/matRad/sequencing/matRad_xiaLeafSequencing.m b/matRad/sequencing/matRad_xiaLeafSequencing.m index c26d618ea..052c918c6 100644 --- a/matRad/sequencing/matRad_xiaLeafSequencing.m +++ b/matRad/sequencing/matRad_xiaLeafSequencing.m @@ -254,14 +254,11 @@ sequencing.beam(i).shapesWeight = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); sequencing.beam(i).bixelIx = 1+offset:numOfRaysPerBeam+offset; sequencing.beam(i).fluence = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); - sequencing.beam(i).sum = zeros(dimOfFluenceMxZ,dimOfFluenceMxX); end - if numOfRaysPerBeam >1 - sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = sequencing.beam(i).sum(indInFluenceMx); - else - sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = wOfCurrBeams(1); - end + + sequencing.w(1+offset:numOfRaysPerBeam+offset,1) = D_0(indInFluenceMx)/numOfLevels*calFac; + offset = offset + numOfRaysPerBeam; end