1. AMDEMD decomposition process
fenjieguocheng.pdf
2. Code of AMDEMD
function [IMF, R]= AMDEMD(x,n,alpha,beta)
%INPUT PARAMETERS:
%x - input image
%n - the number of decomposition
%alpha - The filter size adjustment coefficient
%beta - the threshold value of the inner stopping criterion
%denoise - set as 'bm3d' to perform BM3D denoising (recommended)
[hg,wd]=size(x);
imf=[];
img_lin=x;
count=100;
for cishu=1:n
img_r=img_lin;
sift=1;
max_avg=[];
different_max_avg=[];
while sift<=count
%compute the size of window
I=img_lin;
se=strel('disk',2);
Ridge= imopen(I, se);
Ridge_tophat= I- Ridge;
BW_Ridge_tophat=abs(Ridge_tophat)>eps;
BW_Ridge_roughthin_tophat= bwmorph( BW_Ridge_tophat, 'thin', Inf);
BW_Ridge_roughthin_tophat= bwmorph( BW_Ridge_roughthin_tophat, 'bridge');
BW_EDT_Ridge=bwdist(BW_Ridge_roughthin_tophat);
msk = strel('square',3);
MaxDil_Ridge = imdilate(BW_EDT_Ridge,msk); %MAX filtration
M_n_Ridge = ~(BW_EDT_Ridge - MaxDil_Ridge); %binary map of maxima
% Map of distances between extrema
DBE_n_Ridge = 2*BW_EDT_Ridge.*M_n_Ridge;
nonZeroDBE_n_Ridge= DBE_n_Ridge(DBE_n_Ridge ~= 0);
meanDBE_n_Ridge = mean(nonZeroDBE_n_Ridge);
Valley= imclose(I, se);
Valley_tophat= Valley-I;
BW_Valley_tophat=abs(Valley_tophat)>eps;
BW_Valley_roughthin_tophat= bwmorph( BW_Valley_tophat, 'thin', Inf);
BW_Valley_roughthin_tophat= bwmorph( BW_Valley_roughthin_tophat, 'bridge');
BW_EDT_Valley=bwdist(BW_Valley_roughthin_tophat);
msk = strel('square',3);
MaxDil_Valley = imdilate(BW_EDT_Valley,msk); %MAX filtration
M_n_Valley = ~(BW_EDT_Valley - MaxDil_Valley); %binary map of maxima
% Map of distances between extrema
DBE_n_Valley= 2*BW_EDT_Valley.*M_n_Valley;
nonZeroDBE_n_Valley= DBE_n_Valley(DBE_n_Valley ~= 0);
meanDBE_n_Valley = mean(nonZeroDBE_n_Valley);
win=round(beta*(meanDBE_n_Ridge+meanDBE_n_Valley)/2);
%%
if win<=3
win=5;
end
if mod(win,2)==0
win=win+1;
end
bweightfilter=eightdirection2(win);
linshi=[];
for i=1:hg
for j=1:wd
linshi=img_lin(max(1,i-(win-1)/2):min(i+(win-1)/2,hg),max(1,j-(win-1)/2):min(j+(win-1)/2,wd));
cankaofangxiang(i,j)=CalculateOrientation2(linshi);
end
end
upmax =[];
downmin=[];
upmaxlin =[];
downminlin=[];
%%
%find the Maximum and minimum
m=60;
for i=1:8
upmaxlin(:,:,i)=xingtaimax(img_lin,bweightfilter(:,:,i),m);
downminlin(:,:,i)=xingtaimin(img_lin,bweightfilter(:,:,i),m);
end
%
for i=1:hg
for j=1:wd
upmax(i,j)=upmaxlin(i,j,cankaofangxiang(i,j));
downmin(i,j)=downminlin(i,j,cankaofangxiang(i,j));
end
end
lin_mean=(upmax+downmin)/2;
%%
% smoothing
for i=1:hg
for j=1:wd
linshi1=lin_mean(max(1,i-(win-1)/2):min(i+(win-1)/2,hg),max(1,j-(win-1)/2):min(j+(win-1)/2,wd));
img_avg(i,j)=mean(mean(linshi1));
end
end
img_lin=img_lin-img_avg;
%%
% inner stop
max_avg(sift)=max(max(abs(img_avg)));
if sift >= 2
different_max_avg(sift-1)=abs(max_avg(sift)-max_avg(sift-1));
mean_different=mean(different_max_avg);
mean_different_before=mean(different_max_avg(1:end-1));
ratio_mean_different=abs(mean_different-mean_different_before)/mean_different_before;
if ratio_mean_different < alpha
break;
end
end
sift=sift+1;
end
imf(:,:,cishu)=img_lin;
img_lin=img_r-img_lin;
%%
%outer stop
f=imf(:,:,cishu);
[r,c]=size(f);
matrix_mean=mean2(f);
matrix_std=std2(f);
T=matrix_mean+matrix_std;
magnitude_points=zeros(r,c);
magnitude_points(find(abs(f)>=T))=1;
target_array=abs(f).*magnitude_points;
magnitude_sum=sum(target_array(:));
streak_rate(cishu)=magnitude_sum/(c*r);
if cishu>=3
if streak_rate(cishu-2)>streak_rate(cishu-1) && streak_rate(cishu-1)<streak_rate(cishu)
break;
end
end
end
IMF=imf(:,:,1:cishu-2);
R=imf(:,:,cishu-1)+imf(:,:,cishu)+img_lin;
end
%%
function bbwfilter=eightdirection2(win)
u1=0;
u2=0;
N=(win-1)/2;
sigma1=win*win;
sigma2=N/2;
p=0;
[tpx,tpy]=meshgrid(-N:N,-N:N);
for s=1:8
x=tpx*cos((s-1)*pi/8)-tpy*sin((s-1)*pi/8);
y=tpx*sin((s-1)*pi/8)+tpy*cos((s-1)*pi/8);
f=1/(2*pi*sigma1*sigma2*sqrt(1-p*p))*exp(-1/(2*(1-p*p))*(((x-u1).^2)./(sigma1*sigma1)-2*p*((x-u1)*(y-u2))./(sigma1*sigma2)+((y-u2).^2)./(sigma2*sigma2)));
w=zeros(win,win);
l=find(f>=mean(mean((f))));
w(l)=1;
bbwfilter(:,:,s)=w;
end
end
%%
function eightorientation = CalculateOrientation2(x)
[delx, dely] = gradient(x);
Vx = sum(sum(2*delx.*dely))+eps;
Vy = sum(sum((delx.^2)-(dely.^2))); %I'm sure you can work it out
orientation = atan(Vy/Vx)/pi*180+90;
if (orientation>=0 && orientation<11.25) || (orientation>=168.75 && orientation<=180)
eightorientation=1;
elseif orientation>=11.25 && orientation<33.75
eightorientation=2;
elseif orientation>=33.75 && orientation<56.25
eightorientation=3;
elseif orientation>=56.25 && orientation<78.75
eightorientation=4;
elseif orientation>=78.75 && orientation<101.25
eightorientation=5;
elseif orientation>=101.25 && orientation<123.75
eightorientation=6;
elseif orientation>=123.75 && orientation<146.25
eightorientation=7;
elseif orientation>=146.25 && orientation<168.75
eightorientation=8;
end
end
%%
function getmax=xingtaimax(im,bw,m)
[l,q]=size(bw);
[ll,qq]=size(im);
ime=exp(im*m);
last=convolve2(ime,bw,'same');
getmax=(1/(m))*log(last);
end
%%
function getmin=xingtaimin(im,bw,m)
[l,q]=size(bw);
[ll,qq]=size(im);
ime=exp(-1*im*m);
last=convolve2(ime,bw,'same');
getmin=-(1/(m))*log(last);
end
%%
function y = convolve2(x, m, shape, tol)
%CONVOLVE2 Two dimensional convolution.
% Y = CONVOLVE2(X, M) performs the 2-D convolution of matrices X and
% M. If [mx,nx] = size(X) and [mm,nm] = size(M), then size(Y) =
% [mx+mm-1,nx+nm-1]. Values near the boundaries of the output array are
% calculated as if X was surrounded by a border of zero values.
%
% Y = CONVOLVE2(X, M, SHAPE) where SHAPE is a string returns a
% subsection of the 2-D convolution with size specified by SHAPE:
%
% 'full' - (default) returns the full 2-D convolution
%
% 'valid' - returns only those parts of the convolution
% that can be computed without padding; size(Y) =
% [mx-mm+1,nx-nm+1] when size(X) > size(M)
%
% 'same' - returns the central part of the convolution
% that is the same size as X using zero padding
%
% 'wrap' or
% 'circular' - as for 'same' except that instead of using
% zero-padding the input X is taken to wrap round as
% on a toroid
%
% 'reflect' or
% 'symmetric' - as for 'same' except that instead of using
% zero-padding the input X is taken to be reflected at
% its boundaries
%
% 'replicate' - as for 'same' except that instead of using
% zero-padding the rows at the array boundary are
% replicated
%
% CONVOLVE2 is fastest when mx > mm and nx > nm - i.e. the first
% argument is the input and the second is the mask.
%
% If the rank of the mask M is low, CONVOLVE2 will decompose it into a
% sum of outer product masks, each of which is applied efficiently as
% convolution with a row vector and a column vector, by calling CONV2.
% The function will often be faster than CONV2 or FILTER2 (in some
% cases much faster) and will produce the same results as CONV2 to
% within a small tolerance.
%
% Y = CONVOLVE2(... , TOL) where TOL is a number in the range 0.0 to
% 1.0 computes the convolution using a reduced-rank approximation to
% M, provided this will speed up the computation. TOL limits the
% relative sum-squared error in the effective mask; that is, if the
% effective mask is E, the error is controlled such that
%
% sum(sum( (M-E) .* (M-E) ))
% -------------------------- <= TOL
% sum(sum( M .* M ))
%
% See also CONV2, FILTER2, EXINDEX
% Copyright David Young, Feb 2002, revised Jan 2005, Jan 2009, Apr 2011,
% Feb 2014
% Deal with optional arguments
narginchk(2,4);
if nargin < 3
shape = 'full'; % shape default as for CONV2
tol = 0;
elseif nargin < 4
if isnumeric(shape)
tol = shape;
shape = 'full';
else
tol = 0;
end
end;
% Set up to do the wrap & reflect operations, not handled by conv2
if ismember(shape, {'wrap' 'circular' 'reflect' 'symmetric' 'replicate'})
x = extendarr(x, m, shape);
shape = 'valid';
end
% do the convolution itself
y = doconv(x, m, shape, tol);
end
%-----------------------------------------------------------------------
function y = doconv(x, m, shape, tol)
% Carry out convolution
[mx, nx] = size(x);
[mm, nm] = size(m);
% If the mask is bigger than the input, or it is 1-D already,
% just let CONV2 handle it.
if mm > mx || nm > nx || mm == 1 || nm == 1
y = conv2(x, m, shape);
else
% Get svd of mask
if mm < nm; m = m'; end % svd(..,0) wants m > n
[u,s,v] = svd(m, 0);
s = diag(s);
rank = trank(m, s, tol);
if rank*(mm+nm) < mm*nm % take advantage of low rank
if mm < nm; t = u; u = v; v = t; end % reverse earlier transpose
vp = v';
% For some reason, CONV2(H,C,X) is very slow, so use the normal call
y = conv2(conv2(x, u(:,1)*s(1), shape), vp(1,:), shape);
for r = 2:rank
y = y + conv2(conv2(x, u(:,r)*s(r), shape), vp(r,:), shape);
end
else
if mm < nm; m = m'; end % reverse earlier transpose
y = conv2(x, m, shape);
end
end
end
%-----------------------------------------------------------------------
function r = trank(m, s, tol)
% Approximate rank function - returns rank of matrix that fits given
% matrix to within given relative rms error. Expects original matrix
% and vector of singular values.
if tol < 0 || tol > 1
error('Tolerance must be in range 0 to 1');
end
if tol == 0 % return estimate of actual rank
tol = length(m) * max(s) * eps;
r = sum(s > tol);
else
ss = s .* s;
t = (1 - tol) * sum(ss);
r = 0;
sm = 0;
while sm < t
r = r + 1;
sm = sm + ss(r);
end
end
end
%-----------------------------------------------------------------------
function y = extendarr(x, m, shape)
% Extend x so as to wrap around on both axes, sufficient to allow a
% "valid" convolution with m to return a result the same size as x.
% We assume mask origin near centre of mask for compatibility with
% "same" option.
[mx, nx] = size(x);
[mm, nm] = size(m);
mo = floor((1+mm)/2); no = floor((1+nm)/2); % reflected mask origin
ml = mo-1; nl = no-1; % mask left/above origin
mr = mm-mo; nr = nm-no; % mask right/below origin
% deal with shape option terminology - was inconsistent with exindex
switch shape
case 'wrap'
shape = 'circular';
case 'reflect'
shape = 'symmetric';
end
y = exindex(x, 1-ml:mx+mr, 1-nl:nx+nr, shape);
end