radon变换和逆变换的实现用C＋＋？

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, p /* 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 rIdx /* r value offset from initial array element */

int rLow/* (int) rIdx */

double pixelLow /* amount of pixel's mass to be assigned to */

/* the bin below Idx */

double *yTable, *xTable /* x- and y-coordinate tables */

double *ySinTable, *xCosTable

/* tables for x*cos(angle) and y*sin(angle) */

/* Allocate space for the lookup tables */

yTable = (double *)calloc(2*M, sizeof(double))

xTable = (double *)calloc(2*N, sizeof(double))

xCosTable = (double *)calloc(2*N, sizeof(double))

ySinTable = (double *) calloc(2*M, sizeof(double))

/* x- and y-coordinates are offset from pixel locations by 0.25 */

/* spaced by intervals of 0.5. */

/* We want bottom-to-top to be the positive y direction */

yTable[2*M-1] = -yOrigin - 0.25

for (k = 2*M-2k >=0k--)

yTable[k] = yTable[k+1] + 0.5

xTable[0] = -xOrigin - 0.25

for (k = 1k <2*Nk++)

xTable[k] = xTable[k-1] + 0.5

for (k = 0k <numAnglesk++) {

angle = thetaPtr[k]

pr = pPtr + k*rSize /* pointer to the top of the output column */

cosine = cos(angle)

sine = sin(angle)

/* Radon impulse response locus: R = X*cos(angle) + Y*sin(angle) */

/* Fill the X*cos table and the Y*sin table. Incorporate the */

/* origin offset into the X*cos table to save some adds later. */

for (p = 0p <2*Np++)

xCosTable[p] = xTable[p] * cosine - rFirst

for (p = 0p <2*Mp++)

ySinTable[p] = yTable[p] * sine

/* Remember that n and m will each change twice as fast as the */

/* pixel pointer should change. */

for (n = 0n <2*Nn++) {

pixelPtr = iPtr + (n/2)*M

for (m = 0m <2*Mm++) {

pixel = *pixelPtr

if (pixel) {

pixel *= 0.25

rIdx = (xCosTable[n] + ySinTable[m])

rLow = (int) rIdx

pixelLow = pixel*(1 - rIdx + rLow)

pr[rLow++] += pixelLow

pr[rLow] += pixel - pixelLow

}

if (m%2)

pixelPtr++

}

}

}

free((void *) yTable)

free((void *) xTable)

free((void *) xCosTable)

free((void *) ySinTable)

}

theta=0:180theta是radon变换的角度，可以自由设置，这里是从0°到180°一共181个角度，如果只求等分的60个角度的radon变换，可以改成theta=0:3:177；如果只求一个特殊角度的radon变换可以写成theta=x这里x就是希望的角度。

R = radon(I, theta)

R是存储radon变换的值，它是一个矩阵，列数是theta的个数，表示每一个角度生成一列，行数是被处理的矩阵（I）对角线的长度。

R一般用于inverse radon transform，程序如下：

I1=iradon(R,theta)

imshow(I1)

iradon命令有很多参数，具体使用方法请参考http://www.mathworks.com/access/helpdesk/help/toolbox/images/index.html?/access/helpdesk/help/toolbox/images/iradon.html&http://www.mathworks.com/cgi-bin/texis/webinator/search?pr=Whole_site&db=MSS&prox=page&rorder=750&rprox=750&rdfreq=500&rwfreq=500&rlead=250&sufs=0&order=r&whole=Whole_site&entire_flag=1&is_summary_on=1&ResultCount=10&query=inverse+radon+transform

：）

Radon 变换是平行束对图像的线积分，根据各个角度得到的一系列投影值逆radon重建得到原始图像。变换角度默认是逆时针。r=radon(im,30)得到的是一维数组，平行束与X轴夹角为30度时，距原点不同距离的投影线（平行束）上对图像的线积分。[R,Xp] = RADON(...) XP对应平行束的位置。

---