![[personal profile]](https://www.dreamwidth.org/img/silk/identity/user.png){kind=small}
Поиск изофот. Часть 1: неправильный подход
Классический метод поиска изофот подразумевает то, что вы бинаризуете изображение для каждого уровня, заменяя пиксели с интенсивностью, меньшей Ii (i-й уровень изолиний), нулями, а остальные - единицами. Затем строится сетка, объединяющая в себе несколько пикселей (например, 2х2 или 3х3) и по ней строится матрица для реализации "шагающих квадратов": в зависимости от интенсивности в углах сетки очередной ячейке матрицы присваивается некоторое значение. После, при обработке конкретной изолинии, направление для проверки следующей ячейки выбирается, исходя из значения в предыдущей.
Метод несложный, но очень медленный: если нам надо построить 256 изофот, то придется 256 раз бинаризовать изображение и строить маски для "шагающих квадратов". Распараллелить (как очень хотелось бы) это не получится: если наша картинка имеет размер, скажем, 4000х4000 пикселей, то (помимо самой картинки) понадобится еще 256(количество изолиний) · [16МБ(размер бинаризованного изображения) + 4МБ(размер маски)] = 5120МБ. Вряд ли какой домашний компьютер обладает таким количеством оперативки, не говоря уже о видеокартах (а CUDA очень хороша). Следовательно, в данном случае распараллелить можно лишь нахождение изолиний для конкретного уровня Ii. Но, судя по моим прикидкам, ощутимого прироста производительности этот способ параллелизации не даст, т.к. каждый участок придется "стыковать" с соседним (что в параллели уже не сделать).
Итак, распараллеливание процесса сразу же было отметено как слишком сложное и/или дорогостоящее. Но хотелось же как-то ускорить процесс. И я решил сделать это так:
- предварительное изображение строится не бинаризацией исходного (для каждого изоуровня), а квантованием по заданному алгоритму определения изолиний; т.е. интенсивность в каждом пикселе преобразуется по нужному алгоритму распределения изоуровней (линейному, логарифмическому, степенному и т.п.) и посредством округления получается номер изоуровня;
- маска строится по квадратам 2х2 пикселя но так, что эти квадраты перекрываются, причем маска строится для всех уровней интенсивности в данном квадрате (похоже, именно это и было моей ошибкой);
- для отождествления, что же за изолинии проходят через данных квадрат, строились еще 2 маски: минимумы и максимумы номеров изоуровней в данном квадрате (т.к. через него может проходить несколько изолиний);
- далее в цикле я просматривал все пиксели маски и, как только находил ненулевую ячейку, проводил изолинии, начиная с минимального и заканчивая максимальным уровнем для нее.
#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));
}
}
Сам цикл построения изолиний:
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_:
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;
}
Как только мы при построении контура натыкаемся на пиксель, через который проходит несколько изолиний, снова вызывается функция process_it_ для построения всех прочих изолиний. И вот здесь-то вышеупомянутая ошибка разрывает изолинии. Планирую: провести работу над ошибками и, если ничего не получится, переделать методику (т.е. сделать ячейки маски непересекающимися или даже вообще использовать классический - медленный - алгоритм).
В следующей части повествования я надеюсь рассказать о реализации работающего алгоритма.