Separation of Image Mixture using FastICA

Image separation of mixed and overlapped images is a frequent problem in computer vision (image processing). The following Matlab source code is a demonstration of image separation using FastICA algorithm based on kurtosis.

clc;clear all;close all; 

original_image1=imread ('input1.jpg');  
original_image2=imread ('input2.jpg'); 

figure()
subplot(3,3,1),imshow(original_image1),title('Original Source 1'),
subplot(3,3,2),imshow(original_image2),title('Original Source 2'), 

[x1,y1,z1]=size(original_image1); 
[x2,y2,z2]=size(original_image2); 
 
s1=[reshape(original_image1,[1,x1*y1*z1])]; 
s2=[reshape(original_image2,[1,x2*y2*z2])]; 

             
mixtures=((0.99-0.4).*rand(2,2) + 0.4 )*double([s1;s2]);       

mixture_1=uint8(round(reshape(mixtures(1,:),[x1,y1,z1]))); 
mixture_2=uint8(round(reshape(mixtures(2,:),[x2,y2,z2]))); 

subplot(3,3,4),imshow(mixture_1),title('Mixture 1'), 
subplot(3,3,5),imshow(mixture_2),title('Mixture 2'), 

mixtures_origin=mixtures;                         
mixtures_mean=zeros(2,1); 

for i=1:2 
    mixtures_mean(i)=mean(mixtures(i,:)); 
end                                        
 
for i=1:2 
    for j=1:size(mixtures,2) 
        mixtures(i,j)=mixtures(i,j)-mixtures_mean(i); 
    end 
end

mixtures_cov=cov(mixtures');                 

[E,D]=eig(mixtures_cov);

Q=inv(sqrt(D))*(E)';                      
mixtures_white=Q*mixtures;                    

IsI=cov(mixtures_white');                        
                         
[VariableNum,SampleNum]=size(mixtures_white); 
numofIC=VariableNum;                       
B=zeros(numofIC,VariableNum);

for r=1:numofIC                           
    i=1;maxIterationsNum=250;             
    b=2*(rand(numofIC,1)-.5);                   
    b=b/norm(b);                         
    while i < maxIterationsNum
       if i == maxIterationsNum           
            fprintf('No convergence - %d - %d', r,maxIterationsNum); 
            break; 
       end   
       b_prev=b;                           
       t=mixtures_white'*b; 
       G = 4 * t.^3;
       Gp = 12 * t.^2;
       b=(t'*G*b+mixtures_white*G)/SampleNum-mean(Gp)*b;                                          
       b=b-B*B'*b;                       
       b=b/norm(b);       
       if abs(abs(b'*b_prev)-1) < 1e-10      
           B(:,r)=b;                  
           break; 
       end
       fprintf(1,'%d  |  %d  -  %d\n',r,i,abs(abs(b'*b_prev)-1));      
       i=i+1;
    end
end 

ICAedS=abs(55*(B'*Q*mixtures_origin)); 

original_image1_icaed =uint8 (round(reshape(ICAedS(1,:),[x1,y1,z1]))); 
original_image2_icaed =uint8 (round(reshape(ICAedS(2,:),[x2,y2,z2]))); 

subplot(3,3,7),imshow(original_image1_icaed),title('Restored Image 1'); 
subplot(3,3,8),imshow(original_image2_icaed),title('Restored Image 2');