Ранее я продемонстрировал полностью векторизованный алгоритм СЧЕТЧИК, который понижает дискретизацию помеченных изображений, находя режим 2x2 патчей без подсчета частот значений пикселей. Здесь я обобщаю СЧЕТЧИК на 3D-изображения, используя патч 2x2x2, и на произвольные коэффициенты субдискретизации. COUNTLESS 2D работал, находя основной пиксель без вычисления его общей частоты. Это было определено путем обнаружения наличия или отсутствия двух совпадающих пикселей в фрагменте изображения 2x2. COUNTLESS 3D расширяет эту логику путем поиска среди восьми вокселей (трехмерных пикселей) совпадений четырех пикселей, затем трех и, наконец, двух. COUNTLESS N обобщает этот метод для любого размера ввода, хотя он быстро становится непрактичным. COUNTLESS 2D и COUNTLESS 3D можно рассматривать как специфичные для приложения реализации COUNTLESS 4 и COUNTLESS 8 соответственно. Бесчисленные варианты реализуются с использованием только векторизованных инструкций, что делает их практичными для реализации на языках программирования, таких как Python, которые имеют библиотеки быстрой линейной алгебры, такие как Numpy, но которые имеют медленные операторы цикла и ветвления. Хотя COUNTLESS 2D может превзойти простые алгоритмы, использование COUNTLESS 3D для приложений с высокой пропускной способностью не рекомендуется, поскольку это намного медленнее, чем стандартный подход. Однако в Python это можно реализовать, не прибегая к дополнительным этапам компиляции, и при этом обеспечить приемлемую производительность.

Мотивация

Хотя средний человек не часто встречается, трехмерные изображения (объемные изображения) широко используются в биомедицинской визуализации. Объемное изображение может быть построено из стопки 2D-изображений, полученных с регулярно увеличивающимися интервалами. В аппаратах МРТ используются магниты для неинвазивного получения изображений срезов мозга, а клеточные биологи часто используют микроскопы с лазерным питанием для сканирования образцов на разной глубине. Полученные ими изображения складываются в стопку в отсортированном порядке, чтобы сформировать окончательное изображение.

Моя интересующая ткань - это мозг, и метод получения этих изображений состоит в том, чтобы очень тонко нарезать ее с помощью машины, похожей на нож для резки продуктов, называемой ультрамикротом. Затем полученные срезы визуализируются с помощью электронного микроскопа и собираются в стопку для создания объемного изображения. Срезы большие и варьируются от десятков тысяч до сотен тысяч пикселей на каждую сторону. В отличие от этого, высококачественный телевизор или компьютерный монитор с разрешением 4K отображает изображения только от трех до четырех тысяч пикселей на каждую сторону, а другая сторона часто вдвое меньше.

Чтобы отобразить эти изображения с помощью потребительского оборудования, обычной практикой является их субдискретизация, то есть создание серии более мелких сводных изображений, которые достоверно представляют лежащее в основе изображение с высоким разрешением. Например, на Google Maps мир отображается как одно (или небольшое количество) сшитых изображений, но при увеличении масштаба изображения с более высоким разрешением выбираются и отображаются только для интересующей области. С изображениями с микроскопа мы можем уменьшить их разрешение, усредняя участки 2x2x2, чтобы создать изображение в восемь раз меньше, и делать это итеративно, пока изображение не станет достаточно маленьким. Однако, как только мы генерируем метки для описания того, какой воксель принадлежит какому нейрону, усреднение не является вариантом, поскольку метки являются не аналоговыми сигналами, а дискретными идентификаторами. Самый подходящий способ суммировать их - выбрать наиболее частое значение в патче.

Хранение субдискретов требует дополнительных затрат, хотя они экспоненциально снижаются с каждым дополнительным уровнем mip. При рекурсивном сокращении 2x2 стоимость хранения и передачи стека субдискретизации на 33% превышает стоимость хранения слоя с полным разрешением.

LISTING 1: Size of an Infinite Stack of Downsamples
Let S = The size of all the downsamples
Let r = The reduction factor of a downsample (e.g. 4 for 2x2)
S = 1 + 1/r + 1/r^2 + 1/r^3 + …
S/r = 1/r + 1/r^2 + 1/r^3 + 1/r^4 + …
S — S/r = 1
S = r / (r-1)

Следовательно, стоимость хранения стека субдискретизации 2x2 составляет не более 4/3 или 133% от стоимости только исходного изображения. Стоимость хранения стека субдискретизации 2x2x2 составляет не более 8/7 или 114,3% от полного разрешения. В некоторых случаях это сокращение может быть заманчивым. Для тех, кто так соблазнен, вопрос сводится к тому, как этого можно достичь без счета?

