
clc
clear
imsize=160;
n=imsize;
nhlf=n/2;


%bispect=zeros(imsize,imsize,imsize,imsize);
bispect=zeros(1,15e7);
bispect2=zeros(1,15e7);
dhaka=zeros(imsize,imsize);
dingdong=zeros(imsize,imsize);
dingdong2=zeros(imsize,imsize);
mark=zeros(imsize,imsize);
oops=0
para=1;

   disp( 'Calculating Bispectrum')

frames=24;

for iter=1:frames
    
  disp(strcat(num2str(iter), ' frame'))
    
    
   jacky=num2str(iter);
   image=im2double(imread(strcat('/home/arun/binary/image',jacky,'.png')));
   fim = fft2(image);
   sfim=fftshift(fim);
   dhaka=dhaka+abs(sfim).*abs(sfim);
   sfim=exp(i*angle(sfim));
   sfim2=fliplr(sfim);
   
   mark=zeros(imsize,imsize);
   mark(nhlf,nhlf)=1;
   mark(nhlf+1,nhlf)=1;
   mark(nhlf,nhlf+1)=1;

   
   
   
   
   
   tic 
    
flag=1;

for u=-(n-2)/2:n/2
   for v=-(n-2)/2:n/2
      for u1=-(n-2)/2:n/2
          for v1=-(n-2)/2:n/2
         
            if (u+u1) <= n/2
            if (v+v1) <= n/2
            if (u+u1) >= -(n-2)/2
            if (v+v1) >= -(n-2)/2
                
                 
            if(mark(nhlf+u,nhlf+v)~=0)
            if  mark(nhlf+u1,nhlf+v1)~=0       
            
              
             
             bispect(flag)=bispect(flag)+  sfim(nhlf+u,nhlf+v)*sfim(nhlf+u1,nhlf+v1)*conj(sfim(nhlf+u+u1,nhlf+v+v1));
             mark(nhlf+u+u1,nhlf+v+v1)=mark(nhlf+u+u1,nhlf+v+v1)+1;
             if iter==frames
                  %bispect(flag)= bispect(flag)/abs(bispect(flag));
                  bispect(flag)= bispect(flag)/frames;
             end
             

             
             
             flag=flag+1;
             
             end
             end
            end
            end
            end
            end
              
end            
end
        
end
end

toc

mark=zeros(imsize,imsize);
mark(nhlf,nhlf)=1;
 mark(nhlf+1,nhlf)=1;
 mark(nhlf,nhlf+1)=1;

   
   
   
   
   
   tic 
    
flag=1;

for u=-(n-2)/2:n/2
   for v=-(n-2)/2:n/2
      for u1=-(n-2)/2:n/2
          for v1=-(n-2)/2:n/2
         
            if (u+u1) <= n/2
            if (v+v1) <= n/2
            if (u+u1) >= -(n-2)/2
            if (v+v1) >= -(n-2)/2
                
                 
            if(mark(nhlf+u,nhlf+v)~=0)
            if  mark(nhlf+u1,nhlf+v1)~=0       
            
              
             
             bispect2(flag)=bispect2(flag)+  sfim2(nhlf+u,nhlf+v)*sfim2(nhlf+u1,nhlf+v1)*conj(sfim2(nhlf+u+u1,nhlf+v+v1));
             mark(nhlf+u+u1,nhlf+v+v1)=mark(nhlf+u+u1,nhlf+v+v1)+1;
             if iter==frames
                 % bispect2(flag)= bispect2(flag)/abs(bispect2(flag));
                  bispect2(flag)= bispect2(flag)/frames;
             end
             

             
             
             flag=flag+1;
             
             end
             end
            end
            end
            end
            end
              
end            
end
        
end
end
toc




end



              







dhaka=dhaka/frames;
dhaka=sqrt(dhaka);



dingdong(nhlf,nhlf)=sfim(nhlf,nhlf);
dingdong(nhlf+1,nhlf)=sfim(nhlf+1,nhlf);
dingdong(nhlf,nhlf+1)=sfim(nhlf,nhlf+1);

%dingdong(1,1)=-1;
%dingdong(2,1)=-1;
%dingdong(1,2)=-1;
 mark=zeros(imsize,imsize);
mark(nhlf,nhlf)=1;
mark(nhlf+1,nhlf)=1;
mark(nhlf,nhlf+1)=1;









