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

Category:

Морфологические операции: эрозия и диляция.

Итак, как я уже говорил в предыдущей записи, за время, бесцельно проведенное в Нижнем Новгороде, кое-что полезное я таки сделал. В этой записи расскажу о реализации операций эрозии и диляции.

Бла-бла

Точные определения эрозии (сжатия) и диляции (расширения) бинарных изображений можно прочитать в википедии. Вкратце: эти операции чем-то похожи на свертку изображения с некоей маской, основная же разница от свертки в том, что в случае эрозии сложение заменяется логической операцией "И", а в случае диляции — операцией "ИЛИ". Т.е. эрозия оставляет на изображении лишь те пиксели, которые имеют абсолютно всех соседей по аналогии с маской, а диляция же добавляет каждому пикселю соседей из маски. Понятное дело, операции эти — совершенно необратимые.
Кстати, больший интерес представляют даже не сами по себе эрозия и диляция, а их комбинации: закрытие (эрозия диляции) и открытие (диляция эрозии). Эти операции наглядно показывают необратимость отдельных операций: эрозия диляции позволяет получить вроде бы первоначальный объект, но с закрытыми "дырками"; диляция же эрозии дает вроде бы первоначальный объект, но с удалением небольших выбросов и подчисткой индивидуальных пикселей (и даже небольших объектов, если наша маска достаточно большая).
Маски могут иметь самые разнообразные формы. Мне же интересно было реализовать эрозию и диляцию наиболее часто употребляемой маской — крестообразной (т.е. для 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.
Таким образом, алгоритм получается следующим:
  1. подменяем E маской[E] (эрозия внутренних битов)
  2. делаем &: E = E & B & H
  3. корректируем старший бит E по младшему D и младший бит E по старшему F.
ПРОВЕРИТЬ, ЧТО БЫСТРЕЙ: (A&0X80)>>8 ИЛИ РАЗЫМЕНОВАНИЕ УКАЗАТЕЛЯ (ПО ТАБЛИЦЕ).
Вот эту проверку я было думал сделать, но когда обнаружил, что операция сама по себе так шустро выполняется, что дополнительная оптимизация вряд ли поднимет скорость больше, чем на порядок, понял, что ничего переделывать и не нужно. Я думаю, что если результат тяжелой оптимизации повышает производительность меньше, чем на порядок, такая оптимизация нафиг не нужна! Если же эта оптимизация дается малой кровью (скажем, минут за 5), то почему бы и нет (но этот случай явно не из разряда пятиминутных).

Диляция крестообразной маской:
00000     00110
00110     01111
01010 =>  11111
01100     11110
00000     01100
Здесь противоположная ситуация: вместо "&" делаем "|":
  1. подменяем E маской (диляция внутренних битов),
  2. делаем |: E = E | B | H,
  3. корректируем старший и младший биты операцией |.
  4. Реализация

    Преобразование бинарного изображения в "упакованное" реализуется такой вот элементарщиной:
    /**
     * 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мс стали выполняться. Понятно, что куда будет считать медленнее, т.к. копирование памяти туда-сюда — довольно-таки дорогая операция. С другой стороны, если изображение уже в памяти видеокарты, и с ним надо что-то эдакое сотворить (т.е. помимо всяких эрозий-диляций еще выполнить какие-нибудь хорошо распараллеливаемые операции), можно будет на будущее и для куды написать реализацию. Но не сейчас, т.к., например, операция маркирования связанных областей вообще не параллелится: если даже разбить изображение на кусочки, то затем все равно потребуется делать полную перемаркировку, стыкуя куски на границах. Конечно, зубов бояться — … (вполне возможно, что все-таки пошустрей будет эта операция на очень больших изображениях), но реализацию пока отложу на потом.

    В следующей заметке изложу эпопею поиска оптимального алгоритма выделения связанных областей. Но для начала надо отрихтовать "китайский" вариант, чтобы работал правильно. Ну и подумать насчет параллелизации (мало ли: вдруг на пару порядков быстрей будет).
Tags: c, snippets, обработка изображений
Subscribe

  • Хеши строковых команд для МК

    Долго я к этому шел, но, похоже, пора уже: однобуквенные команды сложно запоминать (особенно если команд толпа, и большая часть с этими буквами…

  • OBS studio

    В общем, надоумили меня попробовать трансляцию в youtube. Напрямую скринкасты он писать не умеет, но может забирать поток с промежуточного…

  • M$ teams…

    Начал с сегодняшнего дня студентам ЮФУ удаленно лекции читать. У них все завязано на различные корпорации зла. И базовая работа - через teams. ОК,…

promo eddy_em september 3, 12:13 8
Buy for 10 tokens
Уже больше полугода занимаюсь разработкой, вот, наконец-то в мастерских взялись за меня и начали выдавать первые детали. Сегодня сделал тестовую сборку (как обычно, местами пришлось "доработать напильником"): Пока прибор без названия (да и как-то не лезет в голову ничего, у меня нет…
  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic
  • 0 comments