function newBeam = matRad_discardApertures(beam,numToKeep)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The sequencing algorithm generates an a priori unkown number of aperture.
% We only want to keep a certain number of them (numToKeep).  These will be
% the ones with the highest intensity-area product.
%
%
% call
%   beam =
%   matRad_discardApertures(beam,numToKeep)
%
% input
%   beam:               beam struct containing original shapes and
%                       intensities
%
%   numToKeep:          number of apertures to keep
%
% output
%   beam:               beam struct containing shapes and re-scaled
%                       intensities of the apertures we are keeping
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


newBeam = beam;
newBeam.shapes = zeros(size(newBeam.shapes,1),size(newBeam.shapes,2),numToKeep);
newBeam.shapesWeight = zeros(numToKeep,1);

%Find the numToKeep apertures having the highest dose-area product
numToKeep = min(numToKeep,beam.numOfShapes);

DAP = zeros(beam.numOfShapes,1);
comPos = zeros(beam.numOfShapes,1);

for shape = 1:beam.numOfShapes
    DAP(shape) = nnz(beam.shapes(:,:,shape)).*beam.shapesWeight(shape);
    x = repmat(1:size(beam.shapes(:,:,shape),2),size(beam.shapes(:,:,shape),1),1);
    comPosRow = sum(beam.shapes(:,:,shape).*x,2)./sum(beam.shapes(:,:,shape),2);
    comPos(shape) = mean(comPosRow(~isnan(comPosRow),1));
end

%{
% This code will sort the aperture shapes in increasing order of centre of mass.
However, some algorithms (in particular, Siochi) already sort shapes in
increasing (left to right) leaf position.  So this is not necessary (and
will in fact some times change the order of the apertures generated by
Siochi).

% sort comPos into ascending order
[comPosSorted,comPosSortInd] = sort(beam.comPos);

if any(comPosSorted ~= beam.comPos)
    tempShapes = zeros(size(beam.shapes));
    tempShapesWeight = zeros(size(beam.shapesWeight));
    tempDAP = zeros(size(beam.DAP));
    
    for shape = 1:beam.numOfShapes
        ind = comPosSortInd(shape);
        
        tempShapes(:,:,shape) = beam.shapes(:,:,ind);
        tempShapesWeight(shape) = beam.shapesWeight(ind);
        tempDAP(shape) = beam.DAP(ind);
    end
    
    % now the shapes should be sorted in increasing center of mass position
    beam.shapes = tempShapes;
    beam.shapesWeight = tempShapesWeight;
    beam.DAP = tempDAP;
    beam.comPos = comPosSorted;
end
%}

[~,comPosToDAPSort] = sort(DAP,'descend');

totDAP_all = sum(DAP(:));
totDAP_keep = sum(DAP(comPosToDAPSort(1:numToKeep)));

segmentKeep = 1;

%Keep only those numToKeep apertures with the highest DAP
%Preserve the shapes of the apertures, but scale the weights so
%that the total DAP is kept
for shape = 1:beam.numOfShapes
    if comPosToDAPSort(shape) <= numToKeep
        newBeam.shapes(:,:,segmentKeep) = beam.shapes(:,:,shape);
        tempNewDAP = totDAP_all*DAP(shape)/totDAP_keep;
        newBeam.shapesWeight(segmentKeep) = tempNewDAP/(nnz(newBeam.shapes(:,:,segmentKeep))); %sequencing.beam.shapesWeight(sequencing.beam.segmentSortedDAP(segment))
        
        segmentKeep = segmentKeep+1;
    else
        continue
    end
end

newBeam.numOfShapes = numToKeep;
