%This script generates the images for figure 3 and performs associated
%analysis.

%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
%% -----set up------
delta=[31]
folder=['../data/flow_phantom/a/UIV/2cm/'];
name=['/phantom_a_128lines_' int2str(delta) 'delta_20000depth7000000freq_out']
load([folder int2str(delta) name '.mat']);
depth=0.02;
ypix=2.3653e-05;
lines=128;
[VxC, VyC]=pulsePIVScaleAndCorrect_f2(delta, depth, Vx, Vy, 1,'glyce',period, split);
if delta==7
    indices=[1:50];
    centre=[1.9,1.1];
elseif delta==15
    indices=[47:50 1:46];
    centre=[1.9,1.05];
elseif delta==31
    indices=[4:50 1:3]; %indices=[5:50 1:4];
    centre=[1.9,1.05];
elseif delta==47
     indices=[25:50 1:24];
    centre=[1.9,1.05];
elseif delta==63
    indices=[1:50];
    centre=[1.9,1.05];
end

%--sort out the time steps--
depth=0.02;
lines=128;
timestep=(lines*2+(delta-1))*(2*depth/1540+23.488e-6); %in s
totalTime=timestep*period;
times=totalTime/split:totalTime/split:totalTime;

%---filter xy---
VxF=zeros(size(VxC));
VyF=zeros(size(VyC));
for i=1:50
    VxF(:,:,i)=medfilt2(VxC(:,:,i),[27 3]);
    VyF(:,:,i)=medfilt2(VyC(:,:,i),[27 3]);
end

VxC=VxF;
VyC=VyF;
clear VxF VyF;

%% ----------------4 vector plots-------------------------
yfrom=100; yto=180;
for t=[1:50]
    VxN(:,:,t)=VxC(:,yfrom:yto,t);
    VyN(:,:,t)=VyC(:,yfrom:yto,t);
end
[a,b,c]=size(VxN);
vmagC=(VxN.^2+VyN.^2).^0.5;
vmagMax=max(vmagC(:));
VxN(a,b,:)=vmagMax;
VyN(a,b,:)=0;
xgxN=xgx(:,yfrom:yto);
ygyN=ygy(:,yfrom:yto);

