eddy_em

Итак, алгоритм китайцев заработал и показал вполне неплохие результаты (но об этом — в следующей записи). Здесь же я просто скопирую то, что написал за время пребывания в НН и внуковском аэропорту.

Вернемся к тому же общему случаю: искомый пиксель нахдится в байте E.

ABC
DEF
GHK

Фильтрация 4-связанных областей схожа с операцией диляции за исключением собственно логики - в отсечении пикселей, не являющихся 4-связанными, необходимо провести операцию:

E4 = E & ( B|H|(E<<1)|(E>>1)|(D<<7)|(F>>7) )

Т.е. мы оставляем пиксель в E лишь в том случае, если хотя бы один из его 4-связанных соседей присутствует.
Соответственно, для крайних пикселей (крайние столбцы/строки) выражение для E4 будет видоизменяться (для крайних строк - наличием/отсутствием B,H; для крайних столбцов - наличием/отсутствием выражений с D, F).
Для упрощения записи можно было бы сделать макрос при помощи boost, но что-то лень мне в дебри буста погружаться, поэтому проблему решаю простым копипастом.
Кстати, здесь нет проверки A,C,G и K, иначе пришлось бы еще больше кода наворачивать (добавились бы особые случаи для углов изображения).

Общее выражение для поиска 8-связанных областей несколько видоизменится:

пусть {x} = (x|(x<<1)|(x>>1)), тогда
E8 = E & ( {B|H} | ((E<<1)|(E>>1)) | ((C|F|K)>>7) | ((A|D|G)<<7) )

Вычислительная сложность сильно возрастает, причем многие операции вычисляются по два раза. Пожалуй, можно попытаться немного ускорить вычисления, если заранее подготовить матрицу {x} (и матрицы (x<<7) и (x>>7) ???).
/* (сомнительно) Аналогично, заранее подготовив матрицы (x<<7) и (x>>7) можно
* попытаться ускорить поиск E4 */

Когда мы уже имеем отфильтрованное изображение, нужно "всего лишь" пробежаться по нему в поисках областей, со всех сторон ограниченных нулями. Грубо говоря, мы маркируем 8-связанные области без проверки на 8-связанность.
Первое, что приходит на ум - двухпроходная маркировка: сначала мы пробегаемся по связанной области, постепенно расширяя границы заключающего ее прямоугольника. Затем мы выделяем память для этого кусочка и заполняем его, пробегаясь по нашей области еще раз (ну или делаем полную копию содержимого бокса, а затем выкидываем все, что не принадлежит искомой связанной области).
Еще вариант - сразу выделять область памяти (вспомогательную) размером со ВСЕ первоначальное изображение (гениально!). В этом случае мы сможем сразу же копировать в эту область памяти нашу х-связанную область, а затем, когда она целиком будет скопирована и будут ясны параметры описывающего ее прямоугольника, выделяем уже память именно под эту область и копируем в нее содержимое бокса.
Логично было бы преобразовать картинку до анализа в бинарный вид, поэтому сначала будем предполагать, что мы имеем дело с бинарным изображением (а не "сжатым" битовым).
Теперь каждый пиксель нашей матрицы соответствует одному пикселю изображения.
На ум приходит использовать метод шагающих квадратов, однако, он применяется для обхода границ -> мало чем поможет (разве что первоначальным определением размеров бокса, но я уже решил делать лишь один проход обнаружения связанных областей). Интернета под рукой нет, чтобы подсмотреть наиболее эффективный способ, поэтому буду изобретать велосипед.

Нулевое приближение

/*
Работаем с копией изображения, т.к. уже пройденные значения будем обнулять.
*/
Итак, выделим память под промежуточный буфер (сразу обнулив ее) и начнем проходить изображение с левого верхнего угла. Как только мы натыкаемся на единицу в E, записываем в соответствующие координаты буфера 1. Далее смотрим, куда нам нужно двигаться. Считаем, что в A,B,C и D нули (т.к. там - уже пройденная область). Нам необходимо проверить пиксели F, G, H и K. Теперь каждый из них становится E. Если кто-то из них установлен в 1, проверяем его соседей и т.д. Однако, на всех шагах, отличных от первого, мы уже обязаны проверить ВСЕХ СОСЕДЕЙ, т.к. не факт, что наша область не представляет собою какой-нибудь извращенной фигуры с отростками и многочисленными полостями.
На ум приходит рекурсия, однако, учитывая то, что изображения могут быть довольно-таки огромными, с горем расстаюсь с этой идеей. А так было бы хорошо: рекурсия сделала бы очень простым этот алгоритм. Но регистры не бесконечные.
Кстати, если не преобразовывать изображение в бинарное, то вполне возможны случаи, когда один и тот же байт будет проходиться вплоть до четырех раз!
После копирования каждого пикселя с 1 в результирующий буфер, "стираем" его с оригинала (сбросом в 0). По окончанию операций обнуляем буфер.
Тут мне опять пришла в голову идея о шагающих квадратах: а что, если в первый проход пробежаться по границе области, определить описывающий ее прямоугольник, выделить буфер по размеру прямоугольника (это - как раз искомое возвращаемое функцией значение); после этого, проходя по объекту, пиксель за пикселем копировать его в буфер, "стирая" с оригинала. Правда, здесь опять та же самая проблема: нет гарантии, что за один проход все будет сделано.

