Как выполнить логическую операцию и логическое индексирование с помощью VIPS в Python?

У меня были следующие коды, использующие Python и OpenCV. Вкратце, у меня есть стопка снимков, сделанных с разной фокусной глубиной. Коды выбирают пиксели в каждой позиции (x, y), которая имеет наибольший лапласианский отклик гуасиана среди всех фокальных глубин (z), создавая таким образом изображение с наложенным фокусом. Функция get_fmap создает двумерный массив, в котором каждый пиксель будет содержать номер фокальной плоскости, имеющей наибольший логарифмический отклик. В следующих кодах закомментированные строки являются моей текущей реализацией VIPS. Они не выглядят совместимыми в определении функции, потому что это лишь частичное решение.

# from gi.repository import Vips

def get_log_kernel(siz, std):
    x = y = np.linspace(-siz, siz, 2*siz+1)
    x, y = np.meshgrid(x, y)
    arg = -(x**2 + y**2) / (2*std**2)
    h = np.exp(arg)
    h[h < sys.float_info.epsilon * h.max()] = 0
    h = h/h.sum() if h.sum() != 0 else h
    h1 = h*(x**2 + y**2 - 2*std**2) / (std**4)
    return h1 - h1.mean()

def get_fmap(img):    # img is a 3-d numpy array.
    log_response = np.zeros_like(img[:, :, 0], dtype='single')
    fmap = np.zeros_like(img[:, :, 0], dtype='uint8')
    log_kernel = get_log_kernel(11, 2)
    # kernel = get_log_kernel(11, 2)
    # kernel = [list(row) for row in kernel]
    # kernel = Vips.Image.new_from_array(kernel)
    # img = Vips.new_from_file("testimg.tif")
    for ii in range(img.shape[2]):           
        # img_filtered = img.conv(kernel)
        img_filtered = cv2.filter2D(img[:, :, ii].astype('single'), -1, log_kernel)
        index = img_filtered > log_response
        log_response[index] = img_filtered[index]
        fmap[index] = ii
    return fmap

а затем fmap будет использоваться для выделения пикселей из разных фокальных плоскостей для создания изображения с фокусировкой.

Это делается на очень большом изображении, и я считаю, что VIPS может с этим справиться лучше, чем OpenCV. Однако официальная документация предоставляет довольно скудную информацию о его привязке к Python. Из информации, которую я могу найти в Интернете, я могу заставить работать свертку изображений (что в моем случае на порядок быстрее, чем OpenCV). Мне интересно, как это реализовать в VIPS, особенно эти строки?

log_response = np.zeros_like(img[:, :, 0], dtype = 'single')

index = img_filtered > log_response

log_response[index] = im_filtered[index]

fmap[index] = ii

person user3667217    schedule 18.10.2015    source источник


Ответы (2)


log_response и fmap инициализируются как трехмерные массивы в коде вопроса, тогда как в тексте вопроса указано, что результат fmap является двумерным массивом. Итак, я предполагаю, что log_response и fmap должны быть инициализированы как 2D-массивы с их формами, такими же, как у каждого изображения. Таким образом, правки будут -

log_response = np.zeros_like(img[:,:,0], dtype='single')
fmap = np.zeros_like(img[:,:,0], dtype='uint8')

Теперь, возвращаясь к теме вопроса, вы выполняете 2D-фильтрацию для каждого изображения одно за другим и получаете максимальный индекс отфильтрованного вывода для всех сложенных изображений. В случае, если вы не знали, согласно документации cv2.filter2D, его также можно использовать с многомерным массивом, давая нам на выходе многомерный массив. Тогда получить максимальный индекс для всех изображений так же просто, как .argmax(2). Таким образом, реализация должна быть чрезвычайно эффективной и простой -

