
# tomografie - Fourier slice theorem

pkg load image

close all


uhlydeg=0:179;
uhly=uhlydeg*pi/180;

dimu=length(uhlydeg);

imagedata=double(imread('mozek_s.png'));

imagedata=imagedata/max(imagedata(:));

size(imagedata)

figure(1)
imshow(imagedata);

[min(imagedata(:)),max(imagedata(:))]


# provedu Radonovu transformaci objektu

sinogram=radon(imagedata,uhlydeg);

sinogram=sinogram/max(sinogram(:));

dimsin=size(sinogram)

figure(2)
imshow(sinogram)

[min(sinogram(:)),max(sinogram(:))]


# inverzni Radonova transformace (zpetna projekce)

figure(3); imshow(iradon(sinogram,uhlydeg),[])


# matice 1D FT projekci

spektr1D=fftshift(fft(ifftshift(sinogram,1)),1);

#spektr1D=spektr1D(39:219,:);  #redukce



dimsp=size(spektr1D);

figure(4)

imshow(abs(spektr1D),[0,max(spektr1D(:))/15])



# vytvareni 2D spektra objektu

dims=dimsp(1);

s=linspace(-1,1,dims);  # souradnice 1D projekce

[xx,yy]=meshgrid(s);  # sit pro rekonstrukci

[uh,ss]=meshgrid(uhly,s); # sit pro sinogram

 

# polarni sit pro sinogram
rr=abs(ss);fi=uh;
idx=find(ss<0);
fi(idx)=fi(idx)+pi;

# konverze na neregularni pravouhlou sit
conver=pol2cart(fi,rr);
conx=conver(:,1);
cony=conver(:,2);
func=spektr1D(:);

size(conx)
size(func)


# interpolace na regularni pravouhlou sit

tic
FTrec=griddata(conx,cony,func,xx,yy);
toc




# NaN ve vysledku (neuspesna interpolace) jsou nahrazeny nulou

idx=isnan(FTrec);

length(idx)

FTrec(idx)=0;


# zpetna 2D FT

rec=abs(fftshift(ifft2(ifftshift(FTrec))));

# otoceni

rec=rec(end:-1:1,:);

figure(5)

imshow(rec,[min(rec(:)),max(rec(:))/1])






