function [qbin,dqbin,Ibin,Idqb]=binning_equalintensity(q,I,nbins)

%this function implements binning of intensity into bins of equal
%intensity. qbin is the middle of the bin, and dqbin is the size of
%the bin (in the same units as q). Idqb is the normalised binned intensity,
%normalised to the bin size.

%reshape q and I to x by 1 vector
q=reshape(q,size(q,1)*size(q,2),1);
I=reshape(I,size(I,1)*size(I,2),1);

%sort q,I to q
[qsort,qsorti]=sort(q);
Isorti=I(qsorti);

%determine amount of intensity per bin
Iperbin=sum(Isorti)/(nbins+1);

%cumsum. This we're going to need a few times
Icum=cumsum(Isorti);

binnumber=ceil(Icum/Iperbin); %now every pixel has a bin number assigned

% newbinmarker=diff(binnumber);
% newbinmarker(length(newbinmarker+1))=[1 newbinmarker];

%first implementation: using a loop function
%find indices
endi=find(binnumber==binnumber(1),1,'last');
bdiv=binnumber(1);
Ibin(1:binnumber(1))=Icum(endi)/bdiv;
dqbin(1:binnumber(1))=(qsort(endi)-qsort(1))/bdiv;
qbin(1:binnumber(1))=(qsort(endi)+qsort(1)).*(1:bdiv)./(bdiv.*2);
leftoff=bdiv+1;
for bini=bdiv+1:max(binnumber)
    %find indices
    starti=find(binnumber==bini,1,'first');
    endi=find(binnumber==bini,1,'last');
    if ~isempty(starti)
        %determine number of bins in this section
        bdiv=binnumber(starti)-binnumber(starti-1);
        Ibin(leftoff:bini)=(Icum(endi)-Icum(starti-1))/bdiv;
        dqbin(leftoff:bini)=(qsort(endi)-qsort(starti-1))/bdiv;
        qbin(leftoff:bini)=qsort(starti)+dqbin(leftoff:bini).*(1:bdiv)./(2);
        leftoff=bini+1;
%         Ibin(bini)=Icum(endi)-Icum(starti-1);
%         dqbin(bini)=qsort(endi)-qsort(starti);
%         qbin(bini)=(qsort(endi)+qsort(starti))/2;
    end %otherwise skip
end
Idqb=Ibin./dqbin;