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

TEXTURE

Типы памяти в 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

Архитектура 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 zВ 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 zЛатентность больше, чем у прямого обращения в память yДополнительные стадии в конвеере: xПреобразование адресов xФильтрация xПреобразование данных zНо зато есть кэш yРазумно использовать, если: xОбъем данных не влезает в shared память xПаттерн доступа хаотичный xДанные переиспользуются разными потоками

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

Texture HW zПреобразование координат: yКоординаты, которые не лежат в диапазоне [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]

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

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

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

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

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

Texture в CUDA (cudaArray) zОсобенности текстур: yОбращение к 1D / 2D / 3D массивам данных по: xЦелочисленным индексам xНормализованным координатам yПреобразование адресов на границах xClamp xWrap yФильтрация данных xPoint xLinear yПреобразование данных Данные могут храниться в формате uchar4 Возвращаемое значение – float4

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

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

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

СВЕРТКА

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

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

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

Примеры zGaussian Blur zEdge Detection

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

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

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

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

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

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

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

EDGE DETECTION

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

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

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

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

Edge Detection zPrewitt mask:

Edge Detection zSobel mask:

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

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

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

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

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

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

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

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

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

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

Фильтры (b)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Image Denoising zШумы в изображении yИмпульсный xSalt & pepper yАддитивный xUniform xGaussian

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

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

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

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

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

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

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

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

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 Оптимизации zBilateral не сепарабельный фильтр yНо можно его разделить zСмешивать исходное изображение с фильтрованным yЕсли в блоке много ненулевых коэф., то с большой вероятностью в этом блоке шум был подавлен успешно yЕсли в блоке много нулевых коэф., то с большой вероятностью в блоке много деталей (границы, текстура и т.д.)

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

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

Сравнение

Сравнение Bilateral

Сравнение NLM

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

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

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

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

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

Сравнение

Lanczos

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

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 …

NEDI (improvement) …

NEDI (improvement) …

NEDI (improvement) …

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

Сравнение

Вопросы

Assignements zПервое и второе задание доступно zСрок сдачи первого –23 Марта zСрок сдачи второго – 30 Марта

ASSIGNMENT 1

Assignment 1.1 zМетод наименьших квадратов yПрограмма принимает размер матрицы в формате MxN – (строки x столбцы) xM > N yПрограмма генерирует случайные числа для матрицы C и для вектора P yПишет в текстовый файл результат xC; (C T C); (C T C) -1 ; α = (C T C) -1 C T P;

Assignment 1.2 z Дан отрезок [a, b] z Ген N отрезков [ai, bi] i={0, N} zДан размер гистограммы nBin yS: сколько отрезков началось слева yE: сколько отрезков закончилось справа a b S:E nBin

Assignment 1.2 z Найти такую точку qi что выполнено оптимальное соотношение: z(pi - ai) * Bin qi () a b S:E nBin

ASSIGNMENT 2

Assignment 2.1 zДан параметр R – относительный радиус ядра (R (0, 1)) zСвертка изображения с использованием FFT zЯдра для свертки – идеальный high pass и low pass

Assignment 2.2 zДан файл test.bmp (RGB) zПодавление высокочастотного шума в исходном изображении. zВычисляются частные производные, используя один из известных шаблонов zВычисляются длина и угол наклона градиента zВ каждой точке по направлению угла делается несколько выборок и определеятся является ли яркость данной точки локальным максимумом функции яркости. yДА : значение пиксела приравнивается 1.0f yНЕТ : значение пиксела приравнивается 0.0f. zПорог, по которому делается дополнительное отсеивание по яркости.

Assignment 2.2

Общие правила по оформлению прорамм zПрограмма должна делать проверки на ошибки: zНаличие девайса? zОткрылся ли нужный файл? zПравильного ли он формата? zПрограмма должна быть скомпилирована в Release и запускаться на Windows XP SP2 с CUDA Toolkit 3.0 zПрограмма должна компилироваться zДля этого должен быть приложен vcproj для VS2005 либо (makefile +.bat)

Общие правила по оформлению программ zЕсли вы используете любые другие инклюды кроме стандартных – не расчитывайте, что они прописаны на проверяющей машине. Пример того, чего не будет на машине: ycutil.h (требует установки CUDA SDK) zПример того, что будет на машине: ycudart.h (ставиться вместе с CUDA toolkit) ystdio.h (стандартная C библиотека)

Вопросы