Работа с граничными условиями / гало регионов в CUDA

Question

Я работаю над обработкой изображений с помощью CUDA и сомневаюсь в обработке пикселей.

Что часто делается с граничными пикселями изображения при применении сверточного фильтра m x m?

В сверточном ядре 3 x 3 легче игнорировать 1 пиксельную границу изображения, особенно когда код улучшается с помощью разделяемой памяти. Действительно, в этом случае не нужно проверять, имеет ли данный пиксель всю доступную окрестность (то есть пиксель с координатой (0, 0) не имеет левых, верхних левых, верхних соседей). Однако удаление границы пикселя 1 исходного изображения может привести к частичным результатам.

В противоположность этому, я хотел бы обработать все пикселей в изображении, также при использовании улучшений общей памяти, например, при загрузке 16 x 16 пикселей, но вычислении внутренней 14 x 14. Также в этом случае игнорирование граничных пикселей создает более четкий код.

Что обычно делается в этом случае?

Кто-нибудь обычно использует мой подход, игнорируя граничные пиксели?

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

Заранее спасибо.

Answer 1

Обычный подход к работе с эффектами границы заключается в добавлении исходного изображения дополнительными строками и столбцами в зависимости от размера фильтра.Некоторые общие варианты для дополненных значений:

• Константа (например, ноль)
• Реплицируйте первую и последнюю строку / столбец столько раз, сколько необходимо
• Отражениеизображение на границах (например, столбец [-1] = столбец [1], столбец [-2] = столбец [2])
• Обтекание значений изображения (например, столбец [-1] = столбец [ширина]-1], столбец [-2] = столбец [ширина-2])

Answer 2

tl; dr: Это зависит от проблемы, которую вы пытаетесь решить - для этого нет решения, которое применимо ко всем проблемам. На самом деле, говоря математически, я подозреваю, что «решения» вообще не может быть, поскольку я считаю, что это некорректная проблема, с которой вам приходится сталкиваться.

(заранее извиняюсь за мое безрассудное злоупотребление математикой)

Чтобы продемонстрировать, давайте рассмотрим ситуацию, когда все компоненты пикселя и значения ядра предполагаются положительными. Чтобы получить представление о том, как некоторые из этих ответов могут ввести нас в заблуждение, давайте еще подумаем о простом усредняющем («рамочном») фильтре. Если мы установим значения за пределами границы изображения на ноль, то это явно приведет к снижению среднего значения для каждого пикселя в пределах ceil (n / 2) (манхэттенское расстояние) от границы. Таким образом, вы получите «темную» границу на отфильтрованном изображении (при условии наличия одного компонента интенсивности или цветового пространства RGB - ваши результаты будут различаться в зависимости от цветового пространства!). Обратите внимание, что аналогичные аргументы могут быть сделаны, если мы установим значения вне границы для любой произвольной константы - среднее значение будет стремиться к этой константе. Постоянная ноль может быть уместной, если края вашего типичного изображения все равно стремятся к 0. Это также верно, если мы рассмотрим более сложный фильтр ядра как гауссову, однако проблема будет менее выраженной, поскольку значения ядра имеют тенденцию быстро уменьшаться с расстоянием от центра.

Теперь предположим, что вместо использования константы мы решили повторить значения ребер. Это то же самое, что сделать рамку вокруг изображения и копировать строки, столбцы или углы достаточно много раз, чтобы фильтр оставался «внутри» нового изображения. Вы также можете думать об этом как о зажиме / насыщении координат выборки. Это имеет проблемы с нашим простым блочным фильтром, потому что он переоценивает значения краевых пикселей. Набор краевых пикселей появится несколько раз, но все они имеют одинаковый вес w=(1/(n*n)). Предположим, что мы отбираем краевой пиксель со значением K 3 раза. Это означает, что его вклад в среднее значение:

K*w + K*w + K*w  = K*3*w

Так эффективно, что один пиксель имеет больший вес в среднем. Обратите внимание, что, поскольку это средний фильтр, вес является постоянным по отношению к ядру. Однако этот аргумент применяется к ядрам с весами, которые также меняются в зависимости от позиции (опять же: подумайте о ядре Гаусса ..).

Предположим, мы обертываем или отражаем координаты выборки, так что мы все еще используем значения из-за границы изображения. Это имеет некоторые ценные преимущества по сравнению с использованием константы, но также не обязательно является «правильным». Например, сколько фотографий вы делаете, когда объекты на верхней границе похожи на объекты внизу? Если вы не снимаете зеркально-гладкие озера, я сомневаюсь, что это правда. Если вы фотографируете камни для использования в качестве текстур в играх, может подойти упаковка или отражение. Я уверен, что здесь следует сделать важные замечания о том, как упаковка и отражение, вероятно, уменьшат любые артефакты, возникающие в результате использования преобразования Фурье. Однако это возвращается к той же идее: у вас есть периодический сигнал, который вы не хотите искажать, вводя паразитные новые частоты или переоценивая амплитуду существующих частот.

