function  [ix,resx,medel]=IEMD_public(image,epsilon,numberimfs,conn)
% Decompose an image into its intrinsic modes IMFs and corresponding
% residues.
%  
%
% INPUT:
% image- a matlab 2D vector;
% epsilon- value of stop criterion, between 0 and 1 but you can try other
% and see what happens
% numberimfs- maximum number of IMF
% conn- Type of neighbourhood in extrema point selection
%  '4m'     two-dimensional four-connected neighborhood, use 'imregionalmax' in Images toolbox
%  '8m'     two-dimensional eight-connected neighborhood, use 'imregionalmax' in Images toolbox
%  '4'     two-dimensional four-connected neighborhood, simple implementation
%  '8'     two-dimensional eight-connected neighborhood, simple implementation
% Note all that these give different results
%
% OUTPUT:
% ix - 3D vector with all IMFs
% resx - 3D vector with all residues
%
%
%
%
%
%
%  
%
% NOTE:
% The interpolation use 'tpaps' in the Matlab 6.5 Splines toolbox. This is
% extremely slow and memory consuming; unless you have a very good computer
% do not use this code on larger images than 128x128. Please make a
% better implementation and tell me about it.
%%


%%

ix=zeros(size(image,1),size(image,2),numberimfs);
resx=zeros(size(image,1),size(image,2),numberimfs);
f=0;

[imf,res,medel]=findimf(image,1,epsilon,conn);

for i=1:numberimfs
   resx(:,:,i)=res;
   ix(:,:,i)=imf;
    
   if i==numberimfs
      f=1;
   end
   [max,min]=maxomin(res,conn);
   [a dummy dummy]=find(max);
   if size(a)<1
      f=1;
   end
  
   [a2 dummy dummy]=find(min);
   if size(a2)<1
      f=1;
   end


   if f==1;
      i  
      break
   end
  
  [imf,res,medel]=findimf(res,1,epsilon,conn);
end
  




%%*******************************************************************
%%

function [imf,res,medel]=findimf(xin,p,epsilon,conn)

sis=size(xin,1);
numpix=sis*sis;
hk=xin-1;
medel=hk;
k=0;
xin1=xin;
for i=1:100
    
  if      sum(sum(abs(medel)))/numpix> epsilon
      
      % k=sum(sum(abs(medel)))/numpix;
       xin=hk;
       [hk, medel,f]=interp(xin,p,conn);
      % b=i;

  end

imf=hk;
res=xin1-hk;

end
%figure;imagesc(imf);colormap(gray);axis off;

%%%%********************************************************************
%%
function [c,medel,f]=interp(x,p,conn)

[y,v]=maxomin(x,conn);
      
[Iy,Jy,Vy] = find(y);
[Iv,Jv,Vv] = find(v);
[Iy,Jy,Vy] = find(y);
  xyy=[Iy';Jy'];
  xyv=[Iv';Jv'];



[m1,n1] = size(y);[m2,n2] = size(v);
[m12,n12,v12] = find(Iy);
[m22,n22,v22] = find(Iv);

  
if size(Iy)<200
  for k=1:2:size(y,1)
    y(1,k)=x(1,k);
    y(size(x,1),k)=x(size(x,1),k);
    y(k,size(x,2))=x(k,size(x,2));
       y(k,1)=x(k,1);
  end
end
[Iv,Jv,Vv] = find(v);
if size(Iv)<200
  for k=1:2:size(v,1)
    v(1,k)=x(1,k);
    v(size(x,1),k)=x(size(x,1),k);
    v(k,size(x,2))=x(k,size(x,2));
       v(k,1)=x(k,1);
  end
end
        
        


if max(size(v12,1),size(v12,2))<3
     medel=zeros(size(x));
     f=1;
c=x-medel;
    return
end
if max(size(v22,1),size(v22,2))<3
   medel=zeros(size(x));
c=x-medel;
f=1;
    return
end



