Лекторы: Боресков А.В. (ВМиК МГУ) Харламов А.А. (NVidia) Иерархия памяти CUDA. Текстуры в CUDA. Цифровая обработка сигналов

TEXTURE

2011 Типы памяти в CUDA Тип памятиДоступУровень выделения Скорость работы Регистры R/WPer-threadВысокая(on-chip) Локальная R/WPer-threadНизкая (DRAM) Shared R/WPer-blockВысокая(on-chip) Глобальная R/WPer-gridНизкая (DRAM) Constant R/OPer-gridВысокая(L1 cache) Texture R/OPer-grid [-] Низкая(DRAM) [+] L1 cache

2011 Архитектура Tesla 10 TPC Interconnection Network ROP L2 ROP L2 ROP L2 ROP L2 ROP L2 ROP L2 ROP L2 ROP L2 DRAM Work Distribution CPUBridgeHost Memory

2011 TEX SM Texture Processing Cluster SM Архитектура Tesla Мультипроцессор Tesla 10 Streaming Multiprocessor Instruction $Constant $ Instruction Fetch Shared Memory SFU SP SFU SP Double Precision Register File

2011 Texture в 3D В CUDA есть доступ к fixed-function HW: Texture Unit 0,01,0 1,10,1 u v 0,0 1,0 0.66,1 0,01,0 0.66,1

2011 Texture HW Латентность больше, чем у прямого обращения в память – Дополнительные стадии в конвеере: Преобразование адресов Фильтрация Преобразование данных Но зато есть кэш – Разумно использовать, если: Объем данных не влезает в shared память Паттерн доступа хаотичный Данные переиспользуются разными потоками

2011 Texture HW Нормализация координат: – Обращение по координатам, которые лежат в диапазоне [0,1] 0,01,0 1,10,1 u v

2011 Texture HW Преобразование координат: – Координаты, которые не лежат в диапазоне [0,1] (или [w, h]) 0,01,0 1,10,1 u v Clamp: -Координата «обрубается» по допустимым границам Wrap - Координата «заворачивается» в допустимый диапозон [1.1, 0.7] [1.1, 0.3] [1.0, 0.7] [0.1, 0.3]

2011 Texture HW Фильтрация: – Если вы используете float координаты, что должно произойти если координата не попадает точно в texel? Point: -Берется ближайший texel Linear: - Билинейная фильтрация Приблизимся вот в эту точку uv 0,01,0 1,10,1

2011 Texture HW Фильтрация Point: -Берется ближайший texel [0.1, 0.3] u v

2011 Texture HW Фильтрация Linear: -Билинейная фильтрация (P 1 *(1-α)+P 2 *(α)) * (1-β) + (P 3 *(1-α)+P 4 *(α)) * (β) [0.1, 0.3] u v P1P2 P3P4 α β

2011 Texture в CUDA Преобразование данных: – cudaReadModeNormalizedFloat : Исходный массив содержит данные в integer, возвращаемое значение во floating point представлении (доступный диапазон значений отображается в интервал [0, 1] или [-1,1]) – cudaReadModeElementType Возвращаемое значение то же, что и во внутреннем представлении

2011 Texture в CUDA (cudaArray) Особый контейнер памяти: cudaArray Черный ящик для приложения Позволяет организовывать данные в 1D/ 2D/3D массивы данных вида: – 1/2/4 компонентные векторы – 8/16/32 bit signed/unsigned integers – 32 bit float – 16 bit float (driver API) Доступ по семейству функций tex1D()/tex2D()/tex3D()

2011 Texture в CUDA (cudaArray) Особенности текстур: – Обращение к 1D / 2D / 3D массивам данных по: Целочисленным индексам Нормализованным координатам – Преобразование адресов на границах Clamp Wrap – Фильтрация данных Point Linear – Преобразование данных Данные могут храниться в формате uchar4 Возвращаемое значение – float4

2011 Texture в CUDA (linear) Можно использовать обычную линейную память Ограничения: – Только для одномерных массивов – Нет фильтрации – Доступ по целочисленным координтам – Обращение по адресу вне допустимого диапазона возвращает ноль

