photo/src/npr.hpp

/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                           License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of the copyright holders may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/

#include "precomp.hpp"
#include "opencv2/photo.hpp"
#include <iostream>
#include <stdlib.h>
#include <limits>
#include "math.h"


using namespace std;
using namespace cv;

double myinf = std::numeric_limits<double>::infinity();

class Domain_Filter
{
    public:
        Mat ct_H, ct_V, horiz, vert, O, O_t, lower_idx, upper_idx;
        void init(const Mat &img, int flags, float sigma_s, float sigma_r);
        void getGradientx( const Mat &img, Mat &gx);
        void getGradienty( const Mat &img, Mat &gy);
        void diffx(const Mat &img, Mat &temp);
        void diffy(const Mat &img, Mat &temp);
        void find_magnitude(Mat &img, Mat &mag);
        void compute_boxfilter(Mat &output, Mat &hz, Mat &psketch, float radius);
        void compute_Rfilter(Mat &O, Mat &horiz, float sigma_h);
        void compute_NCfilter(Mat &O, Mat &horiz, Mat &psketch, float radius);
        void filter(const Mat &img, Mat &res, float sigma_s, float sigma_r, int flags);
        void pencil_sketch(const Mat &img, Mat &sketch, Mat &color_res, float sigma_s, float sigma_r, float shade_factor);
        void Depth_of_field(const Mat &img, Mat &img1, float sigma_s, float sigma_r);
};

void Domain_Filter::diffx(const Mat &img, Mat &temp)
{
    int channel = img.channels();

    for(int i = 0; i < img.size().height; i++)
        for(int j = 0; j < img.size().width-1; j++)
        {
            for(int c =0; c < channel; c++)
            {
                temp.at<float>(i,j*channel+c) =
                    img.at<float>(i,(j+1)*channel+c) - img.at<float>(i,j*channel+c);
            }
        }
}

void Domain_Filter::diffy(const Mat &img, Mat &temp)
{
    int channel = img.channels();

    for(int i = 0; i < img.size().height-1; i++)
        for(int j = 0; j < img.size().width; j++)
        {
            for(int c =0; c < channel; c++)
            {
                temp.at<float>(i,j*channel+c) =
                    img.at<float>((i+1),j*channel+c) - img.at<float>(i,j*channel+c);
            }
        }
}

void Domain_Filter::getGradientx( const Mat &img, Mat &gx)
{
    int w = img.cols;
    int h = img.rows;
    int channel = img.channels();

    for(int i=0;i<h;i++)
        for(int j=0;j<w;j++)
            for(int c=0;c<channel;++c)
            {
                gx.at<float>(i,j*channel+c) =
                    img.at<float>(i,(j+1)*channel+c) - img.at<float>(i,j*channel+c);
            }
}

void Domain_Filter::getGradienty( const Mat &img, Mat &gy)
{
    int w = img.cols;
    int h = img.rows;
    int channel = img.channels();

    for(int i=0;i<h;i++)
        for(int j=0;j<w;j++)
            for(int c=0;c<channel;++c)
            {
                gy.at<float>(i,j*channel+c) =
                    img.at<float>(i+1,j*channel+c) - img.at<float>(i,j*channel+c);

            }
}

void Domain_Filter::find_magnitude(Mat &img, Mat &mag)
{
    int h = img.rows;
    int w = img.cols;

    vector <Mat> planes;
    split(img, planes);

    Mat magXR = Mat(h, w, CV_32FC1);
    Mat magYR = Mat(h, w, CV_32FC1);

    Mat magXG = Mat(h, w, CV_32FC1);
    Mat magYG = Mat(h, w, CV_32FC1);

    Mat magXB = Mat(h, w, CV_32FC1);
    Mat magYB = Mat(h, w, CV_32FC1);

    Sobel(planes[0], magXR, CV_32FC1, 1, 0, 3);
    Sobel(planes[0], magYR, CV_32FC1, 0, 1, 3);

    Sobel(planes[1], magXG, CV_32FC1, 1, 0, 3);
    Sobel(planes[1], magYG, CV_32FC1, 0, 1, 3);

    Sobel(planes[2], magXB, CV_32FC1, 1, 0, 3);
    Sobel(planes[2], magYB, CV_32FC1, 0, 1, 3);

    Mat mag1 = Mat(h,w,CV_32FC1);
    Mat mag2 = Mat(h,w,CV_32FC1);
    Mat mag3 = Mat(h,w,CV_32FC1);

    magnitude(magXR,magYR,mag1);
    magnitude(magXG,magYG,mag2);
    magnitude(magXB,magYB,mag3);

    mag = mag1 + mag2 + mag3;
    mag = 1.0f - mag;
}

