From da6f668748b005f64bb899663c383607655e34bc Mon Sep 17 00:00:00 2001 From: amitantony Date: Mon, 22 Jul 2024 09:59:45 +0200 Subject: [PATCH 01/10] basic world coordinate functions for stf --- matRad/geometry/matRad_cube2worldCoords.m | 42 +++++++++++++++++++++++ matRad/geometry/matRad_getIsoCenter.m | 13 +++---- matRad/geometry/matRad_getWorldAxes.m | 40 +++++++++++++++++++++ matRad/geometry/matRad_world2cubeCoords.m | 38 ++++++++++++++++++++ matRad/matRad_generateStf.m | 28 ++++++++++----- 5 files changed, 146 insertions(+), 15 deletions(-) create mode 100644 matRad/geometry/matRad_cube2worldCoords.m create mode 100644 matRad/geometry/matRad_getWorldAxes.m create mode 100644 matRad/geometry/matRad_world2cubeCoords.m diff --git a/matRad/geometry/matRad_cube2worldCoords.m b/matRad/geometry/matRad_cube2worldCoords.m new file mode 100644 index 000000000..b1cd6e787 --- /dev/null +++ b/matRad/geometry/matRad_cube2worldCoords.m @@ -0,0 +1,42 @@ +function coord = matRad_cube2worldCoords(vCoord, ct) +% matRad function to convert cube coordinates to world coordinates +% +% call +% coord = matRad_worldToCubeCoordinates(vCoord, ct) +% +% ct: matRad ct struct +% vCoord : Voxel Coordinates [vx vy vz] +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 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/LICENSES.txt. 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +coord = []; % [x y z ] mm +if nargin == 2 + if isfield(ct, 'dicomInfo') && isfield(ct.dicomInfo,'ImagePositionPatient') + firstVox = ct.dicomInfo.ImagePositionPatient; + else + firstVox = - (ct.cubeDim./2).*[ct.resolution.x ct.resolution.y ct.resolution.z] ; + end + try + coord(1,:) = firstVox(1) + (vCoord(:,1) - 1 ) *ct.resolution.x; + coord(2,:) = firstVox(2) + (vCoord(:,2) - 1 ) *ct.resolution.y; + coord(3,:) = firstVox(3) + (vCoord(:,3) - 1 ) *ct.resolution.z; + coord = coord'; + catch + + end +end + +end \ No newline at end of file diff --git a/matRad/geometry/matRad_getIsoCenter.m b/matRad/geometry/matRad_getIsoCenter.m index 1b595e02b..bff249216 100644 --- a/matRad/geometry/matRad_getIsoCenter.m +++ b/matRad/geometry/matRad_getIsoCenter.m @@ -65,15 +65,16 @@ end % Transform subcripts from linear indices -[yCoordsV, xCoordsV, zCoordsV] = ind2sub(ct.cubeDim,V); +[coordV(:,2), coordV(:,1), coordV(:,3)] = ind2sub(ct.cubeDim,V); + % Transform to [mm] -xCoordsV = xCoordsV * ct.resolution.x; -yCoordsV = yCoordsV * ct.resolution.y; -zCoordsV = zCoordsV * ct.resolution.z; +coord = matRad_cube2worldCoords(coordV, ct); %idx2worldcoord + % Calculated isocenter. -isoCenter = mean([xCoordsV yCoordsV zCoordsV]); +isoCenter = mean(coord); + % Visualization if visBool @@ -82,7 +83,7 @@ hold on % Plot target - plot3(yCoordsV,xCoordsV,zCoordsV,'kx') + plot3(coord(:,2), coord(:,1), coord(:,3),'kx') % Show isocenter: red point plot3(isoCenter(2),isoCenter(1),isoCenter(3),'r.','MarkerSize',30) diff --git a/matRad/geometry/matRad_getWorldAxes.m b/matRad/geometry/matRad_getWorldAxes.m new file mode 100644 index 000000000..7577f41a5 --- /dev/null +++ b/matRad/geometry/matRad_getWorldAxes.m @@ -0,0 +1,40 @@ +function ct = matRad_getWorldAxes(ct) +% matRad function to compute and store world coordinates int ct.x +% +% call +% ct = matRad_getWorldAxes(ct) +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 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/LICENSES.txt. 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +matRad_cfg = MatRad_Config.instance(); + +if nargin>0 + % + % check if dicominfo exists + if isfield(ct, 'dicomInfo') && isfield(ct.dicomInfo,'ImagePositionPatient') + firstVox = ct.dicomInfo.ImagePositionPatient; + else + firstVox = - (ct.cubeDim./2).*[ct.resolution.x ct.resolution.y ct.resolution.z] ; + end + + ct.x = firstVox(1) + ct.resolution.x*[0:ct.cubeDim(2)-1] ; + ct.y = firstVox(2) + ct.resolution.y*[0:ct.cubeDim(1)-1] ; + ct.z = firstVox(3) + ct.resolution.z*[0:ct.cubeDim(3)-1] ; +else + matRad_cfg.dispWarning('Cannot compute world coordinates without matRad ct structure'); +end + +end \ No newline at end of file diff --git a/matRad/geometry/matRad_world2cubeCoords.m b/matRad/geometry/matRad_world2cubeCoords.m new file mode 100644 index 000000000..a652fb949 --- /dev/null +++ b/matRad/geometry/matRad_world2cubeCoords.m @@ -0,0 +1,38 @@ +function coord = matRad_world2cubeCoords(wCoord, ct) +% matRad function to convert world coordinates to cube coordinates +% +% call +% coord = world2cubeCoords(wCoord, ct) +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 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/LICENSES.txt. 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +coord = []; +if nargin == 2 + try + if ~ (isfield(ct,'x') && isfield(ct,'y') && isfield(ct,'z') ) + ct = matRad_getWorldAxes(ct); + end + + [~,coord(2)]=min(abs(ct.y-wCoord(2))); + [~,coord(1)]=min(abs(ct.x-wCoord(1))); + [~,coord(3)]=min(abs(ct.z-wCoord(3))); + + catch + + end +end + +end diff --git a/matRad/matRad_generateStf.m b/matRad/matRad_generateStf.m index a7c7ddacf..f6742f7b1 100644 --- a/matRad/matRad_generateStf.m +++ b/matRad/matRad_generateStf.m @@ -192,6 +192,11 @@ matRad_cfg.dispError('Could not find target.'); end +% get world coordinate system +if ~isfield(ct,'x') + ct = matRad_getWorldAxes(ct); +end + % Convert linear indices to 3D voxel coordinates [coordsY_vox, coordsX_vox, coordsZ_vox] = ind2sub(ct.cubeDim,V); @@ -218,9 +223,11 @@ % Correct for iso center position. Whit this correction Isocenter is % (0,0,0) [mm] - coordsX = coordsX_vox*ct.resolution.x - pln.propStf.isoCenter(i,1); - coordsY = coordsY_vox*ct.resolution.y - pln.propStf.isoCenter(i,2); - coordsZ = coordsZ_vox*ct.resolution.z - pln.propStf.isoCenter(i,3); + wCoords = matRad_cube2worldCoords([coordsX_vox, coordsY_vox, coordsZ_vox], ct); + + coordsX = wCoords(:,1) - pln.propStf.isoCenter(i,1); + coordsY = wCoords(:,2) - pln.propStf.isoCenter(i,2); + coordsZ = wCoords(:,3) - pln.propStf.isoCenter(i,3); % Save meta information for treatment plan stf(i).gantryAngle = pln.propStf.gantryAngles(i); @@ -325,12 +332,14 @@ % loop over all rays to determine meta information for each ray stf(i).numOfBixelsPerRay = ones(1,stf(i).numOfRays); - + + % mm axes with isocenter at (0,0,0) + mmCubeisoCenter = stf(i).isoCenter - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z]; for j = stf(i).numOfRays:-1:1 - + for ShiftScen = 1:pln.multScen.totNumShiftScen % ray tracing necessary to determine depth of the target - [alphas,l{ShiftScen},rho{ShiftScen},d12,~] = matRad_siddonRayTracer(stf(i).isoCenter + pln.multScen.isoShift(ShiftScen,:), ... + [alphas,l{ShiftScen},rho{ShiftScen},d12,~] = matRad_siddonRayTracer(mmCubeisoCenter + pln.multScen.isoShift(ShiftScen,:), ... ct.resolution, ... stf(i).sourcePoint, ... stf(i).ray(j).targetPoint, ... @@ -591,9 +600,10 @@ % generate a 3D rectangular grid centered at isocenter in % voxel coordinates - [X,Y,Z] = meshgrid((1:ct.cubeDim(2))-stf(i).isoCenter(1)/ct.resolution.x, ... - (1:ct.cubeDim(1))-stf(i).isoCenter(2)/ct.resolution.y, ... - (1:ct.cubeDim(3))-stf(i).isoCenter(3)/ct.resolution.z); + cubeIso = matRad_world2cubeCoords(stf(i).isoCenter,ct); + [X,Y,Z] = meshgrid( (1:ct.cubeDim(2))- cubeIso(1), ... + (1:ct.cubeDim(1))- cubeIso(2), ... + (1:ct.cubeDim(3))- cubeIso(3)); % computes surface patSurfCube = 0*ct.cube{1}; From a47d27deb0e86781772459079750c56c48475f6b Mon Sep 17 00:00:00 2001 From: amitantony Date: Mon, 22 Jul 2024 11:37:35 +0200 Subject: [PATCH 02/10] Change to world coordinate system for the isocenter --- examples/matRad_example10_4DphotonRobust.m | 3 +- examples/matRad_example11_helium.m | 3 +- examples/matRad_example15_brachy.m | 3 +- examples/matRad_example16_photonMC_MLC.m | 3 +- examples/matRad_example1_phantom.m | 3 +- examples/matRad_example2_photons.m | 3 +- examples/matRad_example4_photonsMC.m | 3 +- examples/matRad_example5_protons.m | 3 +- examples/matRad_example6_protonsNoise.m | 3 +- examples/matRad_example7_carbon.m | 6 +- examples/matRad_example8_protonsRobust.m | 3 +- examples/matRad_example9_4DDoseCalcMinimal.m | 3 +- .../matRad_exportDicomCt.m | 9 +- .../matRad_exportDicomRTDoses.m | 9 +- .../matRad_exportDicomRTStruct.m | 9 +- .../@matRad_DoseEngineBase/initDoseCalc.m | 11 +-- .../matRad_ParticleMCsquareEngine.m | 7 +- .../matRad_PencilBeamEngineAbstract.m | 7 +- .../+DoseEngines/matRad_PhotonOmpMCEngine.m | 21 ++-- .../+DoseEngines/matRad_TopasMCEngine.m | 2 +- matRad/geometry/matRad_projectOnComponents.m | 2 +- matRad/gui/widgets/matRad_ViewingWidget.m | 14 +-- matRad/planAnalysis/matRad_compareDose.m | 97 ++++++++++--------- matRad/plotting/matRad_plotAxisLabels.m | 40 +++++--- matRad/plotting/matRad_plotIsoCenterMarker.m | 6 +- .../matRad_plotProjectedGantryAngles.m | 19 ++-- 26 files changed, 165 insertions(+), 127 deletions(-) diff --git a/examples/matRad_example10_4DphotonRobust.m b/examples/matRad_example10_4DphotonRobust.m index 959fa9748..5558f0815 100644 --- a/examples/matRad_example10_4DphotonRobust.m +++ b/examples/matRad_example10_4DphotonRobust.m @@ -248,7 +248,8 @@ %% Visualize results plane = 3; -slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); maxDose = max([max(resultGUI.([quantityOpt])(:,:,slice)) max(resultGUIrobust.([quantityOpt])(:,:,slice))])+1e-4; doseIsoLevels = linspace(0.1 * maxDose,maxDose,10); figure, diff --git a/examples/matRad_example11_helium.m b/examples/matRad_example11_helium.m index e28ded925..b019114a3 100644 --- a/examples/matRad_example11_helium.m +++ b/examples/matRad_example11_helium.m @@ -93,7 +93,8 @@ %% Plot the Resulting Dose Slice % Let's plot the transversal iso-center dose slice -slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); figure subplot(121),imagesc(resultGUI.physicalDose(:,:,slice)),colorbar,colormap(jet),title('physical dose') subplot(122),imagesc(resultGUI.RBExD(:,:,slice)),colorbar,colormap(jet),title('RBExDose') diff --git a/examples/matRad_example15_brachy.m b/examples/matRad_example15_brachy.m index be9c378cf..58a1bfe37 100644 --- a/examples/matRad_example15_brachy.m +++ b/examples/matRad_example15_brachy.m @@ -226,7 +226,8 @@ %% IV.1 Plot the Resulting Dose Slice % Let's plot the transversal iso-center dose slice -slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); figure imagesc(resultGUI.physicalDose(:,:,slice)),colorbar, colormap(jet); diff --git a/examples/matRad_example16_photonMC_MLC.m b/examples/matRad_example16_photonMC_MLC.m index f5e17a6e4..d8658af19 100644 --- a/examples/matRad_example16_photonMC_MLC.m +++ b/examples/matRad_example16_photonMC_MLC.m @@ -80,7 +80,8 @@ matRad_visApertureInfo(resultGUI.apertureInfo) %% Plot the Resulting Dose Slice % Just let's plot the transversal iso-center dose slice -slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); figure, imagesc(resultGUI.physicalDose(:,:,slice)),colorbar, colormap(jet) diff --git a/examples/matRad_example1_phantom.m b/examples/matRad_example1_phantom.m index 374c92330..74f961dc6 100644 --- a/examples/matRad_example1_phantom.m +++ b/examples/matRad_example1_phantom.m @@ -135,7 +135,8 @@ matRadGUI %% Plot the resulting dose slice plane = 3; -slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); doseWindow = [0 max([resultGUI.physicalDose(:)])]; figure,title('phantom plan') diff --git a/examples/matRad_example2_photons.m b/examples/matRad_example2_photons.m index 87746338a..4ca61c384 100644 --- a/examples/matRad_example2_photons.m +++ b/examples/matRad_example2_photons.m @@ -173,7 +173,8 @@ %% Plot the Resulting Dose Slice % Let's plot the transversal iso-center dose slice -slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); figure imagesc(resultGUI.physicalDose(:,:,slice)),colorbar, colormap(jet); diff --git a/examples/matRad_example4_photonsMC.m b/examples/matRad_example4_photonsMC.m index 749dd78af..7c800cd59 100644 --- a/examples/matRad_example4_photonsMC.m +++ b/examples/matRad_example4_photonsMC.m @@ -72,7 +72,8 @@ %% Plot the Resulting Dose Slice % Just let's plot the transversal iso-center dose slice -slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); figure, imagesc(resultGUI.physicalDose(:,:,slice)),colorbar, colormap(jet) diff --git a/examples/matRad_example5_protons.m b/examples/matRad_example5_protons.m index 7f33be717..236610422 100644 --- a/examples/matRad_example5_protons.m +++ b/examples/matRad_example5_protons.m @@ -99,7 +99,8 @@ %% Plot the Resulting Dose Slice % Let's plot the transversal iso-center dose slice -slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); figure imagesc(resultGUI.RBExD(:,:,slice)),colorbar,colormap(jet) diff --git a/examples/matRad_example6_protonsNoise.m b/examples/matRad_example6_protonsNoise.m index 8c880c3c5..9ef7ceabc 100644 --- a/examples/matRad_example6_protonsNoise.m +++ b/examples/matRad_example6_protonsNoise.m @@ -115,7 +115,8 @@ %% Plot the Resulting Dose Slice % Let's plot the transversal iso-center dose slice -slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); figure imagesc(resultGUI.RBExD(:,:,slice)),colorbar, colormap(jet) diff --git a/examples/matRad_example7_carbon.m b/examples/matRad_example7_carbon.m index 0b3756881..417dd15a5 100644 --- a/examples/matRad_example7_carbon.m +++ b/examples/matRad_example7_carbon.m @@ -110,13 +110,15 @@ %% Plot the Resulting Dose Slice % Let's plot the transversal iso-center dose slice -slice = round(pln.propStf.isoCenter(3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); figure, imagesc(resultGUI.RBExD(:,:,slice)),colorbar, colormap(jet); %% Let's check out the LET % Let's plot the transversal iso-center LET slice -slice = round(pln.propStf.isoCenter(3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); figure; imagesc(resultGUI.LET(:,:,slice)),colorbar, colormap(jet); diff --git a/examples/matRad_example8_protonsRobust.m b/examples/matRad_example8_protonsRobust.m index 8a7a6de0f..6395fce54 100644 --- a/examples/matRad_example8_protonsRobust.m +++ b/examples/matRad_example8_protonsRobust.m @@ -134,7 +134,8 @@ %% Visualize results plane = 3; -slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); figure,matRad_plotSliceWrapper(gca,ct,cst,1,resultGUI.RBExD_beam1 ,plane,slice,[],[],colorcube,[],[0 max(resultGUI.RBExD_beam1(:))],[]);title('conventional plan - beam1') figure,matRad_plotSliceWrapper(gca,ct,cst,1,resultGUIrobust.RBExD_beam1,plane,slice,[],[],colorcube,[],[0 max(resultGUIrobust.RBExD_beam1(:))],[]);title('robust plan - beam1') diff --git a/examples/matRad_example9_4DDoseCalcMinimal.m b/examples/matRad_example9_4DDoseCalcMinimal.m index 063bcdfff..827568509 100644 --- a/examples/matRad_example9_4DDoseCalcMinimal.m +++ b/examples/matRad_example9_4DDoseCalcMinimal.m @@ -85,7 +85,8 @@ [resultGUI, timeSequence] = matRad_calc4dDose(ct, pln, dij, stf, cst, resultGUI); % plot the result in comparison to the static dose -slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z); +slice = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +slice = slice(3); figure diff --git a/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomCt.m b/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomCt.m index 03006e436..b8843f0c9 100644 --- a/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomCt.m +++ b/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomCt.m @@ -71,10 +71,11 @@ if ~any(isfield(ct,{'x','y','z'})) %positionOffset = transpose(ct.cubeDim ./ 2); %positionOffset = ct.cubeDim .* [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; - positionOffset = [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; - ct.x = ct.resolution.x*[0:ct.cubeDim(2)-1] - positionOffset(2); - ct.y = ct.resolution.y*[0:ct.cubeDim(1)-1] - positionOffset(1); - ct.z = ct.resolution.z*[0:ct.cubeDim(3)-1] - positionOffset(3); +% positionOffset = [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; +% ct.x = ct.resolution.x*[0:ct.cubeDim(2)-1] - positionOffset(2); +% ct.y = ct.resolution.y*[0:ct.cubeDim(1)-1] - positionOffset(1); +% ct.z = ct.resolution.z*[0:ct.cubeDim(3)-1] - positionOffset(3); + ct = matRad_getWorldAxes(ct); end obj.ct = ct; diff --git a/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTDoses.m b/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTDoses.m index 6d0af05af..1978a64e8 100644 --- a/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTDoses.m +++ b/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTDoses.m @@ -34,10 +34,11 @@ if ~any(isfield(ct,{'x','y','z'})) %positionOffset = transpose(ct.cubeDim ./ 2); %positionOffset = ct.cubeDim .* [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; - positionOffset = [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; - ct.x = ct.resolution.x*[0:ct.cubeDim(2)-1] - positionOffset(2); - ct.y = ct.resolution.y*[0:ct.cubeDim(1)-1] - positionOffset(1); - ct.z = ct.resolution.z*[0:ct.cubeDim(3)-1] - positionOffset(3); +% positionOffset = [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; +% ct.x = ct.resolution.x*[0:ct.cubeDim(2)-1] - positionOffset(2); +% ct.y = ct.resolution.y*[0:ct.cubeDim(1)-1] - positionOffset(1); +% ct.z = ct.resolution.z*[0:ct.cubeDim(3)-1] - positionOffset(3); + ct = matRad_getWorldAxes(ct); end %% Meta data diff --git a/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTStruct.m b/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTStruct.m index 5670a34f5..74214c828 100644 --- a/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTStruct.m +++ b/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTStruct.m @@ -94,10 +94,11 @@ if ~any(isfield(ct,{'x','y','z'})) %positionOffset = transpose(ct.cubeDim ./ 2); %positionOffset = ct.cubeDim .* [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; - positionOffset = [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; - ct.x = ct.resolution.x*[0:ct.cubeDim(2)-1] - positionOffset(2); - ct.y = ct.resolution.y*[0:ct.cubeDim(1)-1] - positionOffset(1); - ct.z = ct.resolution.z*[0:ct.cubeDim(3)-1] - positionOffset(3); +% positionOffset = [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; +% ct.x = ct.resolution.x*[0:ct.cubeDim(2)-1] - positionOffset(2); +% ct.y = ct.resolution.y*[0:ct.cubeDim(1)-1] - positionOffset(1); +% ct.z = ct.resolution.z*[0:ct.cubeDim(3)-1] - positionOffset(3); + ct = matRad_getWorldAxes(ct); end diff --git a/matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m b/matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m index e0be40fce..51eaee04e 100644 --- a/matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m +++ b/matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m @@ -79,14 +79,11 @@ % to guarantee downwards compatibility with data that does not have % ct.x/y/z if ~any(isfield(ct,{'x','y','z'})) - dij.ctGrid.x = ct.resolution.x*[0:ct.cubeDim(2)-1]-ct.resolution.x/2; - dij.ctGrid.y = ct.resolution.y*[0:ct.cubeDim(1)-1]-ct.resolution.y/2; - dij.ctGrid.z = ct.resolution.z*[0:ct.cubeDim(3)-1]-ct.resolution.z/2; -else - dij.ctGrid.x = ct.x; - dij.ctGrid.y = ct.y; - dij.ctGrid.z = ct.z; + ct = matRad_getWorldAxes(ct); end +dij.ctGrid.x = ct.x; +dij.ctGrid.y = ct.y; +dij.ctGrid.z = ct.z; dij.ctGrid.dimensions = [numel(dij.ctGrid.y) numel(dij.ctGrid.x) numel(dij.ctGrid.z)]; dij.ctGrid.numOfVoxels = prod(dij.ctGrid.dimensions); diff --git a/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m b/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m index e030a97f2..77e60ecad 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m @@ -250,6 +250,11 @@ function setDefaults(this) this.config.LET_Sparse_Output = ~this.calcDoseDirect; end + %Create X Y Z vectors if not present + if ~any(isfield(ct,{'x','y','z'})) + ct = matRad_getWorldAxes(ct); + end + for scenarioIx = 1:this.multScen.totNumScen %For direct dose calculation totalWeights = 0; @@ -268,7 +273,7 @@ function setDefaults(this) for i = 1:length(stf) %Create new stf for MCsquare with energy layer ordering and %shifted scenario isocenter - stfMCsquare(i).isoCenter = stf(i).isoCenter + isoCenterShift; + stfMCsquare(i).isoCenter = stf(i).isoCenter - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z] + isoCenterShift; stfMCsquare(i).gantryAngle = mod(180-stf(i).gantryAngle,360); %Different MCsquare geometry stfMCsquare(i).couchAngle = stf(i).couchAngle; stfMCsquare(i).energies = unique([stf(i).ray.energy]); diff --git a/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m b/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m index d1d4bdf1f..794f41320 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m +++ b/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m @@ -85,6 +85,11 @@ function setDefaults(this) % initialize dij = this.initDoseCalc(ct,cst,stf); + %Create X Y Z vectors if not present + if ~any(isfield(ct,{'x','y','z'})) + ct = matRad_getWorldAxes(ct); + end + for shiftScen = 1:this.multScen.totNumShiftScen %Find first instance of the shift to select the shift values @@ -93,7 +98,7 @@ function setDefaults(this) scenStf = stf; % manipulate isocenter for k = 1:numel(scenStf) - scenStf(k).isoCenter = scenStf(k).isoCenter + this.multScen.isoShift(ixShiftScen,:); + scenStf(k).isoCenter = scenStf(k).isoCenter - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z] + this.multScen.isoShift(ixShiftScen,:); end if this.multScen.totNumShiftScen > 1 diff --git a/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m b/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m index d821f125e..60f34fce1 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m @@ -127,6 +127,11 @@ %Now we have to calibrate to the the beamlet width. calibrationFactor = this.absCalibrationFactor * (bixelWidth/50)^2; + + %Create X Y Z vectors if not present + if ~any(isfield(ct,{'x','y','z'})) + ct = matRad_getWorldAxes(ct); + end scenCount = 0; %run over all scenarios @@ -139,7 +144,8 @@ scenCount = scenCount + 1; % manipulate isocenter - shiftedIsoCenter = vertcat(stf(:).isoCenter) + this.multScen.isoShift(scenarioIx,:) + dij.doseGrid.isoCenterOffset; + shiftedIsoCenter = vertcat(stf(:).isoCenter) - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z]... + + this.multScen.isoShift(scenarioIx,:) + dij.doseGrid.isoCenterOffset; this.ompMCgeo.isoCenter = shiftedIsoCenter; tmpStf = stf; @@ -394,7 +400,8 @@ function getOmpMCsource(obj,stf) for i = 1:numOfBeams % loop over all beams % define beam source in physical coordinate system in cm - beamSource(i,:) = (stf(i).sourcePoint + stf(i).isoCenter)/10; + isoCenterBeam = stf(i).isoCenter - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z] + dij.doseGrid.isoCenterOffset; % coord axes start at 0 in mm + beamSource(i,:) = (stf(i).sourcePoint + isoCenterBeam)/10; for j = 1:stf(i).numOfRays % loop over all rays / for photons we only have one bixel per ray! @@ -410,18 +417,18 @@ function getOmpMCsource(obj,stf) % get bixel corner and delimiting vectors. % a) change coordinate system (Isocenter cs-> physical cs) and units mm -> cm - currCorner = (beamletCorners(1,:) + stf(i).isoCenter) ./ obj.scale; + currCorner = (beamletCorners(1,:) + isoCenterBeam) ./ obj.scale; bixelCorner(counter,:) = currCorner; - bixelSide1(counter,:) = (beamletCorners(2,:) + stf(i).isoCenter) ./ obj.scale - currCorner; - bixelSide2(counter,:) = (beamletCorners(4,:) + stf(i).isoCenter) ./ obj.scale - currCorner; + bixelSide1(counter,:) = (beamletCorners(2,:) + isoCenterBeam) ./ obj.scale - currCorner; + bixelSide2(counter,:) = (beamletCorners(4,:) + isoCenterBeam) ./ obj.scale - currCorner; if obj.visBool for k = 1:4 - currCornerVis = (beamletCorners(k,:) + stf(i).isoCenter)/10; + currCornerVis = (beamletCorners(k,:) + isoCenterBeam)/10; % rays connecting source and ray corner plot3([beamSource(i,1) currCornerVis(1)],[beamSource(i,2) currCornerVis(2)],[beamSource(i,3) currCornerVis(3)],'b') % connection between corners - lRayCorner = (beamletCorners(mod(k,4) + 1,:) + stf(i).isoCenter)/10; + lRayCorner = (beamletCorners(mod(k,4) + 1,:) + isoCenterBeam)/10; plot3([lRayCorner(1) currCornerVis(1)],[lRayCorner(2) currCornerVis(2)],[lRayCorner(3) currCornerVis(3)],'r') end end diff --git a/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m b/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m index 40b53edfe..9355f9167 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m @@ -489,7 +489,7 @@ function writeAllFiles(obj,ct,cst,stf,machine,w) % manipulate isocenter for k = 1:numel(stf) - stf(k).isoCenter = stf(k).isoCenter + this.multScen.isoShift(ixShiftScen,:); + stf(k).isoCenter = stf(k).isoCenter - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z] + this.multScen.isoShift(ixShiftScen,:); end % Delete previous topas files so there is no mix-up diff --git a/matRad/geometry/matRad_projectOnComponents.m b/matRad/geometry/matRad_projectOnComponents.m index 2ac55d688..f832e368c 100644 --- a/matRad/geometry/matRad_projectOnComponents.m +++ b/matRad/geometry/matRad_projectOnComponents.m @@ -31,7 +31,7 @@ % % Copyright 2017 the matRad development team. % -% This file is not part of the offical matRad release +% This file is not part of the offical matRad release ? % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/matRad/gui/widgets/matRad_ViewingWidget.m b/matRad/gui/widgets/matRad_ViewingWidget.m index 66840b3b4..55b5ab1b9 100644 --- a/matRad/gui/widgets/matRad_ViewingWidget.m +++ b/matRad/gui/widgets/matRad_ViewingWidget.m @@ -971,9 +971,9 @@ function UpdatePlot(this) cursorText{end+1,1} = ['Cube Index: ' mat2str(vCubeIdx)]; %Space Coordinates coords = zeros(1,3); - coords(1) = cubePos(2)*ct.resolution.y; - coords(2) = cubePos(1)*ct.resolution.x; - coords(3) = cubePos(3)*ct.resolution.z; + coords(1) = ct.y(cubePos(2)); + coords(2) = ct.x(cubePos(1)); + coords(3) = ct.z(cubePos(3)); cursorText{end+1,1} = ['Space Coordinates: ' mat2str(coords,5) ' mm']; ctVal = ct.cubeHU{1}(cubeIx(1),cubeIx(2),cubeIx(3)); @@ -1085,9 +1085,10 @@ function initValues(this) this.VOIPlotFlag(i) = true; end end - + if isfield(pln,'propStf') - planeCenters = ceil(pln.propStf.isoCenter(1,[2 1 3]) ./ [ct.resolution.x ct.resolution.y ct.resolution.z]); + isoCoordinates = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:), ct); + planeCenters = ceil(isoCoordinates); this.numOfBeams=pln.propStf.numOfBeams; else planeCenters = ceil(ct.cubeDim./ 2); @@ -1266,8 +1267,7 @@ function updateValues(this) % set isoCenter values % Note: only defined for the first Isocenter if isfield(pln,'propStf') - uniqueIsoCenters = unique(pln.propStf.isoCenter,'rows'); - this.vIsoCenter = round(uniqueIsoCenters(1,:)./[ct.resolution.x ct.resolution.y ct.resolution.z]); + this.vIsoCenter = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:), ct); else this.plotIsoCenter = false; end diff --git a/matRad/planAnalysis/matRad_compareDose.m b/matRad/planAnalysis/matRad_compareDose.m index 9dc87c27f..f83fb9a70 100644 --- a/matRad/planAnalysis/matRad_compareDose.m +++ b/matRad/planAnalysis/matRad_compareDose.m @@ -109,14 +109,15 @@ [~,s(3)] = max(sum(sum(cube1,1),2)); isoCenter = [ct.resolution.y*s(1) ct.resolution.x*s(2) ct.resolution.z*s(3)]; else - isoCenter = matRad_getIsoCenter(cst,ct,0); + + isoCenter = matRad_world2cubeCoords( matRad_getIsoCenter(cst,ct,0),ct); end resolution = [ct.resolution.x ct.resolution.y ct.resolution.z]; -slicename = {round(isoCenter(2)./resolution(2)),round(isoCenter(1)./resolution(1)),round(isoCenter(3)./resolution(3))}; +sliceName = {isoCenter(1),isoCenter(2),isoCenter(3)}; doseWindow = [0 max([cube1(:); cube2(:)])]; -planename = {'coronal','sagittal','axial'}; +planeName = {'coronal','sagittal','axial'}; %% Integral Energy Output intEnergy1 = matRad_calcIntEnergy(cube1,ct,pln); @@ -156,58 +157,58 @@ end for plane = 1:3 - matRad_cfg.dispInfo('Plotting %s plane...\n',planename{plane}); + matRad_cfg.dispInfo('Plotting %s plane...\n',planeName{plane}); % Initialize Figure - hfig.(planename{plane}).('fig') = figure('Renderer', 'painters', 'Position', [10 50 800 800]); + hfig.(planeName{plane}).('fig') = figure('Renderer', 'painters', 'Position', [10 50 800 800]); set(gcf,'Color',[1 1 1]); % Plot Dose 1 - hfig.(planename{plane}).('cube1').Axes = subplot(2,2,1); - [hfig.(planename{plane}).('cube1').CMap,... - hfig.(planename{plane}).('cube1').Dose,... - hfig.(planename{plane}).('cube1').Ct,... - hfig.(planename{plane}).('cube1').Contour,... - hfig.(planename{plane}).('cube1').IsoDose] = ... - matRad_plotSliceWrapper(gca,ct,cstHandle,1,cube1,plane,slicename{plane},[],[],colorcube,jet,doseWindow,[],100); + hfig.(planeName{plane}).('cube1').Axes = subplot(2,2,1); + [hfig.(planeName{plane}).('cube1').CMap,... + hfig.(planeName{plane}).('cube1').Dose,... + hfig.(planeName{plane}).('cube1').Ct,... + hfig.(planeName{plane}).('cube1').Contour,... + hfig.(planeName{plane}).('cube1').IsoDose] = ... + matRad_plotSliceWrapper(gca,ct,cstHandle,1,cube1,plane,sliceName{plane},[],[],colorcube,jet,doseWindow,[],100); % Plot Dose 2 - hfig.(planename{plane}).('cube2').Axes = subplot(2,2,2); - [hfig.(planename{plane}).('cube2').CMap,... - hfig.(planename{plane}).('cube2').Dose,... - hfig.(planename{plane}).('cube2').Ct,... - hfig.(planename{plane}).('cube2').Contour,... - hfig.(planename{plane}).('cube2').IsoDose] = ... - matRad_plotSliceWrapper(gca,ct,cstHandle,1,cube2,plane,slicename{plane},[],[],colorcube,jet,doseWindow,[],100); + hfig.(planeName{plane}).('cube2').Axes = subplot(2,2,2); + [hfig.(planeName{plane}).('cube2').CMap,... + hfig.(planeName{plane}).('cube2').Dose,... + hfig.(planeName{plane}).('cube2').Ct,... + hfig.(planeName{plane}).('cube2').Contour,... + hfig.(planeName{plane}).('cube2').IsoDose] = ... + matRad_plotSliceWrapper(gca,ct,cstHandle,1,cube2,plane,sliceName{plane},[],[],colorcube,jet,doseWindow,[],100); % Plot absolute difference - hfig.(planename{plane}).('diff').Axes = subplot(2,2,3); - [hfig.(planename{plane}).('diff').CMap,... - hfig.(planename{plane}).('diff').Dose,... - hfig.(planename{plane}).('diff').Ct,... - hfig.(planename{plane}).('diff').Contour,... - hfig.(planename{plane}).('diff').IsoDose] = ... - matRad_plotSliceWrapper(gca,ct,cstHandle,1,differenceCube,plane,slicename{plane},[],[],colorcube,diffCMap,doseDiffWindow,[],100); + hfig.(planeName{plane}).('diff').Axes = subplot(2,2,3); + [hfig.(planeName{plane}).('diff').CMap,... + hfig.(planeName{plane}).('diff').Dose,... + hfig.(planeName{plane}).('diff').Ct,... + hfig.(planeName{plane}).('diff').Contour,... + hfig.(planeName{plane}).('diff').IsoDose] = ... + matRad_plotSliceWrapper(gca,ct,cstHandle,1,differenceCube,plane,sliceName{plane},[],[],colorcube,diffCMap,doseDiffWindow,[],100); % Plot gamma analysis - hfig.(planename{plane}).('gamma').Axes = subplot(2,2,4); + hfig.(planeName{plane}).('gamma').Axes = subplot(2,2,4); gammaCMap = matRad_getColormap('gammaIndex'); - [hfig.(planename{plane}).('gamma').CMap,... - hfig.(planename{plane}).('gamma').Dose,... - hfig.(planename{plane}).('gamma').Ct,... - hfig.(planename{plane}).('gamma').Contour,... - hfig.(planename{plane}).('gamma').IsoDose]=... - matRad_plotSliceWrapper(gca,ct,cstHandle,1,gammaCube,plane,slicename{plane},[],[],colorcube,gammaCMap,doseGammaWindow,[],100); + [hfig.(planeName{plane}).('gamma').CMap,... + hfig.(planeName{plane}).('gamma').Dose,... + hfig.(planeName{plane}).('gamma').Ct,... + hfig.(planeName{plane}).('gamma').Contour,... + hfig.(planeName{plane}).('gamma').IsoDose]=... + matRad_plotSliceWrapper(gca,ct,cstHandle,1,gammaCube,plane,sliceName{plane},[],[],colorcube,gammaCMap,doseGammaWindow,[],100); % Adjusting axes - matRad_plotAxisLabels(hfig.(planename{plane}).('cube1').Axes,ct,plane,slicename{plane},[],100); - set(get(hfig.(planename{plane}).('cube1').Axes, 'title'), 'string', 'Dose 1'); - matRad_plotAxisLabels(hfig.(planename{plane}).('cube2').Axes,ct,plane,slicename{plane},[],100); - set(get(hfig.(planename{plane}).('cube2').Axes, 'title'), 'string', 'Dose 2'); - matRad_plotAxisLabels(hfig.(planename{plane}).('diff').Axes,ct,plane,slicename{plane},[],100); - set(get(hfig.(planename{plane}).('diff').Axes, 'title'), 'string', 'Absolute difference'); - matRad_plotAxisLabels(hfig.(planename{plane}).('gamma').Axes,ct,plane,slicename{plane},[],100); - set(get(hfig.(planename{plane}).('gamma').Axes, 'title'), 'string', {[num2str(gammaPassRate{1,2},5) '% of points > ' num2str(relDoseThreshold) '% pass gamma criterion (' num2str(relDoseThreshold) '% / ' num2str(dist2AgreeMm) 'mm)']; ['with ' num2str(2^n-1) ' interpolation points']}); + matRad_plotAxisLabels(hfig.(planeName{plane}).('cube1').Axes,ct,plane,sliceName{plane},[],100); + set(get(hfig.(planeName{plane}).('cube1').Axes, 'title'), 'string', 'Dose 1'); + matRad_plotAxisLabels(hfig.(planeName{plane}).('cube2').Axes,ct,plane,sliceName{plane},[],100); + set(get(hfig.(planeName{plane}).('cube2').Axes, 'title'), 'string', 'Dose 2'); + matRad_plotAxisLabels(hfig.(planeName{plane}).('diff').Axes,ct,plane,sliceName{plane},[],100); + set(get(hfig.(planeName{plane}).('diff').Axes, 'title'), 'string', 'Absolute difference'); + matRad_plotAxisLabels(hfig.(planeName{plane}).('gamma').Axes,ct,plane,sliceName{plane},[],100); + set(get(hfig.(planeName{plane}).('gamma').Axes, 'title'), 'string', {[num2str(gammaPassRate{1,2},5) '% of points > ' num2str(relDoseThreshold) '% pass gamma criterion (' num2str(relDoseThreshold) '% / ' num2str(dist2AgreeMm) 'mm)']; ['with ' num2str(2^n-1) ' interpolation points']}); end end @@ -217,13 +218,13 @@ if enable(2) == 1 matRad_cfg.dispInfo('Plotting profiles...\n'); fontsize = 12; - profilex{1} = squeeze(cube1(slicename{1},:,slicename{3})); - profiley{1} = squeeze(cube1(:,slicename{2},slicename{3})); - profilez{1} = squeeze(cube1(slicename{1},slicename{2},:)); + profilex{1} = squeeze(cube1(sliceName{1},:,sliceName{3})); + profiley{1} = squeeze(cube1(:,sliceName{2},sliceName{3})); + profilez{1} = squeeze(cube1(sliceName{1},sliceName{2},:)); - profilex{2} = squeeze(cube2(slicename{1},:,slicename{3})); - profiley{2} = squeeze(cube2(:,slicename{2},slicename{3})); - profilez{2} = squeeze(cube2(slicename{1},slicename{2},:)); + profilex{2} = squeeze(cube2(sliceName{1},:,sliceName{3})); + profiley{2} = squeeze(cube2(:,sliceName{2},sliceName{3})); + profilez{2} = squeeze(cube2(sliceName{1},sliceName{2},:)); posX = resolution(1)*(1:length(profilex{1})); posY = resolution(2)*(1:length(profiley{1})); @@ -277,7 +278,7 @@ legend({'Dose 1','Dose 2'},'Location','southeast') legend boxoff - set(hfig.profiles.fig,'name',['Profiles:, x=',num2str(slicename{1}),'mm, y=',num2str(slicename{2}),'mm, z=',num2str(slicename{3}),'mm']); + set(hfig.profiles.fig,'name',['Profiles:, x=',num2str(sliceName{1}),'mm, y=',num2str(sliceName{2}),'mm, z=',num2str(sliceName{3}),'mm']); end diff --git a/matRad/plotting/matRad_plotAxisLabels.m b/matRad/plotting/matRad_plotAxisLabels.m index 2207ca0ec..e3215c2d7 100644 --- a/matRad/plotting/matRad_plotAxisLabels.m +++ b/matRad/plotting/matRad_plotAxisLabels.m @@ -37,19 +37,25 @@ function matRad_plotAxisLabels(axesHandle,ct,plane,slice,defaultFontSize,tickdi defaultFontSize = 12; end +if ~isfield(ct,'x') || isempty(ct.x) + ct = matRad_getWorldAxes(ct); +end + if ~exist('tickdist','var') || isempty(tickdist) - tickdist = 50; + tickdist = abs(ct.x(end)-ct.x(1))/10; end + %% Set axis labels and plot iso center if plane == 3% Axial plane if ~isempty(ct.resolution.x) && ~isempty(ct.resolution.y) - set(axesHandle,'XTick',0:tickdist/ct.resolution.x:1000); - set(axesHandle,'YTick',0:tickdist/ct.resolution.y:1000); - set(axesHandle,'XTickLabel',0:tickdist:1000*ct.resolution.x); - set(axesHandle,'YTickLabel',0:tickdist:1000*ct.resolution.y); + set(axesHandle,'XTick',linspace(0,ct.x(end)-ct.x(1),10)./ct.resolution.x); + set(axesHandle,'YTick',linspace(0,ct.y(end)-ct.y(1),10)./ct.resolution.y); + set(axesHandle,'XTickLabel',round(linspace(ct.x(1),ct.x(end),10))); + set(axesHandle,'YTickLabel',round(linspace(ct.y(1),ct.y(end),10))); xlabel(axesHandle,'x [mm]','FontSize',defaultFontSize) ylabel(axesHandle,'y [mm]','FontSize',defaultFontSize) - title(axesHandle,['axial plane z = ' num2str(ct.resolution.z*slice) ' [mm]'],'FontSize',defaultFontSize) + vcoord = matRad_cube2worldCoords([0,0,slice],ct); + title(axesHandle,['axial plane z = ' num2str(vcoord(3)) ' [mm]'],'FontSize',defaultFontSize); else xlabel(axesHandle,'x [voxels]','FontSize',defaultFontSize) ylabel(axesHandle,'y [voxels]','FontSize',defaultFontSize) @@ -57,13 +63,14 @@ function matRad_plotAxisLabels(axesHandle,ct,plane,slice,defaultFontSize,tickdi end elseif plane == 2 % Sagittal plane if ~isempty(ct.resolution.y) && ~isempty(ct.resolution.z) - set(axesHandle,'XTick',0:tickdist/ct.resolution.z:1000) - set(axesHandle,'YTick',0:tickdist/ct.resolution.y:1000) - set(axesHandle,'XTickLabel',0:tickdist:1000*ct.resolution.z) - set(axesHandle,'YTickLabel',0:tickdist:1000*ct.resolution.y) + set(axesHandle,'XTick',linspace(0,ct.z(end)-ct.z(1),10)./ct.resolution.z); + set(axesHandle,'YTick',linspace(0,ct.y(end)-ct.y(1),10)./ct.resolution.y); + set(axesHandle,'XTickLabel',round(linspace(ct.z(1),ct.z(end),10))); + set(axesHandle,'YTickLabel',round(linspace(ct.y(1),ct.y(end),10))); xlabel(axesHandle,'z [mm]','FontSize',defaultFontSize); ylabel(axesHandle,'y [mm]','FontSize',defaultFontSize); - title(axesHandle,['sagittal plane x = ' num2str(ct.resolution.x*slice) ' [mm]'],'FontSize',defaultFontSize) + vcoord = matRad_cube2worldCoords([slice,0,0],ct); + title(axesHandle,['sagittal plane x = ' num2str(vcoord(1)) ' [mm]'],'FontSize',defaultFontSize); else xlabel(axesHandle,'z [voxels]','FontSize',defaultFontSize) ylabel(axesHandle,'y [voxels]','FontSize',defaultFontSize) @@ -71,13 +78,14 @@ function matRad_plotAxisLabels(axesHandle,ct,plane,slice,defaultFontSize,tickdi end elseif plane == 1 % Coronal plane if ~isempty(ct.resolution.x) && ~isempty(ct.resolution.z) - set(axesHandle,'XTick',0:tickdist/ct.resolution.z:1000) - set(axesHandle,'YTick',0:tickdist/ct.resolution.x:1000) - set(axesHandle,'XTickLabel',0:tickdist:1000*ct.resolution.z) - set(axesHandle,'YTickLabel',0:tickdist:1000*ct.resolution.x) + set(axesHandle,'XTick',linspace(0,ct.z(end)-ct.z(1),10)./ct.resolution.z); + set(axesHandle,'YTick',linspace(0,ct.x(end)-ct.x(1),10)./ct.resolution.x); + set(axesHandle,'XTickLabel',round(linspace(ct.z(1),ct.z(end),10))); + set(axesHandle,'YTickLabel',round(linspace(ct.x(1),ct.x(end),10))); xlabel(axesHandle,'z [mm]','FontSize',defaultFontSize) ylabel(axesHandle,'x [mm]','FontSize',defaultFontSize) - title(axesHandle,['coronal plane y = ' num2str(ct.resolution.y*slice) ' [mm]'],'FontSize',defaultFontSize) + vcoord = matRad_cube2worldCoords([0,slice,0],ct); + title(axesHandle,['coronal plane y = ' num2str(vcoord(2)) ' [mm]'],'FontSize',defaultFontSize) else xlabel(axesHandle,'z [voxels]','FontSize',defaultFontSize) ylabel(axesHandle,'x [voxels]','FontSize',defaultFontSize) diff --git a/matRad/plotting/matRad_plotIsoCenterMarker.m b/matRad/plotting/matRad_plotIsoCenterMarker.m index d63c9eb7e..a96f6240e 100644 --- a/matRad/plotting/matRad_plotIsoCenterMarker.m +++ b/matRad/plotting/matRad_plotIsoCenterMarker.m @@ -40,7 +40,7 @@ for i = 1:size(uniqueIsoCenters,1) - vIsoCenter = round(uniqueIsoCenters(i,:)./[ct.resolution.x ct.resolution.y ct.resolution.z]); + vIsoCenter = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:), ct); if plane == 3% Axial plane vIsoCenterPlot = [vIsoCenter(1) vIsoCenter(2)]; if vIsoCenter(3) == slice @@ -50,7 +50,7 @@ end elseif plane == 2 vIsoCenterPlot = [vIsoCenter(3) vIsoCenter(2)]; - if vIsoCenter(2) == slice + if vIsoCenter(1) == slice isoCenterDirection = 0; else isoCenterDirection = sign(double(vIsoCenter(1) > slice) - 0.5); @@ -58,7 +58,7 @@ elseif plane == 1 vIsoCenterPlot = [vIsoCenter(3) vIsoCenter(1)]; - if vIsoCenter(1) == slice + if vIsoCenter(2) == slice isoCenterDirection = 0; else isoCenterDirection = sign(double(vIsoCenter(2) > slice) - 0.5); diff --git a/matRad/plotting/matRad_plotProjectedGantryAngles.m b/matRad/plotting/matRad_plotProjectedGantryAngles.m index a543f1e66..95eb2a070 100644 --- a/matRad/plotting/matRad_plotProjectedGantryAngles.m +++ b/matRad/plotting/matRad_plotProjectedGantryAngles.m @@ -35,15 +35,14 @@ function matRad_plotProjectedGantryAngles(axesHandle,pln,ct,plane) meanIsoCenter = mean(pln.propStf.isoCenter,1); - xOffset = meanIsoCenter(1)/ct.resolution.x; - yOffset = meanIsoCenter(2)/ct.resolution.y; + cubeIso = matRad_world2cubeCoords(meanIsoCenter,ct); % find radius of inner circle from isocenter - r = 0.8*min([abs([1 ct.cubeDim(1)]-xOffset) abs([1 ct.cubeDim(2)]-yOffset)]); + r = 0.8*min([abs([1 ct.cubeDim(1)]-cubeIso(1)) abs([1 ct.cubeDim(2)]-cubeIso(2))]); % coordinates of circle - x = r*cosd(0:360)+xOffset; - y = r*sind(0:360)+yOffset; + x = r*cosd(0:360)+cubeIso(1); + y = r*sind(0:360)+cubeIso(2); gantryAngleVisColor = 'w'; @@ -52,17 +51,17 @@ function matRad_plotProjectedGantryAngles(axesHandle,pln,ct,plane) % add text txt = '180°'; - text(axesHandle,1.1*r*sind(0)+xOffset,1.1*r*cosd(0)+yOffset,txt,'Color',gantryAngleVisColor) + text(axesHandle,1.1*r*sind(0)+cubeIso(1),1.1*r*cosd(0)+cubeIso(2),txt,'Color',gantryAngleVisColor) txt = '90°'; - text(axesHandle,1.1*r*sind(90)+xOffset,1.1*r*cosd(90)+yOffset,txt,'Color',gantryAngleVisColor) + text(axesHandle,1.1*r*sind(90)+cubeIso(1),1.1*r*cosd(90)+cubeIso(2),txt,'Color',gantryAngleVisColor) txt = '0°'; - text(axesHandle,1.1*r*sind(180)+xOffset,1.1*r*cosd(180)+yOffset,txt,'Color',gantryAngleVisColor) + text(axesHandle,1.1*r*sind(180)+cubeIso(1),1.1*r*cosd(180)+cubeIso(2),txt,'Color',gantryAngleVisColor) txt = '270°'; - text(axesHandle,1.22*r*sind(270)+xOffset,1.22*r*cosd(270)+yOffset,txt,'Color',gantryAngleVisColor) + text(axesHandle,1.22*r*sind(270)+cubeIso(1),1.22*r*cosd(270)+cubeIso(2),txt,'Color',gantryAngleVisColor) % plot gantry angles for i = 1:numel(pln.propStf.gantryAngles) - plot(axesHandle,[0 r*sind(180-pln.propStf.gantryAngles(i))]+xOffset,[0 r*cosd(180-pln.propStf.gantryAngles(i))]+yOffset,'LineWidth',1,'Color',gantryAngleVisColor) + plot(axesHandle,[0 r*sind(180-pln.propStf.gantryAngles(i))]+cubeIso(1),[0 r*cosd(180-pln.propStf.gantryAngles(i))]+cubeIso(2),'LineWidth',1,'Color',gantryAngleVisColor) end end \ No newline at end of file From f55f43b3569a543a952d9c3ff83f0068fda9ace8 Mon Sep 17 00:00:00 2001 From: amitantony Date: Mon, 22 Jul 2024 12:31:54 +0200 Subject: [PATCH 03/10] example 5 bug fix --- examples/matRad_example5_protons.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/matRad_example5_protons.m b/examples/matRad_example5_protons.m index 236610422..8aa22f021 100644 --- a/examples/matRad_example5_protons.m +++ b/examples/matRad_example5_protons.m @@ -141,8 +141,8 @@ matRad_plotSliceWrapper(gca,ct,cst,1,absDiffCube,plane,slice,[],[],colorcube); % Let's plot single profiles that are perpendicular to the beam direction -ixProfileY = round(pln.propStf.isoCenter(1,2)./ct.resolution.y); - +ixProfileY = matRad_world2cubeCoords(pln.propStf.isoCenter(1,:),ct); +ixProfileY = ixProfileY(2) profileOrginal = resultGUI.RBExD(:,ixProfileY,slice); profileShifted = resultGUI_isoShift.RBExD(:,ixProfileY,slice); From 7a1d29cbef1a09e127e31cabbfa6054bc7506839 Mon Sep 17 00:00:00 2001 From: amitantony Date: Mon, 22 Jul 2024 16:39:19 +0200 Subject: [PATCH 04/10] photons OMP MC bug fix --- .../+DoseEngines/matRad_PhotonOmpMCEngine.m | 27 +++++++++++++------ 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m b/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m index 60f34fce1..a79a5a8a7 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m @@ -222,7 +222,19 @@ this.getOmpMCgeometry(dij.doseGrid); %% Create beamlet source - this.getOmpMCsource(stf); + if ~isfield(ct,'x') + ct = matRad_getWorldAxes(ct); + end + + tmpStf = stf; + for k = 1:length(stf) + shiftedIsoCenter = vertcat(stf(:).isoCenter) - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z]... + + dij.doseGrid.isoCenterOffset; + + tmpStf(k).isoCenter = shiftedIsoCenter; + end + + this.getOmpMCsource(tmpStf); end end @@ -400,8 +412,7 @@ function getOmpMCsource(obj,stf) for i = 1:numOfBeams % loop over all beams % define beam source in physical coordinate system in cm - isoCenterBeam = stf(i).isoCenter - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z] + dij.doseGrid.isoCenterOffset; % coord axes start at 0 in mm - beamSource(i,:) = (stf(i).sourcePoint + isoCenterBeam)/10; + beamSource(i,:) = (stf(i).sourcePoint + stf(i).isoCenter)/10; for j = 1:stf(i).numOfRays % loop over all rays / for photons we only have one bixel per ray! @@ -417,18 +428,18 @@ function getOmpMCsource(obj,stf) % get bixel corner and delimiting vectors. % a) change coordinate system (Isocenter cs-> physical cs) and units mm -> cm - currCorner = (beamletCorners(1,:) + isoCenterBeam) ./ obj.scale; + currCorner = (beamletCorners(1,:) + stf(i).isoCenter) ./ obj.scale; bixelCorner(counter,:) = currCorner; - bixelSide1(counter,:) = (beamletCorners(2,:) + isoCenterBeam) ./ obj.scale - currCorner; - bixelSide2(counter,:) = (beamletCorners(4,:) + isoCenterBeam) ./ obj.scale - currCorner; + bixelSide1(counter,:) = (beamletCorners(2,:) + stf(i).isoCenter) ./ obj.scale - currCorner; + bixelSide2(counter,:) = (beamletCorners(4,:) + stf(i).isoCenter) ./ obj.scale - currCorner; if obj.visBool for k = 1:4 - currCornerVis = (beamletCorners(k,:) + isoCenterBeam)/10; + currCornerVis = (beamletCorners(k,:) + stf(i).isoCenter)/10; % rays connecting source and ray corner plot3([beamSource(i,1) currCornerVis(1)],[beamSource(i,2) currCornerVis(2)],[beamSource(i,3) currCornerVis(3)],'b') % connection between corners - lRayCorner = (beamletCorners(mod(k,4) + 1,:) + isoCenterBeam)/10; + lRayCorner = (beamletCorners(mod(k,4) + 1,:) + stf(i).isoCenter)/10; plot3([lRayCorner(1) currCornerVis(1)],[lRayCorner(2) currCornerVis(2)],[lRayCorner(3) currCornerVis(3)],'r') end end From 1524c7bc33897da6a6c17cfb9c4e1726ed9a7c08 Mon Sep 17 00:00:00 2001 From: amitantony Date: Mon, 22 Jul 2024 17:17:05 +0200 Subject: [PATCH 05/10] example 12 bugfix, matrad_generateSingleBixelStf. NOTE: results differ by small degree. is the example result reproducible everytime ? --- matRad/tools/matRad_generateSingleBixelStf.m | 24 +++++++++++--------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/matRad/tools/matRad_generateSingleBixelStf.m b/matRad/tools/matRad_generateSingleBixelStf.m index 71dd0a914..65a223bbf 100644 --- a/matRad/tools/matRad_generateSingleBixelStf.m +++ b/matRad/tools/matRad_generateSingleBixelStf.m @@ -95,18 +95,17 @@ %Get Isocenter voxel as target if ~any(isfield(ct,{'x','y','z'})) - ct.x = ct.resolution.x*[1:ct.cubeDim(2)]-ct.resolution.x/2; - ct.y = ct.resolution.y*[1:ct.cubeDim(1)]-ct.resolution.y/2; - ct.z = ct.resolution.z*[1:ct.cubeDim(3)]-ct.resolution.z/2; + ct = matRad_getWorldAxes(ct); + end %xVox = ct.x + ct.resolution.x/2; %yVox = ct.y + ct.resolution.y/2; %zVox = ct.z + ct.resolution.z/2; -xVox = ct.resolution.x*[1:ct.cubeDim(2)]-ct.resolution.x/2; -yVox = ct.resolution.y*[1:ct.cubeDim(1)]-ct.resolution.y/2; -zVox = ct.resolution.z*[1:ct.cubeDim(3)]-ct.resolution.z/2; +% xVox = ct.resolution.x*[1:ct.cubeDim(2)]-ct.resolution.x/2; +% yVox = ct.resolution.y*[1:ct.cubeDim(1)]-ct.resolution.y/2; +% zVox = ct.resolution.z*[1:ct.cubeDim(3)]-ct.resolution.z/2; % loop over all angles for i = 1:length(pln.propStf.gantryAngles) @@ -123,9 +122,9 @@ stf(i).totalNumOfBixels = 1; stf(i).machine = pln.machine; - x = floor(matRad_interp1(xVox,[1:ct.cubeDim(2)]',stf.isoCenter(1))); - y = floor(matRad_interp1(yVox,[1:ct.cubeDim(1)]',stf.isoCenter(2))); - z = floor(matRad_interp1(zVox,[1:ct.cubeDim(3)]',stf.isoCenter(3))); + x = floor(matRad_interp1(ct.x,[1:ct.cubeDim(2)]',stf.isoCenter(1))); + y = floor(matRad_interp1(ct.y,[1:ct.cubeDim(1)]',stf.isoCenter(2))); + z = floor(matRad_interp1(ct.z,[1:ct.cubeDim(3)]',stf.isoCenter(3))); %Voxel index of Isocenter isoIx = [y x z]; @@ -159,7 +158,10 @@ stf(i).ray.rayPos_bev = [0 0 0]; stf(i).ray.targetPoint_bev = [0 SAD 0]; - stf(i).ray.rayPos = stf.isoCenter; + + shiftedIsoCenter = vertcat(stf(:).isoCenter) - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z] ; + + stf(i).ray.rayPos = shiftedIsoCenter; stf(i).ray.targetPoint = [0 SAD 0] * rotMat_vectors_T; @@ -168,7 +170,7 @@ stf.longitudinalSpotSpacing = longitudinalSpotSpacing; % ray tracing necessary to determine depth of the target - [alphas,l,rho,d12,~] = matRad_siddonRayTracer(stf(i).isoCenter, ... + [alphas,l,rho,d12,~] = matRad_siddonRayTracer(shiftedIsoCenter, ... ct.resolution, ... stf(i).sourcePoint, ... stf(i).ray.targetPoint, ... From d13ffbf91214040ab1cc309ed3902ff69033ecdb Mon Sep 17 00:00:00 2001 From: amitantony Date: Tue, 23 Jul 2024 15:36:21 +0200 Subject: [PATCH 06/10] updates to the functions + tests --- matRad/geometry/matRad_cube2worldCoords.m | 18 +++++++---- matRad/geometry/matRad_getWorldAxes.m | 5 +-- matRad/geometry/matRad_world2cubeCoords.m | 38 +++++++++++++++++----- test/core/test_cube2worldCoords.m | 39 +++++++++++++++++++++++ test/core/test_getWorldAxes.m | 23 +++++++++++++ test/core/test_world2cubeCoords.m | 31 ++++++++++++++++++ 6 files changed, 138 insertions(+), 16 deletions(-) create mode 100644 test/core/test_cube2worldCoords.m create mode 100644 test/core/test_getWorldAxes.m create mode 100644 test/core/test_world2cubeCoords.m diff --git a/matRad/geometry/matRad_cube2worldCoords.m b/matRad/geometry/matRad_cube2worldCoords.m index b1cd6e787..370b3f168 100644 --- a/matRad/geometry/matRad_cube2worldCoords.m +++ b/matRad/geometry/matRad_cube2worldCoords.m @@ -6,7 +6,7 @@ % % ct: matRad ct struct % vCoord : Voxel Coordinates [vx vy vz] -% +% coord: worldCoordinates [x y z ] mm % References % - % @@ -23,20 +23,26 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% coord = []; % [x y z ] mm -if nargin == 2 + +if nargin == 2 && ~isempty(vCoord) if isfield(ct, 'dicomInfo') && isfield(ct.dicomInfo,'ImagePositionPatient') firstVox = ct.dicomInfo.ImagePositionPatient; else - firstVox = - (ct.cubeDim./2).*[ct.resolution.x ct.resolution.y ct.resolution.z] ; + firstVox = -(ct.cubeDim./2).*[ct.resolution.x ct.resolution.y ct.resolution.z] ; end - try + if prod(vCoord<=ct.cubeDim) && prod(vCoord>0) coord(1,:) = firstVox(1) + (vCoord(:,1) - 1 ) *ct.resolution.x; coord(2,:) = firstVox(2) + (vCoord(:,2) - 1 ) *ct.resolution.y; coord(3,:) = firstVox(3) + (vCoord(:,3) - 1 ) *ct.resolution.z; coord = coord'; - catch - + else + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Queried cube coordinate is outside the ct cube'); end + +else + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Cannot compute coordinates without matRad input coordinates or ct structure'); end end \ No newline at end of file diff --git a/matRad/geometry/matRad_getWorldAxes.m b/matRad/geometry/matRad_getWorldAxes.m index 7577f41a5..5ec606cf3 100644 --- a/matRad/geometry/matRad_getWorldAxes.m +++ b/matRad/geometry/matRad_getWorldAxes.m @@ -19,7 +19,7 @@ % LICENSE file. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -matRad_cfg = MatRad_Config.instance(); + if nargin>0 % @@ -34,7 +34,8 @@ ct.y = firstVox(2) + ct.resolution.y*[0:ct.cubeDim(1)-1] ; ct.z = firstVox(3) + ct.resolution.z*[0:ct.cubeDim(3)-1] ; else - matRad_cfg.dispWarning('Cannot compute world coordinates without matRad ct structure'); + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Cannot compute world coordinates without matRad ct structure'); end end \ No newline at end of file diff --git a/matRad/geometry/matRad_world2cubeCoords.m b/matRad/geometry/matRad_world2cubeCoords.m index a652fb949..c10fd5e06 100644 --- a/matRad/geometry/matRad_world2cubeCoords.m +++ b/matRad/geometry/matRad_world2cubeCoords.m @@ -4,6 +4,11 @@ % call % coord = world2cubeCoords(wCoord, ct) % +% +% wCoord: world coordinates (x,y,z)[mm] +% ct : matRad ct structure +% +% coord : cube index (x,y,z) % References % - % @@ -20,19 +25,36 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% coord = []; -if nargin == 2 - try + +if nargin == 2 && ~isempty(wCoord) + if ~ (isfield(ct,'x') && isfield(ct,'y') && isfield(ct,'z') ) ct = matRad_getWorldAxes(ct); end + % world bounds + boundx = [min(ct.x)-ct.resolution.x/2 max(ct.x)+ct.resolution.x/2]; + boundy = [min(ct.y)-ct.resolution.y/2 max(ct.y)+ct.resolution.y/2]; + boundz = [min(ct.z)-ct.resolution.z/2 max(ct.z)+ct.resolution.z/2]; - [~,coord(2)]=min(abs(ct.y-wCoord(2))); - [~,coord(1)]=min(abs(ct.x-wCoord(1))); - [~,coord(3)]=min(abs(ct.z-wCoord(3))); - - catch + if (wCoord(1)>=boundx(1) && wCoord(1)<=boundx(2)) && ... + (wCoord(2)>=boundy(1) && wCoord(2)<=boundy(2)) && ... + (wCoord(3)>=boundz(1) && wCoord(3)<=boundz(2)) + + [~,coord(1)]=min(abs(ct.x-wCoord(1))); + [~,coord(2)]=min(abs(ct.y-wCoord(2))); + [~,coord(3)]=min(abs(ct.z-wCoord(3))); + else + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Queried world coordinate is outside the ct'); + + end - end +% catch +% +% end +else + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Cannot compute cube coordinates without matRad input coordinates or ct structure'); end end diff --git a/test/core/test_cube2worldCoords.m b/test/core/test_cube2worldCoords.m new file mode 100644 index 000000000..451ab2b28 --- /dev/null +++ b/test/core/test_cube2worldCoords.m @@ -0,0 +1,39 @@ +function test_suite = test_cube2worldCoords + +test_functions = localfunctions(); + +initTestSuite; + +function test_empty_cube2worldCoords + ct = get_testCtHelper(); + assertExceptionThrown(@()matRad_cube2worldCoords()); + assertExceptionThrown(@()matRad_cube2worldCoords([],ct)); + assertExceptionThrown(@()matRad_cube2worldCoords([100 100 100],ct)); + + + +%with ct.dicom +function test_Dicom_cube2worldCoords + ct = get_testCtHelper(); + ct.dicomInfo.ImagePositionPatient = [-10 -10 -10]; + vCoord = [31 31 31] ; + expected = [ 80 80 80]; + assertEqual(matRad_cube2worldCoords(vCoord , ct), expected); + +% without ct.dicomInfo +function test_noDicom_cube2worldCoords + ct = get_testCtHelper(); + vCoord = [31 31 31] ; + expected = [ 0 0 0]; + assertEqual(matRad_cube2worldCoords(vCoord , ct), expected); + vCoord = [1 2 3]; + expected = [-90 -87 -84] ; + assertEqual(matRad_cube2worldCoords(vCoord , ct), expected); + + + +function ct = get_testCtHelper() + ct.resolution.x = 3; + ct.resolution.y = 3; + ct.resolution.z = 3; + ct.cubeDim = [60 60 60 ]; diff --git a/test/core/test_getWorldAxes.m b/test/core/test_getWorldAxes.m new file mode 100644 index 000000000..71c7edba3 --- /dev/null +++ b/test/core/test_getWorldAxes.m @@ -0,0 +1,23 @@ +function test_suite = test_getWorldAxes + +test_functions = localfunctions(); + +initTestSuite; + +function test_basic_getWorldAxes + assertExceptionThrown(@()matRad_getWorldAxes()); + + ct = get_testCtHelper(); + expectedX = [-5 -4 -3 -2 -1 0 1 2 3 4]; + expectedY = [-10 -8 -6 -4 -2 0 2 4 6 8]; + expectedZ = [-15 -12 -9 -6 -3 0 3 6 9 12]; + resultCt = matRad_getWorldAxes(ct); + assertEqual(resultCt.x,expectedX); + assertEqual(resultCt.y,expectedY); + assertEqual(resultCt.z,expectedZ); + +function ct = get_testCtHelper() + ct.resolution.x = 1; + ct.resolution.y = 2; + ct.resolution.z = 3; + ct.cubeDim = [10 10 10 ]; \ No newline at end of file diff --git a/test/core/test_world2cubeCoords.m b/test/core/test_world2cubeCoords.m new file mode 100644 index 000000000..455ac7124 --- /dev/null +++ b/test/core/test_world2cubeCoords.m @@ -0,0 +1,31 @@ +function test_suite = test_world2cubeCoords + +test_functions = localfunctions(); + +initTestSuite; + + +function test_empty_world2cubeCoords + ct = get_testCtHelper(); + assertExceptionThrown(@()matRad_world2cubeCoords()); + assertExceptionThrown(@()matRad_world2cubeCoords([],ct)); + + +function test_basic_world2cubeCoords + ct = get_testCtHelper(); + expected = [1 12 15]; + assertEqual(matRad_world2cubeCoords([-30 3 12], ct),expected); +% out of bounds + assertExceptionThrown(@()matRad_world2cubeCoords([-30,4,35],ct)); + + +function ct = get_testCtHelper() + ct.x = -30:3:27; + ct.y = -30:3:27; + ct.z = -30:3:27; + + ct.resolution.x = 3; + ct.resolution.y = 3; + ct.resolution.z = 3; + + ct.cubeDim = [20 20 20 ]; \ No newline at end of file From ba896c1c8eb0b510345dcac98d8e360cdda9649b Mon Sep 17 00:00:00 2001 From: amitantony Date: Fri, 26 Jul 2024 15:25:56 +0200 Subject: [PATCH 07/10] addressing rcomments from the PR review --- .../matRad_exportDicomRTDoses.m | 11 +--- .../matRad_exportDicomRTStruct.m | 11 +--- .../@matRad_DoseEngineBase/initDoseCalc.m | 5 +- .../matRad_ParticleMCsquareEngine.m | 5 +- .../matRad_PencilBeamEngineAbstract.m | 5 +- .../+DoseEngines/matRad_PhotonOmpMCEngine.m | 11 ++-- .../matRad_projectOnComponents.m | 7 ++- matRad/geometry/matRad_cube2worldCoords.m | 24 ++++--- matRad/geometry/matRad_getWorldAxes.m | 63 ++++++++++++------- matRad/geometry/matRad_world2cubeCoords.m | 28 +++++---- matRad/matRad_generateStf.m | 16 ++--- matRad/plotting/matRad_plotAxisLabels.m | 5 +- matRad/tools/matRad_generateSingleBixelStf.m | 12 +--- 13 files changed, 98 insertions(+), 105 deletions(-) rename matRad/{geometry => doseCalc}/matRad_projectOnComponents.m (90%) diff --git a/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTDoses.m b/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTDoses.m index 1978a64e8..78b89e731 100644 --- a/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTDoses.m +++ b/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTDoses.m @@ -31,15 +31,8 @@ %% CT check ct = obj.ct; -if ~any(isfield(ct,{'x','y','z'})) - %positionOffset = transpose(ct.cubeDim ./ 2); - %positionOffset = ct.cubeDim .* [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; -% positionOffset = [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; -% ct.x = ct.resolution.x*[0:ct.cubeDim(2)-1] - positionOffset(2); -% ct.y = ct.resolution.y*[0:ct.cubeDim(1)-1] - positionOffset(1); -% ct.z = ct.resolution.z*[0:ct.cubeDim(3)-1] - positionOffset(3); - ct = matRad_getWorldAxes(ct); -end +ct = matRad_getWorldAxes(ct); + %% Meta data storageClass = obj.rtDoseClassUID; diff --git a/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTStruct.m b/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTStruct.m index 74214c828..ccf7f554b 100644 --- a/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTStruct.m +++ b/matRad/dicom/@matRad_DicomExporter/matRad_exportDicomRTStruct.m @@ -91,15 +91,8 @@ ct = obj.ct; %Create X Y Z vectors if not present -if ~any(isfield(ct,{'x','y','z'})) - %positionOffset = transpose(ct.cubeDim ./ 2); - %positionOffset = ct.cubeDim .* [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; -% positionOffset = [ct.resolution.y, ct.resolution.x, ct.resolution.z] ./ 2; -% ct.x = ct.resolution.x*[0:ct.cubeDim(2)-1] - positionOffset(2); -% ct.y = ct.resolution.y*[0:ct.cubeDim(1)-1] - positionOffset(1); -% ct.z = ct.resolution.z*[0:ct.cubeDim(3)-1] - positionOffset(3); - ct = matRad_getWorldAxes(ct); -end +ct = matRad_getWorldAxes(ct); + %Since we are exporting HU directly --> no rescaling in any case diff --git a/matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m b/matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m index 51eaee04e..937559120 100644 --- a/matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m +++ b/matRad/doseCalc/+DoseEngines/@matRad_DoseEngineBase/initDoseCalc.m @@ -78,9 +78,8 @@ % to guarantee downwards compatibility with data that does not have % ct.x/y/z -if ~any(isfield(ct,{'x','y','z'})) - ct = matRad_getWorldAxes(ct); -end +ct = matRad_getWorldAxes(ct); + dij.ctGrid.x = ct.x; dij.ctGrid.y = ct.y; dij.ctGrid.z = ct.z; diff --git a/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m b/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m index 77e60ecad..db29ccef6 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m @@ -251,9 +251,8 @@ function setDefaults(this) end %Create X Y Z vectors if not present - if ~any(isfield(ct,{'x','y','z'})) - ct = matRad_getWorldAxes(ct); - end + ct = matRad_getWorldAxes(ct); + for scenarioIx = 1:this.multScen.totNumScen %For direct dose calculation diff --git a/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m b/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m index 794f41320..13f91fe95 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m +++ b/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m @@ -86,9 +86,8 @@ function setDefaults(this) dij = this.initDoseCalc(ct,cst,stf); %Create X Y Z vectors if not present - if ~any(isfield(ct,{'x','y','z'})) - ct = matRad_getWorldAxes(ct); - end + ct = matRad_getWorldAxes(ct); + for shiftScen = 1:this.multScen.totNumShiftScen diff --git a/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m b/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m index a79a5a8a7..87324fb27 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m @@ -129,9 +129,8 @@ calibrationFactor = this.absCalibrationFactor * (bixelWidth/50)^2; %Create X Y Z vectors if not present - if ~any(isfield(ct,{'x','y','z'})) - ct = matRad_getWorldAxes(ct); - end + ct = matRad_getWorldAxes(ct); + scenCount = 0; %run over all scenarios @@ -222,10 +221,8 @@ this.getOmpMCgeometry(dij.doseGrid); %% Create beamlet source - if ~isfield(ct,'x') - ct = matRad_getWorldAxes(ct); - end - + ct = matRad_getWorldAxes(ct); + tmpStf = stf; for k = 1:length(stf) shiftedIsoCenter = vertcat(stf(:).isoCenter) - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z]... diff --git a/matRad/geometry/matRad_projectOnComponents.m b/matRad/doseCalc/matRad_projectOnComponents.m similarity index 90% rename from matRad/geometry/matRad_projectOnComponents.m rename to matRad/doseCalc/matRad_projectOnComponents.m index f832e368c..ba3da5a23 100644 --- a/matRad/geometry/matRad_projectOnComponents.m +++ b/matRad/doseCalc/matRad_projectOnComponents.m @@ -31,7 +31,12 @@ % % Copyright 2017 the matRad development team. % -% This file is not part of the offical matRad release ? +% 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/LICENSES.txt. 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. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/matRad/geometry/matRad_cube2worldCoords.m b/matRad/geometry/matRad_cube2worldCoords.m index 370b3f168..4f755dc5a 100644 --- a/matRad/geometry/matRad_cube2worldCoords.m +++ b/matRad/geometry/matRad_cube2worldCoords.m @@ -1,10 +1,10 @@ -function coord = matRad_cube2worldCoords(vCoord, ct) +function coord = matRad_cube2worldCoords(vCoord, gridStruct) % matRad function to convert cube coordinates to world coordinates % % call -% coord = matRad_worldToCubeCoordinates(vCoord, ct) +% coord = matRad_worldToCubeCoordinates(vCoord, gridStruct) % -% ct: matRad ct struct +% gridStruct: matRad ct struct or dij.doseGrid/ctGrid struct % vCoord : Voxel Coordinates [vx vy vz] % coord: worldCoordinates [x y z ] mm % References @@ -25,15 +25,19 @@ coord = []; % [x y z ] mm if nargin == 2 && ~isempty(vCoord) - if isfield(ct, 'dicomInfo') && isfield(ct.dicomInfo,'ImagePositionPatient') - firstVox = ct.dicomInfo.ImagePositionPatient; + if isfield(gridStruct,'cubeDim') + gridStruct.dimensions = gridStruct.cubeDim; + end + + if isfield(gridStruct, 'dicomInfo') && isfield(gridStruct.dicomInfo,'ImagePositionPatient') + firstVox = gridStruct.dicomInfo.ImagePositionPatient; else - firstVox = -(ct.cubeDim./2).*[ct.resolution.x ct.resolution.y ct.resolution.z] ; + firstVox = -(gridStruct.cubeDim./2).*[gridStruct.resolution.x gridStruct.resolution.y gridStruct.resolution.z] ; end - if prod(vCoord<=ct.cubeDim) && prod(vCoord>0) - coord(1,:) = firstVox(1) + (vCoord(:,1) - 1 ) *ct.resolution.x; - coord(2,:) = firstVox(2) + (vCoord(:,2) - 1 ) *ct.resolution.y; - coord(3,:) = firstVox(3) + (vCoord(:,3) - 1 ) *ct.resolution.z; + if prod(prod(vCoord<=gridStruct.dimensions)) && prod(prod(vCoord>=0)) + coord(1,:) = firstVox(1) + (vCoord(:,1) - 1 ) *gridStruct.resolution.x; + coord(2,:) = firstVox(2) + (vCoord(:,2) - 1 ) *gridStruct.resolution.y; + coord(3,:) = firstVox(3) + (vCoord(:,3) - 1 ) *gridStruct.resolution.z; coord = coord'; else matRad_cfg = MatRad_Config.instance(); diff --git a/matRad/geometry/matRad_getWorldAxes.m b/matRad/geometry/matRad_getWorldAxes.m index 5ec606cf3..040e3e934 100644 --- a/matRad/geometry/matRad_getWorldAxes.m +++ b/matRad/geometry/matRad_getWorldAxes.m @@ -1,41 +1,56 @@ -function ct = matRad_getWorldAxes(ct) -% matRad function to compute and store world coordinates int ct.x -% +function gridStruct = matRad_getWorldAxes(gridStruct) +% matRad function to compute and store world coordinates into ct.x +% % call -% ct = matRad_getWorldAxes(ct) +% gridStruct = matRad_getWorldAxes(gridStruct) +% +% gridStruct: can be ct, dij.doseGrid,dij.ctGrid % % References % - % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Copyright 2015 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/LICENSES.txt. No part -% of the matRad project, including this file, may be copied, modified, -% propagated, or distributed except according to the terms contained in the +% Copyright 2015 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/LICENSES.txt. 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. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - +update = false; if nargin>0 - % - % check if dicominfo exists - if isfield(ct, 'dicomInfo') && isfield(ct.dicomInfo,'ImagePositionPatient') - firstVox = ct.dicomInfo.ImagePositionPatient; - else - firstVox = - (ct.cubeDim./2).*[ct.resolution.x ct.resolution.y ct.resolution.z] ; - end + if ~ (isfield(gridStruct,{'x','y','z'})) + update = true; + else + if isempty(gridStruct.x)||isempty(gridStruct.y)||isempty(gridStruct.z) + update = true; + end + end + + if update + % + if isfield(gridStruct,'cubeDim') + gridStruct.dimensions = gridStruct.cubeDim; + end + % check if dicominfo exists + if isfield(gridStruct, 'dicomInfo') && isfield(gridStruct.dicomInfo,'ImagePositionPatient') + firstVox = gridStruct.dicomInfo.ImagePositionPatient; + else + firstVox = - (gridStruct.dimensions./2).*[gridStruct.resolution.x gridStruct.resolution.y gridStruct.resolution.z] ; + end - ct.x = firstVox(1) + ct.resolution.x*[0:ct.cubeDim(2)-1] ; - ct.y = firstVox(2) + ct.resolution.y*[0:ct.cubeDim(1)-1] ; - ct.z = firstVox(3) + ct.resolution.z*[0:ct.cubeDim(3)-1] ; + gridStruct.x = firstVox(1) + gridStruct.resolution.x*[0:gridStruct.dimensions(2)-1] ; + gridStruct.y = firstVox(2) + gridStruct.resolution.y*[0:gridStruct.dimensions(1)-1] ; + gridStruct.z = firstVox(3) + gridStruct.resolution.z*[0:gridStruct.dimensions(3)-1] ; + end else matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Cannot compute world coordinates without matRad ct structure'); + matRad_cfg.dispError('Cannot compute world coordinates without matRad ct or grid structure'); end -end \ No newline at end of file +end diff --git a/matRad/geometry/matRad_world2cubeCoords.m b/matRad/geometry/matRad_world2cubeCoords.m index c10fd5e06..afa80ba16 100644 --- a/matRad/geometry/matRad_world2cubeCoords.m +++ b/matRad/geometry/matRad_world2cubeCoords.m @@ -1,4 +1,4 @@ -function coord = matRad_world2cubeCoords(wCoord, ct) +function coord = matRad_world2cubeCoords(wCoord, gridStruct) % matRad function to convert world coordinates to cube coordinates % % call @@ -6,7 +6,8 @@ % % % wCoord: world coordinates (x,y,z)[mm] -% ct : matRad ct structure +% gridStruct: can be matRad ct, dij.doseGrid, or the ctGrid +% required fields x,y,x,dimensions,resolution % % coord : cube index (x,y,z) % References @@ -25,24 +26,27 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% coord = []; + if nargin == 2 && ~isempty(wCoord) - - if ~ (isfield(ct,'x') && isfield(ct,'y') && isfield(ct,'z') ) - ct = matRad_getWorldAxes(ct); - end + + gridStruct = matRad_getWorldAxes(gridStruct); + % world bounds - boundx = [min(ct.x)-ct.resolution.x/2 max(ct.x)+ct.resolution.x/2]; - boundy = [min(ct.y)-ct.resolution.y/2 max(ct.y)+ct.resolution.y/2]; - boundz = [min(ct.z)-ct.resolution.z/2 max(ct.z)+ct.resolution.z/2]; + boundx = [min(gridStruct.x)-gridStruct.resolution.x/2 max(gridStruct.x)+gridStruct.resolution.x/2]; + boundy = [min(gridStruct.y)-gridStruct.resolution.y/2 max(gridStruct.y)+gridStruct.resolution.y/2]; + boundz = [min(gridStruct.z)-gridStruct.resolution.z/2 max(gridStruct.z)+gridStruct.resolution.z/2]; + % check if queried coordinate is within worldcube if (wCoord(1)>=boundx(1) && wCoord(1)<=boundx(2)) && ... (wCoord(2)>=boundy(1) && wCoord(2)<=boundy(2)) && ... (wCoord(3)>=boundz(1) && wCoord(3)<=boundz(2)) - [~,coord(1)]=min(abs(ct.x-wCoord(1))); - [~,coord(2)]=min(abs(ct.y-wCoord(2))); - [~,coord(3)]=min(abs(ct.z-wCoord(3))); + % find the closest world coord index in gridStruct.x/y/z + % calc absolute differences and locate smallest difference + [~,coord(1)]=min(abs(gridStruct.x-wCoord(1))); + [~,coord(2)]=min(abs(gridStruct.y-wCoord(2))); + [~,coord(3)]=min(abs(gridStruct.z-wCoord(3))); else matRad_cfg = MatRad_Config.instance(); matRad_cfg.dispError('Queried world coordinate is outside the ct'); diff --git a/matRad/matRad_generateStf.m b/matRad/matRad_generateStf.m index f6742f7b1..0cd979710 100644 --- a/matRad/matRad_generateStf.m +++ b/matRad/matRad_generateStf.m @@ -193,9 +193,8 @@ end % get world coordinate system -if ~isfield(ct,'x') - ct = matRad_getWorldAxes(ct); -end +ct = matRad_getWorldAxes(ct); + % Convert linear indices to 3D voxel coordinates [coordsY_vox, coordsX_vox, coordsZ_vox] = ind2sub(ct.cubeDim,V); @@ -224,11 +223,8 @@ % Correct for iso center position. Whit this correction Isocenter is % (0,0,0) [mm] wCoords = matRad_cube2worldCoords([coordsX_vox, coordsY_vox, coordsZ_vox], ct); - - coordsX = wCoords(:,1) - pln.propStf.isoCenter(i,1); - coordsY = wCoords(:,2) - pln.propStf.isoCenter(i,2); - coordsZ = wCoords(:,3) - pln.propStf.isoCenter(i,3); - + Icoords = wCoords - pln.propStf.isoCenter(i,:); + % Save meta information for treatment plan stf(i).gantryAngle = pln.propStf.gantryAngles(i); stf(i).couchAngle = pln.propStf.couchAngles(i); @@ -244,7 +240,7 @@ % rotation matrix are necessary rotMat_system_T = matRad_getRotationMatrix(pln.propStf.gantryAngles(i),pln.propStf.couchAngles(i)); - rot_coords = [coordsX coordsY coordsZ]*rotMat_system_T; + rot_coords = [Icoords]*rotMat_system_T; % project x and z coordinates to isocenter coordsAtIsoCenterPlane(:,1) = (rot_coords(:,1)*SAD)./(SAD + rot_coords(:,2)); @@ -672,7 +668,7 @@ subplot(1,2,2) % Plot target coordinates whitout any rotation - plot3(coordsX,coordsY,coordsZ,'r.') + plot3(Icoords,'r.') hold on; % Rotated projection matrix at isocenter diff --git a/matRad/plotting/matRad_plotAxisLabels.m b/matRad/plotting/matRad_plotAxisLabels.m index e3215c2d7..45967a829 100644 --- a/matRad/plotting/matRad_plotAxisLabels.m +++ b/matRad/plotting/matRad_plotAxisLabels.m @@ -37,9 +37,8 @@ function matRad_plotAxisLabels(axesHandle,ct,plane,slice,defaultFontSize,tickdi defaultFontSize = 12; end -if ~isfield(ct,'x') || isempty(ct.x) - ct = matRad_getWorldAxes(ct); -end +ct = matRad_getWorldAxes(ct); + if ~exist('tickdist','var') || isempty(tickdist) tickdist = abs(ct.x(end)-ct.x(1))/10; diff --git a/matRad/tools/matRad_generateSingleBixelStf.m b/matRad/tools/matRad_generateSingleBixelStf.m index 65a223bbf..04b1d8c72 100644 --- a/matRad/tools/matRad_generateSingleBixelStf.m +++ b/matRad/tools/matRad_generateSingleBixelStf.m @@ -94,18 +94,8 @@ end %Get Isocenter voxel as target -if ~any(isfield(ct,{'x','y','z'})) - ct = matRad_getWorldAxes(ct); +ct = matRad_getWorldAxes(ct); -end - -%xVox = ct.x + ct.resolution.x/2; -%yVox = ct.y + ct.resolution.y/2; -%zVox = ct.z + ct.resolution.z/2; - -% xVox = ct.resolution.x*[1:ct.cubeDim(2)]-ct.resolution.x/2; -% yVox = ct.resolution.y*[1:ct.cubeDim(1)]-ct.resolution.y/2; -% zVox = ct.resolution.z*[1:ct.cubeDim(3)]-ct.resolution.z/2; % loop over all angles for i = 1:length(pln.propStf.gantryAngles) From 4ca359197c6a350107ed9e105019acd8efa6c5d9 Mon Sep 17 00:00:00 2001 From: amitantony Date: Fri, 26 Jul 2024 16:43:12 +0200 Subject: [PATCH 08/10] consistent tranformation to the image coordinate system (isocenter at 0,0 ) --- .../matRad_ParticleMCsquareEngine.m | 2 +- .../matRad_PencilBeamEngineAbstract.m | 2 +- .../+DoseEngines/matRad_PhotonOmpMCEngine.m | 7 +-- .../+DoseEngines/matRad_TopasMCEngine.m | 2 +- matRad/geometry/matRad_world2imageCoords.m | 63 +++++++++++++++++++ matRad/matRad_generateStf.m | 2 +- matRad/tools/matRad_generateSingleBixelStf.m | 2 +- test/core/test_world2imageCoords.m | 31 +++++++++ 8 files changed, 102 insertions(+), 9 deletions(-) create mode 100644 matRad/geometry/matRad_world2imageCoords.m create mode 100644 test/core/test_world2imageCoords.m diff --git a/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m b/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m index db29ccef6..14d736324 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m @@ -272,7 +272,7 @@ function setDefaults(this) for i = 1:length(stf) %Create new stf for MCsquare with energy layer ordering and %shifted scenario isocenter - stfMCsquare(i).isoCenter = stf(i).isoCenter - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z] + isoCenterShift; + stfMCsquare(i).isoCenter = matRad_world2imageCoords(stf(i).isoCenter, ct) + isoCenterShift; stfMCsquare(i).gantryAngle = mod(180-stf(i).gantryAngle,360); %Different MCsquare geometry stfMCsquare(i).couchAngle = stf(i).couchAngle; stfMCsquare(i).energies = unique([stf(i).ray.energy]); diff --git a/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m b/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m index 13f91fe95..924173f5c 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m +++ b/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m @@ -97,7 +97,7 @@ function setDefaults(this) scenStf = stf; % manipulate isocenter for k = 1:numel(scenStf) - scenStf(k).isoCenter = scenStf(k).isoCenter - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z] + this.multScen.isoShift(ixShiftScen,:); + scenStf(k).isoCenter = matRad_world2imageCoords(scenStf(k).isoCenter,ct) + this.multScen.isoShift(ixShiftScen,:); end if this.multScen.totNumShiftScen > 1 diff --git a/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m b/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m index 87324fb27..aa095369b 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m @@ -143,8 +143,8 @@ scenCount = scenCount + 1; % manipulate isocenter - shiftedIsoCenter = vertcat(stf(:).isoCenter) - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z]... - + this.multScen.isoShift(scenarioIx,:) + dij.doseGrid.isoCenterOffset; + shiftedIsoCenter = matRad_world2imageCoords(vertcat(stf(:).isoCenter), ct) + this.multScen.isoShift(scenarioIx,:) + dij.doseGrid.isoCenterOffset; + this.ompMCgeo.isoCenter = shiftedIsoCenter; tmpStf = stf; @@ -225,8 +225,7 @@ tmpStf = stf; for k = 1:length(stf) - shiftedIsoCenter = vertcat(stf(:).isoCenter) - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z]... - + dij.doseGrid.isoCenterOffset; + shiftedIsoCenter = matRad_world2imageCoords(vertcat(stf(:).isoCenter),ct) + dij.doseGrid.isoCenterOffset; tmpStf(k).isoCenter = shiftedIsoCenter; end diff --git a/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m b/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m index 9355f9167..52ec0f86a 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m @@ -489,7 +489,7 @@ function writeAllFiles(obj,ct,cst,stf,machine,w) % manipulate isocenter for k = 1:numel(stf) - stf(k).isoCenter = stf(k).isoCenter - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z] + this.multScen.isoShift(ixShiftScen,:); + stf(k).isoCenter = matRad_world2imageCoords(stf(k).isoCenter,ct) + this.multScen.isoShift(ixShiftScen,:); end % Delete previous topas files so there is no mix-up diff --git a/matRad/geometry/matRad_world2imageCoords.m b/matRad/geometry/matRad_world2imageCoords.m new file mode 100644 index 000000000..15e1fa477 --- /dev/null +++ b/matRad/geometry/matRad_world2imageCoords.m @@ -0,0 +1,63 @@ +function coord = matRad_world2imageCoords(wCoord, gridStruct) +% matRad function to convert world coordinates to Iso Center image coordinates +% voxels step: mm, isocenter at (0,0) +% +% call +% coord = matRad_world2imageCoords(wCoord, ct) +% +% +% wCoord: world coordinates (x,y,z)[mm] +% gridStruct: can be matRad ct, dij.doseGrid, or the ctGrid +% required fields x,y,x,dimensions,resolution +% +% coord : image coordinates (x,y,z)[mm] centered at isocenter +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 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/LICENSES.txt. 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +coord = []; + + +if nargin == 2 && ~isempty(wCoord) + + gridStruct = matRad_getWorldAxes(gridStruct); + + % world bounds + boundx = [min(gridStruct.x)-gridStruct.resolution.x/2 max(gridStruct.x)+gridStruct.resolution.x/2]; + boundy = [min(gridStruct.y)-gridStruct.resolution.y/2 max(gridStruct.y)+gridStruct.resolution.y/2]; + boundz = [min(gridStruct.z)-gridStruct.resolution.z/2 max(gridStruct.z)+gridStruct.resolution.z/2]; + + % check if queried coordinate is within worldcube + if (wCoord(1)>=boundx(1) && wCoord(1)<=boundx(2)) && ... + (wCoord(2)>=boundy(1) && wCoord(2)<=boundy(2)) && ... + (wCoord(3)>=boundz(1) && wCoord(3)<=boundz(2)) + + firstVox = [gridStruct.x(1), gridStruct.y(1), gridStruct.z(1)]; + res = [gridStruct.resolution.x gridStruct.resolution.y gridStruct.resolution.z]; + coord = wCoord - firstVox + res; + else + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Queried world coordinate is outside the ct'); + + end + +% catch +% +% end +else + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Cannot compute image coordinates without matRad input coordinates or ct structure'); +end + +end diff --git a/matRad/matRad_generateStf.m b/matRad/matRad_generateStf.m index 0cd979710..9d75e763d 100644 --- a/matRad/matRad_generateStf.m +++ b/matRad/matRad_generateStf.m @@ -330,7 +330,7 @@ stf(i).numOfBixelsPerRay = ones(1,stf(i).numOfRays); % mm axes with isocenter at (0,0,0) - mmCubeisoCenter = stf(i).isoCenter - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z]; + mmCubeisoCenter = matRad_world2imageCoords(stf(i).isoCenter,ct); for j = stf(i).numOfRays:-1:1 for ShiftScen = 1:pln.multScen.totNumShiftScen diff --git a/matRad/tools/matRad_generateSingleBixelStf.m b/matRad/tools/matRad_generateSingleBixelStf.m index 04b1d8c72..44039542a 100644 --- a/matRad/tools/matRad_generateSingleBixelStf.m +++ b/matRad/tools/matRad_generateSingleBixelStf.m @@ -149,7 +149,7 @@ stf(i).ray.targetPoint_bev = [0 SAD 0]; - shiftedIsoCenter = vertcat(stf(:).isoCenter) - [ct.x(1) ct.y(1) ct.z(1)] + [ct.resolution.x ct.resolution.y ct.resolution.z] ; + shiftedIsoCenter = matRad_world2imageCoords(vertcat(stf(:).isoCenter),ct); stf(i).ray.rayPos = shiftedIsoCenter; stf(i).ray.targetPoint = [0 SAD 0] * rotMat_vectors_T; diff --git a/test/core/test_world2imageCoords.m b/test/core/test_world2imageCoords.m new file mode 100644 index 000000000..e1fdd07fc --- /dev/null +++ b/test/core/test_world2imageCoords.m @@ -0,0 +1,31 @@ +function test_suite = test_world2imageCoords + +test_functions = localfunctions(); + +initTestSuite; + + +function test_empty_world2imageCoords + ct = get_testCtHelper(); + assertExceptionThrown(@()matRad_world2imageCoords()); + assertExceptionThrown(@()matRad_world2imageCoords([],ct)); + + +function test_basic_world2imageCoords + ct = get_testCtHelper(); + expected = [3 36 45]; + assertEqual(matRad_world2imageCoords([-30 3 12], ct),expected); +% out of bounds + assertExceptionThrown(@()matRad_world2imageCoords([-30,4,45],ct)); + + +function ct = get_testCtHelper() + ct.x = -30:3:27; + ct.y = -30:3:27; + ct.z = -30:3:27; + + ct.resolution.x = 3; + ct.resolution.y = 3; + ct.resolution.z = 3; + + ct.cubeDim = [20 20 20 ]; \ No newline at end of file From 483f8d1cb95a975b92f94e3daabd4d46864ecc9b Mon Sep 17 00:00:00 2001 From: amitantony Date: Fri, 26 Jul 2024 17:18:57 +0200 Subject: [PATCH 09/10] renaming it to world2isocentricCoords --- .../matRad_ParticleMCsquareEngine.m | 2 +- .../matRad_PencilBeamEngineAbstract.m | 2 +- .../+DoseEngines/matRad_PhotonOmpMCEngine.m | 4 +-- .../+DoseEngines/matRad_TopasMCEngine.m | 2 +- ...ords.m => matRad_world2isocentricCoords.m} | 4 +-- matRad/matRad_generateStf.m | 2 +- matRad/tools/matRad_generateSingleBixelStf.m | 2 +- test/core/test_world2imageCoords.m | 31 ------------------- test/core/test_world2isocentricCoords.m | 31 +++++++++++++++++++ 9 files changed, 40 insertions(+), 40 deletions(-) rename matRad/geometry/{matRad_world2imageCoords.m => matRad_world2isocentricCoords.m} (95%) delete mode 100644 test/core/test_world2imageCoords.m create mode 100644 test/core/test_world2isocentricCoords.m diff --git a/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m b/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m index 14d736324..b370dc5ad 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_ParticleMCsquareEngine.m @@ -272,7 +272,7 @@ function setDefaults(this) for i = 1:length(stf) %Create new stf for MCsquare with energy layer ordering and %shifted scenario isocenter - stfMCsquare(i).isoCenter = matRad_world2imageCoords(stf(i).isoCenter, ct) + isoCenterShift; + stfMCsquare(i).isoCenter = matRad_world2isocentricCoords(stf(i).isoCenter, ct) + isoCenterShift; stfMCsquare(i).gantryAngle = mod(180-stf(i).gantryAngle,360); %Different MCsquare geometry stfMCsquare(i).couchAngle = stf(i).couchAngle; stfMCsquare(i).energies = unique([stf(i).ray.energy]); diff --git a/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m b/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m index 924173f5c..e72e438e7 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m +++ b/matRad/doseCalc/+DoseEngines/matRad_PencilBeamEngineAbstract.m @@ -97,7 +97,7 @@ function setDefaults(this) scenStf = stf; % manipulate isocenter for k = 1:numel(scenStf) - scenStf(k).isoCenter = matRad_world2imageCoords(scenStf(k).isoCenter,ct) + this.multScen.isoShift(ixShiftScen,:); + scenStf(k).isoCenter = matRad_world2isocentricCoords(scenStf(k).isoCenter,ct) + this.multScen.isoShift(ixShiftScen,:); end if this.multScen.totNumShiftScen > 1 diff --git a/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m b/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m index aa095369b..5504d8b9c 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_PhotonOmpMCEngine.m @@ -143,7 +143,7 @@ scenCount = scenCount + 1; % manipulate isocenter - shiftedIsoCenter = matRad_world2imageCoords(vertcat(stf(:).isoCenter), ct) + this.multScen.isoShift(scenarioIx,:) + dij.doseGrid.isoCenterOffset; + shiftedIsoCenter = matRad_world2isocentricCoords(vertcat(stf(:).isoCenter), ct) + this.multScen.isoShift(scenarioIx,:) + dij.doseGrid.isoCenterOffset; this.ompMCgeo.isoCenter = shiftedIsoCenter; tmpStf = stf; @@ -225,7 +225,7 @@ tmpStf = stf; for k = 1:length(stf) - shiftedIsoCenter = matRad_world2imageCoords(vertcat(stf(:).isoCenter),ct) + dij.doseGrid.isoCenterOffset; + shiftedIsoCenter = matRad_world2isocentricCoords(vertcat(stf(:).isoCenter),ct) + dij.doseGrid.isoCenterOffset; tmpStf(k).isoCenter = shiftedIsoCenter; end diff --git a/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m b/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m index 52ec0f86a..c09ce2a1a 100644 --- a/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m +++ b/matRad/doseCalc/+DoseEngines/matRad_TopasMCEngine.m @@ -489,7 +489,7 @@ function writeAllFiles(obj,ct,cst,stf,machine,w) % manipulate isocenter for k = 1:numel(stf) - stf(k).isoCenter = matRad_world2imageCoords(stf(k).isoCenter,ct) + this.multScen.isoShift(ixShiftScen,:); + stf(k).isoCenter = matRad_world2isocentricCoords(stf(k).isoCenter,ct) + this.multScen.isoShift(ixShiftScen,:); end % Delete previous topas files so there is no mix-up diff --git a/matRad/geometry/matRad_world2imageCoords.m b/matRad/geometry/matRad_world2isocentricCoords.m similarity index 95% rename from matRad/geometry/matRad_world2imageCoords.m rename to matRad/geometry/matRad_world2isocentricCoords.m index 15e1fa477..c920d41c3 100644 --- a/matRad/geometry/matRad_world2imageCoords.m +++ b/matRad/geometry/matRad_world2isocentricCoords.m @@ -1,9 +1,9 @@ -function coord = matRad_world2imageCoords(wCoord, gridStruct) +function coord = matRad_world2isocentricCoords(wCoord, gridStruct) % matRad function to convert world coordinates to Iso Center image coordinates % voxels step: mm, isocenter at (0,0) % % call -% coord = matRad_world2imageCoords(wCoord, ct) +% coord = matRad_world2isocentricCoords(wCoord, ct) % % % wCoord: world coordinates (x,y,z)[mm] diff --git a/matRad/matRad_generateStf.m b/matRad/matRad_generateStf.m index 9d75e763d..16f9ad4cf 100644 --- a/matRad/matRad_generateStf.m +++ b/matRad/matRad_generateStf.m @@ -330,7 +330,7 @@ stf(i).numOfBixelsPerRay = ones(1,stf(i).numOfRays); % mm axes with isocenter at (0,0,0) - mmCubeisoCenter = matRad_world2imageCoords(stf(i).isoCenter,ct); + mmCubeisoCenter = matRad_world2isocentricCoords(stf(i).isoCenter,ct); for j = stf(i).numOfRays:-1:1 for ShiftScen = 1:pln.multScen.totNumShiftScen diff --git a/matRad/tools/matRad_generateSingleBixelStf.m b/matRad/tools/matRad_generateSingleBixelStf.m index 44039542a..f2d876c61 100644 --- a/matRad/tools/matRad_generateSingleBixelStf.m +++ b/matRad/tools/matRad_generateSingleBixelStf.m @@ -149,7 +149,7 @@ stf(i).ray.targetPoint_bev = [0 SAD 0]; - shiftedIsoCenter = matRad_world2imageCoords(vertcat(stf(:).isoCenter),ct); + shiftedIsoCenter = matRad_world2isocentricCoords(vertcat(stf(:).isoCenter),ct); stf(i).ray.rayPos = shiftedIsoCenter; stf(i).ray.targetPoint = [0 SAD 0] * rotMat_vectors_T; diff --git a/test/core/test_world2imageCoords.m b/test/core/test_world2imageCoords.m deleted file mode 100644 index e1fdd07fc..000000000 --- a/test/core/test_world2imageCoords.m +++ /dev/null @@ -1,31 +0,0 @@ -function test_suite = test_world2imageCoords - -test_functions = localfunctions(); - -initTestSuite; - - -function test_empty_world2imageCoords - ct = get_testCtHelper(); - assertExceptionThrown(@()matRad_world2imageCoords()); - assertExceptionThrown(@()matRad_world2imageCoords([],ct)); - - -function test_basic_world2imageCoords - ct = get_testCtHelper(); - expected = [3 36 45]; - assertEqual(matRad_world2imageCoords([-30 3 12], ct),expected); -% out of bounds - assertExceptionThrown(@()matRad_world2imageCoords([-30,4,45],ct)); - - -function ct = get_testCtHelper() - ct.x = -30:3:27; - ct.y = -30:3:27; - ct.z = -30:3:27; - - ct.resolution.x = 3; - ct.resolution.y = 3; - ct.resolution.z = 3; - - ct.cubeDim = [20 20 20 ]; \ No newline at end of file diff --git a/test/core/test_world2isocentricCoords.m b/test/core/test_world2isocentricCoords.m new file mode 100644 index 000000000..b1e3ca702 --- /dev/null +++ b/test/core/test_world2isocentricCoords.m @@ -0,0 +1,31 @@ +function test_suite = test_world2isocentricCoords + +test_functions = localfunctions(); + +initTestSuite; + + +function test_empty_world2isocentricCoords + ct = get_testCtHelper(); + assertExceptionThrown(@()matRad_world2isocentricCoords()); + assertExceptionThrown(@()matRad_world2isocentricCoords([],ct)); + + +function test_basic_world2isocentricCoords + ct = get_testCtHelper(); + expected = [3 36 45]; + assertEqual(matRad_world2isocentricCoords([-30 3 12], ct),expected); +% out of bounds + assertExceptionThrown(@()matRad_world2isocentricCoords([-30,4,45],ct)); + + +function ct = get_testCtHelper() + ct.x = -30:3:27; + ct.y = -30:3:27; + ct.z = -30:3:27; + + ct.resolution.x = 3; + ct.resolution.y = 3; + ct.resolution.z = 3; + + ct.cubeDim = [20 20 20 ]; \ No newline at end of file From 176509ba375760471cd752623f70eea16ace42a1 Mon Sep 17 00:00:00 2001 From: amitantony Date: Tue, 6 Aug 2024 18:14:59 +0200 Subject: [PATCH 10/10] bugfix --- matRad/geometry/matRad_cube2worldCoords.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matRad/geometry/matRad_cube2worldCoords.m b/matRad/geometry/matRad_cube2worldCoords.m index 4f755dc5a..bbabe494a 100644 --- a/matRad/geometry/matRad_cube2worldCoords.m +++ b/matRad/geometry/matRad_cube2worldCoords.m @@ -32,7 +32,7 @@ if isfield(gridStruct, 'dicomInfo') && isfield(gridStruct.dicomInfo,'ImagePositionPatient') firstVox = gridStruct.dicomInfo.ImagePositionPatient; else - firstVox = -(gridStruct.cubeDim./2).*[gridStruct.resolution.x gridStruct.resolution.y gridStruct.resolution.z] ; + firstVox = -(gridStruct.dimensions./2).*[gridStruct.resolution.x gridStruct.resolution.y gridStruct.resolution.z] ; end if prod(prod(vCoord<=gridStruct.dimensions)) && prod(prod(vCoord>=0)) coord(1,:) = firstVox(1) + (vCoord(:,1) - 1 ) *gridStruct.resolution.x;