2011 DEVICEHOST Texture в CUDA cudaArray textureReference tex1D() tex2D() tex3D() cudaBindTextureToArray Kernel2 >> devPtr cudaBindTexture tex1Dfetch() Kernel1 >>

2011 Texture в CUDA (linear) texture g_TexRef; __global__ void kernel1 ( float * data ) { int idx = blockIdx.x * blockDim.x + threadIdx.x; data [idx] = tex1Dfetch(g_TexRef, idx); } int main(int argc, char ** argv) { float *phA = NULL, *phB = NULL, *pdA = NULL, *pdB = NULL; for (int idx = 0; idx < nThreads * nBlocks; idx++) phA[idx] = sinf(idx * 2.0f * PI / (nThreads * nBlocks) ); CUDA_SAFE_CALL( cudaMemcpy ( pdA, phA, nMemSizeInBytes, cudaMemcpyHostToDevice ) ); CUDA_SAFE_CALL( cudaBindTexture(0, g_TexRef, pdA, nMemSizeInBytes) ); dim3 threads = dim3( nThreads ); dim3 blocks = dim3( nBlocks ); kernel1 >> ( pdB ); CUDA_SAFE_CALL( cudaThreadSynchronize() ); CUDA_SAFE_CALL( cudaMemcpy ( phB, pdB, nMemSizeInBytes, cudaMemcpyDeviceToHost ) ); return 0; } ! ! ! ! ! !

2011 Texture в CUDA (cudaArray) texture g_TexRef; __global__ void kernel ( float * data ) { int idx = blockIdx.x * blockDim.x + threadIdx.x; data [idx + blockIdx.y * gridDim.x * blockDim.x] = tex2D(g_TexRef, idx, blockIdx.y); } int main ( int argc, char * argv [] ) { float *phA = NULL, *phB = NULL, *pdA = NULL, *pdB = NULL; // linear memory pointers cudaArray * paA = NULL; // device cudaArray pointer cudaChannelFormatDesc cfDesc = cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat); CUDA_SAFE_CALL( cudaMallocArray(&paA, &cfDesc, nBlocksX * nThreads, nBlocksY) ); for (int idx = 0; idx < nThreads * nBlocksX; idx++) { phA[idx] = sinf(idx * 2.0f * PI / (nThreads * nBlocksX) ); phA[idx + nThreads * nBlocksX] = cosf(idx * 2.0f * PI / (nThreads * nBlocksX) ); } CUDA_SAFE_CALL( cudaMemcpyToArray ( paA, 0, 0, phA, nMemSizeInBytes, cudaMemcpyHostToDevice ) ); CUDA_SAFE_CALL( cudaBindTextureToArray(g_TexRef, paA) ); dim3 threads = dim3( nThreads ); dim3 blocks = dim3( nBlocksX, nBlocksY ); kernel2 >> ( pdB ); CUDA_SAFE_CALL( cudaThreadSynchronize() ); CUDA_SAFE_CALL( cudaMemcpy ( phB, pdB, nMemSizeInBytes, cudaMemcpyDeviceToHost ) ); return 0; } ! ! ! ! ! ! !

СВЕРТКА

2011 Свертка В DSP свертка - это один из основных инструментов Определение свертки: В Дискретном случае: В 2D для изображений:

2011 Свертка Ядро Исходный сигнал Окно Выходной сигнал 31

2011 Свертка Вычислительная сложность: – W x H x N x K – умножений Сепарабельные фильтры Размер входного сигнала Размер ядра Ядро 01 Ядро X Ядро Y

2011 Примеры Gaussian Blur Edge Detection

2011 Gaussian Blur Свертка с ядром:

2011 Gaussian Blur #define SQR(x) ((x) * (x)) texture g_TexRef; __global__ void GaussBlur( float * pFilteredImage, int W, int H, float r) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int idy = blockIdx.y * blockDim.y + threadIdx.y; float wSum = 0.0f; float rResult = 0.0f; for (int ix = -r; ix

2011 Свертка Оптимизации Использовать сепарабельные фильтры – Существенно меньше алгоритмическая сложность Использовать shared память

2011 Свертка Оптимизации Исходное изображение

2011 Эффективно ли такое разбиение изображения

2011 Свертка Оптимизации Исходное изображение

2011 Свертка Smem Оптимизации

2011 Свертка Smem Оптимизации

2011 Свертка Smem Оптимизации

2011 Свертка Smem Оптимизации

2011 Свертка Smem Оптимизации

2011 Свертка Smem Оптимизации

2011 Свертка Smem Оптимизации

EDGE DETECTION

2011 Edge Detection Обнаружение границ – поиск разрывов в яркости изображения Идеальная граница «реальная» граница «шумная» граница

2011 Edge Detection Градиент функции f(x,y) – Это вектор который показывает направление роста – Определяется как G GyGy GxGx θ(x,y) = tan -1 (G y /G x )

2011 Edge Detection Разностная производная: Свертка с ядром:

2011 Edge Detection Разностная производная: Свертка с ядром:

2011 Edge Detection Prewitt mask:

2011 Edge Detection Sobel mask:

2011 Edge Detection Оператор Лапласа:

ШУМОПОДАВЛЕНИЕ

2011 Преобразование Фурье Линейный оператор вида: Обратный оператор:

2011 Преобразование Фурье zУсловие существования 1. 2.Конечное число устранимых разрывов

2011 Преобразование Фурье Пример 1D f(x) x W 0 A

2011 Преобразование Фурье Пример 2D f(x,y) W2W2 A W1W1 x y

2011 Преобразование Фурье Свойства

2011 Преобразование Фурье Свойства

2011 Преобразование Фурье Свойства

2011 Фильтры a)Низкочастотные (low-pass) b)Высокочастотные (high-pass) (b)

