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