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

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

2 TEXTURE

3 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

4 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

5 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

6 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

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

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

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

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

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

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

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

14 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()

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

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

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

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

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

20 СВЕРТКА

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

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

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

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

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

26 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

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

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

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

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

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

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

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

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

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

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

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

38 EDGE DETECTION

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

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

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

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

43 2011 Edge Detection Prewitt mask:

44 2011 Edge Detection Sobel mask:

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

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

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

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

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

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

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

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

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

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

55 2011 Фильтры (b)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

79 2011 Bilateral

80 2011 Bilateral

81 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

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

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

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

85 2011 Сравнение

86 2011 Сравнение Bilateral

87 2011 Сравнение NLM

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

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

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

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

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

93 2011 Сравнение

94 2011 Lanczos

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

96 NEDI (New Edge-Directed Interpolation)

97 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 } ?

98 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 …

99 2011 NEDI (improvement) …

100 2011 NEDI (improvement) …

101 2011 NEDI (improvement) …

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

103 2011 Сравнение

104 2011 Сравнение

105 2011 Вопросы