Так что вы можете сделать, если вы фильтруете фотографии ярко-красных скал под голубым небом? Очевидно, что вы не хотите добавлять оранжевую дымку в голубое небо и голубоватый пух на красных скалах. Отражение выборочной координаты работает, потому что мы ожидаем цвета, аналогичные пикселям, найденным в отраженных координатах ... если только ради аргумента мы не представляем, что ядро ​​фильтра настолько велико, что отраженная координата будет простираться за горизонт.

Давайте вернемся к примеру фильтра коробки.Альтернатива этому фильтру - перестать думать об использовании статического ядра и вспомнить, что должно было делать это ядро.Усредняющий / блочный фильтр предназначен для суммирования компонентов пикселя, а затем деления на количество суммированных пикселей.Идея в том, что это сглаживает шум.Если мы хотим снизить эффективность подавления шума вблизи границы, мы можем просто сложить меньшее количество пикселей и разделить на соответственно меньшее число.Это может быть расширено для фильтров с аналогичными терминами, которые я буду называть «нормализующими» - терминами, которые связаны с площадью или объемом фильтра.Для терминов «площадь» вы подсчитываете количество весов ядра, которые находятся в пределах границы, и игнорируете те весы, которые не являются.Затем используйте этот счет как «площадь» (что может потребовать дополнительного умножения).Для объема (опять же: принимая положительные веса!) Просто сложите веса ядра.Эта идея, вероятно, ужасна для производных фильтров, потому что меньше пикселей, чтобы конкурировать с шумными пикселями, а дифференциалы общеизвестно чувствительны к шуму.Кроме того, некоторые фильтры были получены с помощью числовой оптимизации и / или эмпирических данных, а не с помощью ab-initio / аналитических методов, и, таким образом, может отсутствовать очевидный «нормализующий» фактор.

Answer 3

Ваш вопрос несколько широк, и я считаю, что он смешивает две проблемы:

1. , касающийся граничных условий ;
2. , касающийся областей гало .

Первая проблема ( граничные условия ) встречается, например, при вычислении свертки между изображением и ядром 3 x 3.Когда окно свертки пересекает границу, возникает проблема расширения изображения за пределы его границ.

Вторая проблема ( области гало ) встречается, например, при загрузке16 x 16 мозаика в разделяемой памяти, и для обработки производных второго порядка необходимо обработать внутреннюю 14 x 14 плитку.

Для второго вопроса, мне кажется, полезным является следующий вопрос: Анализ доступа к памятиобъединение моего ядра CUDA .

Что касается расширения сигнала вне его границ, в этом случае полезным инструментом является память текстур благодаря различным предоставленным режимам адресациисм. Различные режимы адресации текстур CUDA .

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

#include <stdio.h>

#include "TimingGPU.cuh"
#include "Utilities.cuh"

texture<float, 1, cudaReadModeElementType> signal_texture;

#define BLOCKSIZE 32

/*************************************************/
/* KERNEL FUNCTION FOR MEDIAN FILTER CALCULATION */
/*************************************************/
__global__ void median_filter_periodic_boundary(float * __restrict__ d_vec, const unsigned int N){

    unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;

    if (tid < N) {

        float signal_center = tex1D(signal_texture, tid - 0);
        float signal_before = tex1D(signal_texture, tid - 1);
        float signal_after  = tex1D(signal_texture, tid + 1);

        printf("%i %f %f %f\n", tid, signal_before, signal_center, signal_after);

        d_vec[tid] = (signal_center + signal_before + signal_after) / 3.f;

    }
}


/********/
/* MAIN */
/********/
int main() {

    const int N = 10;

    // --- Input host array declaration and initialization
    float *h_arr = (float *)malloc(N * sizeof(float));
    for (int i = 0; i < N; i++) h_arr[i] = (float)i;

    // --- Output host and device array vectors
    float *h_vec = (float *)malloc(N * sizeof(float));
    float *d_vec;   gpuErrchk(cudaMalloc(&d_vec, N * sizeof(float)));

    // --- CUDA array declaration and texture memory binding; CUDA array initialization
    cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();
    //Alternatively
    //cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat);

    cudaArray *d_arr;   gpuErrchk(cudaMallocArray(&d_arr, &channelDesc, N, 1));
    gpuErrchk(cudaMemcpyToArray(d_arr, 0, 0, h_arr, N * sizeof(float), cudaMemcpyHostToDevice));

    cudaBindTextureToArray(signal_texture, d_arr); 
    signal_texture.normalized = false; 
    signal_texture.addressMode[0] = cudaAddressModeWrap;

    // --- Kernel execution
    median_filter_periodic_boundary<<<iDivUp(N, BLOCKSIZE), BLOCKSIZE>>>(d_vec, N);
    gpuErrchk(cudaPeekAtLastError());
    gpuErrchk(cudaDeviceSynchronize());

    gpuErrchk(cudaMemcpy(h_vec, d_vec, N * sizeof(float), cudaMemcpyDeviceToHost));

    for (int i=0; i<N; i++) printf("h_vec[%i] = %f\n", i, h_vec[i]);

    printf("Test finished\n");

    return 0;
}