Функциональная обработка изображений в D

в 23:23, , рубрики: Без рубрики

image

Недавно я завершил переработку графического пакета для моей D библиотеки. Вдохновленный модулями std.algorithm и std.range, я пытался достичь следующих целей:

  • Представить все в виде малых комбинируемых компонентов
  • Избежать неявного копирования и предпочтительно использовать ленивые вычисления
  • Использовать шаблоны для производительности

Начиная с первой версии, все компонеты пакета обработки изображений были параметризированы типом цвета. Это не стандартный способ реализации графических библиотек — большинство абстрагируют конкретный тип цвета изображения через ООП интерфейс, или просто конвертируют все изображения в единый формат пикселей, с которыми далее работают в памяти. Однако для большинства случаев это является тратой памяти и времени, обычно разработчики заранее знают в каком конкретно формате будет представлено изображение, за исключением приложений, где графические данные вводятся пользователем (например, граф. редакторы). Вместо этого моя библиотека объявляет все типы изображений как шаблоны с типом-параметром для цвета.

Я весьма доволен результатами работы над библиотекой, поэтому я хочу поделиться несколькими интересными моментами в данном посте.


Библиотека начинается с определения view:

/// Любой тип, у которого есть width, height
/// и может быть проиндексирован для получения
/// цвета конкретного пикселя
enum isView(T) =
    is(typeof(T.init.w) : size_t) && // width
    is(typeof(T.init.h) : size_t) && // height
    is(typeof(T.init[0, 0])     );   // color information

Этот метод объявления статических интерфейсов такой же, который используется в std.range, например, в isInputRange. Вместо того, чтобы объявить интерфейс близкий по смыслу к ООП, статичные интерфейсы в D условно определяются с помощью теста на реализованность определенных фич (прим. пер. утиная типизация). Проверка завершается успешно, если операции, которые должен реализовать тип, компилируются без ошибок или имеют определенный тип. Обычно для этого используется IsExpression или trait compiles.

Аналогично std.range.ElementType определяем шаблон для получения типа, который используется во view для цвета пикселя:

/// Возвращает тип цвета переданного view
alias ViewColor(T) = typeof(T.init[0, 0]);

Далее, мы определяем несколько специализаций для view:

Код

/// Views могут быть неизменяемые
/// или поддерживать запись
enum isWritableView(T) =
    isView!T &&
    is(typeof(T.init[0, 0] = ViewColor!T.init));

/// Опционально view может реализовывать
/// прямой доступ к пикселям. 
/// Такие views называются "direct views"
enum isDirectView(T) =
    isView!T &&
    is(typeof(T.init.scanline(0)) : ViewColor!T[]);

Опять-таки это похоже на определение isForwardRange: проверка на то, что тип реализовывает все базовые фичи, а также некоторые дополнительные возможности, специфичные для этой специализации.

Так как реализация прямого доступа к пикселям может быть определена через scanline direct view, объявляем template mixin, который это реализовывает:

Код

/// Миксин, который реализует примитивы view 
/// на базе существующих примитивов direct view
mixin template DirectView()
{
    alias COLOR = typeof(scanline(0)[0]);
 
    /// Реализация оператора view[x, y]
    ref COLOR opIndex(int x, int y)
    {
        return scanline(y)[x];
    }
 
    /// Реализация оператора view[x, y] = c
    COLOR opIndexAssign(COLOR value, int x, int y)
    {
        return scanline(y)[x] = value;
    }
}

Например, определим шаблон Image, который описывает изображение в памяти:

Код

/// Изображение в памяти
/// Пиксели хранятся в плоском массиве
struct Image(COLOR)
{
    int w, h;
    COLOR[] pixels;
 
    /// Возвращает массив пикселей строчки y
    COLOR[] scanline(int y)
    {
        assert(y>=0 && y<h);
        return pixels[w*y..w*(y+1)];
    }
 
    mixin DirectView;
 
    this(int w, int h)
    {
        size(w, h);
    }
 
    /// Изменение размера без масштабирования
    void size(int w, int h)
    {
        this.w = w;
        this.h = h;
        if (pixels.length < w*h)
            pixels.length = w*h;
    }
}

