I have problem with the following code for fan beam image reconstruction matlab gives the following error ??? Error using ==> interp1 at 165 The interpolation points XI should be real . Please help of you have any idea why its not working.thanks
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
clear all; close all;
%figure(1);
%x=linspace(-1,1);
%y=x;
%[xx,yy]=meshgrid(x,y);
%rr=xx.^2+yy.^2;
%B=(rr<=9/16)-0.75*(rr<=1/4)+0.25*(rr<=1/16);
%imshow(B,[0 1]);
%pause;
%a is half the image diagonal dimension
%a=sqrt(size(B,1)^2 + size(B,2)^2)/2;
%a;
%a=71.4178;
%D is the distance from the fan-beam vertex to the center of rotation
%D should be a few pixels larger than a
D=80;
alpha=1;
gamam=20;
t=-1:0.05:1;
gama=-gamam:alpha:gamam;
rb1=radon_bull(gama);
rb= radon_bull(t);
figure(3); hold on;
plot(t,rb,'k-','Linewidth',3);
%plot(x,rb,'r.','Markersize',10);
title('Radon transform of Bulls Eye');
pause;
r2_vals=r2(t);
figure(4);
plot(t,r2_vals,'k-','Linewidth',3);
pause;
r3_vals=r2(gama);
n=-gamam/alpha:gamam/alpha;
g2=g(n,alpha);
g2 = [g2(gamam+1:2*gamam+1) g2(1:gamam)];
h=convp(r3_vals,g2);%this is a convolution for one beta angle
figure(5)
plot(h,'k-','Linewidth',3);
pause;
s = -gamam:1/50:gamam;
Ih = interp1(n,h,s,'nearest');
figure(6); hold on;
plot(s,Ih,'k-','Linewidth',3);
title('Interpolation of h');
pause;
dx=1/100;
x=0:dx:2;
y=0:dx:2;
[xx yy] = meshgrid(x,y);
l=-gamam:gamam;
B=back3(l,h,xx,yy,D);
figure(7); hold on;
imshow(B);
axis image;
axis xy;
title('weighted back projection')
%and radon_bull.m
function rb = radon_bull(t)
%
t = abs(t);
rb = 0.5*sqrt(9 - 16*t.^2).*(t <= 3/4) - 0.75*sqrt(1 - 4*t.^2).*(t <= 0.5) + 1/8*sqrt(1 - 16*t.^2).*(t <= 0.25);
%r2.m is
function [R]=r(x)
R=0.5*sqrt(9-16*abs(x).^2).*(abs(x)<=3/4)-0.75*sqrt(1-4*abs(x).^2).*(abs(x)<=0.5)+1/8*sqrt(1-16*abs(x).^2).*(abs(x)<=0.25);
end
%g.m
function [G]=g(n,alpha)
%G=(1/(4*pi))*(n/sin(n))^2*rlf(n);
%
%g(gama)=1/(4*pi)*((gama)/sin(gama))^2*rlf(gama)
%rlf is the ram-lak filter for L=10;
for i=1:length(n)
G(i)=1/(4*pi)*(i*alpha/sin(i*alpha))^2*rlf(i);
end
end
%rlf.m
function [RLF]=rlf(n)
for i=1:length(n)
if n(i)==0;
RLF(i)=100./(2*pi);
elseif mod(n(i),2)==0;
RLF(i)=0;
else
RLF(i)=-200./(pi^3*(n(i).^2));
end
end
%l.m
function [L] = l(r,fi,beta,D)
L=sqrt((D+r*sin(beta-fi))^2+(r*cos(beta-fi))^2);
end
back3.m
function B = back3(t,h,x,y,D)
%
%r=sqrt(x^2+y^2);
%fi=acos(x/r);
%B=zeros(M,M);
%=-1:1/100:1;
%y=-1:1/100:1;
%[fi,r] = cart2pol(x,y);
%gamam=20;
%[M M] = size(x);
%B = zeros(M,M);
[M M] = size(x);
B = zeros(M,M);
r=sqrt(x^2+y^2);
fi=acos(x/(sqrt(x^2+y^2)));
%for beta=0:36/sqrt(2):2*180
for beta=0:10:360
v = atan(r*cos(beta-fi)/(D+r*sin(beta-fi)));
B(:) = B(:) + interp1(t,h,v(:),'nearest')*1/(l(r,fi,beta,D)^2);
% B(:) = B(:) + interp1(t,h,v(:),'linear');
%x=r*cos(fi);
%y=r*sin(fi);
end;
B=B*10;
convp.m
function h = convp(f,g)
%CONVP Convolution of discrete periodic functions.
% h = CONVP(f,g) convolves the discrete N-periodic functions f and g.
% The resulting function is a discrete N-periodic function.
N =length(f);
for m=1:N
h(m) = f(1:m)*g(m:-1:1)' + f(m+1:N)*g(N:-1:m+1)';
end;
2 comentarios
Geoff
el 23 de Mayo de 2012
Could you please edit this question and format the code using the 'code' button in the toolbar? It's unreadable right now.
Respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!