2011 Фильтры (b)

2011 Discrete Cosine Transform Широко используется в ЦОС Является основой современных алгоритмов сжатия данных с потерями (JPEG, MPEG) JPEG, 2/10

2011 Discrete Cosine Transform Представитель семейства пространственно-частотных 1D преобразований, задается формулами: Прямое: Обратное: Нормировочные коэффициенты:

2011 Discrete Cosine Transform 8-точечный случай: u=0

2011 Discrete Cosine Transform 8-точечный случай: u=1

2011 Discrete Cosine Transform 8-точечный случай: u=2

2011 Discrete Cosine Transform 8-точечный случай: u=3

2011 Discrete Cosine Transform 8-точечный случай: u=4

2011 Discrete Cosine Transform 8-точечный случай: u=5

2011 Discrete Cosine Transform 8-точечный случай: u=6

2011 Discrete Cosine Transform 8-точечный случай: u=7

2011 Discrete Cosine Transform N-мерное преобразование обладает свойством сепарабельности z2D-визуализация коэффициентов для случая 8х8 (изображение справа) zКоэффициенты A[8x8] преобразования вычисляются один раз z

2011 Discrete Cosine Transform Наивный: 64 нити на блок (8х8) – Загрузка одного пикселя из текстуры – Барьер – Поток вычисляет один коэффициент – Барьер – Запись коэффициента в глобальную память

2011 Насколько это эффективно?

2011 Discrete Cosine Transform Блок потоков обрабатывает несколько блоков 8х8 Один поток обрабатывает вектор 8х1 (1x8)

2011 Image Denoising Шумы в изображении – Импульсный Salt & pepper – Аддитивный Uniform Gaussian

2011 Ранговые фильтры Алгоритм Р.Ф. ранга N: – Для каждого отсчета сигнала i – Выбор окрестности вокруг отсчета i Сортируем по значению Выбираем N-ое значение как результат

