САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
Математико-механический факультет
Кафедра информатики
Кудрин Богдан Константинович
Выпускная квалификационная работа
Разработка и реализация алгоритмов моделирования и
визуализации низкочастотного магнитного поля в устройствах
магнитотерапии
Допущена к защите.
Зав. Кафедрой:
д. ф.-м. н., профессор Новиков Б. А.
Научный руководитель:
к.ф.-м.н., доцент Ампилова Н.Б.
Рецензент:
к.ф.-м.н.,Раба Н.О.
Санкт-Петербург
2016
SAINT-PETERSBURG STATE UNIVERSITY
Mathematics&Mechanics Faculty
Chair of Computer science
Kudrin Bogdan
Graduate qualifying paper
Design and implementation of the algorithms for modelling and visualization of
low-frequency magnetic field in magnetotherapy devices
Admitted for defence.
Head of the chair:
Dr. of Phys. and Math. Sci., professor Novikov B. A.
Scientific supervisor:
Ph.D., associate professor Ampilova N.B.
Reviewer:
Ph.D., Raba N.O.
Saint-Petersburg
2016
2
Оглавлени
ВВЕДЕНИЕ....................................................................................................5
ГЛАВА 1. ДВУМЕРНЫЙ СЛУЧАЙ. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ,
РЕАЛИЗАЦИЯ, ВИЗУАЛИЗАЦИЯ....................................................................... 9
1.1. Математическая модель магнитного поля цилиндрической
катушки индуктивности....................................................................................9
1.2. Двумерный случай. Результаты..................................................13
ГЛАВА 2. ТРЕХМЕРНЫЙ СЛУЧАЙ. ВЫЧИСЛЕНИЕ И
ВИЗУАЛИЗАЦИЯ..................................................................................................16
2.1. Математическая модель магнитного поля нескольких
цилиндрических катушек индуктивности в 3D............................................16
2.2. Описание практической задачи и результат визуализации........19
ГЛАВА 3. Интерполяция данных..............................................................23
3.1 Алгоритм интерполяции вычисленных данных........................23
3.1.1. Математическая модель............................................................24
3.1.2 Алгоритм интерполяции............................................................26
3.1.3. Анализ результатов...................................................................26
3.2 Алгоритм интерполяции экспериментальных данных.............28
3.2.1 Описание задачи.........................................................................28
3.2.2 Математическая модель.............................................................29
3.2.3 Визуализация.............................................................................. 30
3.2.4 Алгоритмы интерполяции и визуализации..............................31
3.2.5 Описание интерфейса................................................................ 31
3.2.6 Пример использования приложения.........................................32
3
ГЛАВА 4. Оптимизация визуализации.....................................................33
Заключение.................................................................................................. 37
Литература................................................................................................... 38
Приложение 1. Описание пользовательского интерфейса вкладки
Вычисление............................................................................................................40
Приложение 2. Описание работы с пакетом ParaView............................41
1.
Описание формата файла .vtk.....................................................41
2.
Описание работы с ParaView.......................................................42
4
ВВЕДЕНИЕ
Объективными природными факторами, окружающими все биообъекты
и непрерывно влияющими на них являются физические поля: магнитные,
электрические, электромагнитные, акустические, тепловые, гравитационные.
Многочисленные экспериментальные факты показывают, что физические
поля очень низкой интенсивности влияют на функционирование живых
систем [1,2]. Изучением этих явлений занимается гелиобиология, а ее
составная часть — магнитобиология — изучает эффекты слабых
низкочастотных переменных магнитных и электрических полей. Как указано
в [19], в настоящее время электронный банк данных по всем аспектам
маг н и тоби ол ог и и н асчитывает более 30000 ссылок. Э фф екты
электромагнитных полей низкой интенсивности достоверно наблюдаются не
только в медицине, но и в биологических экспериментах с бактериями,
насекомыми, млекопитающими.
Магнитотерапия представляет собой лечебный метод, основанный на
использовании низкочастотных магнитных полей. На сегодняшний день в
мире существует несколько десятков успешно функционирующих
магнитотерапевтических аппаратов. Такая терапия хорошо сочетается со
средствами традиционной медицины.
В то же время пока нет единых правил подбора параметров в
различных устройствах магнитотерапии. Это связано прежде всего тем, что
механизм действия малых доз как физических полей так и лекарственных
препаратов на живые системы весьма сложен и недостаточно изучен. Этому
вопросу посвящены работы А.Чижевского, Е.Бурлаковой, В.Самойлова,
А.Коновалова и многих других. Достаточно подробный обзор современного
состояния исследований в этой области дан в работах Л.Н. Галль [4,5].
Один из перспективных подходов к выбору парамет ров,
обеспечивающий наибольшую лечебную эффективность, предложен в работе
[3]. Автор основывает свою методику на фундаментальных положениях
биофизики сенсорных систем, которые описывают отклик биосистемы на
5
внешние воздействия, а также предлагает метод контроля за состоянием
здоровья пациента.
Анализ работ в этой области показывает, что пока определяющим
моментом в изучении действия слабого магнитного поля на живой организм
является эксперимент. Тем не менее, для определенного класса устройств
магнитотерапии можно рассмотреть некоторую
математическую модель,
связанную с генерацией магнитного поля катушками индуктивности. Метод
графического моделирования магнитного поля
соосных катушек
индуктивности в плоскости сечения, проходящей через их общую ось, был
применен в работе [11] и показал хорошее совпадение результатов измерений
и поля, полученного в результате вычислений.
Весьма актуальной
является более общая задача — построение
изображения распределения магнитного поля, генерируемого несколькими
катушками, расположенными в трехмерном пространстве достаточно
произвольным образом. Врач получает графическое изображение поля.
Изменение положения катушек, расстояния между ними, силы тока позволяет
подобрать режим, при котором может быть достигнут наилучший лечебный
эффект при уменьшении продолжительности курса лечения.
Постановка задачи.
Основные цели работы — разработать приложение, решающее задачу
моделирования и визуализации в трехмерном пространстве магнитного поля,
возникающего в устройствах магнитотерапии, и оптимизировать процесс
визуализации полученных результатов.
Для этого были поставлены следующие задачи:
1) разработка и реализация алгоритма интерполяции данных,
полученных в результате моделирования электромагнитного поля с
целью оптимизации времени вычислений
2) разработка и реализация алгоритма интерполяции экспериментально
полученных данных (когда данных недостаточно для трехмерной
интерполяции)
6
3) оптимизация процесса визуализации данных, полученных при
моделировании результатов.
Структура работы.
Первая глава содержит:
1) Описание математической модели магнитного поля, создаваемого
одним витком цилиндрической катушки;
2) описание алгоритма вычисления значений электромагнитного поля
в узлах двумерной сетки для двух случаев:
- осевое сечение одной катушки
- плоское сечение области с несколькими активными катушками
(оси катушек параллельны).
Вторая глава содержит:
1) описание алгоритма вычисления значений электромагнитного поля в
у з л а х т р ех м е р н о й с е т к и д л я н е с кол ь к и х п р о и з вол ь н о
сконфигурированных катушек (базовый алгоритм);
2) пример применения приложения для решения одной практической
задачи, возникающей в магнитотерапии.
Третья глава содержит:
1) описание алгоритма интерполяции данных, полученных в результате
работы базового алгоритма;
2) описание алгоритма интерполяции экспериментально полученных
данных;
В четвертой главе описывается оптимизация процесса визуализации
данных, полученных при моделировании результатов.
Приложение 1 содержит описание работы с графическим интерфейсом
разработанного приложения, а Приложение 2 содержит описание работы с
пакетом визуализации ParaView.
Результаты работы докладывались на конференциях «Технологии
Microsoft в теории и практике программирования 2013» [6] и международной
7
конференции «Герценовские чтения – 2013»[7], а также на конференциях
CEMA 2013[12,13,14], CEMA 2014[15,16,17], CEMA 2015[18].
ГЛАВА 1. ДВУМЕРНЫЙ СЛУЧАЙ. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ,
РЕАЛИЗАЦИЯ, ВИЗУАЛИЗАЦИЯ.
1.1.
Математическая модель магнитного поля цилиндрической
катушки индуктивности.
На рис.1.1
показана схема, используемая при исследовании
пространственной конфигурации магнитного поля, генерируемого
цилиндрической катушкой. Введены следующие обозначения:
8
i - величина тока в рассматриваемом единичном контуре;
R - радиус контура;
О - центр контура и начало цилиндрической системы координат;
М -произвольная точка, в которой вычисляется магнитная индукция
⃗
BM
m - проекция точки М на плоскость контура;
O1 -центр окружности, лежащей в плоскости, параллельной плоскости
контура, и проходящей через точку М;
ρ - радиус этой окружности;
dl
2α
- элементарный участок контура;
- угол между двумя симетрично расположенными элементарными
участками контура;
r - расстояние от концов элементарных участков dl контура до точки
М;
⃗
Bz ,
⃗
B ρ -компоненты
BM
магнитной индукции ⃗
в точке М по
соответсвующим осям z, ρ в цилиндрической координатной системе.
Рис. 1.1 Схема витка цилиндрической катушки индуктивности.
B M в некоторой точке М в
При вычислении магнитной индукция ⃗
пространстве нужно применить общий метод для ее представления. В силу
9
симметрии величина векторного потенциала ⃗
A
не зависит от угл а 2 α
между радиусами, проходящими через два симметрично расположенных по
контуру элемента dl . Для упрощения представления и без умаления
общности точка М выбирается в плоскости xOz, где α=0 . Тогда элементы
dl , на которые можно разложить контур, группируются в пары элементов,
равноотстоящих от М и с координатами α и −α .
Вектор-потенциал ⃗
A
при использовании
в точке М перпендикулярен плоскости xOz и
цилиндрической системы
положить его компоненты равными:
A =A α ,
координат мы можем
A ρ =0
и A z =0 . Используя
формулу для векторного потенциала ⃗
A получаем:
⃗
A =⃗
Aα =
μ0 i ⃗
dlα
∮
4π r r
(1.1.1)
Из рис. 1.1а и 1.1б видно, что:
2
2
2
r = R + ρ −2 ρR cos α+ z
(1.1.2)
2
Подставляя (1.1.2) в (1.1.1) получаем:
π
μi
R cos α dl
A= A α = 0 ∫
2
2 π 0 √ R + ρ2−2 ρR cos α+ z 2
Сделаем подстановку
(1.1.3)
с введением новой переменной β,
α= π +2 β
тогда (1.1.3) примет вид:
π /2
A= A α =
2
μ0 iR
2 sin β−1
dβ
∫
π 0 √ (R+ ρ)2+ 4 ρR sin2 β
(1.1.4)
Знаменатель в подынтегральной функции выражения (1.1.4) можно
представить в виде:
r=2
√ ρR √1−k 2 sin2 β , при k 2= 4 ρR
(R+ ρ)2 + z 2
k
(1.1.5)
Принимая во внимание, что
2
2 sin β−1
√ 1−k2 sin 2 β
1
2−k
2 √ 1−k sin β
−
2
2
2
k √ 1−k sin β
k2
2
=
2
2
(1.1.6)
можно написать:
A=A α =
μ0 i √ R 2−k 2
2
(
K − L)
k
k
2π √ρ
(1.1.7)
10
где:
π /2
K =∫
0
dβ
(1.1.8)
√ 1−k2 sin2 β
π /2
L= ∫ √ 1−k 2 sin2 β dβ
(1.1.9)
0
Формулы (1.1.8) и (1.1.9) представляют собой полные эллиптические
интегралы первого и второго рода как функции величины k. Чтобы измерить
векторный потенциал ⃗
A
необходимо вычислить k при заданных значениях
z и ρ. Затем по таблицам или графикам найти соответствующие значения K и
L и вычислить функцию:
2
F (k )=
2−k
2
K− L
k
k
(1.1.10)
Затем вычисляется вектор-потенциал ⃗
A . В силу симметрии
относительно оси Oz линии вектор-потенциала имеют форму окружностей,
центры которых лежат на данной оси, а плоскости кругов параллельны
плоскости, ограниченной контуром.
Силовые линии магнитного поля лежат в плоскости, проходящей через
о с ь Oz. Вектор магнитной индукция определяется двумя скалярными
со ставляющими B ρ =rot ρ A
и B z =rot z A ,
так
к а к B α=0 .
В
цилиндрической системе координат получаются следующие формулы для
B M в точке М:
компонент вектора магнитной индукции ⃗
Bρ=
2
2
2
μ0i
z
R +ρ +z
(
L−K )
2 π ρ √ ( R+ ρ )2 + z 2 ( R −ρ )2 +z 2
Bz=
2
2
2
μ0 i
1
R − ρ −z
(
L+ K)
2 π √ ( R+ ρ )2+ z 2 ( R− ρ )2 + z 2
(1.1.11)
(1.1.12)
Для вычисления магнитного поля по формулам (1.1.11. и 1.1.12.)
используются значения функций K и L. Так как данные полные
эллиптические интегралы не могут быть выражены в элементарных
функциях, используется интерполяция. В качестве приближенных значений в
данной работе используются линейно интерполированные значения.
11
Табличные значения функций в промежуточных точках и значения конечных
разностей первого порядка взяты из справочника [10].
Также в математическом пакете Maple была разработана программа,
использовавшаяся для проверки значений, полученных с помощью
интерполяции. Сопоставление результатов вычислений показало, что
выбранная линейная интерполяция обеспечивает хорошую точность
получаемых значений.
1.2.
Двумерный случай. Результаты.
Под двумерным случаем подразумевается расчет значений
электромагнитного поля для набора точек,
лежащих в одной плоскости
(плоская сетка).
В качестве прототипа системы была разработана программа на языке
C#, с использованием средств визуализации .NET Framework. Программа
производит расчет значений модулей векторов векторного поля для
прямоугольной участка, являющегося сечением исследуемой области.
Расчет и визуализация производились для двух следующих случаев:
1)
Магнитное поле внутри одной катушки (по полученной карте
цветов проверялась корректность работы программы – карта
цветов визуально сравнивалась с результатами, описанными в
[11]). Значения векторного поля вычислялись по принципу
12
суперпозиции по всем виткам катушки (двойной цикл - по виткам
и по намоткам в каждом витке).
Результат работы программы для данного случая (карта
цветов) изображен на рисунке 1.2.1.
Также на рисунке показан графический интерфейс программы
для данного случая, позволяющий задавать характеристики
катушки.
Увеличение магнитной индукции
Рисунок 1.2.1. Визуализация электромагнитного поля внутри
одной катушки.
2)
Магнитное поле, создаваемое несколькими катушками с
параллельными осями, нижние основания катушек при этом лежат в одной
плоскости. Расчет и визуализация производились для плоскости,
13
перпендикулярной осям одинаковых
катушек, и проходящей через их
середины. Значения скалярного поля вычислялись по принципу
суперпозиции по всем виткам катушки и по всем катушкам.
На рисунке 1.2.2 изображен результат работы программы для
данного случая (карта цветов). Направление тока в крайних
катушках совпадает, направление тока в центральной катушке
противоположно.
Увеличение магнитной индукции
Рисунок 1.2.2. Визуализация плоского сечения электромагнитного поля
(несколько активных катушек)
14
ГЛАВА 2. ТРЕХМЕРНЫЙ СЛУЧАЙ. ВЫЧИСЛЕНИЕ И
ВИЗУАЛИЗАЦИЯ.
2.1.
Математическая модель магнитного поля нескольких
цилиндрических катушек индуктивности в 3D.
Расчет пространственной конфигурации магнитного поля производится
внутри области, имеющей форму прямоугольного параллелепипеда. При этом
ячейки сетки, в узлах которой вычисляются значения, также имеют форму
прямоугольных параллелепипедов. Все ячейки равны между собой.
За начало глобальной системы координат принимается «нижний левый
ближний» угол области.
Рисунок 2.1.1. Схема витка цилиндрической катушки индуктивности.
Катушки позиционируются заданием координат центра нижнего
основания O и заданием координат вектора оси катушки (оси zloc). Катушки
15
могут лежать как вне области, так и внутри нее. Направление вектора оси
катушки выбирается в соответствии с направлением электрического тока в
катушке.
Координатами точки M в локальной цилиндрической системе координат
(ЛЦСК) являются расстояние ρ от точки M до оси катушки и расстояние z до
плоскости основания катушки. Система координат выбирается таким
образом, чтобы угловая координата была равна 0. Таким образом, при
фиксированных катушке и точке M координаты точки в ЛЦСК задаются
однозначно.
В качестве базиса локальной декартовой системы координат (ЛДСК)
каждой катушки берутся 3 взаимно ортогональных вектора: zloc, xloc и yloc. В
связи с тем, что координаты точки М в ГСК и ЛЦСК определены однозначно,
то выбор ЛДСК является произвольным (при переходе ГСК-ЛДСК-ЛЦСК
получаются одни и те же координаты ρ и z точки M, и обратно, при переводе
вектора э/м индукции путем ЛЦСК-ЛДСК-ГСК получаются одни и те же его
координаты в ГСК).
Необходимость введения промежуточной ЛДСК обусловлена удобством
перехода из одной декартовой системы координат в другую (переход
осуществляется сдвигом начала координат и умножением на фиксированную
матрицу перехода). Каждой катушке соответствует прямая и обратная
трансформационные матрицы.
Столбцами прямой трансформационной матрицы являются координаты
базиса ЛДСК катушки в базисе глобальной системы координат. Т.е. если
координаты оси Ox ЛДСК (x1,y1,z1), координаты оси Oy (x2,y2,z2), и оси Oz
(x3,y3,z3), то прямая трансформационная матрица имеет вид:
T=
x1
y1
z1
x2
y2
z2
x3
y3
z3
.
16
Переход из ЛДСК в ЛЦСК и обратно также является простым (координата
по оси z не изменяется, координаты mx и my трансформируются в m и
обратно: вектор Bp раскладывается по базису (mx, my))
Алгоритм вычисления электромагнитного поля в узлах сетки состоит в
следующем:
1. ГСК-ЛДСК
Пусть произвольная точка области имеет глобальные координаты
M = (xM,yM,zM), координаты базиса ЛДСК в ГСК v1 = (x1,y1,z1), v2 =
(x2,y2,z2), v3 = (x3,y3,z3), координаты начала ЛДСК O = (x0,y0,z0). Тогда
вектор координат точки в ЛДСК будет иметь вид:
(mx,my, mz)T = T x (M – O)
(здесь M и O – радиус-векторы соответствующих точек)
2. ЛДСК-ЛЦСК
ρ=√ m 2x +m 2y
z=m z
3. Вычисляется значение вектора электромагнитной индукции в точке как
суперпозиция значений для каждого контура катушки.
В локальной цилиндрической системе координат (ЛЦСК)
используются формулы (1.1.11. и 1.1.12.) для компонент вектора
B M в точке М.
магнитной индукции ⃗
4. Производится перевод вектора электромагнитной индукции в
локальную декартову систему координат (ЛЦСК-ЛДСК).
B z loc=B z cyl
B x loc=B ρ ×
B y loc=B ρ ×
mx
√ m +m
2
x
my
2
y
√m +m
2
x
2
y
Bloc = (Bx loc, By loc,Bz loc)
5. Производится перевод вектора электромагнитной индукции в
глобальную систему координат (локальная декартова - глобальная).
17
Bglob = T-1 x Bloc
6. Пункты 1-5 повторяются для каждой катушки и считается
суперпозиция значений электромагнитной индукции.
Программа, ре а лизующая а лгоритм вычисления значений
электромагнитного поля, создаваемого несколькими несоосными катушками,
разработана на языке программирования Java. Описание пользовательского
интерфейса данной программы приведено в Приложении 1. Визуализация
полученных результатов осуществлена с помощью средств пакета научной
визуализации ParaView.
2.2.
Описание практической задачи и результат визуализации.
Как известно [11], конфигурация катушек в магнитотерапии зависит от
медицинских показаний. На рисунке 2.2.1 изображен один из приборов,
применяемых в магнитотерапии, а на рисунке 2.2.2 изображена конфигурация
катушек, используемая в исследуемом в этой главе случае.
В этом параграфе приводится описание задачи, возникающей в этом
конкретном реальном случае. Также приведены результаты, полученные при
визуализации магнитного поля для данной конфигурации катушек.
18
Рисунок 2.2.1. Один из приборов, используемых в магнитотерапии.
На рисунке 2.2.2 изображены 6 катушек, А 1, A2, B1, B2, C1 и C2. Центры
ОА1 и ОА2 катушек А1 и А2 лежат на малом основании трапеции. Расстояние
между их осями 232 мм. Центры ОС1 и ОС2 катушек С1 и С2 лежат на большом
основании трапеции. Катушки А1 и С1 имеют общую ось, и катушки А2 и С2
тоже имеют общую ось. На рисунке указаны размеры трапеции. Центры О В1 и
ОВ2 катушек В1 и В2 лежат в серединах боковых ребер трапеции. Оси всех
катушек перпендикулярны сторонам трапеции и лежат в плоскости трапеции.
Центр О системы координат находится в нижнем левом углу трапеции.
Ось Ох направлена вдоль большого основания трапеции, ось Оу
перпендикулярна ей. Центры всех катушек в миллиметрах имеют следующие
координаты:
ОА1 (184, 500,0);
ОА2 (416,500,0);
ОВ1 (50, 250, 0);
ОВ2 (550, 250,0);
ОС1 (184,0,0);
ОС2 (416,0,0).
Внешний диаметр каждой катушки равен 116 мм, внутренний – 36 мм.
Высота катушки равняется 34 мм. Каждая катушка имеет 500 намоток.
Электрический ток в каждой катушке имеет величину 3А. Направление тока в
катушках с общими осями противоположное.
19
Рисунок 2.2.2. Схема расположения катушек индуктивности в устройстве.
Требуется выполнить компьютерную симуляцию и визуализацию магнитных
полей, создаваемых следующими 5 парами катушек:
1)
2)
3)
4)
5)
А1 и С1
А1 и С2
В1 и С1
А1 и В1
В1 и В2
Визуализация для каждой пары производится двумя способами:
1. Общий – в области, имеющий размеры 500х600х232 мм, и содержащей
все катушки. Результат данной визуализации приведен на рисунке 2.2.3.
2. Детализированный - в области, имеющей форму прямоугольного
параллелепипеда, имеющего размеры по ширине и высоте 232 мм, а
длина равняется расстоянию между центрами рассматриваемых
катушек. И центры данных катушек лежат на меньших боковых гранях
20
параллелепипеда. Результат данной визуализации приведен на рисунке
2.2.4.
На рисунках 2.2.3 и 2.2.4 синие линии – это силовые линии
электромагнитного поля, создаваемого активными катушками.
Рисунок 2.2.3. Результат визуализации для катушек А1 и С2
(глобальная область)
21
Рисунок 2.2.4. Результат визуализации для катушек А1 и С2
(детализированный).
ГЛАВА 3. Интерполяция данных.
3.1 Алгоритм интерполяции вычисленных данных.
Алгоритм, описанный в параграфе 2.1 является ресурсоемким, так как
используется принцип суперпозиции и общее магнитное поле представляется
в виде суммы магнитных полей, генерируемых отдельными контурами
катушки.
В связи с этим в работе [8] был реализован алгоритм трилинейной
интерполяции электромагнитного поля по базовым значениям, полученным в
результате вычислений по методу, реализованному ранее. Это позволило
повысить скорость вычисления без существенной потери точности.
Результаты оптимизации приведены в параграфе 3.1.3.
В данном параграфе мы рассматриваем три способа интерполяции:
С2
а) линейная по 8 ближайшим соседям (данный способ применим для
любой точки внутри куба, вершинами которого являются 8 соседних точек
базовой сетки);
22
б) квадратичная интерполяция Лагранжа по трем точкам, находящимся
с точкой интерполяции на одной диагонали (точка интерполяции лежит точно
в середине куба с вершинами в соседних точках базовой сетки).
Эксперименты показывают, что любой из методов приводит к
сокращению времени вычислений без существенной потери точности.
3.1.1. Математическая модель.
Трилинейная интерполяция
В качестве интерполяционных узлов используются 8 ближайших узлов
базовой сетки: { (x i , y j , zk ) }. Интерполяционные формулы напоминают
формулы для интерполяции Лагранжа, но фактически являются формулами
для линейной трехмерной интерполяции:
23
l nmp ( x i , y j , z k )=0 , i≠n∨ j≠m∨k≠p
l nmp ( xn , y m , z p )=1
1
1
1
( x−x )( y− y )( z−z )
l nmp ( x , y , z )= ∏ ∏
∏ x −x i y − yj z −zk
i=0 , i≠n j=0 , j≠m k=0 ,k ≠ p ( n
i )( m
j )( p
k)
,
,где lnmp – базисный полином для узла (xn,ym,zp)
И, соответственно, интерполяционный многочлен будет иметь вид:
1
L( x , y , z )= ∑
1
1
∑ ∑ f nmp⋅l nmp ( x , y , z )
n=0 m=0 p=0
, где fnmp – значение, полученное непосредственным вычислением в
точке (xn,ym,zp)
Квадратичная интерполяция Лагранжа
В данном случае, в качестве узлов интерполяции используются 3 точки
базовой сетки, находящиеся с точкой интерполяции на одной диагонали:
24
l n ( xi , yi , z i )=0 ,i≠n
l i ( x i , yi , z i )=1 ,
,где li – базисный полином для узла (xi,yi,zi)
( x−x ) +( y− y ) + ( z− z )
√
∏
√( x −x ) + ( y − y ) +( z −z )
2
2
l n ( x , y , z )=
n=0 ,n≠i
2
i
i
2
n
2
i
i
2
n
i
n
2
i
И, соответственно, интерполяционный многочлен будет иметь вид:
2
L( x , y , z )= ∑ f n⋅ln ( x , y , z )
n=0
, где fn – значение, полученное непосредственным вычислением в точке
(xn,yn,zn)
Таким образом, используется одномерная квадратичная интерполяция
Лагранжа (точка и узлы интерполяции находятся на одной прямой).
3.1.2 Алгоритм интерполяции
Оба алгоритма различаются только одним шагом, зависящим от выбора
способа интерполяции.
1. Выбирается одна из координатных функций векторного поля
электромагнитной индукции, для которой будет производиться интерполяция.
2а. Для заданной точки выбирается 8 ближайших узлов базовой сетки это будут узлы интерполяции.
2б. Для заданной точки выбирается 3 узла базовой сетки, которые
лежат с данной точкой на одной диагонали - это будут узлы интерполяции.
3. Строится интерполяционный многочлен Лагранжа от трех
переменных (в первом случае интерполяция Лагранжа вырождается в
линейную), интерполирующий выбранную координатную функцию
векторного поля электромагнитной индукции.
4. Алгоритм повторяется для каждой точки заданной новой сетки.
25
3.1.3. Анализ результатов.
Анализ скорости работы программ, реализующих данные алгоритмы
проводился следующим образом: сначала производится непосредственное
вычисление трехмерного массива векторов электромагнитной индукции.
Затем строилась сетка, каждый из узлов которой является центром ячейки
базового массива, и для этой сетки производится интерполяция и
производится вычисление с помощью базового алгоритма. Таким образом,
количество узлов одинаковое, что позволяет объективно сравнить скорость и
точность данных алгоритмов.
Анализ погрешности проводился следующим образом: находились
нормы (максимум модуля) каждого из полученных трехмерных массивов:
полученного непосредственным вычислением и двух после интерполяции, а
потом считался коэффициент интерполяции – отношение нормы массива,
п о л у ч е н н о го и н т е р п о л я ц и е й к н о р м е м а с с и в а , п о л у ч е н н о го
непосредственным вычислением. Норма (максимум модуля) была выбрана,
потому что значение придавалось критическому расхождению по всему
массиву, а не в каждом конкретном узле. Как видно из таблицы 1,
погрешность составила не более 8 % для линейной интерполяции, и 12% для
квадратичной. Более высокая точность линейной в данном случае
обусловлена несимметричностью метода одномерной квадратичной
интерполяции (результат зависит от выбора диагонали и направления).
Размер сетки
Коэффициент интерполяции
Время вычисления, секунды
Число
точек
Шаг,
мм
Базовый
алгоритм
Трилинейная
интерполяция
Квадратичная
одномерная
интерполяция
Трилинейная
интерполяция
Квадратичная
одномерная
интерполяция
63602
188160
530145
2495460
10
7
5
3
6
17
49
228
1
2
8
38
<1
1
5
25
0.92
0.93
0.96
0.97
1.03
1.10
1.10
1.12
Таблица 3.1.3.1. Скорость работы алгоритмов.
26
Базовый алгоритм
м
)
то
че
к
(ш
аг
5
м
м
)
Линейная трехмерная интерполяция
53
01
45
63
60
2
18
81
60
то
чк
и
то
че
к
(ш
аг
(ш
аг
7
10
м
м
м
)
сек
50
45
40
35
30
25
20
15
10
5
0
Квадратичная одномерная интерполяция
Диаграмма 3.1.3.1. Скорость работы алгоритмов.
3.2 Алгоритм интерполяции экспериментальных данных.
3.2.1 Описание задачи
Для устройства, описанного в параграфе 2.2, существует возможность
измерять значения электромагнитной индукции в центрах катушек и в точке,
находящейся в середине соединяющей их линии. Полученные результаты
являются скалярными величинами. В каждый момент времени активными
являются ровно две катушки. Порядок включения и выключения катушек
приведен в таблице 3.2.1. Необходимо вычислить интерполированные
значения на линии, соединяющей центры активных катушек.
27
Таблица 3.2.1. Результаты измерений.
Данные для интерполяции записываются в текстовый файл, в котором
13 столбцов и строк столько, сколько пар катушек использовалось для
получения экспериментальных данных.
3.2.2 Математическая модель
Интерполяция проводится для точки, которая лежит с узлами
интерполяции на одной прямой, что позволяет применить алгоритм
одномерной квадратичной интерполяции Лагранжа. В качестве узлов для
интерполяции используются 3 точки:
Центр первой активной катушки.
Координаты этого узла задаются в столбцах 3, 4, 5 входного файла.
Обозначим их ( x 0 , y0 , z 0 ) .
28
Экспериментально полученное значение в этом узле записывается в 10
столбце входного файла. Обозначим его f 0 .
Центр второй активной катушки.
Координаты этого узла задаются в столбцах 7, 8, 9 входного файла.
Обозначим их ( x 2 , y 2 , z 2 ) .
Экспериментально полученное значение в этом узле записывается в 11
столбце входного файла. Обозначим его f 2 .
Середина отрезка, соединяющего точки 1 и 2.
Координаты этого узла вычисляются из координат точек 1 и 2.
Обозначим их ( x 1 , y 1 , z 1 ) . Тогда
( x 1 , y 1 , z 1 )=(
x0 +x 2 y 0 + y 2 z 0 + z 2
,
,
)
2
2
2
.
Экспериментально полученное значение в этом узле записывается в 12
столбце входного файла. Обозначим его f 1 .
Для точки интерполяции (x,y,z) рассмотрим полином
( x−x ) +( y− y ) + ( z− z )
√
∏
√( x − x ) + ( y − y ) +( z − z )
2
2
l n ( x , y , z )=
2
i
n=0 , n≠i
i
2
n
2
i
i
2
n
i
n
2
i
,
l n ( x i , yi , z i )=0 ,i≠n
l i ( xi , yi , z i )=1 ,
где li – базисный полином для узла (xi,yi,zi).
Интерполяционный многочлен будет иметь вид:
2
L( x , y , z )= ∑ f n⋅ln ( x , y , z )
n=0
где fn – значение, полученное измерением в точке (xn,yn,zn).
29
Для получения интерполированного значения в точке необходимо лишь
вычислить значение интерполяционного многочлена в этой точке.
3.2.3 Визуализация
Значение электромагнитной индукции непрерывно зависит от точки.
Поэтому для заданной пары катушек значения интерполируются в точках,
находящихся на расстоянии 1 мм. В специальном поле интерфейса рисуется
диаграмма распределения величины индукции B
в соответствии со
следующей схемой.
B∈ [ 0 ,a ] - голубой
B∈ [ a ,b ] - светло-голубой
B∈ [ b ,c ] - желтый
величине
соответствие цветов
электромагнитной индукции B
B∈ [ c , d ] -оранжевый
B>d
- красный
Параметры a,b,c,d определяет пользователь. В текущей реализации
интервал [0,d] делится на 5 равных частей, где d –максимальное значение
величины индукции в вычисленных точках. Таким образом, диаграмма сразу
определяет распределение интенсивности и согласуется с положением
выбираемой точки интерполяции.
3.2.4 Алгоритмы интерполяции и визуализации
Алгоритм получения интерполированного значения в точке:
1. Выбор строки из текстового файла (выбор активных катушек)
2. Запись данных из строки в массив данных для интерполяции
3. Выбор точки
30
4. Интерполяция
5. Выбрана ли новая пара активных катушек:
Да – переход к шагу 1
Нет – переход к шагу 3
Алгоритм визуализации:
1.
Выбор строки из текстового файла (выбор активных катушек)
2. Запись данных из строки в массив данных для интерполяции
3. Интерполяция в цикле по выбору точки с шагом 1мм
4. Визуализация
Общий алгоритм:
1. Визуализация
2. Выбор точки
3. Интерполяция для выбранной точки
4. Выбрана ли новая пара активных катушек:
Да – переход к шагу 1
Нет – переход к шагу 2
3.2.5 Описание интерфейса
Рисунок 3.2.5.1. Интерфейс вкладки «Интерполяция
экспериментальных данных».
Пользователь приложения указывает полный путь до входного файла
(элемент 1 на Рисунке 3.2.5.1), в котором записаны данные, полученные
экспериментальным путем. Используются файлы в формате .txt, колонки
разделены табуляцией. После указания входного пути к файлу, пользователь
может выбрать пару катушек, для которой необходимо провести вычисления
(элемент 2 на Рисунке 3.2.5.1), и указать расстояние до центра катушки.
Указать расстояние можно двумя способами – ползунком (элемент 3 на
31
Рисунке 3.2.5.1) или введя точное расстояние в предназначенное для этого
поле (элемент 4 на Рисунке 3.2.5.1). При изменении расстояния одним из
приведенных способов в режиме реального времени будет обновляться
результат интерполяции в поле результата. (элемент 5 на Рисунке 3.2.5.1).
Также, сразу после выбора пары катушек отображается диаграмма величин
индукции электромагнитного поля (элемент 6 на Рисунке 3.2.5.1).
Расположение ползунка согласовано с диаграммой под ним.
3.2.6 Пример использования приложения
Шаг 1. Выбор входного файла с экспериментальными данными. В
данном примере (Рисунке 3.2.5.1) это «C:\\magneticFieldFiles\\inputTable.txt»
Шаг 2. Выбор пары катушек, для которых производится интерполяция.
В данном примере (Рисунке 3.2.5.1) это пара 1, ей соответствует 1 строка из
таблицы входных данных.
Шаг 2. Выбор расстояния от центра первой катушки до точки
интерполяции. В данном примере (рис. 2) это расстояние равняется 340 мм и
значение в этой точке равняется 6.915 мТ.
Таким образом, изменяя ползунком расстояние можно в режиме
реального времени наблюдать изменение результатов интерполяции.
ГЛАВА 4. Оптимизация визуализации.
В параграфе 2.2 приведены результаты визуализазии, полученные с
помощью пакета научной визуализации Paraview. На них отсутствует
изображение расположения пациента. В этой главе рассматривается метод
визуализации, который позволяет отображать положение пациента, катушки и
распределение электромагнитного поля в месте взаимодействия с телом
32
пациента. Для этого используются средства Paraview, модель тела и
специальный сценарий - скрипт на языке программирования Python. Данная
схема позволяет нам выбрать пару активных катушек и соответствующий
файл с данными вычислений или интерполяции и выполнить визуализацию.
Также мы можем посмотреть результаты визуализации для каждой из пар
катушек шаг за шагом. Данная визуализация является практическим
инструментом и может позволить физиотерапевту подобрать подходящую
схему лечения.
В пакете Paraview существует возможность схематически отобразить
тело пациента. Для этих целей мы используем шаблон из [21] (файл имеет
расширение .stl) и пакет Blender для редактирования 3D изображений [22].
Для того чтобы повторять похожие действия для различных пар активных
катушек удобно использовать сценарии (скрипты), т. е. макросы, которые
можно запускать с разными параметрами. Сценарии написаны на языке
программирования Python и имеют расширение .py.
Основными действиями при визуализации являются:
1. Загрузка файла
Для визуализации результатов вычислений или результатов
интерполяции необходимо загрузить файл: File – Open… -Select file и нажать
на кнопку Apply
2. Добавление примитивов
Важным пунктом является отображение позиции электромагнитных
катушек в соответствии с масштабом кровати. Данные элементы могут быть
добавлены как примитивы (цилиндр), путем выполнения следующих
действий:
Найти в меню пункт Sources, выбрать object, и затем в свойствах
объекта указать его параметры и нажать Apply.
3. Трансформация объекта
Для установки примитивов в корректное положение, которое в нашем
случае определено размерами кровати и позициями катушек, требуется
33
произвести трансформацию. Для этого нужно выбрать примитив в Pipeline
Browser, и затем выбрать Filters – Alphabetical – Transform. После этого в
свойствах объекта задать параметры и нажать Apply.
4. Ресэмплирование
Для добавления модели в сцену (для отображения позиции тела
пациента) используется шаблон. После задания параметров примитивных
объектов и шаблона, возможно получить финальное изображение. Это
делается следующим образом: Выбрать модель в Pipeline browser, затем
Filters – Alphabetical – Resample with Dataset. Файл с данными вычислений
нужно указать как Input, файл с моделью нужно выбрать как source. После
нажатия кнопки Apply отображается результат влияния вычисленного
магнитного поля на модель. Для выбора цветовой схемы используются
Properties (пункт Colors).
Легко заметить, что часто приходится повторять одни и те же действия
по несколько раз. Для того чтобы оптимизировать взаимодействие с Paraview,
используются сценарии, которые описывают все необходимые действия.
Для записи сценария необходимо выбрать Tools – Start Scenario-name.
Затем необходимо определить необходимые параметры, выполнить
необходимые действия и выбрать Tools – Stop Scenario-name.
Пример записанного сценария может выглядеть следующим образом:
from paraview.simple import *
paraview.simple._DisableFirstRenderCameraReset()
# create a new 'STL Reader'
body30stl = STLReader(FileNames=['C:\\magneticFieldFiles\\body30.stl'])
# get active view
renderView1 = GetActiveViewOrCreate('RenderView')
# show data in view
body30stlDisplay = Show(body30stl, renderView1)
# trace defaults for the display properties.
body30stlDisplay.ColorArrayName = [None, '']
34
# reset view to fit data
renderView1.ResetCamera()
# create a new 'Transform'
transform1 = Transform(Input=body30stl)
transform1.Transform = 'Transform'
# Properties modified on transform1.Transform
transform1.Transform.Translate = [300.0, 170.0, 0.0]
transform1.Transform.Rotate = [0.0, 0.0, -130.0]
transform1.Transform.Scale = [4.0, 4.0, 4.0]
………………………………………………………………
Таким образом, оптимизированный алгоритм визуализации может быть
следующим:
1. Вычисляются магнитные поля для выбранных пар катушек. Файлы
нумеруются в соответствии с порядком выбора данных пар.
2. Полученные файлы помещаются в рабочую директорию, туда же
помещается 3D модель тела (шаблон) и файл, содержащий скрипт. В
скрипте указан путь к данной директории.
3. Скрипт загружается в Paraview.
35
Рисунок 4.1. Пример визуализации с моделью тела для катушек C1 и В2
Рисунок 4.2. Пример визуализации с моделью тела для катушек А2 и C2
36
Заключение
В рамках выпускной квалификационной работы были получены
следующие результаты:
1) Были разработаны и реализованы алгоритмы интерполяции данных,
полученных в результате моделирования электромагнитного поля,
что позволило значительно сократить время вычислений.
2) Был разработан и реализован алгоритм интерполяции значений
элект ромагнитной индукции в заданной области по
экспериментально полученным данным.
3) Реализация данных алгоритмов была включена в ранее
разработанное приложение для моделирования электромагнитного
поля.
4) Реализован более удобный способ взаимодействия разработанного
приложения с пакетом научной визуализации Paraview.
5) Разработанное приложение может быть применено в лечебной
практике, где визуализация результатов вычислений используется
врачом для выбора оптимальных режимов в процессе лечения. В
частности результаты квалификационной работы были
использованы в лаборатории “Применение информационных
технологий в медицине” Технического Университета в Софии.
37
Литература
1. В. Александров. Электромагнитные поля и экология. СПб, 2005.
2. В. Александров. Экологическая роль электромагнетизма. СПб, изд.
СПбГПУ, 2006.
3. В. Ефремов. Метод и лечебно-диагностическая система на основе
низкочастотного магнитного поля малой амплитуды. Дисс. на соиск.
ктн, СПб, 2006.
4. Л.Н.Галль. В мире сверхслабых. Нелинейная квантовая
биоэнергетика. СПб, 2009.
5. Л.Н. Галль. Биоэнергетика – магия жизни. СПб, изд. АСТ, 2010.
6. Куд р и н Б . К . , А м п и л о в а Н . Б . А л г о р и т м п о с т р о е н и я
эл е кт р ома г н и т н о го п ол я н е с о о с н ы х кату ше к . Материалы
межвузовского конкурса- конференции "Технологии Майкрософт в
теории и практике программирования", С.-Петербург, март 2013. С.
98-99.
7. Кудрин Б.К. Компьютерное моделирование электромагнитного поля
нескольких несоосных катушек индуктивности. Материалы научной
конференции «Герценовские чтения – 2013», 15-20 апреля 2013г. СПб.: Изд. РПГУ им. А. И. Герцена, 2013. с. 246-249.
8. Кудрин Б.К. Интерполяция электромагнитного поля в устройствах
магнитотерапии. Материалы научной конференции «Герценовские
чтения – 2014». - СПб.: Изд. РПГУ им. А. И. Герцена, 2014. с. 210212.
9. Кудрин Б.К. Реализация алгоритмов моделирования и визуализации
электромагнитного поля в низкочастотных устройствах
магнитотерапии.
Материалы научной конференции «Герценовские
чтения – 2015» - СПб.: Изд. РПГУ им. А. И. Герцена, 2014.
10. Е. Ямке, Ф. Эмде, Ф. Леш. Специальные функции. Москва, изд.
Наука, 1964.
11. D.Dimitrov. Medical Systems for Influence of Electromagnetic Field on
the Human Body (in Bulgarian), Sophia, Technical University, 2008.
38
12. N. Ampilova, D. Dimitrov, B.Kudrin. Mathematical modeling of low
frequency magnetic field in systems for magnetotherapy . Proc. 8 Int.
Conf. CEMA13, 17-19 Oct. 2013, Sofia, Bulgaria. p.48-51.
13. B. Kudrin, A. Dimitrov. An algorithm for visualization of low-frequency
magnetic signals in systems for magnetotherapy. Proc. 8 Int. Conf.
CEMA13. 17-19 Oct. 2013, Sofia, Bulgaria, p.31-35.
14. B. Kudrin, A. Dimitrov. Computer visualization of low frequency
magnetic signals in systems for magnetotherapy with variable
parameters. Proc. 8 Int. Conf. CEMA13. 17-19 Oct. 2013, Sofia,
Bulgaria, p.36-39.
15. B. Kudrin, I. Soloviev. On interpolation methods of low frequency
magnetic field in systems for magnetotherapy. Proc. 9 Int. Conf.
CEMA14. 16-18 Oct. 2014, Sofia, Bulgaria, p.154-157.
16. B. Kudrin, V.Nikolov, D. Dimitrov. On the mathematical model of
interpolation of low frequency magnetic field using experimental data.
Proc. 9 Int. Conf. CEMA14. 16-18 Oct. 2014, Sofia, Bulgaria, p.80-84.
17. B. Kudrin, V. Nikolov, D. Dimitrov. Algorithms of interpolation and
visualization of low frequency magnetic field by using experimental data.
Proc. 9 Int. Conf. CEMA14. 16-18 Oct. 2014, Sofia, Bulgaria, p.43-45.
18. B.Kudrin, V. Nikolov. On a visualization of low-frequency magnetic field
by using Paraview package. Proc. 10 Int. Conf. CEMA15, 15-17 Oct.
2015, Sofia, Bulgaria. p.18-21.
19. http:\infoventures.com
20. Python programming language: http://www.docs.python.org
21. http://tf3dm.com/3d-model/boy-10833.html
22. https://www.blender.org/
23.
ParaView Tutorial:
http://paraview.org/Wiki/images/d/d5/ParaViewTutorial398.pdf
39
Приложение 1. Описание пользовательского интерфейса
вкладки Вычисление.
Пользовательский интерфейс вкладки Вычисление содержит три
логических блока:
a.
Блок, в котором задаются параметры области (размеры
области и шаг сетки по каждой координате)
b.
Блок, в котором задаются параметры для каждой катушки.
Содержит следующие группы столбцов (номер катушки,
выключатель катушки, координаты начала катушки,
координаты вектора оси катушки, радиус, силу тока, число
витков и число намоток в витке). При выключении катушки
все поля для одной катушки выключаются.
c.
Блок для вычисления вектора электромагнитной индукции в
точке (возможно, не лежащей на сетке).
При нажатии на кнопку «Рассчитать значение поля и сохранить
результаты в файл» происходит расчет электромагнитного поля для
сетки с указанными параметрами и с указанной конфигурацией
активных катушек. Результаты сохраняются в файл формата .vtk,
используемого для последующей визуализации в пакете научной
визуализации ParaView.
40
Приложение 2. Описание работы с пакетом ParaView.
1. Описание формата файла .vtk
Разработанное приложение генерирует файл, предназначенный для
обработки пакетом научной визуализации ParaView. Данный файл имеет
специализированный формат .vtk, и содержит следующие блоки данных:
1) # vtk DataFile Version 2.0 – версия vtk, применяемого для обработки
данного типа файла;
2) Заголовок;
3) ASCII – информация о кодировке;
4) DATASET STRUCTURED_POINTS – указание обработчику, что
используется равномерная сетка, с ячейками, имеющими форму
прямоугольных параллелепипедов;
5) DIMENSIONS 40 100 100
ORIGIN 0 0 0
SPACING 4 4 4
POINT_DATA 400000 блок, содержащий информацию о размерах области, начале
координат, шаге сетки и количестве точек
6) VECTORS vectors double – заголовок блока с данными
И далее в каждой строке по 5 векторов, каждый из которых задан 3
координатами. При генерации файлов векторы перечисляются в
порядке, в котором обходится трехмерный массив с данными
(сначала по 3 координате, потом по второй, потом по первой)
2. Описание работы с ParaView
Для загрузки сгенерированного файла для обработки требуется
выполнить следующие действия:
File-Open-В появившемся меню выбрать файл-ОК (Рисунок 1а и 1б):
41
Рисунок 1а.
Рисунок 1б.
После загрузки файл добавится в Pipeline Browser (Рисунок 2а):
42
Рисунок 2а.
Далее необходимо нажать кнопку Apply на вкладке Properties (Рисунок 2б):
Рисунок 2б.
И файл станет доступным для применения фильтров (в основном окне
появится параллелепипед с белыми ребрами) (Рисунок 2в):
Рисунок 2в.
43
Для того, чтобы посмотреть силовые линии исследуемого участка магнитного
поля, необходимо применить фильтр Stream Tracer, для этого необходимо
кликнуть по иконке
, далее на вкладке Properties нажать Apply.
В основном окне появятся силовые линии (Рисунок 3в):
Рисунок 3а.
Чтобы визуально оценить величину магнитного поля, можно выбрать окраску
силовых линий в соответствии с данной величиной (изменение цвета от
красного к синему означает ослабевание магнитного поля). Для этого при
активном фильтре Stream Tracer необходимо на вкладке Active Variables
Control изменить значение Solid Color на vectors, после чего цвет сразу
поменяется (Рисунок 3б и 3в) :
Рисунок 3б.
44
Рисунок 3в.
Также можно менять место, для которого требуется нарисовать силовые
линии, двигая при активном фильтре Stream Tracer белое перекрестие в
основном окне. После выбора точки для исследования необходимо снова
нажать Apply. Расположение силовых линий на рисунке после этого
изменится. (Рисунок 3г):
Рисунок 3г.
45
В
ParaView и м е е т с я т а к ж е ф и л ьт р Probe
Location, позволяющий
интерактивно определять значение векторного поля по заданной координате.
Чтобы выбрать этот фильтр, нужно выполнить:
Filters – Alphabetical - Probe Location (Рисунок 4а)
Рисунок 4а.
Затем необходимо либо
переместить белое перекрестие в нужную точку,
либо явно указать требуемые координаты на вкладке Properties. (Рисунок 4б)
Рисунок 4б.
На вкладке Information требуется навести курсор на переменную vectors, и
всплывет строка с координатами вектора электромагнитного индукции для
выбранной точки.
Рисунок 4в.
46
Подробное руководство по использованию ParaView находится по ссылке
[23].
47
Отзывы:
Авторизуйтесь, чтобы оставить отзыв