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

Categories:

Поиск изофот. Часть 1: неправильный подход

Разрабатывая свой велосипед, о котором уже неоднократно писал, я наткнулся на необходимость поиска изофот с представлением их для начала хотя бы в виде ломаных (потом можно будет и сглаживание добавить). И вот здесь-то у меня возникли затруднения. Есть множество программ, реализующих поиск и даже сглаживание изофот, но вот внятного изложения алгоритма, реализующего это, на просторах интернета я не нашел. Ковыряться в исходниках GPL-программ, строящих изофоты, мне не хотелось (на это ушло бы гораздо больше времени), поэтому я решил для начала попробовать сделать свою реализацию метода шагающих квадратов. Использовать "классический" алгоритм я не захотел, т.к. он приводит к уменьшению количества найденных узлов (все-таки, я собираюсь использовать именно этот вариант, проводя сетку 2х2 пикселя, но об этом - в следующей публикации). Итак, у меня "родился" свой, неправильный метод.

Классический метод поиска изофот подразумевает то, что вы бинаризуете изображение для каждого уровня, заменяя пиксели с интенсивностью, меньшей Ii (i-й уровень изолиний), нулями, а остальные - единицами. Затем строится сетка, объединяющая в себе несколько пикселей (например, 2х2 или 3х3) и по ней строится матрица для реализации "шагающих квадратов": в зависимости от интенсивности в углах сетки очередной ячейке матрицы присваивается некоторое значение. После, при обработке конкретной изолинии, направление для проверки следующей ячейки выбирается, исходя из значения в предыдущей.

Метод несложный, но очень медленный: если нам надо построить 256 изофот, то придется 256 раз бинаризовать изображение и строить маски для "шагающих квадратов". Распараллелить (как очень хотелось бы) это не получится: если наша картинка имеет размер, скажем, 4000х4000 пикселей, то (помимо самой картинки) понадобится еще 256(количество изолиний) · [16МБ(размер бинаризованного изображения) + 4МБ(размер маски)] = 5120МБ. Вряд ли какой домашний компьютер обладает таким количеством оперативки, не говоря уже о видеокартах (а CUDA очень хороша). Следовательно, в данном случае распараллелить можно лишь нахождение изолиний для конкретного уровня Ii. Но, судя по моим прикидкам, ощутимого прироста производительности этот способ параллелизации не даст, т.к. каждый участок придется "стыковать" с соседним (что в параллели уже не сделать).

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

  1. предварительное изображение строится не бинаризацией исходного (для каждого изоуровня), а квантованием по заданному алгоритму определения изолиний; т.е. интенсивность в каждом пикселе преобразуется по нужному алгоритму распределения изоуровней (линейному, логарифмическому, степенному и т.п.) и посредством округления получается номер изоуровня;
  2. маска строится по квадратам 2х2 пикселя но так, что эти квадраты перекрываются, причем маска строится для всех уровней интенсивности в данном квадрате (похоже, именно это и было моей ошибкой);
  3. для отождествления, что же за изолинии проходят через данных квадрат, строились еще 2 маски: минимумы и максимумы номеров изоуровней в данном квадрате (т.к. через него может проходить несколько изолиний);
  4. далее в цикле я просматривал все пиксели маски и, как только находил ненулевую ячейку, проводил изолинии, начиная с минимального и заканчивая максимальным уровнем для нее.
Построение маски реализовалось простеньким кодом:

	#pragma omp parallel for
	for(y = 0; y < h1; y++){
		int x, idx = y*w1;
		unsigned char *in = levels + y*w, *out = mask + idx, *mi = mins + idx, *ma = maxs + idx;
		unsigned char m, a,b,c,d;
		for(x = 0; x < w1; x++, in++, out++, mi++, ma++){
			a = in[0], b = in[1], c = in[w], d = in[w+1];
			m = MIN(MIN(a,b), MIN(c,d));
			*mi = m;
			*out = ((a!=m)<<3) | ((b!=m)<<2) | ((c!=m)<<1) | (d!=m);
			*ma = MAX(MAX(a,b), MAX(a,b));
		}
	}
