WiscPlan_olderMV/opt/linlsqUtil/calcNormalizedInitialGuess.m

function [w0,dose0] = calcNormalizedInitialGuess(varargin)
% [w0,dose0] = calcNormalizedInitialGuess(B,presc) calculates an initial
% guess for the linear least squares optimizer based on the beamlets, B,
% and the prescription structure, presc.  B is a cell array of Dij cell
% arrays, where each element in Dij is a Matlab sparse matrix.  The
% structure presc has the form:
% presc.siz = [M N Q]  % dimensions of prescription, beamlet must match
% presc.tissue(:)  % structure array of tissue prescription
%                .dMinus  % minimum dose tissue must receive
% 
% [w0,dose] = calcNormalizedInitialGuess(optInputFile)
% Calculates the initial guess for beamlets and a prescription that are not
% loaded into Matlab, but can be accessed using information in
% optInputFile, which is an input file to the linear least squares
% optimizer.
% 
% The goal is for all of the initial beams to have the same weight,
% designed such that the integral dose to the structures that were
% prescribed non-zero dose is equal to the prescribed integral dose.
% 
% RTF 1/19/07

% Check which mode the user has requested
if numel(varargin) == 2
    if ~iscell(varargin{1})
        error('First input must be a cell array of beamlets.');
    elseif ~isstruct(varargin{2})
        error('Second input must be a prescription structure.');
    else
        B = varargin{1};
        presc = varargin{2};
    end
        
    Nbeamlets = 0;
    dose = zeros(presc.siz);
    for k=1:numel(B)
        for m=1:numel(B{k}.Dij)
            if ~isempty(B{k}.Dij{m}.non_zero_indices)
                
                dose(B{k}.Dij{m}.non_zero_indices) = dose(B{k}.Dij{m}.non_zero_indices) ...
                    + (B{k}.Dij{m}.non_zero_values);
                
            end
            
            Nbeamlets = Nbeamlets + 1;
        end
    end
elseif numel(varargin) == 1
    if ~ischar(varargin{1});
        error('Input must be a character string.');
    else
        optInputFile = varargin{1};
    end
    
    optSettings = loadOptSettings(optInputFile);
    presc = optSettings.prescInfo.presc;  % extract the prescription information
    
    % Extract the folder name that contains the input file
    inputFileNameRev = fliplr(optInputFile);  % flip the input filename around
    % pop off the reversed file name
    [fileNameRev,inputFolderRev] = strtok(inputFileNameRev,{'/','\'});
    inputFolder = fliplr(inputFolderRev);
    
    % calculate the initial dose for all unit weights
    dose = zeros(presc.siz,'single');
    Nbeamlets = 0;
    for k=1:optSettings.beamletInfo.Nbeamletbatches
        beamletBatchFile = [inputFolder '/' optSettings.beamletInfo.beamletFolder ...
            '/' optSettings.beamletInfo.beamletBatchBaseName num2str(k-1) '.' ...
            optSettings.beamletInfo.beamletBatchExtension];
        beamlets = open_beamlet_batch(beamletBatchFile);  % load the beamlet batch
        for m=1:length(beamlets)
            if ~isempty(beamlets{m}.non_zero_indices)
                
                    dose(beamlets{m}.non_zero_indices) = dose(beamlets{m}.non_zero_indices) ...
                    + beamlets{m}.non_zero_values;
                
                
            end
            Nbeamlets = Nbeamlets + 1;
            
        end
    end
    if Nbeamlets ~= optSettings.beamletInfo.Nbeamlets
        error(['Beamlets read = ' num2str(Nbeamlets) '. Beamlets expected = ' num2str(optSettings.beamletInfo.Nbeamlets) '.']);
    end
end

dIntTumor = 0;  % integral dose in tumor
dIntTumorPresc = 0;  % prescribed integral tumor dose

% calculate integral dose in the non-zero prescription regions
for k=1:length(presc.tissue)
    if presc.tissue(k).dMinus > 0
        dIntTumorPresc = dIntTumorPresc + sum(presc.tissue(k).dMinus);
        dIntTumor = dIntTumor + sum(dose(presc.tissue(k).ind));
    end
end

if dIntTumor == 0
    error('Zero tumor dose -- respecify initial beamlet weights.');
end

w0 = ones(Nbeamlets,1);   % allocate a vector for the initial beamlet weights
w0 = w0*dIntTumorPresc/dIntTumor;  % normalized initial guess
dose0 = dose*dIntTumorPresc/dIntTumor;