Radon变换源代码疑惑

lucklizhan 2009-04-02 04:28:15

小弟在做关于图像倾斜校正的功能
以前从来没做过图像处理方面的研究

在网上找到的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 *) mxCalloc(2*M, sizeof(double));
xTable = (double *) mxCalloc(2*N, sizeof(double));
xCosTable = (double *) mxCalloc(2*N, sizeof(double));
ySinTable = (double *) mxCalloc(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-2; k >=0; k--)
yTable[k] = yTable[k+1] + 0.5;

xTable[0] = -xOrigin - 0.25;
for (k = 1; k < 2*N; k++)
xTable[k] = xTable[k-1] + 0.5;

for (k = 0; k < numAngles; k++) {
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 = 0; p < 2*N; p++)
xCosTable[p] = xTable[p] * cosine - rFirst;
for (p = 0; p < 2*M; p++)
ySinTable[p] = yTable[p] * sine;

/* Remember that n and m will each change twice as fast as the */
/* pixel pointer should change. */
for (n = 0; n < 2*N; n++) {
pixelPtr = iPtr + (n/2)*M;
for (m = 0; m < 2*M; m++) {
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++;
}
}
}

mxFree((void *) yTable);
mxFree((void *) xTable);
mxFree((void *) xCosTable);
mxFree((void *) ySinTable);
}

...全文