Реализация метода ART

Министерство  образования и науки Российской Федерации

Федеральное государственное  бюджетное образовательное учреждение

высшего профессионального  образования

«Владимирский государственный университет

имени Александра Григорьевича и Николая Григорьевича Столетовых»

(ВлГУ)

Кафедра «Физики и прикладной математики»

 

 

 

Курсовой проект

по дисциплине «Компьютерное моделирование»

на тему «Реализация  метода ART»

 

 

 

 

 

Выполнил:

ст. гр. АИСу-110

Шахов Ю.А.

 

Принял:

Заякин А.А..

 

 

 

 

 

Владимир, 2012

 

СОДЕРЖАНИЕ

1.	Теоретическая часть		2
	1.1	Метод ART (algebraic reconstruction technique)	2
	1.2	Написание расширений MATLAB на языке C	4
2.	Практическая часть		6
3.	Приложение А. Программный  код		8

 

 

 

1. Теоретическая часть

1.1 Метод ART (algebraic reconstruction technique)

Алгебраические алгоритмы  реконструкции образуют большую  группу итерационных  алгоритмов,  в  которых  задача с самого  начала  решается  в дискретной форме.

Накинем на область реконструкции прямоугольную сетку, содержащую J ячеек (элизов) таким образом, чтобы она покрывала реконструируемое изображение. На рисунке 1 показана сетка, содержащая J = 6·6 = 36 элизов.

Рисунок 1. Элемент проекционной матрицы R_ij равен значению j-го элиза, который пересекает i-ый луч.

Пронумеруем элизы слева направо и сверху вниз. Определим вектор изображения X размерностью J, компоненты которого X_j равны среднему значению функции источника в j-м элизе (j = 1…J). Пусть в томографическом эксперименте выполнены измерения луч-сумм Ф ρ_n, θ_m для фиксированного набора прицельных расстояний {ρ_n | n = 1…N} и ракурсов {θ_m | m = 1…M}. Пронумеруем луч-суммы и представим эти данные также в виде вектора Y с размерностью I, равной общему количеству лучей (I = N·M). Тогда интегральные уравнения можно заменить их дискретным аналогом в матричной форме R·X=Y, где R – проекционная матрица, содержащая I строк и J столбцов, каждый элемент которой R_ij равен значению j-го элиза, который пересекает i-ый луч.

Согласно методу ART нужно:

а) выбрать некоторое  начальное приближение вектора дискретной функции источника Х⁽⁰⁾;

б) на k-ой итерации взять одну из луч-сумм из всего набора проекционных данных, рассчитать элементы соответствующей строки проекционной матрицы R и записать их ненулевые значения в вектор SR^(k). Компоненты этого вектора – значения плотности в тех элизах, которые пересекает рассматриваемый луч;

в) если рассматриваемый луч пересекает хотя бы один элиз, то на k-ой итерации произвести коррекцию дискретной функции источника Х⁽⁰⁾ по формуле

;

г) если не достигнута приемлемая точность реконструкции функции  источника, перейти к следующей  итерации. Рекомендуется перебирать луч-суммы в циклической последовательности.

Таким образом, в методе ART на k-й итерации рассматривается  только один луч, а корректируются лишь те компоненты вектора Х(k), которым соответствуют элизы, пересекаемые данным лучом, остальная часть изображения остается без изменений. Согласно формуле изменение значений компонентов вектора X^(k) пропорционально разности , то есть ошибке в одном из линейных уравнений системы перед k-ой итерацией, и пропорционально значению плотности элиза пересекаемого рассматриваемым лучом.

Учет априорной информации

Эффективность итерационных процедур при реконструкции изображений часто можно повысить посредством введения некоторых операций обработки векторов изображений между различными итерациями. Указанные операции в литературе получили название искусственных приемов, а если в них учитывается априорная информация о структуре функции источника, – дескриптивной регуляризации. В целом эти приемы повышают устойчивость решения к случайным изменениям проекционных данных.