Массивы в D реализованы через указатель и длину (и запрашивают память только тогда, когда расширяются или склеиваются), поэтому выражение pixels[w*y… w*(y+1)] не создает копию массива.

Юнит тест проверяет во время компиляции, что Image в самом деле удовлетворяет условиям интерфейса isDirectView:

unittest
{
    static assert(isDirectView!(Image!ubyte));
}

Итак, что мы можем cделать с этой моделью?

Сперва мы можем определить изображения, которые на самом деле не ссылаются на массив пикселей, а вместо этого вычисляют их по требованию:

Код

/// Возвращает view, который расчитывает пиксели
/// по требованию на основе переданной формулы
template procedural(alias formula)
{
    alias fun = binaryFun!(formula, "x", "y");
    alias COLOR = typeof(fun(0, 0));
 
    auto procedural(int w, int h)
    {
        struct Procedural
        {
            int w, h;
 
            auto ref COLOR opIndex(int x, int y)
            {
                return fun(x, y);
            }
        }
        return Procedural(w, h);
    }
}

Данный одноименный шаблон использует std.functional.binaryFun для преобразования строки (которая будет примешана) в предикат, или делегат (лямбду). Так как функция имеет возвращаемый тип auto и возвращает struct, объявленный внутри этой функции, Procedural является примером Voldemort types (прим. пер. есть перевод).

Простейший пример процедурного изображения, заполненного одним цветом:

/// Возвращает view, указанных размеров
/// и заполненное одним цветом
auto solid(COLOR)(COLOR c, int w, int h)
{
    return procedural!((x, y) => c)(w, h);
}

Заметьте, как тип цвета возвращаемого view выводится из типа параметра c, поэтому solid(RGB(1, 2, 3), 10, 10) вернет view из RGB пикселей, даже если оно не имеет полностью определенного имени.


Другая вещь, которая может быть представлена в этой модели — создание views, которые трансформируют другие views различными способами. Определим другой template mixin для часто используемого кода:

Код

/// Миксин, который реализует view на базе
/// другого view, используя функцию преобразования
/// координат
mixin template Warp(V)
{
    V src;
 
    auto ref ViewColor!V opIndex(int x, int y)
    {
        warp(x, y);
        return src[x, y];
    }
 
    static if (isWritableView!V)
    ViewColor!V opIndexAssign(ViewColor!V value, int x, int y)
    {
        warp(x, y);
        return src[x, y] = value;
    }
}

Посмотрим на строку static if (isWritableView!V), которая говорит, что оператор view[x, y] = c должен быть определен только, если нижележащий view его поддерживает. Итого wraped view будет изменяемым только, если нижележащий view тоже можно изменять.

Используя эту функцию, мы можем определить cropping view, который представляет прямоугольный кусок другого view:

Код

/// Обрезанный view по заданному прямоугольнику
auto crop(V)(auto ref V src, int x0, int y0, int x1, int y1)
    if (isView!V)
{
    assert( 0 <=    x0 &&  0 <=    y0);
    assert(x0 <     x1 && y0 <     y1);
    assert(x1 <= src.w && y1 <= src.h);
 
    static struct Crop
    {
        mixin Warp!V;
 
        int x0, y0, x1, y1;
 
        @property int w() { return x1-x0; }
        @property int h() { return y1-y0; }
 
        void warp(ref int x, ref int y)
        {
            x += x0;
            y += y0;
        }
 
        static if (isDirectView!V)
        ViewColor!V[] scanline(int y)
        {
            return src.scanline(y0+y)[x0..x1];
        }
    }
 
    static assert(isDirectView!V == isDirectView!Crop);
 
    return Crop(src, x0, y0, x1, y1);
}

if (isView!V) template constraint проверяет, чтобы первый аргумент соответствовал условиям интерфейса isView.

Как и раньше, crop использует isDirectView, чтобы предоставить прямой доступ к пикселям, если нижележащее изображение это поддерживает. Прямой доступ к пикселям полезен при работе с большим количеством пикселей единовременно, что приводит к увеличению производительности по сравнению с последовательным доступом. Например, при копировании одного изображения в другое гораздо быстрее использовать slice copies (типобезопасная замена memcpy в D), чем присваивать значения каждого пикселя независимо:

Код

/// Копирует один view в другой.
/// Views должны иметь одинаковый размер.
void blitTo(SRC, DST)(auto ref SRC src, auto ref DST dst)
    if (isView!SRC && isWritableView!DST)
{
    assert(src.w == dst.w && src.h == dst.h, "View size mismatch");
    foreach (y; 0..src.h)
    {
        static if (isDirectView!SRC && isDirectView!DST)
            dst.scanline(y)[] = src.scanline(y)[];
        else
        {
            foreach (x; 0..src.w)
                dst[x, y] = src[x, y];
        }
    }
}

Та же идея, что и в crop может быть использована, чтобы реализовать view, который тайлит другой view или масштабирует по nearest-neighbor алгоритму (алгоритмы масштабирования посложнее лучше реализуются в императивном стиле). Код очень похож на crop, поэтому я его сюда не включаю.

Даже если crop берет источник как обычный аргумент, предполагаемое использование этой функции и ей подобных такое, как будто она является методом исходного view: someView.nearestNeighbor(100, 100).tile(1000, 1000).crop(50, 50, 950, 950). Эта возможность обусловлена фичей языка, называемой «Uniform Function Call Syntax» (или просто UFCS), которая позволяет писать a.fun(b...) вместо fun(a, b...). Главным плюсом этой фичи является возможность организации chaining (a.fun1().fun2().fun3() вместо fun3(fun2(fun1(a)))), что используется на всю катушку в Phobos и в данном пакете.


Для простых преобразования, в которых не меняется размер view, мы можем определить вспомогательную функцию, которая упрощает применение формулы, заданной пользователем, к каждой координате пикселя:

Код

/// Возвращает view с преобразованными координатами
/// согласно переданным формулам
template warp(string xExpr, string yExpr)
{
    auto warp(V)(auto ref V src)
        if (isView!V)
    {
        static struct Warped
        {
            mixin Warp!V;
 
            @property int w() { return src.w; }
            @property int h() { return src.h; }
 
            void warp(ref int x, ref int y)
            {
                auto nx = mixin(xExpr);
                auto ny = mixin(yExpr);
                x = nx; y = ny;
            }
 
            private void testWarpY()()
            {
                int y;
                y = mixin(yExpr);
            }
 
            /// Если x координата не изменяется и y не зависит
            /// от x, то мы можем преобразовывать целые scanlines.
            static if (xExpr == "x" &&
                __traits(compiles, testWarpY()) &&
                isDirectView!V)
            ViewColor!V[] scanline(int y)
            {
                return src.scanline(mixin(yExpr));
            }
        }
 
        return Warped(src);
    }
}

warp использует хитрый метод проверки передаваемых формул. Функция testWarpY объявлена как шаблон, тем не менее с нулем шаблонных аргументов. Это заставляет компилятор не проводить семантический анализ тела этой функции до момента ее использования. И так как она не имеет x в области видимости, она может быть успешно инстанирована, только если yExpr не использует x. Выражение __traits(compiles, testWrapY()) как раз это проверяет. Все это позволяет нам определить direct view scanline, только если уверены в том, что можем сделать это безопасно. Пример:

Код

/// Возвращает view с инверитрованной координатой x
alias hflip = warp!(q{w-x-1}, q{y});
 
/// Возвращает view с инверитрованной координатой y
alias vflip = warp!(q{x}, q{h-y-1});
 
/// Возвращает view с инвертированными координатами x и y
alias flip = warp!(q{w-x-1}, q{h-y-1});

Синтаксис q{...} является просто удобным способом определять строковые константы. Эта запись обычно используется для D кода, который потом примешивается где-нибудь. Выражение имеет доступ ко всем символам в месте примешивания — в нашем случае, это функция warp и метод testWarpY структуры Wrapped.

Так как vflip удовлетворяет первым двум условиям, необходимым для объявления метода scanline, someView.vflip() будет являться direct view, если таковым является someView. И это было достигнуто без явной проверки для объявления vflip.

