Image Abstraction の実装

/***********************************************************
* Image Abstraction Implementation
************************************************************
* このコードは画像を抽象化して、イラスト風の画像を生成します。
* 画像の抽象化にはBilateralフィルタを、エッジの抽出には
* Difference-of-Gaussianというモデルを用いています。
*
* 参考文献
* H.Winnemoller et al., "Real-Time Video Abstraction",
* ACM TOG 2005.
*
* 使い方
* usage: > ImageAbstraction.exe [input image]
*
* Copyright:
* This program is coded by tatsy. You can use this code
* for any purpose including commercial usage :-)
************************************************************/

#include <iostream>
#include <vector>
#include <cmath>
using namespace std;

#include <opencv2/opencv.hpp>

const int ksize = 5;
const double sigma_d = 3.0;
const double sigma_r = 4.25;
const double sigma_e2 = 0.5 * 0.5;
const int ne = 2;
const int nb = 4;

const double phi_e = 3.0;
const double phi_q = 3.0;
const double tau = 0.98;
const int q = 10;

void bfilter(cv::Mat& in, cv::Mat& out);

int main(int argc, char** argv) {
    if(argc <= 1) {
        cout << "usage: > ImageAbstraction.exe [input image]" << endl;
        return -1;
    }

    cv::Mat img = cv::imread(argv[1], CV_LOAD_IMAGE_COLOR);
    if(img.empty()) {
        cout << "Failed to load "" << argv[1] << ""." << endl;
        return -1;
    }

    int width = img.cols;
    int height = img.rows;
    img.convertTo(img, CV_32FC3, 1.0 / 255.0);

    // 色空間の変換
    cv::Mat lab;
    cv::cvtColor(img, lab, CV_BGR2Lab);

    // フィルタ処理
    cv::Mat out, tmp;
    lab.convertTo(out, CV_32FC3);
    for(int i=0; i<ne; i++) {
        out.convertTo(tmp, CV_32FC3);
        cv::bilateralFilter(tmp, out, ksize, sigma_r, sigma_d);
    }

    // DoGによるエッジ抽出
    cv::Mat gray;
    cv::cvtColor(out, gray, CV_Lab2BGR);
    cv::cvtColor(gray, gray, CV_BGR2GRAY);
    cv::Mat edge = cv::Mat(height, width, CV_32FC1);
    for(int y=0; y<height; y++) {
        for(int x=0; x<width; x++) {
            double Se = 0.0;
            double Sr = 0.0;
            double weight_e = 0.0;
            double weight_r = 0.0;
            for(int dy=-ksize; dy<=ksize; dy++) {
                for(int dx=-ksize; dx<=ksize; dx++) {
                    int xx = x + dx;
                    int yy = y + dy;
                    if(xx < 0 || yy < 0 || xx>=width || yy>=height) continue;

                    double pn = gray.at<float>(yy, xx) * 255.0;
                    double d2 = dx * dx + dy * dy;
                    double we = exp(-0.5 * d2 / sigma_e2);
                    double wr = exp(-0.5 * d2 / (2.56 * sigma_e2));
                    weight_e += we;
                    weight_r += wr;
                    Se += pn * we;
                    Sr += pn * wr;
                }
            }

            Se = Se / weight_e;
            Sr = Sr / weight_r;
            double diff = Se - tau * Sr;
            if(diff > 0.0) {
                edge.at<float>(y, x) = 1.0f;
            } else {
                edge.at<float>(y, x) = 1.0f + (float)tanh(phi_e * diff);
            }
        }
    }

    // 再度フィルタ処理
    for(int i=0; i<ne; i++) {
        out.convertTo(tmp, CV_32FC3);
        cv::bilateralFilter(tmp, out, ksize, sigma_r, sigma_d);
    }

    // 輝度の量子化
    cv::cvtColor(out, out, CV_Lab2BGR);
    double dq = 1.0 / (double)q;
    for(int c=0; c<3; c++) {
        vector<int> hist(q, 0);
        for(int y=0; y<height; y++) {
            for(int x=0; x<width; x++) {
                double val = out.at<float>(y, x*3+c);
                int bin = (int)(val / dq + 0.5);
                double qnear = bin * dq;
                out.at<float>(y, x*3+c) = (float)(qnear + dq / 2.0 * tanh(phi_q * (val - qnear)));
            }
        }
    }

    // 抽象画像とエッジを組み合わせる  
    for(int y=0; y<height; y++) {
        for(int x=0; x<width; x++) {
            for(int c=0; c<3; c++) {
                out.at<float>(y, x*3+c) = edge.at<float>(y, x) * out.at<float>(y, x*3+c);
            }
        }
    }

    // 結果の出力
    cv::namedWindow("Input");
    cv::namedWindow("Output");
    cv::namedWindow("Edge");
    cv::imshow("Input", img);
    cv::imshow("Output", out);
    cv::imshow("Edge", edge);
    cv::waitKey(0);
    cv::destroyAllWindows();

    // 結果画像の保存
    out.convertTo(out, CV_8UC3, 255.0);
    cv::imwrite("output.jpg", out);
}