Автор работы: Пользователь скрыл имя, 16 Января 2013 в 18:29, реферат
Алгебраические алгоритмы реконструкции образуют большую группу итерационных алгоритмов, в которых задача с самого начала решается в дискретной форме.
Накинем на область реконструкции прямоугольную сетку, содержащую J ячеек (элизов) таким образом, чтобы она покрывала реконструируемое изображение. На рисунке 1 показана сетка, содержащая J = 6·6 = 36 элизов.
Рисунок 1. Элемент проек
1.
Теоретическая часть
2
1.1
Метод ART (algebraic reconstruction technique)
2
1.2
Написание расширений MATLAB на языке C
4
2.
Практическая часть
6
3.
Приложение А. Программный код
mexErrMsgTxt("Input STEP must be in the range from 1 to 179.");
if (!(mxGetScalar(prhs[3])>=0 && mxGetScalar(prhs[3])<=179))
mexErrMsgTxt("Input FINISH must be in the range from 0 to 179.");
if (mxGetScalar(prhs[1])>
mexErrMsgTxt("Input BEGIN must be less than FINISH.");
// считываем исходный массив и его размерность
nrows = mxGetM(prhs[0]);
ncols = mxGetN(prhs[0]);
a = mxGetPr(prhs[0]);
// k = nrows; // т.к. матрица квадратная, кол-во лучей равно кол-ву строк/столбцов
// считываем значения начало, шаг и конец
begin = mxGetScalar(prhs[1]);
step = mxGetScalar(prhs[2]);
finish = mxGetScalar(prhs[3]);
// считаем количество шагов
n = (finish - begin)/step + 1;
// считываем количество лучей
ray = mxGetScalar(prhs[4]); // old k
// создаем выходной массив
plhs[0] = mxCreateDoubleMatrix(n, ray, mxREAL);
b = mxGetPr(plhs[0]);
x=begin;
y=step; // в цикле проходим n шагов каждый раз под разным углом x
run(a,b,nrows,ncols,n,ray,x,y)
return;
}
Листинг кода функции iART:
/*============================
*
* iART.C
*
* The calling syntax is:
*
* C = iart(B, begin, step, finish, count, size)
* B - input array, result of function art()
* begin - start angle
* step - step of angle
* finish - end angle
* count - count of iterations
* size - size of output array
*
*
*=============================
#include <math.h>
#include "mex.h"
#include "math.h"
static int makeid(int i, int j, int nrows) // вычисление ид для матрицы
{
return nrows*i+j; // i - номер столбца, j - номер строки
}
static double delta(int angle, int rows, int cols, int ray) // вычисление шага по оси х
{
double PI = 3.14159265;
double d = 0;
if ((angle < 90) && (angle > 0))
d = (cols + (rows/tan(angle*PI/180)))/(ray + 1);
if (angle == 0)
d = rows*pow((ray+1),(-1));
if (angle == 90)
d = cols*pow((ray+1),(-1));
return d;
}
static void copymatrix(double *y, double *z, int m, int n) // копирование данных из одной матрицы в другую
{
int i,j;
for (i=0; i<n; i++)
{
for (j=0; j<m; j++)
*(z+makeid(i,j, m)) = *(y+makeid(i,j, m));
}
}
static void firstisum0(double *b, double *c, int rows, int cols, int rowsb, int ray)
{
double indata, d;
int i, j, id;
d = delta(0, rows, cols, ray);
for (j=0; j<ray; j++)
{
id = (j + 1) * d;
indata = *(b+makeid(j,0,rowsb)) / ray;
for (i=0; i<cols; i++)
*(c+makeid(i,id,rows)) = indata;
}
}
static void firstisum90(double *b, double *c, int rows, int cols, int rowsb, int ray)
{
double indata, d;
int i, j, id;
d = delta(90, rows, cols, ray);
for (i=0; i<ray; i++)
{
id = (i + 1) * d;
indata = *(b+makeid(i,0,rowsb)) / ray;
for (j=0; j<rows; j++)
*(c+makeid(id,j,rows)) = indata;
}
}
static void firstisum(double *b, double *c, int rows, int cols, int rowsb, int colsb, int angle)
{
double PI = 3.14159265;
int nrow = 0;
double indata, count;
double ox,oy,oy1,d,tg;
int j,j1,z,w,w1;
int ang = angle;
int ray = colsb;
if (angle > 90) ang = 180 - angle;
// y = -tan() * x + oy
d = delta(ang, rows, cols, ray);
ox = d; // начальная точка по оси х
for (j=0; j<ray; j++)
{
tg = tan(ang*PI/180);
oy = tg * ox;
count = 0;
for (z=ox; z>=0; z--)
{
if (z<cols)
{
oy1 = (-1) * tg * z + oy;
w = oy1;
oy1 = (-1) * tg * (z+1) + oy;
if (oy1 < 0) w1 = 0;
else w1 = oy1;
for (j1=w1; j1<=w; j1++)
{
if (j1<rows)
{
count++;
}
}
}
}
indata = *(b+makeid(j,0,rowsb)) / count;
for (z=ox; z>=0; z--)
{
if (z<cols)
{
oy1 = (-1) * tg * z + oy;
w = oy1;
oy1 = (-1) * tg * (z+1) + oy;
if (oy1 < 0) w1 = 0;
else w1 = oy1;
for (j1=w1; j1<=w; j1++)
{
if (j1<rows)
{
if (angle < 90) *(c+makeid(z,(rows-j1-1),rows)
else *(c+makeid(z,j1,rows)) = indata;
}
}
}
}
ox = ox + d; //приращение шага по оси х
}
}
static void isum0(double *b, double *c, int rows, int cols, int rowsb, int colsb, int nrow)
{
double sum, delt, d;
int i, j;
int ray = colsb, id;
d = delta(0, rows, cols, ray);
for (j=0; j<ray; j++)
{
id = (j + 1) * d;
sum=0;
for (i=0; i<cols; i++)
sum = sum + *(c+makeid(i,id,rows));
delt = (*(b+makeid(j,nrow,rowsb)) - sum) / ray;
for (i=0; i<cols; i++)
*(c+makeid(i,id,rows)) = *(c+makeid(i,id,rows)) + delt;
}
}
static void isum90(double *b, double *c, int rows, int cols, int rowsb, int colsb, int nrow)
{
double sum, delt, d;
int i, j;
int ray = colsb, id;
d = delta(90, rows, cols, ray);
for (i=0; i<ray; i++)
{
id = (i + 1) * d;
sum=0;
for (j=0; j<rows; j++)
sum = sum + *(c+makeid(id,j,rows));
delt = (*(b+makeid(i,nrow,rowsb)) - sum) / rows;
for (j=0; j<rows; j++)
*(c+makeid(id,j,rows)) = *(c+makeid(id,j,rows)) + delt;
}
}
static void isum(double *b, double *c, int rows, int cols, int rowsb, int colsb, int nrow, int angle, double *matrix)
{
double PI = 3.14159265;
double sum,ox,oy,oy1,d,tg,delt;
int j,j1,z,w,w1,count;
int ang = angle;
int ray = colsb;
if (angle > 90) ang = 180 - angle;
// y = -tan() * x + oy
d = delta(ang, rows, cols, ray);
ox = d; // начальная точка по оси х
for (j=0; j<ray; j++)
{
tg = tan(ang*PI/180);
oy = tg * ox;
sum = 0;
count = 0;
for (z=ox; z>=0; z--)
{
if (z<cols)
{
oy1 = (-1) * tg * z + oy;
w = oy1;
oy1 = (-1) * tg * (z+1) + oy;
if (oy1 < 0) w1 = 0;
else w1 = oy1;
for (j1=w1; j1<=w; j1++)
{
if (j1<rows)
{
if (angle < 90) sum = sum + *(matrix+makeid(z,(rows-j1-1),
else sum = sum + *(matrix+makeid(z,j1,rows));
count++;
}
}
}
}
for (z=ox; z>=0; z--)
{
if (z<cols)
{
oy1 = (-1) * tg * z + oy;
w = oy1;
oy1 = (-1) * tg * (z+1) + oy;
if (oy1 < 0) w1 = 0;
else w1 = oy1;
for (j1=w1; j1<=w; j1++)
{
if (j1<rows)
{
delt = (*(b+makeid(j,nrow,rowsb)) - sum) / count;
if (angle < 90) *(c+makeid(z,(rows-j1-1),rows)
else *(c+makeid(z,j1,rows)) = *(c+makeid(z,j1,rows)) + delt;
}
}
}
}
sum = 0;
ox = ox + d; //приращение шага по оси х
}
}
static void firstrun(double *b, double *c, int nrows, int ncols, int rowsc, int colsc, int x)
{
if (x == 0)
firstisum0(b, c, rowsc, colsc, nrows, ncols);
if (x == 90)
firstisum90(b, c, rowsc, colsc, nrows, ncols);
if ((x != 0)&&(x != 90))
firstisum(b, c, rowsc, colsc, nrows, ncols, x);
}
static void run(double *b, double *c, int nrows, int ncols, int rowsc, int colsc, int n, int x, int y, int count, int ray)
{
mxArray *array;
double *matrix;
int begin, step;
int i, j, i1, x1;
double n1;
n1 = n * pow(2,(-1));
begin = x;
step = y;
array = mxCreateDoubleMatrix(rowsc, colsc, mxREAL);
matrix = mxGetPr(array);
firstrun(b,c,nrows,ncols,
for (j=1; j<=count; j++)
{
x = begin;
y = step;
i1 = n1;
x1 = begin + i1 * step;
for (i=0; i<n1; i++)
{
copymatrix(c,matrix,rowsc,
if (x == 0)
isum0(b, c, rowsc, colsc, nrows, ncols, i);
if (x == 90)
isum90(b, c, rowsc, colsc, nrows, ncols, i);
if ((x != 0)&&(x != 90))
isum(b, c, rowsc, colsc, nrows, ncols, i, x, matrix);
copymatrix(c,matrix,rowsc,
if (x1 == 0)
isum0(b, c, rowsc, colsc, nrows, ncols, i1);
if (x1 == 90)
isum90(b, c, rowsc, colsc, nrows, ncols, i1);
if ((x1 != 0)&&(x1 != 90))
isum(b, c, rowsc, colsc, nrows, ncols, i1, x1, matrix);
x=x+y;
x1=x1+y;
i1++;
}
}
/*
x = begin;
y = step;
i = 0;
for (j=0; j<count; j++)
{
if (j == n) i = 0;
copymatrix(c,matrix,rowsc,
if (x == 0)
isum0(b, c, rowsc, colsc, nrows, ncols, i);
if (x == 90)
isum90(b, c, rowsc, colsc, nrows, ncols, i);
if ((x != 0)&&(x != 90))
isum(b, c, rowsc, colsc, nrows, ncols, i, x, matrix);
x=x+y;
if (x>179) x = begin;
i++;
} */
}
void mexFunction( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
{
double *b, *c; // b - исходный массив (art, луч суммы), c - выходной массив (типа phantom)
double nrows, ncols; // и его размерность
int begin, step, finish; // угол, начало шаг и конец
int n; // количество шагов
int ray; // количество лучей
int x,y; // копия для begin и step
int count; // количество итераций
int rowsc, colsc; // размерность матрицы C
// проверяем корректность полученных данных
if (nrhs!=6)
mexErrMsgTxt("Six inputs required.");
if (nlhs!=1)
mexErrMsgTxt("One output required.");
if (!(mxIsDouble(prhs[0])))
mexErrMsgTxt("Input array must be of type double.");
if (!(mxGetScalar(prhs[1])>=0 && mxGetScalar(prhs[1])<=179))
mexErrMsgTxt("Input BEGIN must be in the range from 0 to 179.");
if (!(mxGetScalar(prhs[2])>0 && mxGetScalar(prhs[2])<=179))
mexErrMsgTxt("Input STEP must be in the range from 1 to 179.");
if (!(mxGetScalar(prhs[3])>=0 && mxGetScalar(prhs[3])<=179))
mexErrMsgTxt("Input FINISH must be in the range from 0 to 179.");
if (mxGetScalar(prhs[1])>
mexErrMsgTxt("Input BEGIN must be less than FINISH.");
if (mxGetScalar(prhs[4])<=0)
mexErrMsgTxt("Input COUNT must be positive.");
// считываем исходный массив и его размерность
nrows = mxGetM(prhs[0]);
ncols = mxGetN(prhs[0]);
b = mxGetPr(prhs[0]);
ray = ncols; // кол-во лучей равно кол-ву столбцов
// считываем значения начало, шаг и конец
begin = mxGetScalar(prhs[1]);
step = mxGetScalar(prhs[2]);
finish = mxGetScalar(prhs[3]);
// считаем количество шагов
n = (finish - begin)/step + 1;
// n должно быть ровно количеству строк nrows если не так то нас обманывают
if (n!=mxGetM(prhs[0]))
mexErrMsgTxt("Inputs BEGIN, STEP and FINISH must be correct for input array A.");
// считываем количество итераций пересчета
count = mxGetScalar(prhs[4]);
// считываем размерность выходного массива
rowsc = mxGetScalar(prhs[5]);
colsc = rowsc;
// создаем выходной массив
plhs[0] = mxCreateDoubleMatrix(rowsc, colsc, mxREAL);
c = mxGetPr(plhs[0]);
x=begin;
y=step;
run(b, c, nrows, ncols, rowsc, colsc, n, x, y, count, ray);
return;
}