eddy_em | Выделение 4-связных компонент на изображении

Я уже давным-давно писал об этом алгоритме, но когда понадобилось его однозначно и надежно применить, оказалось, что на некоторых тестовых изображениях он давал сбой. В конце-концов, я нашел, где "собака порылась": я неправильно пересчитывал данные в массиве связей, что оставляло некоторые компоненты "сиротами" (хотя фактически они были связаны с кем-то еще).

Вот такая функция в итоге получилась:

/**
 * label 4-connected components on image
 * (slow algorythm, but easy to parallel)
 *
 * @param I (i)    - image ("packed")
 * @param W,H      - size of the image (W - width in pixels)
 * @param CC (o)   - connected components boxes
 * @return an array of labeled components
 */
size_t *cclabel4(uint8_t *Img, int W, int H, ConnComps **CC){
    size_t *assoc;
    // check table and rename all "oldval" into "newval"
    inline void remark(size_t newval, size_t oldval){
        // find the least values
        do{newval = assoc[newval];}while(assoc[newval] != newval);
        do{oldval = assoc[oldval];}while(assoc[oldval] != oldval);
        // now change larger value to smaller
        if(newval > oldval){
            assoc[newval] = oldval;
        }else{
            assoc[oldval] = newval;
        }
    }
    if(W < MINWIDTH || H < MINHEIGHT) return NULL;
    uint8_t *f = filter4(Img, W, H); // remove all non 4-connected pixels
    size_t *labels = bin2ST(f, W, H);
    FREE(f);
    size_t Nmax = W*H/4; // max number of 4-connected labels
    assoc = MALLOC(size_t, Nmax); // allocate memory for "remark" array
    size_t last_assoc_idx = 1; // last index filled in assoc array
    for(int y = 0; y < H; ++y){
        bool found = false;
        size_t *ptr = &labels[y*W];
        size_t curmark = 0; // mark of pixel to the left
        for(int x = 0; x < W; ++x, ++ptr){
            if(!*ptr){found = false; continue;} // empty pixel
            size_t U = (y) ? ptr[-W] : 0; // upper mark
            if(found){ // there's a pixel to the left
                if(U && U != curmark){ // meet old mark -> remark one of them in assoc[]
                    remark(U, curmark);
                    curmark = U; // change curmark to upper mark (to reduce further checks)
                }
            }else{ // new mark -> change curmark
                found = true;
                if(U) curmark = U; // current mark will copy upper value
                else{ // current mark is new value
                    curmark = last_assoc_idx++;
                    assoc[curmark] = curmark;
                }
            }
            *ptr = curmark;
        }
    }
    size_t *indexes = MALLOC(size_t, last_assoc_idx); // new indexes
    size_t cidx = 1;
    for(size_t i = 1; i < last_assoc_idx; ++i){
        // search new index
        register size_t realnew = i, newval = 0;
        do{
            realnew = assoc[realnew];
            if(indexes[realnew]){ // find least index
                newval = indexes[realnew];
                break;
            }
        }while(assoc[realnew] != realnew);
        if(newval){
            indexes[i] = newval;
            continue;
        }
        // enter next label
        indexes[i] = cidx++;
    }
    Box *boxes = MALLOC(Box, cidx);
    OMP_FOR()
    for(size_t i = 1; i < cidx; ++i){ // init borders
        boxes[i].xmin = W;
        boxes[i].ymin = H;
    }
// parallel:
#pragma omp parallel shared(boxes)
    {
        Box *l_boxes = MALLOC(Box, cidx);
        for(size_t i = 1; i < cidx; ++i){ // init borders
            l_boxes[i].xmin = W;
            l_boxes[i].ymin = H;
        }
        #pragma omp for nowait
        for(int y = 0; y < H; ++y){
            size_t *lptr = &labels[y*W];
            for(int x = 0; x < W; ++x, ++lptr){
                if(!*lptr) continue;
                register size_t mark = indexes[*lptr];
                *lptr = mark;
                Box *b = &l_boxes[mark];
                ++b->area;
                if(b->xmax < x) b->xmax = x;
                if(b->xmin > x) b->xmin = x;
                if(b->ymax < y) b->ymax = y;
                if(b->ymin > y) b->ymin = y;
            }
        }
        #pragma omp critical
        for(size_t i = 1; i < cidx; ++i){
            Box *ob = &boxes[i], *ib = &l_boxes[i];
            if(ob->xmax < ib->xmax) ob->xmax = ib->xmax;
            if(ob->xmin > ib->xmin) ob->xmin = ib->xmin;
            if(ob->ymax < ib->ymax) ob->ymax = ib->ymax;
            if(ob->ymin > ib->ymin) ob->ymin = ib->ymin;
            ob->area += ib->area;
        }
        FREE(l_boxes);
    }
    FREE(assoc);
    FREE(indexes);
    if(CC){
        *CC = MALLOC(ConnComps, 1);
        (*CC)->Nobj = cidx; (*CC)->boxes = boxes;
    }else{
        FREE(boxes);
    }
    return labels;
}