Так как используемая абстракция не полагается на динамический полиморфизм, компилятор свободен проводить встраивание вызовов всех слоев преобразований. Инвертирование изображения дважды не порождает операцию, и, в самом деле, i[5, 5] и i.hflip().hflip()[5, 5] генерируют одинаковый машинный код. Компиляторы D с более совершенным бекэндом могут проводоить еще более агрессивные оптимизации: например, если определили функцию flipXY, которая инвертирует оси X и Y, и rotateCW (для поворота изображения на 90° против часовой стрелки) как src.flipXY().hflip(), тогда четыре удачных вызова rotateCW вырезаются при оптимизации.


Давайте перейдем к операциям над самими пикселями. Главная функция в std.algorithmmap, которая возвращает range, лениво применяющий выражение над другим range. Наша colorMap использует эту идею для цветов:

Код

/// Возвращает view, который применяет предикат
/// к каждому пикселю нижележащего view. 
template colorMap(alias pred)
{
    alias fun = unaryFun!(pred, false, "c");
 
    auto colorMap(V)(auto ref V src)
        if (isView!V)
    {
        alias OLDCOLOR = ViewColor!V;
        alias NEWCOLOR = typeof(fun(OLDCOLOR.init));
 
        struct Map
        {
            V src;
 
            @property int w() { return src.w; }
            @property int h() { return src.h; }
 
            auto ref NEWCOLOR opIndex(int x, int y)
            {
                return fun(src[x, y]);
            }
        }
 
        return Map(src);
    }
}

Используя colorMap, определить функцию, которая инвертирует цвета изображения, также просто как:

alias invert = colorMap!q{~c};

colorMap не требует совпадения типов цвета источника и результата. Это позволяет использовать ее для преобразования цвета: read(«image.bmp»).parseBMP!RGB().colorMap!(c => BGRX(c.b, c.g, c.r)) вернет RGB битмап как BGRX view.


Обработка изображений очень часто хорошо поддается параллелизации. std.parallelism помогает сделать задачу параллельной обработки изображений тривиальной:

Код

/// Разделяет view на сегменты и
/// применяет fun к каждому сегменту
/// параллельно.
/// Возвращает массив сегментов,
/// который может быть объединен,
/// используя vjoin или vjoiner.
template parallel(alias fun)
{
    auto parallel(V)(auto ref V src, size_t chunkSize = 0)
        if (isView!V)
    {
        auto processSegment(R)(R rows)
        {
            auto y0 = rows[0];
            auto y1 = y0 + rows.length;
            auto segment = src.crop(0, y0, src.w, y1);
            return fun(segment);
        }
 
        import std.range : iota, chunks;
        import std.parallelism : taskPool, parallel;
 
        if (!chunkSize)
            chunkSize = taskPool.defaultWorkUnitSize(src.h);
 
        auto range = src.h.iota.chunks(chunkSize);
        alias Result = typeof(processSegment(range.front));
        auto result = new Result[range.length];
        foreach (n; range.length.iota.parallel(1))
            result[n] = processSegment(range[n]);
        return result;
    }
}

Даже если parallel разделяет свое имя с теми функциями, что присутствуют в std.parallelism, конфликта не возникает, так как они имеют различные сигнатуры и работают на различных типах.

Вместе с этим, операцию можно разделить между несколькими потоками заменой image.process() на image.parallel!(segment => segment.process()).vjoin().


