Этот пост, по-хорошему, заслуживает место во flame.comp, т.к. побудительным мотивом изначально являлся вопрос alex_public, заданный в недрах флейма в том форуме:
Но могу для разнообразия подкинуть ещё один примерчик для понимания. Простейшая, но очень нужная во многих областях задачка: имеем двухмерный массив чисел и надо получить из него новый массив, в котором каждая точка является средним от своих соседей (4-ёх или 8-ми, не суть). Что об этом скажет linq? )))
Но, поскольку речь идёт об пределах возмжностей кода на C#, я всё же запощу это сюда — в интересах тех, кого тоже волнует вопрос о производительности "математического" кода на дотнете, а перечитывать топики глубиной за сотню постов терпения нет.
Итак, что об этом скажет Linq? Вопрос, понятное дело, с подковыркой — типа, да Linq же — это убогая надстройка над IEnumerable, которую потом и соплями прикрутили к RDBMS. А тут у нас — принципиально другой материал, неперечислимый и, в отличие от, скажем, графов, к перечислениям несводимый.
Лично моими устами об этом Linq скажет следующее:
int[,] FourNeighborAverage(int[,] data) => from d in data select (d[-1, 0] + d[1, 0] + d[0, 1] + d[0, -1]) / 4;
Как это работает?
Ну, для начала у нас есть public static класс Array2d, где определён вот такой вот extension метод:
public static T[,] Select<T>(this T[,] source, Kernel<T> kernel)
То есть выражение d => (d[-1, 0] + d[1, 0] + d[0, 1] + d[0, -1]) / 4 трактуется как Kernel<int>:
public delegate T Kernel<T>(IRelativeAccess2d<T> point);
В свою очередь, интерфейс IRelativeAccess2d<T> определён весьма незатейливо:
public interface IRelativeAccess2d<T>
{
T this[int dx, int dy] { get; }
}
Теперь всё становится понятным: в наш Kernel передаётся интерфейс для доступа к "соседним" относительно текущего элементам массива, что позволяет нам описывать выражения типа "среднее от своих четырёх соседей".
Реализовать наш метод Select не очень трудно. Наивная реализация выглядит так:
public static T[,] Select<T>(this T[,] source, Kernel<T> kernel)
{
var h = source.GetLength(0);
var w = source.GetLength(1);
var result = new T[h, w];
for (int i = 0; i < h; i++)
for (int j = 0; j < w; j++)
result[i, j] = kernel(new RelativeArrayAccess2d<T>(source, i, j));
return result;
}
структура RelativeArrayAccess2d<T> слишком тривиальна, чтобы её приводить — она хранит ссылку на массив source и координаты текущего элемента; при обращении по this[dx, dy] мы прибавляем относительные координаты к текущим и возвращаем элемент.
Этого кода уже достаточно для того, чтобы реализовать копирование массива в стиле linq:
from d in data select d[0, 0];
К сожалению, первый же запуск с ядром усреднения вылетит с ArrayIndexOutOfBoundsException. Логично — мы пробуем обращаться за пределы исходного массива. Суровые программисты на плюсах на такие мелочи внимания не обращают — ну, или требуют от программиста ручного указания размеров "каёмки", в которой фильтр не работает.
Правильное решение этой задачи — довольно сложно: нам придётся учиться подпиливать фильтр на ходу таким образом, чтобы он корректно работал на границах массива (у элемента рядом со "стороной" прямоугольника не 4, а 3 соседа, а у углового элемента — только 2. Строго говоря, именно их и нужно усреднять).
Но пока что мы обойдёмся решением в стиле С++ — ограничением области вывода фильтра.
Сделать это не очень сложно — мы просто один раз вызовем kernel со специальной реализацией IRelativeAccess, которая запоминает, куда идёт обращение:
private sealed class KernelMeasure<T> : IRelativeAccess2d<T>
{
...
T IRelativeAccess2d<T>.this[int x, int y]
{
get
{
xmin = Math.Min(xmin, x);
xmax = Math.Max(xmax, x);
ymin = Math.Min(ymin, y);
ymax = Math.Max(ymax, y);
return default(T);
}
}
}
Вообще говоря, это будет работать не для всех типов Kernel — например, мы могли использовать Math.Random() для выбора смещения. Но для практически интересных случаев, когда смещения всегда постоянны, этот подход сработает.
Итого, у нас получается более-менее работоспособный код, который позволяет удобно записывать формулы фильтрации:
public static T[,] Select<T>(this T[,] source, Kernel<T> kernel)
{
var h = source.GetLength(0);
var w = source.GetLength(1);
var km = new KernelMeasure<T>();
kernel(km); // провели измеренияvar t = 0 - km.xmin;
var b = h - km.xmax;
var l = 0 - km.ymin;
var r = w - km.ymax;
var result = new T[h, w];
for (int i = t; i < b; i++)
for (int j = l; j < r; j++)
result[i, j] = kernel(new RelativeArrayAccess2d<T>(source, i, j));
return result;
}
На этом можно было бы и остановиться, если бы не ехидный комментарий в том же флейме на 30 уровней выше
Ну как же, если мы говорим про работу с коллекциями (а не про СУБД), то обычный императивный код (типа for() ... ) никуда не делся. И при этом он по прежнему является как наиболее обобщённым, так и наиболее производительным.
CHALLENGE ACCEPTED!
По поводу "обобщённости" спорить смысла нет — очевидно, что приведённый from d in data select... можно без изменений заставить работать не только с 2d-массивом, но и с произвольной структурой, которая представляет из себя двумерную сетку, индексированную целочисленными координатами. Например, мы можем работать с "изображением" из double с размерами миллион на миллион, то есть имея 10^12 точек или около 7 терабайт. Перебор их циклом for упрётся даже не во время, а в размер памяти и диска локальной машины. Так что мы пилим всё на тайлы, раскидываем по нескольким тысячам серверов, отправляем код для исполнения, дожидаемся, мёрджим. Ну, то есть получаем опять хэндл чудовищного "холста" для дальнейших операций — например, edge-detect, а потом ещё свёртка, а потом понижение размерности, и уже в конце получаем что-нибудь приемлемое для памяти локальной машины. И всё — в компактном и удобном для чтения (и написания) виде.
Вопрос у нас к тому пункту, где про "наиболее производительным".
Давайте напишем тот идеальный императивный код, который бы написал на C# программист, не доверяющий всем этми новомодным линкам, с их "динамикой и рефлексией":
static int[,] PerfectFourNeighborAverage(this int[,] source)
{
var h = source.GetLength(0);
var w = source.GetLength(1);
var result = new int[h, w];
for (int i = 1; i < h-1; i++)
for (int j = 1; j < w-1; j++)
result[i, j] = (source[i, j-1] + source[i, j+1] + source[i-1, j] + source[i+1,j])/4;
return result;
}
Возьмём случайный массив размером 1000*1000, и попробуем сравнить производительность нашего "обобщённого" кода со "специфичным":
BenchmarkDotNet=v0.10.14, OS=Windows 10.0.15063.1088 (1703/CreatorsUpdate/Redstone2)
Intel Core i7-6600U CPU 2.60GHz (Skylake), 1 CPU, 4 logical and 2 physical cores
Frequency=2742190 Hz, Resolution=364.6720 ns, Timer=TSC
[Host] : .NET Framework 4.6.1 (CLR 4.0.30319.42000), 32bit LegacyJIT-v4.7.2650.0
DefaultJob : .NET Framework 4.6.1 (CLR 4.0.30319.42000), 32bit LegacyJIT-v4.7.2650.0
Method | Mean | Error | StdDev | Scaled | ScaledSD |
----------- |----------:|----------:|----------:|-------:|---------:|
Perfect | 15.008 ms | 0.1706 ms | 0.1512 ms | 1.00 | 0.00 |
Linq | 33.423 ms | 0.6599 ms | 0.8345 ms | 2.23 | 0.06 |
Ну, на первый взгляд, linq справился не так уж и плохо — всего лишь в 2.23 раза медленнее.
Тем не менее, эта незначительная цена за удобную абстракцию легко может отпугнуть фанатов тактовыжимания.
Поэтому в следующей серии я попробую показать, что можно с этим сделать, ничего не меняя в коде бенчмарка. То есть те самые преимущества, которые бесплатны для прикладного программиста — от него потребуется разве что перекомпилировать свой проект. Никаких хинтов в коде — ни императивного развёртывания циклов, ни ручной векторизации, ни даже умелой расстановки "декларативных" прагм.
Только хардкор, только голый from d in data select (d[-1, 0] + d[1, 0] + d[0, 1] + d[0, -1]) / 4.
Уйдемте отсюда, Румата! У вас слишком богатые погреба.
Если начинаешь заниматься битовыжиманием, то гоняться нужно вот с этой штукой http://halide-lang.org/
Всё остальное не серьёзно.
В принципе эту штуку можно и на C# натянуть с почти тем же синтаксисом.
И вычислительная модель тоже позволяет раскидать вычисления по кластеру.
... << RSDN@Home 1.0.0 alpha 5 rev. 0>>
Пусть это будет просто:
просто, как только можно,
но не проще.
(C) А. Эйнштейн
Здравствуйте, Sinclair, Вы писали:
S>надо получить из него новый массив, в котором каждая точка является средним от своих соседей (4-ёх или 8-ми, не суть). Что об этом скажет linq? )))
Зачем одевать носок через голову?
S>Поэтому в следующей серии я попробую показать, что можно с этим сделать
Надеюсь, что тогда станет понятно зачем. Прямо сейчас я в сомнениях.
Возможно говоря linq вы понимаете работу с IEnumerable<T>, который не выгодно загружать целиком в виде массива? Можно использовать тип данных, поддерживающий виртуализацию, c public indexer и каким-нибудь кешем внутри, в этом случае память не переполнится и PerfectFourNeighborAverage() станет нормальным решением.
Здравствуйте, Sinatr, Вы писали: S>Зачем одевать носок через голову?
Ок. Ваши предложения для решения задачи (на C#)? S>>Поэтому в следующей серии я попробую показать, что можно с этим сделать S>Надеюсь, что тогда станет понятно зачем. Прямо сейчас я в сомнениях.
Пока что — задача уйти от кода PerfectFourNeighborAverage(), в котором слишком много лишних деталей, которых не было в формулировке исходной задачи.
Это делает его (очевидно) плохо читаемым (например, легко ли обнаружить one-off ошибки?) и (спорно) плохо оптимизируемым. S>Возможно говоря linq вы понимаете работу с IEnumerable<T>, который не выгодно загружать целиком в виде массива?
int[,] не реализует IEnumerable<T>. S>Можно использовать тип данных, поддерживающий виртуализацию, c public indexer и каким-нибудь кешем внутри, в этом случае память не переполнится и PerfectFourNeighborAverage() станет нормальным решением.
Очень вряд ли PerfectFourNeighborAverage() станет нормальным решением для семитерабайтного массива, даже если использовать внутри какой-нибудь кеш.
Уйдемте отсюда, Румата! У вас слишком богатые погреба.
Здравствуйте, Sinclair, Вы писали:
S>>Зачем одевать носок через голову? S>Ок. Ваши предложения для решения задачи (на C#)?
PerfectFourNeighborAverage() или аналогичное решение, которое не использует linq от слова вообще.
S>Это делает его (очевидно) плохо читаемым (например, легко ли обнаружить one-off ошибки?) и (спорно) плохо оптимизируемым.
Любые достаточно сложные алгоритмы плохочитаемые и это нормально. Использовать linq, чтобы легче понять алгоритм или не использовать, чтобы он был оптимальным?
S>Очень вряд ли PerfectFourNeighborAverage() станет нормальным решением для семитерабайтного массива, даже если использовать внутри какой-нибудь кеш.
Вы имели в виду размер базы данных или именно "массив"? Возможно вы меня не поняли, я предлагал использовать другой тип данных + ваш метод с использованием indexer (indexer != array). А если речь именно о размере, то неважно какой алгоритм, все равно придется все данные обработать, все 7 терабайт.
Здравствуйте, Sinatr, Вы писали:
S>PerfectFourNeighborAverage() или аналогичное решение, которое не использует linq от слова вообще.
Это — худший возможный вариант. Его дорого как писать, так и исполнять.
S>Любые достаточно сложные алгоритмы плохочитаемые и это нормально. Использовать linq, чтобы легче понять алгоритм или не использовать, чтобы он был оптимальным?
В первую очередь алгоритм должен быть понятным. Потому что понятный алгоритм можно сделать корректным. Нас принципиально не интересуют алгоритмы, которые быстро выдают неверный результат. S>>Очень вряд ли PerfectFourNeighborAverage() станет нормальным решением для семитерабайтного массива, даже если использовать внутри какой-нибудь кеш. S>Вы имели в виду размер базы данных или именно "массив"?
Я имел в виду "массив данных". S>Возможно вы меня не поняли, я предлагал использовать другой тип данных + ваш метод с использованием indexer (indexer != array). А если речь именно о размере, то неважно какой алгоритм, все равно придется все данные обработать, все 7 терабайт.
Ну, тут как раз важно, какой алгоритм. Например, можно ли обрабатывать данные в несколько потоков. Или, лучше, на нескольких машинах.
Можно ли векторизовать алгоритм.
Алгоритм с явным выписыванием for for крайне трудно поддаётся каким-либо трансформациям — нам нужен очень умный компилятор, чтобы восстановить замысел программиста, и реализовать его по-новому, а не так, как написано.
Уйдемте отсюда, Румата! У вас слишком богатые погреба.
Ты потратил (скажи сам , сколько) времени, на то, чтобы доказать, что этот простейший фильтр реализовать с помощью LinQ таки можно, написав в итоге 2 страницы текста с делегатом, интерфейсом, kernel и т.д. там, где в императивном стиле пишется простейший двойной цикл с одним оператором внутри. Пишется "на автомате", почти не думая, так как подобное любой опытный программист писал много раз.
А результат — просадка скорости в 2 с лишним раза.
И это все при том, что стоит серьезно поменять задачу, и ты будешь опять придумывать решение, и хорошо, если придумаешь.
Здравствуйте, Pavel Dvorkin, Вы писали:
PD>И это все при том, что стоит серьезно поменять задачу, и ты будешь опять придумывать решение, и хорошо, если придумаешь.
Всё как раз наоборот. Стоит поменять задачу — скажем, увеличить радиус гаусс-блюра в фильтре — и код, написанный программистом "почти не думая", окажется не у дел. А если он вдруг потратит время на оптимизацию этого "недумательного" кода, то всё станет ещё хуже: там же объём кода уже раз в 10 вырос — все вот эти loop unrolling, SIMD intrincics, это же всё не просто так. Теперь понять, куда воткнуть дополнительные коэффициенты, крайне тяжело.
А ещё нужно не забыть поменять пределы в циклах, и т.д. и т.п. PD>Зачем ?
Затем, чтобы отучить людей бояться performance penalty.
Уйдемте отсюда, Румата! У вас слишком богатые погреба.
S>Всё как раз наоборот. Стоит поменять задачу — скажем, увеличить радиус гаусс-блюра в фильтре — и код, написанный программистом "почти не думая", окажется не у дел.
Я под "поменять задачу" несколько иное имел в виду. Ты вторгся в матричные операции, а фильтры такого рода — это простейшее из того, что там есть. Для обработки изображения, может, этого фильтра и хватит, а вот для чего-то большего ?
Вот сейчас я со своими дипломниками как раз делал задачу, в которой хоть и есть простенький фильтр, но есть еще
Все это 99%, а простенький фильтр — 1%. Это тоже предложишь на LinQ переводить ? А если нет — зачем мне 1% кода переводить на LinQ ?
>А если он вдруг потратит время на оптимизацию этого "недумательного" кода, то всё станет ещё хуже: там же объём кода уже раз в 10 вырос — все вот эти loop unrolling, SIMD intrincics, это же всё не просто так. Теперь понять, куда воткнуть дополнительные коэффициенты, крайне тяжело.
Эээ, не передергивай. Он как раз потратил время на оптимизацию, потому что в алгоритмах обработки изображений O(N^2) — это очень хорошо, а зачастую бывает существенно хуже. А ты этой оптимизацией и не занимался, да и не ясно, сможешь ли вообще ей заняться.
PD>>Зачем ? S>Затем, чтобы отучить людей бояться performance penalty.
В целом же скажу вот что. LinQ заточен под последовательную обработку данных. Ты сумел показать, что и этот фильтр можно подвести под механизм последовательной обработки, пусть понимаемой несколько расширенно. Однако когда дело дойдет до тех алгоритмов, где обработка ведется не последовательным, а более хитрым способом, ничего путного у тебя не получится
Это как с последовательным доступом к файлам vs. прямым доступом. Если надо писать или читать в/из поток(а) — классы последовательного доступа замечательно работают. Если же надо файл модифицировать сложным образом — применять классы последовательного доступа, наверное, все же можно , но лучше не стоит.
Здравствуйте, Pavel Dvorkin, Вы писали: PD>Я под "поменять задачу" несколько иное имел в виду. Ты вторгся в матричные операции, а фильтры такого рода — это простейшее из того, что там есть. Для обработки изображения, может, этого фильтра и хватит, а вот для чего-то большего ?
А для большего — придумаем своё решение. Linq не позицинируется как замена вообще всему.
Но ему можно найти достаточно неожиданные применения, что я и хочу продемонстрировать.
Не нужно сводить Linq к "медленному способу обработки перечислений".
PD>Вот сейчас я со своими дипломниками как раз делал задачу, в которой хоть и есть простенький фильтр, но есть еще
PD>Бинаризация по Sauvola (алгоритм со скользящим окном с целью как раз уменьшить число операций) PD>http://www.ee.oulu.fi/mvg/files/pdf/pdf_24.pdf
Сходу тяжко врубиться в 20 страниц научного текста. Как выглядит код бинаризации по Sauvola? Какая часть алгоритма является "переменной" (зависит от специфики конкретной задачи), а какая — общей?
PD>Нахождение линий по Хафу PD>https://ru.wikipedia.org/wiki/%D0%9F%D1%80%D0%B5%D0%BE%D0%B1%D1%80%D0%B0%D0%B7%D0%BE%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5_%D0%A5%D0%B0%D1%84%D0%B0
Опять же — интересно было бы посмотреть на то, как выглядит типичный код на С++, который ищет какие-то примитивы при помощи преобразования Хафа или PD>и кое-что еще
PD>Все это 99%, а простенький фильтр — 1%. Это тоже предложишь на LinQ переводить ? А если нет — зачем мне 1% кода переводить на LinQ ?
Нет. Я предложу писать код в максимально читаемой форме, где сказано "что", но не сказано "как".
PD>Эээ, не передергивай. Он как раз потратил время на оптимизацию, потому что в алгоритмах обработки изображений O(N^2) — это очень хорошо, а зачастую бывает существенно хуже. А ты этой оптимизацией и не занимался, да и не ясно, сможешь ли вообще ей заняться.
Я-то как раз смогу. Причём эта оптимизация пригодится не только в этой строчке этой программы, а во всех местах, где у меня вычисляется фильтр по изображению.
PD>И научить вот такому PD>http://rsdn.org/forum/flame.comp/4098976.1
PD>В целом же скажу вот что. LinQ заточен под последовательную обработку данных.
Неа. Как раз про последовательность в Linq нету совсем-совсем ничего.
PD>Ты сумел показать, что и этот фильтр можно подвести под механизм последовательной обработки, пусть понимаемой несколько расширенно. Однако когда дело дойдет до тех алгоритмов, где обработка ведется не последовательным, а более хитрым способом, ничего путного у тебя не получится
Ждём пример такого алгоритма. Как обычно — интересует алгоритм, в котором есть "переменная часть". В случае с фильтром у нас было задание ядра; в случае с СУБД — задание предиката для поиска.
PD>Это как с последовательным доступом к файлам vs. прямым доступом. Если надо писать или читать в/из поток(а) — классы последовательного доступа замечательно работают. Если же надо файл модифицировать сложным образом — применять классы последовательного доступа, наверное, все же можно , но лучше не стоит.
Вот как раз "императивный for" — это алгоритм последовательного доступа. А Linq — он весь про "дай мне результат". Что потенциально позволяет выбирать "последовательность" так, как нам это удобно.
Уйдемте отсюда, Румата! У вас слишком богатые погреба.
Здравствуйте, Sinclair, Вы писали:
S>Сходу тяжко врубиться в 20 страниц научного текста. Как выглядит код бинаризации по Sauvola? Какая часть алгоритма является "переменной" (зависит от специфики конкретной задачи), а какая — общей?
Держи
void binarize(int image_width, int image_height, unsigned char *bin_image, unsigned char *gray_image,
int w, double k){
int whalf = w>>1;
// Calculate the integral image, and integral of the squared imagedouble **integral_image,**rowsum_image;
double **integral_sqimg,**rowsum_sqimg;
unsigned char** gray_image_ptr = new unsigned char*[image_height];
gray_image_ptr[0] = gray_image;
for(int i = 1; i < image_height; i++)
gray_image_ptr[i] = gray_image_ptr[i-1] + image_width;
unsigned char** bin_image_ptr = new unsigned char*[image_height];
bin_image_ptr[0] = bin_image;
for(int i = 1; i < image_height; i++)
bin_image_ptr[i] = bin_image_ptr[i-1] + image_width;
integral_image = new double*[image_width];
for(int i = 0; i < image_width; i++)
integral_image[i] = new double[image_height];
rowsum_image = new double*[image_width];
for(int i = 0; i < image_width; i++)
rowsum_image[i] = new double[image_height];
integral_sqimg = new double*[image_width];
for(int i = 0; i < image_width; i++)
integral_sqimg[i] = new double[image_height];
rowsum_sqimg = new double*[image_width];
for(int i = 0; i < image_width; i++)
rowsum_sqimg[i] = new double[image_height];
int xmin,ymin,xmax,ymax;
double diagsum,idiagsum,diff,sqdiagsum,sqidiagsum,sqdiff,area;
double mean,std,threshold;
// START for(int j=0; j<image_height; j++){
rowsum_image[0][j] = gray_image_ptr[j][0];
rowsum_sqimg[0][j] = gray_image_ptr[j][0]*gray_image_ptr[j][0];
}
for(int i=1; i<image_width; i++){
for(int j=0; j<image_height; j++){
rowsum_image[i][j] = rowsum_image[i-1][j] + gray_image_ptr[j][i];
rowsum_sqimg[i][j] = rowsum_sqimg[i-1][j] + gray_image_ptr[j][i]*gray_image_ptr[j][i];
}
}
for(int i=0; i<image_width; i++){
integral_image[i][0] = rowsum_image[i][0];
integral_sqimg[i][0] = rowsum_sqimg[i][0];
}
for(int i=0; i<image_width; i++){
for(int j=1; j<image_height; j++){
integral_image[i][j] = integral_image[i][j-1] + rowsum_image[i][j];
integral_sqimg[i][j] = integral_sqimg[i][j-1] + rowsum_sqimg[i][j];
}
}
// END
//Calculate the mean and standard deviation using the integral imagefor(int i=0; i<image_width; i++){
for(int j=0; j<image_height; j++){
xmin = __max(0,i-whalf);
ymin = __max(0,j-whalf);
xmax = __min(image_width-1,i+whalf);
ymax = __min(image_height-1,j+whalf);
area = (xmax-xmin+1)*(ymax-ymin+1);
// area can't be 0 here
// proof (assuming whalf >= 0):
// we'll prove that (xmax-xmin+1) > 0,
// (ymax-ymin+1) is analogous
// It's the same as to prove: xmax >= xmin
// image_width - 1 >= 0 since image_width > i >= 0
// i + whalf >= 0 since i >= 0, whalf >= 0
// i + whalf >= i - whalf since whalf >= 0
// image_width - 1 >= i - whalf since image_width > i
// --IM
// ASSERT(area);if(!xmin && !ymin){ // Point at origin
diff = integral_image[xmax][ymax];
sqdiff = integral_sqimg[xmax][ymax];
}
else if(!xmin && ymin){ // first column
diff = integral_image[xmax][ymax] - integral_image[xmax][ymin-1];
sqdiff = integral_sqimg[xmax][ymax] - integral_sqimg[xmax][ymin-1];
}
else if(xmin && !ymin){ // first row
diff = integral_image[xmax][ymax] - integral_image[xmin-1][ymax];
sqdiff = integral_sqimg[xmax][ymax] - integral_sqimg[xmin-1][ymax];
}
else{ // rest of the image
diagsum = integral_image[xmax][ymax] + integral_image[xmin-1][ymin-1];
idiagsum = integral_image[xmax][ymin-1] + integral_image[xmin-1][ymax];
diff = diagsum - idiagsum;
sqdiagsum = integral_sqimg[xmax][ymax] + integral_sqimg[xmin-1][ymin-1];
sqidiagsum = integral_sqimg[xmax][ymin-1] + integral_sqimg[xmin-1][ymax];
sqdiff = sqdiagsum - sqidiagsum;
}
mean = diff/area;
std = sqrt((sqdiff - diff*diff/area)/(area-1));
threshold = mean*(1+k*((std/128)-1));
if(gray_image_ptr[j][i] <= threshold)
bin_image_ptr[j][i] = 0;
else
bin_image_ptr[j][i] = 255;
}
}
delete[] gray_image_ptr;
delete[] bin_image_ptr;
for(int i = 0; i < image_width; i++)
delete[] integral_image[i];
delete[] integral_image;
for(int i = 0; i < image_width; i++)
delete[] rowsum_image[i];
delete[] rowsum_image;
for(int i = 0; i < image_width; i++)
delete[] integral_sqimg[i];
delete[] integral_sqimg;
for(int i = 0; i < image_width; i++)
delete[] rowsum_sqimg[i];
delete []rowsum_sqimg;
}
Сама идея достаточно проста — вычисляются суммы и суммы квадратов для каждого пикселя , считая его центром окна размером в w (обычно 15-25), потом простенькая формула
Проблема в том, что если это делать в лоб, то будет O(N*M*w^2)
Мы его сами не писали, но вот тебе код, и даже на C#
PD>>Все это 99%, а простенький фильтр — 1%. Это тоже предложишь на LinQ переводить ? А если нет — зачем мне 1% кода переводить на LinQ ? S>Нет. Я предложу писать код в максимально читаемой форме, где сказано "что", но не сказано "как".
Ну напишем 1% кода "в максимально читаемой форме", а 99% останутся в нечитаемой.
PD>>Эээ, не передергивай. Он как раз потратил время на оптимизацию, потому что в алгоритмах обработки изображений O(N^2) — это очень хорошо, а зачастую бывает существенно хуже. А ты этой оптимизацией и не занимался, да и не ясно, сможешь ли вообще ей заняться. S>Я-то как раз смогу. Причём эта оптимизация пригодится не только в этой строчке этой программы, а во всех местах, где у меня вычисляется фильтр по изображению.
PD>>Ты сумел показать, что и этот фильтр можно подвести под механизм последовательной обработки, пусть понимаемой несколько расширенно. Однако когда дело дойдет до тех алгоритмов, где обработка ведется не последовательным, а более хитрым способом, ничего путного у тебя не получится S>Ждём пример такого алгоритма. Как обычно — интересует алгоритм, в котором есть "переменная часть". В случае с фильтром у нас было задание ядра; в случае с СУБД — задание предиката для поиска.
Ну вот тебе Савола. В общем-то, тот же фильтр , так что бери его в качестве переменной части. Правда, окно там не размера 2, а обычно 15-25.
Ну и для Хафа тоже есть переменная часть. Например, можно окружности искать
PD>>Это как с последовательным доступом к файлам vs. прямым доступом. Если надо писать или читать в/из поток(а) — классы последовательного доступа замечательно работают. Если же надо файл модифицировать сложным образом — применять классы последовательного доступа, наверное, все же можно , но лучше не стоит. S>Вот как раз "императивный for" — это алгоритм последовательного доступа.
Вообще-то, кроме for, есть еще и while
Кстати, как насчет binary search с помощью LinQ ? Разумеется, без написания цикла.
Здравствуйте, Pavel Dvorkin, Вы писали:
PD>Здравствуйте, Sinclair, Вы писали:
S>>Сходу тяжко врубиться в 20 страниц научного текста. Как выглядит код бинаризации по Sauvola? Какая часть алгоритма является "переменной" (зависит от специфики конкретной задачи), а какая — общей?
PD>Держи
PD>[ccode] PD>void binarize(int image_width, int image_height, unsigned char *bin_image, unsigned char *gray_image, PD> int w, double k){ PD> int whalf = w>>1;
Я правильно понимаю, что в integral_image[i][j] в итоге лежит сумма всех пикселов "левее и выше" текущего пиксела + он сам?
То есть integral_image[0][0] == gray_image_ptr[0][0]; а integral_image[image_width-1][image_height-1] == сумме всех пикселов исходного изображения?
Уйдемте отсюда, Румата! У вас слишком богатые погреба.
Здравствуйте, Pavel Dvorkin, Вы писали:
PD>Здравствуйте, Sinclair, Вы писали:
S>>Сходу тяжко врубиться в 20 страниц научного текста. Как выглядит код бинаризации по Sauvola? Какая часть алгоритма является "переменной" (зависит от специфики конкретной задачи), а какая — общей?
PD>Держи
PD>
PD>void binarize(int image_width, int image_height, unsigned char *bin_image, unsigned char *gray_image,
PD> int w, double k){
PD> int whalf = w>>1;
PD> // Calculate the integral image, and integral of the squared image
PD> double **integral_image,**rowsum_image;
PD> double **integral_sqimg,**rowsum_sqimg;
PD> unsigned char** gray_image_ptr = new unsigned char*[image_height];
PD> gray_image_ptr[0] = gray_image;
PD> for(int i = 1; i < image_height; i++)
PD> gray_image_ptr[i] = gray_image_ptr[i-1] + image_width;
PD> unsigned char** bin_image_ptr = new unsigned char*[image_height];
PD> bin_image_ptr[0] = bin_image;
PD> for(int i = 1; i < image_height; i++)
PD> bin_image_ptr[i] = bin_image_ptr[i-1] + image_width;
PD> integral_image = new double*[image_width];
PD> for(int i = 0; i < image_width; i++)
PD> integral_image[i] = new double[image_height];
PD> rowsum_image = new double*[image_width];
PD> for(int i = 0; i < image_width; i++)
PD> rowsum_image[i] = new double[image_height];
PD> integral_sqimg = new double*[image_width];
PD> for(int i = 0; i < image_width; i++)
PD> integral_sqimg[i] = new double[image_height];
PD> rowsum_sqimg = new double*[image_width];
PD> for(int i = 0; i < image_width; i++)
PD> rowsum_sqimg[i] = new double[image_height];
PD> int xmin,ymin,xmax,ymax;
PD> double diagsum,idiagsum,diff,sqdiagsum,sqidiagsum,sqdiff,area;
PD> double mean,std,threshold;
PD> // START
PD> for(int j=0; j<image_height; j++){
PD> rowsum_image[0][j] = gray_image_ptr[j][0];
PD> rowsum_sqimg[0][j] = gray_image_ptr[j][0]*gray_image_ptr[j][0];
PD> }
PD> for(int i=1; i<image_width; i++){
PD> for(int j=0; j<image_height; j++){
PD> rowsum_image[i][j] = rowsum_image[i-1][j] + gray_image_ptr[j][i];
PD> rowsum_sqimg[i][j] = rowsum_sqimg[i-1][j] + gray_image_ptr[j][i]*gray_image_ptr[j][i];
PD> }
PD> }
PD> for(int i=0; i<image_width; i++){
PD> integral_image[i][0] = rowsum_image[i][0];
PD> integral_sqimg[i][0] = rowsum_sqimg[i][0];
PD> }
PD> for(int i=0; i<image_width; i++){
PD> for(int j=1; j<image_height; j++){
PD> integral_image[i][j] = integral_image[i][j-1] + rowsum_image[i][j];
PD> integral_sqimg[i][j] = integral_sqimg[i][j-1] + rowsum_sqimg[i][j];
PD> }
PD> }
PD> // END
PD> //Calculate the mean and standard deviation using the integral image
PD> for(int i=0; i<image_width; i++){
PD> for(int j=0; j<image_height; j++){
PD> xmin = __max(0,i-whalf);
PD> ymin = __max(0,j-whalf);
PD> xmax = __min(image_width-1,i+whalf);
PD> ymax = __min(image_height-1,j+whalf);
PD> area = (xmax-xmin+1)*(ymax-ymin+1);
PD> // area can't be 0 here
PD> // proof (assuming whalf >= 0):
PD> // we'll prove that (xmax-xmin+1) > 0,
PD> // (ymax-ymin+1) is analogous
PD> // It's the same as to prove: xmax >= xmin
PD> // image_width - 1 >= 0 since image_width > i >= 0
PD> // i + whalf >= 0 since i >= 0, whalf >= 0
PD> // i + whalf >= i - whalf since whalf >= 0
PD> // image_width - 1 >= i - whalf since image_width > i
PD> // --IM
PD> // ASSERT(area);
PD> if(!xmin && !ymin){ // Point at origin
PD> diff = integral_image[xmax][ymax];
PD> sqdiff = integral_sqimg[xmax][ymax];
PD> }
PD> else if(!xmin && ymin){ // first column
PD> diff = integral_image[xmax][ymax] - integral_image[xmax][ymin-1];
PD> sqdiff = integral_sqimg[xmax][ymax] - integral_sqimg[xmax][ymin-1];
PD> }
PD> else if(xmin && !ymin){ // first row
PD> diff = integral_image[xmax][ymax] - integral_image[xmin-1][ymax];
PD> sqdiff = integral_sqimg[xmax][ymax] - integral_sqimg[xmin-1][ymax];
PD> }
PD> else{ // rest of the image
PD> diagsum = integral_image[xmax][ymax] + integral_image[xmin-1][ymin-1];
PD> idiagsum = integral_image[xmax][ymin-1] + integral_image[xmin-1][ymax];
PD> diff = diagsum - idiagsum;
PD> sqdiagsum = integral_sqimg[xmax][ymax] + integral_sqimg[xmin-1][ymin-1];
PD> sqidiagsum = integral_sqimg[xmax][ymin-1] + integral_sqimg[xmin-1][ymax];
PD> sqdiff = sqdiagsum - sqidiagsum;
PD> }
PD> mean = diff/area;
PD> std = sqrt((sqdiff - diff*diff/area)/(area-1));
PD> threshold = mean*(1+k*((std/128)-1));
PD> if(gray_image_ptr[j][i] <= threshold)
PD> bin_image_ptr[j][i] = 0;
PD> else
PD> bin_image_ptr[j][i] = 255;
PD> }
PD> }
PD> delete[] gray_image_ptr;
PD> delete[] bin_image_ptr;
PD> for(int i = 0; i < image_width; i++)
PD> delete[] integral_image[i];
PD> delete[] integral_image;
PD> for(int i = 0; i < image_width; i++)
PD> delete[] rowsum_image[i];
PD> delete[] rowsum_image;
PD> for(int i = 0; i < image_width; i++)
PD> delete[] integral_sqimg[i];
PD> delete[] integral_sqimg;
PD> for(int i = 0; i < image_width; i++)
PD> delete[] rowsum_sqimg[i];
PD> delete []rowsum_sqimg;
PD> }
PD>
Уфф. Короче, на первый взгляд, эквивалентом этой порнографии на C# будет примерно вот такое:
public static byte[,] Binarize(byte[,] grayImage, int w, double k)
{
int image_height = grayImage.GetLength(0);
int image_width = grayImage.GetLength(1);
var integral_image = (double[,])Array.CreateInstance(typeof(double), new int[] { image_height + 1, image_width + 1 }, new int[] { -1, -1 });
var integral_sqimg = (double[,])Array.CreateInstance(typeof(double), new int[] { image_height + 1, image_width + 1 }, new int[] { -1, -1 });
var result = new byte[image_height, image_width];
for (int i = 0; i < image_height; i++)
for (int j = 0; j < image_width; j++)
{
integral_image[i, j] = grayImage[i, j] + integral_image[i - 1, j] + integral_image[i, j - 1] - integral_image[i - 1, j - 1]; ;
integral_sqimg[i, j] = grayImage[i, j] * grayImage[i, j] + integral_sqimg[i - 1, j] + integral_sqimg[i, j - 1] - integral_sqimg[i - 1, j - 1];
}
var whalf = w / 2;
for (int i = 0; i < image_height; i++)
for (int j = 0; j < image_width; j++)
{
var xmin = Math.Max(0, i - whalf) - 1;
var ymin = Math.Max(0, j - whalf) - 1;
var xmax = Math.Min(image_width - 1, i + whalf);
var ymax = Math.Min(image_height - 1, j + whalf);
var area = (xmax - xmin) * (ymax - ymin);
var diff = integral_image[xmax, ymax] + integral_image[xmin, ymin] - integral_image[xmax, ymin] - integral_image[xmin, ymax];
var sqdiff = integral_sqimg[xmax, ymax] + integral_sqimg[xmin, ymin] - integral_sqimg[xmax, ymin] - integral_sqimg[xmin, ymax];
var mean = diff / area;
var std = Math.Sqrt((sqdiff - diff * mean) / (area - 1));
var threshold = mean * (1 + k * ((std / 128) - 1));
result[i, j] = grayImage[i, j] > threshold ? (byte)255 : (byte)0;
}
return result;
}
Здесь ровно два прохода по N*M вместо трёх в вашем коде, памяти выделяется в два раза меньше (экономией на наркомании с gray_image_ptr, который должен быть функтором, а не массивом, можно пренебречь), ветвлений (и возможностей напороть с ошибками) значительно меньше.
В общем, пока что просто исправлены очевидные косяки исходного С++ кода.
Теперь можно подумать и о том, как сделать этот С#-код понятнее
Уйдемте отсюда, Румата! У вас слишком богатые погреба.
Здравствуйте, Sinclair, Вы писали:
S>Уфф. Короче, на первый взгляд, эквивалентом этой порнографии на C# будет примерно вот такое:
Насчет порнографии — к Sauvola. Это тот код, который я у него взял, я лишь поменял тот класс Matrix, который у него был, на простой массив и расставил указатели.
Код правилен или нет — судить не берусь, надо тестировать.
S>Здесь ровно два прохода по N*M вместо трёх в вашем коде,
В моем коде, сделанном из этого, не только число проходов изменено, но и многое другое, в частности, код сделан многопоточным, потому что основное требование было — скорость.
S>Теперь можно подумать и о том, как сделать этот С#-код понятнее
Здравствуйте, Sinclair, Вы писали:
S>Павел, огромное спасибо за красивую задачу!
S>Ответ на неё я дам в части 4 своего повествования.
S>Скептики будут посрамлены!
Успехов! Говорю без всякой иронии — мне самому будет интересно посмотреть.
Ну и конечно, обязательное условие — не увеличивать объем требуемой памяти.
А все же, как насчет бинарного поиска с использованием LinQ ? Дело, конечно, не в нем, просто много алгоритмов на векторах и матрицах, в которых проход отнюдь не последовательный, даже в твоем расширенном понимании. Если не хочешь бинарный поиск, а тебе нужно, чтобы обязательно было что-то вроде
Здравствуйте, Pavel Dvorkin, Вы писали: PD>Насчет порнографии — к Sauvola. Это тот код, который я у него взял, я лишь поменял тот класс Matrix, который у него был, на простой массив и расставил указатели.
Похоже, что Sauvola хороший математик, но плохой программист
PD>Код правилен или нет — судить не берусь, надо тестировать.
Скорее всего правилен, т.к. там меньше места, где сделать ошибку. Плюс-минус опечатки вроде j вместо i. PD>В моем коде, сделанном из этого, не только число проходов изменено, но и многое другое, в частности, код сделан многопоточным, потому что основное требование было — скорость.
А можно, чисто для ориентирования, прикинуть, сколько строк в вашей финальной версии метода binarize()? S>>Теперь можно подумать и о том, как сделать этот С#-код понятнее
Уйдемте отсюда, Румата! У вас слишком богатые погреба.
Здравствуйте, Sinclair, Вы писали:
S>Похоже, что Sauvola хороший математик, но плохой программист
Совершенно верно. Точнее, он математик, который запрограммировал свой алгоритм.
PD>>Код правилен или нет — судить не берусь, надо тестировать. S>Скорее всего правилен, т.к. там меньше места, где сделать ошибку. Плюс-минус опечатки вроде j вместо i.
Протестируй Берешь gray-scaled image и вперед.
S>А можно, чисто для ориентирования, прикинуть, сколько строк в вашей финальной версии метода binarize()?
Думаю, раза в 3 примерно больше. Во-первых, потоки (на чистом WinAPI, с ивентами) , во-вторых, я учел тот факт, что будет массовая обработка, поэтому все эти массивы не выделяются каждый раз, а вместо этого VirtualAlloc одного большого блока, и если в следующий раз больше не нужно, то он переиспользуется — это тоже дало некоторый выигрыш. Ну и ряд других хитростей.
Но в конце концов все это переписали под CUDA. Вот там я и занимался примерно тем же, чем ты сейчас — придумывал реализации для такого алгоритма, который обычно пишется не думая.
Здравствуйте, Pavel Dvorkin, Вы писали:
PD>Протестируй Берешь gray-scaled image и вперед.
Для этого нужно иметь образец, и результат его применения.
Рандомное-то изображение этот код прожёвывает на ура — но без образца непонятно, правильный ли получается вариант.
PD>Думаю, раза в 3 примерно больше. Во-первых, потоки (на чистом WinAPI, с ивентами) , во-вторых, я учел тот факт, что будет массовая обработка, поэтому все эти массивы не выделяются каждый раз, а вместо этого VirtualAlloc одного большого блока, и если в следующий раз больше не нужно, то он переиспользуется — это тоже дало некоторый выигрыш. Ну и ряд других хитростей.
Ужас какой. PD>Но в конце концов все это переписали под CUDA. Вот там я и занимался примерно тем же, чем ты сейчас — придумывал реализации для такого алгоритма, который обычно пишется не думая.
В райском мире светлого будущего переписывать прикладную программу не понадобится.
Уйдемте отсюда, Румата! У вас слишком богатые погреба.
Здравствуйте, Sinclair, Вы писали:
PD>>Протестируй Берешь gray-scaled image и вперед. S>Для этого нужно иметь образец, и результат его применения. S>Рандомное-то изображение этот код прожёвывает на ура — но без образца непонятно, правильный ли получается вариант.
Попробовал залить на RSDN — 500. Видимо, размер его не устраивает — оригинальный (8bpp BMP) — 21 Мб, сжатый до JPG — 6.5 Мб. Примерно 5000 * 4000 пикселей.
Время обработки на моем стареньком Phenom II X4 955.
Оригинальный Sauvola (тот, что дал тебе) 2.5 сек
Оптимизированный и с потоками — 200-250 мсек.
Уменьшать не буду — неинтересно.
Как тебе его передать ?
PD>>Думаю, раза в 3 примерно больше. Во-первых, потоки (на чистом WinAPI, с ивентами) , во-вторых, я учел тот факт, что будет массовая обработка, поэтому все эти массивы не выделяются каждый раз, а вместо этого VirtualAlloc одного большого блока, и если в следующий раз больше не нужно, то он переиспользуется — это тоже дало некоторый выигрыш. Ну и ряд других хитростей. S>Ужас какой.
Может, конечно, и ужас, но когда стоит задача выжать по скорости все, что можно (потому что если будет недостаточно быстро,то это все равно, как если бы и вообще не было — пойдешь на любой ужас).
PD>>Но в конце концов все это переписали под CUDA. Вот там я и занимался примерно тем же, чем ты сейчас — придумывал реализации для такого алгоритма, который обычно пишется не думая. S>В райском мире светлого будущего переписывать прикладную программу не понадобится.
Может быть, но только не раньше чем в райском мире светлого будущего.
А пока что, когда я начал разбираться с CUDA , то я сразу понял одно — почти все, что я знаю о программировании, здесь не работает. Дело в том, что там аппаратные потоки, и написать надо так, чтобы они могли работать параллельно. А там очень жесткие (по крайней мере тогда) были требования к тому, чтобы работало параллельно. Шаг вправо, шаг влево — и все работать будет, но только последовательно, и получим хороший penalty. Вот мы его и изгоняли. В итоге это чудо занимало примерно 500 строк, зато время обработки было 50 мсек на дешевой NVidia карте того времени (2008 год).
, но лучше не стоит.