Помимо маркировки областей она еще и для ускорения процесса распознавания возвращает массив вспомогательных данных:

// simple box with given borders
typedef struct{
    uint16_t xmin;
    uint16_t xmax;
    uint16_t ymin;
    uint16_t ymax;
    uint32_t area; // total amount of object pixels inside the box
} Box;

typedef struct{
    size_t Nobj;
    Box *boxes;
} ConnComps;

Структура Box отражает данные конкретной области (границы и количество пикселей в объекте). А структура ConnComps характеризует изображение: по сути она представляет собой массив Box длины Nobj (но стоит иметь в виду, что Nobj — количество объектов без единицы, т.к. массив с индексом нуль пустой, нумерация ведется от единицы, ведь нуль — фон).
Теперь значительно проще выделять полезные объекты: сначала выбрасываем те, у которых поле area слишком маленькое (скажем, меньше пяти). Далее можно форм-фактор вычислить, равный (xmax-xmin+1)/(ymax-ymin+1): если он не лежит вблизи единицы, то это — либо сильно искаженные комой звезды, либо не звезды. А если area меньше, чем 3/4 площади бокса, то это — тоже явно не кружочек (хотя, может быть сильно расфокусированным изображением, здесь уже не все так однозначно).
Выделение 14233 объектов на тестовом снимке 4k×4k с 0.5-м телескопа выполнялось в среднем 32мс (мне кажется, что это ужасно много) на моем рабочем компьютере (восемь потоков на i7-6700, 32ГБ оперативки). Если убрать параллелизацию при окончательной маркировке и вычислении боксов, получится проигрыш в 3мс. Было бы интересно сравнить с другими алгоритмами, но мне уже откровенно некогда этим заниматься.
Что до ошибок в предыдущей версии, можно сравнить функцию remark здесь и там. Оказалось совершенно без надобности бродить по массиву assoc и заменять все индексы на новый: достаточно пройтись по списку связей до конца (когда индекс массива совпадает со значением) и заменить его на новый. Встречаем связь — ищем "корни", меняем бóльший на меньший.
Следующей моей ошибкой было пытаться прямо в этом массиве (без вспомогательного) пересмотреть все индексы, чтобы по порядку перенумеровать блоки. Я завел массив indexes, в котором и хранил окончательные связи (indexes[n] дает номер блока для значения n с изображения). Здесь уже поиск "корневого" значения в массиве assoc не обязательно проводить, если в соответствующей ячейке indexes есть отличное от нуля число (оно и равно "корневому" для этого значения). В конце-концов величина cidx и равна размеру результирующего массива (1 + количество найденных связных областей).
В последнем цикле нельзя, к сожалению, для подсчета объема каждой области написать

#pragma omp parallel shared(boxes) reduction(+: boxes->area)

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

Было бы интересно, если кто-нибудь сравнил бы быстродействие этого с чем-нибудь общеупотребляемым (скажем, OpenCV).