COUNTLESS 5 - Небольшое расширение 2D-задачи

В основе COUNTLESS 2D лежит идея о том, что если две метки в наборе из четырех совпадают, нет необходимости учитывать какую-либо дополнительную информацию для объявления режима этого патча. Естественно спросить, можно ли сформулировать аналогичный принцип для набора из пяти ярлыков - наименьшего возможного расширения этой проблемы.

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

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

LISTING 2: The Relevant Combinations of Five Labels
Let the five integer labels be denoted A, B, C, D, E
Combinations of three: 5C3 = 5! / 3!2! = 10 combinations
    ABC, ABD, ABE, 
         ACD, ACE, 
              ADE,
         BCD, BCE, 
              BDE,
              CDE
Combinations of two: 5C2 = 5! / 2!3! = 10 combinations
    AB, AC, AD, AE
        BC, BD, BE
            CD, CE
                DE
Combinations of one: 5C1 = 5 combinations
    A, B, C, D, E

Чтобы оценить каждую комбинацию, мы можем обобщить оператор PICK (A, B) из COUNTLESS 2D, чтобы он мог принимать любое количество аргументов:

PICK(A,B,C, ..., N) := A if A == B == C == ... == N else 0    EQN. 1

Это может быть реализовано в псевдокоде Python / numpy («&» - это побитовый оператор AND):

PICK(A,B,C ..., M,N) := A * (A == B & B == C & ... & M == N) EQN. 2

Каждое применение PICK к данной комбинации даст метку или ноль. Если любая из наших комбинаций из трех дает ненулевой результат, это квалифицируется как режим-кандидат. Если есть более одного совпадения, например, если на самом деле есть четыре совпадающих метки, мы можем выбрать любой из возможных режимов произвольно. У короткозамкнутого логического оператора or есть правильная семантика для выбора метки, если таковая существует.

MODE(A,B,C,D,E) := PICK(A,B,C) || PICK(A,C,D) || PICK(A,D,E)
    || PICK(B,C,D) || PICK(B,D,E) || PICK(C,D,E) || PICK(A,B)
    || PICK(A,C) || PICK(A,D) || PICK(A,E) || PICK(B,C) 
    || PICK(B,D) || PICK(B,E) || PICK(C,D) || PICK(C,E)
    || PICK(D,E) || E                                         EQN. 3

Как описано ранее, || оператор может быть реализован так:

LOGICAL_OR(X,Y) := X + (X == 0) * Y                           EQN. 4

Уравнение 3 состоит из семнадцати элементов, которые необходимо вычислить и объединить логическим оператором ИЛИ. Есть ли способ уменьшить количество необходимых элементов? Невозможно существенно снизить временную сложность общего алгоритма, но есть небольшая экономия, которая имеет значение для небольшого количества этикеток. Обратите внимание, что для случая совпадения двух мы можем расширить трюк, используемый в COUNTLESS 2D, чтобы избежать вычисления совпадений последнего элемента. Используя обозначение PQ для обозначения PICK (P, Q), если ни AB, AC, AD, BC, BD, ни CD не совпадают, тогда мы будем вынуждены выбрать один из AE, BE, CE, DE или E , все они идентичны E. Следовательно, мы можем опустить AE, BE, CE и DE, оставив для расчета тринадцать элементов, что дает значительную экономию.

MODE(A,B,C,D,E) := PICK(A,B,C) || PICK(A,C,D) || PICK(A,D,E)
    || PICK(B,C,D) || PICK(B,D,E) || PICK(C,D,E) 
    || PICK(A,B) || PICK(A,C) || PICK(A,D) || PICK(B,C) 
    || PICK(B,D) || PICK(C,D) || E                            EQN. 5

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

Объединение всего этого дает некоторый рабочий код Python для СЧЕТЧИКА 5. Обратите внимание, что это всего лишь ступенька к БЕСЧИСЛЕННОМУ 3D, возможно, лучше названному СЧЁТНОМУ 8, поскольку он решает режим вокселей 2 x 2 x 2.

БЕССЧЕТНОЕ 3D (также известное как БЕСЧИСЛЕННОЕ 8)

БЕССЧЕТНЫЙ 3D, практическое применение этой концепции теперь в пределах досягаемости. Понижение частоты дискретизации фрагмента объемного изображения 2x2x2 эквивалентно нахождению режима восьми целых чисел. Простое требование для большинства или равного количества кандидатов - четыре одинаковых лейбла. Если имеется ничья 4–4, мы можем выбрать одно из совпадений произвольно. Чтобы оценить СЧЕТЧИК 8, мы должны рассмотреть совпадения длины четыре, затем три, затем два.

