eddy_em: (hram nauki)
[personal profile] eddy_em
Итак, как я уже говорил в предыдущей записи, алгоритм китайцев у меня заработал. Вот — сравнение с моим:
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);
This account has disabled anonymous posting.
If you don't have an account you can create one now.
HTML doesn't work in the subject.
More info about formatting

If you are unable to use this captcha for any reason, please contact us by email at support@dreamwidth.org

May 2025

S M T W T F S
    123
45678910
11121314151617
1819202122 2324
25262728293031

Most Popular Tags

Style Credit

Expand Cut Tags

No cut tags
Page generated May. 29th, 2025 12:31 pm
Powered by Dreamwidth Studios