fmap = cv2.filter2D(img,-1,log_kernel).argmax(2)
person Divakar    schedule 18.10.2015
comment
Спасибо. Я допустил ошибку, приводя этот минимальный пример. Я отредактировал свой пост. Действительно хорошо знать, что cv2.filter2d имеет такое использование. Я обязательно попробую. Однако VIPS по умолчанию использует несколько ядер, а с OpenCV мне придется реализовать это сам. А добавление многопроцессорности в OpenCV также добавит накладных расходов. Я не думаю, что модуль multiprocessing плюс OpenCV будет работать близко к VIPS. Я тестировал фильтрацию по Гауссу на изображении 3000x4000 с ядром 11x11. VIPS завершил задачу за 0,0025 секунды. OpenCV занял 0,1334 секунды. - person user3667217; 18.10.2015
comment
@ user3667217 Итак, вы уже реализовали решение в VIPS? Si так, поделитесь в качестве ответа здесь? Что ж, я не слышал об этом до сегодняшнего дня от вас, так что спасибо! Сообщите мне любые другие мысли по этому поводу! Вот ВИПС, мне бы непременно интересно смотрелся. - person Divakar; 18.10.2015
comment
У меня нет полного решения в VIPS, поэтому я спрашиваю, как выполнять логические операции и индексацию в VIPS. Судя по всему, VIPS активно развивается (на данный момент уже 8-я версия), но гораздо менее известна. То, что мне удалось сделать, не слишком далеко от существующего руководства по питону. Я получил его только для выполнения гауссовой фильтрации и использования моего собственного ядра для фильтрации. Я все еще изучаю это. Я буду коды в моем OP. VIPS действительно потрясающий. Например, он потратил 0,0011 секунды на чтение изображения, в то время как openCV потратил 0,043 секунды. Это просто невероятно быстро. Определенно стоит попробовать! - person user3667217; 18.10.2015
comment
@ user3667217 Сумасшедшие вещи, которыми вы со мной делитесь! Я имею в виду, что использую OpenCV довольно давно. Итак, судя по тому, что вы говорите, этот VIPS кажется новым и определенно захватывающим! Я обязательно попробую когда-нибудь. Но в этом вопросе я не думаю, что смогу помочь с VIPS. Стоит подождать, чтобы увидеть, что у кого-то есть больше идей по этому поводу. - person Divakar; 18.10.2015
comment
Вы когда-нибудь пробовали OpenCV с TBB? Мои приведенные выше результаты получены из однопоточного OpenCV. OpenCV + TBB - одна из альтернатив, но я не могу найти много примеров. Я до сих пор не думаю, что они вместе превосходят VIPS, thuogh. Но, конечно, VIPS гораздо более ограничен с точки зрения количества выполняемых функций. - person user3667217; 18.10.2015
comment
@ user3667217 Нет, пока я не работал ни с каким распараллеливанием на OpenCV, но определенно над картами в ближайшем будущем. Я думаю, что многопроцессорность - это то, с чего я бы начал. Распараллеливание - это одна из областей, которые мне действительно нужно изучить на Python! До сих пор все мои распараллеливания были ограничены NumPy. - person Divakar; 18.10.2015
comment
Главное, что делает vips, - это потоковая передача изображений. Он может загружать, обрабатывать и сохранять все параллельно, передавая изображение через ваш набор ядер ЦП. Это дает огромное ускорение и экономию памяти для больших изображений. Ваш тест gaussblur может быть несправедливым: vips будет использовать здесь разделяемую свертку, и, конечно же, это намного быстрее, чем полная двумерная свертка. - person jcupitt; 19.10.2015
comment
@ user894763 Очень полезная информация, спасибо! Определенно на картах для опробования после вдохновения от вас и OP! - person Divakar; 19.10.2015

После ознакомления с руководством Python VIPS и некоторых проб и ошибок я придумал свой собственный ответ. Моя реализация numpy и OpenCV, о которой идет речь, может быть переведена на VIPS следующим образом:

import pyvips

