eddy_em

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

Бла-бла

Точные определения эрозии (сжатия) и диляции (расширения) бинарных изображений можно прочитать в википедии. Вкратце: эти операции чем-то похожи на свертку изображения с некоей маской, основная же разница от свертки в том, что в случае эрозии сложение заменяется логической операцией "И", а в случае диляции — операцией "ИЛИ". Т.е. эрозия оставляет на изображении лишь те пиксели, которые имеют абсолютно всех соседей по аналогии с маской, а диляция же добавляет каждому пикселю соседей из маски. Понятное дело, операции эти — совершенно необратимые.
Кстати, больший интерес представляют даже не сами по себе эрозия и диляция, а их комбинации: закрытие (эрозия диляции) и открытие (диляция эрозии). Эти операции наглядно показывают необратимость отдельных операций: эрозия диляции позволяет получить вроде бы первоначальный объект, но с закрытыми "дырками"; диляция же эрозии дает вроде бы первоначальный объект, но с удалением небольших выбросов и подчисткой индивидуальных пикселей (и даже небольших объектов, если наша маска достаточно большая).
Маски могут иметь самые разнообразные формы. Мне же интересно было реализовать эрозию и диляцию наиболее часто употребляемой маской — крестообразной (т.е. для 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.

ПРОВЕРИТЬ, ЧТО БЫСТРЕЙ: (A&0X80)>>8 ИЛИ РАЗЫМЕНОВАНИЕ УКАЗАТЕЛЯ (ПО ТАБЛИЦЕ).

Вот эту проверку я было думал сделать, но когда обнаружил, что операция сама по себе так шустро выполняется, что дополнительная оптимизация вряд ли поднимет скорость больше, чем на порядок, понял, что ничего переделывать и не нужно. Я думаю, что если результат тяжелой оптимизации повышает производительность меньше, чем на порядок, такая оптимизация нафиг не нужна! Если же эта оптимизация дается малой кровью (скажем, минут за 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;
}

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

/**
 * 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++;
}

-O3