%Bispectrum Radon Transform Reconstruction for 2D Image
%copyright 2010 Arun Surya
%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 1 of the License. 
%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.
%You should have received a copy of the GNU General Public License along with %this program.  If not, see <http://www.gnu.org/licenses/>.
close all
clear all
x=imread('C:\camera.png');
x=im2double(x);
rad=radon(x,1:10:180);
rad=rad(1:366,:);
for p=1:18
x=im2double(x);
n=366;
nhlf=n/2;
nhlfp=nhlf+1;
nhlfm=nhlf-1;
nbi=n;
kill=1:256;
nbihlfp=nbi/2 +1;
im=rad(1:366,p)';
%im(5)=25;
%im(45)=5;
%im(60)=2;
%im(80)=10;
%im=rectpuls(kill-60,20) + 3*rectpuls(kill-100,10)  ;
%im=rectpuls(kill-60,95)
%im=rectpuls(kill-5,2)
%im=rectpuls(kill-60,10) + 3*rectpuls(kill-100,10);
figure(1);
subplot(2,1,1)
plot(im)
title(strcat('Projection  ',num2str(p)));
pause(3);
fim = fft(im);
sfim=fftshift(fim);

bispect=zeros(nbi,nbi);
count=1;
for v=1:nhlf
  for u=1:nhlf
       if (nhlf+u+v) <= n 
         if (nhlf+u+v) >= 1
           bispect(nhlf+v, nhlf+u)= sfim(nhlf+u).*sfim(nhlf+v) .*conj(sfim(nhlf+u+v));
           
end            
end
end
end
dingdong=ones(1,n);
sfimangle=angle(sfim);
abssfim=exp(i*sfimangle);
%dingdong(1)=sfim(1)
dingdong(1)=abssfim(1);
%fultose=angle(bispect);
%bispect=exp(i.*fultose);
count=0
for v=1:nhlf
  for u=1:nhlf
      
       if (nhlf+u+v) <= n
           if (nhlf+u+v) >= 1
               
               
           dingdong(nhlf+u+v)=dingdong(nhlf+u+v)+dingdong(nhlf+u)*dingdong(nhlf+v)*conj(bispect(nhlf+v, nhlf+u));
           dingdong(nhlf+u+v)=dingdong(nhlf+u+v)/abs(dingdong(nhlf+u+v));
           count=count+1;
           
end            
end
end
end
for j=2:nhlf
    dingdong(j)=conj(dingdong(n-j+2));
end
pukasa=angle(dingdong);
%pukasa=angle(dingdong);
koo=(ifft(abs(sfim).*exp(i.*angle(dingdong))));
%koo=fftshift(fft(abs(sfim).*exp(i.*angle(sfim))));
pukasi=angle(sfim);
%koo=koo/max(max(koo));
jacky(1:366,p)=abs(koo);

subplot(2,1,2)
plot(im);
plot(abs(koo));
title(strcat('Reconstructed  ',num2str(p)));
pause(3);

end

I2 = iradon(jacky,1:10:180);
figure(2);
imshow(I2);