void Domain_Filter::compute_Rfilter(Mat &output, Mat &hz, float sigma_h)
{
    int h = output.rows;
    int w = output.cols;
    int channel = output.channels();

    float a = (float) exp((-1.0 * sqrt(2.0)) / sigma_h);

    Mat temp = Mat(h,w,CV_32FC3);

    output.copyTo(temp);
    Mat V = Mat(h,w,CV_32FC1);

    for(int i=0;i<h;i++)
        for(int j=0;j<w;j++)
            V.at<float>(i,j) = pow(a,hz.at<float>(i,j));

   for(int i=0; i<h; i++)
    {
        for(int j =1; j < w; j++)
        {
            for(int c = 0; c<channel; c++)
            {
                temp.at<float>(i,j*channel+c) = temp.at<float>(i,j*channel+c) +
                    (temp.at<float>(i,(j-1)*channel+c) - temp.at<float>(i,j*channel+c)) * V.at<float>(i,j);
            }
        }
    }

    for(int i=0; i<h; i++)
    {
        for(int j =w-2; j >= 0; j--)
        {
            for(int c = 0; c<channel; c++)
            {
                temp.at<float>(i,j*channel+c) = temp.at<float>(i,j*channel+c) +
                    (temp.at<float>(i,(j+1)*channel+c) - temp.at<float>(i,j*channel+c))*V.at<float>(i,j+1);
            }
        }
    }

    temp.copyTo(output);
}

