1 /*
2 * By downloading, copying, installing or using the software you agree to this license.
3 * If you do not agree to this license, do not download, install,
4 * copy or use the software.
5 *
6 *
7 * License Agreement
8 * For Open Source Computer Vision Library
9 * (3 - clause BSD License)
10 *
11 * Redistribution and use in source and binary forms, with or without modification,
12 * are permitted provided that the following conditions are met :
13 *
14 * * Redistributions of source code must retain the above copyright notice,
15 * this list of conditions and the following disclaimer.
16 *
17 * * Redistributions in binary form must reproduce the above copyright notice,
18 * this list of conditions and the following disclaimer in the documentation
19 * and / or other materials provided with the distribution.
20 *
21 * * Neither the names of the copyright holders nor the names of the contributors
22 * may be used to endorse or promote products derived from this software
23 * without specific prior written permission.
24 *
25 * This software is provided by the copyright holders and contributors "as is" and
26 * any express or implied warranties, including, but not limited to, the implied
27 * warranties of merchantability and fitness for a particular purpose are disclaimed.
28 * In no event shall copyright holders or contributors be liable for any direct,
29 * indirect, incidental, special, exemplary, or consequential damages
30 * (including, but not limited to, procurement of substitute goods or services;
31 * loss of use, data, or profits; or business interruption) however caused
32 * and on any theory of liability, whether in contract, strict liability,
33 * or tort(including negligence or otherwise) arising in any way out of
34 * the use of this software, even if advised of the possibility of such damage.
35 */
36
37 #include "precomp.hpp"
38
39 #include <cmath>
40 #include <vector>
41 #include <memory>
42 #include <stdlib.h>
43 #include <iostream>
44 #include <iterator>
45 #include <algorithm>
46
47
48 #ifdef HAVE_EIGEN
49 # if defined __GNUC__ && defined __APPLE__
50 # pragma GCC diagnostic ignored "-Wshadow"
51 # endif
52 # include <Eigen/Dense>
53 # include <Eigen/SparseCore>
54 # include <Eigen/SparseCholesky>
55 # include <Eigen/IterativeLinearSolvers>
56 # include <Eigen/Sparse>
57
58
59
60
61 #if __cplusplus <= 199711L
62 #include <map>
63 typedef std::map<long long /* hash */, int /* vert id */> mapId;
64 #else
65 #include <unordered_map>
66 typedef std::unordered_map<long long /* hash */, int /* vert id */> mapId;
67 #endif
68
69 namespace cv
70 {
71 namespace ximgproc
72 {
73
74 class FastBilateralSolverFilterImpl : public FastBilateralSolverFilter
75 {
76 public:
77
create(InputArray guide,double sigma_spatial,double sigma_luma,double sigma_chroma,double lambda,int num_iter,double max_tol)78 static Ptr<FastBilateralSolverFilterImpl> create(InputArray guide, double sigma_spatial, double sigma_luma, double sigma_chroma, double lambda, int num_iter, double max_tol)
79 {
80 CV_Assert(guide.type() == CV_8UC1 || guide.type() == CV_8UC3);
81 FastBilateralSolverFilterImpl *fbs = new FastBilateralSolverFilterImpl();
82 Mat gui = guide.getMat();
83 fbs->init(gui,sigma_spatial,sigma_luma,sigma_chroma,lambda,num_iter,max_tol);
84 return Ptr<FastBilateralSolverFilterImpl>(fbs);
85 }
86
filter(InputArray src,InputArray confidence,OutputArray dst)87 void filter(InputArray src, InputArray confidence, OutputArray dst) CV_OVERRIDE
88 {
89
90 CV_Assert(!src.empty() && (src.depth() == CV_8U || src.depth() == CV_16S || src.depth() == CV_16U || src.depth() == CV_32F) && src.channels()<=4);
91 CV_Assert(!confidence.empty() && (confidence.depth() == CV_8U || confidence.depth() == CV_32F) && confidence.channels()==1);
92 if (src.rows() != rows || src.cols() != cols)
93 {
94 CV_Error(Error::StsBadSize, "Size of the filtered image must be equal to the size of the guide image");
95 return;
96 }
97 if (confidence.rows() != rows || confidence.cols() != cols)
98 {
99 CV_Error(Error::StsBadSize, "Size of the confidence image must be equal to the size of the guide image");
100 return;
101 }
102
103 std::vector<Mat> src_channels;
104 std::vector<Mat> dst_channels;
105 if(src.channels()==1)
106 src_channels.push_back(src.getMat());
107 else
108 split(src,src_channels);
109
110 Mat conf = confidence.getMat();
111
112 for(int i=0;i<src.channels();i++)
113 {
114 Mat cur_res = src_channels[i].clone();
115
116 solve(cur_res,conf,cur_res);
117 cur_res.convertTo(cur_res, src.type());
118 dst_channels.push_back(cur_res);
119 }
120
121 dst.create(src.size(),src_channels[0].type());
122 if(src.channels()==1)
123 {
124 Mat& dstMat = dst.getMatRef();
125 dstMat = dst_channels[0];
126 }
127 else
128 merge(dst_channels,dst);
129 CV_Assert(src.type() == dst.type() && src.size() == dst.size());
130 }
131
132 // protected:
133 void solve(cv::Mat& src, cv::Mat& confidence, cv::Mat& dst);
134 void init(cv::Mat& reference, double sigma_spatial, double sigma_luma, double sigma_chroma, double lambda, int num_iter, double max_tol);
135
136 void Splat(Eigen::VectorXf& input, Eigen::VectorXf& dst);
137 void Blur(Eigen::VectorXf& input, Eigen::VectorXf& dst);
138 void Slice(Eigen::VectorXf& input, Eigen::VectorXf& dst);
139
diagonal(Eigen::VectorXf & v,Eigen::SparseMatrix<float> & mat)140 void diagonal(Eigen::VectorXf& v,Eigen::SparseMatrix<float>& mat)
141 {
142 mat = Eigen::SparseMatrix<float>(v.size(),v.size());
143 for (int i = 0; i < int(v.size()); i++)
144 {
145 mat.insert(i,i) = v(i);
146 }
147 }
148
149
150
151 private:
152
153 int npixels;
154 int nvertices;
155 int dim;
156 int cols;
157 int rows;
158 std::vector<int> splat_idx;
159 std::vector<std::pair<int, int> > blur_idx;
160 Eigen::VectorXf m;
161 Eigen::VectorXf n;
162 Eigen::SparseMatrix<float, Eigen::ColMajor> blurs;
163 Eigen::SparseMatrix<float, Eigen::ColMajor> S;
164 Eigen::SparseMatrix<float, Eigen::ColMajor> Dn;
165 Eigen::SparseMatrix<float, Eigen::ColMajor> Dm;
166
167 struct grid_params
168 {
169 float spatialSigma;
170 float lumaSigma;
171 float chromaSigma;
grid_paramscv::ximgproc::FastBilateralSolverFilterImpl::grid_params172 grid_params()
173 {
174 spatialSigma = 8.0;
175 lumaSigma = 8.0;
176 chromaSigma = 8.0;
177 }
178 };
179
180 struct bs_params
181 {
182 float lam;
183 float A_diag_min;
184 float cg_tol;
185 int cg_maxiter;
bs_paramscv::ximgproc::FastBilateralSolverFilterImpl::bs_params186 bs_params()
187 {
188 lam = 128.0;
189 A_diag_min = 1e-5f;
190 cg_tol = 1e-5f;
191 cg_maxiter = 25;
192 }
193 };
194
195 grid_params grid_param;
196 bs_params bs_param;
197
198 };
199
200
201
init(cv::Mat & reference,double sigma_spatial,double sigma_luma,double sigma_chroma,double lambda,int num_iter,double max_tol)202 void FastBilateralSolverFilterImpl::init(cv::Mat& reference, double sigma_spatial, double sigma_luma, double sigma_chroma, double lambda, int num_iter, double max_tol)
203 {
204
205 bs_param.lam = lambda;
206 bs_param.cg_maxiter = num_iter;
207 bs_param.cg_tol = max_tol;
208
209 if(reference.channels()==1)
210 {
211 dim = 3;
212 cols = reference.cols;
213 rows = reference.rows;
214 npixels = cols*rows;
215 long long hash_vec[3];
216 for (int i = 0; i < 3; ++i)
217 hash_vec[i] = static_cast<long long>(std::pow(255, i));
218
219 mapId hashed_coords;
220 #if __cplusplus <= 199711L
221 #else
222 hashed_coords.reserve(cols*rows);
223 #endif
224
225 const unsigned char* pref = (const unsigned char*)reference.data;
226 int vert_idx = 0;
227 int pix_idx = 0;
228
229 // construct Splat(Slice) matrices
230 splat_idx.resize(npixels);
231 for (int y = 0; y < rows; ++y)
232 {
233 for (int x = 0; x < cols; ++x)
234 {
235 long long coord[3];
236 coord[0] = int(x / sigma_spatial);
237 coord[1] = int(y / sigma_spatial);
238 coord[2] = int(pref[0] / sigma_luma);
239
240 // convert the coordinate to a hash value
241 long long hash_coord = 0;
242 for (int i = 0; i < 3; ++i)
243 hash_coord += coord[i] * hash_vec[i];
244
245 // pixels whom are alike will have the same hash value.
246 // We only want to keep a unique list of hash values, therefore make sure we only insert
247 // unique hash values.
248 mapId::iterator it = hashed_coords.find(hash_coord);
249 if (it == hashed_coords.end())
250 {
251 hashed_coords.insert(std::pair<long long, int>(hash_coord, vert_idx));
252 splat_idx[pix_idx] = vert_idx;
253 ++vert_idx;
254 }
255 else
256 {
257 splat_idx[pix_idx] = it->second;
258 }
259
260 pref += 1; // skip 1 bytes (y)
261 ++pix_idx;
262 }
263 }
264 nvertices = static_cast<int>(hashed_coords.size());
265
266 // construct Blur matrices
267 Eigen::VectorXf ones_nvertices = Eigen::VectorXf::Ones(nvertices);
268 diagonal(ones_nvertices,blurs);
269 blurs *= 10;
270 for(int offset = -1; offset <= 1;++offset)
271 {
272 if(offset == 0) continue;
273 for (int i = 0; i < dim; ++i)
274 {
275 Eigen::SparseMatrix<float, Eigen::ColMajor> blur_temp(hashed_coords.size(), hashed_coords.size());
276 blur_temp.reserve(Eigen::VectorXi::Constant(nvertices,6));
277 long long offset_hash_coord = offset * hash_vec[i];
278 for (mapId::iterator it = hashed_coords.begin(); it != hashed_coords.end(); ++it)
279 {
280 long long neighb_coord = it->first + offset_hash_coord;
281 mapId::iterator it_neighb = hashed_coords.find(neighb_coord);
282 if (it_neighb != hashed_coords.end())
283 {
284 blur_temp.insert(it->second,it_neighb->second) = 1.0f;
285 blur_idx.push_back(std::pair<int,int>(it->second, it_neighb->second));
286 }
287 }
288 blurs += blur_temp;
289 }
290 }
291 blurs.finalize();
292
293 //bistochastize
294 int maxiter = 10;
295 n = ones_nvertices;
296 m = Eigen::VectorXf::Zero(nvertices);
297 for (int i = 0; i < int(splat_idx.size()); i++)
298 {
299 m(splat_idx[i]) += 1.0f;
300 }
301
302 Eigen::VectorXf bluredn(nvertices);
303
304 for (int i = 0; i < maxiter; i++)
305 {
306 Blur(n,bluredn);
307 n = ((n.array()*m.array()).array()/bluredn.array()).array().sqrt();
308 }
309 Blur(n,bluredn);
310
311 m = n.array() * (bluredn).array();
312 diagonal(m,Dm);
313 diagonal(n,Dn);
314
315 }
316 else
317 {
318 dim = 5;
319 cv::Mat reference_yuv;
320 cv::cvtColor(reference, reference_yuv, COLOR_BGR2YCrCb);
321
322 cols = reference_yuv.cols;
323 rows = reference_yuv.rows;
324 npixels = cols*rows;
325 long long hash_vec[5];
326 for (int i = 0; i < 5; ++i)
327 hash_vec[i] = static_cast<long long>(std::pow(255, i));
328
329 mapId hashed_coords;
330 #if __cplusplus <= 199711L
331 #else
332 hashed_coords.reserve(cols*rows);
333 #endif
334
335 const unsigned char* pref = (const unsigned char*)reference_yuv.data;
336 int vert_idx = 0;
337 int pix_idx = 0;
338
339 // construct Splat(Slice) matrices
340 splat_idx.resize(npixels);
341 for (int y = 0; y < rows; ++y)
342 {
343 for (int x = 0; x < cols; ++x)
344 {
345 long long coord[5];
346 coord[0] = int(x / sigma_spatial);
347 coord[1] = int(y / sigma_spatial);
348 coord[2] = int(pref[0] / sigma_luma);
349 coord[3] = int(pref[1] / sigma_chroma);
350 coord[4] = int(pref[2] / sigma_chroma);
351
352 // convert the coordinate to a hash value
353 long long hash_coord = 0;
354 for (int i = 0; i < 5; ++i)
355 hash_coord += coord[i] * hash_vec[i];
356
357 // pixels whom are alike will have the same hash value.
358 // We only want to keep a unique list of hash values, therefore make sure we only insert
359 // unique hash values.
360 mapId::iterator it = hashed_coords.find(hash_coord);
361 if (it == hashed_coords.end())
362 {
363 hashed_coords.insert(std::pair<long long, int>(hash_coord, vert_idx));
364 splat_idx[pix_idx] = vert_idx;
365 ++vert_idx;
366 }
367 else
368 {
369 splat_idx[pix_idx] = it->second;
370 }
371
372 pref += 3; // skip 3 bytes (y u v)
373 ++pix_idx;
374 }
375 }
376 nvertices = static_cast<int>(hashed_coords.size());
377
378 // construct Blur matrices
379 Eigen::VectorXf ones_nvertices = Eigen::VectorXf::Ones(nvertices);
380 diagonal(ones_nvertices,blurs);
381 blurs *= 10;
382 for(int offset = -1; offset <= 1;++offset)
383 {
384 if(offset == 0) continue;
385 for (int i = 0; i < dim; ++i)
386 {
387 Eigen::SparseMatrix<float, Eigen::ColMajor> blur_temp(hashed_coords.size(), hashed_coords.size());
388 blur_temp.reserve(Eigen::VectorXi::Constant(nvertices,6));
389 long long offset_hash_coord = offset * hash_vec[i];
390 for (mapId::iterator it = hashed_coords.begin(); it != hashed_coords.end(); ++it)
391 {
392 long long neighb_coord = it->first + offset_hash_coord;
393 mapId::iterator it_neighb = hashed_coords.find(neighb_coord);
394 if (it_neighb != hashed_coords.end())
395 {
396 blur_temp.insert(it->second,it_neighb->second) = 1.0f;
397 blur_idx.push_back(std::pair<int,int>(it->second, it_neighb->second));
398 }
399 }
400 blurs += blur_temp;
401 }
402 }
403 blurs.finalize();
404
405
406 //bistochastize
407 int maxiter = 10;
408 n = ones_nvertices;
409 m = Eigen::VectorXf::Zero(nvertices);
410 for (int i = 0; i < int(splat_idx.size()); i++)
411 {
412 m(splat_idx[i]) += 1.0f;
413 }
414
415 Eigen::VectorXf bluredn(nvertices);
416
417 for (int i = 0; i < maxiter; i++)
418 {
419 Blur(n,bluredn);
420 n = ((n.array()*m.array()).array()/bluredn.array()).array().sqrt();
421 }
422 Blur(n,bluredn);
423
424 m = n.array() * (bluredn).array();
425 diagonal(m,Dm);
426 diagonal(n,Dn);
427 }
428 }
429
Splat(Eigen::VectorXf & input,Eigen::VectorXf & output)430 void FastBilateralSolverFilterImpl::Splat(Eigen::VectorXf& input, Eigen::VectorXf& output)
431 {
432 output.setZero();
433 for (int i = 0; i < int(splat_idx.size()); i++)
434 {
435 output(splat_idx[i]) += input(i);
436 }
437
438 }
439
Blur(Eigen::VectorXf & input,Eigen::VectorXf & output)440 void FastBilateralSolverFilterImpl::Blur(Eigen::VectorXf& input, Eigen::VectorXf& output)
441 {
442 output.setZero();
443 output = input * 10;
444 for (int i = 0; i < int(blur_idx.size()); i++)
445 {
446 output(blur_idx[i].first) += input(blur_idx[i].second);
447 }
448 }
449
450
Slice(Eigen::VectorXf & input,Eigen::VectorXf & output)451 void FastBilateralSolverFilterImpl::Slice(Eigen::VectorXf& input, Eigen::VectorXf& output)
452 {
453 output.setZero();
454 for (int i = 0; i < int(splat_idx.size()); i++)
455 {
456 output(i) = input(splat_idx[i]);
457 }
458 }
459
460
solve(cv::Mat & target,cv::Mat & confidence,cv::Mat & output)461 void FastBilateralSolverFilterImpl::solve(cv::Mat& target,
462 cv::Mat& confidence,
463 cv::Mat& output)
464 {
465
466 Eigen::SparseMatrix<float, Eigen::ColMajor> M(nvertices,nvertices);
467 Eigen::SparseMatrix<float, Eigen::ColMajor> A_data(nvertices,nvertices);
468 Eigen::SparseMatrix<float, Eigen::ColMajor> A(nvertices,nvertices);
469 Eigen::VectorXf b(nvertices);
470 Eigen::VectorXf y(nvertices);
471 Eigen::VectorXf y0(nvertices);
472 Eigen::VectorXf y1(nvertices);
473 Eigen::VectorXf w_splat(nvertices);
474
475 Eigen::VectorXf x(npixels);
476 Eigen::VectorXf w(npixels);
477
478 if(target.depth() == CV_16S)
479 {
480 const int16_t *pft = reinterpret_cast<const int16_t*>(target.data);
481 for (int i = 0; i < npixels; i++)
482 {
483 x(i) = (cv::saturate_cast<float>(pft[i])+32768.0f)/65535.0f;
484 }
485 }
486 else if(target.depth() == CV_16U)
487 {
488 const uint16_t *pft = reinterpret_cast<const uint16_t*>(target.data);
489 for (int i = 0; i < npixels; i++)
490 {
491 x(i) = cv::saturate_cast<float>(pft[i])/65535.0f;
492 }
493 }
494 else if(target.depth() == CV_8U)
495 {
496 const uchar *pft = reinterpret_cast<const uchar*>(target.data);
497 for (int i = 0; i < npixels; i++)
498 {
499 x(i) = cv::saturate_cast<float>(pft[i])/255.0f;
500 }
501 }
502 else if(target.depth() == CV_32F)
503 {
504 const float *pft = reinterpret_cast<const float*>(target.data);
505 for (int i = 0; i < npixels; i++)
506 {
507 x(i) = pft[i];
508 }
509 }
510
511 if(confidence.depth() == CV_8U)
512 {
513 const uchar *pfc = reinterpret_cast<const uchar*>(confidence.data);
514 for (int i = 0; i < npixels; i++)
515 {
516 w(i) = cv::saturate_cast<float>(pfc[i])/255.0f;
517 }
518 }
519 else if(confidence.depth() == CV_32F)
520 {
521 const float *pfc = reinterpret_cast<const float*>(confidence.data);
522 for (int i = 0; i < npixels; i++)
523 {
524 w(i) = pfc[i];
525 }
526 }
527
528 //construct A
529 Splat(w,w_splat);
530
531 diagonal(w_splat,A_data);
532 A = bs_param.lam * (Dm - Dn * (blurs*Dn)) + A_data ;
533
534 //construct b
535 b.setZero();
536 for (int i = 0; i < int(splat_idx.size()); i++)
537 {
538 b(splat_idx[i]) += x(i) * w(i);
539 }
540
541 //construct guess for y
542 y0.setZero();
543 for (int i = 0; i < int(splat_idx.size()); i++)
544 {
545 y0(splat_idx[i]) += x(i);
546 }
547 y1.setZero();
548 for (int i = 0; i < int(splat_idx.size()); i++)
549 {
550 y1(splat_idx[i]) += 1.0f;
551 }
552 for (int i = 0; i < nvertices; i++)
553 {
554 y0(i) = y0(i)/y1(i);
555 }
556
557
558 // solve Ay = b
559 Eigen::ConjugateGradient<Eigen::SparseMatrix<float>, Eigen::Lower|Eigen::Upper> cg;
560 cg.compute(A);
561 cg.setMaxIterations(bs_param.cg_maxiter);
562 cg.setTolerance(bs_param.cg_tol);
563 // y = cg.solve(b);
564 y = cg.solveWithGuess(b,y0);
565 std::cout << "#iterations: " << cg.iterations() << std::endl;
566 std::cout << "estimated error: " << cg.error() << std::endl;
567
568 //slice
569 if(target.depth() == CV_16S)
570 {
571 int16_t *pftar = (int16_t*) output.data;
572 for (int i = 0; i < int(splat_idx.size()); i++)
573 {
574 pftar[i] = cv::saturate_cast<short>(y(splat_idx[i]) * 65535.0f - 32768.0f);
575 }
576 }
577 else if(target.depth() == CV_16U)
578 {
579 uint16_t *pftar = (uint16_t*) output.data;
580 for (int i = 0; i < int(splat_idx.size()); i++)
581 {
582 pftar[i] = cv::saturate_cast<ushort>(y(splat_idx[i]) * 65535.0f);
583 }
584 }
585 else if (target.depth() == CV_8U)
586 {
587 uchar *pftar = (uchar*) output.data;
588 for (int i = 0; i < int(splat_idx.size()); i++)
589 {
590 pftar[i] = cv::saturate_cast<uchar>(y(splat_idx[i]) * 255.0f);
591 }
592 }
593 else
594 {
595 float *pftar = (float*)(output.data);
596 for (int i = 0; i < int(splat_idx.size()); i++)
597 {
598 pftar[i] = y(splat_idx[i]);
599 }
600 }
601
602
603 }
604
605
606 ////////////////////////////////////////////////////////////////////////////
607 ////////////////////////////////////////////////////////////////////////////
createFastBilateralSolverFilter(InputArray guide,double sigma_spatial,double sigma_luma,double sigma_chroma,double lambda,int num_iter,double max_tol)608 Ptr<FastBilateralSolverFilter> createFastBilateralSolverFilter(InputArray guide, double sigma_spatial, double sigma_luma, double sigma_chroma, double lambda, int num_iter, double max_tol)
609 {
610 return Ptr<FastBilateralSolverFilter>(FastBilateralSolverFilterImpl::create(guide, sigma_spatial, sigma_luma, sigma_chroma, lambda, num_iter, max_tol));
611 }
612
fastBilateralSolverFilter(InputArray guide,InputArray src,InputArray confidence,OutputArray dst,double sigma_spatial,double sigma_luma,double sigma_chroma,double lambda,int num_iter,double max_tol)613 void fastBilateralSolverFilter(InputArray guide, InputArray src, InputArray confidence, OutputArray dst, double sigma_spatial, double sigma_luma, double sigma_chroma, double lambda, int num_iter, double max_tol)
614 {
615 Ptr<FastBilateralSolverFilter> fbs = createFastBilateralSolverFilter(guide, sigma_spatial, sigma_luma, sigma_chroma, lambda, num_iter, max_tol);
616 fbs->filter(src, confidence, dst);
617 }
618
619 }
620
621 }
622
623 #else
624
625 namespace cv
626 {
627 namespace ximgproc
628 {
629
createFastBilateralSolverFilter(InputArray,double,double,double,double,int,double)630 Ptr<FastBilateralSolverFilter> createFastBilateralSolverFilter(InputArray, double, double, double, double, int, double)
631 {
632 CV_Error(Error::StsNotImplemented, "createFastBilateralSolverFilter : needs to be compiled with EIGEN");
633 }
634
fastBilateralSolverFilter(InputArray,InputArray,InputArray,OutputArray,double,double,double,double,int,double)635 void fastBilateralSolverFilter(InputArray, InputArray, InputArray, OutputArray, double, double, double, double, int, double)
636 {
637 CV_Error(Error::StsNotImplemented, "fastBilateralSolverFilter : needs to be compiled with EIGEN");
638 }
639
640
641 }
642
643 }
644
645 #endif // HAVE_EIGEN
646