Выбор  небольшого  значения  релаксационного  параметра  оказывает положительное  влияние  на  качество  реконструкции,  уменьшает  влияние погрешностей в уравнениях и предупреждает появление ложных фрагментов в изображении.

В реализуемом методе добавляется различное значение релаксационного параметра на каждой итерации. Оно зависит от разницы  исходной луч-суммы и текущей луч-суммы для одного и того же луча. Таким образом формула принимает следующий вид

.

 

1.2 Написание расширений MATLAB на языке C

В документации по системе MATLAB для подобных расширений употребляется  термин MEX-файл (Matlab EXtension), и сама MATLAB по этому расширению имени файла может определить, что данный модуль является ее расширением.

В модулях-расширениях MATLAB (MEX-файлах) для обмена параметрами  всех типов с вычислительной средой MATLAB используется ровно одна структура, хотя и чрезвычайно гибкая. Называется эта структура - mxArray (Matlab Extension Array).

Все типы данных MATLAB - массивы, скаляры, строки, клеточные и многомерные  массивы, объекты и т.п. выражаются при помощи этой единственной структуры.

Программный интерфейс создаваемой динамической библиотеки тоже достаточно прост. Должна быть обязательно экспортирована единственная функция с двумя параметрами - входным массивом структур mxArray и выходным массивом структур mxArray. Соответственно, задача модуля расширения заключается в том, чтобы на основе входной информации создать выходной массив.

Прототип интерфейсной функции объявлен в файле следующим  образом:

void mexFunction(    

int nlhs, /* количество выходных параметров */    

mxArray *plhs[],/*массив указателей на вых.параметры*/    

int nrhs, /* количество входных параметров */              

/**/    

const mxArray *prhs[]/*массив указ-ей на вход.пар-ры */

);

Отметим, что в английской компьютерной литературе входные и  выходные параметры называются "left hand side parameters" (параметры слева от знака присваивания) и "right hand side parameters" (параметры справа от знака присваивания). Сокращения этих названий вошли в имена параметров интерфейсной функции.

Матрицы (двумерные массивы) являются подмножеством mxArray, для чего в этой структуре предусмотрены поля pr (вещественная часть) и pi (мнимая часть). Каждое из этих полей представляет собой одномерный массив, содержащий элементы матриц (double-числа) поколоночно - сначала все элементы первого столбца матрицы, затем - второго и т.д.

Так как извлечение требуемого типа данных для обработки напрямую из общей структуры может быть сложным (из-за наличия разного рода флажков, типов, размерностей и т.п.), то MATLAB содержит специальную библиотеку вспомогательных функций, облегчающую работу с тем или иным типом данных внутри расширения MATLAB. Все функции этой библиотеки, работающие с массивом mxArray, имеют префикс mx (напр. mxGetPr()).

После извлечения данных из массива типа  mxArray все операции над переменными пишутся на языке C.

 

2. Практическая часть

Метод ART реализован двумя  функциями: art и iart.

В качестве входных данных функция ART принимает четыре параметра:

• массив A, представляющий из себя результат функции phantom;
• начальный угол, с которого начинается просвечивание объекта лучами;
• шаг, на который изменяется угол просвечивания объекта;
• конечный угол, которым заканчивают просвечивание объекта;
• количество лучей.

Выходным значением  функции ART является матрица лучсумм.

Функция iART предназначена построение исходного изображения по полученным лучсуммам. Входные данные для функции iART:

• массив B, представляющий из себя результат функции ART;
• начальный угол, с которого начинается просвечивание объекта лучами;
• шаг, на который изменяется угол просвечивания объекта;
• конечный угол, которым заканчивают просвечивание объекта;
• количество итераций для просчета результата;
• размерность выходного массива.

Выходным значением  является матрица с восстановленным  изображением.

Для запуска функций  необходимо выполнить следующие шаги:

1. скопировать файлы art.c и iart.c в папку с MATLAB;
2. выполнить команды “mex art.c” и “mex iart.c”;