Практические примеры:

  • Создание короткой анимации:
    Код

    import std.algorithm;
    import std.parallelism;
    import std.range;
    import std.stdio;
    import std.string;
     
    import ae.utils.geometry;
    import ae.utils.graphics.color;
    import ae.utils.graphics.draw;
    import ae.utils.graphics.gamma;
    import ae.utils.graphics.image;
     
    void main()
    {
        enum W = 4096;
     
        const FG = L16(0);
        const BG = L16(ushort.max);
     
        auto image = Image!L16(W, W);
        image.fill(BG);
     
        enum OUTER = W/2 * 16/16;
        enum INNER = W/2 * 13/16;
        enum THICK = W/2 *  3/16;
     
        image.fillCircle(W/2, W/2, OUTER, FG);
        image.fillCircle(W/2, W/2, INNER, BG);
        image.fillRect(0, W/2-INNER, W/2, W/2+INNER, BG);
        image.fillRect(W/2-THICK/2, W/2-INNER, W/2+THICK/2, W/2+INNER, FG);
     
        enum frames = 32;
        foreach (n; frames.iota.parallel)
            image
            .rotate(TAU * n / frames, BG)
            .copy
            .downscale!(W/16)
            .lum2pix(gammaRamp!(ushort, ubyte, ColorSpace.sRGB))
            .toPNG
            .toFile("loading-%02d.png".format(n++));
    }
    

    Эта программа рисует изначальное изображение в большем разрешении, используя 16-битнyю яркость, которе потом преобразовывает в 8-битное sRGB изображение после уменьшения. Уменьшение масштаба позволяет избежать aliasing, и цветовое преобразование необходимо для точного изменения размера.

    Вот что получилось после конвертации PNG кадров в GIF:
    image

  • Отрисовка серо-белого изображения множества Мандельброта:
    Код

    auto mandelbrot(int w, int h)
    {
        import std.algorithm, std.range;
        import ae.utils.graphics.view;
        return procedural!((x, y)
        {
            auto c = (2.0*x/w - 1.5) + (2.0*y/h - 1.0)*1i;
            return cast(ubyte)(1+
                recurrence!((a, n) => c + a[n-1]^^2)(0+0i)
                .take(ubyte.max)
                .countUntil!(z => z.re^^2 + z.im^^2 > 4));
        })(w, h);
    }
     
    void main()
    {
        import std.stdio;
        import ae.utils.graphics.image;
        mandelbrot(500, 500).toPNG().toFile("mandel.png");
    }
    

    Эта программа использует recurrence, take and countUntil примитивы над D ranges, также как и встроенную поддержку комплексных чисел, позволяя лаконично реализовать алгоритм, который обычно требует множество строк для реализации (Однако встроенные комплексные числа сейчас в процессе устаревания в пользу std.complex).

    Результат работы:

    Картинка

    image


Подход с использованием шаблонов обещает большой выигрыш в производительности. Как простой бенчмарк, данная программа уменьшает масштаб всех изображений в директории на 25%:

Код

import std.file;
import std.path;
import std.stdio;
 
import ae.utils.graphics.color;
import ae.utils.graphics.gamma;
import ae.utils.graphics.image;
 
void main()
{
    alias BGR16 = Color!(ushort, "b", "g", "r");
    auto gamma = GammaRamp!(ushort, ubyte)(2.2);
    foreach (de; dirEntries("input", "*.bmp", SpanMode.shallow))
    {
        static Image!BGR scratch1;
        static Image!BGR16 scratch2, scratch3;
 
        de
        .read
        .parseBMP!BGR(scratch1)
        .parallel!(segment =>
            segment
            .pix2lum(gamma)
            .copy(scratch2)
            .downscale!4(scratch3)
            .lum2pix(gamma)
            .copy
        )(4)
        .vjoin
        .toBMP
        .toFile(buildPath("output-d", de.baseName));
    }
}

Я сравнил результат работы с эквивалентной командой ImageMagick:

convert           
 input/*.bmp      
 -depth 16        
 -gamma 0.454545  
 -filter box      
 -resize 25%      
 -gamma 2.2       
 -depth 8         
 output-im/%02d.bmp

Версия на D выполняется в 4-5 разра быстрее. Конечно, это нечестное сравнение: даже если обе используют 16-битную глубину цвета, гамма коррекцию, многопоточность и оптимизированы под одинаковые архитектуры, D программа содержит код специально оптимизированный для этой задачи. Если не брать в расчет различные JIT техники, пакет нельзя сравнивать с библиотеками обработки изображений общего назначения.


Графический пакет доступен на GitHub. Спасибо David Ellsworth за его вклад в эту статью.

Автор: NCrashed

Источник

Поделиться

* - обязательные к заполнению поля