CC BY
57
2015 ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА Сер. 10 Вып. 2
ИНФОРМАТИКА
УДК 004.942 А. М. Батюков
АНАЛИЗ ЦИФРОВЫХ ИЗОБРАЖЕНИЙ, ОСНОВАННЫЙ НА ПОСТРОЕНИИ СТАЦИОНАРНОГО ПОТОКА НА ГРАФЕ*)
Санкт-Петербургский государственный университет, Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7/9
В работе рассмотрена модификация алгоритма построения стационарного потока на графе, применяемого для классификации изображений. Указан способ адаптации алгоритма к использованию на многоядерных вычислительных системах путем разбиения исходного изображения на смежные области. Дана сравнительная оценка сложности базового и модифицированного алгоритмов. Исследованы вопросы зависимости от числа областей скорости вычислений и точности классификации изображений. Приведены результаты численных экспериментов. Библиогр. 7 назв. Табл. 2.
Ключевые слова: анализ изображений, стационарный поток на графе, максимизация взвешенной энтропии, параллельное программирование.
A. M. Batyukov
DIGITAL IMAGES PROCESSING BASED ON GRAPH STATIONARY FLOW CONSTRUCTION
St. Petersburg State University, 7/9, Universitetskaya embankment, St. Petersburg, 199034, Russian Federation
A modification of the algorithm to construct a stationary flow on a graph, which is applied to classification of digital images, is considered. An adaptation of the algorithm for multi core systems by the image partitioning is discussed. The dependence of calculation rate and classification accuracy on the number of partition elements is investigated. The numerical results are given. Bibliogr. 7. Table 2.
Keywords: image analysis, stationary flow on a graph, maximization of weighted entropy, parallel programming.
Введение. Одной из задач, решением которых занимаются при анализе цифровых изображений, является получение характеристик, которые в дальнейшем могут быть использованы для классификации. К таким относятся текстурные статистические признаки Харалика, фрактальные размерности, мультифрактальные спектры, некоторые морфологические признаки. Применимость и эффективность определенных методов анализа к тем или иным классам изображений могут быть определены экспериментально.
Батюков Александр Михайлович — аспирант; e-mail: batsun@gmail.com Batyukov Alexander Michailovich — post-graduent student; e-mail: batsun@gmail.com
*) Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (грант № 13-01-00782).
В настоящей работе будем рассматривать изображения, связанные с процессами распространения вещества. Примерами могут служить изображения разнообразных диффузионных процессов, изображения распространения электрического разряда в среде. Часто подобные изображения свойственны для медицины: например, подобные процессы наблюдаются на снимках различных тканей при воздействии на них лекарственными препаратами в малых дозах, а также при распространении в тканях разного рода инородных клеток, появляющихся в результате заболеваний. Для таких изображений оказывается удобным описывать процесс как движение по ориентированному графу. При этом каждому пикселю изображения ставится в соответствие вершина графа. Каждой вершине приписывается вес (мера), равная интенсивности соответствующего пикселя. Каждая дуга графа выражает отношение соседства двух вершин, при этом для каждых двух соседних вершин существуют две дуги, ведущие из одной в другую. Количество соседей обусловливается выбором окрестности, как правило, изучают окрестность из четырех ближайших вершин. Каждой дуге приписывается вес, определяемый как мера являющейся его началом вершины, деленная на число выходящих из нее ребер. Дуги с приписанными им весами можно рассматривать как поток. Данный поток нормируется так, чтобы сумма весов на дугах была равна 1. Можно считать, что полученный таким образом граф описывает начальное состояние некоторого процесса с помощью марковской цепи на нем.
Интерес представляют стационарные состояния этого процесса, а именно такие, когда построенная цепь является стационарной. Стационарной называют марковскую цепь, для каждой вершины которой сумма весов на входящих в нее ребрах равна сумме весов на исходящих. Построение стационарных процессов на графе - хорошо изученная область математики. Она очень близка к так называемой транспортной задаче линейного программирования с линейными ограничениями.
Постановка задачи. Опишем изображение как решетку пикселей ). Каждой ее точке сопоставляется значение интенсивности пикселя 0 < I(г,з) < 256. Рассмотрим окрестность ближайших соседей пикселя, состоящую из 4 пикселей. По изображению строим ориентированный граф, каждая вершина которого отвечает пикселю решетки. Всем вершинам приписываем вес, равный интенсивности пикселя. Вершины графа соединены двумя разнонаправленными дугами, если соответствующие этим вершинам пиксели - соседние. Всем выходящим из вершины дугам приписываем интенсивность узла, деленную на число соседей. Значения в вершинах и дугах графа нормируем таким образом, чтобы сумма мер вершин была равна 1. Обозначим полученный поток на графе рц. По построению мера каждой вершины равна сумме мер выходящих дуг. Нас интересует такое распределение пц значений на дугах, при котором построенная цепь будет стационарной, а это означает, что для каждой вершины графа сумма весов по всем входящим в нее дугам равна сумме весов по всем выходящим.
Взвешенная энтропия определяется в терминах начального распределения весов
как
(1)
Рассмотрим экстремальную задачу поиска максимума функции
при условиях
"^Щк Ukj = 0,
i j
J2uij = 1
Примем, что начальный поток pj - приближенное решение этой задачи (несмотря на то, что он не удовлетворяет ограничениям). Точное решение данной задачи даст построенный на графе стационарный поток. Величина взвешенной энтропии определяется начальным состоянием и полученным стационарным [1].
Решение такой задачи существует, если на графе есть замкнутые циклы. Интересна история возникновения метода решения. В 1946 г. ленинградский архитектор Г. В. Шелейховский решал задачу прогнозирования транспортных потоков в городе -транспортную задачу с линейными ограничениями. Он предложил [2] метод балансирования, заключающийся в последовательном преобразовании строк и столбцов матрицы смежности графа транспортных потоков умножением и делением на подходящие множители. В терминах весов на дугах графа это означает домножение на данные множители весов на исходящих и входящих в некоторую вершину дугах. И. В. Романовский обнаружил, что в случае, когда предложенный метод сходится, его пределом является сбалансированный план транспортных потоков, на котором достигается максимум некоторой специальной нелинейной целевой функции, называемой взвешенной энтропией. Чуть позже Л. М. Брэгман нашел доказательство сходимости метода балансировки [3] и обобщил [4] метод Шелейховского для решения произвольных задач линейного программирования. В настоящее время для транспортных задач принято называть такой метод именем Шелейховского, а во всех остальных случаях -методом Брэгмана.
Описанный метод балансировки был применен для анализа цифровых изображений путем построения стационарного распределения на ориентированном графе, связанном с изображением, в работах [5, 6]. Было показано [6], что для изображений растворов серебра различной концентрации взвешенная энтропия может рассматриваться как классификационный признак. Однако авторы отмечают, что время работы метода на изображениях большого размера (или высокого разрешения) велико.
В настоящей работе описаны две реализации алгоритма метода балансировки -базовый (предложенный в указанных выше статьях) и модифицированный, основанный на разбиении изображения на части и параллельного построения стационарного потока на каждой из частей. Приведены оценки скорости работы и результаты экспериментов для двух классов изображений патологий печени.
Метод решения — базовый алгоритм. Алгоритм построения стационарного потока на графе основан на методе Брэгмана. Как было показано в [2], сходимость метода не зависит от порядка перебора столбцов и строк матрицы смежности. В работе [5] Н. Б. Ампиловой и Е. И. Петренко для организации перебора вершин было предложено использовать приоритетную очередь, в которой приоритетом является дисбаланс вершины q(n) = \j0ut(n) — jin(n)|, где j0ut(i) - сумма мер исходящих из вершины дуг, а jin(i) - сумма мер входящих. На каждом шаге для вершины n с максимальным приоритетом (она находится в начале приоритетной очереди) вычисляется коэффициент
(4)
Меры всех дуг, исходящих из п, делятся на 7, меры всех дуг, входящих в п, умножаются на 7. После чего производятся нормировка потока и пересчет значений в затронутых вершинах. Вершины с изменившимся приоритетом при этом меняют свое место в очереди. Критерием остановки алгоритма является то, что дисбаланс вершины в начале очереди становится меньше некоторого заданного е. При условии е = 0 мы получили бы стационарный поток, но вычисления ведутся с некоторой точностью е > 0, поэтому в результате имеем некоторое е-приближенное решение, по которому вычисляют приближенное значение взвешенной энтропии.
Указанный алгоритм был применен для расчета взвешенной энтропии изображений растворов серебра различной концентрации [6]. Было показано, что использование значения взвешенной энтропии в качестве классификационного признака достаточно эффективно.
Метод решения — модифицированный алгоритм. Основная идея модификации алгоритма состоит в том, чтобы разделить исходные данные на некоторое число подмножеств, вычисления для которых можно осуществлять независимо друг от друга. Изображение рассматривается не как единая решетка пикселей, а как набор из таких решеток. Простейший случай разбиения множества на четыре решетки - поделить решетку на две равные части по каждому из измерений. Обратим внимание, что пиксели, находящиеся на границе двух решеток (в дальнейшем будем называть их граничными точками), входят в обе из них.
На каждой из таких решеток строится свой собственный граф, на котором все дальнейшие вычисления производятся независимо от проводимых на других графах. При этом на каждом графе при построении потока сумму весов во всех вершинах нормируют не 1, а значением, равным отношению числа вершин в данном графе к общему числу вершин в графах всех частей. Полученные е-стационарные распределения на графах объединяются в одно общее: в качестве значения дугам, входящим и выходящим из граничных точек и принадлежащим нескольким частям, присваивается вес, равный сумме весов соответственно входящих и исходящих ребер каждой из частей. Заметим, что сумма весов во всех вершинах объединенного графа равна 1.
Утверждение 1. В результате вычислений по модифицированному алгоритму из е-приближенных к стационарному распределений получается не более чем 4е-приближенное к стационарному.
Доказательство. Действительно, для неграничных точек каждого из объединяемых множеств неравенство д(п) < е выполнено (это критерий остановки алгоритма балансировки на соответствующем подмножестве). Каждая граничная точка может принадлежать не более чем 4 частям, потому для нее выполнено соотношение
Как видно, распределение осталось е-стационарным во всех точках, кроме граничных. В граничных точках распределение не более чем 4е-стационарно. I
Оценим временную сложность базового алгоритма. Пусть размеры исходного изображения М на N пикселей. Тогда в графе ММ вершин. Для построения приоритетной очереди необходимо вычислить для каждой вершины ее дисбаланс, после чего среди всех вершин найти вершину с максимальным приоритетом, число требуемых операций О(ММ). На каждом шаге алгоритма для такой вершины подсчитывается коэффициент 7, после чего у некоторого числа дуг (и соответственно вершин) в графе
4
(5)
1
изменяется вес. Пересчет новых дисбалансов вершин и поиск вершины с максимальным дисбалансом требуют O(MN) операций.
Модифицированный алгоритм представляет собой K запусков базовых алгоритмов на графах с числом вершин Соответственно временная сложность вычислений на каждом шаге составляет О Каждый шаг модифицированного алгоритма требует меньше операций, чем шаг базового, в K раз.
Число шагов алгоритма, которое потребуется для построения е-стационарного распределения с нужной точностью, зависит напрямую не от M и N, а от характера изначального распределения. При этом можно ожидать, что общее число шагов, требуемых для построения K е-стационарных распределений на всех подграфах, не превосходит числа шагов, требуемых для построения е-стационарного распределения на полном графе. Это означает, что общая временная сложность вычислений по модифицированному алгоритму меньше временной сложности вычислений по базовому.
Можно заметить, что в случае использования такой модификации алгоритма все вычисления, кроме связанных с разбиением графа на подграфы и заключительного объединения результатов, - независимы, а потому могут проводиться параллельно друг с другом на многоядерном вычислителе. Это дает еще большие возможности по уменьшению времени вычислений.
Результаты численных экспериментов. При выборе алгоритма, который будет применен в конкретном случае классификации, основными факторами являются точность классификации и скорость работы алгоритма. Для оценки этих характеристик классифицирующего алгоритма его запускают на предварительно разделенных на классы по необходимым признакам изображениях и получают численную характеристику каждого из них. Алгоритм считают эффективно классифицирующим, когда полученные численные характеристики кластеризуемы на числовой прямой и разные кластеры соответствуют различным классам изображений. Может быть рассчитан некий показатель, который можно назвать процентом точности (accuracy) алгоритма. В качестве этого показателя, как правило [7], используется отношение числа верно отнесенных в соответствии с численной характеристикой к тому или иному кластеру точек к их общему количеству. Кроме того, запуск различных алгоритмов на одном и том же наборе исходных изображений позволяет получить усредненные данные о скорости работы разных алгоритмов.
Для проверки утверждения был проведен численный эксперимент по вычислению взвешенной энтропии для изображений патологий тканей печени. Изображения имеют размер 2584 пикселя в ширину, 1936 пикселей в высоту. Они предварительно были разделены на два класса согласно изображенным на них заболеваниям печени. Первый класс составляли изображения тканей печени с жировой дистрофией, второй - с выраженным полнокровием. Первый класс изображений состоял из 11 представителей, второй - из 14.
Для каждого из изображений взвешенная энтропия была определена двумя способами - по базовому и модифицированному алгоритмам с точностью е = 0.0001. Для модифицированного алгоритма были рассмотрены его версии с K = 4 и K = 16, разбиение на части производилось путем деления каждой из сторон изображения на 2 или 4 равных части. После вычисления взвешенной энтропии для каждого эксперимента полученные числа были кластеризованы на числовой прямой на два кластера. Кластеризующим алгоритмом служила тривиальная вариация k-means, в качестве исходных центров кластеров брался максимум и минимум значений точек на числовой прямой. Попадание результата вычислений в первый кластер считалось показателем
того, что алгоритм классифицировал изображение как изображение тканей печени с жировой дистрофией. Попадание изображения во второй кластер считалось показателем того, что алгоритм классифицировал изображение как изображение тканей печени с выраженным полнокровием.
Эксперименты проводились на многоядерном вычислителе с числом ядер 4. Базовый алгоритм не допускает распараллеливания, поэтому скорость его работы не зависит от числа ядер системы. Как было показано выше, модифицированный алгоритм допускает распараллеливание. Поэтому для него вычисления производились трижды. В первом эксперименте вычислениям было отдано 1 ядро процессора, во втором -2 ядра, в третьем - все 4 имеющихся ядра.
Параллелизм в программе обеспечивался на уровне потоков, работающих в составе одного процесса. Вначале в рамках основного потока выполнения программы изображение разбивалось на необходимое число частей. После этого для каждой части создавался собственный вспомогательный поток (pthread). В каждом из вновь созданных потоков по соответствующей ему части изображения строился свой собственный ориентированный граф, и на нем запускался независимо от других частей алгоритм балансировки. По окончании работы всех вспомогательных потоков в рамках основного потока графы объединялись в общий граф. Логика работы основного потока была запрограммирована на скриптовом языке Tcl, логика работы вспомогательных потоков (основные вычисления) - на языке C. Вычисления производились на операционной системе с вытесняющей многозадачностью (Linux 3.2). Ограничения на использование вычислительных ядер для процесса выставлялись средствами операционной системы.
В результате вычислений по базовому алгоритму были получены значения центров кластеров 1.37360 и 1.34875 соответственно для первого и второго классов изображений, а по модифицированному алгоритму при K = 4 - значения центров кластеров 1.37326 и 1.34804, при K =16 - значения центров кластеров 1.35829 и 1.34103.
Алгоритмы сравнивались по двум факторам - проценту точности и скорости работы. Процент точности рассчитывался как отношение числа верно отнесенных к классам изображений к общему их числу. Результаты представлены в табл. 1 и 2.
Таблица 1. Сравнение алгоритмов классификации. Процент точности
Алгоритм Результат классификации
Верно Ошибочно Точность, %
Базовый 24 1 96
Модифицированный К = 4 23 2 92
К = 16 19 6 76
Таблица 2. Сравнение алгоритмов классификации. Время вычислений
Алгоритм Число ядер Время вычислений, мин
Базовый 1 11:45
1 4:49
К = 4 2 2:37
Модифицированный 4 1:22
1 3:15
К = 16 2 1:59
4 1:08
В результате численного эксперимента было установлено, что при вычислениях по модифицированному алгоритму точность падает с 96 до 92% при K = 4 и до 76% при K = 16. При этом скорость вычислений по модифицированному алгоритму c K = 4 даже при использовании 1 ядра процессора выше в 2.4 раза, а 4 ядер - выше в 8.5 раз. При K =16 скорость работы алгоритма на 4 ядрах более чем в 10 раз превосходит скорость работы базового алгоритма. Стоит отметить, что для K =16 необходимая степень параллельности вычислений не могла была быть достигнута вследствие технических ограничений используемого компьютера (число ядер которого равно 4). Есть все основания полагать, что при вычислении на компьютере с большим числом ядер (8 или 16) скорость могла бы еще возрасти.
Заключение. Эксперимент показывает, что модифицированный алгоритм при K = 4 можно применять для классификации исследуемых изображений практически в той же мере, что и немодифицированный. При этом точность в использованной выборке падает всего на 4%. При этом при K =16 процент точности падает значительно. Можно предположить, что это связано с тем, что число граничных точек, в которых свойство стационарности построенного потока выполняется с меньшей точностью по отношению к общему числу точек (в них условие стационарности выполняется с точностью е), возросло. На указанной выборке падение точности составило 20%. Это означает, что алгоритм с большими значениями K для исследуемых изображений тканей печени можно использовать только как предварительно оценивающий степень близости картинки к тому или иному классу.
Для изображений больших размеров (или большего разрешения) можно ожидать, что для задач классификации могут быть применены и большие значения K, вследствие того, что число граничных точек растет пропорционально линейным размерам изображения, в то время как общее число точек изображения растет квадратично.
Особое внимание нужно уделить тому факту, что скорость работы модифицированного алгоритма выше не только за счет распараллеливания, но и потому, что, как показано выше, он эффективнее алгоритмически.
Литература
1. Marcus Lind. An Introduction to symbolic dynamics and coding. Cambridge: Cambridge University Press, 1995. 495 с.
2. Шелейховский Г. В. Композиция городского плана как проблема транспорта. Минск: Гипро-гор, 1946. 129 с.
3. Брэгман Л. М. Доказательство сходимости метода Г. В. Шелейховского для задачи с транспортными ограничениями // Журн. вычисл. математики и матем. физики. 1967. Т. 7, № 1. C. 147—156.
4. Брэгман Л. М. Релаксационный метод нахождения общей точки выпуклых множеств и его применение для решения задач выпуклого программирования //Журн. вычисл. математики и матем. физики. 1967. Т. 7, № 3. С. 620-631.
5. Ампилова Н. Б, Романовский И. В., Петренко Е. И. О максимизации энтропии при линейных ограничениях // Труды Междунар. науч. конференции «Космос, астрономия и программирование (Лавровские чтения)». СПб.: Математико-механический факультет С.-Петерб. ун-та, 2008. С. 181-185.
6. Ампилова Н. Б. Стационарные процессы на графах и анализ изображений //Компьютерные инструменты в образовании. 2013. № 2. С. 24-29.
7. Fawcett T. An introduction to ROC analysis // Pattern Recognition Letters. 2006. N 27. P. 861874.
References
1. Marcus Lind. An Introduction to symbolic dynamics and coding. Cambridge, Cambridge University Press, 1995, 495 p.
2. Sheleihovsky G. V. Kompozicija gorodskogo plana как problema transporta [City plan composition as a transport problem]. Minsk, Gyprogor Publ., 1946, 129 p. (in Russ.)
3. Bregman L. M. Dokazatel'stvo shodimosti metoda G. V. Shelejhovskogo dlja zadachi s transportnymi ogranichenijami [Proof of Sheleihovsky's method convergence for transport problem]. Appl. math. and math. phys. journ., 1967, vol. 7, no. 1, pp. 147—156. (in Russ.)
4. Bregman L. M. Relaksacionnyj metod nahozhdenija obshhej tochki vypuklyh mnozhestv i ego primenenie dlja reshenija zadach vypuklogo programmirovanija [Convex sets common point findind relaxation method and it's applications for convex programming problems solving]. Appl. math. and math. phys. journ., 1967, vol. 7, no. 3, pp. 620—631. (in Russ.)
5. Ampilova N. B, Romanovskiy I. V., Petrenko E. I. O maksimizacii jentropii pri linejnyh ogranichenijah [About entropy maximization in case of linear restrictions]. Intern. science conference "Space, astronomy and programming". St. Petersburg, Mathematics and mechenics faculty St. Petersburg State University, 2008, pp. 181-185. (in Russ.)
6. Ampilova N. B. Stacionarnye processy na grafah i analiz izobrazhenij [Image analysis and stationary flows on graphs]. Comp. tools in education, 2013, no. 2, pp. 24-29. (in Russ.)
7. Fawcett T. An introduction to ROC analysis. Pattern Recognition Letters, 2006, no. 27, pp. 861874.
Статья рекомендована к печати проф. А. П. Жабко. Статья поступила в редакцию 17 февраля 2015 г.
