Емельянов Эдуард Владимирович (eddy_em) wrote,
Емельянов Эдуард Владимирович
eddy_em

Category:

[C][pthreads] Использование преобразований Хафа (Hough transform) для выделения прямых на изображени

Об алгоритме преобразования Хафа можно найти много источников в интернете, поэтому касаться алгоритма я не буду, а сразу перейду к программной реализации метода и демонстрации примера.

Строить массив [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 — изображение, в которое и будет заноситься наше преобразование Хафа
THREAD_NUMBER — макрос, определяемый на момент компиляции (это — количество ядер вашего процессора), именно столько потоков мы и будем создавать.
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.





Tags: c, hough transform, велосипедостроение
Subscribe

  • Дурацкий перекресток

    Был на днях в Пятигорске. Ну и движение там! Просто жесть!!! Вечные пробки, куча "кругов" и грохотящие трамваи… А когда выезжал оттуда, на углу пр.…

  • А что, в С так нельзя?

    Пытаюсь передать в функцию цвет как массив. Функция такая: void Pattern_draw3(Img3 *img, Pattern *p, int xul, int yul, uint8_t colr[3]); И…

  • Документация! Никто не любит документацию…

    Хочу использовать занятную header-only библиотеку STB для чтения-записи изображений (останется лишь фитсы добавить), а также добавления к ним…

promo eddy_em август 17, 2019 12:33 3
Buy for 10 tokens
Юра намедни напечатал корпус для хронометра. Для первого блина получилось неплохо: И еще немного фотографий:
  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your reply will be screened

  • 2 comments