Уравнение 6 ниже показывает количество вызовов PICK, которые необходимо сделать. Обратите внимание, что термин 7C2, приведенный ниже, является результатом нашего упрощения 8C2 с использованием метода, показанного в разделе БЕСЧИСЛЕННЫЙ 5.

Total PICKs = 8C4 + 8C3 + 7C2
            =  70 +  56 +  21 
            = 147                                             EQN. 6

Обратите внимание, что PICK (A, B, C, D) в 1,5 раза дороже, чем PICK (A, B, C) (шесть операций против четырех) и в 3 раза дороже, чем PICK (A, B) (шесть операций против два). Этот факт будет важен позже.

Ниже представлена ​​реализация COUNTLESS 8. Обратите внимание, что создание 148 PICKS приведет к увеличению требуемой памяти, если мы не будем осторожны. Используя выражения генератора Python 3, мы можем существенно снизить требования к памяти для этой программы, создавая только несколько дополнительных производных изображений за раз.

Как и в случае с COUNTLESS 2D, вывод оператора PICK бессмысленен, если совпадающие метки равны нулю (он возвращает 0 независимо от того, совпадают они или нет), поэтому мы сдвигаем данные вверх на единицу, чтобы учесть нулевые метки, и смещаем вниз в конце.

Динамическое программирование БЕСПЛАТНОЕ 3D

Пока все хорошо, но можно еще немного ускорить процесс. Как показано в уравнении 6, производительность СЧЕТЧИКА 3D примерно пропорциональна количеству операций PICK, которые должны быть вычислены, которые, в свою очередь, равны количеству операций, которые должны быть вычислены. PICK (A, B, C, D) требует шести операций, в то время как PICK (A, B, C) требует четырех, а PICK (A, B) требует двух. Нам также необходимо вычислить три операции на логическое или.

LISTING 3: Number of Operations Required for COUNTLESS 3D
# of COUNTLESS_3D operations 
    = 8C4 PICK(A,B,C,D) + 8C3 PICK(A,B,C) + 7C2 PICK(A,B)
      + (8C4 + 8C3 + 7C2 + 1 - 1) logical ors
    = 8C4 x 6 + 8C3 x 4 + 7C2 x 2
      + (8C4 + 8C3 + 7C2) x 3
    = 686 + 441     // picks + logical or
    = 1,127 operations

Или в более общем плане:

До сих пор, поскольку решение проблемы основывается на более крупных совпадениях перед меньшими, мы сначала вычисляли наибольшие совпадения. Однако более крупные комбинации могут быть построены из более мелких с использованием методов динамического программирования. За счет дополнительной памяти мы можем рассматривать совпадения двух как подзадачи совпадений трех, а совпадения трех как подзадачи совпадений четырех. Поскольку каждый новый слой строится поверх старых за счет добавления одного элемента, этот метод позволит нам снизить коэффициент (r-1) в члене PICK, исключив дублирующуюся работу.

Вытягивая множитель два, мы можем упростить это уравнение до:

Однако это несколько утомляет вычислять и сообщать, поэтому давайте упростим это с помощью некоторых обозначений Ландау. Сумма всех комбинаций от 0 до N равна 2^N. Поскольку комбинации симметричны, и мы суммируем до N / 2, пропуская небольшой объем вычислений, получается экспоненциальная временная сложность 2 ^ (N-1).

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

LISTING 4: Not All Subproblems are Used
Let A, B, C, D, E represent five labels.
Let the notation AB := PICK(A,B) and ABC := PICK(A,B,C) and so on.
Combinations of Length Two: 
    AB, AC, AD, AE
        BC, BD, BE
            CD, CE
                DE
Combinations of Length Three:
   
   ABC, ABD, ABE, ACD, ACE, ADE
                  BCD, BCE, BDE
                            CDE
Notice that all length three combinations are prefixed with length two problems. For example, PICK(A,B,C) can be written as PICK(PICK(A,B), C). However, certain length two subproblems are not used in length three. AE, BE, CE, and DE can never be used as a prefix because E always occurs at the end.
Therefore, substantial memory can be saved by not memoizing subproblems that include the terminal element.

Собирая все это вместе, мы получаем версию COUNTLESS 3D для динамического программирования, которая значительно быстрее, хотя требует больше памяти. Это требует сохранения результатов предыдущей итерации. Это приводит к пиковому потреблению памяти:

