Как ускорить вложенный цикл операций поиска?

Question

Я программирую полутонирование изображений для лазерной гравировки . При заданной настройке лазер только включается или выключается, поэтому я могу дать ему двоичные изображения с глубиной 1 бит. Поэтому я преобразую изображения в градациях серого с глубиной 8 бит (от 0 до 255) в двоичные изображения с глубиной 1 бит (от 0 до 1).

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

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

import math
import time

import numpy as np

TONES = [[0, 0,
          0, 0],
         [0, 1,
          0, 0],
         [1, 1,
          0, 0],
         [1, 1,
          0, 1],
         [1, 1,
          1, 1]]

def process_tones():
    """Converts the tones above to the right shape."""
    tones_dict = dict()

    for t in TONES:
        brightness = sum(t)
        bitmap_tone = np.reshape(t, (2, 2)) * 255
        tones_dict[brightness] = bitmap_tone
    return(tones_dict)

def halftone(gray, tones_dict):
    """Generate a new image where each pixel is replaced by one with the values in tones_dict.
    """

    num_rows = gray.shape[0]
    num_cols = gray.shape[1]
    num_tones = len(tones_dict)
    tone_width = int(math.sqrt(num_tones - 1))

    output = np.zeros((num_rows * tone_width, num_cols * tone_width),
                         dtype = np.uint8)

    # Go through each pixel
    for i in range(num_rows):
        i_output = range(i * tone_width, (i + 1)* tone_width)

        for j in range(num_cols):
            j_output = range(j * tone_width, (j + 1)* tone_width)

            pixel = gray[i, j]
            brightness = int(round((num_tones - 1) * pixel / 255))

            output[np.ix_(i_output, j_output)] = tones_dict[brightness]

    return output

def generate_gray_image(width = 100, height = 100):
    """Generates a random grayscale image.
    """

    return (np.random.rand(width, height) * 256).astype(np.uint8)

gray = generate_gray_image()
tones_dict = process_tones()

start = time.time()
for i in range(10):
    binary = halftone(gray, tones_dict = tones_dict)
duration = time.time() - start
print("Average loop time: " + str(duration))

Результат:

Average loop time: 3.228989839553833

Средний циклзанимает 3 секунды для изображения размером 100x100, что кажется длинным по сравнению с функциями OpenCV

Я проверил Как ускорить вложенный цикл Python? и Цикл по пикселям в изображении и я не сразу вижу, как векторизовать эту операцию.

Как я могу ускорить этот вложенный цикл операций поиска?

Answer 1

Хитрость заключается в том, чтобы не выполнять итерацию с такой низкой степенью детализации, как у вас, а вместо этого перекладывать большую часть работы на оптимизированные функции numpy. меньших изображений (назовите их «каналами»), каждое из которых содержит данные для одной из позиций в полутоновой сетке.

Отдельный каналЗатем изображения могут быть сгенерированы простым поиском, который в Numpy мы можем просто сделать путем индексации таблицы поиска с изображением в градациях серого (т.е. LUT[image]).

Таблицы поиска

Допустим, мы определяем «размер плитки» (размер одного полутонового рисунка) и отдельные тоновые плитки следующим образом:

TILE_SIZE = (2, 2) # Rows, Cols

TONES = np.array(
    [[0, 0,
      0, 0],
     [0, 1,
      0, 0],
     [1, 1,
      0, 0],
     [1, 1,
      0, 1],
     [1, 1,
      1, 1]]
    , dtype=np.uint8) * 255

Сначала используем np.linspace для расчета соответствия между оттенками серого и индексами тонов. Затем для каждой позиции мы создаем таблицу поиска из определения тонов (для этого используем технику поиска).

def generate_LUTs(tones, tile_size):
    num_tones, num_tiles = tones.shape
    tile_rows, tile_cols = tile_size
    assert(num_tiles == (tile_rows * tile_cols))

    # Generate map between grayscale value and tone index
    gray_level = np.linspace(0, (num_tones - 1), 256, dtype=np.float32)
    tone_map = np.uint8(np.round(gray_level))

    # Generate lookup tables for each tile
    LUTs = []
    for tile in range(num_tiles):
        LUTs.append(tones[:,tile][tone_map])

    return LUTs

Объединение каналов

Теперь, чтобы объединить каналы вполное выходное изображение.

Первым шагом является reshape каждого изображения канала, так что оно имеет только один столбец.

Затем мы можем объединить все каналыизображения, которые совместно используют одну и ту же строку полутонового рисунка, используя np.hstack.

Далее мы изменим результаты, напримерчто они имеют то же количество строк, что и входное изображение (то есть теперь у них будет вдвое больше столбцов).

Мы снова объединяем все измененные изображения, используя np.hstack.

Наконец, мы изменим результат так, чтобы он имел правильное количество строк (в соответствии с размером плитки), и мы закончили.

В коде (обобщенно для любого размера плитки):

def halftone(image, LUTs, tile_size):
    tiles = []
    for tile in range(len(LUTs)):
        tiles.append(LUTs[tile][image])

    image_rows, _ = image.shape
    tile_rows, tile_cols = tile_size

    merged_rows = []
    for row in range(tile_rows):
        row_tiles = tiles[row * tile_cols:(row + 1) * tile_cols]
        merged_row = np.hstack([row_tile.reshape(-1, 1) for row_tile in row_tiles])
        merged_rows.append(merged_row.reshape(image_rows, -1))

    return np.hstack(merged_rows).reshape(image_rows * tile_rows, -1)

---

Пример использования:

LUTs = generate_LUTs(TONES, TILE_SIZE)
binary = halftone(gray, LUTs, TILE_SIZE)

---

Пример вывода:

Ас плитками 3х3:

Answer 2

