% Loads the results of dose painting optimizations into the results
% structure.

resultsFile = 'resultsTomo.mat';

% folder containing optimization results
resultsFolder = 'F:\rtflynn\projects\dosePainting\dataAndAnalysis\protonsTomoContrastPhantom\tomo\optimization';

inputFileBaseName = 'inputTomo';
inputFileExtension = 'txt';

% phantom radii
Rp = [8 10 12 15];

% radii for centers of boost regions
Rb = [5 4 3 2];  

% contrasts
contrasts = [-0.1 -0.2 -0.5 -1 0.1 0.2 0.5 1 2];

% parameters for the dose painting phantom
pProp.dx = 0.1;       % pixel side length (cm)
pProp.Rt = 6;       % background radius (cm)
pProp.Nboost = 6;   % number of inserts, equi-angularly spaced
pProp.rHi = 1;        % highest insert radius
pProp.rLo = 0.15;      % lowest insert radius
pProp.rhoPhant = 1;   % density of the phantom (g/cc)

results = [];  % results structure

for m=1 %1:length(Rp)
    for n=1:length(Rb)
        for q=1:length(contrasts)
            pProp.Rp = Rp(m);   % phantom size (cm)
            pProp.Rb = Rb(n);   % distance between phantom center and insert centers (cm)
            pProp.contrast = contrasts(q);   % contrast between background and each insert
            pProp.Npix = ceil((2*pProp.Rp + 1)/pProp.dx);  % number of pixels on each side of the grid
            results{m,n,q}.pProp = pProp;
            inputFileName = [resultsFolder '/' inputFileBaseName num2str((q-1) + n*length(contrasts) + m*length(contrasts)*length(Rb)) '.txt'];
            [optSettings,optResults] = loadOptResults(inputFileName,'last');
            
            % save the dose distribution, beamlet weights
            results{m,n,q}.dose = optResults.dose;
            results{m,n,q}.weights = optResults.weights;
            
            % save other useful information
            results{m,n,q}.objFunc = optResults.objFunc;
            results
            
            
            % loop through all tissues, extracting names and dvhs
            for k=1:length(optResults.presc.tissue)
                results{m,n,q}.tissue(k).name = optResults.presc.tissue(k).name;
                results{m,n,q}.tissue(k).dvhbins = optResults.presc.tissue(k).dvhbins;
                results{m,n,q}.tissue(k).dvh = optResults.presc.tissue(k).dvh;
                
                % create the differential DVH, calculate some dose statistics
                dvhDiff = [diff(results{m,n,q}.tissue(k).dvh) 0];
                dbins = results{m,n,q}.tissue(k).dvhbins;
				ED = sum(dvhDiff.*dbins)/sum(dvhDiff);  % expectation value of dose
                ED2 = sum(dvhDiff.*dbins.^2)/sum(dvhDiff);  % expectation of dose squared
                Dvar = ED2 - ED^2;  % dose variance
                Dstd = sqrt(Dvar);  % dose standard deviation
                
                results{m,n,q}.tissue(k).Davg = ED;
                results{m,n,q}.tissue(k).Dstd = Dstd;
            end
        end
    end
end

% 
% 
% for m=1:M
%     for n=1:N
%         for q=1:Q
%             siz = results{m,n,q}.optSettings.presc.siz;
%             Nbeamlets = numel(results{m,n,q}.optSettings.w0);
%             Niterations = results{m,n,q}.optSettings.Niterations;
%             
%             % load the dose results file
%             doseFolder = [resultsFolder '/' results{m,n,q}.optSettings.outputFolder];
%             doseFile = [doseFolder '/dosebatch' num2str(Niterations) '.img'];
%             fid = fopen(doseFile,'rb');
%             dose = fread(fid,'float=>float');
%             fclose(fid);
%             dose = reshape(dose,siz);
%             
%             % load the beamlet weights
%             weightsFolder = [resultsFolder '/' results{m,n,q}.optSettings.outputFolder];
%             weightsFile = [doseFolder '/weightbatch' num2str(Niterations)
%             '.img'];
%             fid = fopen(weightsFile,'rb');
%             weights = fread(fid,'float=>float');
%             fclose(fid);
%             weights = reshape(weights,[1 Nbeamlets]);
%             
%             % load the objective function
%             objFolder = [resultsFolder '/' results{m,n,q}.optSettings.outputFolder];
%             objFile = [doseFolder '/' results{1,1,1}.optSettings.objFuncFileName];
%             fid = fopen(objFile,'rb');
%             objFunc = fread(fid,'float=>float');
%             fclose(fid);
%             objFunc = reshape(objFunc,[1 Niterations+1]);
%             
%             % save results
%             results{m,n,q}.weights = weights;
%             results{m,n,q}.dose = dose;
%             results{m,n,q}.objFunc = objFunc;
%             
%             % calculate the dvh for each tissue
%             for k=1:length(results{m,n,q}.optSettings.presc.tissue)
%                 mask = zeros(siz);
%                 mask(results{m,n,q}.optSettings.presc.tissue(k).ind) = 1;
%                 results{m,n,q}.optSettings.presc.tissue(k).dvhbins = dbins;
%                 results{m,n,q}.optSettings.presc.tissue(k).dvh = dvh(results{m,n,q}.dose,mask,dbins);
%                 dvhDiff = [diff(results{m,n,q}.optSettings.presc.tissue(k).dvh) 0]; % differential DVH
% 				ED = sum(dvhDiff.*dbins)/sum(dvhDiff);  % expectation value of dose
%                 ED2 = sum(dvhDiff.*dbins.^2)/sum(dvhDiff);  % expectation of dose squared
%                 Dvar = ED2 - ED^2;  % dose variance
%                 Dstd = sqrt(Dvar);  % dose standard deviation
%                 
%                 results{m,n,q}.optSettings.presc.tissue(k).Davg = ED;
%                 results{m,n,q}.optSettings.presc.tissue(k).Dstd = Dstd;
%             end
%         end
%     end
%     fprintf('Finished m = %g of %g\n',m,M);
% end
% 
% results = results;
% 
% save(resultsFile,'results');