Added the guts

ignored the bin files

.gitignore

@@ -0,0 +1,2 @@
+
+*.bin

TRT_TPS_srb.asv

@@ -0,0 +1,160 @@
+function [BED_Map,DoseMap] =  TRT_TPS()
+cd('C:\Documents and Settings\Joe\My Documents\Research\Phd Project\Multimodal Radiation Therapy')
+I124half = 100.32; %[hrs]
+I131half = 192.96; %[hrs]
+
+lambda_I124 = log(2)/I124half;
+lambda_I131 = log(2)/I131half;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%The user locates the respective files to be used
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+[file, pathname] = uigetfile({'*.mat'},'Please locate the DVK file.','Multiselect','off');
+I131_DVK_3x3x3mm=load(fullfile(pathname,file));
+I131_DVK_3x3x3mm=I131_DVK_3x3x3mm.I131_DVK_3x3x3mm;
+
+[file, pathname] = uigetfile({'*.am'},'Please locate the Amira planning CT file.','Multiselect','off');
+temp_CT = am2geom(fullfile(pathname,file));
+CT = temp_CT.data;
+
+[file, pathname] = uigetfile({'*.am'},'Please locate the Amira ROI files.','Multiselect','on');
+
+Contours = {};
+
+for t=1:numel(file)
+    temp_ROI = am2geom(fullfile(pathname,file{t}));
+    Contours{t}=[];
+    [tmppath, tmpname] = fileparts(file{t});
+    Contours{t}.name = tmpname;
+    [tmppath, tmpname] = fileparts(file{t});
+    Contours{t}.voxel_size = temp_ROI.voxel_size;
+    Contours{t}.start = temp_ROI.start;
+    Contours{t}.data = temp_ROI.data;
+    Contours{t}.alpha_beta = input(['What is the alpha/beta ratio for ',tmpname,'(in Gy)? ']);
+    Contours{t}.repair_constant = input(['What is the repair constant for ',tmpname,'(in hr^-1)? ']);
+    clear temp_ROI;
+end
+
+[file, pathname] = uigetfile({'*.am'},'Please locate the Amira PET files.','Multiselect','on');
+
+for t=1:numel(file)
+    temp_PET = am2geom(fullfile(pathname,file{t}));
+    PET(:,:,:,t) = temp_PET.data;
+end
+voxel_size = temp_PET.voxel_size(1)* temp_PET.voxel_size(2)* temp_PET.voxel_size(3)/1000; %voxel volume [ml]
+start_point = temp_PET.start;
+
+clear temp_PET;
+
+for r = 1:size(PET,4) 
+    temp_time(r) = input(['What is post injection time of PET dataset ',num2str(r),'(in hrs)?']);
+end
+
+nbins = input('Into how many time bins should the data be broken? ');
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Need to adjust the time vector to include t=0 and t=infinite
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+for s = 1:size(PET,4)+2 %Need to add 2 time points, one before and one after the initial PET timepoints
+    if s==1
+        time(s) = 0; %time(1)=0 for the extrapolation step
+    elseif s < size(PET,4)+1 || s == size(PET,4)+1
+        time(s) = temp_time(s-1);
+    else
+        time(s) = temp_time(s-2)*4; %time(n+2) is going to be a decayed data set
+    end
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Finds dose rate maps for I-131 map (PET needs to be converted to I-131 maps first)
+%Need to extrapolate to t=0 and t=infinite - [MBq]*[mGy/MBq-sec] = [mGy/sec]
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+I131_MBq = zeros(size(PET,1),size(PET,2),size(PET,3),size(PET,4)+2);
+DoseRateMap = zeros(size(I131_MBq));
+
+% temp_time
+% time
+% size(PET,4)
+
+for p = 1:numel(time)
+    if p == 1 %Defines time(1)=0 hrs, need extrapolate back to t=0 hrs for each voxel
+        for i=1:size(I131_MBq,1)
+            for j=1:size(I131_MBq,2)
+                for k=1:size(I131_MBq,3)
+                    PET_extrap = interp1(temp_time,squeeze(PET(i,j,k,:)*voxel_size),time,'linear','extrap');
+                    I131_MBq(i,j,k,p) = PET_extrap(1); %Assign the first vector element to create extrapolated t=0
+                end
+            end
+        end
+%         disp('pre-convolution')
+%         size(I131_MBq(:,:,:,p))
+%         size(I131_DVK_3x3x3mm)
+        DoseRateMap(:,:,:,p) = convfft(I131_MBq(:,:,:,p),I131_DVK_3x3x3mm); %Computes the initial dose-rate for t=0
+    elseif p < numel(time)
+        I131_MBq(:,:,:,p) = PET(:,:,:,p-1)*exp(-lambda_I131*time(p))/exp(-lambda_I124*time(p))*voxel_size; %Need to correct for difference in isotope decay
+        DoseRateMap(:,:,:,p) = convfft(I131_MBq(:,:,:,p),I131_DVK_3x3x3mm);
+    else
+        I131_MBq(:,:,:,p) = PET(:,:,:,p-2)*exp(-lambda_I131*time(p))/exp(-lambda_I124*time(p))*exp(-lambda_I131*(time(p)-time(p-1)))*voxel_size; %Simulates decayed data
+        DoseRateMap(:,:,:,p) = convfft(I131_MBq(:,:,:,p),I131_DVK_3x3x3mm);
+    end
+end
+
+% explicit convolutions must be performed between full matrices since
+% spatial dependence needs to be preserved, so this can't be converted to
+% sparse unless you take the Fourier transforms of each and reduce the
+% convolution to a multiplication
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Integrates the DoseRateMap over time to compute DoseMap - [mGy/sec]*[sec] = [mGy]
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+DoseRate = [];
+DoseRate.time = time;
+DoseRate.map = {};
+
+for i=1:size(DoseRateMap,4)
+    DoseRate.map{i} = sparse3D(DoseRateMap(:,:,:,i));
+end
+
+DoseMap = zeros(size(I131_MBq,1),size(I131_MBq,2),size(I131_MBq,3));
+size(DoseMap)
+
+
+time_in = time.*3600; %need to convert the time into seconds for proper integration [secs]
+
+for i=1:numel(time_in);
+    data_in(i) = ; %integration is accomplished for one voxel at a time [mGy/sec] vector
+last_data = DoseRateMap(i,j,k,numel(time)); %[mGy/sec]
+DoseMap(i,j,k)= trapzf(time_in,data_in) + last_data/lambda_I131*3600; %Need to make sure any residual dose is included [mGy/sec]*[sec]
+clear data_in;
+clear last_data;
+
+
+% simple trapezoidal rieman sum integration can be performed on sparse
+% matrices
+
+%size(DoseMap)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Calls PieceWiseBED to the G factor calculation and then creates the BED
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+BED_Map = zeros(size(DoseMap));
+G_factor = zeros(size(DoseMap));
+
+mu = ones(size(DoseMap));
+alpha_beta = ones(size(DoseMap));
+
+for i=1:size(DoseMap,1)
+    for j=1:size(DoseMap,2)
+        for k=1:size(DoseMap,3)
+            G_factor = PieceWiseBED(mu(i,j,k),squeeze(DoseRateMap(i,j,k,:)),DoseMap(i,j,k),nbins,time,p); %Calls PieceWiseBED
+            BED_Map(i,j,k) = DoseMap(i,j,k)*(1 + G_factor./alpha_beta(i,j,k)*DoseMap(i,j,k)); %Uses G to compute BED
+            clear G_factor;
+        end
+    end
+end
+
+save('BED_Map','BED_Map');
+save('DoseMap','DoseMap');
+
+% fid = fopen(['DoseMap' '.img'],'w');
+% fwrite(fid,DoseMap, 'float'); %Pinnacle requires the data to be in 'float' format
+% fclose(fid);

WiscPlanMenu.m

@@ -0,0 +1,141 @@
+%% WiscPlan Menu
+
+% Stephen R Bowen April 2009
+%
+% simple GUI to organize scripts and direct users to appropriate parts of
+% WiscPlan treatment planning system
+
+clear all;
+
+addpath('./photonDoseCalc/scripts_m');
+addpath('./opt/linlsqUtil');
+
+modules = menu('Select WiscPlan Module','Geometry Definition','Tomotherapy Dose Calculator','IMRT Optimization','Prescription Definition','Exit');
+    
+    if modules == 1
+        geomFiles = menu('Select Geometry Definition Function','Import Planning CT & ROIs from Pinnacle','Import Planning CT & ROIs from Amira','Crop & Resample Geometry Matrix','Exit');
+            
+            if geomFiles == 1               
+                pinn2mat = menu(['Script reads in planning CT .img file and plan.roi file, converts '...
+                                  'them to a Matlab structure called Geometry and saves the structure to a .mat file'],'Edit','Run','Exit');
+
+                    if pinn2mat == 1
+                        edit('geometry/pinn2matlab.m');
+                    
+                    elseif pinn2mat == 2
+                        run('geometry/pinn2matlab.m');
+                                               
+                    elseif pinn2mat == 3                        
+                        close(figure);
+                        
+                    end
+
+            elseif geomFiles == 2
+                am2mat = menu(['Script reads in Amira Mesh File .am for planning CT and individual ROI masks,'...
+                                ' converts them to a Matlab structure called Geometry, and saves the structure to a .mat file'],'Edit','Run','Exit');
+                
+                    if am2mat == 1
+                        edit('geometry/GeometryGenerator.m');
+                    
+                    elseif am2mat == 2
+                        run('geometry/GeometryGenerator.m');
+                                               
+                    elseif am2mat == 3                        
+                        close(figure);
+                        
+                    end
+                
+            elseif geomFiles == 3
+                matResamp = menu(['Script reads in Geometry structure from .mat file, crops & resamples '...
+                                  'the matrix and saves the new structure as .mat file'],'Edit','Run','Exit');       
+                                       
+                    if matResamp == 1
+                        edit('geometry/matlabPatientResample.m')
+                    
+                    elseif matResamp == 2
+                        run('geometry/matlabPatientResample.m')
+                        
+                    elseif matResamp == 3
+                        close(figure);
+                        
+                    end
+                    
+            elseif geomFiles == 4
+                close(figure);
+                       
+            end
+
+        
+    elseif modules == 2
+        tomoDoseCalc = menu(['Script reads in Geometry structure from .mat file and either calculates dose transfer matrix from photon beamlets '...
+                             'in Matlab, or writes setup files necessary to run dose calc on Condor or Windows XP'],'Edit','Run','Exit');
+      
+            if tomoDoseCalc == 1
+                edit('photonDoseCalc/helicalDoseCalcSetup.m');
+                
+            elseif tomoDoseCalc == 2
+                run('photonDoseCalc/helicalDoseCalcSetup.m');
+                
+            elseif tomoDoseCalc == 3
+                close(figure);
+                
+            end
+        
+    elseif modules == 3
+
+        optSetup = menu(['Script reads in Geometry structure from .mat file and either runs dose optimization '...
+                                 'in Matlab, or writes setup files necessary to run opt on Condor or Windows XP'],'Edit','Run','Exit');
+                    
+            if optSetup == 1
+                    edit('opt/linlsqOptSetup.m');
+                        
+            elseif optSetup == 2
+                    run('opt/linlsqOptSetup.m');
+                        
+            elseif optSetup == 3
+                    close(figure);                     
+                        
+            end
+    
+        
+    elseif modules == 4
+        prescFiles = menu('Select Prescription Setup Option','Dose Painting Presc Setup','TRT Presc Setup','Exit');
+                
+            if prescFiles == 1
+                dosePaint = menu(['Script reads in  prescription .bin files and replaces uniform dose to PTVs '...
+                                 'with non-uniform PET-based prescriptions'],'Edit','Run','Exit');
+                    
+                    if dosePaint == 1
+                        edit('opt/linlsqOpt/input/dosePaintSetup.m');
+                        
+                    elseif dosePaint == 2
+                        run('opt/linlsqOpt/input/dosePaintSetup.m');
+                        
+                    elseif dosePaint == 3
+                        close(figure);                      
+                        
+                    end
+                
+            elseif prescFiles == 2
+                TRT = menu('Script reads in PTV prescription .bin files and subtracts non-uniform TRT dose','Edit','Run','Exit');
+                    
+                    if TRT == 1
+                        edit('opt/linlsqOpt/input/TRTprescSetup.m');
+                        
+                    elseif TRT == 2
+                        run('opt/linlsqOpt/input/TRTprescSetup.m');
+                        
+                    elseif TRT == 3
+                        close(figure);
+                        
+                    end
+                    
+            elseif prescFiles == 3
+                close(figure);
+                
+            end
+            
+    elseif modules == 5
+        close(figure);
+                                
+    end

