Matlab程序 关于FBP

%%
clc;
clear;
close all;
tic; %计时开始
imageN=128;       %图像大小imageN*imageN，即64*64
I= phantom（imageN）;%读取Shepp-Logan头像
[imageNimageN]=size（I）;
subplot（221）%最后显示的图有2行2列，现在显示的是第一个
imshow（I）;
title（‘64*64原始图像‘）;%显示原始图像
%%
%4096*4096的矩阵
nx=imageN;       %图像的宽64
ny=imageN;       %图像的高64

na=128;          %投影角度64
nb=128;          %每个角度下的射线数64

mask=true（nxny）;  %预先生成64*64的单位矩阵，加快速度
ray_pix = 1;       % 像素大小

x = [0:nx-1] - （nx-1）/2;
y = [0:ny-1] - （ny-1）/2;

[xy] = ndgrid（x y）; %相对于中心的坐标转换
x = x（mask（:））;
y = y（mask（:））;
np = length（x）;		 

angle = [0:na-1]‘/na * pi;
tau = cos（angle） * x‘ + sin（angle） * y‘;	
tau = tau / ray_pix;		
tau = tau + （nb+1）/2;		% 射线编号从1开始

ibl = floor（tau）;		    % 左半部

% floor 为向后取整1.6-1;  1.2-1;

val = 1 - （tau-ibl）;		% 左半部的加权