<!DOCTYPE article PUBLIC "-//NLM//DTD JATS (Z39.96) Journal Archiving and Interchange DTD v1.0 20120330//EN" "JATS-archivearticle1.dtd">
<article xmlns:xlink="http://www.w3.org/1999/xlink">
  <front>
    <journal-meta />
    <article-meta>
      <title-group>
        <article-title>Исследование параллельных CUDA- реализаций алгоритма построения карты диспарантности по разноракурсным изображениям*</article-title>
      </title-group>
      <pub-date>
        <year>2016</year>
      </pub-date>
      <fpage>574</fpage>
      <lpage>581</lpage>
      <abstract>
        <p>В статье решается задача повышения быстродействия алгоритма формирования по разноракурсным изображениям карты диспарантности, по которой затем строится модель трехмерной сцены. Основным наиболее затратным этапом алгоритма является определение относительных сдвигов точек на разноракурсных изображениях. Ранее авторами был разработан эффективный алгоритм построения карт диспарантности [1], в котором для повышения быстродействия и надежности используется процедура формирования пирамиды изображений. Настоящая работа посвящена исследованию быстродействия и эффективности реализации соответствующего параллельного алгоритма в CUDA-среде. Ключевые слова: цифровая обработка изображений, реконструкция 3D-сцен по разноракурсным изображениям, сопоставление изображений, CUDA-технология.</p>
      </abstract>
    </article-meta>
  </front>
  <body>
    <sec id="sec-1">
      <title>-</title>
      <p>схема основных этапов предложеной технологии реконструкции 3D-сцены по разноракурсным
изображениям приведена на рисунке 1.
(1)
(2)
(3)
Рис. 1. Схема этапов технологии
Введение дополнительного этапа начального приближенного совмещения разноракурсных
изображений позволяет существенно повысить быстродействие этапа сопоставления с
использованием «пирамиды изображений». Алгоритм пирамиды изображений строится в виде
иерархической схемы вычислений. Пирамида изображений формируется в виде набора изображений,
получаемых уменьшением разрешения в два раза по обеим координатам так, что на N -м уровне
пирамиды формируется изображение, разрешение которого в 2N раз меньше исходного
разрешения.</p>
      <p>Число уровней пирамиды задается в виде параметра. На первом шаге данного этапа
обрабатывается изображение с наименьшим разрешением и начальным сдвигом равным нулю. На
следующем шаге используется информация о сдвиге, найденном на предшествующем шаге так, что
на каждом следующем шаге значения координат удваиваются. Найденные таким образом
соответствующие точки используются для определения фундаментальной матрицы.</p>
      <p>На следующем этапе осуществляется поиск соответствующих точек для построения карты
диспарантности. Используется метод, основанный на использовании в качестве штрафного
коэффициента в минимизируемой функции расстояния до эпиполярной линии [4]. Существо
метода сводится к следующему.</p>
      <p>Пусть координаты точек на первом изображении u, v , а координаты соответствующих им
точек на втором – u  u, v  v , где u , v – относительные сдвиги координат u , v
соответственно, а</p>
      <p>I u, v  , I u  u,v  v
где
au,v  wc  wd  wf ,
wd  exp u0 , v0   u,v 2 ,
– функции распределения яркости отсчетов на этих изображениях.</p>
      <p>Задача состоит в поиске для каждой точки u, v на первом изображении соответствующей
точки u  u, v  v на втором изображении посредством минимизации критерия сходства:
E u0 ,v0 , u, v 
</p>
      <p>
u,vDu0 ,v0 </p>
      <p>
        a u, v I u, v  I u  u, v  v ,
где Du0 ,v0  – заданная область вокруг точки u0,v0  , а a u, v – весовая функция, задаваемая
в указанной области в виде произведения трех коэффициентов:
wd  exp I u0 , v0   I u, v 2 ,
(4)
 au  bv  c 
wd  exp   (5)
 a2  b2 
Коэффициент w f имеет смысл функции «штрафа» при удалении точки u, v от
эпиполярной линии au  bv  c , определяемой по координатам точки u0,v0  с использованием
фундаментальной матрицы F [
        <xref ref-type="bibr" rid="ref12 ref13">11,12</xref>
        ]:
a  u0F11  v0F12  F13 ,
b  u0F21  v0F22  F23 ,
c  u0F31  v0F32  F33.
С использованием полученных значений относительных сдвигов соответствующих точек на
разноракурсных изображениях формируется карта диспарантности (этап 5).
      </p>
      <p>Наибольшие вычислительные затраты в приведенной общей схеме алгоритма имеют место
на этапе 2 – предварительный поиск сдвигов для нахождения фундаментальной матрицы и на
этапе 4 – нахождение итоговых относительных сдвигов. Эти этапы обладают большим
внутренним параллелизмом. В частности, нахождение относительных сдвигов для каждой точки  x, y
можно выполнять независимо. При этом количество параллельных процессов равно
произведению числа пикселей изображения и числа всех возможных сдвигов u, v в области поиска D2 .</p>
      <p>В настоящей работе проводится исследование двух вариантов распараллеливания,
различающихся числом дескрипторов, вычисляемых в каждой отдельной нити и, как следствие,
накладными расходами на пересылки.
2. Описание вариантов параллельного алгоритма</p>
      <p>Варианты параллельного алгоритма реализуются на втором и четвертом этапах сквозной
технология построения карты диспарантности, показанной на рис.1. Отличие четвертого этапа
от второго заключается в использовании штрафного коэффициента для уточнения
относительных сдвигов. Использование штрафного коэффициента практически не влияет на
вычислительную сложность алгоритма. Поэтому, построение параллельного алгоритма рассмотрено на
примере второго этапа.</p>
      <p>Под первым вариантом параллельного алгоритма будем понимать параллельную
CUDA-реализацию алгоритма, которая использовалась в статье [1]. На рисунке 2 приведена структурная
схема этой CUDA-реализации.</p>
      <p>Перед запуском CUDA-ядра копируются необходимые данные из оперативной памяти в
память (global memory) графической видеокарты. Для результатов также выделяется память на
видеокарте. Каждая нить рассчитывает одну евклидову норму для двух выбранных дескрипторов.
Под дескриптором данной точки понимается вектор признаков фрагмента изображения с
центром в искомой соответствующей точке. Для нахождения соответствующей точки, необходимо
рассчитать евклидовы нормы для всех дескрипторов точек, выбранных в области поиска в
качестве возможных соответствующих. При этом число создаваемых нитей равно произведению
размера изображения (в пикселах) и размера области поиска (также в пикселах).
Копирование данных на GPU ( , ). Создание на GPU массива всех евклидовых</p>
      <p>норм дескрипторов ([размер окна поиска]*[размер ])
Блок вычисления</p>
      <p>Thread(0)
Блок вычисления</p>
      <p>Thread(1)</p>
      <p>Блок вычисления
Thread([размер окна
по</p>
      <p>иска]*[размер ])
Копирование массива all_distances на CPU
Нахождение наименьшего значения евклидовой нормы для каждого дескриптора
пикселя на и дескрипторов точек в области поиска на .</p>
      <p>Получение относительного сдвига для каждого пикселя
Рис. 2. Структурная схема первого варианта параллельного алгоритма
На рисунке 3 приведена укрупненная схема вычисления, выполняемых одной нитью.
Данный блок вычислений обозначен на рисунке 2 обозначен как «Блок вычислений Thread». Каждая
нить выполняет вычисления независимо от остальных нитей. Данные, которые используются для
вычислений, находятся в памяти видеокарты (global memory). Таким образом, результатом
работы нити является евклидова норма для двух выбранных дескрипторов.
Недостатком первого варианта параллельного алгоритма является добавление этапа выбора
минимальной евклидовой нормы вектора дескрипторов. Данный этап выполняется
последовательно на CPU. При этом предварительно копируются данные массива с вычисленными
евклидовыми нормами из видеопамяти. После определения минимальной нормы на CPU определяется
относительный сдвиг.</p>
      <p>Структурная схема второго варианта параллельного алгоритма представлена на рисунке 4.