Насколько быстрее должна быть динамическая версия COUNTLESS 3D? Используя уравнения 8 и 10:

LISTING 5: Comparing the Ratio of Eqn. 10 and Eqn. 8
Non-Dynamic = 1,127 operations ; from LISTING 3
Dynamic     = 5 * ( 8C4 + 8C3 + 7C2 ) ; from Eqn 9
            = 735 operations
Dynamic / Non-Dynamic = 0.65 or 1.53x faster

Поскольку идеализированное решение динамического программирования требует только 65% работы, мы должны примерно ожидать, что решение динамического программирования будет примерно в 1,5 раза быстрее. На практике я считаю, что это в 1,3–2,3 раза быстрее, вероятно, из-за разницы в скорости операторов.

БЕСЧИСЛЕННЫЙ N - ОБЩИЙ БЕСЧИСЛЕННЫЙ

Легко обобщить то, что мы узнали, формулируя БЕСЧИСЛЕННОЕ 3D. COUNTLESS N - это очень требовательный к памяти и вычислительным ресурсам алгоритм. Я не ожидаю, что он будет эффективен за пределами узкого диапазона параметров. Однако приведенный ниже алгоритм, когда задан коэффициент уменьшения 2x2x2, по производительности соответствует СЧЕТЧИКУ 3D.

Представление

Показатели производительности COUNTLESS 2D сосредоточены на чистой скорости выполнения. Однако для СЧЕТЧИКА 3D и СЧЕТЧИКА ND из-за количества комбинаций также требуются измерения использования памяти.

Чтобы определить, насколько быстр этот алгоритм, был разработан набор для сравнения, который работал на двухъядерном компьютере с частотой 2,8 ГГц, i7 Macbook Pro (примерно 2014 г.), кэш-памятью L2 256 КБ, кэш-памятью третьего уровня 4 МБ. Максимальное объединение, усреднение и шаг были включены для сравнения скорости, хотя они не подходят для этой задачи.

Python 3.6.2 с numpy-1.13.3 и clang-900.0.39.2 использовались для следующих экспериментов.

Протестированные алгоритмы:

  • шаг: выберите каждый второй пиксель.
  • countless3d: реализация countless3d с низким объемом памяти.
  • dynamic_countless3d: реализация функции countless3d для динамического программирования.
  • countless3d_generalized: реализация COUNTLESS N с низким объемом памяти с коэффициентом субдискретизации 2,2,2 для соответствия COUNTLESS 3D
  • countless3d_dynamic_generalized: версия динамического программирования COUNTLESS N с коэффициентом субдискретизации 2,2,2 для соответствия COUNTLESS 3d

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

  • downsample_with_averaging: среднее количество пикселей в окне 2x2.
  • downsample_with_max_pooling: выберите пиксель с максимальным значением в окне 2x2.

Код, используемый для тестирования алгоритмов, можно найти здесь. Приведенные ниже числа были протестированы на ревизии 15c0077.

Пропускная способность

Следующие числа получены в результате пятикратного выполнения алгоритмов, указанных выше, для воксельного куба 512 x 512 x 512, когда куб выделяется с использованием uint8, uint16, uint32 и uint64. Результаты в таблице 1 и на рис. 3 указаны в мегавокселях в секунду (MVx / sec), поскольку обычно требуется обработать объем определенного «физического» размера. Однако таблица 2, рис. 4 и рис. 5 указаны в мегабайтах в секунду (МБ / с), поскольку различия между запусками в основном заключаются в размере памяти.

Четыре варианта COUNTLESS 3D состоят из трехмерной и N-мерной реализации базового алгоритма вместе с их аналогами динамического программирования. Производительность конкретных реализаций очень похожа на характеристики обобщенных реализаций как в базовых, так и в динамических алгоритмах. Все эти сравнения находятся в пределах 5% друг от друга. Обработка этого куба 128 MVx занимает около пяти секунд с использованием алгоритма динамического программирования.

Что касается MVx / сек, для массивов uint8 динамическое программирование дает увеличение скорости примерно в 2,2 раза. Для uint16 и uint32 он ближе к 1,5x. Для uint64 это примерно 1,27x. Таким образом, наблюдается очевидная неутешительная тенденция к снижению производительности по мере увеличения объема памяти, однако это несколько опровергается производительностью, измеренной в МБ / с.