void Domain_Filter::compute_boxfilter(Mat &output, Mat &hz, Mat &psketch, float radius)
{
    int h = output.rows;
    int w = output.cols;
    Mat lower_pos = Mat(h,w,CV_32FC1);
    Mat upper_pos = Mat(h,w,CV_32FC1);

    lower_pos = hz - radius;
    upper_pos = hz + radius;

    lower_idx = Mat::zeros(h,w,CV_32FC1);
    upper_idx = Mat::zeros(h,w,CV_32FC1);

    Mat domain_row = Mat::zeros(1,w+1,CV_32FC1);

    for(int i=0;i<h;i++)
    {
        for(int j=0;j<w;j++)
            domain_row.at<float>(0,j) = hz.at<float>(i,j);
        domain_row.at<float>(0,w) = (float) myinf;

        Mat lower_pos_row = Mat::zeros(1,w,CV_32FC1);
        Mat upper_pos_row = Mat::zeros(1,w,CV_32FC1);

        for(int j=0;j<w;j++)
        {
            lower_pos_row.at<float>(0,j) = lower_pos.at<float>(i,j);
            upper_pos_row.at<float>(0,j) = upper_pos.at<float>(i,j);
        }

        Mat temp_lower_idx = Mat::zeros(1,w,CV_32FC1);
        Mat temp_upper_idx = Mat::zeros(1,w,CV_32FC1);

        for(int j=0;j<w;j++)
        {
            if(domain_row.at<float>(0,j) > lower_pos_row.at<float>(0,0))
            {
                temp_lower_idx.at<float>(0,0) = (float) j;
                break;
            }
        }
        for(int j=0;j<w;j++)
        {
            if(domain_row.at<float>(0,j) > upper_pos_row.at<float>(0,0))
            {
                temp_upper_idx.at<float>(0,0) = (float) j;
                break;
            }
        }

        int temp = 0;
        for(int j=1;j<w;j++)
        {
            int count=0;
            for(int k=(int) temp_lower_idx.at<float>(0,j-1);k<w+1;k++)
            {
                if(domain_row.at<float>(0,k) > lower_pos_row.at<float>(0,j))
                {
                    temp = count;
                    break;
                }
                count++;
            }

            temp_lower_idx.at<float>(0,j) = temp_lower_idx.at<float>(0,j-1) + temp;

            count = 0;
            for(int k=(int) temp_upper_idx.at<float>(0,j-1);k<w+1;k++)
            {


                if(domain_row.at<float>(0,k) > upper_pos_row.at<float>(0,j))
                {
                    temp = count;
                    break;
                }
                count++;
            }

            temp_upper_idx.at<float>(0,j) = temp_upper_idx.at<float>(0,j-1) + temp;
        }

        for(int j=0;j<w;j++)
        {
            lower_idx.at<float>(i,j) = temp_lower_idx.at<float>(0,j) + 1;
            upper_idx.at<float>(i,j) = temp_upper_idx.at<float>(0,j) + 1;
        }

    }
    psketch = upper_idx - lower_idx;
}
void Domain_Filter::compute_NCfilter(Mat &output, Mat &hz, Mat &psketch, float radius)
{
    int h = output.rows;
    int w = output.cols;
    int channel = output.channels();

    compute_boxfilter(output,hz,psketch,radius);

    Mat box_filter = Mat::zeros(h,w+1,CV_32FC3);

    for(int i = 0; i < h; i++)
    {
        box_filter.at<float>(i,1*channel+0) = output.at<float>(i,0*channel+0);
        box_filter.at<float>(i,1*channel+1) = output.at<float>(i,0*channel+1);
        box_filter.at<float>(i,1*channel+2) = output.at<float>(i,0*channel+2);
        for(int j = 2; j < w+1; j++)
        {
            for(int c=0;c<channel;c++)
                box_filter.at<float>(i,j*channel+c) = output.at<float>(i,(j-1)*channel+c) + box_filter.at<float>(i,(j-1)*channel+c);
        }
    }

    Mat indices = Mat::zeros(h,w,CV_32FC1);
    Mat final =   Mat::zeros(h,w,CV_32FC3);

    for(int i=0;i<h;i++)
        for(int j=0;j<w;j++)
            indices.at<float>(i,j) = (float) i+1;

    Mat a = Mat::zeros(h,w,CV_32FC1);
    Mat b = Mat::zeros(h,w,CV_32FC1);

    // Compute the box filter using a summed area table.
    for(int c=0;c<channel;c++)
    {
        Mat flag = Mat::ones(h,w,CV_32FC1);
        multiply(flag,c+1,flag);

        Mat temp1, temp2;
        multiply(flag - 1,h*(w+1),temp1);
        multiply(lower_idx - 1,h,temp2);
        a = temp1 + temp2 + indices;

        multiply(flag - 1,h*(w+1),temp1);
        multiply(upper_idx - 1,h,temp2);
        b = temp1 + temp2 + indices;

        int p,q,r,rem;
        int p1,q1,r1,rem1;

        // Calculating indices
        for(int i=0;i<h;i++)
        {
            for(int j=0;j<w;j++)
            {

                r = (int) b.at<float>(i,j)/(h*(w+1));
                rem = (int) b.at<float>(i,j) - r*h*(w+1);
                q = rem/h;
                p = rem - q*h;
                if(q==0)
                {
                    p=h;
                    q=w;
                    r=r-1;
                }
                if(p==0)
                {
                    p=h;
                    q=q-1;
                }

                r1 = (int) a.at<float>(i,j)/(h*(w+1));
                rem1 = (int) a.at<float>(i,j) - r1*h*(w+1);
                q1 = rem1/h;
                p1 = rem1 - q1*h;
                if(p1==0)
                {
                    p1=h;
                    q1=q1-1;
                }

                final.at<float>(i,j*channel+2-c) = (box_filter.at<float>(p-1,q*channel+(2-r)) - box_filter.at<float>(p1-1,q1*channel+(2-r1)))
                    /(upper_idx.at<float>(i,j) - lower_idx.at<float>(i,j));
            }
        }
    }

    final.copyTo(output);
}
void Domain_Filter::init(const Mat &img, int flags, float sigma_s, float sigma_r)
{
    int h = img.size().height;
    int w = img.size().width;
    int channel = img.channels();

    ////////////////////////////////////     horizontal and vertical partial derivatives /////////////////////////////////

    Mat derivx = Mat::zeros(h,w-1,CV_32FC3);
    Mat derivy = Mat::zeros(h-1,w,CV_32FC3);

    diffx(img,derivx);
    diffy(img,derivy);

    Mat distx = Mat::zeros(h,w,CV_32FC1);
    Mat disty = Mat::zeros(h,w,CV_32FC1);

    //////////////////////// Compute the l1-norm distance of neighbor pixels ////////////////////////////////////////////////

    for(int i = 0; i < h; i++)
        for(int j = 0,k=1; j < w-1; j++,k++)
            for(int c = 0; c < channel; c++)
            {
                distx.at<float>(i,k) =
                    distx.at<float>(i,k) + abs(derivx.at<float>(i,j*channel+c));
            }

    for(int i = 0,k=1; i < h-1; i++,k++)
        for(int j = 0; j < w; j++)
            for(int c = 0; c < channel; c++)
            {
                disty.at<float>(k,j) =
                    disty.at<float>(k,j) + abs(derivy.at<float>(i,j*channel+c));
            }

    ////////////////////// Compute the derivatives of the horizontal and vertical domain transforms. /////////////////////////////

    horiz = Mat(h,w,CV_32FC1);
    vert = Mat(h,w,CV_32FC1);

    Mat final = Mat(h,w,CV_32FC3);

    Mat tempx,tempy;
    multiply(distx,sigma_s/sigma_r,tempx);
    multiply(disty,sigma_s/sigma_r,tempy);

    horiz = 1.0f + tempx;
    vert = 1.0f + tempy;

    O = Mat(h,w,CV_32FC3);
    img.copyTo(O);

    O_t = Mat(w,h,CV_32FC3);

    if(flags == 2)
    {

        ct_H = Mat(h,w,CV_32FC1);
        ct_V = Mat(h,w,CV_32FC1);

        for(int i = 0; i < h; i++)
        {
            ct_H.at<float>(i,0) = horiz.at<float>(i,0);
            for(int j = 1; j < w; j++)
            {
                ct_H.at<float>(i,j) = horiz.at<float>(i,j) + ct_H.at<float>(i,j-1);
            }
        }

        for(int j = 0; j < w; j++)
        {
            ct_V.at<float>(0,j) = vert.at<float>(0,j);
            for(int i = 1; i < h; i++)
            {
                ct_V.at<float>(i,j) = vert.at<float>(i,j) + ct_V.at<float>(i-1,j);
            }
        }
    }

}