Копирование данных на GPU ( , ). Создание на GPU массива всех
евкли</p>
      <p>довых норм дескрипторов ([размер ])
Блок вычисления</p>
      <p>Thread(0)
Блок вычисления</p>
      <p>Thread(1)
Блок вычисления</p>
      <p>Thread([размер ])
Копирование массива сдвигов на CPU из видеопамяти
Получение относительного сдвига для каждого пикселя наCPU
Рис. 4. Структурная схема второго варианта параллельного алгоритма
В этом варианте (в отличие от предыдущего) удается избежать формирования массива
избыточных данных на GPU реализации. Также отсутствует последовательный этап сравнений
евклидовых норм на CPU.
Предлагаемые параллельные реализации отличаются степенью параллелизма. Во второй
CUDA-реализации используется более высокий уровень распараллеливания. Нить, во втором
варианте параллельного алгоритма, осуществляет больший объем вычислений, чем нить в первом
варианте. На рисунке 5 приведена укрупненная схема блока вычислений, которые выполняются
одной нитью.
3. Результаты экспериментов</p>
      <p>При проведении экспериментов использовались стереоизображения из набора изображений
«Tsukuba», которые часто используются в качестве тестовых в задаче сопоставления
изображений. Этот выбор был продиктован тем, что указанная база изображений содержит также
эталонные карты диспарантности, по которым возможно сопоставление результатов. Исходные
изображения представлены на рисунке 6.</p>
      <p>а)
Рис. 6. Исходные изображения «Tsukuba»
С помощью предлагаемой технологии для указанных изображений была сформирована
карта диспарантности (рис. 7, а). На рисунке 7, б) для сравнения приведена эталонная карта
диспарантности, рассчитанная по тем же изображениям с использованием априорной информации
о параметрах камер и хранящаяся в наборе изображений «Tsukuba».
В качестве исходных данных использовались изображения размером 640×480 пикселей. В
таблице 1 приведены результаты для различных размеров изображений: 20×15, 40×30, 80×60,
160×120, 320×240 и 640×480 пикселей.</p>
      <p>Результаты для последовательного алгоритма и первого варианта параллельного алгоритма
отличаются от результатов приведенного в работе [1]. Связано это с использованием другой
вычислительной системы. Кроме того, в статье [1] не учитывался этап сравнения евклидовых норм
дескрипторов.</p>
      <p>Время выполнения второго параллельного алгоритма больше, чем время выполнения
последовательного для изображений размером 20×15 пикселей. Однако, для остальных изображений
пирамиды получено ускорение второго параллельного алгоритма относительно
последовательного. Было установлено, что время выполнения второго варианта параллельного алгоритма
изображений размером 20×15, 40×30 и 80×60 пикселей примерно одинаковое. Это означает, что
ресурсы видеокарты используются неэффективно. Однако, для изображений размером 160×120,
320×240 и 640×480 пикселей получены примерно одинаковые значения ускорения второго
варианта параллельного алгоритма относительно последовательного, более чем в 15 раз. Тот факт,
В таблице 1 приведены результаты сравнительных исследований времени реализации
алгоритма на CPU и GPU при различном числе уровней пирамиды, задаваемых на первом этапе
устранения больших относительных сдвигов (данные получены при реализации с
использованием GeForce GTX 750 Ti, и Intel Core i7-6700K, 16GB DDR4, ОС Windows 10).
Таблица 1. Полученные результаты исследования</p>
      <p>1 (20×15) 2 (40×30) 3 (80×60) 4(160×120)5(320×240) 6(640×480)
Уровень пирамиды
(разрешение изображения в пикселях)
Время выполнения на CPU (мс)
Время выполнения на GPU
параллельной реализации №1 (мс)
Время выполнения этапа сравнения (мс)
Время выполнения параллельной
реализации №1 и этапа сравнения (мс)
Ускорение (с учетом этапа сравнения)
Время выполнения на GPU
параллельной реализации №2 (мс)
Ускорение
5.3
2.3
0.02
2.32
2.28
9.2
что ускорение для разных изображений остается практически неизменным означает
эффективное использование ресурсов-мощностей видеокарты (CUDA occupancy).</p>
      <p>Время выполнения первого варианта параллельного алгоритма меньше времени выполнения
