1 /*M/////////////////////////////////////////////////////////////////////////////////////// 2 // 3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. 4 // 5 // By downloading, copying, installing or using the software you agree to this license. 6 // If you do not agree to this license, do not download, install, 7 // copy or use the software. 8 // 9 // 10 // License Agreement 11 // For Open Source Computer Vision Library 12 // 13 // Copyright (C) 2013, OpenCV Foundation, all rights reserved. 14 // Third party copyrights are property of their respective owners. 15 // 16 // Redistribution and use in source and binary forms, with or without modification, 17 // are permitted provided that the following conditions are met: 18 // 19 // * Redistribution's of source code must retain the above copyright notice, 20 // this list of conditions and the following disclaimer. 21 // 22 // * Redistribution's in binary form must reproduce the above copyright notice, 23 // this list of conditions and the following disclaimer in the documentation 24 // and/or other materials provided with the distribution. 25 // 26 // * The name of the copyright holders may not be used to endorse or promote products 27 // derived from this software without specific prior written permission. 28 // 29 // This software is provided by the copyright holders and contributors "as is" and 30 // any express or implied warranties, including, but not limited to, the implied 31 // warranties of merchantability and fitness for a particular purpose are disclaimed. 32 // In no event shall the Intel Corporation or contributors be liable for any direct, 33 // indirect, incidental, special, exemplary, or consequential damages 34 // (including, but not limited to, procurement of substitute goods or services; 35 // loss of use, data, or profits; or business interruption) however caused 36 // and on any theory of liability, whether in contract, strict liability, 37 // or tort (including negligence or otherwise) arising in any way out of 38 // the use of this software, even if advised of the possibility of such damage. 39 // 40 //M*/ 41 42 #include "precomp.hpp" 43 #include "opencv2/photo.hpp" 44 #include <iostream> 45 #include <stdlib.h> 46 #include <limits> 47 #include "math.h" 48 49 50 using namespace std; 51 using namespace cv; 52 53 double myinf = std::numeric_limits<double>::infinity(); 54 55 class Domain_Filter 56 { 57 public: 58 Mat ct_H, ct_V, horiz, vert, O, O_t, lower_idx, upper_idx; 59 void init(const Mat &img, int flags, float sigma_s, float sigma_r); 60 void getGradientx( const Mat &img, Mat &gx); 61 void getGradienty( const Mat &img, Mat &gy); 62 void diffx(const Mat &img, Mat &temp); 63 void diffy(const Mat &img, Mat &temp); 64 void find_magnitude(Mat &img, Mat &mag); 65 void compute_boxfilter(Mat &output, Mat &hz, Mat &psketch, float radius); 66 void compute_Rfilter(Mat &O, Mat &horiz, float sigma_h); 67 void compute_NCfilter(Mat &O, Mat &horiz, Mat &psketch, float radius); 68 void filter(const Mat &img, Mat &res, float sigma_s, float sigma_r, int flags); 69 void pencil_sketch(const Mat &img, Mat &sketch, Mat &color_res, float sigma_s, float sigma_r, float shade_factor); 70 void Depth_of_field(const Mat &img, Mat &img1, float sigma_s, float sigma_r); 71 }; 72 73 void Domain_Filter::diffx(const Mat &img, Mat &temp) 74 { 75 int channel = img.channels(); 76 77 for(int i = 0; i < img.size().height; i++) 78 for(int j = 0; j < img.size().width-1; j++) 79 { 80 for(int c =0; c < channel; c++) 81 { 82 temp.at<float>(i,j*channel+c) = 83 img.at<float>(i,(j+1)*channel+c) - img.at<float>(i,j*channel+c); 84 } 85 } 86 } 87 88 void Domain_Filter::diffy(const Mat &img, Mat &temp) 89 { 90 int channel = img.channels(); 91 92 for(int i = 0; i < img.size().height-1; i++) 93 for(int j = 0; j < img.size().width; j++) 94 { 95 for(int c =0; c < channel; c++) 96 { 97 temp.at<float>(i,j*channel+c) = 98 img.at<float>((i+1),j*channel+c) - img.at<float>(i,j*channel+c); 99 } 100 } 101 } 102 103 void Domain_Filter::getGradientx( const Mat &img, Mat &gx) 104 { 105 int w = img.cols; 106 int h = img.rows; 107 int channel = img.channels(); 108 109 for(int i=0;i<h;i++) 110 for(int j=0;j<w;j++) 111 for(int c=0;c<channel;++c) 112 { 113 gx.at<float>(i,j*channel+c) = 114 img.at<float>(i,(j+1)*channel+c) - img.at<float>(i,j*channel+c); 115 } 116 } 117 118 void Domain_Filter::getGradienty( const Mat &img, Mat &gy) 119 { 120 int w = img.cols; 121 int h = img.rows; 122 int channel = img.channels(); 123 124 for(int i=0;i<h;i++) 125 for(int j=0;j<w;j++) 126 for(int c=0;c<channel;++c) 127 { 128 gy.at<float>(i,j*channel+c) = 129 img.at<float>(i+1,j*channel+c) - img.at<float>(i,j*channel+c); 130 131 } 132 } 133 134 void Domain_Filter::find_magnitude(Mat &img, Mat &mag) 135 { 136 int h = img.rows; 137 int w = img.cols; 138 139 vector <Mat> planes; 140 split(img, planes); 141 142 Mat magXR = Mat(h, w, CV_32FC1); 143 Mat magYR = Mat(h, w, CV_32FC1); 144 145 Mat magXG = Mat(h, w, CV_32FC1); 146 Mat magYG = Mat(h, w, CV_32FC1); 147 148 Mat magXB = Mat(h, w, CV_32FC1); 149 Mat magYB = Mat(h, w, CV_32FC1); 150 151 Sobel(planes[0], magXR, CV_32FC1, 1, 0, 3); 152 Sobel(planes[0], magYR, CV_32FC1, 0, 1, 3); 153 154 Sobel(planes[1], magXG, CV_32FC1, 1, 0, 3); 155 Sobel(planes[1], magYG, CV_32FC1, 0, 1, 3); 156 157 Sobel(planes[2], magXB, CV_32FC1, 1, 0, 3); 158 Sobel(planes[2], magYB, CV_32FC1, 0, 1, 3); 159 160 Mat mag1 = Mat(h,w,CV_32FC1); 161 Mat mag2 = Mat(h,w,CV_32FC1); 162 Mat mag3 = Mat(h,w,CV_32FC1); 163 164 magnitude(magXR,magYR,mag1); 165 magnitude(magXG,magYG,mag2); 166 magnitude(magXB,magYB,mag3); 167 168 mag = mag1 + mag2 + mag3; 169 mag = 1.0f - mag; 170 } 171 172 void Domain_Filter::compute_Rfilter(Mat &output, Mat &hz, float sigma_h) 173 { 174 int h = output.rows; 175 int w = output.cols; 176 int channel = output.channels(); 177 178 float a = (float) exp((-1.0 * sqrt(2.0)) / sigma_h); 179 180 Mat temp = Mat(h,w,CV_32FC3); 181 182 output.copyTo(temp); 183 Mat V = Mat(h,w,CV_32FC1); 184 185 for(int i=0;i<h;i++) 186 for(int j=0;j<w;j++) 187 V.at<float>(i,j) = pow(a,hz.at<float>(i,j)); 188 189 for(int i=0; i<h; i++) 190 { 191 for(int j =1; j < w; j++) 192 { 193 for(int c = 0; c<channel; c++) 194 { 195 temp.at<float>(i,j*channel+c) = temp.at<float>(i,j*channel+c) + 196 (temp.at<float>(i,(j-1)*channel+c) - temp.at<float>(i,j*channel+c)) * V.at<float>(i,j); 197 } 198 } 199 } 200 201 for(int i=0; i<h; i++) 202 { 203 for(int j =w-2; j >= 0; j--) 204 { 205 for(int c = 0; c<channel; c++) 206 { 207 temp.at<float>(i,j*channel+c) = temp.at<float>(i,j*channel+c) + 208 (temp.at<float>(i,(j+1)*channel+c) - temp.at<float>(i,j*channel+c))*V.at<float>(i,j+1); 209 } 210 } 211 } 212 213 temp.copyTo(output); 214 } 215 216 void Domain_Filter::compute_boxfilter(Mat &output, Mat &hz, Mat &psketch, float radius) 217 { 218 int h = output.rows; 219 int w = output.cols; 220 Mat lower_pos = Mat(h,w,CV_32FC1); 221 Mat upper_pos = Mat(h,w,CV_32FC1); 222 223 lower_pos = hz - radius; 224 upper_pos = hz + radius; 225 226 lower_idx = Mat::zeros(h,w,CV_32FC1); 227 upper_idx = Mat::zeros(h,w,CV_32FC1); 228 229 Mat domain_row = Mat::zeros(1,w+1,CV_32FC1); 230 231 for(int i=0;i<h;i++) 232 { 233 for(int j=0;j<w;j++) 234 domain_row.at<float>(0,j) = hz.at<float>(i,j); 235 domain_row.at<float>(0,w) = (float) myinf; 236 237 Mat lower_pos_row = Mat::zeros(1,w,CV_32FC1); 238 Mat upper_pos_row = Mat::zeros(1,w,CV_32FC1); 239 240 for(int j=0;j<w;j++) 241 { 242 lower_pos_row.at<float>(0,j) = lower_pos.at<float>(i,j); 243 upper_pos_row.at<float>(0,j) = upper_pos.at<float>(i,j); 244 } 245 246 Mat temp_lower_idx = Mat::zeros(1,w,CV_32FC1); 247 Mat temp_upper_idx = Mat::zeros(1,w,CV_32FC1); 248 249 for(int j=0;j<w;j++) 250 { 251 if(domain_row.at<float>(0,j) > lower_pos_row.at<float>(0,0)) 252 { 253 temp_lower_idx.at<float>(0,0) = (float) j; 254 break; 255 } 256 } 257 for(int j=0;j<w;j++) 258 { 259 if(domain_row.at<float>(0,j) > upper_pos_row.at<float>(0,0)) 260 { 261 temp_upper_idx.at<float>(0,0) = (float) j; 262 break; 263 } 264 } 265 266 int temp = 0; 267 for(int j=1;j<w;j++) 268 { 269 int count=0; 270 for(int k=(int) temp_lower_idx.at<float>(0,j-1);k<w+1;k++) 271 { 272 if(domain_row.at<float>(0,k) > lower_pos_row.at<float>(0,j)) 273 { 274 temp = count; 275 break; 276 } 277 count++; 278 } 279 280 temp_lower_idx.at<float>(0,j) = temp_lower_idx.at<float>(0,j-1) + temp; 281 282 count = 0; 283 for(int k=(int) temp_upper_idx.at<float>(0,j-1);k<w+1;k++) 284 { 285 286 287 if(domain_row.at<float>(0,k) > upper_pos_row.at<float>(0,j)) 288 { 289 temp = count; 290 break; 291 } 292 count++; 293 } 294 295 temp_upper_idx.at<float>(0,j) = temp_upper_idx.at<float>(0,j-1) + temp; 296 } 297 298 for(int j=0;j<w;j++) 299 { 300 lower_idx.at<float>(i,j) = temp_lower_idx.at<float>(0,j) + 1; 301 upper_idx.at<float>(i,j) = temp_upper_idx.at<float>(0,j) + 1; 302 } 303 304 } 305 psketch = upper_idx - lower_idx; 306 } 307 void Domain_Filter::compute_NCfilter(Mat &output, Mat &hz, Mat &psketch, float radius) 308 { 309 int h = output.rows; 310 int w = output.cols; 311 int channel = output.channels(); 312 313 compute_boxfilter(output,hz,psketch,radius); 314 315 Mat box_filter = Mat::zeros(h,w+1,CV_32FC3); 316 317 for(int i = 0; i < h; i++) 318 { 319 box_filter.at<float>(i,1*channel+0) = output.at<float>(i,0*channel+0); 320 box_filter.at<float>(i,1*channel+1) = output.at<float>(i,0*channel+1); 321 box_filter.at<float>(i,1*channel+2) = output.at<float>(i,0*channel+2); 322 for(int j = 2; j < w+1; j++) 323 { 324 for(int c=0;c<channel;c++) 325 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); 326 } 327 } 328 329 Mat indices = Mat::zeros(h,w,CV_32FC1); 330 Mat final = Mat::zeros(h,w,CV_32FC3); 331 332 for(int i=0;i<h;i++) 333 for(int j=0;j<w;j++) 334 indices.at<float>(i,j) = (float) i+1; 335 336 Mat a = Mat::zeros(h,w,CV_32FC1); 337 Mat b = Mat::zeros(h,w,CV_32FC1); 338 339 // Compute the box filter using a summed area table. 340 for(int c=0;c<channel;c++) 341 { 342 Mat flag = Mat::ones(h,w,CV_32FC1); 343 multiply(flag,c+1,flag); 344 345 Mat temp1, temp2; 346 multiply(flag - 1,h*(w+1),temp1); 347 multiply(lower_idx - 1,h,temp2); 348 a = temp1 + temp2 + indices; 349 350 multiply(flag - 1,h*(w+1),temp1); 351 multiply(upper_idx - 1,h,temp2); 352 b = temp1 + temp2 + indices; 353 354 int p,q,r,rem; 355 int p1,q1,r1,rem1; 356 357 // Calculating indices 358 for(int i=0;i<h;i++) 359 { 360 for(int j=0;j<w;j++) 361 { 362 363 r = (int) b.at<float>(i,j)/(h*(w+1)); 364 rem = (int) b.at<float>(i,j) - r*h*(w+1); 365 q = rem/h; 366 p = rem - q*h; 367 if(q==0) 368 { 369 p=h; 370 q=w; 371 r=r-1; 372 } 373 if(p==0) 374 { 375 p=h; 376 q=q-1; 377 } 378 379 r1 = (int) a.at<float>(i,j)/(h*(w+1)); 380 rem1 = (int) a.at<float>(i,j) - r1*h*(w+1); 381 q1 = rem1/h; 382 p1 = rem1 - q1*h; 383 if(p1==0) 384 { 385 p1=h; 386 q1=q1-1; 387 } 388 389 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))) 390 /(upper_idx.at<float>(i,j) - lower_idx.at<float>(i,j)); 391 } 392 } 393 } 394 395 final.copyTo(output); 396 } 397 void Domain_Filter::init(const Mat &img, int flags, float sigma_s, float sigma_r) 398 { 399 int h = img.size().height; 400 int w = img.size().width; 401 int channel = img.channels(); 402 403 //////////////////////////////////// horizontal and vertical partial derivatives ///////////////////////////////// 404 405 Mat derivx = Mat::zeros(h,w-1,CV_32FC3); 406 Mat derivy = Mat::zeros(h-1,w,CV_32FC3); 407 408 diffx(img,derivx); 409 diffy(img,derivy); 410 411 Mat distx = Mat::zeros(h,w,CV_32FC1); 412 Mat disty = Mat::zeros(h,w,CV_32FC1); 413 414 //////////////////////// Compute the l1-norm distance of neighbor pixels //////////////////////////////////////////////// 415 416 for(int i = 0; i < h; i++) 417 for(int j = 0,k=1; j < w-1; j++,k++) 418 for(int c = 0; c < channel; c++) 419 { 420 distx.at<float>(i,k) = 421 distx.at<float>(i,k) + abs(derivx.at<float>(i,j*channel+c)); 422 } 423 424 for(int i = 0,k=1; i < h-1; i++,k++) 425 for(int j = 0; j < w; j++) 426 for(int c = 0; c < channel; c++) 427 { 428 disty.at<float>(k,j) = 429 disty.at<float>(k,j) + abs(derivy.at<float>(i,j*channel+c)); 430 } 431 432 ////////////////////// Compute the derivatives of the horizontal and vertical domain transforms. ///////////////////////////// 433 434 horiz = Mat(h,w,CV_32FC1); 435 vert = Mat(h,w,CV_32FC1); 436 437 Mat final = Mat(h,w,CV_32FC3); 438 439 Mat tempx,tempy; 440 multiply(distx,sigma_s/sigma_r,tempx); 441 multiply(disty,sigma_s/sigma_r,tempy); 442 443 horiz = 1.0f + tempx; 444 vert = 1.0f + tempy; 445 446 O = Mat(h,w,CV_32FC3); 447 img.copyTo(O); 448 449 O_t = Mat(w,h,CV_32FC3); 450 451 if(flags == 2) 452 { 453 454 ct_H = Mat(h,w,CV_32FC1); 455 ct_V = Mat(h,w,CV_32FC1); 456 457 for(int i = 0; i < h; i++) 458 { 459 ct_H.at<float>(i,0) = horiz.at<float>(i,0); 460 for(int j = 1; j < w; j++) 461 { 462 ct_H.at<float>(i,j) = horiz.at<float>(i,j) + ct_H.at<float>(i,j-1); 463 } 464 } 465 466 for(int j = 0; j < w; j++) 467 { 468 ct_V.at<float>(0,j) = vert.at<float>(0,j); 469 for(int i = 1; i < h; i++) 470 { 471 ct_V.at<float>(i,j) = vert.at<float>(i,j) + ct_V.at<float>(i-1,j); 472 } 473 } 474 } 475 476 } 477 478 void Domain_Filter::filter(const Mat &img, Mat &res, float sigma_s = 60, float sigma_r = 0.4, int flags = 1) 479 { 480 int no_of_iter = 3; 481 int h = img.size().height; 482 int w = img.size().width; 483 float sigma_h = sigma_s; 484 485 init(img,flags,sigma_s,sigma_r); 486 487 if(flags == 1) 488 { 489 Mat vert_t = vert.t(); 490 491 for(int i=0;i<no_of_iter;i++) 492 { 493 sigma_h = (float) (sigma_s * sqrt(3.0) * pow(2.0,(no_of_iter - (i+1))) / sqrt(pow(4.0,no_of_iter) -1)); 494 495 compute_Rfilter(O, horiz, sigma_h); 496 497 O_t = O.t(); 498 499 compute_Rfilter(O_t, vert_t, sigma_h); 500 501 O = O_t.t(); 502 503 } 504 } 505 else if(flags == 2) 506 { 507 508 Mat vert_t = ct_V.t(); 509 Mat temp = Mat(h,w,CV_32FC1); 510 Mat temp1 = Mat(w,h,CV_32FC1); 511 512 float radius; 513 514 for(int i=0;i<no_of_iter;i++) 515 { 516 sigma_h = (float) (sigma_s * sqrt(3.0) * pow(2.0,(no_of_iter - (i+1))) / sqrt(pow(4.0,no_of_iter) -1)); 517 518 radius = (float) sqrt(3.0) * sigma_h; 519 520 compute_NCfilter(O, ct_H, temp,radius); 521 522 O_t = O.t(); 523 524 compute_NCfilter(O_t, vert_t, temp1, radius); 525 526 O = O_t.t(); 527 } 528 } 529 530 res = O.clone(); 531 } 532 533 void Domain_Filter::pencil_sketch(const Mat &img, Mat &sketch, Mat &color_res, float sigma_s, float sigma_r, float shade_factor) 534 { 535 536 int no_of_iter = 3; 537 init(img,2,sigma_s,sigma_r); 538 int h = img.size().height; 539 int w = img.size().width; 540 541 /////////////////////// convert to YCBCR model for color pencil drawing ////////////////////////////////////////////////////// 542 543 Mat color_sketch = Mat(h,w,CV_32FC3); 544 545 cvtColor(img,color_sketch,COLOR_BGR2YCrCb); 546 547 vector <Mat> YUV_channel; 548 Mat vert_t = ct_V.t(); 549 550 float sigma_h = sigma_s; 551 552 Mat penx = Mat(h,w,CV_32FC1); 553 554 Mat pen_res = Mat::zeros(h,w,CV_32FC1); 555 Mat peny = Mat(w,h,CV_32FC1); 556 557 Mat peny_t; 558 559 float radius; 560 561 for(int i=0;i<no_of_iter;i++) 562 { 563 sigma_h = (float) (sigma_s * sqrt(3.0) * pow(2.0,(no_of_iter - (i+1))) / sqrt(pow(4.0,no_of_iter) -1)); 564 565 radius = (float) sqrt(3.0) * sigma_h; 566 567 compute_boxfilter(O, ct_H, penx, radius); 568 569 O_t = O.t(); 570 571 compute_boxfilter(O_t, vert_t, peny, radius); 572 573 O = O_t.t(); 574 575 peny_t = peny.t(); 576 577 for(int k=0;k<h;k++) 578 for(int j=0;j<w;j++) 579 pen_res.at<float>(k,j) = (shade_factor * (penx.at<float>(k,j) + peny_t.at<float>(k,j))); 580 581 if(i==0) 582 { 583 sketch = pen_res.clone(); 584 split(color_sketch,YUV_channel); 585 pen_res.copyTo(YUV_channel[0]); 586 merge(YUV_channel,color_sketch); 587 cvtColor(color_sketch,color_res,COLOR_YCrCb2BGR); 588 } 589 590 } 591 } 592
