%This script generates the images for figure 9 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]
name=['../data/in_vivo/UIV/aorta_128lines_31delta_20000depth7000000freq']
load([name '.mat']);
indices=[13:25 1:12];
depth=0.02;
ypix=2.3653e-05;
lines=128;
[VxC, VyC]=pulsePIVScaleAndCorrect_f2(delta, depth, Vx, Vy, 1,'glyce',period, split);

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


%---filter t---
%use for high pulse rate
VxC2(:,:,1)=VxC(:,:,25);
VxC2(:,:,2:26)=VxC(:,:,1:25);
VxC2(:,:,27)=VxC(:,:,1);
VyC2(:,:,1)=VyC(:,:,25);
VyC2(:,:,2:26)=VyC(:,:,1:25);
VyC2(:,:,27)=VyC(:,:,1);
[a,b,c]=size(VxC2);
for l=1:a
    for r=1:b
        VxC3(l,r,:)=medfilt1(squeeze(VxC2(l,r,:)),[3]);
        VyC3(l,r,:)=medfilt1(squeeze(VyC2(l,r,:)),[3]);
    end
end
VxC=VxC3(:,:,2:26);
VyC=VyC3(:,:,2:26);

%% ----------------4 vector plots-------------------------
yfrom=1; yto=259;
for t=[1:25]
    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(1,1,:)=vmagMax;
VyN(1,1,:)=0;
xgxN=xgx(:,yfrom:yto);
ygyN=ygy(:,yfrom:yto);

idx=[3 5 7 9 20]
for n=[1:5];
    t=indices(idx(n));
    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]-70).*ypix*1000,i);
    set(gca,'Position',[0.14 0.17 0.7 0.75]);
    hold on;

    h=quiverc4ePIV(xgxN(1:4:end,1:4:end)*xpix*1000,(ygyN(1:4:end,1:4:end)-70).*ypix*1000,VxN(1:4:end,1:4:end,t),VyN(1:4:end,1:4:end,t),2,1);
    set(gca,'FontSize',8);
    title(['time = ' num2str(idx(n)/25.*times(end),2) ' s'],'FontSize',8);
    xlabel('axial position /mm','FontSize',8);
    ylabel('depth /mm','FontSize',8);

    ylim([0 20]);
    xlim([0 38]);
    hold off;
    m=m+1;
    set(gcf, 'PaperUnits', 'centimeters');
    set(gcf, 'PaperSize', [8.4 5]);
    set(gcf, 'PaperPosition', [0 0 8.4 5]);
    figfilename=['aorta_Delta' int2str(delta) '_vector_plots_' int2str(idx(n)) '.png']
    print(gcf, '-painters', '-dpng', '-r600', figfilename);
end
a=1;
%% ----------------4 vector plots-zoom in-------------------------
yfrom=1; yto=259;
for t=[1:25]
    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(1,1,:)=vmagMax;
VyN(1,1,:)=0;
xgxN=xgx(:,yfrom:yto);
ygyN=ygy(:,yfrom:yto);
idx=[3 5 7 9 20]
for n=[1:4];
    t=indices(idx(n));
    figure;
    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]-70).*ypix*1000,i);
    set(gca,'Position',[0.2 0.17 0.8 0.75]);
    hold on;
    h=quiverc4ePIV(xgxN(1:4:end, 1:2:end)*xpix*1000,(ygyN(1:4:end, 1:2:end)-70).*ypix*1000,VxN(1:4:end, 1:2:end,t),VyN(1:4:end, 1:2:end,t),2,1);
    set(gca,'FontSize',8);
    title(['time = ' num2str(idx(n)/25.*times(end),2) ' s'],'FontSize',8);
    xlabel('axial position /mm','FontSize',8);
    ylabel('depth /mm','FontSize',8);
    ylim([2 11]);
    xlim([16 32]);
    hold off;
    m=m+1;
    set(gcf, 'PaperUnits', 'centimeters');
    set(gcf, 'PaperSize', [5.5 5]);
    set(gcf, 'PaperPosition', [0 0 5.5 5]);
    figfilename=['aorta_Delta' int2str(delta) '_vector_plots_zoom_' int2str(idx(n)) '.png']
    print(gcf, '-painters', '-dpng', '-r600', figfilename);
end

%% --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.12;
widthBox=0.2;
centre=[2.32,0.54];
[x y]=rectROI(widthBox,heightBox,centre,15);
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:25]
    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', 'alias_sep', 'time'};
load('../data/in_vivo/Doppler/aorta_Doppler.mat', vars{:});


[v_long,t_long]=size(average);
velocitiesDS=([0:127]-(64-alias_sep)).*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',2*flipud(intensity(:,indices)).^0.5);
set(g,'AlphaDataMapping','scaled');
set(gca,'YDir','normal');
f=get(gca,'position');
f2=[f(1)*1.05 f(2)*0.78 f(3)*1.1 f(4)*1.1];
set(gca, 'position', f2,'FontSize',8,'FontWeight','normal');
hold off;
ylim([-0.1 1]);
xlim([0 0.35]);
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/in_vivo/Doppler/locator.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);
%title('Doppler','FontSize',8,'FontWeight','normal');

%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.03 h(3)*0.6 h(4)*0.6];
set(gca, 'position', h2);
plot(x(1:2),y(1:2),'c','linewidth',0.1);
plot(x(2:3),y(2:3),'c','linewidth',0.1);
plot(x(3:4),y(3:4),'c','linewidth',0.1);
plot([x(4) x(1)],[y(4) y(1)],'c','linewidth',0.1);
ax=gca;
ax.TitleFontSizeMultiplier = 1;
%title('UIV','FontSize',8,'FontWeight','normal');
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', [8.4 5]);
set(gcf, 'PaperPosition', [0 0 8.4 5]);
figfilename=['aorta_Delta' int2str(delta) '_Doppler_comp.png'];
print( '-painters', '-dpng', '-r600', figfilename);

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

peak=max(average);
peak=background*1.5;
[a,b]=size(average);
for t=1:b
    threshDS=find(average(:,t)>(peak));%(t)/2
    maxDS(t)=velocitiesDS(threshDS(end));
end
maxDS2=interp1(timesDS,maxDS,times);
meanDiffmaxDS=mean(abs(averageUIV(indices(1:end-1))-maxDS2(1:end-1)))
percentDiffmaxDS=meanDiffmaxDS/max(averageUIV).*100

figure; plot(averageUIV(indices))
hold on; plot(maxDS2)