disp('Recovering')
pause(1);
flag=0;
tic

for iter=1:1
flag=1;
for u=-(n-2)/2:n/2
   for v=-(n-2)/2:n/2
      for u1=-(n-2)/2:n/2
          for v1=-(n-2)/2:n/2
         
            if (u+u1) <= n/2
            if (v+v1) <= n/2
            if (u+u1) >= -(n-2)/2
            if (v+v1) >= -(n-2)/2
                
       
                
          
                
                
             if(mark(nhlf+u,nhlf+v)~=0)
             if  mark(nhlf+u1,nhlf+v1)~=0
               
               
                      
             
            dingdong(nhlf+u+u1,nhlf+v+v1)= dingdong(nhlf+u+u1,nhlf+v+v1)+dingdong(nhlf+u,nhlf+v)*(dingdong(nhlf+u1,nhlf+v1))*conj(bispect(flag));            
            mark(nhlf+u+u1,nhlf+v+v1)=mark(nhlf+u+u1,nhlf+v+v1)+1;  
             
            if dingdong(nhlf+u+u1,nhlf+v+v1) ~= 0

                 
           
                 if mark(nhlf+u+u1,nhlf+v+v1)>1
                %dingdong(nhlf+u+u1,nhlf+v+v1)=dingdong(nhlf+u+u1,nhlf+v+v1)/2;
                dingdong(nhlf+u+u1,nhlf+v+v1)=dingdong(nhlf+u+u1,nhlf+v+v1)/abs(dingdong(nhlf+u+u1,nhlf+v+v1));
                 end  
                            
            end
           
            
                
           
          flag=flag+1;
      
            end
            end
            end
            end
            end
            end
              
end            
end
        
end
end
toc
end

dingdong2(nhlf,nhlf)=sfim2(nhlf,nhlf);
dingdong2(nhlf+1,nhlf)=sfim2(nhlf+1,nhlf);
dingdong2(nhlf,nhlf+1)=sfim2(nhlf,nhlf+1);

%dingdong(1,1)=-1;
%dingdong(2,1)=-1;
%dingdong(1,2)=-1;
 mark=zeros(imsize,imsize);
mark(nhlf,nhlf)=1;
mark(nhlf+1,nhlf)=1;
mark(nhlf,nhlf+1)=1;
tic

flag=1;
for u=-(n-2)/2:n/2
   for v=-(n-2)/2:n/2
      for u1=-(n-2)/2:n/2
          for v1=-(n-2)/2:n/2
         
            if (u+u1) <= n/2
            if (v+v1) <= n/2
            if (u+u1) >= -(n-2)/2
            if (v+v1) >= -(n-2)/2
                
       
                
          
                
                
             if(mark(nhlf+u,nhlf+v)~=0)
             if  mark(nhlf+u1,nhlf+v1)~=0
               
               
                      
             
            dingdong2(nhlf+u+u1,nhlf+v+v1)= dingdong2(nhlf+u+u1,nhlf+v+v1)+dingdong2(nhlf+u,nhlf+v)*(dingdong2(nhlf+u1,nhlf+v1))*conj(bispect2(flag));            
            mark(nhlf+u+u1,nhlf+v+v1)=mark(nhlf+u+u1,nhlf+v+v1)+1;  
             
            if dingdong2(nhlf+u+u1,nhlf+v+v1) ~= 0

                 
           
                 if mark(nhlf+u+u1,nhlf+v+v1)>1
               % dingdong2(nhlf+u+u1,nhlf+v+v1)=dingdong2(nhlf+u+u1,nhlf+v+v1)/2;
                dingdong2(nhlf+u+u1,nhlf+v+v1)=dingdong2(nhlf+u+u1,nhlf+v+v1)/abs(dingdong2(nhlf+u+u1,nhlf+v+v1));
                 end  
                            
            end
           
            
                
           
          flag=flag+1;
      
            end
            end
            end
            end
            end
            end
              
end            
end
        
end
end
toc

dingdong=dingdong+fliplr(dingdong2);
dingdong=dingdong+(flipud(fliplr(conj(dingdong))));














subplot(2,2,1)
imshow(angle(sfim));
subplot(2,2,2)
imshow(angle(dingdong));
subplot(2,2,3)
imshow(image);
subplot(2,2,4)
imshow(abs(ifft2((dhaka.*exp(i*angle(dingdong))))))
figure(2)
imshow(mark)