Деформация цветовых пространств методом Moving Least Squares

— Меньше оваций, больше ассигнаций!
(из фильма «O Brother, Where Art Thou?»)

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

Сегодня предлагаю отдохнуть от чисто программистского трюкачества и заняться настоящими хардкорным факин-балавством, но уже не в терминах простых и понятных объектов, а будем искривлять и портить то, что так любим: пространства.  Цветовые. Например, так:

Слева вверху Lab - плоскость a-b, слева внизу HSV - плоскость H-V

Видео для привлечения внимания (дальше следует безумный лонгрид — убегайте отсюда пока живы!):

Еще картинка для привлечения внимания. Зачем нам портить цветовые пространства, например. Ну, на мой вкус — это более натуральный способ цветокоррекции, скажем так.

mls-example

Приложение для экспериментов

Выдаем сразу, что-бы не пугать слабых духом. Как всегда можно свободно скачать для сборки из проекта ImageMetalling.

Вспомнить всё

Начнем с легкой разминки перед подходами. Если вы впервые на этом ресурсе и не читали некоторые из постов и, например, не совсем в теме, то храните про запас ссылки на поясняющие посты: рассказ про LUT-ы; геометрия искажений; и еще про геометрию; как получить пиксел в GPU для использования его в расчетах дальше; как показать цветовое пространство в виде 3D фигурки. Все это нужно для понимания, вспомогательного, но избыточного кода приложения этого поста.

Основаная идея

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

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

А теперь представьте, все это можно заменить одним единственным и единственно правильным механизмом? Единственным, потому что он заменяет все инструменты разом. Правильным — в силу перцептуального воздействия в соответствии с природой восприятия человека, а не прямолинейно в лоб, как это решает известный всем Великий и Могучий Фоторедактор.

Вот сегодня про такой механизм и поговорим. И поговорим как можно искривлять не только само изображение, но представление цвета в нем. Простым механическим способом доступным любому креативному стилисту прошедшему курс компьютерной грамотности.

Деформация объектов методом MLS

В нескольких постах ранее были описаны методы деформации изображений в терминах аффинных преобразований и в терминах нелинейной трансформации на B-Spline. Это прекрасные, фундаментальные инструменты. Основная проблема с ними — либо ограниченное количество контрольных точек, которыми мы можем управлять, либо слишком плавная экстраполяция целевых искажений. С цветовыми пространствами нам может понадобится, как можно более жесткая связь контрольной точки и её окружения с её деформацией. Более того, набор точек может быть произвольным. Такая задача, достаточно хорошо решена для анимационных задач в компьютерной графике, но почему-то до сих пор неактивно используется в инструментах больших фото-редакторов, да и вообще не сильно используется как инструмент коррекции цвета. В общем совершим акт гражданского самосожжения — подарим миру новый взгляд на космос.

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

Решение задачи

И так представим, что у нас есть N-контрольных точек плоскости, для пущей умозрительности можете представить себе плоскость с разбросанными шариками на абсолютно не прогибаемой резиновой поверхности, но в плоскости поверхности эта резина свободно может растягиваться (лежит на твердом столе, к примеру). За каждый шарик можно тянуть и фиксировать в целевой точке. Теперь представьте, что это не просто поверхность, а раскрашена цветом какого-то среза какого-то Color Space, например RGB, или HSV, или для наглядности и плавности путь будет Lab.

Lab
L=80

Тогда в терминах аффинных трансформаций мы должны найти решение для уравнения:

{F = {\sum_{i=0}^{N-1}w_i{|f_x(p_i)-q_i|}^2 \rightarrow 0}}  (1)

где: {w_i=\frac{1}{{|p_i-x|}^{2\alpha}}} — среднеквадратичный вес, {f_x} — искомая функция трансформации точки {x}{p_i} — исходное положение контрольных точек, {q_i} — целевое положение контрольных точек.

