逆拉东变换（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, int N,
	int xOrigin, int yOrigin, int numAngles, int rFirst,
	int rSize);


/**
**pr：进行radon变换后输出矩阵的对于一个特定theta角的列的首地址
**pixel：要进行radon变换的像素值乘以0.25以后的值（由于每一个像素点取了相邻四个点提高精度，故在计算时pixel也要相应乘以0.25，类似于一个点占0.25的比例，然后四个点刚好凑足1的份额）
**r：进行radon变换的该点与初始的r值——rFirst之间的差
**/
void incrementRadon(double *pr, double pixel, double r)
{
	int r1;
	double delta;

	r1 = (int)r;   //对于每一个点，r值不同，所以，通过这种方式，可以把这一列中相应行的元素的值给赋上
	delta = r - r1;
	pr[r1] += pixel * (1.0 - delta); //radon变换本来就是通过记录目标平面上某一点的被映射后点的积累厚度来反推原平面的直线的存在性的，故为+=  
	pr[r1 + 1] += pixel * delta;  //两个点互相配合，提高精度 
}
/***
**参数解释：
**pPtr：经过radon变换后输出的一维数组，该一维数组是其实要还原成一个rSize*numAngles的矩阵
**iPtr：需要进行radon变换的矩阵的一维数组存储形势
**thetaPtr：指定进行radon变换的弧度的数组，该角度就是极坐标中偏离正方向的角度
**M：要进行radon变换的矩阵的行数
**N：要进行radon变换的矩阵的列数
**xOrigin：要进行radon变换的矩阵的的中心的横坐标
**yOrigin：要进行radon变换的矩阵的中心的纵坐标
**numAngles：thetaPtr数组中元素的个数
**rFist：极坐标中初始点与变换原点的距离
**rSize：整个radon变换中极坐标的点之间的最远距离
***/
static void radon(double *pPtr, double *iPtr, double *thetaPtr, int M, int N,
	int xOrigin, int yOrigin, int numAngles, int rFirst, int rSize)
{
	int k, m, n;              /* loop counters */
	double angle;             /* radian angle value */
	double cosine, sine;      /* cosine and sine of current angle */
	double *pr;               /* points inside output array */
	double *pixelPtr;         /* points inside input array */
	double pixel;             /* current pixel value */
	double *ySinTable, *xCosTable;
	/* tables for x*cos(angle) and y*sin(angle) */
	double x, y;
	double r, delta;
	int r1;

	/* Allocate space for the lookup tables */
	xCosTable = (double *)mxCalloc(2 * N, sizeof(double));  //MATLAB的内存申请函数，对应C语言可以换成calloc函数 
	ySinTable = (double *)mxCalloc(2 * M, sizeof(double));

	for (k = 0; k < numAngles; k++)
	{
		//每一个theta角，经过radon变化后，就会产生一列数据，这一列数据中，共有rSize个数据
		angle = thetaPtr[k];
		pr = pPtr + k*rSize;  /* pointer to the top of the output column */
		cosine = cos(angle);
		sine = sin(angle);

		/*
		**radon 变换中，极坐标下，沿r轴的theta角和每一个像素点的分布都是非线性的，而此处采用的是线性radon变换，
		**为了提高精度，把每一个像素点分成其四周四个相邻的像素点来进行计算！x、y坐标的误差是正负0.25
		*/
		for (n = 0; n < N; n++)
		{
			x = n - xOrigin;
			xCosTable[2 * n] = (x - 0.25)*cosine;   //由极坐标的知识知道，相对于变换的原点，这个就是得到了该点的横坐标
			xCosTable[2 * n + 1] = (x + 0.25)*cosine;
		}
		for (m = 0; m < M; m++)
		{
			y = yOrigin - m;
			ySinTable[2 * m] = (y - 0.25)*sine;   //同理，相对于变换的原点，得到了纵坐标
			ySinTable[2 * m + 1] = (y + 0.25)*sine;
		}

		pixelPtr = iPtr;
		for (n = 0; n < N; n++)
		{
			for (m = 0; m < M; m++)   //便利原矩阵中的每一个像素点
			{
				pixel = *pixelPtr++;
				if (pixel != 0.0)   //如果该点像素值不为0，也即图像不连续
				{
					pixel *= 0.25;

					//一个像素点分解成四个临近的像素点进行计算，提高精确度
					r = xCosTable[2 * n] + ySinTable[2 * m] - rFirst;
					incrementRadon(pr, pixel, r);

					r = xCosTable[2 * n + 1] + ySinTable[2 * m] - rFirst;
					incrementRadon(pr, pixel, r);

					r = xCosTable[2 * n] + ySinTable[2 * m + 1] - rFirst;
					incrementRadon(pr, pixel, r);

					r = xCosTable[2 * n + 1] + ySinTable[2 * m + 1] - rFirst;
					incrementRadon(pr, pixel, r);
				}
			}
		}
	}

	mxFree((void *)xCosTable);   //MATLAB的内存释放函数，对应C语言可以换成free函数
	mxFree((void *)ySinTable);
}