eddy_em: (hram nauki)
eddy_em ([personal profile] eddy_em) wrote2013-04-21 07:06 pm

Маркировка связанных областей, итог

Итак, как я уже говорил в предыдущей записи, алгоритм китайцев у меня заработал. Вот — сравнение с моим:
graph
Производительность, по вертикали — логарифм времени выполнения операции, по горизонтали — корень из количества пикселей на картинке. Синее, красное, розовое — мое; зеленое, голубое и желтое — китайцев. Верхняя пара кривых — для 50% соотношения пикселей изображения и фона; средняя — для 10%; нижняя — для 90%.
В ходе тестов, кстати, выяснилось, что мой алгоритм глючит: периодически он пропускает переименование меток (видимо, какой-то косяк в функции перемаркировки, но мне лень его искать: мой алгоритм настолько тормозной, что я его использовать не буду, выкину нафиг).
Пока что я не стал добавлять поиск 4-связанных областей по алгоритму китайцев, сделаю это позже. Ну и еще не сделано выделение объектов (это — примитивная операция, сделать несложно). В общем же все достаточно просто:
/**
 * label 4-connected components on image
 * (slow algorythm, but easy to parallel)
 *
 * @param I (i)    - image ("packed")
 * @param W,H,W_0  - size of the image (W - width in pixels, W_0 - width in octets)
 * @param Nobj (o) - number of objects found
 * @return an array of labeled components
 */
CCbox *cclabel4(unsigned char *Img, int W, int H, int W_0, size_t *Nobj){
	unsigned char *I = FC_filter(Img, W_0, H);
	#include "cclabling.h"
	FREE(I);
	return ret;
}
// label 8-connected components, look cclabel4
CCbox *cclabel8(unsigned char *I, int W, int H, int W_0, size_t *Nobj){
	#define LABEL_8
	#include "cclabling.h"
	#undef LABEL_8
	return ret;
}

И вот реализация алгоритма китайцев:
	CCbox *ret = NULL;
	size_t N = 0, // current label
		Ntot = 0, // number of found objects
		Nmax = W*H/4; // max number of labels
	size_t *labels = char2st(I, W, H, W_0);
	int w = W - 1, h = H - 1;
	int y;
	size_t *assoc = Malloc(Nmax, sizeof(size_t)); // allocate memory for "remark" array
	inline void remark(size_t old, size_t *newv){
		size_t new = *newv;
		if(old == new || assoc[old] == assoc[new])
			return;
		if(!assoc[old]){ // try to remark non-marked value
			assoc[old] = new;
			return;
		}
		// value was remarked -> we have to remark "new" to assoc[old]
		// and decrement all marks, that are greater than "new"
		size_t ao = assoc[old], an = assoc[new];
		// now we must check assoc[old]: if it is < new -> change *newv, else remark assoc[old]
		if(ao < an)
			*newv = old;
		else{
			size_t x = ao; ao = an; an = x; // swap values
		}
		OMP_FOR()
		for(int i = 1; i <= N; i++){
			size_t m = assoc[i];
			if(m < an) continue;
			if(m == an)
				assoc[i] = ao;
			else
				assoc[i]--;
		}
		Ntot--;
	}
	size_t *ptr = labels;
	for(y = 0; y < H; y++){  // FIXME!!!! throw out that fucking "if" checking coordinates!!!
		bool found = false;
		size_t curmark = 0; // mark of pixel to the left
		for(int x = 0; x < W; x++, ptr++){
			size_t curval = *ptr;
			if(!curval){ found = false; continue;}
			size_t *U = (y) ? &ptr[-W] : NULL;
			size_t upmark = 0; // mark from line above
			if(!found){ // new pixel, check neighbours above
				found = true;
				// now check neighbours in upper string:
				if(U){
					if(x && U[-1]){ // there is something in upper left corner -> use its value
						upmark = U[-1];
					}else // check point above only if there's nothing in left up
						if(U[0]) upmark = U[0];
					if(x < w && U[1]){ // there's something in upper right
						if(upmark){ // to the left of it was pixels
							remark(U[1], &upmark);
						}else
							upmark = U[1];
					}
				}
				if(!upmark){ // there's nothing above - set current pixel to incremented counter
					size_t *D = (y < h) ? &ptr[W] : NULL;
					if(  !(x && ((D && D[-1]) || ptr[-1]))   // no left neighbours
					  && !(x < w && ((D && D[1]) || ptr[1])) // no right neighbours
					  && !(D && D[0])){ // no neighbour down
						*ptr = 0; // erase this hermit!
						continue;
					}
					upmark = ++N;
					assoc[upmark] = ++Ntot; // refresh "assoc"
				}
				*ptr = upmark;
				curmark = upmark;
			}else{ // there was something to the left -> we must chek only U[1]
				if(U){
					if(x < w && U[1]){ // there's something in upper right
						remark(U[1], &curmark);
					}
				}
				*ptr = curmark;
			}
		}
	}
	// Step 2: rename markers
	// first correct complex assotiations in assoc
	OMP_FOR()
	for(y = 0; y < H; y++){
		size_t *ptr = &labels[y*W];
		for(int x = 0; x < W; x++, ptr++){
			size_t p = *ptr;
			if(!p){continue;}
			*ptr = assoc[p];
		}
	}
	FREE(assoc);
	FREE(labels);

