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

План Работа с Текстурами Свертка Шумоподавление Масштабирование изображений Сжатие изображений

РАБОТА С ТЕКСТУРАМИ

Типы памяти в CUDA Тип памяти ДоступУровень выделения Скорость работы Расположение Регистры R/WPer-threadВысокаяSM Локальная R/WPer-threadНизкаяDRAM Shared R/WPer-blockВысокаяSM Глобальная R/WPer-gridНизкаяDRAM Constant R/OPer-gridВысокаяDRAM Texture R/OPer-gridВысокаяDRAM Легенда: -интерфейсы доступа Легенда: -интерфейсы доступа

Архитектура 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

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

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

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

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

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

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

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

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

cudaArray Особый контейнер памяти Черный ящик для приложения Позволяет организовывать данные в 1D/ 2D/3D массивы данных вида: – 1/2/4 компонентные векторы – 8/16/32 bit signed/unsigned integers – 32 bit float – 16 bit float (driver API) На CC 2.x возможна запись из ядра

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

Surface в CUDA (cudaArray) Чтение и запись – cudaArraySurfaceLoadStore Обращение к 1D/2D массивам данных по: – Целочисленным индексам – Byte-адресация по x Преобразование адресов на границах – Clamp, Zero, Trap

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

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; }

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 cfD=cudaCreateChannelDesc(32,0,0,0, cudaChannelFormatKindFloat); CUDA_SAFE_CALL( cudaMallocArray(&paA, &cfD, 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; }

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

СВЕРТКА

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

i Свертка Ядро Исходный сигнал Окно Выходной сигнал Легенда: -сигнал -ядро свертки Легенда: -сигнал -ядро свертки j

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

Примеры Gaussian Blur Edge Detection

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

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

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

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

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

Свертка Smem Оптимизации Легенда: -сигнал -ядро свертки Легенда: -сигнал -ядро свертки

Свертка Smem Оптимизации Легенда: -сигнал -ядро свертки Легенда: -сигнал -ядро свертки

Свертка Smem Оптимизации Легенда: -сигнал -ядро свертки Легенда: -сигнал -ядро свертки

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

EDGE DETECTION

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

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

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

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

Edge Detection Prewitt mask:

Edge Detection Sobel mask:

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

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

Discrete Cosine Transform Является основой современных алгоритмов сжатия данных с потерями (JPEG, MPEG) JPEG, 2/10

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

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

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

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

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

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

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

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

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

Discrete Cosine Transform N-мерное преобразование обладает свойством сепарабельности Коэффициенты A[8x8] преобразования вычисляются один раз

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

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

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

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

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

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

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

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

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

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

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

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

Bilateral

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

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

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

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

Сравнение

Сравнение Bilateral

Сравнение NLM

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

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

Простые методы Билинейная интерполяция P1 P7 P3 P9 P4P4 P2 P5P6 P8 Легенда: -исходные пиксели -искомые пиксели Легенда: -исходные пиксели -искомые пиксели

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

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

Сравнение

Lanczos

Gradient interpolation P1 P7 P3 P9 P7 P2 P5P6 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 (New Edge-Directed Interpolation) Легенда: -исходные пиксели -искомые пиксели Легенда: -исходные пиксели -искомые пиксели

NEDI (New Edge-Directed Interpolation) Легенда: -исходные пиксели -искомые пиксели Легенда: -исходные пиксели -искомые пиксели

NEDI (New Edge-Directed Interpolation) Легенда: -исходные пиксели -искомые пиксели Легенда: -исходные пиксели -искомые пиксели

NEDI α1α1 α3α3 α2α2 α4α4 x 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 } ? Легенда: -исходные пиксели -искомые пиксели Легенда: -исходные пиксели -искомые пиксели

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, … Легенда: -исходные пиксели -искомые пиксели Легенда: -исходные пиксели -искомые пиксели

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, … 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 Легенда: -исходные пиксели -искомые пиксели Легенда: -исходные пиксели -искомые пиксели

NEDI (improvement) … … Легенда: -исходные пиксели -искомые пиксели Легенда: -исходные пиксели -искомые пиксели

NEDI (improvement) … … Легенда: -исходные пиксели -искомые пиксели - вспомогательные пиксели Легенда: -исходные пиксели -искомые пиксели - вспомогательные пиксели

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

Сравнение

ФРАКТАЛЬНОЕ СЖАТИЕ АЛГОРИТМА

Мат.часть Сжимающее отображение f, X X: Если f сжимающее, то

Фрактальное сжатие РангДомен

Фрактальное сжатие РангДомен

Фрактальное сжатие РангДомен

Фрактальное сжатие РангДомен

... Фрактальное сжатие РангДомен

Фрактальное сжатие РангДомен

Фрактальное сжатие РангДомен

Фрактальное сжатие

Блок пикселей заменяется преобразованием – Сдвиг, поворот / отражение – s,t – преобразование яркости float2 x rgb s,t – сжимают во время компресси – Например s:5,t:3 – в интервале [0, 1]

Декомпрессия Сжимающее отображение f, X X: Если f сжимающее, то Берем любое изображение Строим Доменное изображение Применим полученный набор преобразований к блокам

Сколько раз необходимо применять преобразования?

Декомпрессия

Артифакты

Артефакты

Вопросы

Задание 1 Вычислить число pi методом монтекарло – Генерация точки (x, y) в [-1, 1] x [-1, 1] – Расчет сколько точек попало в ед. круг Построение Kd дерева для точек – Задние на thrust

Задание 2 Изменение времени звуковой дорожки без изменения тона – Задание на CUFFT Обучение фильтра – Разбить изображение на блоки – Классифицировать каждый блок – Обучить фильтр имея пару (R, T) R – «правильное» изображение T – «входной» сигнал

Ресурсы нашего курса Steps3d.Narod.Ru Google Site CUDA.CS.MSU.SU Google Group CUDA.CS.MSU.SU Google Mail CS.MSU.SU Google SVN Tesla.Parallel.Ru Twirpx.Com Nvidia.Ru

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

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

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

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

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

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

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

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

Фильтры (b)