Dev varRBE robopt world coordinate system

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,

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')

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);
 

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)
 

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')

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);
 

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)
 

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);
 

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)
 

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);
 

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')

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 
 

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;

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

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
 
 

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);

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]);

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

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;
 
@@ -216,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
 

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

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

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)