geometry/GeometryGenerator.m

@@ -0,0 +1,176 @@
+function Geometry = GeometryGenerator()
+
+%% % % GeometryGenerator % % % % %
+%  This function builds the geometry files needed for the WiscPlan dose 
+%  optimization.  Requires that you have my amiraStructReader function. You 
+%  need contours in any amira format. Contours should ALL BE CROPPED TO THE
+%  SAME DIMENSION.
+%
+%  DLB - 6/5/09
+
+s = pwd;
+
+% Get the planning CT
+[planningCT, CT_pathname, filterindex] = uigetfile({'*.am','Amira file (*.am)'},'Select the planning CT.','MultiSelect','off');
+cd(CT_pathname);
+
+% Generate a list of contours to include in Geometry.ROIS
+[files, pathname, filterindex] = uigetfile({'*.am','Amira file (*.am)'; '*.*','All files (*.*)'},'Select all the contours for the Geometry file (use Ctrl).','MultiSelect','on');
+cd(pathname);
+
+% Set the initial Geometry structure members
+data = amiraStructReader([CT_pathname planningCT]);
+Geometry.voxel_size = [data.xpix data.ypix data.zpix];
+Geometry.start = [data.boundingbox(1) data.boundingbox(3) data.boundingbox(5)];
+Geometry.data = single(data.image);
+
+%% Adding ROIS to Geometry
+for i = 1:length(files)
+    data = amiraStructReader([pathname files{i}]);
+    % remove the .am file extension
+    name = strrep(files{i},'.am','');
+    % remove everything before the first underscore
+    name = name(findstr(name,'_')+1:length(name));
+    % anything that remotely has the letters GTV is renamed GTV
+    if ~isempty(findstr(name,'GTV'))
+            name = 'GTV';
+    end
+    
+    % assign parameters
+    Geometry.ROIS{i}.name = name;
+    Geometry.ROIS{i}.ind = int32(find(data.image));
+    Geometry.ROIS{i}.dim = [data.dims(1) data.dims(2) data.dims(3)];
+    Geometry.ROIS{i}.pix_dim = [data.xpix data.ypix data.zpix];
+    Geometry.ROIS{i}.start = [data.boundingbox(1) data.boundingbox(3) data.boundingbox(5)];
+    Geometry.ROIS{i}.start_ind = '(1,1,1)';
+end
+
+
+cd(pathname);
+
+% save the Geometry file in the same location as your contours.
+save('Geometry.mat','Geometry');
+
+% return to the original directory
+cd(s);
+
+
+%% This is the amiraStructReader in case you don't have it. %%
+function data = amiraStructReader(varargin)
+
+%% % % amiraStructReader % % % % %
+%  Built for when you might need additional amira header info in matlab.
+%  Same functionality as amiraReader, but adds stuff in structure format 
+%
+%  DLB - 5/5/09  Cinco de Mayo!!!  But another 12 hour day :*(
+%
+
+warning off all;  % you'll get truncation warnings because you read inf in line 24.
+
+cd('C:\Data\');
+
+switch nargin 
+    case 1
+        filename = varargin{1};
+    case 0
+        [file, pathname] = uigetfile({'*.am'},'Please locate Amira file.','Multiselect','off');
+        filename = [pathname file];
+end
+
+fid = fopen(filename,'r','ieee-be');
+datachar = char(fread(fid,inf,'char'))'; % this will be the char data in the header
+switch isempty(findstr(datachar,'BINARY-LITTLE-ENDIAN'))
+    case 1
+        endianness = 'be';
+    case 0
+        endianness = 'le';
+end
+lattice = findstr(datachar,'define Lattice');
+parameter = findstr(datachar,'Parameters');
+start_data_position = findstr(datachar,'@1');  % note there will be two @1, so take the second one!
+
+fseek(fid,lattice+14,'bof');
+data.dims = str2num(fgetl(fid));
+
+bounding = findstr(datachar,'Bounding');
+fseek(fid,bounding + 11,'bof');
+data.boundingbox = str2num(fgetl(fid));
+
+
+data.x = (data.boundingbox(2)-data.boundingbox(1));
+data.y = (data.boundingbox(4)-data.boundingbox(3));
+data.z = (data.boundingbox(6)-data.boundingbox(5))*(1+1/(data.dims(3)-1));
+
+data.xpix = data.x/(data.dims(1)-1);
+data.ypix = data.y/(data.dims(2)-1);
+data.zpix = data.z/data.dims(3);
+
+%% Determine if we're dealing with a labelfield or regular data
+switch isempty(findstr(datachar,'HxByteRLE'))
+    case 1  % regular data case
+        data_type = char(datachar(findstr(datachar,'Lattice { ')+10:findstr(datachar,'Data } ')-2));
+        switch data_type
+            case 'byte'
+                out_type = 'int8';
+            case 'short'
+                out_type = 'int16';
+            case 'ushort'
+                out_type = 'uint16'; 
+            case 'int32'
+                out_type = 'int32';
+            case 'float'
+                out_type = 'float';
+            case 'double' 
+                out_type = 'double';
+        end
+        fid = fopen(filename,'r',['ieee-' endianness]);
+        fseek(fid,start_data_position(2)+2,'bof');
+        data.image = reshape(fread(fid,data.dims(1)*data.dims(2)*data.dims(3),out_type),data.dims(1),data.dims(2),data.dims(3));
+ 
+    case 0 % labelfield case
+        ByteRLE = findstr(datachar,'HxByteRLE,');
+        DataSection = findstr(datachar,'# Data section follows');
+        
+        fseek(fid,ByteRLE+9,'bof');
+        byte_length = str2num(char(fread(fid,DataSection - ByteRLE-13,'char'))');
+
+        switch endianness
+            case 'le'
+                fclose(fid);
+                fid = fopen(filename,'r',['ieee-' endianness]);
+                fseek(fid,start_data_position(2)+2,'bof');
+                data.image = fread(fid,byte_length,'int8');
+            case 'be'
+                fseek(fid,start_data_position(2)+2,'bof');
+                data.image = fread(fid,byte_length,'int8');
+        end
+
+        out_line = zeros(data.dims(1)*data.dims(2)*data.dims(3),1);
+
+        final_offset = sum(data.image(find(data.image<0))+127);
+
+        net_counter = 0;  % this sifts through the out_line variable
+        offset = 0;  % this adds an offset to data_counter when writing VERBATIM
+
+        for i = 1:floor((byte_length-final_offset)/2) % 
+            data_counter = 2*i+offset;  % this will sift through the data variable
+            if ((data.image(data_counter) == 0) && (data.image(data_counter-1) > 0)) % check if the value is 0 AND that we're not writing VERBATIM
+                net_counter = data.image(data_counter-1)+net_counter;
+            else
+                if data.image(data_counter-1) < 0 % check if we're writing VERBATIM
+                    for j = 1:128+data.image(data_counter-1)
+                        out_line(j+net_counter) = data.image(data_counter+j-1);
+                    end
+                    net_counter = 128+data.image(data_counter-1)+net_counter;
+                    offset = offset + 127+ data.image(data_counter-1);
+                else  % assume that we're copying value data(data_counter), data(data_counter-1) times.
+                    for j = 1:data.image(data_counter-1)
+                        out_line(j+net_counter) = data.image(data_counter);
+                    end
+                    net_counter = data.image(data_counter-1)+net_counter;
+                end
+            end
+        end 
+
+        data.image = reshape(out_line,data.dims(1),data.dims(2),data.dims(3));
+end

geometry/matlabPatientResample.asv

@@ -0,0 +1,149 @@
+% Stephen R Bowen April 2009
+% 
+% Crops and resamples the planning CT data set and associated ROI masks,
+% typically imported from Pinnacle using pinn2matlab.m.  CT data is
+% linearly resampled and ROI masks are nearest-neighbor resampled so as to
+% preserve binary values.  Bounding box coordinates are specified in the
+% crop ranges user inputs [xmin xmax], [ymin ymax], [zmin zmax].  Be aware
+% of half-voxel shifts in the boundary box used by Amira and remain
+% consistent in defining your crop ranges.
+
+%%%%% Start of user supplied inputs %%%%%
+
+% 
+infile = '../plan/HN003/HN003.mat';
+outfile = '../plan/HN003/HN003final.mat';
+
+% downsamples patient Geometry files, including their tissue masks
+
+% crop ranges (if empty, no cropping is done)
+xCrop = [-10.695 12.289];
+yCrop = [-66.202 -39.682];
+zCrop = [-6.801 8.669];
+
+% resolution after resampling (if empty, resolution is unchanged
+dxNew = 0.442;
+dyNew = 0.442;
+dzNew = 0.442;
+
+x0New = -10.695;
+y0New = -66.202;
+z0New = -6.801;
+
+%%%%% End of user supplied inputs %%%%%
+
+load(infile);
+
+% get dimensions of original data
+[M,N,Q] = size(Geometry.data);
+
+dx = Geometry.voxel_size(1);
+dy = Geometry.voxel_size(2);
+dz = Geometry.voxel_size(3);
+
+x0 = Geometry.start(1);
+y0 = Geometry.start(2);
+z0 = Geometry.start(3);
+
+% x,y,z coordinates of original data
+x = single([0:M-1]*dx + x0);
+y = single([0:N-1]*dy + y0);
+z = single([0:Q-1]*dz + z0);
+
+% If no cropping is desired, the resampling will be performed such that the
+% bounding planes of the 3D grid data are unchanged.
+
+% bounding plane coordinates:
+xp1 = x(1) - dx/2;
+xpM = x(end) + dx/2;
+yp1 = y(1) - dy/2;
+ypN = y(end) + dy/2;
+zp1 = z(1) - dz/2;
+zpQ = z(end) + dz/2;
+
+% create the new coordinate system
+if ~isempty(xCrop)
+    xp1 = xCrop(1);
+    xpM = xCrop(2);
+end
+
+if ~isempty(yCrop)
+    yp1 = yCrop(1);
+    ypN = yCrop(2);
+end
+
+if ~isempty(zCrop)
+    zp1 = zCrop(1);
+    zpQ = zCrop(2);
+end
+
+if isempty(dxNew)
+    dxNew = dx;
+end
+
+if isempty(dyNew)
+    dyNew = dy;
+end
+
+if isempty(dzNew)
+    dzNew = dz;
+end
+
+% create new coordinate systems
+if isempty(x0New)
+    xNew = single(xp1+dxNew/2:dxNew:xpM-dxNew/2);
+else
+    xNew = x0New:dxNew:xpM;
+    xNew = xNew(xNew >= xp1 & xNew <= xpM);
+end
+
+if isempty(y0New)
+    yNew = single(yp1+dyNew/2:dyNew:ypN-dyNew/2);
+else
+    yNew = y0New:dyNew:ypN;
+    yNew = yNew(yNew >= yp1 & yNew <= ypN);
+end
+
+if isempty(z0New)
+    zNew = single(zp1+dzNew/2:dzNew:zpQ-dzNew/2);
+else
+    zNew = z0New:dzNew:zpQ;
+    zNew = zNew(zNew >= zp1 & zNew <= zpQ);
+end
+
+clear Geometry_ds;
+Geometry_ds.voxel_size = [dxNew dyNew dzNew];
+Geometry_ds.start = [xNew(1) yNew(1) zNew(1)];
+
+fprintf('Resampling CT data:\n');
+Geometry_ds.data = gridResample3D(x,y,z,Geometry.data,xNew,yNew,zNew,'linear');
+fprintf('Finished resampling CT data.\n');
+
+Geometry_ds.ROIS = Geometry.ROIS;
+
+% resample the ROIs
+for k=1:length(Geometry.ROIS)    
+    % expand the original ROI mask
+    roimask = zeros(M,N,Q,'int8');  % create the original mask
+    roimask(Geometry.ROIS{k}.ind) = 1;
+
+    % resample the ROI mask
+    fprintf('Resampling ROI mask: %s',Geometry.ROIS{k}.name);
+    roimaskNew = gridResample3D(x,y,z,roimask,xNew,yNew,zNew,'nearest');
+    fprintf('  Finished resampling %s.\n',Geometry.ROIS{k}.name);
+    
+    % save the new roi mask data
+    Geometry_ds.ROIS{k}.ind = int32(find(roimaskNew));
+    Geometry_ds.ROIS{k}.dim = size(Geometry_ds.data);
+    Geometry_ds.ROIS{k}.pixdim = [dxNew dyNew dzNew];
+    Geometry_ds.ROIS{k}.start = [xNew(1) yNew(1) zNew(1)];
+    clear roimask roimaskNew;
+end
+
+if isfield(Geometry,'Trial');
+    Geometry_ds.Trial = Geometry.Trial;
+end
+
+Geometry = Geometry_ds;
+
+save(outfile,'Geometry');

geometry/matlabPatientResample.m

@@ -0,0 +1,156 @@
+%% matlabPatientResample.m
+
+% Stephen R Bowen April 2009
+% 
+% Crops and resamples the planning CT data set and associated ROI masks,
+% typically imported from Pinnacle using pinn2matlab.m.  CT data is
+% linearly resampled and ROI masks are nearest-neighbor resampled so as to
+% preserve binary values.  Bounding box coordinates are specified in the
+% crop ranges user inputs [xmin xmax], [ymin ymax], [zmin zmax].  Be aware
+% of half-voxel shifts in the boundary box used by Amira defined by the
+% centers of the 8 corner voxels of the image volume and remain
+% consistent in defining your crop ranges.
+
+%%
+
+%%%%% Start of user supplied inputs %%%%%
+
+% location of original Geometry structure and resampled/cropped Geometry
+% structure
+infile = '../plan/HN003/HN003.mat';
+outfile = '../plan/HN003/HN003final.mat';
+
+% crop ranges (if empty, no cropping is done)
+xCrop = [-10.695 12.289];
+yCrop = [-66.202 -39.682];
+zCrop = [-6.801 8.669];
+
+% resolution after resampling (if empty, resolution is unchanged)
+dxNew = 0.442;
+dyNew = 0.442;
+dzNew = 0.442;
+
+% new start coordinates if necessary
+x0New = -10.695;
+y0New = -66.202;
+z0New = -6.801;
+
+%%%%% End of user supplied inputs %%%%%
+
+%%
+
+load(infile);
+
+% get dimensions of original data
+[M,N,Q] = size(Geometry.data);
+
+dx = Geometry.voxel_size(1);
+dy = Geometry.voxel_size(2);
+dz = Geometry.voxel_size(3);
+
+x0 = Geometry.start(1);
+y0 = Geometry.start(2);
+z0 = Geometry.start(3);
+
+% x,y,z coordinates of original data
+x = single([0:M-1]*dx + x0);
+y = single([0:N-1]*dy + y0);
+z = single([0:Q-1]*dz + z0);
+
+% If no cropping is desired, the resampling will be performed such that the
+% bounding planes of the 3D grid data are unchanged.
+
+% bounding plane coordinates:
+xp1 = x(1) - dx/2;
+xpM = x(end) + dx/2;
+yp1 = y(1) - dy/2;
+ypN = y(end) + dy/2;
+zp1 = z(1) - dz/2;
+zpQ = z(end) + dz/2;
+
+% create the new coordinate system
+if ~isempty(xCrop)
+    xp1 = xCrop(1);
+    xpM = xCrop(2);
+end
+
+if ~isempty(yCrop)
+    yp1 = yCrop(1);
+    ypN = yCrop(2);
+end
+
+if ~isempty(zCrop)
+    zp1 = zCrop(1);
+    zpQ = zCrop(2);
+end
+
+if isempty(dxNew)
+    dxNew = dx;
+end
+
+if isempty(dyNew)
+    dyNew = dy;
+end
+
+if isempty(dzNew)
+    dzNew = dz;
+end
+
+% create new coordinate systems
+if isempty(x0New)
+    xNew = single(xp1+dxNew/2:dxNew:xpM-dxNew/2);
+else
+    xNew = x0New:dxNew:xpM;
+    xNew = xNew(xNew >= xp1 & xNew <= xpM);
+end
+
+if isempty(y0New)
+    yNew = single(yp1+dyNew/2:dyNew:ypN-dyNew/2);
+else
+    yNew = y0New:dyNew:ypN;
+    yNew = yNew(yNew >= yp1 & yNew <= ypN);
+end
+
+if isempty(z0New)
+    zNew = single(zp1+dzNew/2:dzNew:zpQ-dzNew/2);
+else
+    zNew = z0New:dzNew:zpQ;
+    zNew = zNew(zNew >= zp1 & zNew <= zpQ);
+end
+
+clear Geometry_ds;
+Geometry_ds.voxel_size = [dxNew dyNew dzNew];
+Geometry_ds.start = [xNew(1) yNew(1) zNew(1)];
+
+fprintf('Resampling CT data:\n');
+Geometry_ds.data = gridResample3D(x,y,z,Geometry.data,xNew,yNew,zNew,'linear');
+fprintf('Finished resampling CT data.\n');
+
+Geometry_ds.ROIS = Geometry.ROIS;
+
+% resample the ROIs
+for k=1:length(Geometry.ROIS)    
+    % expand the original ROI mask
+    roimask = zeros(M,N,Q,'int8');  % create the original mask
+    roimask(Geometry.ROIS{k}.ind) = 1;
+
+    % resample the ROI mask
+    fprintf('Resampling ROI mask: %s',Geometry.ROIS{k}.name);
+    roimaskNew = gridResample3D(x,y,z,roimask,xNew,yNew,zNew,'nearest');
+    fprintf('  Finished resampling %s.\n',Geometry.ROIS{k}.name);
+    
+    % save the new roi mask data
+    Geometry_ds.ROIS{k}.ind = int32(find(roimaskNew));
+    Geometry_ds.ROIS{k}.dim = size(Geometry_ds.data);
+    Geometry_ds.ROIS{k}.pixdim = [dxNew dyNew dzNew];
+    Geometry_ds.ROIS{k}.start = [xNew(1) yNew(1) zNew(1)];
+    clear roimask roimaskNew;
+end
+
+if isfield(Geometry,'Trial');
+    Geometry_ds.Trial = Geometry.Trial;
+end
+
+Geometry = Geometry_ds;
+
+save(outfile,'Geometry');

geometry/pinn2matlab.asv

@@ -0,0 +1,492 @@
+% pinn2matlab.m
+%
+%  Ryan T Flynn 23 August 2007
+%
+% Reads in raw Pinnacle data to a MATLAB data structure. The patient CT,
+% ROIs, and Trial files can all be read in if available and desired. The
+% result is a structure called 'Geometry', which contains ROIS and Trial
+% fields.
+%
+% The CT image and all of the ROIs are flipped onto a coordinate system in
+% which the position of each voxel index (i,j,k) are given by:
+%
+% x(i) = x0 + i*delx; (i=0,...,M-1)
+% y(j) = y0 + j*delx; (j=0,...,N-1)
+% z(k) = z0 + k*delx; (k=0,...,Q-1)
+% 
+% As opposed to the Pinnacle coordinate system, which
+% 
+% x(i) = x0 + i*delx;        (i=0,...,M-1)
+% y(j) = y0 + (N-1-j)*dely;  (j=0,...,N-1)
+% z(k) = z0 + k*delz;        (k=0,...,Q-1)
+% 
+% Stephen R Bowen April 2009
+% 
+% The MATLAB data structure is called Geometry with the following format:
+% 
+% Geometry.voxel_size  (dx dy dz)
+%         .start       (x y z)
+%         .data        (M x N x Q single float matrix)
+%         .ROIS        {1 x Nroi cell array} if readROIS flag on
+% 
+% Within each index of the ROIS cell array, a data structure exists
+% containing the vector information from a Pinnacle mesh contour and indices
+% of the ROI mask used by WiscPlan:
+% 
+% Geometry.ROIS{1}.name         'string'
+%                 .num_curves   int
+%                 .curves       {1 x num_curves cell array}
+%                 .ind          (int32 data array)
+%                 .dim          (M N Q)
+%                 .pix_dim      (dx dy dz)
+%                 .start        (x y z)
+%                 .start_ind    'string' 
+%
+% This enables seamless translation between Pinnacle and MATLAB
+
+%%%%% Start of user supplied inputs %%%%%
+
+% Read-in flags
+readTrial = 'no';  % either 'yes' or 'no'
+readROIS = 'yes';   % either 'yes' or 'no'
+
+% locations of the Pinnacle files
+roifile = '../plan/HN003/plan.roi';
+imageheaderfile = '../plan/HN003/HN003_planningCT.header';
+imagefile = '../plan/HN003/HN003_planningCT.img';  % CT file
+save_file = '../plan/HN003/HN003.mat';
+
+%%%%% End of user supplied inputs %%%%%
+
+start_ind = '(1,1,1)';  % this tag will be added to the ROIS so it is known where the start voxel is
+
+clear roi ROIS Trial Geometry CTimage;  % start with a clean slate of structures
+
+% extract geometric information from the header file
+fid_imageheaderfile = fopen(imageheaderfile,'r');
+
+tline = fgets(fid_imageheaderfile);
+while tline ~= -1
+    % check the line for key words
+    if length(findstr(tline,'byte_order'))
+        eval(tline);
+    if length(findstr(tline,'x_dim')) & ~length(findstr(tline,'fname_index_start'))
+        eval(tline);  % run the line to get x_dim
+    elseif length(findstr(tline,'y_dim'))
+        eval(tline);  % run the line to get y_dim
+    elseif length(findstr(tline,'z_dim'))
+        eval(tline);  % run the line to get z_dim
+    elseif length(findstr(tline,'x_pixdim'))
+        eval(tline);  % run the line to get x_pixdim
+    elseif length(findstr(tline,'y_pixdim'))
+        eval(tline);  % run the line to get y_pixdim
+    elseif length(findstr(tline,'z_pixdim'))
+        eval(tline);  % run the line to get z_pixdim
+    elseif length(findstr(tline,'x_start'))
+        eval(tline);  % run the line to get x_start
+    elseif length(findstr(tline,'y_start'))
+        eval(tline);  % run the line to get x_start
+    elseif length(findstr(tline,'z_start'))
+        eval(tline);  % run the line to get x_start
+    end
+    tline = fgets(fid_imageheaderfile);
+end
+
+fclose(fid_imageheaderfile);
+
+% Read in the CT data
+fid_image = fopen(imagefile,'rb');
+if byte_order == 0
+    CTimage = fread(fid_image,'short=>int16');
+elseif byte_order == 1
+    CTimage = fread(fid_image,'short=>int16','ieee-be');
+end
+CTimage = reshape(CTimage,x_dim,y_dim,z_dim);
+fclose(fid_image);
+
+% flip the CT image onto the non-Pinnacle coordinate system
+CTimage = flipdim(CTimage,2);
+
+% save the results to the Geometry structure
+Geometry.start = single([x_start y_start z_start]);
+Geometry.voxel_size = single([x_pixdim y_pixdim z_pixdim]);
+Geometry.data = single(CTimage)/1024; % convert the CT data to density data, Pinnacle style
+Geometry.data(Geometry.data < 0) = 0; % truncate anything less than zero
+
+clear CTimage;
+
+if strcmp(readROIS,'yes')
+    % read in the roi file
+    fid_roifile = fopen(roifile,'r');
+
+    roinames = {};  % start a cell array of the roi names
+    Nroi = 0;  % number of rois read in
+
+    % Flags to indicate which sets of angled brackets in the roi file tline is
+    % inside.
+    inroi = 0;
+    incurve = 0;
+    inpoints = 0;
+
+    roi_num = 0;  % current roi number
+
+    tline = fgets(fid_roifile);
+    while tline ~= -1
+        % check the line for key words
+        if length(findstr(tline,'roi={'))
+            inroi = 1;   % mark that we are now currently inside of an roi
+            roi_num = roi_num + 1;
+            % next line contains the roi name
+            tline = fgets(fid_roifile);
+            % pop off first token in line, the remainder of the line is the roi name
+            [T,R] = strtok(tline);
+            roi.name = strtrim(R);
+            % pop off lines until we get to the number of curves in this roi
+            while ~length(findstr(tline,'num_curve'))
+                tline = fgets(fid_roifile);
+            end
+            % pop off the num_curve string
+            [T,R] = strtok(tline);
+            % pop off the equals sign
+            [T,R] = strtok(R);
+            % pop off the number of curves in this roi
+            T = strtok(R,';');
+            % save the number of curves to the roi stucture
+            eval(['roi.num_curves = ' num2str(T) ';']);
+            roi.curves = {};  % get the curves structure started
+            % read in the next curve structure
+            curve_num = 0;  % number of the current curve
+            while roi.num_curves > 0 & curve_num < roi.num_curves
+                while ~length(findstr(tline,'curve={'));
+                    tline = fgets(fid_roifile);
+                end
+                curve_num = curve_num + 1;
+                incurve = 1;  % inside the curve structure now
+                % find the number of points in this structure
+                while ~length(findstr(tline,'num_points'))
+                    tline = fgets(fid_roifile);
+                end
+                % pop off the num_points string
+                [T,R] = strtok(tline);
+                % pop off the equals sign
+                [T,R] = strtok(R);
+                % pop off the number of points in this curve
+                T = strtok(R,';');
+                eval(['num_points = ' num2str(T) ';']);
+                % find the points identifier
+                while ~length(findstr(tline,'points={'))
+                    tline = fgets(fid_roifile);
+                end
+                inpoints = 1;  % inside the points structure now
+                % read in the block of points data
+                block = fscanf(fid_roifile,'%g',[3 num_points]);
+                % save the block of points to the roi stucture
+                roi.curves{curve_num} = block';
+                % read in the right parantheses for the points and curve
+                % structures
+                while ~length(findstr(tline,'};'))
+                    tline = fgets(fid_roifile);
+                end
+                inpoints = 0;   % stepped outside of the points structure
+                while ~length(findstr(tline,'};'))
+                    tline = fgets(fid_roifile);
+                end
+                incurve = 0;  % stepped outside of the curves structure
+            end
+            % read until the roi closing bracket appears
+            while ~length(findstr(tline,'};'))
+                tline = fgets(fid_roifile);
+            end
+            inroi = 0;  % we are now outside of the roi
+            fprintf('read in %s roi\n',roi.name);
+            ROIS{roi_num} = roi;
+        end
+        tline = fgets(fid_roifile);
+    end
+
+    fclose(fid_roifile);
+
+    % ROIS have all been read-in, now just have to convert them to mask
+    % matrices.  We'll use roipoly for this.
+    % Recall that we must use the Pinnacle coordinate system for this to work,
+    % that is, (x,y,z) coordinates are given by:
+    % x(i) = x_start + i*x_pixdim, i=0..x_dim-1
+    % y(j) = y_start + (y_dim-1)*y_pixdim - j*y_pixdim, j=0..y_dim-1
+    % z(k) = z_start + k*z_pixdim, k=0..z_dim-1
+    %
+    % Not all treatment planning systems use this type of coordinate system
+    % definition, so it is very important to get them straight.
+    %
+    % To get around this we will manipulate the data such that we'll have:
+    %
+    % x(i) = x_start + i*x_pixdim, i=0..x_dim-1
+    % y(j) = y_start + j*y_pixdim, j=0..y_dim-1
+    % z(k) = z_start + k*z_pixdim, k=0..z_dim-1
+    %
+    % This can be done by a simple fliplr operation on all of the CT slices
+
+    % define the coordinate system
+    x = x_start:x_pixdim:x_start + (x_dim-1)*x_pixdim;
+    y = y_start:y_pixdim:y_start + (y_dim-1)*y_pixdim;
+    z = z_start:z_pixdim:z_start + (z_dim-1)*z_pixdim;
+
+    % define the locations of the corners of the pixels in each slice
+    xCorners = [x - x_pixdim/2 x(end) + x_pixdim/2];
+    yCorners = [y - y_pixdim/2 y(end) + y_pixdim/2];
+    
+    % loop through all rois
+    for roi_num=1:length(ROIS)
+        % set up the roi mask
+        roimask = zeros(x_dim,y_dim,z_dim,'int8');
+        roimaskCorners = zeros(x_dim+1,y_dim+1,z_dim,'int8');
+        % loop through all of the curves in the roi
+        for curve_num=1:length(ROIS{roi_num}.curves)
+            % make a copy of the curve for easy access
+            current_curve = ROIS{roi_num}.curves{curve_num};
+            % find the z-index (slice number) of this curve
+            q = round((current_curve(1,3)-z_start)/z_pixdim) + 1;
+            
+            % put these index vectors into roipoly
+            if q >= 1 & q <= z_dim
+                roisliceCorners = double(roimaskCorners(:,:,q));
+                % find which corners are inside the contour
+                BWcorners = roipoly(yCorners,xCorners,roisliceCorners,current_curve(:,2),current_curve(:,1)); 
+                % Mark all all pixels bordering corners that are inside the
+                % contour
+                roi_vox = sum(BWcorners(:));  % number of voxels in this ROI
+                % find the voxel overlap between the current roi and BW:
+                overlap_vox = sum(sum(BWcorners.*roisliceCorners));
+                if overlap_vox == roi_vox
+                    roisliceCorners = roisliceCorners - BWcorners;
+                else  % if there is not perfect overlap, add the rois
+                    roisliceCorners = roisliceCorners + BWcorners;
+                end
+                roisliceCorners(roisliceCorners > 0) = 1;  % make sure all mask values are unity
+                roimaskCorners(:,:,q) = int8(roisliceCorners);  % update the overall mask
+            end
+        end
+        
+        % save time be resampling only a subregion
+        ind = find(roimaskCorners);
+        [I,J,K] = ind2sub([x_dim+1 y_dim+1 z_dim],ind);
+        indx = min(I)-3:max(I)+3;
+        indy = min(J)-3:max(J)+3;
+        indz = min(K)-3:max(K)+3;
+        
+        indx = indx(indx >= 1 & indx <= x_dim);
+        indy = indy(indy >= 1 & indy <= y_dim);
+        indz = indz(indz >= 1 & indz <= z_dim);
+        
+        % convert the corners to a 3-D roi mask
+        roimask(indx,indy,indz) = ceil(gridResample3D(xCorners,yCorners,z,roimaskCorners,x(indx),y(indy),z(indz)));
+        
+        % save the indices of the roi mask
+        ROIS{roi_num}.ind = int32(find(roimask ~= 0));
+        ROIS{roi_num}.dim = [x_dim y_dim z_dim];
+        ROIS{roi_num}.pixdim = [x_pixdim y_pixdim z_pixdim];
+        ROIS{roi_num}.start = [x_start y_start z_start];
+        ROIS{roi_num}.start_ind = start_ind;
+        fprintf('Converted %s vectors to an roi mask.\n',ROIS{roi_num}.name);
+    end
+    
+    % save ROIS to the Geometry structure
+    Geometry.ROIS = ROIS;
+    clear ROIS;
+end
+
+if strcmp(readTrial,'yes')
+    % read in the Trial information
+    fid_trialfile = fopen(trialfile,'r');
+    tline = fgets(fid_trialfile);
+
+    structstack = {};   % stack for determining which field we're inside
+    structind = [];    % structure index values for each field
+    lines = 1;
+
+    while tline ~= -1   % read through the trial file
+        tline = regexprep(tline,'"','''');  % replace " with '
+        tline = regexprep(tline,'[','');  % bad characters to remove
+        tline = regexprep(tline,']','');  % bad characters to remove
+        tline = regexprep(tline,'#','num');
+        tline = regexprep(tline,'\\','''');
+        % fprintf(tline);
+        if length(findstr(tline,'{'))
+            % mark the structure we just entered, removing any blank space
+            structstack{length(structstack) + 1} = regexprep(strtok(tline,'='),' ','');
+            if length(structstack) == 1
+                structstack{1} = 'Trial';  % make sure structure name is always the same
+            end
+            if strcmp(structstack{end},'Beam')  % found a beam structure
+                % find the number of existing beams at this location
+                eval_line = ['beams = getfield(' structstack{1}];
+                for k=2:length(structstack)
+                    eval_line = [eval_line ',''' structstack{k} ''''];
+                end
+                eval_line = [eval_line ');'];
+                try
+                    eval(eval_line);
+                    num_beams = length(beams) + 1;  % this is the current beam number
+                catch  % there was an error, so there are no beams yet
+                    num_beams = 1;  % we're on the first beam now
+                end
+            end
+        elseif length(findstr(tline,'}'))
+            structstack(end) = [];  % step outside of the current field
+        else  % the line must be a subfield assignment
+            if strcmp(structstack{end},'Points')  % this case treated specially
+                % read-in values until the next '}' is reached
+                lines = lines + 1;
+                mlc_vec = ['['];
+                while ~length(findstr(tline,'}'))
+                    % build up the vector of MLC positions
+                    mlc_vec = [mlc_vec regexprep(tline,' ','')];
+                    tline = fgets(fid_trialfile);
+                    % fprintf(tline);
+                    lines = lines + 1;
+                end
+                mlc_vec = [mlc_vec ']'];
+                % evaluate the points vector
+                eval_statement = structstack{1};
+                for k=2:length(structstack)
+                    if strcmp(structstack{k},'Beam')
+                        % include the beam number
+                        eval_statement = [eval_statement '.' structstack{k} '(' num2str(num_beams) ')'];
+                    else  % non-beam structure, so don't worry about it
+                        eval_statement = [eval_statement '.' structstack{k}];
+                    end
+                end
+                eval_statement = [eval_statement ' = ' mlc_vec ';'];
+                eval(eval_statement);
+                structstack(end) = [];  % step outside of the current Points[] field
+            else
+                % see if there are any subfields in the expression
+                % evaluate the statement
+                eval_statement = structstack{1};
+                for k=2:length(structstack)
+                    if strcmp(structstack{k},'Beam')
+                        % include the beam number
+                        eval_statement = [eval_statement '.' structstack{k} '(' num2str(num_beams) ')'];
+                    else  % non-beam structure, so don't worry about it
+                        eval_statement = [eval_statement '.' structstack{k}];
+                    end
+                end
+                eval_statement = [eval_statement '.' tline];
+                eval(eval_statement);
+            end
+        end
+        tline = fgets(fid_trialfile);
+        lines = lines + 1;
+        if isempty(structstack)
+            break;   % we're done if we're no longer inside any sets of curly brackets
+        end
+    end
+    fclose(fid_trialfile);
+
+    % dose grid dimensions
+    xdim = Trial.DoseGrid.Dimension.X;
+    ydim = Trial.DoseGrid.Dimension.Y;
+    zdim = Trial.DoseGrid.Dimension.Z;
+    siz = [xdim ydim zdim];  % dosegrid size vector
+
+    % vector describing the start location of the dose grids
+    xstart = Trial.DoseGrid.Origin.X;
+    ystart = Trial.DoseGrid.Origin.Y;
+    zstart = Trial.DoseGrid.Origin.Z;
+    start = [xstart ystart zstart];  % dosegrid start vector
+
+    % vector describing the size of each voxel
+    dx = Trial.DoseGrid.VoxelSize.X;
+    dy = Trial.DoseGrid.VoxelSize.Y;
+    dz = Trial.DoseGrid.VoxelSize.Z;
+    voxel_size = [dx dy dz];
+
+    % Read in the beamlet dose distributions, DRRs, and fluences and store them in
+    % Trial.BeamList.Beam(k), where k is the beam number
+    for k=0:num_beams-1
+        % read-in the dose distribution
+        num_vec = num2str(k*5+4);
+        % ensure that num_vec has a length of 3
+        while length(num_vec) < 3
+            num_vec = ['0' num_vec];
+        end
+        dose_file = [dose_file_base_name '.' num_vec];
+        fid = fopen(dose_file,'rb','ieee-be');
+        if fid ~= -1
+            dose = fread(fid,'float=>float');
+            fclose(fid);
+            try
+                dose = reshape(dose,siz);
+                % flip the dose distribution in accordance with the CT image
+                for j=1:siz(3)
+                    dose(:,:,j) = single(fliplr(double(dose(:,:,j))));
+                end
+                % save the dose
+                Trial = setfield(Trial,'BeamList','Beam',{k+1},'dose',dose);
+            catch
+                fprintf('Dose for beam %g did not meet specifications so it was not read-in.\n',k+1);
+            end
+        end
+
+        % read-in the DRR
+        % find the size of this drr
+        drr_x = Trial.BeamList.Beam(1).FilmImageList.FilmImage.Image.Dimension.X;
+        drr_y = Trial.BeamList.Beam(1).FilmImageList.FilmImage.Image.Dimension.Y;
+        siz_drr = [drr_x drr_y];
+
+        num_vec = num2str(k*5+3);
+        % ensure that num_vec has a length of 3
+        while length(num_vec) < 3
+            num_vec = ['0' num_vec];
+        end
+        drr_file = [dose_file_base_name '.' num_vec];
+        fid = fopen(drr_file,'rb','ieee-be');
+        if fid ~= -1
+            drr = fread(fid,'float');
+            fclose(fid);
+            try  % proceed only if the drr was read-in properly
+                drr = reshape(drr,siz_drr);
+                % save the DRR
+                Trial = setfield(Trial,'BeamList','Beam',{k+1},'FilmImageList','FilmImage','Image','drr',drr);
+            catch
+                fprintf('DRR for beam %g did not meet specifications so it was not read-in.\n',k+1);
+            end
+        end
+
+        % read in the fluence map
+        num_vec = num2str(k*5+2);
+        fluence_siz = [101 101 3];  % all fluence must have this size otherwise and error will be thrown
+        % ensure that num_vec has a length of 3
+        while length(num_vec) < 3
+            num_vec = ['0' num_vec];
+        end
+        fluence_file = [dose_file_base_name '.' num_vec];
+        fid = fopen(fluence_file,'rb','ieee-be');
+        if fid ~= -1
+            fluence = fread(fid,'float');
+            fclose(fid);
+            % check to see if fluence matches dose grid size
+            if Trial.FluenceGridMatchesDoseGrid ~= 1
+                error('Fluence grid does not match dose grid for beam %g, don''t know what to do here.',k+1);
+            end
+            try  % proceed only if the fluence is read-in properly
+                fluence = reshape(fluence,fluence_siz);
+                % save the fluence
+                Trial = setfield(Trial,'BeamList','Beam',{k+1},'FluenceMap',fluence);
+            catch
+                fprintf('Fluence map for beam %g did not meet specifications so it was not read-in.\n',k+1);
+            end
+        end
+    end
+
+    % add up all of the doses
+    dose_tot = zeros(siz);
+    for k=1:num_beams
+        if isfield(Trial.BeamList.Beam(k),'dose') & ~isempty(Trial.BeamList.Beam(k).dose)
+            dose_tot = dose_tot + double(Trial.BeamList.Beam(k).dose);
+        end
+    end
+    Geometry.Trial = Trial;
+    clear Trial;
+end
+
+save(save_file,'Geometry');

geometry/pinn2matlab.m

@@ -0,0 +1,497 @@
+%% pinn2matlab.m
+
+%  Ryan T Flynn 23 August 2007
+%
+% Reads in raw Pinnacle data to a MATLAB data structure. The patient CT,
+% ROIs, and Trial files can all be read in if available and desired. The
+% result is a structure called 'Geometry', which contains ROIS and Trial
+% fields.
+%
+% The CT image and all of the ROIs are flipped onto a coordinate system in
+% which the position of each voxel index (i,j,k) are given by:
+%
+% x(i) = x0 + i*delx; (i=0,...,M-1)
+% y(j) = y0 + j*delx; (j=0,...,N-1)
+% z(k) = z0 + k*delx; (k=0,...,Q-1)
+% 
+% As opposed to the Pinnacle coordinate system, which
+% 
+% x(i) = x0 + i*delx;        (i=0,...,M-1)
+% y(j) = y0 + (N-1-j)*dely;  (j=0,...,N-1)
+% z(k) = z0 + k*delz;        (k=0,...,Q-1)
+% 
+% Stephen R Bowen April 2009
+% 
+% The MATLAB structure Geometry has the following format:
+% 
+% Geometry.voxel_size  (dx dy dz)
+%         .start       (x y z)
+%         .data        (M x N x Q single float matrix)
+%         .ROIS        {1 x Nroi cell array} if readROIS flag on
+% 
+% Within each index of the ROIS cell array, a data structure exists
+% containing the vector information from each Pinnacle mesh contour and indices
+% of the corresponding ROI mask used by WiscPlan:
+% 
+% Geometry.ROIS{1}.name         'string'
+%                 .num_curves   int
+%                 .curves       {1 x num_curves cell array}
+%                 .ind          (int32 data array)
+%                 .dim          (M N Q)
+%                 .pix_dim      (dx dy dz)
+%                 .start        (x y z)
+%                 .start_ind    'string' 
+%
+% This enables seamless translation between Pinnacle and MATLAB.
+
+%%
+
+%%%%% Start of user supplied inputs %%%%%
+
+% Read-in flags
+readTrial = 'no';  % either 'yes' or 'no'
+readROIS = 'yes';   % either 'yes' or 'no'
+
+% locations of the Pinnacle files
+roifile = '../plan/HN003/plan.roi';
+imageheaderfile = '../plan/HN003/HN003_planningCT.header';
+imagefile = '../plan/HN003/HN003_planningCT.img';  % CT file
+save_file = '../plan/HN003/HN003.mat';
+
+%%%%% End of user supplied inputs %%%%%
+
+%%
+
+start_ind = '(1,1,1)';  % this tag will be added to the ROIS so it is known where the start voxel is
+
+clear roi ROIS Trial Geometry CTimage;  % start with a clean slate of structures
+
+% extract geometric information from the header file
+fid_imageheaderfile = fopen(imageheaderfile,'r');
+
+tline = fgets(fid_imageheaderfile);
+while tline ~= -1
+    % check the line for key words
+    if length(findstr(tline,'byte_order'))
+        eval(tline);
+    elseif length(findstr(tline,'x_dim')) & ~length(findstr(tline,'fname_index_start'))
+        eval(tline);  % run the line to get x_dim
+    elseif length(findstr(tline,'y_dim'))
+        eval(tline);  % run the line to get y_dim
+    elseif length(findstr(tline,'z_dim'))
+        eval(tline);  % run the line to get z_dim
+    elseif length(findstr(tline,'x_pixdim'))
+        eval(tline);  % run the line to get x_pixdim
+    elseif length(findstr(tline,'y_pixdim'))
+        eval(tline);  % run the line to get y_pixdim
+    elseif length(findstr(tline,'z_pixdim'))
+        eval(tline);  % run the line to get z_pixdim
+    elseif length(findstr(tline,'x_start'))
+        eval(tline);  % run the line to get x_start
+    elseif length(findstr(tline,'y_start'))
+        eval(tline);  % run the line to get x_start
+    elseif length(findstr(tline,'z_start'))
+        eval(tline);  % run the line to get x_start
+    end
+    tline = fgets(fid_imageheaderfile);
+end
+
+fclose(fid_imageheaderfile);
+
+% Read in the CT data
+fid_image = fopen(imagefile,'rb');
+if byte_order == 0
+    CTimage = fread(fid_image,'short=>int16');
+elseif byte_order == 1
+    CTimage = fread(fid_image,'short=>int16','ieee-be');
+end
+
+CTimage = reshape(CTimage,x_dim,y_dim,z_dim);
+fclose(fid_image);
+
+% flip the CT image onto the non-Pinnacle coordinate system
+CTimage = flipdim(CTimage,2);
+
+% save the results to the Geometry structure
+Geometry.start = single([x_start y_start z_start]);
+Geometry.voxel_size = single([x_pixdim y_pixdim z_pixdim]);
+Geometry.data = single(CTimage)/1024; % convert the CT data to density data, Pinnacle style
+Geometry.data(Geometry.data < 0) = 0; % truncate anything less than zero
+
+clear CTimage;
+
+if strcmp(readROIS,'yes')
+    % read in the roi file
+    fid_roifile = fopen(roifile,'r');
+
+    roinames = {};  % start a cell array of the roi names
+    Nroi = 0;  % number of rois read in
+
+    % Flags to indicate which sets of angled brackets in the roi file tline is
+    % inside.
+    inroi = 0;
+    incurve = 0;
+    inpoints = 0;
+
+    roi_num = 0;  % current roi number
+
+    tline = fgets(fid_roifile);
+    while tline ~= -1
+        % check the line for key words
+        if length(findstr(tline,'roi={'))
+            inroi = 1;   % mark that we are now currently inside of an roi
+            roi_num = roi_num + 1;
+            % next line contains the roi name
+            tline = fgets(fid_roifile);
+            % pop off first token in line, the remainder of the line is the roi name
+            [T,R] = strtok(tline);
+            roi.name = strtrim(R);
+            % pop off lines until we get to the number of curves in this roi
+            while ~length(findstr(tline,'num_curve'))
+                tline = fgets(fid_roifile);
+            end
+            % pop off the num_curve string
+            [T,R] = strtok(tline);
+            % pop off the equals sign
+            [T,R] = strtok(R);
+            % pop off the number of curves in this roi
+            T = strtok(R,';');
+            % save the number of curves to the roi stucture
+            eval(['roi.num_curves = ' num2str(T) ';']);
+            roi.curves = {};  % get the curves structure started
+            % read in the next curve structure
+            curve_num = 0;  % number of the current curve
+            while roi.num_curves > 0 & curve_num < roi.num_curves
+                while ~length(findstr(tline,'curve={'));
+                    tline = fgets(fid_roifile);
+                end
+                curve_num = curve_num + 1;
+                incurve = 1;  % inside the curve structure now
+                % find the number of points in this structure
+                while ~length(findstr(tline,'num_points'))
+                    tline = fgets(fid_roifile);
+                end
+                % pop off the num_points string
+                [T,R] = strtok(tline);
+                % pop off the equals sign
+                [T,R] = strtok(R);
+                % pop off the number of points in this curve
+                T = strtok(R,';');
+                eval(['num_points = ' num2str(T) ';']);
+                % find the points identifier
+                while ~length(findstr(tline,'points={'))
+                    tline = fgets(fid_roifile);
+                end
+                inpoints = 1;  % inside the points structure now
+                % read in the block of points data
+                block = fscanf(fid_roifile,'%g',[3 num_points]);
+                % save the block of points to the roi stucture
+                roi.curves{curve_num} = block';
+                % read in the right parantheses for the points and curve
+                % structures
+                while ~length(findstr(tline,'};'))
+                    tline = fgets(fid_roifile);
+                end
+                inpoints = 0;   % stepped outside of the points structure
+                while ~length(findstr(tline,'};'))
+                    tline = fgets(fid_roifile);
+                end
+                incurve = 0;  % stepped outside of the curves structure
+            end
+            % read until the roi closing bracket appears
+            while ~length(findstr(tline,'};'))
+                tline = fgets(fid_roifile);
+            end
+            inroi = 0;  % we are now outside of the roi
+            fprintf('read in %s roi\n',roi.name);
+            ROIS{roi_num} = roi;
+        end
+        tline = fgets(fid_roifile);
+    end
+
+    fclose(fid_roifile);
+
+    % ROIS have all been read-in, now just have to convert them to mask
+    % matrices.  We'll use roipoly for this.
+    % Recall that we must use the Pinnacle coordinate system for this to work,
+    % that is, (x,y,z) coordinates are given by:
+    % x(i) = x_start + i*x_pixdim, i=0..x_dim-1
+    % y(j) = y_start + (y_dim-1)*y_pixdim - j*y_pixdim, j=0..y_dim-1
+    % z(k) = z_start + k*z_pixdim, k=0..z_dim-1
+    %
+    % Not all treatment planning systems use this type of coordinate system
+    % definition, so it is very important to get them straight.
+    %
+    % To get around this we will manipulate the data such that we'll have:
+    %
+    % x(i) = x_start + i*x_pixdim, i=0..x_dim-1
+    % y(j) = y_start + j*y_pixdim, j=0..y_dim-1
+    % z(k) = z_start + k*z_pixdim, k=0..z_dim-1
+    %
+    % This can be done by a simple fliplr operation on all of the CT slices
+
+    % define the coordinate system
+    x = x_start:x_pixdim:x_start + (x_dim-1)*x_pixdim;
+    y = y_start:y_pixdim:y_start + (y_dim-1)*y_pixdim;
+    z = z_start:z_pixdim:z_start + (z_dim-1)*z_pixdim;
+
+    % define the locations of the corners of the pixels in each slice
+    xCorners = [x - x_pixdim/2 x(end) + x_pixdim/2];
+    yCorners = [y - y_pixdim/2 y(end) + y_pixdim/2];
+    
+    % loop through all rois
+    for roi_num=1:length(ROIS)
+        % set up the roi mask
+        roimask = zeros(x_dim,y_dim,z_dim,'int8');
+        roimaskCorners = zeros(x_dim+1,y_dim+1,z_dim,'int8');
+        % loop through all of the curves in the roi
+        for curve_num=1:length(ROIS{roi_num}.curves)
+            % make a copy of the curve for easy access
+            current_curve = ROIS{roi_num}.curves{curve_num};
+            % find the z-index (slice number) of this curve
+            q = round((current_curve(1,3)-z_start)/z_pixdim) + 1;
+            
+            % put these index vectors into roipoly
+            if q >= 1 & q <= z_dim
+                roisliceCorners = double(roimaskCorners(:,:,q));
+                % find which corners are inside the contour
+                BWcorners = roipoly(yCorners,xCorners,roisliceCorners,current_curve(:,2),current_curve(:,1)); 
+                % Mark all all pixels bordering corners that are inside the
+                % contour
+                roi_vox = sum(BWcorners(:));  % number of voxels in this ROI
+                % find the voxel overlap between the current roi and BW:
+                overlap_vox = sum(sum(BWcorners.*roisliceCorners));
+                if overlap_vox == roi_vox
+                    roisliceCorners = roisliceCorners - BWcorners;
+                else  % if there is not perfect overlap, add the rois
+                    roisliceCorners = roisliceCorners + BWcorners;
+                end
+                roisliceCorners(roisliceCorners > 0) = 1;  % make sure all mask values are unity
+                roimaskCorners(:,:,q) = int8(roisliceCorners);  % update the overall mask
+            end
+        end
+        
+        % save time be resampling only a subregion
+        ind = find(roimaskCorners);
+        [I,J,K] = ind2sub([x_dim+1 y_dim+1 z_dim],ind);
+        indx = min(I)-3:max(I)+3;
+        indy = min(J)-3:max(J)+3;
+        indz = min(K)-3:max(K)+3;
+        
+        indx = indx(indx >= 1 & indx <= x_dim);
+        indy = indy(indy >= 1 & indy <= y_dim);
+        indz = indz(indz >= 1 & indz <= z_dim);
+        
+        % convert the corners to a 3-D roi mask
+        roimask(indx,indy,indz) = ceil(gridResample3D(xCorners,yCorners,z,roimaskCorners,x(indx),y(indy),z(indz)));
+        
+        % save the indices of the roi mask
+        ROIS{roi_num}.ind = int32(find(roimask ~= 0));
+        ROIS{roi_num}.dim = [x_dim y_dim z_dim];
+        ROIS{roi_num}.pixdim = [x_pixdim y_pixdim z_pixdim];
+        ROIS{roi_num}.start = [x_start y_start z_start];
+        ROIS{roi_num}.start_ind = start_ind;
+        fprintf('Converted %s vectors to an roi mask.\n',ROIS{roi_num}.name);
+    end
+    
+    % save ROIS to the Geometry structure
+    Geometry.ROIS = ROIS;
+    clear ROIS;
+end
+
+if strcmp(readTrial,'yes')
+    % read in the Trial information
+    fid_trialfile = fopen(trialfile,'r');
+    tline = fgets(fid_trialfile);
+
+    structstack = {};   % stack for determining which field we're inside
+    structind = [];    % structure index values for each field
+    lines = 1;
+
+    while tline ~= -1   % read through the trial file
+        tline = regexprep(tline,'"','''');  % replace " with '
+        tline = regexprep(tline,'[','');  % bad characters to remove
+        tline = regexprep(tline,']','');  % bad characters to remove
+        tline = regexprep(tline,'#','num');
+        tline = regexprep(tline,'\\','''');
+        % fprintf(tline);
+        if length(findstr(tline,'{'))
+            % mark the structure we just entered, removing any blank space
+            structstack{length(structstack) + 1} = regexprep(strtok(tline,'='),' ','');
+            if length(structstack) == 1
+                structstack{1} = 'Trial';  % make sure structure name is always the same
+            end
+            if strcmp(structstack{end},'Beam')  % found a beam structure
+                % find the number of existing beams at this location
+                eval_line = ['beams = getfield(' structstack{1}];
+                for k=2:length(structstack)
+                    eval_line = [eval_line ',''' structstack{k} ''''];
+                end
+                eval_line = [eval_line ');'];
+                try
+                    eval(eval_line);
+                    num_beams = length(beams) + 1;  % this is the current beam number
+                catch  % there was an error, so there are no beams yet
+                    num_beams = 1;  % we're on the first beam now
+                end
+            end
+        elseif length(findstr(tline,'}'))
+            structstack(end) = [];  % step outside of the current field
+        else  % the line must be a subfield assignment
+            if strcmp(structstack{end},'Points')  % this case treated specially
+                % read-in values until the next '}' is reached
+                lines = lines + 1;
+                mlc_vec = ['['];
+                while ~length(findstr(tline,'}'))
+                    % build up the vector of MLC positions
+                    mlc_vec = [mlc_vec regexprep(tline,' ','')];
+                    tline = fgets(fid_trialfile);
+                    % fprintf(tline);
+                    lines = lines + 1;
+                end
+                mlc_vec = [mlc_vec ']'];
+                % evaluate the points vector
+                eval_statement = structstack{1};
+                for k=2:length(structstack)
+                    if strcmp(structstack{k},'Beam')
+                        % include the beam number
+                        eval_statement = [eval_statement '.' structstack{k} '(' num2str(num_beams) ')'];
+                    else  % non-beam structure, so don't worry about it
+                        eval_statement = [eval_statement '.' structstack{k}];
+                    end
+                end
+                eval_statement = [eval_statement ' = ' mlc_vec ';'];
+                eval(eval_statement);
+                structstack(end) = [];  % step outside of the current Points[] field
+            else
+                % see if there are any subfields in the expression
+                % evaluate the statement
+                eval_statement = structstack{1};
+                for k=2:length(structstack)
+                    if strcmp(structstack{k},'Beam')
+                        % include the beam number
+                        eval_statement = [eval_statement '.' structstack{k} '(' num2str(num_beams) ')'];
+                    else  % non-beam structure, so don't worry about it
+                        eval_statement = [eval_statement '.' structstack{k}];
+                    end
+                end
+                eval_statement = [eval_statement '.' tline];
+                eval(eval_statement);
+            end
+        end
+        tline = fgets(fid_trialfile);
+        lines = lines + 1;
+        if isempty(structstack)
+            break;   % we're done if we're no longer inside any sets of curly brackets
+        end
+    end
+    fclose(fid_trialfile);
+
+    % dose grid dimensions
+    xdim = Trial.DoseGrid.Dimension.X;
+    ydim = Trial.DoseGrid.Dimension.Y;
+    zdim = Trial.DoseGrid.Dimension.Z;
+    siz = [xdim ydim zdim];  % dosegrid size vector
+
+    % vector describing the start location of the dose grids
+    xstart = Trial.DoseGrid.Origin.X;
+    ystart = Trial.DoseGrid.Origin.Y;
+    zstart = Trial.DoseGrid.Origin.Z;
+    start = [xstart ystart zstart];  % dosegrid start vector
+
+    % vector describing the size of each voxel
+    dx = Trial.DoseGrid.VoxelSize.X;
+    dy = Trial.DoseGrid.VoxelSize.Y;
+    dz = Trial.DoseGrid.VoxelSize.Z;
+    voxel_size = [dx dy dz];
+
+    % Read in the beamlet dose distributions, DRRs, and fluences and store them in
+    % Trial.BeamList.Beam(k), where k is the beam number
+    for k=0:num_beams-1
+        % read-in the dose distribution
+        num_vec = num2str(k*5+4);
+        % ensure that num_vec has a length of 3
+        while length(num_vec) < 3
+            num_vec = ['0' num_vec];
+        end
+        dose_file = [dose_file_base_name '.' num_vec];
+        fid = fopen(dose_file,'rb','ieee-be');
+        if fid ~= -1
+            dose = fread(fid,'float=>float');
+            fclose(fid);
+            try
+                dose = reshape(dose,siz);
+                % flip the dose distribution in accordance with the CT image
+                for j=1:siz(3)
+                    dose(:,:,j) = single(fliplr(double(dose(:,:,j))));
+                end
+                % save the dose
+                Trial = setfield(Trial,'BeamList','Beam',{k+1},'dose',dose);
+            catch
+                fprintf('Dose for beam %g did not meet specifications so it was not read-in.\n',k+1);
+            end
+        end
+
+        % read-in the DRR
+        % find the size of this drr
+        drr_x = Trial.BeamList.Beam(1).FilmImageList.FilmImage.Image.Dimension.X;
+        drr_y = Trial.BeamList.Beam(1).FilmImageList.FilmImage.Image.Dimension.Y;
+        siz_drr = [drr_x drr_y];
+
+        num_vec = num2str(k*5+3);
+        % ensure that num_vec has a length of 3
+        while length(num_vec) < 3
+            num_vec = ['0' num_vec];
+        end
+        drr_file = [dose_file_base_name '.' num_vec];
+        fid = fopen(drr_file,'rb','ieee-be');
+        if fid ~= -1
+            drr = fread(fid,'float');
+            fclose(fid);
+            try  % proceed only if the drr was read-in properly
+                drr = reshape(drr,siz_drr);
+                % save the DRR
+                Trial = setfield(Trial,'BeamList','Beam',{k+1},'FilmImageList','FilmImage','Image','drr',drr);
+            catch
+                fprintf('DRR for beam %g did not meet specifications so it was not read-in.\n',k+1);
+            end
+        end
+
+        % read in the fluence map
+        num_vec = num2str(k*5+2);
+        fluence_siz = [101 101 3];  % all fluence must have this size otherwise and error will be thrown
+        % ensure that num_vec has a length of 3
+        while length(num_vec) < 3
+            num_vec = ['0' num_vec];
+        end
+        fluence_file = [dose_file_base_name '.' num_vec];
+        fid = fopen(fluence_file,'rb','ieee-be');
+        if fid ~= -1
+            fluence = fread(fid,'float');
+            fclose(fid);
+            % check to see if fluence matches dose grid size
+            if Trial.FluenceGridMatchesDoseGrid ~= 1
+                error('Fluence grid does not match dose grid for beam %g, don''t know what to do here.',k+1);
+            end
+            try  % proceed only if the fluence is read-in properly
+                fluence = reshape(fluence,fluence_siz);
+                % save the fluence
+                Trial = setfield(Trial,'BeamList','Beam',{k+1},'FluenceMap',fluence);
+            catch
+                fprintf('Fluence map for beam %g did not meet specifications so it was not read-in.\n',k+1);
+            end
+        end
+    end
+
+    % add up all of the doses
+    dose_tot = zeros(siz);
+    for k=1:num_beams
+        if isfield(Trial.BeamList.Beam(k),'dose') & ~isempty(Trial.BeamList.Beam(k).dose)
+            dose_tot = dose_tot + double(Trial.BeamList.Beam(k).dose);
+        end
+    end
+    Geometry.Trial = Trial;
+    clear Trial;
+end
+
+save(save_file,'Geometry');

linlsqOptSetup.asv

@@ -0,0 +1,153 @@
+%% linlsqOptSetup.m
+
+% Stephen R Bowen April 2009
+%
+% Runs in MATLAB or writes setup files necessary to carry out an
+% optimization of photon beamlet or proton spot weight distributions by
+% iteratitively minimizing a linear least squares objective function.  The
+% objective function has the following components:
+% 
+% -Tissue relative importance alpha/Ntissue (alpha binary)
+% -Overdose and underdose relative penalties in each tissue Bplus & Bminus
+% -Maximum and minimum doses at each voxel within a tissue dplus & dminus
+% -Overdose-volume and underdose-volume relative penalties BVplus & BVminus
+% -Maximum and minimum doses dVplus & dVminus evaluated at each relative
+%  tissue maximum and minimum volume percentatges Vplus & Vminus 
+
+%%
+
+%%%%% Start of user supplied inputs %%%%%
+
+% specify relative location of Geometry file including ROI masks
+geometryFile = '../plan/HN003/HN003final.mat';
+
+% define optimization settings
+optSettings = [];
+optSettings.optInfo.Niterations = 1000;                                     % total number of iterations of beamlet weight optimization
+optSettings.optInfo.Nperbatch = 500;                                        % iteration interval between writing planned dose and weights to file
+optSettings.optInfo.optRun = 'no';                                          % flag to run optimization in MATLAB
+optSettings.optInfo.waitForResults = 'no';                                  % flag to load final dose distribution and DVH
+optSettings.optInfo.saveOptInfo = 'yes';                                    % save opt settings in optinput.txt file located in opt folder
+optSettings.beamletInfo.saveBeamlets = 'no';                                % save beamlets when dose calc performed in MATLAB and not Condor
+optSettings.beamletInfo.saveBeamletHeader = 'yes';                          % save beamlet info in beamlet_header.txt file located in beamlet folder
+optSettings.beamletInfo.beamletFolder = 'HN003beamletbatches';              % specify location of beamlet folder relative to opt folder
+optSettings.beamletInfo.Nbeamletbatches = 492;                              % total number of beamlet batches
+optSettings.prescInfo.savePresc = 'yes';                                    % save presc info in prescription.txt file located in input folder
+optSettings.optInfo.optFolder = 'C:\data_temp\steveTPS\opt\linlsqOpt';      % directory where opt executable is located must be absolute path
+optSettings.optInfo.modFactor = 3.3;                                        % mod factor truncating maximum beamlet weight relative to average weight
+
+%   Prescription table
+%   (code name)     (Pinnacle tissue name) (alpha)  (Bplus) (dplus)  (Bminus)  (dminus) (BVplus) (dVplus) (Vplus) (BVminus)  (dVminus)  (Vminus)        
+prescTable = {      
+    'GTV'                'GTV'               0         0        0        0         0        0        0        0       0         0         0        
+    'CTV60'              'CTV 60'            0         0        0        0         0        0        0        0       0         0         0                
+    'CTV54'              'CTV 54'            0         0        0        0         0        0        0        0       0         0         0          
+    'CTV-rt-scv-5'       'CTV rt scv 54'     0         0        0        0         0        0        0        0       0         0         0                  
+    'CTV-lt-scv-50'      'CTV lt scv 50'     0         0        0        0         0        0        0        0       0         0         0                  
+    'larynx'             'larynx'            1         1        0        0         0        0        0        0       0         0         0                
+    'opticApparatus'     'optic apparatus'   1         1        0        0         0        0        0        0       0         0         0             
+    'oralCavity'         'oral cavity'       1         1        0        0         0        0        0        0       0         0         0        
+    'RtParotid'          'rt parotid'        1         1        0        0         0        0        0        0       0         0         0         
+    'LtParotid'          'lt parotid'        1         1        0        0         0        0        0        0       0         0         0         
+    'Cord'               'Cord'              1         1        0        0         0        0        0        0       0         0         0         
+    'Brainstem'          'Brainstem'         1         1        0        0         0        0        0        0       0         0         0         
+    'eyes'               'eyes'              1         1        0        0         0        0        0        0       0         0         0        
+    'couch'              'couch'             0         0        0        0         0        0        0        0       0         0         0        
+    'external'           'external'          0         0        0        0         0        0        0        0       0         0         0        
+    'ext+1cm'            'ext + 1 cm'        0         0        0        0         0        0        0        0       0         0         0        
+    'ring'               'ring'              0         1        0        0         0        0        0        0       0         0         0 
+    'PTV70'              'PTV  70'           1       200       70      200        70        0        0        0       0         0         0 
+    'PTV60'              'PTV  60'           0         0        0        0         0        0        0        0       0         0         0 
+    'PTV54'              'PTV 54'            1        10        0        0         0        0        0        0       0         0         0 
+    'PTV50'              'PTV 50'            1        10        0        0         0        0        0        0       0         0         0 
+    'LResidParotid'      'L resid parotid'   1         1        0        0         0        0        0        0       0         0         0    
+    'normal'             'normal'            1        10        0        0         0        0        0        0       0         0         0 
+    'PTV60min70'         'PTV 60 - PTV 70'   1        10       60      200        60        0        0        0       0         0         0 
+    }; 
+
+%%%%% End of user supplied inputs %%%%%
+
+%%
+
+% total number of initial non-zero weighted beamlets from dose calc
+optSettings.beamletInfo.Nbeamlets = beamletNum(optSettings.beamletInfo.beamletFolder,optSettings.beamletInfo.Nbeamletbatches);
+
+load(geometryFile);  % load geometry file
+
+% match isocenter of coordinate system used to calculate photon beamlets or
+% proton spots, defined by the patient geometry
+x_iso = Geometry.start(1) + (Geometry.voxel_size(1).*double(size(Geometry.data,1)))./2;
+y_iso = Geometry.start(2) + (Geometry.voxel_size(2).*double(size(Geometry.data,2)))./2;
+z_iso = Geometry.start(3);
+
+iso = [x_iso y_iso z_iso];
+
+% size of CT grid
+[M,N,Q] = size(Geometry.data);
+
+% Create the prescription structure
+presc = [];   % erase the old prescription structure, if it exists
+presc.siz = [M N Q];     % dimensions of the prescription and tissue masks
+
+ROIFlags = zeros(1,size(prescTable,1));  % flags to ensure all ROIs are found
+
+% assign the parameters to the prescription structure
+for j=1:size(prescTable,1)
+    % pull out the prescription indices for each tissue
+    tissName = prescTable{j,2};  % current tissue name
+
+    % find the current tissue in the Geometry.ROIS cells
+    for k=1:length(Geometry.ROIS)
+        if strcmp(tissName,Geometry.ROIS{k}.name)
+            ind = Geometry.ROIS{k}.ind;  % extract ROI indices
+            break;
+        end
+    end
+    
+    presc.tissue(j).name = prescTable{j,1};
+    presc.tissue(j).ind = ind;
+    presc.tissue(j).alpha = prescTable{j,3};
+    presc.tissue(j).betaPlus = prescTable{j,4};
+    presc.tissue(j).dPlus = single(zeros(size(ind)) + prescTable{j,5});
+    presc.tissue(j).betaMinus = prescTable{j,6};
+    presc.tissue(j).dMinus = single(zeros(size(ind)) + prescTable{j,7});
+    presc.tissue(j).betaVPlus = prescTable{j,8};
+    presc.tissue(j).dVPlus = prescTable{j,9};
+    presc.tissue(j).vPlus = prescTable{j,10};
+    presc.tissue(j).betaVMinus = prescTable{j,11};
+    presc.tissue(j).dVMinus = prescTable{j,12};
+    presc.tissue(j).vMinus = prescTable{j,13};
+    ROIFlags(j) = 1;
+end
+
+% ensure that all of the ROIS were found
+if sum(ROIFlags) ~= length(ROIFlags)
+    error('Missed one or more ROIS.\n');
+else
+    fprintf('All ROIs located successfully for downsampled prescription.\n');
+end
+
+optSettings.prescInfo.presc = presc;  % set prescription
+
+% Save prescription and beamlet information, which is necessary for initial
+% guess calculation.
+linlsqOptimization(optSettings);  
+
+% calculate the initial guess if flagged
+[w0,dose0] = calcNormalizedInitialGuess([optSettings.optInfo.optFolder '/optInput.txt']);
+optSettings.initialGuessInfo.w0 = w0;
+optSettings.initialGuessInfo.saveInitialGuess = 'yes';                     
+linlsqOptimization(optSettings);
+
+if (optSettings.optInfo.optRun == 'no')
+    fprintf('Optimization setup complete.\n');
+else
+    [weights, dose] = linlsqOptimization(optSettings);
+    fprintf('Optimization run complete.\n');
+    
+    if (optSettings.optInfo.waitForResults == 'yes')
+    [optSettings,optResults] = loadOptResults([optSettings.optInfo.optFolder '/optInput.txt'],'last');
+    frpintf('Final dose distribution and DVHs loaded.\n');
+    end
+    
+end

opt/linlsqOpt/AV001cbeamlets/beamlet_header.txt

@@ -0,0 +1,12 @@
+beamletDim
+113 113 19
+Nbeamlets
+11200
+Nbeamletbatches
+175
+beamletFolder
+AV001cbeamlets
+beamletBatchBaseName
+beamletbatch
+beamletBatchExtension
+bin

opt/linlsqOpt/AV001fbeamlets/beamlet_header.txt

@@ -0,0 +1,12 @@
+beamletDim
+120 120 21
+Nbeamlets
+11322
+Nbeamletbatches
+177
+beamletFolder
+AV001fbeamlets
+beamletBatchBaseName
+beamletbatch
+beamletBatchExtension
+bin

opt/linlsqOpt/AV003fbeamlets/beamlet_header.txt

@@ -0,0 +1,12 @@
+beamletDim
+155 109 20
+Nbeamlets
+9600
+Nbeamletbatches
+150
+beamletFolder
+AV003fbeamlets
+beamletBatchBaseName
+beamletbatch
+beamletBatchExtension
+bin

opt/linlsqOpt/AV004cbeamlets/beamlet_header.txt

@@ -0,0 +1,12 @@
+beamletDim
+87 107 19
+Nbeamlets
+12592
+Nbeamletbatches
+197
+beamletFolder
+AV004cbeamlets
+beamletBatchBaseName
+beamletbatch
+beamletBatchExtension
+bin

opt/linlsqOpt/HN002bbeamlets/beamlet_header.txt

@@ -0,0 +1,12 @@
+beamletDim
+84 104 15
+Nbeamlets
+7872
+Nbeamletbatches
+123
+beamletFolder
+HN002bbeamlets
+beamletBatchBaseName
+beamletbatch
+beamletBatchExtension
+bin