% Copyright (C) 2019 Emanuele Ortu, Eleonora Sulas, Danilo Pani.
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
% To contact the authors:
% sulaseleonora1992@gmail.com, danilo.pani@unica.it,
% Address: Department of Electrical and Electronic Engineering,
% University of Cagliari,Via Marengo 3 - 09123 Cagliari, Italy
% Cagliari, 15 Jan 2020
% Please cite the following article when using this code:
% TBA
% This matlab code extract the upper and the lower envelope of the Doppler
% image.
% input:
% filename: the filename as a string and bmp format
% output:
% x_up: upper envelope
% x_down: lower envelope
function [x_up,x_down]=envelope_extraction(filename)
if nargin==0
[file,path] = uigetfile({'*.bmp'});
filename=fullfile(path,file);
end
% loading
x = imread(filename);
x = rgb2gray(x);
% delete x axis line(from white to black)
ascissa = find(sum(x')==255*length(x));
x(ascissa,:)= 0;
% median otsu 2D threshold
T = median_otsu2d(x,5);
BW = imbinarize(x,T/255);
BW = bwareaopen(BW,70,4);
x1 = BW;
x_down = [];
x_up = [];
up = x1(1:ascissa-1,:);
down = x1(ascissa+1:end,:);
for j = 1:length(x1)
ind = find(down(:,j)==1);
if isempty(ind) == 1
x_down(j) = 0;
else
x_down(j)= ind(end);
end
ind = find(up(:,j)==1);
if isempty(ind) == 1
x_up(j) = ascissa;
else
x_up(j)= ind(1);
end
end
x_down = - x_down;
x_up = ascissa - (x_up);
end
function th = median_otsu2d(img1,k)
dim = size(img1);
if length(dim)==3
img1 = rgb2gray(img1);
end
k_2=fix(k/2);
[m,n]=size(img1);
f=img1(1+k_2:m-k_2,1+k_2:n-k_2);
g=zeros(m-k_2,n-k_2);
f=double(f);
c=zeros(256,256);
for i=1+k_2:m-k_2
for j=1+k_2:n-k_2
g(i-k_2,j-k_2)=median(median((img1(i-k_2:i+k_2,j-k_2:j+k_2))));
c(f(i-k_2,j-k_2)+1,g(i-k_2,j-k_2)+1)=c(f(i-k_2,j-k_2)+1,g(i-k_2,j-k_2)+1)+1;
end
end
[m1,n1]=size(f);
p = c./(m1*n1);
figure
surf(p)
for s=0:255
for t=0:255
w0=sum(sum(p(1:s,1:t)));
w1=sum(sum(p(s+1:end,t+1:end)));
u0=[sum(sum(p(1:s,1:t)'*(1:s)'))/w0, sum(sum(p(1:s,1:t)*(1:t)'))/w0 ];
u1=[sum(sum(p(s+1:end,t+1:end)'*(s+1:256)'))/w1, sum(sum(p(s+1:end,t+1:end)*(t+1:256)'))/w1];
uT=[sum(sum(p(1:end,1:end)'*(1:256)')), sum(sum(p(1:end,1:end)*(1:256)'))];
sigma=(w0*((u0-uT)*(u0-uT)'))+(w1*((u1-uT)*(u1-uT)'));
trace(s+1,t+1)=(w0*(((u0(1)-uT(1))^2)+(u0(2)-uT(2))^2))+(w1*(((u1(1)-uT(1))^2)+(u1(2)-uT(2))^2));
end
end
[mas,ind] = max(trace);
[mas2, ind2] = max(mas);
th = ind2;
end