"Гюйгенс".

Один из вариантов прохождения по изображению в поисках связанных областей (который, собственно, первым пришел мне в голову еще вчера, но лучше всего, конечно, он бы выглядел в рекурсивном исполнении).

лирика по поводу рекурсии
Еще раз, чтобы убедиться лишний раз в том, что рекурсия невозможна, я набросал элементарнейшую проверялку:

#include <stdio.h>

int incr(int a){
    fprintf(stderr, "%d ", a++);
    incr(a);
    return a;
}

main(){
    incr(1);
    return 0;
}

Проверялочка эта должна вылетать с сегфолтом как только скончаются все ресурсы на рекурсию. Получилось вот что:
... 261504 Ошибка сегментирования (core dumped)
... 261350 Ошибка сегментирования (core dumped)
... 261333 Ошибка сегментирования (core dumped)

ЕМНИП, что-то подобное было на ЛОРе. И, возможно, предел этот даже как-то регулируется. Число это подозрительно близко к 256к (262144): возможно, оно так и есть (просто еще кучу места в стеке занимают всякие вспомогательные буферы из stdio).
Я немного видоизменил программку: она теперь еще и по 1к выделяет:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int incr(int a){
    char *x = malloc(1024);
    memset(x, 1, 1024);
    fprintf(stderr, "%d ", a++);
    incr(a);
    return a;
}

main(){
    incr(1);
    return 0;
}

Теперь цифры "более другие":
... 174258 Ошибка сегментирования (core dumped)
... 174372 Ошибка сегментирования (core dumped)
... 174325 Ошибка сегментирования (core dumped)

Кто сожрал 85к? Почему именно 85к? Заменяю выделение памяти на char x = 2; - получаю то же самое. Добавление переменных ничего не меняет. В любом случае, все-таки про рекурсию придется забыть: если искомая область будет иметь площадь, превышающую эдак 170 тысяч пикселей (а это не так уж и много), то я получу большой фигвам.

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

"Гюйгенс модифицированный"

Для ограничения одной-единственной вспомогательной матрицей, можно попробовать оставить одну "волну": от первоначальной точки проводим расходящуюся во все стороны "волну". Пусть для простоты она будет "квадратной". Пробегаясь по всем точкам очередного фронта "волны", кроме особых (о них - позже, их тоже надо контролировать), "стираем" единицы с изображения, перенося их в маску. Как только натыкаемся на нуль, если до этого тоже был нуль, заносим единичку в матрицу спец.маски в соответствующий пиксель, а также в пиксель, соседний с текущим, но расположеный на "волне" с радиусом на 1 больше (т.е. для левой поверхности переносим пиксель влево на 1, для правой - вправо на 1 и т.п.). Если натыкаемся на 1, делаем обнуление соответствующего пикселя на "волне" с радиусом на 1 больше + его предыдущего соседа в спец.маске. Таким образом, "тени" получаются постепенно сходящимися. Но все равно, остается возможность того, что эти "тени" что-нибудь да сокроют, не говоря уже о том, что эта "волна" только в очень редких случаях сразу выявит хотя бы внешний контур искомой фигуры. Именно для этого и нужна наша спец.маска: как только "кончится" искомая фигура по "большой волне", нужно будет просмотреть все пиксели маски (внутри уже получившегося прямоугольника + 1 по радиусу, т.к. снаружи все равно нули), пробегая по их внешней границе: как только обнаруживаем 1 в изображении, "выпускаем" новую "волну". Пока бежим первый раз, пропускаем все единички в маске, копируя их на "волну" маски с большим радиусом.
Метод получился опять каким-то уж очень извращенным. Глаза уже слипаются, поэтому пока - отбой.

первое приближение к решению

