1 // This file is part of OpenCV project.
2 // It is subject to the license terms in the LICENSE file found in the top-level directory
3 // of this distribution and at http://opencv.org/license.html.
4 
5 // This is not a standalone header, see inpainting.cpp
6 
7 namespace cv
8 {
9 namespace xphoto
10 {
11 
12 struct fsr_parameters
13 {
14     // default variables
15     int block_size = 16;
16     double conc_weighting = 0.5;
17     double rhos[4] = { 0.80, 0.70, 0.66, 0.64 };
18     double threshold_stddev_Y[3] = { 0.014, 0.030, 0.090 };
19     double threshold_stddev_Cx[3] = { 0.006, 0.010, 0.028 };
20     // quality profile dependent variables
21     int block_size_min, fft_size, max_iter, min_iter, iter_const;
22     double orthogonality_correction;
fsr_parameterscv::xphoto::fsr_parameters23     fsr_parameters(const int quality)
24     {
25         if (quality == xphoto::INPAINT_FSR_BEST)
26         {
27             block_size_min = 2;
28             fft_size = 64;
29             max_iter = 400;
30             min_iter = 50;
31             iter_const = 2000;
32             orthogonality_correction = 0.2;
33         }
34         else if (quality == xphoto::INPAINT_FSR_FAST)
35         {
36             block_size_min = 4;
37             fft_size = 32;
38             max_iter = 100;
39             min_iter = 20;
40             iter_const = 1000;
41             orthogonality_correction = 0.5;
42         }
43         else
44         {
45             CV_Error(Error::StsBadArg, "Unknown quality level set, supported: FAST, BEST");
46 
47         }
48     }
49 };
50 
51 
52 static void
icvSgnMat(const Mat & src,Mat & dst)53 icvSgnMat(const Mat& src, Mat& dst) {
54     dst = Mat::zeros(src.size(), CV_64F);
55     for (int y = 0; y < src.rows; ++y)
56     {
57         for (int x = 0; x < src.cols; ++x)
58         {
59             double curr_val = src.at<double>(y,x);
60             if (curr_val > 0)
61             {
62                 dst.at<double>(y,x) = 1;
63             }
64             else if (curr_val)
65             {
66                 dst.at<double>(y,x) = -1;
67             }
68         }
69     }
70 }
71 
72 
73 static double
icvStandardDeviation(const Mat & distorted_block_2d,const Mat & error_mask_2d)74 icvStandardDeviation(const Mat& distorted_block_2d, const Mat& error_mask_2d) {
75     if (countNonZero(error_mask_2d) < 1)
76     {
77         return NAN; // block with no undistorted pixels shouldn't be chosen for processing (only if block_size_min is reached)
78     }
79     Scalar tmp_stddev, tmp_mean;
80     Mat mask8u;
81     error_mask_2d.convertTo(mask8u, CV_8U, 2.0);
82     meanStdDev(distorted_block_2d, tmp_mean, tmp_stddev, mask8u);
83     double sigma_n = tmp_stddev[0] / 255;
84     if (sigma_n < 0)
85     {
86         sigma_n = 0;
87     }
88     else if (sigma_n > 1)
89     {
90         sigma_n = 1;
91     }
92     return sigma_n;
93 }
94 
95 static void
icvExtrapolateBlock(Mat & distorted_block,Mat & error_mask,fsr_parameters & fsr_params,double rho,double normedStdDev,Mat & extrapolated_block)96 icvExtrapolateBlock(Mat& distorted_block, Mat& error_mask, fsr_parameters& fsr_params, double rho, double normedStdDev, Mat& extrapolated_block)
97 {
98     double fft_size = fsr_params.fft_size;
99     double orthogonality_correction = fsr_params.orthogonality_correction;
100     int M = distorted_block.rows;
101     int N = distorted_block.cols;
102     int fft_x_offset = cvFloor((fft_size - N) / 2);
103     int fft_y_offset = cvFloor((fft_size - M) / 2);
104 
105     // weighting function
106     Mat w = Mat::zeros(fsr_params.fft_size, fsr_params.fft_size, CV_64F);
107     error_mask.copyTo(w(Range(fft_y_offset, fft_y_offset + M), Range(fft_x_offset, fft_x_offset + N)));
108     for (int u = 0; u < fft_size; ++u)
109     {
110         for (int v = 0; v < fft_size; ++v)
111         {
112             w.at<double>(u, v) *= std::pow(rho, std::sqrt(std::pow(u + 0.5 - (fft_y_offset + M / 2), 2) + std::pow(v + 0.5 - (fft_x_offset + N / 2), 2)));
113         }
114     }
115     Mat W;
116     dft(w, W, DFT_COMPLEX_OUTPUT);
117     Mat W_padded;
118     hconcat(W, W, W_padded);
119     vconcat(W_padded, W_padded, W_padded);
120 
121     // frequency weighting
122     Mat frequency_weighting = Mat::ones(fsr_params.fft_size, fsr_params.fft_size / 2 + 1, CV_64F);
123     for (int y = 0; y < fft_size; ++y)
124     {
125         for (int x = 0; x < (fft_size / 2 + 1); ++x)
126         {
127             double y2 = fft_size / 2 - std::abs(y - fft_size / 2);
128             double x2 = fft_size / 2 - std::abs(x - fft_size / 2);
129             frequency_weighting.at<double>(y, x) = 1 - std::sqrt(x2*x2 + y2 * y2)*std::sqrt(2) / fft_size;
130         }
131     }
132     // pad image to fft window size
133     Mat f(Size(fsr_params.fft_size, fsr_params.fft_size), CV_64F, Scalar::all(0));
134     distorted_block.copyTo(f(Range(fft_y_offset, fft_y_offset + M), Range(fft_x_offset, fft_x_offset + N)));
135 
136     // create initial model
137     Mat G = Mat::zeros(fsr_params.fft_size, fsr_params.fft_size, CV_64FC2); // complex
138 
139     // calculate initial residual
140     Mat Rw_tmp, Rw;
141     dft(f.mul(w), Rw_tmp, DFT_COMPLEX_OUTPUT);
142     Rw = Rw_tmp(Range(0, fsr_params.fft_size), Range(0, fsr_params.fft_size / 2 + 1));
143 
144     // estimate ideal number of iterations (GenserIWSSIP2017)
145     // calculate stddev if not available (e.g., for smallest block size)
146     if (normedStdDev == 0) {
147         normedStdDev = icvStandardDeviation(distorted_block, error_mask);
148     }
149     int num_iters = cvRound(fsr_params.iter_const * normedStdDev);
150     if (num_iters < fsr_params.min_iter) {
151         num_iters = fsr_params.min_iter;
152     }
153     else if (num_iters > fsr_params.max_iter) {
154         num_iters = fsr_params.max_iter;
155     }
156 
157     int iter_counter = 0;
158     while (iter_counter < num_iters)
159     { // Spectral Constrained FSE (GenserIWSSIP2018)
160         Mat projection_distances(Rw.size(), CV_64F);
161         Mat Rw_mag = Mat(Rw.size(), CV_64F);
162         std::vector<Mat> channels(2);
163         split(Rw, channels);
164         magnitude(channels[0], channels[1], Rw_mag);
165         projection_distances = Rw_mag.mul(frequency_weighting);
166 
167         double minVal, maxVal;
168         int maxLocx = -1;
169         int maxLocy = -1;
170         minMaxLoc(projection_distances, &minVal, &maxVal);
171 
172         for (int y = 0; y < projection_distances.rows; ++y)
173         { // assure that first appearance of max Value is selected
174             for (int x = 0; x < projection_distances.cols; ++x)
175             {
176                 if (std::abs(projection_distances.at<double>(y, x) - maxVal) < 0.001)
177                 {
178                     maxLocy = y;
179                     maxLocx = x;
180                     break;
181                 }
182             }
183             if (maxLocy != -1)
184             {
185                 break;
186             }
187         }
188         int bf2select = maxLocy + maxLocx * projection_distances.rows;
189         int v = static_cast<int>(std::max(0.0, std::floor(bf2select / fft_size)));
190         int u = static_cast<int>(std::max(0, bf2select % fsr_params.fft_size));
191 
192 
193         // exclude second half of first and middle col
194         if ((v == 0 && u > fft_size / 2) || (v == fft_size / 2 && u > fft_size / 2))
195         {
196             int u_prev = u;
197             u = fsr_params.fft_size - u;
198             Rw.at<std::complex<double> >(u, v) = std::conj(Rw.at<std::complex<double> >(u_prev, v));
199         }
200 
201         // calculate complex conjugate solution
202         int u_cj = -1;
203         int v_cj = -1;
204         // fill first lower col (copy from first upper col)
205         if (u >= 1 && u < fft_size / 2 && v == 0)
206         {
207             u_cj = fsr_params.fft_size - u;
208             v_cj = v;
209         }
210         // fill middle lower col (copy from first middle col)
211         if (u >= 1 && u < fft_size / 2 && v == fft_size / 2)
212         {
213             u_cj = fsr_params.fft_size - u;
214             v_cj = v;
215         }
216         // fill first row right (copy from first row left)
217         if (u == 0 && v >= 1 && v < fft_size / 2)
218         {
219             u_cj = u;
220             v_cj = fsr_params.fft_size - v;
221         }
222         // fill middle row right (copy from middle row left)
223         if (u == fft_size / 2 && v >= 1 && v < fft_size / 2)
224         {
225             u_cj = u;
226             v_cj = fsr_params.fft_size - v;
227         }
228         // fill cell upper right (copy from lower cell left)
229         if (u >= fft_size / 2 + 1 && v >= 1 && v < fft_size / 2)
230         {
231             u_cj = fsr_params.fft_size - u;
232             v_cj = fsr_params.fft_size - v;
233         }
234         // fill cell lower right (copy from upper cell left)
235         if (u >= 1 && u < fft_size / 2 && v >= 1 && v < fft_size / 2)
236         {
237             u_cj = fsr_params.fft_size - u;
238             v_cj = fsr_params.fft_size - v;
239         }
240 
241         /// add coef to model and update residual
242         if (u_cj != -1 && v_cj != -1)
243         {
244             std::complex< double> expansion_coefficient = orthogonality_correction * Rw.at< std::complex<double> >(u, v) / W.at<std::complex<double> >(0, 0);
245             G.at< std::complex<double> >(u, v) += fft_size * fft_size * expansion_coefficient;
246             G.at< std::complex<double> >(u_cj, v_cj) = std::conj(G.at< std::complex<double> >(u, v));
247 
248             Mat expansion_mat(Rw.size(), CV_64FC2, Scalar(expansion_coefficient.real(), expansion_coefficient.imag()));
249             Mat W_tmp1 = W_padded(Range(fsr_params.fft_size - u, fsr_params.fft_size - u + Rw.rows), Range(fsr_params.fft_size - v, fsr_params.fft_size - v + Rw.cols));
250             Mat W_tmp2 = W_padded(Range(fsr_params.fft_size - u_cj, fsr_params.fft_size - u_cj + Rw.rows), Range(fsr_params.fft_size - v_cj, fsr_params.fft_size - v_cj + Rw.cols));
251             Mat res_1(W_tmp1.size(), W_tmp1.type());
252             mulSpectrums(expansion_mat, W_tmp1, res_1, 0);
253             expansion_mat.setTo(Scalar(expansion_coefficient.real(), -expansion_coefficient.imag()));
254             Mat res_2(W_tmp1.size(), W_tmp1.type());
255             mulSpectrums(expansion_mat, W_tmp2, res_2, 0);
256             Rw -= res_1 + res_2;
257 
258             ++iter_counter; // ... as two basis functions were added
259         }
260         else
261         {
262             std::complex<double> expansion_coefficient = orthogonality_correction * Rw.at< std::complex<double> >(u, v) / W.at< std::complex<double> >(0, 0);
263             G.at< std::complex<double> >(u, v) += fft_size * fft_size * expansion_coefficient;
264             Mat expansion_mat(Rw.size(), CV_64FC2, Scalar(expansion_coefficient.real(), expansion_coefficient.imag()));
265             Mat W_tmp = W_padded(Range(fsr_params.fft_size - u, fsr_params.fft_size - u + Rw.rows), Range(fsr_params.fft_size - v, fsr_params.fft_size - v + Rw.cols));
266             Mat res_tmp(W_tmp.size(), W_tmp.type());
267             mulSpectrums(expansion_mat, W_tmp, res_tmp, 0);
268             Rw -= res_tmp;
269 
270         }
271         ++iter_counter;
272     }
273 
274     // get pixels from model
275     Mat g;
276     idft(G, g, DFT_SCALE);
277 
278     // extract reconstructed pixels
279     Mat g_real(M, N, CV_64F);
280     for (int x = 0; x < M; ++x)
281     {
282         for (int y = 0; y < N; ++y)
283         {
284             g_real.at<double>(x, y) = g.at< std::complex<double> >(fft_y_offset + x, fft_x_offset + y).real();
285         }
286     }
287     g_real.copyTo(extrapolated_block);
288     Mat orig_samples;
289     error_mask.convertTo(orig_samples, CV_8U);
290     distorted_block.copyTo(extrapolated_block, orig_samples); // copy where orig_samples is nonzero
291 }
292 
293 
294 static void
icvGetTodoBlocks(Mat & sampled_img,Mat & sampling_mask,std::vector<std::tuple<int,int>> & set_todo,int block_size,int block_size_min,int border_width,double homo_threshold,Mat & set_process_this_block_size,std::vector<std::tuple<int,int>> & set_later,Mat & sigma_n_array)295 icvGetTodoBlocks(Mat& sampled_img, Mat& sampling_mask, std::vector< std::tuple< int, int > >& set_todo, int block_size, int block_size_min, int border_width, double homo_threshold, Mat& set_process_this_block_size, std::vector< std::tuple< int, int > >& set_later, Mat& sigma_n_array)
296 {
297     std::vector< std::tuple< int, int > > set_now;
298     set_later.clear();
299     size_t list_length = set_todo.size();
300     int img_height = sampled_img.rows;
301     int img_width = sampled_img.cols;
302     Mat reconstructed_img;
303     sampled_img.copyTo(reconstructed_img);
304 
305     // calculate block lists
306     for (size_t entry = 0; entry < list_length; ++entry)
307     {
308         int xblock_counter = std::get<0>(set_todo[entry]);
309         int yblock_counter = std::get<1>(set_todo[entry]);
310 
311         int left_border = std::min(xblock_counter*block_size, border_width);
312         int top_border = std::min(yblock_counter*block_size, border_width);
313         int right_border = std::max(0, std::min(img_width - (xblock_counter + 1)*block_size, border_width));
314         int bottom_border = std::max(0, std::min(img_height - (yblock_counter + 1)*block_size, border_width));
315 
316         // extract blocks from images
317         Mat distorted_block_2d = reconstructed_img(Range(yblock_counter*block_size - top_border, std::min(img_height, (yblock_counter*block_size + block_size + bottom_border))), Range(xblock_counter*block_size - left_border, std::min(img_width, (xblock_counter*block_size + block_size + right_border))));
318         Mat error_mask_2d = sampling_mask(Range(yblock_counter*block_size - top_border, std::min(img_height, (yblock_counter*block_size + block_size + bottom_border))), Range(xblock_counter*block_size - left_border, std::min(img_width, (xblock_counter*block_size + block_size + right_border))));
319 
320         // determine normalized and weighted standard deviation
321         if (block_size > block_size_min)
322         {
323             double sigma_n = icvStandardDeviation(distorted_block_2d, error_mask_2d);
324             sigma_n_array.at<double>( yblock_counter, xblock_counter) = sigma_n;
325 
326             // homogeneous case
327             if (sigma_n < homo_threshold)
328             {
329                 set_now.emplace_back(xblock_counter, yblock_counter);
330                 set_process_this_block_size.at<double>(yblock_counter, xblock_counter) = 255;
331 
332             }
333             else
334             {
335                 int yblock_counter_quadernary = yblock_counter * 2;
336                 int xblock_counter_quadernary = xblock_counter * 2;
337                 int yblock_offset = 0;
338                 int xblock_offset = 0;
339 
340                 for (int quader_counter = 0; quader_counter < 4; ++quader_counter)
341                 {
342                     if (quader_counter == 0)
343                     {
344                         yblock_offset = 0;
345                         xblock_offset = 0;
346                     }
347                     else if (quader_counter == 1)
348                     {
349                         yblock_offset = 0;
350                         xblock_offset = 1;
351                     }
352                     else if (quader_counter == 2)
353                     {
354                         yblock_offset = 1;
355                         xblock_offset = 0;
356                     }
357                     else if (quader_counter == 3)
358                     {
359                         yblock_offset = 1;
360                         xblock_offset = 1;
361                     }
362 
363                     set_later.emplace_back(xblock_counter_quadernary + xblock_offset, yblock_counter_quadernary + yblock_offset);
364                 }
365 
366             }
367         }
368 
369     }
370 }
371 
372 
373 static void
icvDetermineProcessingOrder(const Mat & _sampled_img,const Mat & _sampling_mask,const int quality,const std::string & channel,Mat & reconstructed_img)374 icvDetermineProcessingOrder(
375     const Mat& _sampled_img, const Mat& _sampling_mask,
376     const int quality, const std::string& channel, Mat& reconstructed_img
377 )
378 {
379     fsr_parameters fsr_params(quality);
380     int block_size = fsr_params.block_size;
381     int block_size_max = fsr_params.block_size;
382     int block_size_min = fsr_params.block_size_min;
383     double conc_weighting = fsr_params.conc_weighting;
384     int fft_size = fsr_params.fft_size;
385     double rho = fsr_params.rhos[0];
386     Mat sampled_img, sampling_mask;
387     _sampled_img.convertTo(sampled_img, CV_64F);
388     reconstructed_img = sampled_img.clone();
389 
390     _sampling_mask.convertTo(sampling_mask, CV_64F);
391 
392     double threshold_stddev_LUT[3];
393     if (channel == "Y")
394     {
395         std::copy(fsr_params.threshold_stddev_Y, fsr_params.threshold_stddev_Y + 3, threshold_stddev_LUT);
396     }
397     else if (channel == "Cx")
398     {
399         std::copy(fsr_params.threshold_stddev_Cx, fsr_params.threshold_stddev_Cx + 3, threshold_stddev_LUT);
400     }
401     else
402     {
403         CV_Error(Error::StsBadArg, "channel type unsupported!");
404     }
405 
406 
407     double threshold_stddev = threshold_stddev_LUT[0];
408 
409     std::vector< std::tuple< int, int > > set_later;
410     int img_height = sampled_img.rows;
411     int img_width = sampled_img.cols;
412 
413     // initial scan of distorted blocks
414     std::vector< std::tuple< int, int > > set_todo;
415     int blocks_column = divUp(img_height, block_size);
416     int blocks_line = divUp(img_width, block_size);
417     for (int y = 0; y < blocks_column; ++y)
418     {
419         for (int x = 0; x < blocks_line; ++x)
420         {
421             Mat curr_block = sampling_mask(Range(y*block_size, std::min(img_height, (y + 1)*block_size)), Range(x*block_size, std::min(img_width, (x + 1)*block_size)));
422             double min_block, max_block;
423             minMaxLoc(curr_block, &min_block, &max_block);
424             if (min_block == 0)
425             {
426                 set_todo.emplace_back(x, y);
427             }
428         }
429     }
430 
431     // loop over all distorted blocks and extrapolate them depending on
432     // their block size
433     int border_width = 0;
434     while (block_size >= block_size_min)
435     {
436         int blocks_per_column = cvCeil(img_height / block_size);
437         int blocks_per_line = cvCeil(img_width / block_size);
438         Mat nen_array = Mat::zeros(blocks_per_column, blocks_per_line, CV_64F);
439         Mat proc_array = Mat::zeros(blocks_per_column, blocks_per_line, CV_64F);
440         Mat sigma_n_array = Mat::zeros(blocks_per_column, blocks_per_line, CV_64F);
441         Mat set_process_this_block_size = Mat::zeros(blocks_per_column, blocks_per_line, CV_64F);
442         if (block_size > block_size_min)
443         {
444             if (block_size < block_size_max)
445             {
446                 set_todo = set_later;
447             }
448             border_width = cvFloor(fft_size - block_size) / 2;
449             icvGetTodoBlocks(sampled_img, sampling_mask, set_todo, block_size, block_size_min, border_width, threshold_stddev, set_process_this_block_size, set_later, sigma_n_array);
450         }
451         else
452         {
453             set_process_this_block_size.setTo(Scalar(255));
454         }
455 
456         // if block to be extrapolated, increase nen of neighboring pixels
457         for (int yblock_counter = 0; yblock_counter < blocks_per_column; ++yblock_counter)
458         {
459             for (int xblock_counter = 0; xblock_counter < blocks_per_line; ++xblock_counter)
460             {
461                 Mat curr_block = sampling_mask(Range(yblock_counter*block_size, std::min(img_height, (yblock_counter + 1)*block_size)), Range(xblock_counter*block_size, std::min(img_width, (xblock_counter + 1)*block_size)));
462                 double min_block, max_block;
463                 minMaxLoc(curr_block, &min_block, &max_block);
464                 if (min_block == 0)
465                 {
466                     if (yblock_counter > 0 && xblock_counter > 0)
467                     {
468                         nen_array.at<double>(yblock_counter - 1, xblock_counter - 1)++;
469                     }
470                     if (yblock_counter > 0)
471                     {
472                         nen_array.at<double>(yblock_counter - 1, xblock_counter)++;
473                     }
474                     if (yblock_counter > 0 && xblock_counter < (blocks_per_line - 1))
475                     {
476                         nen_array.at<double>(yblock_counter - 1, xblock_counter + 1)++;
477                     }
478                     if (xblock_counter > 0)
479                     {
480                         nen_array.at<double>(yblock_counter, xblock_counter - 1)++;
481                     }
482                     if (xblock_counter < (blocks_per_line - 1))
483                     {
484                         nen_array.at<double>(yblock_counter, xblock_counter + 1)++;
485                     }
486                     if (yblock_counter < (blocks_per_column - 1) && xblock_counter>0)
487                     {
488                         nen_array.at<double>(yblock_counter + 1, xblock_counter - 1)++;
489                     }
490                     if (yblock_counter < (blocks_per_column - 1))
491                     {
492                         nen_array.at<double>(yblock_counter + 1, xblock_counter)++;
493                     }
494                     if (yblock_counter < (blocks_per_column - 1) && xblock_counter < (blocks_per_line - 1))
495                     {
496                         nen_array.at<double>(yblock_counter + 1, xblock_counter + 1)++;
497                     }
498                 }
499             }
500         }
501 
502         // determine if block itself has to be extrapolated
503         for (int yblock_counter = 0; yblock_counter < blocks_per_column; ++yblock_counter)
504         {
505             for (int xblock_counter = 0; xblock_counter < blocks_per_line; ++xblock_counter)
506             {
507                 Mat curr_block = sampling_mask(Range(yblock_counter*block_size, std::min(img_height, (yblock_counter + 1)*block_size)), Range(xblock_counter*block_size, std::min(img_width, (xblock_counter + 1)*block_size)));
508                 double min_block, max_block;
509                 minMaxLoc(curr_block, &min_block, &max_block);
510                 if (min_block != 0)
511                 {
512                     nen_array.at<double>(yblock_counter, xblock_counter) = -1;
513                 }
514                 else
515                 {
516                 // if border block, increase nen respectively
517                     if (yblock_counter == 0 && xblock_counter == 0)
518                     {
519                         nen_array.at<double>(yblock_counter, xblock_counter) = nen_array.at<double>(yblock_counter, xblock_counter) + 5;
520                     }
521                     if (yblock_counter == 0 && xblock_counter == (blocks_per_line - 1))
522                     {
523                         nen_array.at<double>(yblock_counter, xblock_counter) = nen_array.at<double>(yblock_counter, xblock_counter) + 5;
524                     }
525                     if (yblock_counter == (blocks_per_column - 1) && xblock_counter == 0)
526                     {
527                         nen_array.at<double>(yblock_counter, xblock_counter) = nen_array.at<double>(yblock_counter, xblock_counter) + 5;
528                     }
529                     if (yblock_counter == (blocks_per_column - 1) && xblock_counter == (blocks_per_line - 1))
530                     {
531                         nen_array.at<double>(yblock_counter, xblock_counter) = nen_array.at<double>(yblock_counter, xblock_counter) + 5;
532                     }
533                     if (yblock_counter == 0 && xblock_counter != 0 && xblock_counter != (blocks_per_line - 1))
534                     {
535                         nen_array.at<double>(yblock_counter, xblock_counter) = nen_array.at<double>(yblock_counter, xblock_counter) + 3;
536                     }
537                     if (yblock_counter == (blocks_per_column - 1) && xblock_counter != 0 && xblock_counter != (blocks_per_line - 1))
538                     {
539                         nen_array.at<double>(yblock_counter, xblock_counter) = nen_array.at<double>(yblock_counter, xblock_counter) + 3;
540                     }
541                     if (yblock_counter != 0 && yblock_counter != (blocks_per_column - 1) && xblock_counter == 0)
542                     {
543                         nen_array.at<double>(yblock_counter, xblock_counter) = nen_array.at<double>(yblock_counter, xblock_counter) + 3;
544                     }
545                     if (yblock_counter != 0 && yblock_counter != (blocks_per_column - 1) && xblock_counter == (blocks_per_line - 1))
546                     {
547                         nen_array.at<double>(yblock_counter, xblock_counter) = nen_array.at<double>(yblock_counter, xblock_counter) + 3;
548                     }
549                 }
550             }
551         }
552 
553         // if all blocks have 8 not extrapolated neighbors, penalize nen of blocks without any known samples by one
554         double min_nen_tmp, max_nen_tmp;
555         minMaxLoc(nen_array, &min_nen_tmp, &max_nen_tmp);
556         if (min_nen_tmp == 8) {
557             for (int yblock_counter = 0; yblock_counter < blocks_per_column; ++yblock_counter)
558             {
559                 for (int xblock_counter = 0; xblock_counter < blocks_per_line; ++xblock_counter)
560                 {
561                     Mat curr_block = sampling_mask(Range(yblock_counter*block_size, std::min(img_height, (yblock_counter + 1)*block_size)), Range(xblock_counter*block_size, std::min(img_width, (xblock_counter + 1)*block_size)));
562                     double min_block, max_block;
563                     minMaxLoc(curr_block, &min_block, &max_block);
564                     if (max_block == 0)
565                     {
566                         nen_array.at<double>(yblock_counter, xblock_counter)++;
567                     }
568                 }
569             }
570         }
571 
572         // do actual processing per block
573         int all_blocks_finished = 0;
574         while (all_blocks_finished == 0) {
575             // clear proc_array
576             proc_array.setTo(Scalar(1));
577 
578             // determine blocks to extrapolate
579             double min_nen = 99;
580             int bl_counter = 0;
581             // add all homogeneous blocks that shall be processed to list
582             // using same priority
583             // begins with highest prioroty or lowest nen array value
584             std::vector< std::tuple< int, int > > block_list;
585             for (int yblock_counter = 0; yblock_counter < blocks_per_column; ++yblock_counter)
586             {
587                 for (int xblock_counter = 0; xblock_counter < blocks_per_line; ++xblock_counter)
588                 {
589             // decision if block contains errors
590                     double tmp_val = nen_array.at<double>(yblock_counter, xblock_counter);
591                     if (tmp_val >= 0 && tmp_val < min_nen && set_process_this_block_size.at<double>(yblock_counter, xblock_counter) == 255) {
592                         bl_counter = 0;
593                         block_list.clear();
594                         min_nen = tmp_val;
595                         proc_array.setTo(Scalar(1));
596                     }
597                     if (tmp_val == min_nen && proc_array.at<double>(yblock_counter, xblock_counter) != 0 && set_process_this_block_size.at<double>(yblock_counter, xblock_counter) == 0) {
598                         nen_array.at<double>(yblock_counter, xblock_counter) = -1;
599                     }
600                     if (tmp_val == min_nen && proc_array.at<double>(yblock_counter, xblock_counter) != 0 && set_process_this_block_size.at<double>(yblock_counter, xblock_counter) != 0) {
601                         block_list.emplace_back(yblock_counter, xblock_counter);
602                         bl_counter++;
603                         // block neighboring blocks from processing
604                         if (yblock_counter > 0 && xblock_counter > 0)
605                         {
606                             proc_array.at<double>(yblock_counter - 1, xblock_counter - 1) = 0;
607                         }
608                         if (yblock_counter > 0)
609                         {
610                             proc_array.at<double>(yblock_counter - 1, xblock_counter) = 0;
611                         }
612                          if (yblock_counter > 0 && xblock_counter > 0)
613                         {
614                             proc_array.at<double>(yblock_counter - 1, xblock_counter - 1) = 0;
615                         }
616                         if (yblock_counter > 0)
617                         {
618                             proc_array.at<double>(yblock_counter - 1, xblock_counter) = 0;
619                         }
620                         if (yblock_counter > 0 && xblock_counter < (blocks_per_line - 1))
621                         {
622                             proc_array.at<double>(yblock_counter - 1, xblock_counter + 1) = 0;
623                         }
624                         if (xblock_counter > 0)
625                         {
626                             proc_array.at<double>(yblock_counter, xblock_counter - 1) = 0;
627                         }
628                         if (xblock_counter < (blocks_per_line - 1))
629                         {
630                             proc_array.at<double>(yblock_counter, xblock_counter + 1) = 0;
631                         }
632                         if (yblock_counter < (blocks_per_column - 1) && xblock_counter > 0)
633                         {
634                             proc_array.at<double>(yblock_counter + 1, xblock_counter - 1) = 0;
635                         }
636                         if (yblock_counter < (blocks_per_column - 1))
637                         {
638                             proc_array.at<double>(yblock_counter + 1, xblock_counter) = 0;
639                         }
640                         if (yblock_counter < (blocks_per_column - 1) && xblock_counter < (blocks_per_line - 1))
641                         {
642                             proc_array.at<double>(yblock_counter + 1, xblock_counter + 1) = 0;
643                         }
644                     }
645                 }
646             }
647             int max_bl_counter = bl_counter;
648             block_list.emplace_back(-1, -1);
649             if (bl_counter == 0)
650             {
651                 all_blocks_finished = 1;
652             }
653             // blockwise extrapolation of all blocks that can be processed in parallel
654             for (bl_counter = 0; bl_counter < max_bl_counter; ++bl_counter)
655             {
656                 int yblock_counter = std::get<0>(block_list[bl_counter]);
657                 int xblock_counter = std::get<1>(block_list[bl_counter]);
658 
659                 // calculation of the extrapolation area's borders
660                 int left_border = std::min(xblock_counter*block_size, border_width);
661                 int top_border = std::min(yblock_counter*block_size, border_width);
662                 int right_border = std::max(0, std::min(img_width - (xblock_counter + 1)*block_size, border_width));
663                 int bottom_border = std::max(0, std::min(img_height - (yblock_counter + 1)*block_size, border_width));
664 
665                 // extract blocks from images
666                 Mat distorted_block_2d = reconstructed_img(Range(yblock_counter*block_size - top_border, std::min(img_height, (yblock_counter*block_size + block_size + bottom_border))), Range(xblock_counter*block_size - left_border, std::min(img_width, (xblock_counter*block_size + block_size + right_border))));
667                 Mat error_mask_2d = sampling_mask(Range(yblock_counter*block_size - top_border, std::min(img_height, (yblock_counter*block_size + block_size + bottom_border))), Range(xblock_counter*block_size - left_border, std::min(img_width, xblock_counter*block_size + block_size + right_border)));
668                 // get actual stddev value as it is needed to estimate the
669                 // best number of iterations
670                 double sigma_n_a = sigma_n_array.at<double>(yblock_counter, xblock_counter);
671 
672                 // actual extrapolation
673                 Mat extrapolated_block_2d;
674                 icvExtrapolateBlock(distorted_block_2d, error_mask_2d, fsr_params, rho, sigma_n_a, extrapolated_block_2d);
675 
676                 // update image and mask
677                 extrapolated_block_2d(Range(top_border, extrapolated_block_2d.rows - bottom_border), Range(left_border, extrapolated_block_2d.cols - right_border)).copyTo(reconstructed_img(Range(yblock_counter*block_size, std::min(img_height, (yblock_counter + 1)*block_size)), Range(xblock_counter*block_size, std::min(img_width, (xblock_counter + 1)*block_size))));
678 
679                 Mat signs;
680                 icvSgnMat(error_mask_2d(Range(top_border, error_mask_2d.rows - bottom_border), Range(left_border, error_mask_2d.cols - right_border)), signs);
681                 Mat tmp_mask = error_mask_2d(Range(top_border, error_mask_2d.rows - bottom_border), Range(left_border, error_mask_2d.cols - right_border)) + (1 - signs) *conc_weighting;
682                 tmp_mask.copyTo(sampling_mask(Range(yblock_counter*block_size, std::min(img_height, (yblock_counter + 1)*block_size)), Range(xblock_counter*block_size, std::min(img_width, (xblock_counter + 1)*block_size))));
683 
684                 // update nen-array
685                 nen_array.at<double>(yblock_counter, xblock_counter) = -1;
686                 if (yblock_counter > 0 && xblock_counter > 0)
687                 {
688                     nen_array.at<double>(yblock_counter - 1, xblock_counter - 1)--;
689                 }
690                 if (yblock_counter > 0)
691                 {
692                     nen_array.at<double>(yblock_counter - 1, xblock_counter)--;
693                 }
694                 if (yblock_counter > 0 && xblock_counter < blocks_per_line - 1)
695                 {
696                     nen_array.at<double>(yblock_counter - 1, xblock_counter + 1)--;
697                 }
698                 if (xblock_counter > 0)
699                 {
700                     nen_array.at<double>(yblock_counter, xblock_counter - 1)--;
701                 }
702                 if (xblock_counter < blocks_per_line - 1)
703                 {
704                     nen_array.at<double>(yblock_counter, xblock_counter + 1)--;
705                 }
706                 if (yblock_counter < blocks_per_column - 1 && xblock_counter>0)
707                 {
708                     nen_array.at<double>(yblock_counter + 1, xblock_counter - 1)--;
709                 }
710                 if (yblock_counter < blocks_per_column - 1)
711                 {
712                     nen_array.at<double>(yblock_counter + 1, xblock_counter)--;
713                 }
714                 if (yblock_counter < blocks_per_column - 1 && xblock_counter < blocks_per_line - 1)
715                 {
716                     nen_array.at<double>(yblock_counter + 1, xblock_counter + 1)--;
717                 }
718 
719             }
720 
721         }
722 
723         // set parameters for next extrapolation tasks (higher texture)
724         block_size = block_size / 2;
725         border_width = (fft_size - block_size) / 2;
726         if (block_size == 8)
727         {
728             threshold_stddev = threshold_stddev_LUT[1];
729             rho = fsr_params.rhos[1];
730         }
731         if (block_size == 4)
732         {
733             threshold_stddev = threshold_stddev_LUT[2];
734             rho = fsr_params.rhos[2];
735         }
736         if (block_size == 2)
737         {
738             rho = fsr_params.rhos[3];
739         }
740 
741         // terminate function - no heterogeneous blocks left
742         if (set_later.empty())
743         {
744             break;
745         }
746     }
747 }
748 
749 
750 static
inpaint_fsr(const Mat & src,const Mat & mask,Mat & dst,const int algorithmType)751 void inpaint_fsr(const Mat &src, const Mat &mask, Mat &dst, const int algorithmType)
752 {
753     CV_Assert(algorithmType == xphoto::INPAINT_FSR_BEST || algorithmType == xphoto::INPAINT_FSR_FAST);
754     CV_Check(src.channels(), src.channels() == 1 || src.channels() == 3, "");
755     switch (src.type())
756     {
757         case CV_8UC1:
758         case CV_8UC3:
759             break;
760         case CV_16UC1:
761         case CV_16UC3:
762         {
763             double minRange, maxRange;
764             minMaxLoc(src, &minRange, &maxRange);
765             if (minRange < 0 || maxRange > 65535)
766             {
767                 CV_Error(Error::StsUnsupportedFormat, "Unsupported source image format!");
768                 break;
769             }
770             src.convertTo(src, CV_8U, 1/257.0);
771             break;
772         }
773         case CV_32FC1:
774         case CV_64FC1:
775         case CV_32FC3:
776         case CV_64FC3:
777         {
778             double minRange, maxRange;
779             minMaxLoc(src, &minRange, &maxRange);
780             if (minRange < -FLT_EPSILON || maxRange > (1.0 + FLT_EPSILON))
781             {
782                 CV_Error(Error::StsUnsupportedFormat, "Unsupported source image format!");
783                 break;
784             }
785             src.convertTo(src, CV_8U, 255.0);
786             break;
787         }
788         default:
789             CV_Error(Error::StsUnsupportedFormat, "Unsupported source image format!");
790             break;
791     }
792     dst.create(src.size(), src.type());
793     Mat mask_01;
794     threshold(mask, mask_01, 0.0, 1.0, THRESH_BINARY);
795     if (src.channels() == 1)
796     { // grayscale image
797         Mat y_reconstructed;
798         icvDetermineProcessingOrder(src, mask_01, algorithmType, "Y", y_reconstructed);
799         y_reconstructed.convertTo(dst, CV_8U);
800     }
801     else if (src.channels() == 3)
802     { // RGB image
803         Mat ycrcb;
804         cvtColor(src, ycrcb, COLOR_BGR2YCrCb);
805         std::vector<Mat> channels(3);
806         split(ycrcb, channels);
807         Mat y = channels[0];
808         Mat cb = channels[2];
809         Mat cr = channels[1];
810         Mat y_reconstructed, cb_reconstructed, cr_reconstructed;
811         y = y.mul(mask_01);
812         cb = cb.mul(mask_01);
813         cr = cr.mul(mask_01);
814         icvDetermineProcessingOrder(y, mask_01, algorithmType, "Y", y_reconstructed);
815         icvDetermineProcessingOrder(cb, mask_01, algorithmType, "Cx", cb_reconstructed);
816         icvDetermineProcessingOrder(cr, mask_01, algorithmType, "Cx", cr_reconstructed);
817         Mat ycrcb_reconstructed;
818         y_reconstructed.convertTo(channels[0], CV_8U);
819         cr_reconstructed.convertTo(channels[1], CV_8U);
820         cb_reconstructed.convertTo(channels[2], CV_8U);
821         merge(channels, ycrcb_reconstructed);
822         cvtColor(ycrcb_reconstructed, dst, COLOR_YCrCb2BGR);
823     }
824 }
825 
826 }}  // namespace
827