Здесь levels - массив квантованных уровней интенсивности; mask - маска, которую мы строим; mins и maxs - массивы граничных значений уровней сигнала в данной ячейке маски.

Сам цикл построения изолиний:


	contours = calloc(nLevl, sizeof(Contour)); // all-picture contours array
	for(y = 0; y < h1; y++)
		for(int x = 0; x < w1; x++)
			if(!process_it_(x, y, 255, ima->data)) goto endOF;
Функция process_it_ принимает аргументами координаты ячейки маски, верхний граничный уровень (в данном случае используется unsigned char, поэтому больше 255 изолиний сделать не получится) и указатель на массив - изображение. Возвращает она всегда TRUE, кроме ошибок работы с памятью - тогда возвращается FALSE.

Функция process_it_:


int process_it_(int x, int y, unsigned char restr, float *imdata){
	Contour *cCur;
	unsigned char lvl, *pt0;
	do{
		int i = y*w1+x;
		pt0 = mask + i;
		if(*pt0 == 0 || *pt0 > 14) break; // empty - go further
		// don't empty - go by the contour
		lvl = mins[i];
		if(lvl == restr || lvl > maxs[i]) break;
		if(!contours[lvl]){ // countour wasn't created - create it
			contours[lvl] = new_clist(lvl);
			if(!contours[lvl]){
				g_err(_("Can't allocate memory"));
				return FALSE;
			}
		}
		cCur = new_contour();
		if(!cCur){
			g_err(_("Can't allocate memory"));
			return FALSE;
		}
		cList_add(contours[lvl], cCur);
		if(!proc_contour(x, y, 0, cCur, imdata, lvl)) return FALSE;
	int on = cCur->N;
		if(!proc_contour(x, y, 1, cCur, imdata, lvl)) return FALSE;
	if(cCur->N != on){
	//	DBG("was %d pts, become %d pts, x1=%g", on, cCur->N,cCur->first->x);
		mask[y*w1+x] = 32;
	}
		if( (fabs(cCur->first->x - cCur->last->x) +
			fabs(cCur->first->y - cCur->last->y)) < 2.) // closed contour
				cCur->closed = TRUE;
	}while(*pt0 != 0 && *pt0 < 15); // cycle by all contours pass through given point
	return TRUE;
}
проверяет значение в данном пикселе маски (кстати, эту проверку надо бы в цикл перенести, чтобы пореже дергать вызов функции). Если значение говорит о том, что здесь проходит изолиния, идет проверка значения в массиве минимумов на соответствие граничному условию, от этой точки начинает формироваться контур со значением интенсивности, равным минимуму для данной ячейки. Контур строится как по, так и против часовой стрелки (если контур разомкнут и мы наткнулись на выпуклую часть). Далее производится проверка следующих уровней (если они есть). Если начальная и конечная точки контура близки друг к другу, контур помечается как замкнутый.

И вот обнаружена еще одна ошибка: из-за строчки mask[y*w1+x] = 32; другие контуры, если они проходят через данную точку, обрабатываться не будут!

Функция proc_contour:


int proc_contour(int x1, int y1, char flag,
				Contour *cCur, float *imdata,
				unsigned char lvl){
	int idx1 = y1*w1+x1;
	gboolean frst = flag;
	float Lvl = isoscale[lvl]; // intensity at isoline level lvl
	float *point;
	do{
		int iid = y1*w+x1;
		unsigned char m =   ((levels[iid]>lvl)<<3)   |
							((levels[iid+1]>lvl)<<2) |
							((levels[iid+w]>lvl)<<1) |
							(levels[iid+w+1]>lvl);
		if(!frst){
			if(m == 0 || m > 14) break;
			cPoint *p = new_point();
			if(!p) return FALSE;
			c_add(cCur, p, flag);
			// get coordinates p->x, p->y (by simplified linear approximation)
			point = imdata + iid; // LU pixel for a square
			// find coordinates
			getXY[m](p,point[0],point[1],point[w],point[w+1],Lvl);
			if(p->x < -0.5 || p->x > 1.5) p->x = 0.5;
			if(p->y < -0.5 || p->y > 1.5) p->y = 0.5;
			p->x += (float)x1;
			p->y += (float)y1;
		}
		// zero mask there's only 1 isoline pass through the square
		if(mins[idx1] >= maxs[idx1]-1) mask[idx1] = 15;
		else{
			if(lvl == mins[idx1]) mins[idx1]++;
			//else if(lvl == maxs[idx1]) maxs[idx1]--;
			else{
				process_it_(x1, y1, lvl, imdata);
				mins[idx1]++;
				if(mins[idx1] >= maxs[idx1]-1) mask[idx1] = 0;
			}
		}
		gboolean found = FALSE; // whether next isoline point was found
		int x2, y2;
		for(int i = 0; i < SQSIDES; i++){
		//for(int i = 0; i < 4; i++){
			unsigned char pp;
			// check right direction
			x2 = x1 + dxdy[i][0];
			y2 = y1 + dxdy[i][1];
			//x2 = x1 + DxDy[m][i][0];
			//y2 = y1 + DxDy[m][i][1];
			// is next square out of an picture?
			if(x2 < 0 || x2 >= w1 || y2 < 0 || y2 >= h1) continue;
			int ii = y2*w1+x2;
			pp = mask[ii];
			// has the next square a contour points?
			if(pp == 0 || pp > 14) continue;
			// is the contour alien?
			if(mins[ii] > lvl || maxs[ii] < lvl) continue;
			found = TRUE;
			break;
		}
		x1 = x2; y1 = y2;
		if(!found){
			break; // end of isoline
		}
		idx1 = y1*w1 + x1;
		frst = FALSE;
	}while(1); // cycle through the contour
	return TRUE;
}
В зависимости от значения аргумента flag она обходит контур по или против часовой стрелки. Координаты очередной точки контура находятся простой линейной интерполяцией, поэтому изломов контура не избежать (но это можно будет компенсировать сглаживанием). Линейной интерполяцией занимается массив функций getXY, вычисляющий координаты в зависимости от значения данной ячейки маски.

Как только мы при построении контура натыкаемся на пиксель, через который проходит несколько изолиний, снова вызывается функция process_it_ для построения всех прочих изолиний. И вот здесь-то вышеупомянутая ошибка разрывает изолинии. Планирую: провести работу над ошибками и, если ничего не получится, переделать методику (т.е. сделать ячейки маски непересекающимися или даже вообще использовать классический - медленный - алгоритм).

В следующей части повествования я надеюсь рассказать о реализации работающего алгоритма.

Tags: c, велосипедостроение, изолинии, изофоты, так делать не надо, шагающие квадраты
Subscribe

  • DHT22/DHT11 на STM32F103

    Добил шайтана! Сначала ожидал, что нужно будет полноценным захватом ШИМ пользоваться, но т.к. в протоколе неинформативная часть имеет постоянную…

  • Свеженькие железячки

    Получил сегодня с али ожидаемые железки, в т.ч. для восстановления моей файлопомойки. Во-первых, это блок питания на 16.8В, который я брал для…

  • BMP180 на STM32F103

    Добавил еще один сниппет — работа с BMP180 (датчик температуры и давления). Опять в даташите формула для вычисления "компенсированных" значений…

promo eddy_em august 17, 2019 12:33 3
Buy for 10 tokens
Юра намедни напечатал корпус для хронометра. Для первого блина получилось неплохо: И еще немного фотографий:
  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your reply will be screened

  • 0 comments