for t=indices([8 13 23 33]);
    figure;
    
    subplot(1,2,2);
    colours2=(([0:vmagMax/100:vmagMax])./(2*vmagMax)+0.51)*128;
    velocitiesC=0:vmagMax/100:vmagMax;
    w=1;
    g=image(w,velocitiesC,colours2');
    set(gca,'Position',[0.95 0.17 0.04 0.75]);
    set(gca,'FontSize',6);   
    set(gca,'XTick',[]);
    xlabel({'velocity      ','/ms^{-1}      '},'FontSize',6);
    
    subplot(1,2,1);
    colours=[gray(64); jet(64)];
    colormap(colours);
    c=round(yfrom/length(VxC)*971)+1;
    d=round(yto/length(VxC)*971)-1
    i=log(1+abs(hilbert(Iavg(:,c:d,t))))';
    i=64*i./max(i(:));
    [a,b]=size(i);
    j=image([1:b].*xpix*1000,[c:d].*ypix*1000,i);
    set(gca,'Position',[0.14 0.17 0.7 0.75]);
    hold on;
    h=quiverc4ePIV(xgxN(4:8:end,4:4:end)*xpix*1000,ygyN(4:8:end,4:4:end).*ypix*1000,VxN(4:8:end,4:4:end,t),VyN(4:8:end,4:4:end,t),2,1);
    set(gca,'FontSize',8);
    title(['time = ' num2str(t/50.*times(end),2) ' s'],'FontSize',8);
    xlabel('axial position /mm','FontSize',8);
    ylabel('depth /mm','FontSize',8);
    hold off;
    m=m+1;
    set(gcf, 'PaperUnits', 'centimeters');
    set(gcf, 'PaperSize', [6.3 5]);
    set(gcf, 'PaperPosition', [0 0 6.3 5]);
    figfilename=['phantom_a_Delta' int2str(delta) '_vector_plots_' int2str(t) '.png']
    print(gcf, '-painters', '-dpng', '-r600', figfilename);
end
a=1;

%% -----velocity profile---------
for t=indices
    vprofile(:,t)=mean(-VxC([1:22 36:63],yfrom:yto,t),1);
end
figure;
colors=colormap(jet(50));
set(gca,'Color','none','Position',[0.27 0.08 0.7 0.15]);
r=([1:length(vprofile(:,t))]).*iastep2.*ypix.*-1;
r=(r-r(end)/2)*1000;
c=1;
hold on; plot(r,vprofile(:,indices(1)),'Color',colors(c,:),'LineWidth',0.5);
set(gca,'FontSize',8);
xlabel('radial position /mm','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
xlim([-3 3]);
ylim([-0.5 1.5]);
idx=[2:2:50];
for t=indices(idx)
    pos=get(gca,'position');
    pos=pos+[0 +0.0350 0 0];
    axes('position',pos,'Color', 'none');
    hold on; plot(r,vprofile(:,t),'Color',colors(idx(c),:),'LineWidth',0.5);
    xlim([-3 3]);
    ylim([-0.5 1.5]);
    axis off;
    c=c+1;
end
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [5 16]);
set(gcf, 'PaperPosition', [0 0 4 16]);
figfilename=['phantom_a_Delta' int2str(delta) '_vel_profiles.png'];
print( '-painters', '-dpng', '-r600', figfilename);

%% ----flow rate-comparison----
% ----flow rate-transonics----
M=dlmread('../data/flow_phantom/a/flow_rate/phantom_a_flow.csv',',',[4 0 936 14]);
timeT=M(:,1);
timeT=timeT./1000;
flowT=M(:,12);
flowTMax=M(:,14);
flowTMin=M(:,15);

%scale by bucket flow rate to eliminate speed of sound error with
%Transonics
scaleFactor=383.8./mean(flowT);
flowTMax=flowTMax.*scaleFactor;
flowTMin=flowTMin.*scaleFactor;
flowT=flowT.*scaleFactor;

%---flow rate UIV---
[flow2, Flow,upperEdge, lowerEdge]=flowRatePulse_f(xpix, ypix, times, 50, depth, lines, VxC, VyC, Iavg, iastep1, iastep2, 390, 440, 670, 730);
flow=flow2(indices);
figure;plot(times,flow,'r','LineWidth',0.1);
set(gca,'Position',[0.25 0.18 0.73 0.78]);
Tindices=[30:length(flowT) 1:29];
hold on; plot(timeT, flowT(Tindices),'k','LineWidth',0.1);
hold on; plot(timeT,flowTMax(Tindices),'k:','LineWidth',0.1);
hold on; plot(timeT,flowTMin(Tindices),'k:','LineWidth',0.1);
strings={'UIV',['flow meter' 10 'mean+/-std']};
legend(strings,'Position',[0.66 0.8 0.2 0.1],'FontSize',6);
legend('boxoff');
xlabel('time /s','FontSize',8);
ylabel('flow rate /ml min^{-1}','FontSize',8);
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [6.3 5]);
set(gcf, 'PaperPosition', [0 0 6.3 5]);
figfilename=['phantom_a_Delta' int2str(delta) '_flow_rate.png'];
print( '-painters', '-dpng', '-r600', figfilename);

