逆拉东变换（inverse Radontransform）【MATLAB函数】

iradon函数

I=iradon(R,theta,interp,filter,frequency_scaling,output_size)

iradon函数是基于R-L滤波器的滤波反投影法。实现重建图像的过程如下：

1. 把投影矩阵R转换到频域，生成fft（R）；
2. fft（R）和滤波函数H相乘，得到滤波后的频域投影矩阵fft（R）*H；
3. 把fft（R）*H 转换到空域，得到空域中的滤波后的投影矩阵R’=ifft（fft（R）*H ）;
4. R’插值后，得到处理好的投影矩阵R’’;
5. 直接反投影得到重建图像矩阵I。

参数介绍

• R是投影矩阵。

• theta可以是一个包含所有扫描角度的向量，这时每两个相邻角度间隔相等；也可以是一个标量值，等于相邻两个扫描角度的间隔。

• interp是插值函数，有以下几种差值方式可以选择：
  nearest：最邻近插值方法（nearest neighbor interpolation）。这种插值方法在已知数据的最邻近点设置插值点，对插值点的数值进行四舍五入，对超出范围的数据点返回NaN。
  linear：线性插值（linear interpolation），这是interp中的默认数值。该方法采用直线将相邻的数据点相连，对超出数据范围的数据点返回NaN。（执行速度较快，有足够的精度，最为常用。）
  spline：三次样条插值（cubic spline interpolation），该方法采用三次样条函数获取插值数据点，在已知点为端点的情况下，插值函数至少具有相同的一阶和二阶导数。（执行速度最慢，精度高，最平滑。）
  pchip：分段三次厄米多项式差值（piecewise cubic Hermite interpolation）。
  cubic：三次多项式插值，与分段三次厄米多项式插值方法相同。
  v5cubic：MATLAB5中使用的三次多项式插值。
  nearest 最快的插值方法，但是数据平滑方面最差，数据是不连续的。
  cubic 较慢，精度高，平滑度好，当希望得到平滑的曲线时可以使用该选项。

• filter是滤波函数，有以下几种滤波器可以选择：

滤波函数H	含义
none	没有滤波
ram-lak	R-L滤波函数的傅里叶函数和频域中每个角度的投影相乘，实现滤波
Shepp-Logan	sinc函数*R-L函数——
Cosine	cosine函数*R-L函数——
Hamming	Hamming函数*R-L函数——
Hann	hann函数*R-L函数——
其他	设定的函数*R-L函数——

• frequency_scaling是一个标量值，取值范围[0,1]，通过缩放滤波函数的频率修改滤波函数。
  默认值为1。如果取值小于1，则滤波函数的频率在归一化后被压缩到适合[0，frequency_scaling]的范围。

比如默认值为1时，滤波函数的窗口为[0,10Hz],frequency_scaling=0.5时，滤波函数的窗口为[0,20Hz]。（待确定）在频域中，频率若大于frequency_scaling对应的频率，则该频率处的函数值为0。

• output_size是一个标量，用来规定重建图像的行数和列数。
  默认等于2 * floor(size(R,1)/(2*sqrt(2)))。改变output_size会改变重建图像的大小，但是不会改变像素点的个数。

radon函数

[R,Xp]=radon（I，theta，N）

参数介绍

如果theta是标量，返回R是列向量 ，表示theta角度下图像I的radon变换。
如果theta是向量，返回R是矩阵，每一列表示某一theta角度下图像I的radon变换。
忽略theta，则默认为是0:179的向量。
返回行向量Xp，投影数据R的每一行的射线坐标。

已设N，则用N个点计算投影， 且R有N行。
未设N，则用 2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3 个点计算投影，即使射线通过图像对角线，这个值也足够。

另： 调用c程序radonc( )实现radon变换
mfilename函数中，I是第一个变量，theta是第二个变量，N是第三个变量。如果生成投影的行维度r不等N，则进行N等分，线性插值，变成行维度为N的投影。

附上用c++实现的代码：

#include <math.h>   
#include "mex.h"   

static void radon(double *pPtr, double *iPtr, double *thetaPtr, int M<