[C][pthreads] Использование преобразований Хафа (Hough transform) для выделения прямых на изображени
Об алгоритме преобразования Хафа можно найти много источников в интернете, поэтому касаться алгоритма я не буду, а сразу перейду к программной реализации метода и демонстрации примера.
Строить массив [R, Θ] можно двумя способами. Первый - последовательно пройтись по всем точкам этого массива, занося туда суммарную интенсивность пикселей, расположенных на соответствующей прямой. Сложность этого способа пропорциональна произведению площади изображения-преобразования Хафа на характерный размер оригинального изображения.
Второй способ - последовательно пройтись по всем точкам изображения, инкрементируя на величину интенсивности в данном пикселе все ячейки преобразования Хафа, соответствующие возможным прямым, проходящим через данный пиксель. Сложность этого способа пропорциональна произведению площади изображения на характерный размер изображения-преобразования Хафа.
Какой из способов использовать - выбор наиболее оптимальной производительности. В моем случае второй способ был удобнее (т.к. изображение-преобразование Хафа имеет меньшую площадь).
Еще один вопрос: что заносить в ячейки преобразования Хафа. Мы можем заносить туда единицу для каждого пикселя, чья интенсивность выше некоего порогового значения, а можем заносить саму эту интенсивность. Я выбрал второй способ, т.к. он лучше отражает наиболее интенсивные линии на фоне более слабых.
Так как это преобразование - задача довольно ресурсоемкая, ее можно попытаться распараллелить. Единственным узким местом является точка занесения данных в результирующий массив. Поэтому желательно либо синхронизовать эту точку (при помощи мьютексов, например), либо же использовать атомарные операции.
Для начала рассмотрим реализацию преобразования Хафа на CPU, распараллеливание будем осуществлять потоками (при помощи библиотеки pthreads). Можно использовать и openMP, но данная реализация ближе к следующему шагу - переносу алгоритма на GPU.
Преобразование Хафа будет иметь размер Rmax x angles и дискретизовано по R с шагом 1, от нуля до Rmax. По Θ оно дискретизуется с шагом angles/270 (начало координат для изображения мы считаем лежащим в нижнем левом углу, поэтому углы [-180°, -90°) не будут содержать никакой полезной информации).
Как изображение, так и преобразование Хафа нормированы на единицу (т.е. данные в них лежат в диапазоне [0,1]).
Приведенная ниже функция и является основной: инициализирует необходимые данные и запускает потоки.
Аргументы функции:
int fill_hough_lines(float *ima, float min, float max, int imW, int imH, int Rmax, int angles, float *hough){
int i, Y0, Y1, dY;
float *hptr = hough, hmax;
Hough_kernel_args HD[THREAD_NUMBER];
pthread_t threads[THREAD_NUMBER];
pthread_mutex_t mutex;
pthread_mutex_init(&mutex, NULL);
init_sincos(angles);
dY = (imH + THREAD_NUMBER/2) / THREAD_NUMBER;
Y0 = 0; Y1 = dY;
for(i = 0; i < THREAD_NUMBER; i++, Y0+=dY, Y1+=dY){
if(Y0 >= imW) break;
if(Y1 > imW) Y1 = imW;
HD[i].min = min; HD[i].max = max;
HD[i].image = ima; HD[i].hough = hough;
HD[i].Y0 = Y0; HD[i].Y1 = Y1; HD[i].w = imW;
HD[i].Rmax = Rmax; HD[i].angles = angles;
HD[i].mutex = &mutex;
pthread_create(&threads[i], NULL, fill_lin_hough_array, &HD[i]);
}
for(i = 0; i < THREAD_NUMBER; i++)
pthread_join(threads[i], NULL);
pthread_mutex_destroy(&mutex);
hmax = *hough++;
Y1 = Rmax * angles;
for(i = 1; i < Y1; i++, hough++)
if(*hough > hmax) hmax = *hough;
hough = hptr;
for(i = 0; i < Y1; i++, hough++)
*hough /= hmax;
return 1;
}
В конце мы нормируем полученное преобразование Хафа (т.к. операция эта достаточно простая, я ее не распараллеливал - кстати, при распараллеливании она выполняется дольше).
Структура Hough_kernel_args - вспомогательная, чтобы не передавать в поток уйму аргументов:
typedef struct{
float *image;
float min;
float max;
float *hough;
int w;
int Y0;
int Y1;
int Rmax;
int angles;
pthread_mutex_t *mutex;
}Hough_kernel_args;
Функция init_sincos подготавливает массив синусов и косинусов для дальнейших вычислений. Так как нам потом придется в каждой итерации для каждого угла вычислять синусы и косинусы, то быстрее будет сразу затабулировать их значения.
static float *Sin_d = NULL, *Cos_d = NULL;
static int sincosize = 0;
void init_sincos(int angles){
float anglestep = (float)angles / 270.;
float conv = M_PI / 180.;
float *cp, *sp;
int k;
if(!Sin_d || !Cos_d || angles != sincosize){
free(Sin_d);
free(Cos_d);
Sin_d = malloc(angles * sizeof(float));
Cos_d = malloc(angles * sizeof(float));
cp = Cos_d; sp = Sin_d;
for(k = 0; k < angles; k++){
float theta = ((float)k / anglestep - 90.) * conv;
*sp++ = sinf(theta); // или sincosf(theta, sp++, cp++);
*cp++ = cosf(theta);
}
sincosize = angles;
}
}
Для хранения массивов и их размера используются глобальные переменные. Инициализация массивов производится лишь в том случае, если они еще не были созданы, либо их размер не совпадает с требуемым. Функцию sincosf я здесь не использовал, т.к. она у меня почему-то не работала.
Ну и наконец ядро преобразования Хафа:
const float ImTreshold = 0.1;
void *fill_lin_hough_array(void *data){
Hough_kernel_args *HD = (Hough_kernel_args*) data;
float *ima = HD->image, *hough = HD->hough;
float min = HD->min, wd = HD->max - min;
pthread_mutex_t *mutex = HD->mutex;
int Y0 = HD->Y0, Y1 = HD->Y1, imW = HD->w, Rmax = HD->Rmax, angles = HD->angles;
int i, j, k, Y, R;
for(j = Y0; j < Y1; j++){
for(i = 0.; i < imW; i++, ima++){
float imdata = (*ima - min) / wd;
if(imdata > ImTreshold){
Y = 0;
pthread_mutex_lock(mutex);
for(k = 0; k < angles; k++, Y+=Rmax){
R = (int)(0.5 + i*Cos_d[k] + j*Sin_d[k]);
if(R > 0 && R < Rmax) hough[R + Y] += imdata;
}
pthread_mutex_unlock(mutex);
}
}
}
return NULL;
}
Здесь каждый поток просматривает строки изображения в диапазоне [Y0, Y1) и заполняет итоговый массив, пробегаясь по всем углам и вычисляя соответствующие радиусы.
Работает этот алгоритм следующим образом. Пусть у нас есть вот такое изображение:
и мы хотим распознать на нем характерные линии. Для этого построим преобразование Хафа. В качестве максимального значения по R выберем длину диагонали изображения, а углы будем вычислять с шагом в полградуса (для первого приближения этого достаточно - далее можно сузить область поиска и уменьшить шаг как по R, так и по Θ).
Итак, преобразование Хафа для этого изображения выглядит следующим образом:
Чтобы лучше разглядеть относительные интенсивности, покажу еще в 3D:
Если теперь мы выберем центры пятен с максимальной интенсивностью на преобразовании Хафа, получим характерные линии:
(здесь я не вычислял центры точно, а лишь примерно пощелкал мышкой в характерных точках).
В следующий раз я приведу пример преобразования Хафа с использованием CUDA.
Строить массив [R, Θ] можно двумя способами. Первый - последовательно пройтись по всем точкам этого массива, занося туда суммарную интенсивность пикселей, расположенных на соответствующей прямой. Сложность этого способа пропорциональна произведению площади изображения-преобразования Хафа на характерный размер оригинального изображения.
Второй способ - последовательно пройтись по всем точкам изображения, инкрементируя на величину интенсивности в данном пикселе все ячейки преобразования Хафа, соответствующие возможным прямым, проходящим через данный пиксель. Сложность этого способа пропорциональна произведению площади изображения на характерный размер изображения-преобразования Хафа.
Какой из способов использовать - выбор наиболее оптимальной производительности. В моем случае второй способ был удобнее (т.к. изображение-преобразование Хафа имеет меньшую площадь).
Еще один вопрос: что заносить в ячейки преобразования Хафа. Мы можем заносить туда единицу для каждого пикселя, чья интенсивность выше некоего порогового значения, а можем заносить саму эту интенсивность. Я выбрал второй способ, т.к. он лучше отражает наиболее интенсивные линии на фоне более слабых.
Так как это преобразование - задача довольно ресурсоемкая, ее можно попытаться распараллелить. Единственным узким местом является точка занесения данных в результирующий массив. Поэтому желательно либо синхронизовать эту точку (при помощи мьютексов, например), либо же использовать атомарные операции.
Для начала рассмотрим реализацию преобразования Хафа на CPU, распараллеливание будем осуществлять потоками (при помощи библиотеки pthreads). Можно использовать и openMP, но данная реализация ближе к следующему шагу - переносу алгоритма на GPU.
Преобразование Хафа будет иметь размер Rmax x angles и дискретизовано по R с шагом 1, от нуля до Rmax. По Θ оно дискретизуется с шагом angles/270 (начало координат для изображения мы считаем лежащим в нижнем левом углу, поэтому углы [-180°, -90°) не будут содержать никакой полезной информации).
Как изображение, так и преобразование Хафа нормированы на единицу (т.е. данные в них лежат в диапазоне [0,1]).
Приведенная ниже функция и является основной: инициализирует необходимые данные и запускает потоки.
Аргументы функции:
- ima - входное изображение
- min, max - экстремальные значения интенсивности по изображению
- imW, imH - размеры изображения (ширина и высота)
- Rmax, angles - размеры выходного преобразования Хафа (ширина равна Rmax, высота равна angles)
- hough - изображение, в которое и будет заноситься наше преобразование Хафа
int fill_hough_lines(float *ima, float min, float max, int imW, int imH, int Rmax, int angles, float *hough){
int i, Y0, Y1, dY;
float *hptr = hough, hmax;
Hough_kernel_args HD[THREAD_NUMBER];
pthread_t threads[THREAD_NUMBER];
pthread_mutex_t mutex;
pthread_mutex_init(&mutex, NULL);
init_sincos(angles);
dY = (imH + THREAD_NUMBER/2) / THREAD_NUMBER;
Y0 = 0; Y1 = dY;
for(i = 0; i < THREAD_NUMBER; i++, Y0+=dY, Y1+=dY){
if(Y0 >= imW) break;
if(Y1 > imW) Y1 = imW;
HD[i].min = min; HD[i].max = max;
HD[i].image = ima; HD[i].hough = hough;
HD[i].Y0 = Y0; HD[i].Y1 = Y1; HD[i].w = imW;
HD[i].Rmax = Rmax; HD[i].angles = angles;
HD[i].mutex = &mutex;
pthread_create(&threads[i], NULL, fill_lin_hough_array, &HD[i]);
}
for(i = 0; i < THREAD_NUMBER; i++)
pthread_join(threads[i], NULL);
pthread_mutex_destroy(&mutex);
hmax = *hough++;
Y1 = Rmax * angles;
for(i = 1; i < Y1; i++, hough++)
if(*hough > hmax) hmax = *hough;
hough = hptr;
for(i = 0; i < Y1; i++, hough++)
*hough /= hmax;
return 1;
}
В конце мы нормируем полученное преобразование Хафа (т.к. операция эта достаточно простая, я ее не распараллеливал - кстати, при распараллеливании она выполняется дольше).
Структура Hough_kernel_args - вспомогательная, чтобы не передавать в поток уйму аргументов:
typedef struct{
float *image;
float min;
float max;
float *hough;
int w;
int Y0;
int Y1;
int Rmax;
int angles;
pthread_mutex_t *mutex;
}Hough_kernel_args;
Функция init_sincos подготавливает массив синусов и косинусов для дальнейших вычислений. Так как нам потом придется в каждой итерации для каждого угла вычислять синусы и косинусы, то быстрее будет сразу затабулировать их значения.
static float *Sin_d = NULL, *Cos_d = NULL;
static int sincosize = 0;
void init_sincos(int angles){
float anglestep = (float)angles / 270.;
float conv = M_PI / 180.;
float *cp, *sp;
int k;
if(!Sin_d || !Cos_d || angles != sincosize){
free(Sin_d);
free(Cos_d);
Sin_d = malloc(angles * sizeof(float));
Cos_d = malloc(angles * sizeof(float));
cp = Cos_d; sp = Sin_d;
for(k = 0; k < angles; k++){
float theta = ((float)k / anglestep - 90.) * conv;
*sp++ = sinf(theta); // или sincosf(theta, sp++, cp++);
*cp++ = cosf(theta);
}
sincosize = angles;
}
}
Для хранения массивов и их размера используются глобальные переменные. Инициализация массивов производится лишь в том случае, если они еще не были созданы, либо их размер не совпадает с требуемым. Функцию sincosf я здесь не использовал, т.к. она у меня почему-то не работала.
Ну и наконец ядро преобразования Хафа:
const float ImTreshold = 0.1;
void *fill_lin_hough_array(void *data){
Hough_kernel_args *HD = (Hough_kernel_args*) data;
float *ima = HD->image, *hough = HD->hough;
float min = HD->min, wd = HD->max - min;
pthread_mutex_t *mutex = HD->mutex;
int Y0 = HD->Y0, Y1 = HD->Y1, imW = HD->w, Rmax = HD->Rmax, angles = HD->angles;
int i, j, k, Y, R;
for(j = Y0; j < Y1; j++){
for(i = 0.; i < imW; i++, ima++){
float imdata = (*ima - min) / wd;
if(imdata > ImTreshold){
Y = 0;
pthread_mutex_lock(mutex);
for(k = 0; k < angles; k++, Y+=Rmax){
R = (int)(0.5 + i*Cos_d[k] + j*Sin_d[k]);
if(R > 0 && R < Rmax) hough[R + Y] += imdata;
}
pthread_mutex_unlock(mutex);
}
}
}
return NULL;
}
Здесь каждый поток просматривает строки изображения в диапазоне [Y0, Y1) и заполняет итоговый массив, пробегаясь по всем углам и вычисляя соответствующие радиусы.
Работает этот алгоритм следующим образом. Пусть у нас есть вот такое изображение:
и мы хотим распознать на нем характерные линии. Для этого построим преобразование Хафа. В качестве максимального значения по R выберем длину диагонали изображения, а углы будем вычислять с шагом в полградуса (для первого приближения этого достаточно - далее можно сузить область поиска и уменьшить шаг как по R, так и по Θ).
Итак, преобразование Хафа для этого изображения выглядит следующим образом:
Чтобы лучше разглядеть относительные интенсивности, покажу еще в 3D:
Если теперь мы выберем центры пятен с максимальной интенсивностью на преобразовании Хафа, получим характерные линии:
(здесь я не вычислял центры точно, а лишь примерно пощелкал мышкой в характерных точках).
В следующий раз я приведу пример преобразования Хафа с использованием CUDA.