А вот — так сказать, на память — мое нерабочее поделие:
	CCbox *ret = NULL;
	size_t N _U_ = 0, Ncur = 0;
	size_t *labels = char2st(I, W, H, W_0);
	int y;
	// Step 1: mark neighbours by strings
	size_t *ptr = labels;
	for(y = 0; y < H; y++){
		bool found = false;
		for(int x = 0; x < W; x++, ptr++){
			if(!*ptr){ found = false; continue;}
			if(!found){
				found = true;
				++Ncur;
			}
			*ptr = Ncur;
		}
	}
	// Step 2: fill associative array with remarking
	N = Ncur+1; // size of markers array (starting from 1)
	size_t i, *assoc = Malloc(N, sizeof(size_t));
	for(i = 0; i < N; i++) assoc[i] = i; // primary initialisation
	// now we should scan image again to rebuild enumeration
	// THIS PROCESS AVOID PARALLELISATION???
	Ncur = 0;
	size_t h = H-1
	#ifdef LABEL_8
		,w = W-1;
	#endif
	;
	inline void remark(size_t old, size_t *newv){
		size_t new = *newv;
		if(assoc[old] == new)
			DBG("(%zd)\n", new);
		if(assoc[old] == old){ // remark non-remarked value
			assoc[old] = new;
			return;
		}
		// value was remarked -> we have to remark "new" to assoc[old]
		// and decrement all marks, that are greater than "new"
		size_t ao = assoc[old];
		// now we must check assoc[old]: if it is < new -> change *newv, else remark assoc[old]
		if(ao < new)
			*newv = ao;
		else{
			size_t x = ao; ao = new; new = x; // swap values
		}
		int xx = old + W / 2;
		if(xx > N) xx = N;
		OMP_FOR()
		for(int i = 1; i <= xx; i++){
			size_t m = assoc[i];
			if(m < new) continue;
			if(m == new)
				assoc[i] = ao;
			else if(m <= Ncur)
				assoc[i]--;
		}
		Ncur--;
	}
	for(y = 0; y < H; y++){  // FIXME!!!! throw out that fucking "if" checking coordinates!!!
		size_t *ptr = &labels[y*W];
		bool found = false;
		for(int x = 0; x < W; x++, ptr++){
			size_t curval = *ptr;
			if(!curval){ found = false; continue;}
			if(found) continue; // we go through remarked pixels
			found = true;
			// now check neighbours in upper and lower string:
			size_t *U = (y) ? &ptr[-W] : NULL;
			size_t *D = (y < h) ? &ptr[W] : NULL;
			size_t upmark = 0; // mark from line above
			if(U){
				#ifdef LABEL_8
				if(x && U[-1]){ // there is something in upper left corner -> make remark by its value
					upmark = assoc[U[-1]];
				}else // check point above only if there's nothing in left up
				#endif
					if(U[0]) upmark = assoc[U[0]];
				#ifdef LABEL_8
				if(x < w) if(U[1]){ // there's something in upper right
					if(upmark){ // to the left of it was pixels
						remark(U[1], &upmark);
					}else
						upmark = assoc[U[1]];
				}
				#endif
			}
			if(!upmark){ // there's nothing above - set current pixel to incremented counter
				#ifdef LABEL_8  // check, whether pixel is not single
				if(  !(x && ((D && D[-1]) || ptr[-1]))   // no left neighbours
				  && !(x < w && ((D && D[1]) || ptr[1])) // no right neighbours
				  && !(D && D[0])){ // no neighbour down
					*ptr = 0; // erase this hermit!
					continue;
				}
				#endif
				upmark = ++Ncur;
			}
			// now remark DL and DR corners (bottom pixel will be checked on next line)
			if(D){
				#ifdef LABEL_8
				if(x && D[-1]){
					remark(D[-1], &upmark);
				}
				if(x < w && D[1]){
					remark(D[1], &upmark);
				}
				#else
				if(D[0]) remark(D[0], &upmark);
				#endif
			}
			 // remark current
			remark(curval, &upmark);
		}
	}
	// Step 3: rename markers
	// first correct complex assotiations in assoc
	OMP_FOR()
	for(y = 0; y < H; y++){
		size_t *ptr = &labels[y*W];
		for(int x = 0; x < W; x++, ptr++){
			size_t p = *ptr;
			if(!p){continue;}
			*ptr = assoc[p];
		}
	}
	FREE(assoc);
	FREE(labels);