Морфологические операции: эрозия и диляция.
Итак, как я уже говорил в предыдущей записи, за время, бесцельно проведенное в Нижнем Новгороде, кое-что полезное я таки сделал. В этой записи расскажу о реализации операций эрозии и диляции.
Бла-бла
Точные определения эрозии (сжатия) и диляции (расширения) бинарных изображений можно прочитать в википедии. Вкратце: эти операции чем-то похожи на свертку изображения с некоей маской, основная же разница от свертки в том, что в случае эрозии сложение заменяется логической операцией "И", а в случае диляции - операцией "ИЛИ". Т.е. эрозия оставляет на изображении лишь те пиксели, которые имеют абсолютно всех соседей по аналогии с маской, а диляция же добавляет каждому пикселю соседей из маски. Понятное дело, операции эти - совершенно необратимые.
Кстати, больший интерес представляют даже не сами по себе эрозия и диляция, а их комбинации: закрытие (эрозия диляции) и открытие (диляция эрозии). Эти операции наглядно показывают необратимость отдельных операций: эрозия диляции позволяет получить вроде бы первоначальный объект, но с закрытыми "дырками"; диляция же эрозии дает вроде бы первоначальный объект, но с удалением небольших выбросов и подчисткой индивидуальных пикселей (и даже небольших объектов, если наша маска достаточно большая).
Маски могут иметь самые разнообразные формы. Мне же интересно было реализовать эрозию и диляцию наиболее часто употребляемой маской - крестообразной (т.е. для 4-связанных пикселей). Т.к. интернета у меня не было, исходников лептоники я скачать не мог, чтобы подсмотреть, как же там эти операции реализованы. Я помню лишь, что в лептонике бинарные изображения хранились в "упакованном" виде: по 8 пикселей на байт. Я подумал, что это - действительно разумный подход, т.к. многие операции можно сделать значительно быстрей. На первый взгляд более логичным кажется хранить пиксели в "идеальном" виде: по 64 пикселя на машинное слово, однако, дальше я покажу, почему так делать не получится.
Итак, сначала - немного соображений по поводу эрозии и диляции (то, что я набросал, будучи в НН).
Эрозия крестообразной маской, как я уже говорил выше, представляет собой вот что:
00000 00000 01110 00000 01110 => 00100 01110 00000 00000 00000 Бит равен единице лишь тогда, когда присутствуют все его 4 связи. Т.е. при работе с "упакованным" изображением для срединных битов байта выполняем логическое "&" текущего байта с байтами верхней и нижней строк. Для крайних битов необходимо проверить еще и 4-связные байты. Обозначим текущий байт буквой "E", а соседние с ним проименуем так:
ABC DEF GHK Младший бит E проверяется по старшему биту F, младшим битам B и H и второму биту E. Старший бит E проверяется по младшему биту D, старшим битам B и E и седьмому биту E.
Таким образом, алгоритм получается следующим:
Вот эту проверку я было думал сделать, но когда обнаружил, что операция сама по себе так шустро выполняется, что дополнительная оптимизация вряд ли поднимет скорость больше, чем на порядок, понял, что ничего переделывать и не нужно. Я думаю, что если результат тяжелой оптимизации повышает производительность меньше, чем на порядок, такая оптимизация нафиг не нужна! Если же эта оптимизация дается малой кровью (скажем, минут за 5), то почему бы и нет (но этот случай явно не из разряда пятиминутных).
Диляция крестообразной маской:
00000 00110 00110 01111 01010 => 11111 01100 11110 00000 01100 Здесь противоположная ситуация: вместо "&" делаем "|":
Бла-бла
Точные определения эрозии (сжатия) и диляции (расширения) бинарных изображений можно прочитать в википедии. Вкратце: эти операции чем-то похожи на свертку изображения с некоей маской, основная же разница от свертки в том, что в случае эрозии сложение заменяется логической операцией "И", а в случае диляции - операцией "ИЛИ". Т.е. эрозия оставляет на изображении лишь те пиксели, которые имеют абсолютно всех соседей по аналогии с маской, а диляция же добавляет каждому пикселю соседей из маски. Понятное дело, операции эти - совершенно необратимые.
Кстати, больший интерес представляют даже не сами по себе эрозия и диляция, а их комбинации: закрытие (эрозия диляции) и открытие (диляция эрозии). Эти операции наглядно показывают необратимость отдельных операций: эрозия диляции позволяет получить вроде бы первоначальный объект, но с закрытыми "дырками"; диляция же эрозии дает вроде бы первоначальный объект, но с удалением небольших выбросов и подчисткой индивидуальных пикселей (и даже небольших объектов, если наша маска достаточно большая).
Маски могут иметь самые разнообразные формы. Мне же интересно было реализовать эрозию и диляцию наиболее часто употребляемой маской - крестообразной (т.е. для 4-связанных пикселей). Т.к. интернета у меня не было, исходников лептоники я скачать не мог, чтобы подсмотреть, как же там эти операции реализованы. Я помню лишь, что в лептонике бинарные изображения хранились в "упакованном" виде: по 8 пикселей на байт. Я подумал, что это - действительно разумный подход, т.к. многие операции можно сделать значительно быстрей. На первый взгляд более логичным кажется хранить пиксели в "идеальном" виде: по 64 пикселя на машинное слово, однако, дальше я покажу, почему так делать не получится.
Итак, сначала - немного соображений по поводу эрозии и диляции (то, что я набросал, будучи в НН).
Эрозия крестообразной маской, как я уже говорил выше, представляет собой вот что:
00000 00000 01110 00000 01110 => 00100 01110 00000 00000 00000 Бит равен единице лишь тогда, когда присутствуют все его 4 связи. Т.е. при работе с "упакованным" изображением для срединных битов байта выполняем логическое "&" текущего байта с байтами верхней и нижней строк. Для крайних битов необходимо проверить еще и 4-связные байты. Обозначим текущий байт буквой "E", а соседние с ним проименуем так:
ABC DEF GHK Младший бит E проверяется по старшему биту F, младшим битам B и H и второму биту E. Старший бит E проверяется по младшему биту D, старшим битам B и E и седьмому биту E.
Таким образом, алгоритм получается следующим:
- подменяем E маской[E] (эрозия внутренних битов)
- делаем &: E = E & B & H
- корректируем старший бит E по младшему D и младший бит E по старшему F.
Вот эту проверку я было думал сделать, но когда обнаружил, что операция сама по себе так шустро выполняется, что дополнительная оптимизация вряд ли поднимет скорость больше, чем на порядок, понял, что ничего переделывать и не нужно. Я думаю, что если результат тяжелой оптимизации повышает производительность меньше, чем на порядок, такая оптимизация нафиг не нужна! Если же эта оптимизация дается малой кровью (скажем, минут за 5), то почему бы и нет (но этот случай явно не из разряда пятиминутных).
Диляция крестообразной маской:
00000 00110 00110 01111 01010 => 11111 01100 11110 00000 01100 Здесь противоположная ситуация: вместо "&" делаем "|":
- подменяем E маской (диляция внутренних битов),
- делаем |: E = E | B | H,
- корректируем старший и младший биты операцией |. Реализация
Преобразование бинарного изображения в "упакованное" реализуется такой вот элементарщиной:
/** * Convert boolean image into pseudo-packed (1 char == 8 pixels) * @param im (i) - image to convert * @param W, H - size of image im (must be larger than 1) * @param W_0 (o) - (stride) new width of image * @return allocated memory area with "packed" image */ unsigned char *bool2char(bool *im, int W, int H, int *W_0){ if(W < 2 || H < 2) errx(1, "bool2char: image size too small"); int y, W0 = (W + 7) / 8; unsigned char *ret = Malloc(W0 * H, 1); OMP_FOR() for(y = 0; y < H; y++){ int x, i, X; bool *ptr = &im[y*W]; unsigned char *rptr = &ret[y*W0]; for(x = 0, X = 0; x < W0; x++, rptr++){ for(i = 7; i > -1 && X < W; i--, X++, ptr++){ *rptr |= *ptr << i; } } } if(W_0) *W_0 = W0; return ret; } Просто пробегаемся в цикле по оригинальному изображению и пакуем по 8 бит в один байт. Параметр W_0, возвращаемый функцией, по сути является тем самым stride из всяких кудовских функций.
Единственное, что я делаю для ускорения вычислений - так это заполнение матриц a&(a<<1)&(a>>1) и a|(a<<1)|(a>>1), т.к. эти вычисления выполняются ну очень часто. Заполнение происходит при первом же обращении к морфологическим функциям и реализуется так:
/** * This function inits masks arrays for erosion and dilation * You may call it yourself or it will be called when one of * `erosion` or `dilation` functions will be ran first time */ void morph_init(){ if(__Init_done) return; int i; ER = Malloc(256, 1); DIL = Malloc(256, 1); for(i = 0; i < 256; i++){ ER[i] = i & ((i << 1) | 1) & ((i >> 1) | (0x80)); // don't forget that << and >> set borders to zero DIL[i] = i | (i << 1) | (i >> 1); } __Init_done = true; } Здесь небольшая коррекция сделана для матрицы ER: мы не можем просто написать
ER[i] = i & (i << 1) & (i >> 1); так как две операции сдвига сбросят крайние пиксели в i, поэтому-то приходится хитрить.
Чтобы не выполнять на каждом пикселе уймы лишних if'ов, для таких элементарных операций я сделал просто: разбил цикл на кусочки: "внутренность" изображения обрабатывается в цикле, а вот крайние столбцы и строки - по-отдельности.
Операция эрозии - самая простая, т.к. здесь все пиксели на краю изображения просто обнуляются:
/** * Make morphological operation of erosion * @param image (i) - input image * @param W, H - size of image * @return allocated memory area with erosion of input image */ unsigned char *erosion(unsigned char *image, int W, int H){ if(W < 2 || H < 2) errx(1, "Erosion: image size too small"); if(!__Init_done) morph_init(); unsigned char *ret = Malloc(W*H, 1); int y, w = W-1, h = H-1; OMP_FOR() for(y = 1; y < h; y++){ // reset first & last rows of image unsigned char *iptr = &image[W*y]; unsigned char *optr = &ret[W*y]; unsigned char p = ER[*iptr] & 0x7f & iptr[-W] & iptr[W]; int x; if(!(iptr[1] & 0x80)) p &= 0xfe; *optr++ = p; iptr++; for(x = 1; x < w; x++, iptr++, optr++){ p = ER[*iptr] & iptr[-W] & iptr[W]; if(!(iptr[-1] & 1)) p &= 0x7f; if(!(iptr[1] & 0x80)) p &= 0xfe; *optr = p; } p = ER[*iptr] & 0xfe & iptr[-W] & iptr[W]; if(!(iptr[-1] & 1)) p &= 0x7f; *optr++ = p; iptr++; } return ret; }
Операция диляции чуть похуже: т.к. здесь нужно проверять и крайние строки, и крайние столбцы, пришлось хитрить. Я знаю, что сложные макросы можно написать при помощи буста, но без интернета мне как-то не сподручно было документацию достать. Поэтому я пошел простейшим путем: внутренность функции вынес в отдельный файл (точно так же я поступлю и для маркировки связанных областей), а потом многократно включал его с разными #define'ами.
Вот сама функция, реализующая диляцию:
/** * Make morphological operation of dilation * @param image (i) - input image * @param W, H - size of image * @return allocated memory area with dilation of input image */ unsigned char *dilation(unsigned char *image, int W, int H){ if(W < 2 || H < 2) errx(1, "Dilation: image size too small"); if(!__Init_done) morph_init(); unsigned char *ret = Malloc(W*H, 1); int y = 0, w = W-1, h = H-1; // top of image, y = 0 #define IM_UP #include "dilation.h" #undef IM_UP // mid of image, y = 1..h-1 #include "dilation.h" // image bottom, y = h y = h; #define IM_DOWN #include "dilation.h" #undef IM_DOWN return ret; } А вот - ее внутренности:
/* * HERE'S NO ANY "FILE-GUARDS" BECAUSE FILE IS MULTIPLY INCLUDED! * I use this because don't want to dive into depths of BOOST * * multi-including with different defines before - is a simplest way to * modify simple code for avoiding extra if-then-else */ #if ! defined IM_UP && ! defined IM_DOWN // image without upper and lower borders OMP_FOR() for(y = 1; y < h; y++) #endif { unsigned char *iptr = &image[W*y]; unsigned char *optr = &ret[W*y]; unsigned char p = DIL[*iptr] #ifndef IM_UP | iptr[-W] #endif #ifndef IM_DOWN | iptr[W] #endif ; int x; if(iptr[1] & 0x80) p |= 1; *optr = p; optr++; iptr++; for(x = 1; x < w; x++, iptr++, optr++){ p = DIL[*iptr] #ifndef IM_UP | iptr[-W] #endif #ifndef IM_DOWN | iptr[W] #endif ; if(iptr[1] & 0x80) p |= 1; if(iptr[-1] & 1) p |= 0x80; *optr = p; } p = DIL[*iptr] #ifndef IM_UP | iptr[-W] #endif #ifndef IM_DOWN | iptr[W] #endif ; if(iptr[-1] & 1) p |= 0x80; *optr = p; optr++; iptr++; } Все совершенно просто и понятно.
Как видно, эти функции отлично распараллеливаются, т.е. в теории их можно было бы и кудой считать. Однако, так делать не нужно: на моем ноутбуке (i5-3210M CPU @ 2.50GHz x2cores x2hyper; 6GB RAM) без OpenMP эрозия/диляция изображения 2000х2000 выполнялась около трех миллисекунд; с OpenMP, раскидывая вычисления на все 4 ядра, скорость повысилась до ~1.7мс на операцию; а когда я еще -O3 указал, операции вообще за 0.7мс стали выполняться. Понятно, что куда будет считать медленнее, т.к. копирование памяти туда-сюда - довольно-таки дорогая операция. С другой стороны, если изображение уже в памяти видеокарты, и с ним надо что-то эдакое сотворить (т.е. помимо всяких эрозий-диляций еще выполнить какие-нибудь хорошо распараллеливаемые операции), можно будет на будущее и для куды написать реализацию. Но не сейчас, т.к., например, операция маркирования связанных областей вообще не параллелится: если даже разбить изображение на кусочки, то затем все равно потребуется делать полную перемаркировку, стыкуя куски на границах. Конечно, зубов бояться - … (вполне возможно, что все-таки пошустрей будет эта операция на очень больших изображениях), но реализацию пока отложу на потом.
В следующей заметке изложу эпопею поиска оптимального алгоритма выделения связанных областей. Но для начала надо отрихтовать "китайский" вариант, чтобы работал правильно. Ну и подумать насчет параллелизации (мало ли: вдруг на пару порядков быстрей будет).