Рисунок 4 показывает, что производительность между uint8 и uint64 на самом деле немного увеличивается, хотя uint64 немного проседает по сравнению с uint32 в динамической реализации. Результаты uint16 кажутся необычными. Они сравнимы по производительности с uint8 и почти вдвое выше для uint32. Я подозреваю, что либо внутри библиотеки, либо на моем чипе uint16s обрабатываются как uint32s. Улучшение, связанное с динамическим программированием, является существенным, но в следующем разделе мы увидим, что за него приходится платить.

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

Использование памяти

Память, необходимая для хранения подзадач в реализации динамического программирования COUNTLESS 3D, очевидно, значительна. Ниже я пять раз измерил потребность в памяти для обработки куба размером 128 МБ 512x512x512. На рисунке 6 показан действующий базовый алгоритм использования, а на рисунке 7 показан алгоритм динамического программирования. Следующие ниже графики были созданы с использованием Python 3.4 в Ubuntu 14.04 Linux с использованием инструмента mprof.

Базовый алгоритм довольно последовательно использует около 350 МБ ОЗУ, что примерно в 2,7 раза больше памяти, необходимой для хранения куба. Алгоритм динамического программирования требует около 1150 МБ, или примерно в 9 раз больше объема памяти куба.

Обсуждение

СЧЁТНОЕ 3D - это «медленный» алгоритм. Вариант более быстрого динамического программирования (DP) имеет экспоненциальную временную сложность O(2^(N-1)), в то время как базовый вариант несколько медленнее этого. Вариант DP имеет комбинаторно растущую пространственную сложность. Однако, по исключительной удаче, частные случаи 2x2 (N = 4) и 2x2x2 (N = 8) обычно используются и достаточно малы, чтобы быть практичными даже в этих ужасных условиях. Было показано, что случай 2x2 в некоторых случаях способен превзойти стандартный алгоритм подсчета, который был реализован как O (N²), хотя его можно было уменьшить до O (N) за счет использования хэш-карт (хотя есть вероятность, что он будет медленнее для малых N из-за увеличения постоянных накладных расходов).

Увеличение скорости алгоритма DP в 1,5–2,2 раза достигается за счет высокой стоимости памяти. Базовый вариант COUNTLESS 3D потребляет менее трети памяти по сравнению с вариантом DP. Заплатите волынщику в 3 раза больше памяти за удвоение скорости. К счастью, даже в случае uint64 (~ 8–10 ГБ) он все еще находится в разумном диапазоне для установок потребительского уровня для одного процесса.

Можно создать код C или Cython, который работает во много раз быстрее, чем COUNTLESS 3D, почти без использования дополнительного места. Тем не менее, COUNTLESS 3D может служить своеобразной нише: запуск 3D даунсэмплинга для небольших наборов данных с использованием обычных установок Python / numpy без дополнительной компиляции. Этот 3D-алгоритм может помочь исследователям и программистам, которым нужна библиотека, которая будет надежно развернута везде, где это делает Python / numpy.

Один из вариантов использования, при котором такие исследователи могут получить максимальную производительность, - это понижающая дискретизация двоичных масок интересующих областей. Некоторые виды семантической сегментации могут использовать небольшое количество меток. В этом случае массивы uint8 предоставляют более чем достаточно меток и выигрывают от скорости обработки ~ 20 MVx / сек. В коннектомике более типичными требованиями являются uint16 для небольших объемов и uint32 или uint64 для крупных проектов. В таких ситуациях я рекомендую скомпилированный код.

Основная причина, по которой существует эта статья, заключается в том, что мне было интересно исследовать полное обобщение БЕССЧЕТНОГО 2D. Спасибо за чтение. 😇

Будущая работа

В будущем я планирую показать сравнения между реализациями C. Я ожидаю, что он будет работать с частотой ›100 МВx / сек. Некоторый код C существует на Github, но я еще не уверен, что он правильный.

Благодарности

Несколько человек предоставили полезные советы и помощь в разработке бесчисленного количества 3D и в написании этой статьи. Кисук Ли предоставил визуализацию трехмерных данных, показанных на рисунках 1 и 2. Крис Джордан предоставил начальную часть реализации C counting и бесчисленное множество. Д-р Джордж Надь предложил обобщить БЕСПЛАТНОЕ 2D. Джереми Мэйтин-Шепард из Google изначально разработал код Python для striding и downsample_with_averaging для использования с нейрогланером. Особая благодарность Seung Lab за предоставленные метки нейронной сегментации.

Код

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

11 марта 2019 г .: теперь доступен в составе tinybrain пакета PyPI.

Эта статья является продолжением статьи БЕССЧЕТНОЕ - 2-кратное высокопроизводительное даунсэмплинг помеченных изображений с использованием Python и Numpy (2017 г.).