sty = tpaps(xyy,Vy',p); %figure(1);fnplt(sty)
valy=fnval(sty,{1:m1,1:n1});
stv = tpaps(xyv,Vv',p); %figure(2);fnplt(stv)
valv=fnval(stv,{1:m2,1:n2});
f=0;
medel=(valy+valv)/2;
%figure(1);imagesc(valy);colormap(gray);axis off;
%figure(2);imagesc(valv);colormap(gray);axis off;
%figure(3);imagesc(medel);colormap(gray);axis off;
c=x-medel;
%figure(4);imagesc(c);colormap(gray);axis off;





%**************************
%%
%*************************************************************************


function [max,min]=maxomin(x,conn)

y=zeros(size(x));

v=zeros(size(x));



switch conn
    case '4'

  
for j=2:size(x,2)-1
for i=2:size(x,1)-1
  
%     if x(i,j)>x(i+1,j)
%         if x(i,j)>x(i-1,j)
%             if x(i,j)>x(i,j-1)
%                 if x(i,j)>x(i,j+1)
%                    y(i,j)=x(i,j);
%                end
%            end
%        end
%    end

                      if x(i,j)>x(i+1,j) &  x(i,j)>x(i-1,j) & x(i,j)>x(i,j-1)& x(i,j)>x(i,j+1)
               y(i,j)=x(i,j);
                      end

%       if x(i,j)<x(i+1,j)
%         if x(i,j)<x(i-1,j)
%             if x(i,j)<x(i,j-1)
%                if x(i,j)<x(i,j+1)
%                   v(i,j)=x(i,j);
%                end
%             end
%         end
%   end

                    if x(i,j)<x(i+1,j) & x(i,j)<x(i-1,j) & x(i,j)<x(i,j-1) & x(i,j)<x(i,j+1)
                    v(i,j)=x(i,j);
                    end
                    
    end  
end  


case '8'
  
for j=2:size(x,2)-1
for i=2:size(x,1)-1
  
%     if x(i,j)>x(i+1,j)
%       if x(i,j)>x(i+1,j+1)
%         if x(i,j)>x(i+1,j-1)
%           if x(i,j)>x(i,j+1)
%             if x(i,j)>x(i,j-1)
%               if x(i,j)>x(i-1,j)
%                 if x(i,j)>x(i-1,j+1)
%                   if x(i,j)>x(i-1,j-1)
%                    y(i,j)=x(i,j);
%                   end
%                 end
%               end
%             end
%           end
%         end
%       end
%     end
%

             if x(i,j)>x(i+1,j) & x(i,j)>x(i+1,j+1)  & x(i,j)>x(i+1,j-1) & x(i,j)>x(i,j+1) & x(i,j)>x(i,j-1) & x(i,j)>x(i-1,j) &  x(i,j)>x(i-1,j+1) &  x(i,j)>x(i-1,j-1)
             y(i,j)=x(i,j);
        end

                            
%     if x(i,j)<x(i+1,j)
%       if x(i,j)<x(i+1,j+1)
%         if x(i,j)<x(i+1,j-1)
%           if x(i,j)<x(i,j+1)
%             if x(i,j)<x(i,j-1)
%               if x(i,j)<x(i-1,j)
%                 if x(i,j)<x(i-1,j+1)
%                   if x(i,j)<x(i-1,j-1)
%                    v(i,j)=x(i,j);
%                   end
%                 end
%               end
%             end
%           end
%         end
%       end
%     end

                 if x(i,j)<x(i+1,j) & x(i,j)<x(i+1,j+1)& x(i,j)<x(i+1,j-1) & x(i,j)<x(i,j+1) & x(i,j)<x(i,j-1) & x(i,j)<x(i-1,j) & x(i,j)<x(i-1,j+1) & x(i,j)<x(i-1,j-1)
                v(i,j)=x(i,j);
            end
     end    
end  
case '4m'
    
    P = imregionalmax(x,4);
    y = double(P) .* double(x);
    P = imregionalmin(x,4);
    v = double(P) .* double(x);
    
    case '8m'
    
    P = imregionalmax(x,8);
    y = double(P) .* double(x);
    P = imregionalmin(x,8);
    v = double(P) .* double(x);
end

max=y;
min=v;


%**************************************************************************
%*****************************

clc
clear all


disp('PROCESAMIENTO DIGITAL DE  LA SEÑAL PROYECTO BEMD')
disp('IMAGE COMPRESSION, PRESENTED BY LUTIO')


disp('WELCOME TO BEMD FOR IMAGE COMPRESSION')

disp('USE THE FOLLOWING INPUT PARAMETERS:')

X= input(' introduce and specify the name and format of the image you want to compress = ')

figure (1)
imshow(X)
whos
% colormap(map);
% colorbar;
X = double(X(:,:,1));

epsilon = input('epsilon : value of stop criterion, it is more usefull if you use values between 0 and 1 = ')
numberimfs = input('numberimfs :  maximum number of IMF = ')
conn = input('conn : Type of neighbourhood in extrema point selection use any of the following values 4m,8m,4 and 8 = ')

[ix,resx,medel]=IEMD_public(X,epsilon,numberimfs,conn);
figure (2)
imshow(ix)
figre (3)
imshow(resx), title('REMAINDERS')

X=imread('ludacris.gif'); %/aca la imagen en el mismo directorio.
epsilon=0.7;
numberimfs=2;
conn='4'; %aca puede ser '8' o '4m' o '8m'