e0404 · wahln · Aug 7, 2024 · Jul 22, 2024 · Jul 22, 2024 · Jul 22, 2024
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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)
 
@@ -140,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);
 

diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -31,14 +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);
-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
@@ -91,14 +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);
-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
@@ -78,15 +78,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;
-end
+ct = matRad_getWorldAxes(ct);
+
+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
@@ -250,6 +250,10 @@ function setDefaults(this)
                 this.config.LET_Sparse_Output	 = ~this.calcDoseDirect;
             end
 
+            %Create X Y Z vectors if not present
+            ct = matRad_getWorldAxes(ct);
+
+
             for scenarioIx = 1:this.multScen.totNumScen
                 %For direct dose calculation
                 totalWeights = 0;
@@ -268,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 + 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
@@ -85,6 +85,10 @@ function setDefaults(this)
             % initialize
             dij = this.initDoseCalc(ct,cst,stf);
 
+            %Create X Y Z vectors if not present
+            ct = matRad_getWorldAxes(ct);
+
+
             for shiftScen = 1:this.multScen.totNumShiftScen
 
                 %Find first instance of the shift to select the shift values
@@ -93,7 +97,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 = 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
@@ -127,6 +127,10 @@
 
             %Now we have to calibrate to the the beamlet width.
             calibrationFactor = this.absCalibrationFactor * (bixelWidth/50)^2;
+
+            %Create X Y Z vectors if not present
+            ct = matRad_getWorldAxes(ct);
+
 
             scenCount = 0;
             %run over all scenarios
@@ -139,7 +143,8 @@
                     scenCount = scenCount + 1;
 
                     % manipulate isocenter
-                    shiftedIsoCenter = vertcat(stf(:).isoCenter) + 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;
 
@@ -216,7 +221,16 @@
             this.getOmpMCgeometry(dij.doseGrid);
 
             %% Create beamlet source
-            this.getOmpMCsource(stf);
+            ct = matRad_getWorldAxes(ct);
+
+            tmpStf = stf;
+             for k = 1:length(stf)
+                 shiftedIsoCenter = matRad_world2isocentricCoords(vertcat(stf(:).isoCenter),ct) + dij.doseGrid.isoCenterOffset;
+
+                 tmpStf(k).isoCenter = shiftedIsoCenter;
+            end
+
+            this.getOmpMCsource(tmpStf);
         end
     end
 

diff --git 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 = 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_projectOnComponents.m → matRad/doseCalc/matRad_projectOnComponents.m b/matRad/geometry/matRad_projectOnComponents.m → 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
@@ -0,0 +1,52 @@
+function coord = matRad_cube2worldCoords(vCoord, gridStruct)
+% matRad function to convert cube coordinates to world coordinates
+% 
+% call
+%   coord = matRad_worldToCubeCoordinates(vCoord, gridStruct)
+%   
+%   gridStruct: matRad ct struct or dij.doseGrid/ctGrid struct
+%   vCoord : Voxel Coordinates [vx vy vz]
+%   coord: worldCoordinates [x y z ] mm
+% 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 && ~isempty(vCoord)
+    if isfield(gridStruct,'cubeDim')
+        gridStruct.dimensions = gridStruct.cubeDim;
+    end
+
+    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 
+    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();
+        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
diff --git 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)