Hi All,
I am new to Julia, coming from Matlab, recently started working on OPT (Optical Projection Tomography) project where one of requirements is open source.
I tried to implement something similar to Matlab radon/iradon pair, using as a base the C++ code developed previously by T. Knopp (https://github.com/tknopp/numcpp).
For debugging, i mirrored Julia implementation of fbp in Matlab.
Quality-wise, sinograms and reconstructions of phantom are in good agreement with each other and with original.
(there is one minor scale issue on iradon, as I didn’t find suitable “imresize”, but that can be fixed later)
The main drawback is that the “naïve” Julia functions are much slower than optimized Matlab’s.
For 1000x1000 phantom image, the ratios of processing times were
radon - 16.7
iradon - 20.5
iradon translated to Matlab from Julia - 23.5
That is not practical for the project.
The main issue looks like simply slow brackets.
I paste the code and test program below, hope to hear any suggestion on speed-up.
(It’s amateurish but sort of does the job).
I think sooner or later the question of radon/iradon will appear in Julia anyway.
# tomography.jl
########################################################
function radon(A,angles)
angles = angles + pi/2;
w,h = size(A);
Nr_div2 = floor(sqrt(w*w+h*h)/2);
Nr = Int64(Nr_div2*2+1);
L = length(angles);
sinogram = zeros(Nr,L);
SIN = sin(angles);
COS = cos(angles);
R = linspace(-0.5,0.5,Nr)*Nr;
RSIN = zeros(Nr,L);
RCOS = zeros(Nr,L);
w2 = w/2;
h2 = h/2;
for a=1:L
for k=1:Nr
RSIN[k,a]=R[k]*SIN[a];
RCOS[k,a]=R[k]*COS[a];
end
end
for a=1:L
for k=1:Nr
for l=1:Nr
x = RCOS[k,a] - RSIN[l,a] + w2;
y = RSIN[k,a] + RCOS[l,a] + h2;
x1 = floor(x);
x2 = ceil(x);
y1 = floor(y);
y2 = ceil(y);
if !(x1 <= 0 || x2 > w || y1 <= 0 || y2 > h)
# BILINEAR INTERPOLATION - STARTS
# x1 y1 - btm left
# x1 y2 - top left
# x2 y1 - btm rght
# x2 y2 - top rght
f11 = A[Int64(x1),Int64(y1)];
f21 = A[Int64(x1),Int64(y2)]; # that looks like, correct way
f12 = A[Int64(x2),Int64(y1)];
f22 = A[Int64(x2),Int64(y2)];
value = 0;
if x2==x1 && y2==y1
value = f11;
elseif x2==x1
value = 1/(y2-y1)*( f11*(y2-y) + f22*(y-y1) );
elseif y2==y1
value = 1/(x2-x1)*( f11*(x2-x) + f22*(x-x1) );
else
value = 1/(x2-x1)/(y2-y1)*( f11*(x2-x)*(y2-y) +
f21*(x-x1)*(y2-y) +
f12*(x2-x)*(y-y1) +
f22*(x-x1)*(y-y1) );
end
# BILINEAR INTERPOLATION - ENDS
sinogram[k,a] += value;
end
end
end
end
return sinogram;
end
#########################################
function iradon(A,angles) # input angles [rad]
angles = angles + pi/2;
N, nAngles = size(A);
I = zeros(N,N); # reconstruction
x = linspace(-0.5,0.5,N);
filter = abs(linspace(-1, 1, N));
# FT domain filtering
for t=1:length(angles)
fhat = fftshift(fft(slice(A,:,t)));
A[:,t] = real(ifft(ifftshift(fhat.*filter)));
end
XCOS = zeros(N,length(angles));
XSIN = zeros(N,length(angles));
for k=1:N
for a=1:length(angles)
XCOS[k,a]=x[k]*cos(angles[a]);
XSIN[k,a]=x[k]*sin(angles[a]);
end
end
for m=1:N
for n=1:N
for t=1:nAngles
r = XCOS[m,t] + XSIN[n,t];
index = Int64(round((r/sqrt(2.0) + 0.5)*N)); # sqrt(2) magnified
# index = Int64(round((r + 0.5)*N)); # to keep pixel resolution
if index>0 && index<=N
I[m,n] += A[index,t];
end
end
end
end
return I;
end
test program:
include("tomography.jl")
using Images, TestImages, ImageView
angles = (0:1:359)/360*2*pi; # [rad]
z = shepp_logan(1000; highContrast=true);
ImageView.imshow(z)
z = float(z);
tic();
sinogram = radon(z,angles);
display(toc());
v1 = minimum(sinogram);
v2 = maximum(sinogram);
todisplay = (sinogram-v1)/(v2-v1);
ImageView.imshow(todisplay);
tic();
reconstruction = iradon(sinogram,angles);
display(toc());
v1 = minimum(reconstruction);
v2 = maximum(reconstruction);
todisplay = (reconstruction-v1)/(v2-v1);
ImageView.imshow(todisplay);
it would be really good if this were published in a package!
