Integration of APM / probabilistic optimization

APM/MEX/matRad_calcFourthRangeMom.m

@@ -0,0 +1,126 @@
+function [fourthMom] = matRad_calcFourthRangeMom(numOfSpots,numComp,w,mMeanA,mMeanB,mMeanC,mMeanD,...
+                                                                      mWidthA,mWidthB,mWidthC,mWidthD,...
+                                                                      mWeightA,mWeightB,mWeightC,mWeightD,...
+                                                                      voxelPos,mSysCovRadDept) %#codegen
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% matRad function to calculate the fourth raw moment. 
+% 
+% call
+%   [fourthMom] = matRad_calcFourthRangeMom(numOfSpots,numComp,w,mMeanA,mMeanB,mMeanC,mMeanD, ...
+%                                                                mWidthA,mWidthB,mWidthC,mWidthD,...
+%                                                                mWeightA,mWeightB,mWeightC,mWeightD,...
+%                                                                voxelPos,mSysCovRadDept)
+%
+% input
+%   numOfSpots           total number of spots
+%   numComp              total number of Gaussian components e.g. 10 or 13
+%   mMeanA-mMeanD:      four mean vectors of spot components j o q n
+%   mWidthA-mWidthD:     four width vectors of spot components j o q n
+%   mWeightA-mWeightD:   four weight vectors of spot components j o q n
+%   voxelPos:            voxel position 
+%   mSysCovRadDept:      full covariance matrix of range error
+%
+% output
+%   fourthMom:           fourth central moment
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Copyright 2018 Hans-Peter Wieser
+% 
+% 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.
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+assert(isa(numOfSpots,'double'));
+assert(isa(numComp,'double'));
+assert(isa(w,'double'));
+assert(isa(mMeanA,'double'));
+assert(isa(mMeanB,'double'));
+assert(isa(mMeanC,'double'));
+assert(isa(mMeanD,'double'));
+assert(isa(mWidthA,'double'));
+assert(isa(mWidthB,'double'));
+assert(isa(mWidthC,'double'));
+assert(isa(mWidthD,'double'));
+assert(isa(mWeightA,'double'));
+assert(isa(mWeightB,'double'));
+assert(isa(mWeightC,'double'));
+assert(isa(mWeightD,'double'));
+assert(isa(voxelPos,'double'));
+assert(isa(mSysCovRadDept,'double'));
+
+coder.varsize('numOfSpots',      [1 1], [0 0]);
+coder.varsize('numComp',         [1 1], [0 0]);
+coder.varsize('w',               [60 1],[1 0]);
+coder.varsize('mMeanA',          [15 60],[1 1]);
+coder.varsize('mMeanB',          [15 60],[1 1]);
+coder.varsize('mMeanC',          [15 60],[1 1]);
+coder.varsize('mMeanD',          [15 60],[1 1]);
+coder.varsize('mWidthA',         [15 60],[1 1]);
+coder.varsize('mWidthB',         [15 60],[1 1]);
+coder.varsize('mWidthD',         [15 60],[1 1]);
+coder.varsize('mWidthD',         [15 60],[1 1]);
+coder.varsize('mWeightA',        [15 60],[1 1]);
+coder.varsize('mWeightB',        [15 60],[1 1]);
+coder.varsize('mWeightC',        [15 60],[1 1]);
+coder.varsize('mWeightD',        [15 60],[1 1]);
+coder.varsize('voxelPos',        [1 1],[0 0]);
+coder.varsize('mSysCovRadDept',  [60 60],[1 1]);
+
+
+mPSI_joqn = zeros(numOfSpots,numOfSpots,numOfSpots,numOfSpots);
+
+for j = 1:numOfSpots
+    for o = 1:numOfSpots
+        for q = 1:numOfSpots
+           for n = 1:numOfSpots
+
+            sigma = [mSysCovRadDept(j,j)  mSysCovRadDept(j,o)    mSysCovRadDept(j,q)   mSysCovRadDept(j,n);...
+                     mSysCovRadDept(o,j)  mSysCovRadDept(o,o)    mSysCovRadDept(o,q)   mSysCovRadDept(o,n);...
+                     mSysCovRadDept(q,j)  mSysCovRadDept(q,o)    mSysCovRadDept(q,q)   mSysCovRadDept(q,n);...
+                     mSysCovRadDept(n,j)  mSysCovRadDept(n,o)    mSysCovRadDept(n,q)   mSysCovRadDept(n,n)];
+
+            dev_j   = voxelPos - mMeanA(:,j);
+            dev_o   = voxelPos - mMeanB(:,o);
+            dev_q   = voxelPos - mMeanC(:,q);
+            dev_n   = voxelPos - mMeanD(:,n);
+
+            vW = reshape(reshape(reshape(mWeightA(:,j) * mWeightB(:,o)',[],1) * mWeightC(:,q)',[],1) * mWeightD(:,n)',[numComp,numComp,numComp,numComp]);
+
+            Y = 0;
+            for J = 1:numComp
+                 for O = 1:numComp
+                      for Q = 1:numComp  
+                          for N = 1:numComp
+
+                            lambda = diag([mWidthA(J,j) mWidthB(O,o) mWidthC(Q,q) mWidthD(N,n)]);  
+
+                            LaSi    = (lambda + sigma);
+                            InvLaSi = inv(LaSi);
+
+                             Y = Y + (vW(J,O,Q,N)/((2*pi)^(4/2)*sqrt(det(LaSi)))) * ...
+                                 exp(-0.5*([dev_j(J) dev_o(O) dev_q(Q) dev_n(N)]) * (InvLaSi) * ([dev_j(J) dev_o(O) dev_q(Q) dev_n(N)])');
+
+                          end
+                      end
+                 end
+            end
+
+            mPSI_joqn(j,o,q,n) = Y;
+
+           end
+        end
+    end
+end
+
+
+fourthMom =  w'* reshape(reshape(reshape(mPSI_joqn,[numOfSpots^3 numOfSpots]) * w, [numOfSpots^2 numOfSpots]) * w , [numOfSpots numOfSpots]) * w;    
+
+
+end

APM/MEX/matRad_calcGeoDistsFast.m

@@ -0,0 +1,116 @@
+function [isoLatDistsX,isoLatDistsZ] = matRad_calcGeoDistsFast(rot_coords_bev, ...
+                                                                sourcePoint_bev, ...
+                                                                targetPoint_bev, ...
+                                                                SAD, ...
+                                                                radDepthIx, ...
+                                                                lateralCutOff) %#codegen                        
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% matRad calculation of lateral distances from central ray used for
+% dose calcultion
+% 
+% call
+%   [ix,x_latDists,z_latDists] = ...
+%           matRad_calcGeoDists(rot_coords_bev, ...
+%                               sourcePoint_bev, ...
+%                               targetPoint_bev, ...
+%                               SAD, ...
+%                               radDepthIx, ...
+%                               lateralCutOff)
+%
+% input
+%   rot_coords_bev:     coordinates in bev of the voxels with index V,
+%                       where also ray tracing results are availabe 
+%   sourcePoint_bev:    source point in voxel coordinates in beam's eye view
+%   targetPoint_bev:    target point in voxel coordinated in beam's eye view
+%   SAD:                source-to-axis distance
+%   radDepthIx:         sub set of voxels for which radiological depth
+%                       calculations are available
+%   lateralCutOff:      lateral cutoff specifying the neighbourhood for
+%                       which dose calculations will actually be performed
+%
+% output
+%   ix:                 indices of voxels where we want to compute dose
+%                       influence data
+%   isoLatDistsX:       lateral x-distance to the central ray projected to
+%                       iso center plane
+%   isoLatDistsZ:       lateral z-distance to the central ray projected to
+%                       iso center plane
+%   radialDist_sq:      squared radial distance to the central ray (where the
+%                       actual computation of the radiological depth takes place)
+%
+% 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.
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% ROTATE A SINGLE BEAMLET AND ALIGN WITH BEAMLET WHO PASSES THROUGH
+% ISOCENTER
+
+
+assert(isa(rot_coords_bev,'double'));
+assert(isa(sourcePoint_bev,'double'));  
+assert(isa(targetPoint_bev,'double'));  
+assert(isa(SAD,'double'));  
+assert(isa(radDepthIx,'double'));  
+assert(isa(lateralCutOff,'double'));  
+
+coder.varsize('rot_coords_bev',[1 3]);
+coder.varsize('sourcePoint_bev',[1 3]);
+coder.varsize('targetPoint_bev',[1 3]);
+coder.varsize('SAD',1);
+coder.varsize('radDepthIx',1);
+coder.varsize('lateralCutOff',1);    
+
+% Put [0 0 0] position in the source point for beamlet who passes through
+% isocenter
+a = -sourcePoint_bev';
+
+% Normalize the vector
+a = a/norm(a);
+
+% Put [0 0 0] position in the source point for a single beamlet
+b = (targetPoint_bev - sourcePoint_bev)';
+
+% Normalize the vector
+b = b/norm(b);
+
+% Define function for obtain rotation matrix.
+if all(a==b) % rotation matrix corresponds to eye matrix if the vectors are the same
+    rot_coords_temp = rot_coords_bev;
+else
+    % Define rotation matrix
+    v   = cross(a,b);
+    ssc =  [0 -v(3) v(2); v(3) 0 -v(1); -v(2) v(1) 0];
+    R   = eye(3) + ssc + ssc^2*(1-dot(a,b))/(norm(cross(a,b))^2);
+
+    % Rotate every CT voxel 
+    rot_coords_temp = rot_coords_bev*R;
+end
+
+% Put [0 0 0] position CT in center of the beamlet.
+latDistsX = rot_coords_temp(:,1) + sourcePoint_bev(1);
+latDistsZ = rot_coords_temp(:,3) + sourcePoint_bev(3);
+
+% check of radial distance exceeds lateral cutoff (projected to iso center)
+rad_distancesSq = latDistsX.^2 + latDistsZ.^2;
+subsetMask = rad_distancesSq ./ rot_coords_temp(:,2).^2 <= lateralCutOff^2 /SAD^2;
+
+isoLatDistsX = latDistsX(subsetMask)./rot_coords_temp(subsetMask,2)*SAD;
+isoLatDistsZ = latDistsZ(subsetMask)./rot_coords_temp(subsetMask,2)*SAD; 
+
+
+

APM/MEX/matRad_calcSecLatMom.m

@@ -0,0 +1,58 @@
+function PSI_ijlm =  matRad_calcSecLatMom(vLaSi11,vLaSi22,vLaSi12,vLaSi21,Dev_j,Dev_m) %#codegen
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% matRad function to calculate the second raw moment. Please note that i
+% and l depict voxel indices and j and m pencil beam indices. This function
+% allows to calculate the the correlation between spot j and multiple
+% spots m simultaniously. 
+% 
+% call
+%   PSI_ijlm =  matRad_calcSecLatMom(vLaSi11,vLaSi22,vLaSi12,vLaSi21,Dev_j,Dev_m)
+%
+% input
+%   vLaSi11:        matRads resultGUI struct
+%   vLaSi22:        matRad plan meta information struct
+%   vLaSi12:        matRad dij dose influence struct
+%   vLaSi21:        matRad steering information struct
+%   Dev_j:          matRad critical structure struct cst
+%   Dev_m
+
+% output
+%   PSI_ijlm:         matRad critical structure struct cst
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Copyright 2018 the matRad development team. 
+% 
+% This file is part of the matRad project. It is subject to the license 
+% terms in the LICENSE file found in the top-level directory of this 
+% distribution and at https://github.com/e0404/matRad/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.
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+assert(isa(vLaSi11,'double'));
+assert(isa(vLaSi22,'double'));  
+assert(isa(vLaSi12,'double'));  
+assert(isa(vLaSi21,'double'));  
+assert(isa(Dev_j,'double'));  
+assert(isa(Dev_m,'double'));  
+
+coder.varsize('vLaSi11',[1 1]);
+coder.varsize('vLaSi22',[Inf 1],[1 0]);
+coder.varsize('vLaSi12',[Inf 1],[1 0]);
+coder.varsize('vLaSi21',[Inf 1],[1 0]);
+coder.varsize('Dev_j',[1 1]);
+coder.varsize('Dev_m',[Inf 1],[1 0]);
+
+Det           = abs((vLaSi11 .* vLaSi22) - (vLaSi12 .* vLaSi21));
+FracDet       = (1./(2*pi.*real(sqrt(Det))));
+ExpTerm       = -.5 .* ((Dev_j.*(vLaSi22./Det) + Dev_m.*(-vLaSi12./Det)).* Dev_j +...
+                        (Dev_j.*(-vLaSi21./Det) + Dev_m.*(vLaSi11./Det)).* Dev_m);
+PSI_ijlm      = FracDet .* exp(ExpTerm);
+
+end
+