второго для изображений размером 20×15 и 40×30 пикселей. Для изображений размером 80×60
пикселей время выполнения двух вариантов параллельного алгоритма практически одинаково.
Однако, время обработки с использованием второго варианта для изображений размером
160×120, 320×240 и 640×480 пикселей, оказалось меньше, чем с использованием первого
варианта.
4. Заключение</p>
      <p>Результаты исследований показали, что использование первой параллельной
CUDA-реализации алгоритма [1] позволяет получить ускорение по сравнению с последовательной
реализацией более чем в 12 раз. Однако, из-за избыточного хранения данных на GPU, применение
данного подхода на больших изображениях, например, космических снимках, имеет ограничения,
связанные с размером памяти видеокарты. Кроме того, приходится дополнительно хранить и
обрабатывать больший объем информации в оперативной памяти. Использование второго варианта
параллельного алгоритма (рисунок 4) позволило достичь ускорения в 17 раз по сравнению с
последовательным алгоритмом для изображений размером 640×480, и в 1,5 раза по сравнению с
первым вариантом параллельного алгоритма.
Литература
Research of parallel CUDA implementations of the algorithm for
constructing the disparity map from stereo images
S.P. Korolyov Samara State Aerospace University1, Image Processing Systems Institute of the</p>
      <p>RAS2