img = []
for ii in range(num_z_levels):
    img.append(pyvips.Image.new_from_file("testimg_z" + str(ii) + ".tif")

def get_fmap(img)
    log_kernel = get_log_kernel(11,2)  # get_log_kernel is my own function, which generates a 2-d numpy array.
    log_kernel = [list(row) for row in log_kernel]  # pyvips.Image.new_from_array takes 1-d list array.
    log_kernel = pyvips.Image.new_from_array(log_kernel)  # Turn the kernel into Vips array so it can be used by Vips.
    log_response = img[0].conv(log_kernel)

    for ii in range(len(img)):
        img_filtered = img[ii+1].conv(log_kernel)
        log_response = (img_filtered > log_response).ifthenelse(img_filtered, log_response)
        fmap = (img_filtered > log_response).ifthenelse(ii+1, 0)

Логическая индексация достигается ifthenelse методом:

result_img = (test_condition).ifthenelse(value_if_true, value_if_false)

Синтаксис достаточно гибкий. Условием тестирования может быть сравнение двух изображений одинакового размера или изображения и значения, например img1 > img2 или img > 5. Аналогично, value_if_true может быть одиночным значением или изображением Vips.

person user3667217    schedule 19.10.2015
comment
Я сопровождаю vips. Ваш код выглядит хорошо. Вы видели функцию vips logmat? Это может сработать для вас. Создает маску журнала vips.ecs.soton.ac.uk/supported/current/doc/html/libvips/ использовать из Python как (например) Vips.Image.logmat(2, 0.1). Он также будет генерировать приблизительные целочисленные маски и разделяемые маски, что может дать полезное ускорение. Вы можете использовать строку "float" вместо Vips.Precision.FLOAT, если хотите. - person jcupitt; 19.10.2015
comment
(часть 2) Вам нужно изображение index? Я думаю, вы могли бы переписать это как fmap = (img_filtered > log_response).ifthenelse(ii + 1, 0), я, наверное, что-то упускаю. Мне было бы любопытно, как скорость и использование памяти сравнивается с opencv, если у вас есть какие-то тесты, которыми вы могли бы поделиться. - person jcupitt; 19.10.2015
comment
Попробуйте: x = Vips.Image.logmat(2, 0.1); x.matrixprint() увидеть коврик для журнала, который делает vips. Добавьте precision = "float", чтобы получить версию с плавающей запятой. - person jcupitt; 19.10.2015
comment
О, ты сопровождающий !! Спасибо за ваш комментарий. Я только начал использовать VIPS на прошлой неделе и до сих пор не могу найти, где что находится. Я не мог полностью понять уравнение относительно ядра журнала. Он такой же, как у меня? (Я только что разместил это в своем вопросе). Мне не нужно index изображение, поэтому то, что вы предложили, должно быть лучше. - person user3667217; 19.10.2015
comment
В руководстве есть страница со всеми операциями vips: vips.ecs.soton.ac.uk/supported/current/doc/html/libvips/ Я обычно ^ F в этом поиске, не знаю, поможет ли это. - person jcupitt; 19.10.2015
comment
Как только я смогу полностью преобразовать части, покрываемые OpenCV, в VIPS, я определенно буду рад предоставить тесты. Однако здесь есть один быстрый вопрос. Я читал, что VIPS может обрабатывать только изображения размером до 2 ^ 31 пикселей? Это что правильно? - person user3667217; 19.10.2015
comment
Оси изображения ограничены 2 ^ 31, пиксели - квадратом, поэтому 2 ^ 62. Я регулярно обрабатываю слайды размером 200 000 x 200 000 пикселей на своем ноутбуке. Он также может делать эти огромные изображения на 32-битных машинах, хотя в наши дни это менее актуально, слава богу. - person jcupitt; 19.10.2015
comment
Ядро журнала выглядит похоже, но я бы проверил то, что делает vips (очевидно). - person jcupitt; 19.10.2015