Т.е. как обычно, в каждой непонятной ситуации — минимизируй это! В нашем случае нам нужно найти такую функцию {f_x}, при которой сумма всех разностей точек {q_i} — т.е. то куда мы пришли и значения функции от исходной точки {f_x(p_i)}, т.е. откуда условно тянем — было минимально, или в терминах математики стремилась к нулю.

С другой стороны любое аффинное преобразование можно записать через уравнение:

{f(x) = xM + T}  (2)

Где {M} — матрица трансформации (почти как CCM) и  {T} — матрица трансляции, хотя какая к черту разница. Можно предположить, что и в этом случае для каждой переменной из (1) мы можем найти минимальное значение для (2).

Так как для (2)  производные по каждой из переменных в {f(x)} равны нулю, мы можем решить (2) непосредственно для {T} через матрицу {M}. Частные производные по свободным переменным в {T} создают линейную систему уравнений.

Теперь, как видно, мы можем выразить {T} через {M}:

{{T = f(x) - xM} \Rightarrow {T = q_* - p_*M}}

Где {p_*} и {q_*} — взвешенные центроиды:

{p_* = \frac{\sum_{i=0}^{N-1}w_ip_i}{\sum_{i=0}^{N-1}w_i}}

{q_* = \frac{\sum_{i=0}^{N-1}w_iq_i}{\sum_{i=0}^{N-1}w_i}}

Что дает нам возможность записать (2) как{f_x} от {M}:

{{f(x) = xM +{q_* - p_*M}} \Rightarrow {f(x) = (x - p_*)M + q_* }}  (3)

При таком раскладе, понятно, что нам все равно для каких точек искать минимум уравнения (1). Т.е. искать для точек {f_x(p_i)} или расстояний между взвешенными центроидами и контрольными точками. Поэтому мы будем решать (1) не для самих контрольных точек, а для новых: {\hat{p_i} = p_i-p_*} и {\hat{q_i} = q_i-q_*}. Тогда решение задачи, с учетом (3) сводится к нахождению {M} из :

{F = {\sum_{i=0}^{N-1}w_i{|\hat{p_i}M-\hat{q_i}|}^2 \rightarrow 0}}   (4)

Ура! дальше все просто, задача сводится к дифференцированию (4) и нахождению {M} как решения системы линейных уравнений через обратные матрицы полученной системы. Для того, что бы мы могли найти обратную матрицу, нам нужно свести систему к решению с квадратными матрицами в обоих частях уравнения. И так из (4) предварительно получаем:

