%This function finds the top and bottom edges of the tube and uses these to
%find the tube axis. The velocity is then integrated radially and
%circumferentially to find the flow rate.

%This version uses specific lines 'reg_lines' and is for phantom b which
%has an image artefact

%For full details see "Ultrasound Imaging Velocimetry with interleaved
%images for improved pulsatile arterial flow measurements: a new correction
%method, experimental and in vivo validation", J Royal Soc Interface, 2017

%K H Fraser, C Poelma 2016


function [flow2, Flow,upperEdge,lowerEdge]=flowRatePulse_f(xPix, yPix, times, timesteps, depth, lines, VxC, VyC, Iavg, iastep1, iastep2, minTop, maxTop, minBot, maxBot, reg_lines);

for t=1:timesteps
    meanData(:,:,t)=Iavg(:,:,t)'; %interlaced
end

for t=1:timesteps
    %-- this section is where you need to adjust the limits according to your images
    A(:,:)=log(1+abs(hilbert(squeeze(meanData(:,:,t)')))); %find the inner edges of the tube
    for n=1:lines
         % figure; plot(A(n,:,t)); %uncomment this the first time through
        windowSize = 10;
        B=filtfilt(ones(1,windowSize)/windowSize,1,A(n,:));
        windowSize = 17;
        C=filtfilt(ones(1,windowSize)/windowSize,1,A(n,:));
         %   plot(B,'g'); hold on; plot(C,'r') %uncomment this the first time through
        D=B-C;
         %    hold on; plot(D); %uncomment this the first time through
        maximum=max(D(minTop:maxTop));
        temp=find(D==maximum);
        pos1(n,t)=temp(1);
        maximum=max(D(minBot:maxBot));
        temp=find(D==maximum);
        pos2(n,t)=temp(end);
    end %XX
end

%--end of limits section
for t=1:timesteps
    A(:,:)=log(1+abs(hilbert(squeeze(meanData(:,:,t)'))));
    fit1{t} = fit([1:128]',pos1(:,t),'poly1'); %linear fit to upper artery edge in pixels
    fit2{t} = fit([1:128]',pos2(:,t),'poly1'); %linear fit to lower artery edge in pixels
    %this figure is for you to check the edges are in the right places
%     if t./5==round(t./5)
%         figure; imagesc(A'); colormap gray;
%         hold on; plot(fit1{t});
%         hold on; plot(fit2{t},'b');
%         legend('upper wall','lower wall');
%     end
end
for t=1:timesteps
    m=1;
    for n=1:iastep1:lines %convert these edges to coordinate system
        %y coordinates in m
        yCoords(m,:,t)=([1:iastep2:length(meanData)]-(fit2{t}(n)+fit1{t}(n))/2)*yPix;
        upperEdge(m,t)=fit1{t}(n); %upper and lower edge positions in pixels
        lowerEdge(m,t)=fit2{2}(n);
        diameter(m,t)=fit2{2}(n)-fit1{2}(n); %diameter in pixels
        m=m+1;
    end
end

%find magnitude from velocity components
Vmag=(VxC.^2+VyC.^2).^0.5;
Vmag=Vmag.*sign(VxC);

%Calculate flow
for t=1:timesteps
    for n=1:lines/iastep1-1
        r=1;
        for R=round(upperEdge(n,t)/iastep2)+6:round(lowerEdge(n,t)/iastep2)-6 %the -5 here means we remove the contribution from close to the lower wall
            %flow for each semi-cylinder, for each column of interrogation regions, in m^3/s
           dFlow(n,r,t)=VxC(n,R,t)*abs(yCoords(n,R,t))*iastep2*yPix*pi;
            r=r+1;
        end
    end
end
flow=sum(-dFlow,2); %flow for each column of interrogation regions in m^3/s
region=round(reg_lines./iastep1);
flow2=squeeze(mean(flow(region,:,:)*1e6*60,1)); %flow in ml/min

Flow=mean(flow2) %mean flow rate, in ml/min