Итак, сегодня я накачал всяких статеек. Почитал. Обычный алгоритм выделения связанных областей (как я мог о нем не вспомнить?) жутко тормозной: он требует как минимум четыре прохода (1 - маркировка соседних областей, 2 - построение списка эквивалентности номеров, 3 - переименование областей, 4 - выделение объектов). При этом на первом шаге мы должны проверять, были ли до пикселя нули, на втором - верхнего соседа (и его боковых, если нужно заполнить 8-связанную область), на третьем - сверяться с таблицей, на четвертом - многократно сканировать изображение, пока все объекты не будут выделены. Конечно, четвертый шаг можно укоротить: на третьем же шаге выделить массив прямоугольников и обновлять размеры соответствующего прямоугольника. Тогда четвертый шаг сведется к сканированию внутренностей каждого прямоугольника с заполнением его маски единицами там, где в массиве номеров пиксель == N.
Процедура перемаркировки такова: проходя матрицу с маркерами, проверяем каждый пиксель. Если пиксель нулевой, устанавливаем флаг обнаружения в false и продолжаем. Если флаг == true, это значит, что данный номер мы уже перемаркировали - продолжаем. Если же пиксель - "свеженайденный", надо проверить все 6 пикселей (для 8-связных областей) в предыдущем и последующем рядах (снизу - т.к. возможен вариант коротенького кусочка, а то и единичного пикселя над длинной линией - из-за флага обнаружения мы пролетим под ним, а вот на этапе сканирования его, изучая то, что под ним, сделаем правильную ремаркировку). Вполне возможно, что нет необходимости проверять сразу и пиксель над, и под нашим: пиксель точно под ним по-любому будет проверен (он - либо часть более длинной линии, которую мы обнаружим в DL или DR углах, либо часть вертикальной, которая обнаружится при сканировании нижней линии).
А вообще, полезно рассмотреть, какие варианты расположения соседей могут быть, чтобы не пропустить мелочей (нумерация не важна, поэтому просто 1/0; кроме того, т.к вариантов аж 256, рассмотрим лишь некоторые): (X - что угодно, U - единица в одном из пикселей, V - единица в одном и более пикселей)

 A   B   C   D   E   F   G   H   I   J   K   L   M   N   O   P   Q   R
XXX 000 UUU 000 101 101 11X 000 010 000 000 000 000 000 000 000 000 000
11X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X 01X
XXX 000 000 VVV 000 VVV VVV 010 000 000 000 000 000 000 000 000 000 000

A - просто проходим мимо
B - инкрементируем счетчик уникальных объектов и перемаркируем текущий номер
C - ремаркируем текущее значение на то, которое содержится в массиве с индексом,
	имеющим значение счетчика в пикселе с единичкой
D - то же, что в B + ремаркировка значений в нижней строке
E - ремаркировка текущего + правого верхнего счетчика по массиву с индексом,
	равным значению счетчика в левом верхнем пикселе
F - ремаркировка всех по левому верхнему
G - для упрощения можно сделать эквивалентом F: нет нужды усложнять алгоритм ради
	пары лишних операций
H - нет нужды проверять этот вариант, т.к. на стадии I он ремаркируется
I - заменяет H

Итак, общие тенденции:

Если в верхней строке ничего нет, проверяем, не является ли наше текущее значение счетчика уже перемаркированным: если является, оставляем как есть, иначе - счетчик уникальных объектов инкрементируется и текущее значение перемаркируется им; объекты внизу (если есть) тоже перемаркируются им;
если в левом верхнем углу !=0, все номера перемаркируются им, а пиксель над нами проверять нет нужды;
если обнаружен !=0 ровно над нами, правый верхний проверять нет нужды, а оставшиеся перемаркируются им.

При организации процедуры переименования я натолкнулся на следующее: если ряд X переименовывается в нижележащий Y, а потом Y переименовывается в Z, то ряд X так и остается переименованным в Y. Следовательно, при переименовании в массиве ассоциаций на стадии перенумерования надо проверять, выполняется ли равенство assoc[assoc[X]] == assoc[X], т.е. зафиксировать сложные связи.
Оказалось не все так просто с переименованием: надо будет нарисовать на бумажке, что-то засыпаю уже. На сегодняшний день ничего не планирую (спать надо лечь рано, т.к. завтра в 5 утра вставать), так что, пожалуй, закруглюсь и буду спать. Сегодня еще вопросы с криостатом утрясти надо будет.
Пока засыпал, придумал, как правильно скорректировать алгоритм переименования: скажем, переименовываем мы пиксель с изначальной меткой "X", имеющий новую метку "Y" (т.е. assoc[X] != X) в "Z". В этом случае надо наоборот переименовать Z в Y, а величины, больше Z, но <= Ncur декрементировать, затем декрементировать и Ncur. Все довольно просто, а очевидным плюсом является то, что не надо многократно сканировать изображение для переименования "по месту": это будет сделано на шаге 3.
Конечно, все познается в сравнении, поэтому сразу обливать дерьмом этот метод нельзя: стоит попытаться его реализовать. А учитывая то, что все шаги, кроме второго, превосходно параллелизуются (1, 3 - хоть по строкам, хоть по прямоугольным областям, 4 - по номерам объектов), на куде, пожалуй, он будет работать просто превосходно!

еще мысли

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

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