void Domain_Filter::filter(const Mat &img, Mat &res, float sigma_s = 60, float sigma_r = 0.4, int flags = 1)
{
    int no_of_iter = 3;
    int h = img.size().height;
    int w = img.size().width;
    float sigma_h = sigma_s;

    init(img,flags,sigma_s,sigma_r);

    if(flags == 1)
    {
        Mat vert_t = vert.t();

        for(int i=0;i<no_of_iter;i++)
        {
            sigma_h = (float) (sigma_s * sqrt(3.0) * pow(2.0,(no_of_iter - (i+1))) / sqrt(pow(4.0,no_of_iter) -1));

            compute_Rfilter(O, horiz, sigma_h);

            O_t = O.t();

            compute_Rfilter(O_t, vert_t, sigma_h);

            O = O_t.t();

        }
    }
    else if(flags == 2)
    {

        Mat vert_t = ct_V.t();
        Mat temp = Mat(h,w,CV_32FC1);
        Mat temp1 = Mat(w,h,CV_32FC1);

        float radius;

        for(int i=0;i<no_of_iter;i++)
        {
            sigma_h = (float) (sigma_s * sqrt(3.0) * pow(2.0,(no_of_iter - (i+1))) / sqrt(pow(4.0,no_of_iter) -1));

            radius = (float) sqrt(3.0) * sigma_h;

            compute_NCfilter(O, ct_H, temp,radius);

            O_t = O.t();

            compute_NCfilter(O_t, vert_t, temp1, radius);

            O = O_t.t();
        }
    }

    res = O.clone();
}

void Domain_Filter::pencil_sketch(const Mat &img, Mat &sketch, Mat &color_res, float sigma_s, float sigma_r, float shade_factor)
{

    int no_of_iter = 3;
    init(img,2,sigma_s,sigma_r);
    int h = img.size().height;
    int w = img.size().width;

    /////////////////////// convert to YCBCR model for color pencil drawing //////////////////////////////////////////////////////

    Mat color_sketch = Mat(h,w,CV_32FC3);

    cvtColor(img,color_sketch,COLOR_BGR2YCrCb);

    vector <Mat> YUV_channel;
    Mat vert_t = ct_V.t();

    float sigma_h = sigma_s;

    Mat penx = Mat(h,w,CV_32FC1);

    Mat pen_res = Mat::zeros(h,w,CV_32FC1);
    Mat peny = Mat(w,h,CV_32FC1);

    Mat peny_t;

    float radius;

    for(int i=0;i<no_of_iter;i++)
    {
        sigma_h = (float) (sigma_s * sqrt(3.0) * pow(2.0,(no_of_iter - (i+1))) / sqrt(pow(4.0,no_of_iter) -1));

        radius = (float) sqrt(3.0) * sigma_h;

        compute_boxfilter(O, ct_H, penx, radius);

        O_t = O.t();

        compute_boxfilter(O_t, vert_t, peny, radius);

        O = O_t.t();

        peny_t = peny.t();

        for(int k=0;k<h;k++)
            for(int j=0;j<w;j++)
                pen_res.at<float>(k,j) = (shade_factor * (penx.at<float>(k,j) + peny_t.at<float>(k,j)));

        if(i==0)
        {
            sketch = pen_res.clone();
            split(color_sketch,YUV_channel);
            pen_res.copyTo(YUV_channel[0]);
            merge(YUV_channel,color_sketch);
            cvtColor(color_sketch,color_res,COLOR_YCrCb2BGR);
        }

    }
}