CC BY
26УДК 004.932:532.5
Обработка изображений векторных полей с применением линейной интегральной свертки
Е. В. Ворожцов
Институт теоретической и прикладной механики им. С. А. Христиановича СО РАН, 630090, Новосибирск, Россия
Обсуждаются недостатки таких известных способов получения линий тока из векторных полей, как использование функций тока и интегрирование дифференциальных уравнений, справедливых вдоль линий тока. Описывается сравнительно новый метод визуализации линий тока - метод линейной интегральной свертки, который вытекает из заданного в цифровой форме векторного поля и существенно использует алгоритмы цифровой обработки изображений. Показано, что метод реализуется в три этапа: генерация изображения типа "белый шум" во всей области; сглаживание данного изображения вдоль линий тока на основе заданного векторного поля; улучшение полученного сглаженного цифрового изображения, для чего, в частности, предлагается использовать морфологическую операцию утончения сегментов линий тока. Дается количественная оценка точности метода линейной интегральной свертки с применением нескольких критериев. Приводятся некоторые результаты вычислительных экспериментов на двумерных векторных полях, точные линии тока которых известны. Показано, что метод линейной интегральной свертки является достаточно общим и универсальным и не требует для своей реализации какой-либо априорной информации об исследуемом явлении.
Ключевые слова: научная визуализация, векторные поля, линии тока, цифровые фильтры, линейная интегральная свертка, гидродинамика.
The shortcomings of such well-known techniques for obtaining the streamlines as the use of stream functions and the integration of differential equations, which are valid along the streamlines, are discussed. A comparatively new method for visualizing the streamlines - the method of linear integral convolution - is described. The method departs from the vector field given in digital form and employs substantially the algorithms of digital image processing. It is implemented in three stages: (1) generation of the image of the "white noise" type in the entire region; (2) the smoothing of the given image along the streamlines based on the given vector field; (3) improvement of the obtained smoothed digital image. To this end, it is proposed to use, in particular, the morphological operation of the thinning of the streamline segments. Quantitative estimates are given for the accuracy of the linear integral convolution method by using several criteria. Some results of computer experiments on two-dimensional vector fields, for which the exact streamlines are known, are presented. The method of linear integral convolution is shown to be very general and universal, it does not require any a priori information about the phenomenon under study for its realization.
Key words: scientific visualization, vector fields, streamlines, digital filters, linear integral convolution, hydrodynamics.
Введение. Научные исследования разнообразных динамических процессов часто связаны с необходимостью детального анализа этих процессов на основе имеющихся векторных полей. В свою очередь, этот анализ опирается на детальное изучение линий тока, получаемых из векторных полей. При этом к алгоритмам компьютерной визуализации линий тока предъявляются требования точности, скорости и высокой детальности картины линий тока [1].
Одним из известных способов построения линий тока по заданному векторному полю является численное интегрирование обыкновенных дифференциальных уравнений вдоль линии тока:
dx/dt = u, dy/dt = v, dz/dt = w, (1)
где x,y,z - декартовы прямоугольные координаты точки в исследуемой пространственной области; t - параметр вдоль линии тока; u,v,w - проекции вектора скорости среды на оси координат Ox, Oy, Oz. Известен также способ построения линий тока как линий уровня функции тока. Оба эти способа построения линий тока реализованы в системе инженерных и научных расчетов MATLAB. Выбор данного пакета программ в настоящей работе обусловлен тем, что в его составе имеется пакет прикладных программ "Image processing toolbox", который включает более 100 функций, реализующих наиболее распространенные методы и алгоритмы обработки цифровых изображений [2]. В системе MATLAB имеется встроенная функция streamline(X,Y,Z,U,V,W, startx, starty, startz), в которой численно интегрируются либо три уравнения (1), либо, в случае линий тока на плоскости, какие-либо два уравнения системы (1). Массивы X, Y, Z определяют координаты для U, V, W; startx, starty и startz определяют стартовые положения линий тока. Основная трудность при использовании функции streamline(...) связана с выбором начальных точек линий тока. При неудачном выборе этих точек можно не заметить (пропустить) некоторые важные особенности изучаемого векторного поля, такие как области завихрения или рециркуляции, линии растекания, критические точки и т. п.
Известный способ построения линий тока как линий уровня функции тока имеет существенное ограничение: он применим только к двумерным векторным полям. Когда число пространственных переменных равно трем, функцию тока удается построить лишь в нескольких частных случаях [3].
2. Метод LIC. Рассмотрим метод LIC ("line integral convolution" - линейная интегральная свертка).
2.1. Описание метода. Метод LIC, предложенный в работе [1], вызвал большой интерес различных исследователей, нуждающихся в картинах линий тока. Достоинства LIC-метода оказались столь значительными, что в настоящее время в мире насчитываются тысячи научных работ, в которых либо LIC-метод применяется в прикладных задачах, либо предлагаются его модификации.
В экспериментальной аэродинамике широко известен метод нанесения пленки из смеси масла и краски (или сажи) на поверхность обтекаемого тела. Когда поток газа обдувает поверхность тела, смесь масла и краски преобразуется в линии тока на поверхности. Данный способ визуализации течения наиболее эффективен при наблюдении отрывов потока газа и его повторных присоединений [4]. Однако следует отметить наличие элементов случайности:
1) так как указанная смесь наносится на поверхность вручную, толщина ее пленки различна на разных участках поверхности;
2) на наклонных участках поверхности толщина пленки меняется еще и под действием силы тяжести (пленка стекает);
3) на различных участках пленки локальная концентрация сажи в масле также неодинакова.
Тем не менее метод масляно-сажевой пленки позволяет получать очень информативные картины линий тока на поверхностях обтекаемых тел (см., например, работу [5]).
Метод LIC включает следующие основные этапы:
- генерация цифрового изображения типа "белый шум" во всей пространственной области;
- сглаживание данного изображения вдоль линий тока на основе заданного поля скоростей;
- улучшение полученного сглаженного изображения: контрастирование, выравнивание гистограммы яркостей и т. п.
Разработанная головная "MATLAB''-программа, реализующая LIC-метод, имеет имя StreamLIC.m. В этой программе метод LIC реализован для случая двух пространственных переменных: LIC-изображения линий тока строятся в сечении y = const, т. е. действующими пространственными переменными являются переменные x и z.
Для улучшения качества итогового LIC-изображения сначала осуществлялась интерполяция исходного поля скоростей (u,w) на более мелкую равномерную пространственную сетку в плоскости (x, z). Это небходимо для того, чтобы не потерять информацию о мелких деталях векторного поля в областях сгущения исходной сетки [6].
Построение мелкой сетки выполнялось по формуле
hf = min(hb ИзУ/тап (2)
где hi, h3 - среднеарифметические значения шагов исходной неравномерной сетки вдоль осей x и z, соответственно; fmagn - задаваемый пользователем коэффициент увеличения детальности сетки при измельчении исходной сетки, fmagn > i. Реальное ограничение сверху на fmagn диктуется лишь быстродействием используемого компьютера. При этом чем детальнее исходная пространственная сетка, тем меньшее значение fmagn можно брать. Заметим, что машинное время, необходимое для получения LIC-изображения, пропорционально величине fm2agn.
Формула (2) показывает, что новая, более мелкая сетка состоит из квадратных ячеек, длина стороны ячейки равна hf. Таково требование исходного LIC-метода. Для того чтобы ослабить данное требование и рассматривать ячейки мелкой сетки с неравными сторонами, требуется довольно существенная переработка программы StreamLIC.m. Однако в этом нет необходимости. Следует отметить, что в работах [7, 8] LIC-метод обобщен на случай криволинейной пространственной расчетной сетки.
Величины u, w интерполировались из центров исходных ячеек в центры ячеек более мелкой квадратной сетки с помощью встроенной "MATLAB''-функции interp2(x,z,U,xf,zf,'linear'), где U - одна из величин u, w; (x,z) - координаты точек исходной сетки, в которых заданы значения U; (xf,zf) - координаты центров ячеек более мелкой сетки; при задании аргумента 'linear' реализуется билинейная интерполяция, при задании аргумента 'cubic' - бикубическая сплайн-интерполяция.
Обозначим через Wik числовое значение "белого шума", генерируемого в центре ячейки мелкой сетки. "Белый шум" { Wk} генерировался по формуле Wik = randn(N3j,Nif) (randn - "MATLAB''-функция, формирующая массив размера NfxNy, элементами которого являются случайные величины, распределенные по нормальному закону с математическим ожиданием 0 и среднеквадратичным отклонением i; Nlf, N3f - числа ячеек мелкой сетки в направлениях осей x и z соответственно). Элемент случайности, присутствующий в алгоритме LIC-метода - единственный элемент, который является общим для LIC-метода и масло-сажевого метода из экспериментальной аэродинамики. Однако метод линейной интегральной свертки имеет существенное преимущество по сравнению с масло-сажевым методом. LIC-метод позволяет определять линии тока не только на поверхностях обтекаемых тел, но и в пространстве между телами, заполненном движущейся средой.
В дальнейших вычислениях с использованием LIC-метода осуществлялось деление длины сторон мелкой сетки на величину hf, для того чтобы получить квадратные ячейки со стороной, равной единице.
Для улучшения "белого шума" в работе [i] применялась дополнительная фильтрация нижних частот с помощью фильтра "Butterworth" четвертого порядка с частотами обрезания 128, 84, 64 и 32. Подробное описание указанного фильтра приведено в [2, 9]. В разработанной программе StreamLIC.m применение фильтра "Butterworth" считалось нецелесообразным по следующим причинам:
- численная реализация этого сглаживающего фильтра требует дополнительных расходов машинного времени;
- с уменьшением частоты обрезания "белый шум" становится более размытым, что ухудшает качество LIC-изо-бражения (см. рис. 6 в работе [i]).
Локальное поведение векторного поля можно аппроксимировать, вычисляя локальную линию тока, начинающуюся в центре ячейки (i, k) мелкой сетки. Пусть P0 - стартовая точка линии тока, P0 = (xk zk). Адвекция точки P0 вдоль линии тока вычисляется по формулам
Xm = Xm-i + AXm-i, (3)
Zm = Zm-i + ^m-U m=i,2v.. .
Алгоритм вычисления приращений координат Axm-i, Azm-i при движении вдоль линии тока со скоростью v = (u, w) сводится к рассмотрению следующих восьми случаев.
Случай i. Пусть (im-i, km-i) - центр ячейки, в которой находится точка Pm-i, пусть также
u;m-i,km-i >a w.^ >a (w/u)im-ikm-i<i
т. е. вектор скорости, выходящий из центра ячейки (im-i, km-i), пересекает правую сторону ячейки (im-i, km-i), (рис. i,a). Тогда, проведя элементарное геометрическое вычисление, получаем:
AXm-i = 0,5 Azm-i =AXm-i (w / u), _ t ■
Рис. 1. Случаи расположения вектора скорости относительно сторон ячейки Случай 2. Пусть и1 к ] > 0,w¡ к > 0, но ^ / и) , к > 1 (рис. 1,6). Тогда
^т-1 = 0,5, Агт_1 = /(^/и^^. Случай 3. и( ^ ] > 0,w¡ к < 0, /и), к ] |< 1 (рис. 1,в). Тогда
Ах„_1 = °Д Агт_1 = Ах„-1(w/ и)^^ < Случай 4. и, к ] > 0,w, к < 0, |(w/и), к ] | >1 (рис. 1,г). Тогда
Агт-1 = 0, 5, АХт-1 = Агт-1 / (W / и> 0.
Случай 5. и, к < 0,w,. к >0, /и), к ] |< 1 (рис. 1,д). Тогда
АХт-1 = 0,5, А^т-1 = АХт-1 (W / «),„_,,^ > 0.
Случай 6. и, ^ < 0,w, к ] |(w/и), к ] | >1 (рис. 1,е). Тогда
Агт-1 = 0, 5, АХт-1 = Агт-1 / (W / М> 0.
Случай 7. и, ^ < 0,w, к < 0, /и), к ] |< 1 (рис. 1,ж). Тогда
АХт-1 = -0,5, А^т-1 = АХт-1 (W / «< 0.
Случай 8. и1 < 0, w, к < 0, |(w / и), к ] | >1 (рис. 1,з). Тогда
ДТт-1 =-0,5, Ахт-1 = АТт-1 / (W / и ^^ < 0.
Сравнивая случаи 1, 3, 5, 7, нетрудно заметить, что соответствующие формулы для вычисления приращений Ахт-1 и Агт-1 можно объединить в следующие две формулы:
АХт-1 = 0,5 ^ёПК-.А,-, X А2т-1 = АХт-1 (W / М • (4)
Аналогичным образом сравнивая случаи 2, 4, 6, 8, эти случаи также можно объединить в две формулы
^т-1 = 0,5 ^ёПКт-.А,-, ); АХт-1 = А2т-1 / (W / и ) ^^ - (5)
Вычислив приращения Ахт -1 и Агт -1 по одной из формул (4), (5), по формулам (3) легко найти координаты хт, 2т следующей точки Рт. Однако при этом важно обеспечить, чтобы точка Рт вошла в соседнюю ячейку сетки [1]. Для этого в соответствии с [1] в формулах (4), (5) вместо множителя 0,5 использовался множитель е1 = 0,5(1 + е)
6
а
в
г
(е- малое положительное число; в приведенных ниже примерах расчетов значение е = 10-6). Итак, формулы для вычисления Ахт _1 и Агт _1 имеют окончательный вид
- если \ ^/и), к ] \ < 1, то
Лхт_1 = е1 51§П(и,т_1,кт_1 ) А2т_1 = Ахт_1 (™ / и^^ ; (6)
- если \ (^/и) к ] \ >1, то
А2т_1 = е1 51§П(^,т_1,кт_1 ) АХт_1 = А^т_1 / (^ / и^^ • (7)
Локальная линия тока продолжается из центра ячейки и в отрицательном направлении, т. е. с компонентами (-и,-м>). Это обеспечивает симметрию линии тока относительно центра ячейки. Для того чтобы найти соответствующие приращения Ахт_1, Ах'т_х, в формулы (6), (7) подставлялось значение -и вместо и и значение -V вместо V Теперь подсчитываются соответствующие длины дуг локальной линии тока:
А? , = (Ах2 , +Аг2 ,)°'\ Ал' , =(Ах'2, + Аг'2,)°'\
Сглаживание значения яркости изображения в ячейке (,, к) мелкой сетки осуществляется по формуле ЫС-ме-тода:
I г Т^т. кт Ьт + 1 К к'т * .к = —-1-тГ-, (8)
I Ьт + 1 К
т=0 т=0
где
Г rSm-l+^Sm - 1
hm = J k(w)dw, h'm = J' " k(w)dw, (9)
s0 = 0, sm = sm-1 + Asm-1. Величина l в соотношении (8) равна такому значению i, при котором si < L < si+1, (L - длина локальной линии тока, задаваемая в терминах единичной длины стороны ячейки мелкой сетки). В работе [1] обсуждалось оптимальное значение параметра L, которое должно задаваться пользователем. С учетом рекомендаций [1] определено, что с точки зрения качества LIC-изображения оптимальными являются значения L в интервале 4 < L < 10. Следует отметить, что при расчете выходного сигнала Sik по формуле (8) суммарный расход времени пропорционален L2.
Функция k(w) в (9) - это ядро свертки. В работе [1] приводились аргументы в пользу применения функции k(w) вида
k (w) = 1 + cos(cw) 1 + cos(dw +e) (10)
2 2
Здесь первый множитель (1/2)(1 + cos(cw)) соответствует ядру, применяемому в фильтре Hanning, обсуждение которого приведено, например, в работе [9]. Множитель (1/2)(1 + cos(dw + Р)) - это функция ряби (ripple function); c и d - константы наращивания (dilation), р - сдвиг по фазе функции ряби, задаваемый в радианах.
Вследствие простоты ядра k(w) (10), интегралы в (9) легко вычисляются в аналитическом виде
Jk(w)dw = 4[b - а + sin(bc)-sin(ac) + sin(bd+e)-sin(«d+e)] +
1 sin(b(c-d)-e)-sin(a(c-d)-e sin(b(c+d)+e)-sin(a(c+d)+e)
+4[ 2(c-d) + 2jc+d) ].
(11)
Из (П) следует, что необходимо задавать различные значения с и й, чтобы избежать обращения знаменателя в нуль.
Множители (i/2)(i + cos(cw)) и (i/2)(i + cos(dw + Р)) имеют периоды Ti = 2п/с и T2 = 2n/d соответственно, т. е. периоды зависят от констант c,d. Если в функции ряби задавать низкую частоту (т. е. малые значения d), то с изменением сдвига по фазе в эта функция будет существенно меняться. Это проявляется в LIC-изображении как периодическое размывание и укручение (от англ. steepening - увеличение крутизны профиля) профилей яркости. Использование в функции ряби более высокой частоты приводит к оконному фильтру с почти постоянной выходной частотой, так как общая форма фильтра меняется незначительно. Но тогда отрезки линий тока будут более короткими. Если частота ряби превосходит предел Найквиста [9], определяемый размером пикселя, то длины отрезков линий тока становятся равными нулю. Как указывалось в [i], за счет увеличения пространственного разрешения (увеличения знаменателя fmagn в формуле (2)) всегда можно добиться и оптимальной выходной частоты, и увеличения длин сегментов линий тока, но это, как отмечалось выше, приводит к увеличению объема вычислений, необходимых для получения LIC-изображения.
Очевидно, что оптимальный подбор констант с, d, р для каждого исследуемого течения дает дополнительные возможности для улучшения качества линий тока, генерируемых по LIC-методу.
Полученное по формуле (8) LIC-изображение S = {Sik} может непосредственно использоваться для визуализации линий тока. Однако имеются возможности для дальнейшего улучшения качества LIC-изображения, предложенные в работах [6-8]. В программе StreamLIC.m реализованы два способа улучшения LIC-изображения, предложенные в [6].
Первый способ состоит в том, чтобы использовать LIC-изображение S в качестве исходного для следующей итерации по LIC-методу. Иными словами, на второй итерации вместо исходного "белого шума" W в формуле (8) используется LIC-изображение S. В результате получаются линии тока, которые различаются в большей степени, чем при выполнении лишь одной итерации по LIC-методу. Как указывалось в [6], дальнейшие итерации (третья и т. д.) не приводят к значительным различиям в результатах по сравнению со случаем, когда LIC-алгоритм выполняется дважды. Кроме того, суммарное потребляемое машинное время прямо пропорционально числу итераций по LIC-методу. Поэтому рекомендуется выполнять не более двух итераций.
Другой способ улучшения LIC-изображения, предложенный в работе [6], заключается в приближенном выравнивании гистограммы яркости исходного "белого шума". Фильтр выравнивания гистограммы яркости сдвигает значения Wik к крайним значениям диапазона. Для приближенного выравнивания гистограммы в [6] предлагалось перевычислять значения Wik по формуле
W, Wkk > 0,
W . » (12)
1- W r, W k < о.
В (12) величина r задается пользователем; в работе [6] для r рекомендован диапазон 5 < r < 7.
Так как вычисления по формуле (8) - это операция сглаживания исходного LIC-изображения, то в результате происходит потеря контрастности получаемого изображения. В работе [10] предлагалось применять контрастирование LIC-изображения по формуле
4
S'k = SP, (13)
где Sik - изображение, получаемое отображением LIC-изображения Sik на интервал [0,1]. Заметим, что если 0 < S k < 1, то показатель степени в (13) варьируется от 1/8 до 4.
Теперь приведем некоторые результаты численных экспериментов с использованием LIC-метода для течения среды с круглым вихрем, функция тока которого имеет вид
y (x,z) = (1/2) [(x -X) + (z - Zc)2] (14)
(Xc, Zc) - координаты центра вихря). Тогда с помощью формул u = dy/dz, w = -dy/dx находим: u(x,z) = z — Zc, w(x,z) = -(x - Xc). Вектор вихря скорости rot v = (0, 2, 0) имеет лишь одну ненулевую компоненту, так что |rot v| = 2.
а б в
Рис. 2. Линии тока круглого вихря: а - линии тока, полученные как линии уровня функции тока y(x,z); б - LIC-изображение lW, цвет C = lW; в - LIC-изображение lW, цвет C вычислен по формуле (15) при у = 0,5
На рис. 2,а показаны линии тока рассматриваемого течения, полученные как линии уровня функции тока (14) с помощью MATLAB-функции contour.m. Видно, что линии тока - концентрические окружности с центром в точке (XC,ZC).
Следуя работе [9], обозначим через lW результат однократного применения метода LIC к "белому шуму" W, через llW = IS - результат двукратного применения метода LIC. В примерах расчетов по LIC-методу, приведенных на рис. 2,б, в, в соотношении (10) использовались значения d = 0,1, ß= 0,15. Для того чтобы выполнялось неравенство c Ф d, было выбрано значение c = 0,05. Это позволяло получить LIC-изображения, сопоставимые по качеству (прежде всего, по длинам сегментов линий тока, ухватываемым LIC-методом) с приведенными в работе [1]. На рис. 2,б показано цветное LIC-изображение lW, полученное при вызове MATLAB-функции mesh(xx,zz, lW), где окраска пропорциональна высоте поверхности lW. Если в список параметров функции mesh ввести четвертый параметр C, управляющий локальным цветом, то можно усилить цветовое различие между соседними линиями тока на LIC-изображении. С помощью нескольких численных экспериментов эмпирически было найдено, что задание цвета по формуле
Ca = {SY, S* > 0 (15)
л [-1 Sjt Г, Sjt < 0
с показателем у = 1/а приводит к существенному усилению цветового контраста LIC-изображения S = lW (ср. рис. 2,б и 2,в). Автором настоящей работы проводились вычислительные эксперименты, в которых не применялось приближенное выравнивание гистограммы по формуле (12), тем не менее выполнялась операция вычисления цвета согласно (15). В таких случаях целесообразность операции (15) вытекает из того, что положительную постоянную у можно подобрать эмпирически из требования повышения контраста итогового LIC-изображения.
В случаях, показанных на рис. 2,б, в, в (8) использовалось значение L = 10 длины окна осреднения как в положительном, так и в отрицательном направлениях вектора скорости. В формуле (2) при переходе к более мелкой сетке коэффициент увеличения fmagn = 5. Это означает, что на исходной сетке окно сглаживания имело сравнительно небольшой размер примерно 4*4 ячейки.
На рис. 3,а показано LIC-изображение, полученное в результате двукратного применения LIC-метода. Сравнивая рис. 3,а и рис. 2,в, видим, что повторное применение LIC приводит к улучшению связности сегментов линий тока; иными словами, средняя длина указанных сегментов увеличилась.
Пусть, как и в работе [6], e - оператор приближенного выравнивания гистограммы яркости по формуле (12). На рис. 3,б показано LIC-изображение, полученное в результате предварительного выравнивания гистограммы "белого шума" и последующего двукратного применения метода LIC.
Поскольку во многих профильных журналах цветные рисунки не публикуются, целесообразно рассмотреть возможность инвертирования цветного LIC-изображения в черно-белую форму. Данная операция осуществлялась в два этапа. На первом этапе LIC-изображение S отображалось на интервал яркостей [0, 1] по формуле S' = aS + b, где
а = 1,0/(Smax — Smin), b = — aSmirb Smax = max S¡k
Smm = min Sik.
Рис. 3. Линии тока круглого вихря, полученные ЫС-методом: а - ЫС-изображение 1Ж, цвет С вычислен по (!5) при у = 0,5; б - ЫС-изображение 11еШ при а = 5, у = 0,2
На втором этапе черно-белое ЫС-изображение генерировалось путем вызова МАТЬАВ-функции B=im2bw( * ', "ЬЬ), где "ЬЬ - заданное значение порога, 0 < "ЬЬ < 1; в - получаемое бинарное изображение. Чем выше значение порога, тем больше черных пикселей в черно-белом Ь1С-изобра-жении, поскольку функция im2bw.ni полагает значение яркости изображения в ячейке (,, к) равным нулю (черный цвет), если Б'к < "ЬЬ. Эта операция называется бинаризацией по порогу. В разработанных автором МАТЬАВ-программах ЫСТеБ^т и Б"геатЬ1С.т наряду с генерацией цветного ЫС-изображения осуществляется также бинаризация. На рис. 4 приведены результаты бинаризации двух различных ЫС-изображений с различными значениями порога "ЬЬ.
На бинарных изображениях линии тока показаны только двумя цветами (черным и белым) вместо большого количества цветов в цветных изображениях, приведенных на рис. 2, 3. На рис. 4 видно, что толщины линий тока переменные вдоль линий тока. Этого можно избежать, применяя морфологическую операцию утончения объектов на бинарном изображении. При этом сегментами линий тока будем считать те связные области, значения пикселей которых равны единице (белый цвет), а области изображения, значения пикселей которых равны нулю, будем
Рис. 4. Черно-белые картины линий тока, полученные бинаризацией LIC-изображений l2 W: th = 0,53 (а) и l2eW, th = 0,56 (б)
Рис. 5. Применение морфологической операции утончения: а - изображение bl2W; б - изображение itbl2W
б
а
б
а
в
г
считать фоном. Операция В = ЬмтогрЬ(В,'-ЬЫп') выполняет эрозию объекта с учетом ряда условий для сохранения 8-связности участков остова [2, 11]. В результате сегменты линий тока, имевшие переменную толщину, превращаются в сегменты линий толщиной в один пиксель.
Обозначим через Ь оператор бинаризации изображения с помощью функции im2bw.m, через t - оператор морфологического утончения объектов в бинарном изображении. На изображении линии тока показаны белым цветом на черном фоне. Так как обычно линии тока изображаются черным цветом на белом фоне, можно осуществить операцию инвертирования изображения, при которой белые пиксели становятся черными, а черные - белыми. Обозначим этот оператор инвертирования символом /', т. е. = 1 - М^. На рис. 5,а показано исходное бинарное изображение Ьl2W, /^„=20 в (2), размер изображения 390x780 пикселей. На рис. 5,6 видно, что подавляющее большинство сегментов линий тока теперь имеет толщину в один пиксель.
На рис. 6,а показаны линии тока, полученные с помощью МЛТЬЛБ-функции Б-ЬгеатМпе(...) для задачи о безвихревом обтекании сферы потоком идеальной несжимаемой жидкости. В сечении у = 0, проходящем через центр сферы, значения составляющих скорости ы,м> задавались в центрах ячеек равномерной сетки в соответствии с точным решением [1]. Результат применения ЫС-метода, приведенный на рис. 6,6, получен при следующих константах в (10): ё = 0,05, с = 2ё, в= 3ё. В (8) использовалось значение Ь = 10 длины окна осреднения как в положительном, так и в отрицательном направлениях вектора скорости.
2.2. Точность Ь1С-метода. В работе [10] осуществлена сравнительная количественная оценка шести методов визуализации линий тока включая метод ЫС на основе векторных полей. Было сгенерировано 500 векторных полей с различным количеством критических точек различных типов (например, точек типа седла, притягивающего узла и т. д.). Реализованы следующие количественные критерии оценки точности различных алгоритмов построения линий тока.
1. Определение количества критических точек в векторном поле и вычисление на этой основе процентной доли количества векторных полей (из 500 полей), в которых конкретный метод дал неправильное количество критических точек. По данному критерию ЫС-метод оказался наиболее точным из шести рассмотренных методов.
2. Определение типов выявленных критических точек. Здесь ЫС-метод по точности оказался на третьем месте среди шести методов.
3. Определение средней ошибки локализации критической точки на основе изображения линий тока, построенного одним из шести методов. Положение критической точки определялось пользователем визуально на экране монитора компьютера и сразу заносилось компьютером в память, где хранились и точные координаты соответствующих критических точек. После этого вычислялась абсолютная ошибка в определении положения критической точки. В результате установлено, что метод ЫС обеспечивает наиболее высокую точность определения положения критических точек по сравнению с остальными пятью методами.
4. Определение средней абсолютной угловой ошибки в направлении касательной к линии тока вблизи места ее вхождения в критическую точку. В этом случае метод ЫС уступил по точности остальным пяти методам: указанная ошибка для ЫС превышала 10о, в то время как для наилучшего по точности метода она лежала в интервале 0,5-1,0 о.
В дополнение к описанным выше количественным критериям оценки методов визуализации линий тока на примере круглого вихря было проведено исследование влияния измельчения сетки, т. е. коэффициента/„¡¡¡„ в (2), на точность определения локальных углов наклона касательных к линиям тока в ЫС-изображении tЬl2W. Пусть
а б
Рис. 6. Линии тока в задаче об обтекании сферы: а - линии тока, полученные с помощью функции 81геат1те(...); б - ЫС-изображение 12Ш, цвет С вычислен по (15) при у = 0,5
Ш W = {Ьк1}, к = 1,., Ы3у, 1 = 1,., Ыу. В пределах окна изображения размером 3x3 с центром в пикселе Ьк1 сегмент линии тока аппроксимируется прямолинейным отрезком, описываемым уравнением г = а0 + а\х. Тогда угол ак1 между этим отрезком и положительным направлением оси х находится из уравнения 1£ ак1 = (ёг/ёх)^ = аь В соответствии с методом наименьших квадратов составляется функция
N
Р(а0, а1) = Е Рт (а0 + а1Хт - )2,
т=1
где (х„,гт) - координаты всех белых пикселей в пределах рассматриваемого окна изображения; рт - заданные положительные весовые параметры (в описываемых ниже расчетах использовались значения рт = 1 V т). Из условий экстремума дР/да0 = 0, дР/дах = 0 легко найти выражения для а0 и аь в частности а! = (АА - ^АУС^сА - ^12), где
N NN N N
А0 = Е Рт , = Е РтХт , А2 = Е Рт2т , А3 = Е РтХт2т, А4 = Е РтХ1.
т=1 т=1 т=1 т=1 т=1
Вычисления по методу наименьших квадратов выполнялись только в тех окнах размером 3x3, в которых количество N белых пикселей удовлетворяло неравенствам 2 < N < 6. Пусть а™ - точный угол наклона касательной к линии тока; в случае задачи о круглом вихре
акХ = агс1^[-(хк,1--Хс)/( г^ - 2с) ]. С целью уменьшения величины |ак1- - а™ | выполнялись также следующие операции: если а < 0, то а = а + 2п. Далее,
если ак1< 1,6л(4,7< аке* <6,29), то аке* =2п- аке* , если 4,7 <ак1<6,29л ак* <1,6, то ак1 = 2п - ак1, если (2,4< ак>1<2,7)л(4,7<акех<6,29), то аке* = аке* - п, если (2,4< ак* <2,7)л(4,7< ак>1- < 6,29), то ам = а*,; - п. После этой коррекции углов акл и а к* среднеквадратичная угловая ошибка Да находилась по формуле
да =
(1/N. )Е (а,-а;*.)
(16)
где N. - количество точек линий тока в ЫС-изображении tЬl2W, в которых вычислялись величины ак1 и а™ . На
рис. 7 показана поверхность Да = Да(х, г) (углы измеряются в радианах) для случая изображения tЬl2W, инвертированный вариант которого показан на рис. 5,б. Видно, что локальные погрешности Да выше в следующих секторах углов а: -п/4 < а < п/4, 3п / 4 < а < 5п / 4. На рис. 8 синим и красным цветом показаны линии Да = Да(/та„
для случаев применения соответственно билинейной интерполяции и бикубической сплайн-интерполяции при переходе от исходной сетки к более мелкой квадратной сетке. Поскольку "белый шум" W меняется при каждом запуске МЛТЬЛБ-программы генерации ЫС-изобра-жения, на рис. 8 показаны результаты первого, второго и третьего запусков программы при фиксированном значении /та„
На рис. 8 видно, что при использовании бикубической сплайн-интерполяции абсолютный максимум погрешности (16) ниже, чем в случае билинейной интерполяции. Отсюда следует предпочтительность использования бикубической сплайн-интерполяции при формировании
Рис. 7. Поверхность Да = Да(х, г)
0,5
LIC-изображений. Конечно, этот метод интерполяции требует больших затрат машинного времени, чем билинейная интерполяция, но это время меньше времени, необходимого для вычислений по формулам LIC-метода.
Необходимо также отметить, что при увеличении fmagn погрешность Да не уменьшается.
Важно отметить следующее достоинство LIC-метода: поскольку вычисления по формуле (8) проводятся в каждой ячейке независимо от других, то метод LIC весьма удобен для распараллеливания (см. также [1]).
Заключение. Метод линейной интегральной свертки является универсальным методом построения линий тока по заданным локальным векторам скорости. Реализация метода проста. В отличие от известного метода интегрирования обыкновенных дифференциальных уравнений вдоль линии тока LIC-метод не требует задания начальных координат точек линий тока, иными словами, не требует никакой априорной информации об исследуемом течении. Поскольку в методе LIC сканируются все ячейки пространственной сетки, с его помощью можно успешно выявлять разнообразные локальные и нелокальные нетривиальные особенности векторного поля, некоторые из которых перечислены во введении. Вместе с тем имеется проблема сходимости метода при измельчении сетки.
Список литературы
1. Cabral В., Leedom L. Imaging vector fields using line integral convolution // Proc. of the SIGGRAPH'93 (Special Interest Group on Computer
Graphics and Interactive Techniques), Anaheim (USA). 1-6 Aug. 1993. N.Y.: ACM Press, 1993, P. 263-270.
2. Рудаков П. И. Обработка сигналов и изображений. MATLAB 5.x. / П. И. Рудаков, В. И. Сафонов. М.: Диалог-МИФИ, 2000.
3. Лойцянский Л. Г. Механика жидкости и газа. М.: Дрофа, 2003.
4. Merzkirch W. Flow Visualization. N.Y.: Acad. Press, 1987.
5. Желтоводов А. А., Максимов А. И., Шевченко А. М. Топология пространственного отрыва в условиях симметричного взаимодейст-
вия пересекающихся скачков уплотнения и волн разрежения с турбулентным пограничным слоем // Теплофизика и аэромеханика. 1998. Т. 5, № 3. С. 319-340.
6. Okada A., Lane D. Enhanced line integral convolution with flow feature detection // Techn. Rep. NAS-96-007. NAS, June 1996.
7. Forssell L. K. Visualizing flow over curvilinear grid surfaces using line integral convolution // Proc. of the Visualization'94. IEEE Computer
Society, 17-21 Oct, 1994. Washington: D.C., 1994. P. 240-247.
8. Forssell L. K., Cohen S. D. Using line integral convolution for flow visualization: Curvilinear grids, variable-speed animation, and unsteady
flows // IEEE Transactions on Visualization and Computer Graphics. 1995. V. 1, N 2. P. 133-141.
9. Strampp W. Mathematische Methoden Der Signalverarbeitung // Strampp W., Vorozhtsov E. V. München, Wien: Oldenbourg Verlag, 2004.
10. Laidlaw D. H., Kirby R. M., Davidson J. S., et al. Quantitative comparative evaluation of 2D vector field visualization methods // Proc. of the IEEE Visualization 2001, Oct. 24-26, 2001. IEEE Computer Society / San-Diego: Electronic Edition, 2001.
11. Гонсалес Р. Цифровая обработка изображений // Р. Гонсалес, Р. Вудс. М.: Техносфера, 2005.
Ворожцов Евгений Васильевич - д-р физ.-мат. наук, проф., вед. науч. сотр. Ин-та теорет. и прикл. механики им. С. А. Христиановича СО РАН;
тел. (383) 330-38-04, e-mail: vorozh@itam.nsc.ru
018
0 ;ч-1-1-1---1---а---^---1-1-
г 4 е б 10 12 14 16 15 20
/гладл
Рис. 8. Ошибка Да = Да(/та^п): синий цвет - билинейная интерполяция, красный -бикубическая сплайн-интерполяция, сплошные линии -
первый запуск МАТЬАВ-программы, штриховые -второй запуск, штрихпунктирные линии - третий запуск
Дата поступления - 09.10.2009 г.