This paper solves the problem of increasing the efficiency of the algorithm for the disparity
map formation from stereo images. A three-dimensional model of the scene is reconstructed
using this disparity map. The most computationally complex stage of the algorithm is
determination of the relative points’ shifts in stereo images. In previous paper we proposed an
efficient algorithm for disparity map formation in which a pyramid of images was formed to
improve the efficiency and reliability. This paper is dedicated to the study of the efficiency
of the corresponding parallel algorithm implementation in CUDA environment.</p>
      <p>Keywords: digital image processing, 3D-scene reconstruction, image matching, CUDA.</p>
    </sec>
  </body>
  <back>
    <ref-list>
      <ref id="ref1">
        <mixed-citation>
          <string-name>
            <surname>Kotov</surname>
            <given-names>AP</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Fursov</surname>
            <given-names>VA</given-names>
          </string-name>
          ,
          <article-title>Goshin YeV Technology for fast 3d-scene reconstruction from stereo images // Computer optics</article-title>
          .
          <source>2014</source>
          . Vol
          <volume>39</volume>
          №4 P.
          <fpage>600</fpage>
          -
          <lpage>605</lpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref2">
        <mixed-citation>
          <string-name>
            <surname>Hartley R.I. Theory</surname>
          </string-name>
          and Practice of Projective Rectification //
          <source>International Journal of Computer Vision</source>
          .
          <year>1999</year>
          . Vol.
          <volume>35</volume>
          . P.
          <volume>115</volume>
          -
          <fpage>127</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref3">
        <mixed-citation>
          3.
          <string-name>
            <surname>Fursov</surname>
            <given-names>VA</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Goshin</surname>
            <given-names>YeV</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Bibikov</surname>
            <given-names>SA</given-names>
          </string-name>
          .
          <article-title>3D-scene stereo reconstruction on sheaves of epipolar planes</article-title>
          // Mechatronics automation control
          <year>2013</year>
          . №.
          <volume>9</volume>
          (
          <issue>150</issue>
          ). P.
          <volume>19</volume>
          -
          <fpage>24</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref4">
        <mixed-citation>
          <string-name>
            <surname>Fursov</surname>
            <given-names>VA</given-names>
          </string-name>
          , Goshin,
          <string-name>
            <surname>YeV.</surname>
          </string-name>
          <article-title>Information technology for digital terrain model reconstruction from stereo images // Computer optics</article-title>
          .
          <source>2014</source>
          . Vol
          <volume>38</volume>
          №2 P.
          <fpage>335</fpage>
          -
          <lpage>342</lpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref5">
        <mixed-citation>
          5.
          <string-name>
            <surname>Lowe</surname>
            <given-names>DG</given-names>
          </string-name>
          .
          <article-title>Object recognition from local scale-invariant features //Computer vision</article-title>
          ,
          <year>1999</year>
          .
          <source>The proceedings of the seventh IEEE international conference on. Ieee</source>
          ,
          <year>1999</year>
          . Т. 2. С.
          <volume>1150</volume>
          -
          <fpage>1157</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref6">
        <mixed-citation>
          6.
          <string-name>
            <surname>Bay</surname>
            <given-names>H</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Ess</surname>
            <given-names>A</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Tuytelaars</surname>
            <given-names>T</given-names>
          </string-name>
          , Van Gool L.
          <article-title>Speeded-up robust features (SURF)</article-title>
          .
          <source>Computer vision and image understanding</source>
          .
          <source>2008</source>
          . Vol.
          <volume>110</volume>
          . №. 3. P.
          <volume>346</volume>
          -
          <fpage>359</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref7">
        <mixed-citation>
          7.
          <string-name>
            <surname>Fursov</surname>
            <given-names>VA.</given-names>
          </string-name>
          , Goshin,
          <string-name>
            <surname>YeV.</surname>
          </string-name>
          <article-title>Conformed identification in corresponding points detection problem [In Russian]</article-title>
          .
          <source>Computer optics 2012;</source>
          Vol
          <volume>36</volume>
          №1 p.
          <fpage>131</fpage>
          -
          <lpage>135</lpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref8">
        <mixed-citation>
          8.
          <string-name>
            <surname>Fursov</surname>
            <given-names>VA.</given-names>
          </string-name>
          , Goshin, YeV.
          <article-title>Solving a camera autocalibration problem with a conformed identification method [In Russian]</article-title>
          .
          <source>Computer optics 2014</source>
          . Vol
          <volume>36</volume>
          №4 P.
          <fpage>605</fpage>
          -
          <lpage>610</lpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref9">
        <mixed-citation>
          <string-name>
            <given-names>Harris C.</given-names>
            ,
            <surname>Stephens</surname>
          </string-name>
          <string-name>
            <surname>M.</surname>
          </string-name>
          <article-title>A combined corner and edge detector</article-title>
          .
          <source>Alvey vision conference</source>
          .
          <year>1988</year>
          .
        </mixed-citation>
      </ref>
      <ref id="ref10">
        <mixed-citation>
          Vol.
          <volume>15</volume>
          . С.
          <volume>50</volume>
          .
        </mixed-citation>
      </ref>
      <ref id="ref11">
        <mixed-citation>
          10.
          <string-name>
            <surname>Tao</surname>
            <given-names>M</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Bai</surname>
            <given-names>J</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Kohli</surname>
            <given-names>P</given-names>
          </string-name>
          , Paris S. SimpleFlow:
          <article-title>A Non-iterative, Sublinear Optical Flow Algorithm</article-title>
          . Computer Graphics Forum.
          <year>2012</year>
          . Vol.
          <volume>31</volume>
          . No. 2. P.
          <volume>345</volume>
          -
          <fpage>353</fpage>
          .
        </mixed-citation>
      </ref>
      <ref id="ref12">
        <mixed-citation>
          11.
          <string-name>
            <surname>Forsyth</surname>
            <given-names>D</given-names>
          </string-name>
          ,
          <string-name>
            <surname>Ponce J. Computer</surname>
          </string-name>
          <article-title>Vision: A Modern Approach</article-title>
          . Moscow.: “Williams” Publisher,
          <year>2004</year>
          . 928 p.
        </mixed-citation>
      </ref>
      <ref id="ref13">
        <mixed-citation>
          12.
          <string-name>
            <surname>Gruzman IS</surname>
          </string-name>
          . et al.
          <article-title>Digital image processing in information systems</article-title>
          . NGTU, Novosibirsk.
          <year>2002</year>
          . 352 p.
        </mixed-citation>
      </ref>
    </ref-list>
  </back>
</article>