Эту проблему можно решить очень быстро с помощью чистого numpy.

• Первое вычисление brightness векторным способом.
• Следующий индекс tones с яркостью для преобразования grayдля 4d массива формы HxWx2x2
• используйте np.transpose для реорганизации массива, чтобы чередовать введенные размеры из tones с исходными из gray. Изображение преобразуется в Hx2xWx2
• вертикальные размеры "сглаживание / слияние" (H-от gray и 2 от tone), то же самое для горизонтальных размеров (W от gray, 2 от tone). Эта операция выполняется путем изменения формы (H * 2) x (W * 2)

. Пожалуйста, вставьте следующий код под кодом из вопроса и запустите его.

def process_tones2():
    tones = np.array(TONES, dtype='u1')
    size = int(np.sqrt(tones.shape[-1]))
    tones = 255 * tones.reshape(-1, size, size)
    bins = tones.sum(axis=(-2,-1), dtype=int) // size ** 2
    iperm = np.argsort(bins)
    return bins[iperm], tones[iperm]

def halftone_fast(gray, bins, tones):
    height, width = gray.shape
    tone_height, tone_width = tones.shape[-2:]
    brightness = np.round(gray / 255 * (len(tones) - 1)).astype('u1')
    binary4d = tones[brightness]
    binary4d = binary4d.transpose((0,2,1,3))
    binary = binary4d.reshape(height * tone_height, width * tone_width)
    return binary

bins, tones = process_tones2()
start = time.time()
for i in range(10):
    binary2 = halftone_fast(gray, bins, tones)
duration = time.time() - start
print("Average loop time: " + str(duration))
print("Error:", np.linalg.norm(binary.astype(float) - binary2))

На моей машине я получил следующие результаты:

Average loop time: 2.3393328189849854
Average loop time: 0.0032405853271484375
Error: 0.0

Ускорение примерно в 1000 раз.

Обратите внимание, что аргумент bins не используется в halftone_fast(). Причина в том, что он не нужен для полутонов. Код вопроса работает только в том случае, если TONES образуют линейное пространство уровней яркости, начиная с 0 и заканчивая на всех. Поэтому brightness работает как индекс для отсортированного списка tones.

Если отображение не является линейным, то для вычисления правильных индексов в массиве tones нужно использовать np.digitize(gray, bins).

Answer 3

Ваш алгоритм состоит из двух частей: вычисление «яркости» каждого пикселя и замена пикселей полутоновыми точками.

Сначала предположим, что входное изображение имеет форму ( ч , ш ).

grayscale = np.array(...)
h, w = grayscale.shape

Уровни яркости

Вычисление яркости выполняется в два этапа:

1. Определите границы для каждого уровня яркости. Это может быть достигнуто с помощью np.linspace для разделения диапазона [0, 256) на num_tones равных по размеру кусков.

   bins = np.linspace(0, 256, num_tones + 1)
   # e.g. with 4 tones: [0, 64, 128, 192, 256]
   

2. Определите, какиеуровень каждого пикселя падает. Это может быть достигнуто с помощью np.digitize.

   # (subtract 1 because digitize counts from 1)
   levels = np.digitize(grayscale, bins) - 1  # shape (h, w)
   

   Тогда levels[i, j] это уровень яркости grayscale[i,j] (от 0 до num_tones,включительно).

Полутона

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

tones = np.array(...)  # shape(num_tones, x, y)
x, y = tones.shape[1:]

Используя уровни яркости изображения в качестве индексного массива ¹ для tones, вы получаете полутоновую матрицу каждого пикселя.

halftones = tones[levels]  # shape (h, w, x, y)
# halftones[i, j] is the halftone for grayscale[i, j]

Тогда нужно просто упорядочить элементы в правильном порядке и сгладить массив.

# Reorder axes so halftone rows are before image columns
ordered = halftones.swapaxes(1, 2)  # shape (h, x, w, y)

# Make it 2-dimensional
result = ordered.reshape(h * x, w * y)

Скорость

Я написал скрипт для сравнения скоростей исходного кода, моего ответа и ответа tstanisl . Результаты:

Best times
halftone:      0.346237126000000
np_halftone:   0.000565907715000
halftone_fast: 0.000437084295000

Оба ответа выполняются в несколько сотен (600 для моего, 800 для tstanisl) раз быстрее, чем исходный код, при этом tstanisl работает лучше моего примерно на 30%.

В обмен на эту скорость, моя функция имеет одно незначительное преимущество от tstanisl's и оригинала: если вы хотите использовать пользовательские тоны, у которых нет общих значений, соответствующих их яркости, этот алгоритм все равно будет работать (например, если вы хотите инвертироватьцвета в полутонах). В противном случае tstanisl более эффективен.

---

¹ Последний пример в связанном разделе руководства пользователя Numpy на самом деле очень похож на это - он говорит о сопоставлении значений цвета изображения сRGB утраивается.

Answer 4

Яркость рассчитывается непосредственно векторизованным способом: (gray * ((num_tones - 1) / 255)).round(). Эти значения затем сопоставляются с tones_dict, в результате чего mapped_brightness. Затем эту матрицу необходимо изменить, объединив каждую строку по горизонтали, а затем изменив окончательный результат.

def halftone(gray, tones_dict):
    """Generate a new image where each pixel is replaced by one with the values in tones_dict.
    """

    num_tones = len(tones_dict)
    tone_width = int(math.sqrt(num_tones - 1))

    mapped_brightness = np.array(
        [list(map(tones_dict.__getitem__, row)) 
         for row in (gray * ((num_tones - 1) / 255)).round()],
        dtype=np.uint8)

    output = np.array(
        [np.concatenate(row, axis=1) for row in mapped_brightness]
    ).reshape(*(n * tone_width for n in gray.shape))

    return output