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

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

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

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

4 Типы памяти в 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 Легенда: -интерфейсы доступа Легенда: -интерфейсы доступа

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

6 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

7 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

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

9 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]

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

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

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

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

14 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 возможна запись из ядра

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

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

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

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

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

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

21 СВЕРТКА

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

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

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

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

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

27 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

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

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

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

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

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

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

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

35 EDGE DETECTION

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

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

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

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

40 Edge Detection Prewitt mask:

41 Edge Detection Sobel mask:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

67 Bilateral

68

69 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

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

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

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

73 Сравнение

74 Сравнение Bilateral

75 Сравнение NLM

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

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

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

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

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

81 Сравнение

82 Lanczos

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

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

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

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

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

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

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

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

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

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

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

94 Сравнение

95

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

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

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

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

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

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

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

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

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

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

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

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

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

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

110 Артифакты

111 Артефакты

112 Вопросы

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

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

115 Ресурсы нашего курса 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

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

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

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

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

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

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

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

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

124 Фильтры (b)