Синтаксис запуска функций:

B = art(A, begin, step, finish, ray) – функция ART.

C = iart(B, begin, step, finish, count, size) – функция iART.

Результат выполнения функций представлен на рисунке.

A = phantom();

B = art(A,0,1,179,256);

C = iart(B,0,1,179,1,256);

B1 = art(A,0,1,179,512);

C1 = iart(B1,0,1,179,1,256);

C2 = iart(B1,0,1,179,50,256);

C3 = iart(B1,0,1,179,300,256);

subplot(2,4,1), imshow(A), title('Original')

subplot(2,4,2), imshow(A), title('Original')

subplot(2,4,3), imshow(A), title('Original')

subplot(2,4,4), imshow(A), title('Original')

subplot(2,4,5), imshow(C), title('iART 256')

subplot(2,4,6), imshow(C1), title('iART 512')

subplot(2,4,7), imshow(C2), title('iART 512 iteration 50')

subplot(2,4,8), imshow(C3), title('iART 512 iteration 300')

 

Приложение А

Листинг программного кода.

Листинг кода функции ART:

/*=================================================================

*

* ART.C   

 *

* The calling syntax is:

*

*      B = art(A, begin, step, finish, ray)

*         

*          A - input array, result of function phantom()

*          begin - start angle

*          step - step of angle

*          finish - end angle

*          ray - count of rays

*

*

*=================================================================*/

 

#include <math.h>

#include "mex.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 sum0(double *a, double *b, int rows, int cols, int rowsb, int ray, int nrow)

{

    double sum, d = 0;

    int i, j, id = 0;

    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 + *(a+makeid(i,id,rows));

        *(b+makeid(j,nrow,rowsb)) = sum;

    }

}

 

static void sum90(double *a, double *b, int rows, int cols, int rowsb, int ray, int nrow)

{

    double sum, d;

    int i, j, 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 + *(a+makeid(id,j,rows));

        *(b+makeid(i,nrow,rowsb)) = sum;

    }

}

 

static void sum(double *a, double *b, int rows, int cols, int rowsb, int ray, int nrow, int angle)

{

    double PI = 3.14159265;

    double sum,ox,oy,oy1,d,tg;

    int j,j1,z,w,w1;

    int ang = angle;

    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;

        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 + *(a+makeid(z,(rows-j1-1),rows));

                        else sum = sum + *(a+makeid(z,j1,rows));

                    }

            }

        }

        *(b+makeid(j,nrow,rowsb)) = sum;

        sum = 0;

        ox = ox + d; //приращение шага по оси х

    }

}

 

static void run(double *a, double *b, int nrows, int ncols, int n, int ray, int x, int y)

{

    int i;

    for (i=0; i<n; i++)

    {

        if (x == 0)

            sum0(a, b, nrows, ncols, n, ray, i);

        if (x == 90)

            sum90(a, b, nrows, ncols, n, ray, i);

        if ((x != 0)&&(x != 90))

            sum(a, b, nrows, ncols, n, ray, i, x);

        x=x+y;

    }

}

 

void mexFunction( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )   

{

    double *a, *b; // a - исходный массив (phantom), b - выходной массив (луч суммы)

    double nrows, ncols; // и его размерность

    //double d, ox, oy, oy1, tg; // d - в нее сохраняем дельту, ox - точка по оси х, oy - аналогично

    int begin, step, finish; // угол, начало шаг и конец

    int n; // количество шагов

    int x,y; // копия для begin и step

    int ray; // количество лучей

   

    // проверяем корректность полученных данных

    if (nrhs!=5)

        mexErrMsgTxt("Five inputs required.");

    if (nlhs!=1)

        mexErrMsgTxt("One output required.");

    if (!(mxIsDouble(prhs[0])))

        mexErrMsgTxt("Input array must be of type double.");

    if (mxGetM(prhs[0])!=mxGetN(prhs[0]))

        mexErrMsgTxt("Input array must be square.");

    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))