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 次查看(过去 30 天)
显示 更早的评论
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 个评论
Geoff
2012-5-23
Could you please edit this question and format the code using the 'code' button in the toolbar? It's unreadable right now.
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!