flowT2 = interp1([timeT], [flowT(Tindices)],times(1:end-1));
hold on; plot(times(1:end-1), flowT2,'b');
rmsdifference=(mean((flow(1:end-1)-flowT2').^2))^0.5;
percenterror=mean(abs((flow(1:end-1)-flowT2')./flowT2'))*100;
meanFlowerror=mean(abs(flow(1:end-1)-flowT2'))
percentFlowerror=meanFlowerror./(max(flow)).*100
%% ----womersley profile----

for t=1:50
    R(t)=ypix*mean(lowerEdge(:,t)-upperEdge(:,t))/2;
end
Rmean=mean(R);
flow_f=flowT(Tindices);
pos=round(length(timeT)/50:length(timeT)/50:length(timeT));
nu=3.5e-6;

vmean=flowT(Tindices)./(pi*Rmean^2);
vmean2=vmean./(60*1000*1000);
pos=round(length(timeT)/50:length(timeT)/50:length(timeT));
vel3= womer_Vmean(timeT(pos), 20, vmean2(pos), Rmean,nu);

vprofile_average=(vprofile(round(length(vprofile)/2):end,:)+vprofile(fliplr(1:round(length(vprofile)/2)),:))/2;

figure;  plot([1:50], real(vel3(:,1)),'k');
hold on; plot([1:50], vprofile(round(length(vprofile)/2), indices),'r');
legend('Womersley at R=0', 'UIV at R=0');
figfilename=['phantom_a_Delta' int2str(delta) '_WomersleyComparison_R=0b.png'];
print( '-painters', '-dpng', '-r600', figfilename);

womerpro=real(vel3);
UIVpro=interp1(r(fliplr(1:round(length(vprofile)/2))),vprofile_average(:,indices),[0:0.01:1].*R(indices(t))*1000);
rmsdifference=(mean(mean((womerpro-UIVpro').^2)))^0.5;
percenterror=mean(mean(abs((womerpro-UIVpro')./womerpro)))*100;
meanWomerror=mean(mean(abs(womerpro-UIVpro')))
percenterror=meanWomerror./(max(max(womerpro))).*100
percenterrorAxis=abs(womerpro(:,1)-UIVpro(1,:)')./abs(womerpro(:,1));

%% -----velocity profile comparison with Womersley---------

figure;
colors=colormap(jet(50));
set(gca,'Color','none','Position',[0.27 0.08 0.7 0.2]);
r=([1:length(vprofile(:,1))]).*iastep2.*ypix.*-1;
r=(r-r(end)/2)*1000;
c=1;

hold on; plot([0:0.01:1].*Rmean*1000, real(vel3(2,:)),'k','LineWidth',0.5);
hold on; plot(r(fliplr(1:round(length(vprofile)/2))),vprofile_average(:,indices(2)),'Color',colors(c,:),'LineWidth',1.5);

set(gca,'FontSize',8);
xlabel('radial position /mm','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
xlim([0 Rmean*1000]);
ylim([-0.5 1.5]);
idx=[4:2:50];
for t=indices(idx)
    pos=get(gca,'position');
    pos=pos+[0 +0.0350 0 0];
    axes('position',pos,'Color', 'none');
    hold on; plot([0:0.01:1].*Rmean*1000, real(vel3(t,:)),'k','LineWidth',0.5);
    hold on; plot(r(fliplr(1:round(length(vprofile)/2))),vprofile_average(:,indices(t)),'Color',colors(idx(c),:),'LineWidth',1.5);
    xlim([0 Rmean*1000]);
    ylim([-0.5 1.5]);
    axis off;
    c=c+1;
end
set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [5 16]);
set(gcf, 'PaperPosition', [0 0 4 16]);
figfilename=['phantom_a_Delta' int2str(delta) '_vel_profiles.png'];
print( '-painters', '-dpng', '-r600', figfilename);

%% --simulate the Doppler spectrum--
%extract a box of data
figure;
depth=-0.1:0.1:2;
width=0:0.1:3.8;
j=imagesc(width,depth,log(1+abs(hilbert(squeeze(Iavg(iastep1:iastep1:end,iastep2:iastep2:end,1)))))');
axis equal;  axis tight;
heightBox=0.2;
widthBox=3;
[x y]=rectROI(widthBox,heightBox,centre,0);
i = impoly(gca, [x; y]');
BW=createMask(i,j);
velocities=-50:150;
intensity=zeros(length(velocities),50);
%simulate the spectrum within that box
T=1;
for t=[1:50]
    vmag=(VxC(:,:,t)'.^2+VyC(:,:,t)'.^2).^0.5;
    vmag=-1.*vmag.*sign(VxC(:,:,t)');
    pos=find([BW(1:end-1,1:end-1)]);
    vmag_small=vmag(pos);
    pos=find(isfinite(vmag_small));
    vmag_smaller=vmag_small(pos);
    for n=1:length(vmag_smaller)
        v=round(vmag_smaller(n)*100);
        pos=find(velocities==v);
        intensity(pos,T)=intensity(pos,T)+1;
    end
    averageUIV(t)=mean(vmag_smaller);
    T=T+1;
end
intensity=flipud(intensity);

%overlay spectra
vars={'average', 'ourVpix', 'unwrap', 'time', 'scaledMaxV', 'scaledMinV'};
load('../data/flow_phantom/a/Doppler/a_Doppler.mat', vars{:});
[v_long,t_long]=size(average);
velocitiesDS=([0:v_long]-unwrap).*ourVpix.*0.01;
timesDS=time;

%spectra
fig=figure;
subplot(2,1,1); imagesc(timesDS, velocitiesDS, average); colormap gray;
intensity2= ind2rgb(flipud(intensity),jet);
hold on; g=imagesc(times, velocities./100,intensity2(:,indices,:));
set(g,'AlphaData',flipud(intensity(:,indices)).^0.5);
set(g,'AlphaDataMapping','scaled');
set(gca,'YDir','normal');
f=get(gca,'position');
f2=[f(1)*1.45 f(2)*0.8 f(3) f(4)*1.1];
set(gca, 'position', f2,'FontSize',8,'FontWeight','normal');
hold off;
ylim([-0.5 1.5]);
xlim([0 0.96]);
xlabel('time /s','FontSize',8);
ylabel('velocity /m s^{-1}','FontSize',8);
ax=gca;
ax.TitleFontSizeMultiplier = 1;
title({'Doppler spectrum (grayscale)','overlayed with UIV spectrum (colour)'},'FontSize',8,'FontWeight','normal');

%locator DS
locator = imread('../data/flow_phantom/a/Doppler/a_Doppler_location.png');
subplot(2,2,3); imagesc(locator);
axis equal; axis tight; axis off;
h=get(gca,'position');
h2=[0.35 0.01 h(3)*0.6 h(4)*0.6];
ax=gca;
ax.TitleFontSizeMultiplier = 1;
set(ax, 'position', h2);

%locator UIV
subplot(2,2,4);
i=log(1+abs(hilbert(Iavg(:,:,t))))';
i=64*i./max(i(:));
RGB = ind2rgb8(i,gray);
depth=-0.1:0.1:2;
width=0:0.1:3.8;
j=image(width,depth,RGB);
axis equal;  axis tight; axis off;
hold on;
ylim([0 2]);
h=get(gca,'position');
h2=[0.7 0.05 h(3)*0.6 h(4)*0.6];
set(gca, 'position', h2);
plot(x(1:2),y(1:2),'c');
plot(x(2:3),y(2:3),'c');
plot(x(3:4),y(3:4),'c');
plot([x(4) x(1)],[y(4) y(1)],'c');
ax=gca;
ax.TitleFontSizeMultiplier = 1;
sometext={'sample volumes:'};
ax1 = axes('Position',[0 0 1 1],'Visible','off');
text(.025,0.25,sometext,'FontSize',8);
sometext={'Doppler'};
text(.15,0.15,sometext,'FontSize',8);
sometext={'UIV'};
text(.6,0.15,sometext,'FontSize',8);

set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [6.3 5]);
set(gcf, 'PaperPosition', [0 0 6.3 5]);
figfilename=['phantom_a_Delta' int2str(delta) '_Doppler_comp.png'];
print( '-painters', '-dpng', '-r600', figfilename);

%error
background=mean(average(1,:));
sumAverageDS=(velocitiesDS(1:end-1)'*ones(1,length(timesDS))).*(average-background);
weightedAverageDS=sum(sumAverageDS)./sum(average-background);
weightedAverageDS2=interp1(timesDS,weightedAverageDS,times);
meanDiffaveDS=mean(abs(averageUIV(indices)-weightedAverageDS2))
percentDiffaveDS=meanDiffaveDS/max(averageUIV).*100;

peak=max(average);
for t=1:length(average)
    threshDS=find(average(:,t)>(peak(t)/2));
    maxDS(t)=velocitiesDS(threshDS(end));
end
maxDS2=interp1(timesDS,maxDS,times);
meanDiffmaxDS=mean(abs(averageUIV(indices)-maxDS2))
percentDiffmaxDS=meanDiffmaxDS/max(averageUIV).*100


%% ------write out errors file
filename=['phantom_a_Delta' int2str(delta) 'errors.txt'];

fid=fopen(filename,'w');
formatSpec = 'mean velocity error compared to ave Doppler is %5.3f m/s or %5.3f \n';
fprintf(fid,formatSpec,meanDiffaveDS,percentDiffaveDS);
formatSpec = 'mean velocity error compared to max Doppler is %5.3f m/s or %5.3f \n';
fprintf(fid,formatSpec,meanDiffmaxDS,percentDiffmaxDS);
formatSpec = 'mean flow rate error compared to transit time is %5.3f ml/min or %5.3f \n';
fprintf(fid,formatSpec,meanFlowerror,percentFlowerror);
formatSpec = 'mean velocity error compared to Womersley is %5.3f m/s or %5.3f \n';
fprintf(fid,formatSpec,meanWomerror,percenterror);
fclose('all');