2011 Медиана Ранговый фильтр N=0.5 Сортировка не обязательна для 8bit значений – Строим гистограмму for(int i

2011 Медиана Построение гистограммы сигнал Сигнал + фартук в Smem

2011 Что с банк-конфликтами?

2011 Фильтрация (Аддитивный шум) Размытие – это low-pass фильтр Каким должен быть фильтр? – Подавлять шум? – Сохранять детальность? noise E(noise) = 0

2011 Gaussian Blur Blur (размытие) изображение Свертка с ядром:

2011 Gaussian Blur Blur (размытие) изображение Свертка с ядром:

2011 Адаптивное размытие Свертка с ядром: ClrSpaceDist – это фотометрическая близость

2011 Bilateral

2011 Bilateral

2011 Bilateral Kernel #define SQR(x) ((x) * (x)) texture g_TexRef; __global__ void BilateralBlur( float * pFilteredImage, int W, int H, float r) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int idy = blockIdx.y * blockDim.y + threadIdx.y; float wSum = 0.0f; float rResult = 0.0f; float c = tex2D(g_TexRef, idx, idy); for (int ix = -r; ix

2011 Bilateral Оптимизации Bilateral не сепарабельный фильтр – Но можно его разделить Смешивать исходное изображение с фильтрованным – Если в блоке много ненулевых коэф., то с большой вероятностью в этом блоке шум был подавлен успешно – Если в блоке много нулевых коэф., то с большой вероятностью в блоке много деталей (границы, текстура и т.д.)

2011 Non Local Means ClrSpaceDist – оценивать по блокам пикселей

2011 Non Local Means На вычисление одного веса: – N b xN b вычислений, N размер блока На фильтрацую одного пиксела: – N b xN b x RxR, R размер окна

2011 Сравнение

2011 Сравнение Bilateral

2011 Сравнение NLM

МАСШТАБИРОВАНИЕ ИЗОБРАЖЕНИЙ

Артифакты Алиасинг – При увелечении – ступенчатость – При уменьшении - муар Ringing Потеря четкости Субпиксельный сдвиг – Влияет на формальные метрики

2011 Простые методы Билинейная интерполяция P1 P7 P3 P9 P4?P4? P2? P5?P6? P8?

2011 Простые методы Билинейная интерполяция – Сепарабельная – Очень быстрая – Поддерживается в HW Точность фильтрации

2011 Простые методы Бикубическая интерполяция – Сепарабельная – Лучше качество t=0 t=1 t=0.5

2011 Сравнение

2011 Lanczos

2011 Gradient interpolation P1 P7 P3 P9 P7? P2? P5?P6? P8? Dxd = abs(P3 – P5) Dyd = abs(P1 – P9) If (Dxd > Dyd) //граница P1P5P9 P5 = (P1 + P9) * 0.5f; If (Dyd > Dxd) //граница P3P5P7 P5 = (P1 + P9) * 0.5f; If (Dyd ~= Dxd) //граница не определена P5 = (P1 + P3 + P7 + P9) * 0.25f;

NEDI (New Edge-Directed Interpolation)

NEDI α1α1 α2α2 α3α3 α4α4 χ 2i 2j X = F(2i+1,2j+1) = α 1 *F(2i,2j)+ α 2 *F(2i+2,2j)+ α 3 *F(2i,2j+2)+ α 4 *F(2i+2,2j+2); α = {α 1, α 2, α 3, α 4 } ?

2011 NEDI α1α1 P 2,2 α4α4 α2α2 P 4,2 α3α3 P 2,4 P 6,2 P 8,2 P 6,4 P 8,4 P 2,6 P 4,8 P 4,6 P 2,8 P 6,6 P 8,6 P 6,8 P 8,8 X i,j = α 1 *F(i-2,j-2) + α 2 *F(i+2,j-2) + α 3 *F(i-2,j+2) + α 4 *F(i+2,j+2); For i = 2,4,6,8 | For j = 2,4,6,8 E i,j = P j,j – X i,j – Approximation Error SSE = ΣΣ sqr( E i,j ); α = Arg min(SSE); α P = {P ij } Cα = P (SSE) / α = 0.0 α = (C T C) -1 C T P …

2011 NEDI (improvement) …

2011 NEDI (improvement) …

2011 NEDI (improvement) …

2011 NEDI: Pros and Cons Pros: NEDI – Четкие тонкие края Cons: Очень медленно на CPU – Умножение матриц 4х16 – Обращение матрицы – Рингинг CUDA: – Большой объем данных на тред – Много регистров – Сложно использовать Smem – Много ветвлений

2011 Сравнение

2011 Сравнение

2011 Вопросы