{F' = {\sum_{i=0}^{N-1}w_i{|\hat{p_i}M-\hat{q_i}|} = 0 \Rightarrow \sum_{i=0}^{N-1}\hat{p}^Tw_i{\hat{p_i}M-\sum_{i=0}^{N-1}w_i\hat{p}^T\hat{q_i}}=0}}

Бум!:

{M= {\left(\begin{array}{c}{\sum_{i=0}^{N-1}w_i\hat{p_i}^T{\hat{p_i}}}\end{array}\right)}^{-1}\sum_{i=0}^{N-1}w_i\hat{p_i}^T\hat{q_i}}   (5)

Всё, искомая функция найдена, используем (5) для подстановки в (3):

{f(x) = q_* + (x-p_*){\left(\begin{array}{c}{\sum_{i=0}^{N-1}w_i\hat{p_i}^T{\hat{p_i}}}\end{array}\right)}^{-1}\sum_{i=0}^{N-1}w_i\hat{p_i}^T\hat{q_i}}    (6)

В итоге нам остается запрограммировать решение (6), что мы и сделали в в примере ImageMetalling-16. Поскольку солвер системы уравнений использовать не нужно, а матрички всего лишь 2×2, то, теоретически весь механизм должен быть в меру шустрым.

Расширение условий оптимизации

Если присмотреться к {M}, то мы обнаружим, что есть некоторые особенности с ней связанные, поскольку мы решали задачу как задачу аффинного преобразования. Т.е. задачу: масштабирования и сдвига. Понятно, что в реальном мире деформации не ограничиваются только этими двумя историями, поэтому {M}, может быть выражена более сложным вариантом. Например: преобразованием подобия, которые допускают неравномерность трансляции исходного поля. И так представим {M} как {M = |M'M''|} при этом  {|M'M''| = \lambda|M_1M_2|}, где {M'}{M''} связываются через коэффициент подобия  {\lambda} с другими векторами {M_1} и {M_2}. Поскольку все {M_{xx}} подобны другу другу как вектора плоскости, то можно записать равенство подобия как уравнение для двумерного случая: {M_1^TM_1 = M_2^TM_2 = {\lambda}^2 \Rightarrow M_1^TM_2 = 0}. Или переписать как: {M_2 = M_1^{\bot}}. Из чего следует, что (1) мы теперь можем записать как:

{\sum_{i}^{N-1}w_i|\left(\begin{array}{c}{\hat{p_i}}\\ {-\hat{p_i}^{\bot}}\end{array}\right)M_1-\hat{q_i}^T|^{2}\rightarrow 0}   (7)

Поскольку {M} мы выразили через {M_1} и {M_2} оба их нашли, то (5) можем переписать решив (7), как сделали на предыдущем шаге:

{M = \frac{1}{\mu_s}\sum_{i}^{N-1}w_i\left(\begin{array}{c}{\hat{p_i}}\\ {-\hat{p_i}^{\bot}}\end{array}\right)M_1(\hat{q_i}^T-\hat{q_i}^{\bot^T})}   (8)

Где {\mu_s=\sum_{i=0}^{N-1}w_i\hat{p_i}\hat{p_i}^T}.

Наконец, в исходной публикации (а тут можно сказать её вольное переложение) показано как связать преобразование подобия с ковариационной матрицей и найти {M} для так называемого rigid-преобразования, позволяющего деформировать свойства пространства как можно более жестко. На первый взгляд, для графики — это выглядит достаточно натурально, что можно предположить, хорошо скажется на деформации некоторых цветовых пространств для наших целей. Но это еще не точно и только предстоит проверить на практике, что вы можете сделать скачав пример и покрутив деформации в разных пространствах и с разными настройками поиска {M}.

 

Реализация

Основной код солвера напишем на C++ с учетом того, что код будет исполняться как на GPU в ядрах Metal, так и на CPU, для этого его придется обернуть в Objective-C++ расширение. Очевидно, что запускать код солвера на GPU станет выгодно при большом размере массива контрольных точек и пространства, которое мы будем модифицировать.

kernel void kernel_mlsSolver(
                             constant float2  *input_points  [[buffer(0)]],
                             device float2    *output_points [[buffer(1)]],
                             constant float2  *p      [[buffer(2)]],
                             constant float2  *q      [[buffer(3)]],
                             constant int    &count   [[buffer(4)]],
                             constant MLSSolverKind  &kind [[buffer(5)]],
                             constant float    &alpha [[buffer(6)]],
                             uint gid [[thread_position_in_grid]]
                             ){
    ///
    /// Ядро солвера расчитываем новые координаты для каждой точки   
    ///
    float2 point = input_points[gid];
    IMPMLSSolver solver = IMPMLSSolver(point,p,q,count,kind,alpha);
    output_points[gid] = solver.value(point);
}
class IMPMLSSolver{
    
public:
        
    /**
     Создать солвер для МНК
     
     @param point текущая точка относительно которой рассчитываем деформацию 
     @param p массив исходный контрольных точек
     @param q массив контрольных целевых точек
     @param count количество контрольных точек
     @param kind тип солвера: афинное преобразование, преобразования подобия, преобразование с максимально возможной жосткостью 
     @param alpha - степень деформации
     */
    IMPMLSSolver(const float2 point, 
#ifndef __METAL_VERSION__              
              const float2 *p, 
              const float2 *q,
#else
//
// Важно понимать, что это не просто другой тип данных, а другой механизм расположения 
// данных в памяти. В конкретном решении для ядна на Metal мы передаем массивы в 
// константную память GPU.                  
//                 
              constant float2 *p, 
              constant float2 *q,
#endif
              const int count, 
              const MLSSolverKind kind = MLSSolverKindAffine, 
              const float alpha = 1.0): 
    point_(point),
    kind_(kind), 
    alpha_(alpha), 
    count_(count),
    weight_(0),
    p_(p), q_(q)
    {
        
#ifndef __METAL_VERSION__
        //
        // Для CPU сборки выделяем временные буферы для расчета промежуточных значений
        // солвера
        //
        w_ = new float[count_];   
        pHat_ = new float2[count_];   
        qHat_ = new float2[count_];   
#endif
        
        // Находим веса  
        solveW();
        
        // Находим вектора с "шапкой"
        solveHat();
        
        // Решаем основное оптимизационное уравнение
        solveM();  
        
#ifndef __METAL_VERSION__
        delete pHat_;
        delete qHat_;
        pHat_=qHat_=0;
#endif
        
    }    
        
    /**
     Возвращает по координате точку деформации 

     @param point исходная позиция
     @return новая позиция 
     */
    float2 value(float2 point) {
        if (count_ <= 0) return point;   
        return (point - pStar_) * M + qStar_;
    }
    
    ~IMPMLSSolver() {
#ifndef __METAL_VERSION__
        delete w_;
        if (pHat_) delete pHat_;
        if (qHat_) delete qHat_;
#endif
    }
    
private:
    
    MLSSolverKind kind_;
    float   alpha_;
    int     count_;
    float2  point_;
#ifndef __METAL_VERSION__              
    const float2 *p_;
    const float2 *q_;
    float *w_;
    float2 *pHat_;
    float2 *qHat_;    
#else
    constant float2 *p_;
    constant float2 *q_;
    thread float w_[MLS_MAXIMUM_POINTS];
    thread float2 pHat_[MLS_MAXIMUM_POINTS];
    thread float2 qHat_[MLS_MAXIMUM_POINTS];
#endif
    
    float weight_;
    float2 pStar_;
    float2 qStar_;
    float mu_;
    
    float2x2 M;
    
    void solveW() {
        
        weight_ = 0;
        pStar_ = float2(0);
        qStar_ = float2(0);

        for (int i=0; i<count_ && i<MLS_MAXIMUM_POINTS; i++) {
            
            float d =  pow(distance(p_[i], point_), 2*alpha_);
            
            if (d < FLT_EPSILON)  d = FLT_EPSILON; 
            
            w_[i] = 1.0 / d;
            weight_ = weight_ + w_[i];
            
            pStar_ += w_[i] * p_[i];
            qStar_ += w_[i] * q_[i];

        }
        
        if (weight_ < FLT_EPSILON)  weight_ = FLT_EPSILON;
        
        pStar_ = pStar_ / weight_;                
        qStar_ = qStar_ / weight_;
    }  
    
    void solveHat(){        
        
        mu_ = 0;
        
        float _rmu1 = 0;
        float _rmu2 = 0;
        
        for (int i=0; i<count_ && i<MLS_MAXIMUM_POINTS; i++) {
            
            pHat_[i] = p_[i] - pStar_;                        
            qHat_[i] = q_[i] - qStar_;
            
            switch (kind_) {            
                case MLSSolverKindSimilarity:
                    mu_ += similarityMu(i);
                    break;                
                case MLSSolverKindRigid:
                    _rmu1 += rigidMu1(i); 
                    _rmu2 += rigidMu2(i);
                    break;
                default:
                    break;
            }
        }
        
        switch (kind_) {            
            case MLSSolverKindRigid:
                mu_ = sqrt(_rmu1*_rmu1 + _rmu2*_rmu2);
                break;
            default:
                break;
        }
        
        if (mu_ < FLT_EPSILON)  mu_ = FLT_EPSILON; 
        
        mu_ = 1/mu_;
    }
    
    void solveM() {
        switch (kind_) {
            case MLSSolverKindAffine:
                M = affineM();
                break;
            case MLSSolverKindSimilarity:
            case MLSSolverKindRigid:
                M = similarityM(point_);
                break;
        }
    }  
                   
    float2x2 affineM() {
        float2x2 mi = (float2x2){(float2){0,0},(float2){0,0}};
        float2x2 mj = (float2x2){(float2){0,0},(float2){0,0}};

        for (int i=0; i<count_ && i<MLS_MAXIMUM_POINTS; i++) {
            
            float2x2 pti({
                pHat_[i], 
                float2(0)                
            });
            
            float2x2 ppi({
                (float2){w_[i] * pHat_[i].x, 0.0}, 
                (float2){w_[i] * pHat_[i].y, 0.0}                
            });
            
            mi = mi + (float2x2)(pti * ppi);
            
            
            float2x2 ptj({
                w_[i] * pHat_[i], 
                float2(0)                
            });
            
            float2x2 qpj({
                (float2){qHat_[i].x, 0.0}, 
                (float2){qHat_[i].y, 0.0}                
            });
            
            mj = mj + (float2x2)(ptj * qpj);
        }        
        return inverse(mi) * mj;
    }
    
    
    float similarityMu(int i)  {
        return w_[i]*dot(pHat_[i], pHat_[i]);
    }
    
    float2x2 similarityM(float2 value) {
        
        float2x2 m = (float2x2){(float2){0,0},(float2){0,0}};
        
        for (int i=0; i<count_ && i<MLS_MAXIMUM_POINTS; i++) {
            
            float2 sp = -1 * __slashReflect(pHat_[i]);
            float2 sq = -1 * __slashReflect(qHat_[i]);
            
            float2x2 _p({
                (float2){w_[i] * pHat_[i].x, w_[i] * sp.x}, 
                (float2){w_[i] * pHat_[i].y, w_[i] * sp.y}                
            });
            
            float2x2 _q({qHat_[i], sq});
                        
            m = m + (float2x2)(_p * _q);
        }
        return  mu_ * m; 
    }
    
    float rigidMu1(int i) {
        return w_[i]*dot(qHat_[i], pHat_[i]);        
    }
    
    float rigidMu2(int i) {
        return w_[i]*dot(qHat_[i],  __slashReflect(pHat_[i]));        
    }
};

 

Примеры

В приложении можно тыкнуть на цвета картинки, цвета отметятся на плоскости цветового пространства Lab: a-b (в коде можете его заменить на любое другое). Цвета, которые хочется закрепить и не изменять также стоит отметить. А затем потаскать за точку-цвет, который хочется изменить на другой. На этом примере, я выделил цвет бутылки и хаотично потыкал по другим. Как видно сильно изменяя один цвет, остальные практически не затрагиваются, что может быть весьма удобно для всяких коммерческих фотографов занимающихся предметной съемкой на цифру, или для творческих личностей вечно мечтающих подогнать цвет цифры к плёнке.

mls-bottle-diff

А вот так выглядит RGB кубик после наших издевательств над Lab-пространством. Для референса аппликуха рисует еще hald-LUT, и HV плоскость из HSV пространства.

mls-bottle-lut

На этом примере проведена легкая цветокоррекция на вкус автора. Убран холодный оттенок исходного слайда.

mls-correction-diff

Вот так вот примерно трансформируется LUT исходного RGB-кубика.

mls-correction-lut

 

Выводы

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

Всем пис! И да, надеюсь, что этот пост не взорвал ваши лобные доли.

 


 

Еще раз напомню: авторы блога не преследуют задачи быть предельно корректным, но если заметили явную ашипку, если написали явную глупость, если что-то не понятно: комментируйте или пишите на: imagemetalling [*] gmail.com .

Реклама

Добавить комментарий

Заполните поля или щелкните по значку, чтобы оставить свой комментарий:

Логотип WordPress.com

Для комментария используется ваша учётная запись WordPress.com. Выход /  Изменить )

Google+ photo

Для комментария используется ваша учётная запись Google+. Выход /  Изменить )

Фотография Twitter

Для комментария используется ваша учётная запись Twitter. Выход /  Изменить )

Фотография Facebook

Для комментария используется ваша учётная запись Facebook. Выход /  Изменить )

Connecting to %s