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// Intel License Agreement
11// For Open Source Computer Vision Library
12//
13// Copyright (C) 2000, Intel Corporation, all rights reserved.
14// Third party copyrights are property of their respective icvers.
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 Intel Corporation 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/* ////////////////////////////////////////////////////////////////////
43//
44// Geometrical transforms on images and matrices: rotation, zoom etc.
45//
46// */
47
48#include <queue>
49#include <type_traits>
50
51#include "precomp.hpp"
52
53using namespace cv;
54
55#undef CV_MAT_ELEM_PTR_FAST
56#define CV_MAT_ELEM_PTR_FAST( mat, row, col, pix_size ) \
57 ((mat).data.ptr + (size_t)(mat).step*(row) + (pix_size)*(col))
58
59template<typename T>
60typename std::enable_if<std::is_floating_point<T>::value, T>::type round_cast(float val) {
61 return cv::saturate_cast<T>(val);
62}
63
64template<typename T>
65typename std::enable_if<!std::is_floating_point<T>::value, T>::type round_cast(float val) {
66 return cv::saturate_cast<T>(val + 0.5);
67}
68
69inline float
70min4( float a, float b, float c, float d )
71{
72 a = MIN(a,b);
73 c = MIN(c,d);
74 return MIN(a,c);
75}
76
77#define KNOWN 0 //known outside narrow band
78#define BAND 1 //narrow band (known)
79#define INSIDE 2 //unknown
80#define CHANGE 3 //servise
81
82typedef struct CvHeapElem
83{
84 float T;
85 int i,j;
86 int order; // to keep insertion order
87
88 bool operator > (const CvHeapElem& rhs) const {
89 if (T > rhs.T) {
90 return true;
91 } else if (T < rhs.T) {
92 return false;
93 }
94 return order > rhs.order;
95 }
96}
97CvHeapElem;
98
99
100class CvPriorityQueueFloat
101{
102private:
103 CvPriorityQueueFloat(const CvPriorityQueueFloat & ); // copy disabled
104 CvPriorityQueueFloat& operator=(const CvPriorityQueueFloat &); // assign disabled
105
106protected:
107 std::priority_queue<CvHeapElem, std::vector<CvHeapElem>,std::greater<CvHeapElem> > queue;
108 int next_order;
109
110public:
111 bool Add(const Mat &f) {
112 int i,j;
113 for (i=0; i<f.rows; i++) {
114 for (j=0; j<f.cols; j++) {
115 if (f.at<uchar>(i0: i, i1: j)!=0) {
116 if (!Push(i,j,T: 0)) return false;
117 }
118 }
119 }
120 return true;
121 }
122
123 bool Push(int i, int j, float T) {
124 queue.push(x: {.T: T, .i: i, .j: j, .order: next_order});
125 ++next_order;
126 return true;
127 }
128
129 bool Pop(int *i, int *j) {
130 if (queue.empty()) {
131 return false;
132 }
133 *i = queue.top().i;
134 *j = queue.top().j;
135 queue.pop();
136 return true;
137 }
138
139 bool Pop(int *i, int *j, float *T) {
140 if (queue.empty()) {
141 return false;
142 }
143 *i = queue.top().i;
144 *j = queue.top().j;
145 *T = queue.top().T;
146 queue.pop();
147 return true;
148 }
149
150 CvPriorityQueueFloat(void) : queue(), next_order() {
151 }
152};
153
154static inline float VectorScalMult(const cv::Point2f& v1, const cv::Point2f& v2)
155{
156 return v1.x*v2.x+v1.y*v2.y;
157}
158
159static inline float VectorLength(const cv::Point2f& v1)
160{
161 return v1.x*v1.x+v1.y*v1.y;
162}
163
164///////////////////////////////////////////////////////////////////////////////////////////
165//HEAP::iterator Heap_Iterator;
166//HEAP Heap;
167
168static float FastMarching_solve(int i1,int j1,int i2,int j2, const Mat &f, const Mat &t)
169{
170 double sol, a11, a22, m12;
171 a11=t.at<float>(i0: i1,i1: j1);
172 a22=t.at<float>(i0: i2,i1: j2);
173 m12=MIN(a11,a22);
174
175 if( f.at<uchar>(i0: i1,i1: j1) != INSIDE )
176 if( f.at<uchar>(i0: i2,i1: j2) != INSIDE )
177 if( fabs(x: a11-a22) >= 1.0 )
178 sol = 1+m12;
179 else
180 sol = (a11+a22+sqrt(x: (double)(2-(a11-a22)*(a11-a22))))*0.5;
181 else
182 sol = 1+a11;
183 else if( f.at<uchar>(i0: i2,i1: j2) != INSIDE )
184 sol = 1+a22;
185 else
186 sol = 1+m12;
187
188 return (float)sol;
189}
190
191/////////////////////////////////////////////////////////////////////////////////////
192
193
194static void
195icvCalcFMM(Mat &f, Mat &t, CvPriorityQueueFloat *Heap, bool negate) {
196 int i, j, ii = 0, jj = 0, q;
197 float dist;
198
199 while (Heap->Pop(i: &ii,j: &jj)) {
200
201 unsigned known=(negate)?CHANGE:KNOWN;
202 f.at<uchar>(i0: ii,i1: jj) = (uchar)known;
203
204 for (q=0; q<4; q++) {
205 i=0; j=0;
206 if (q==0) {i=ii-1; j=jj;}
207 else if(q==1) {i=ii; j=jj-1;}
208 else if(q==2) {i=ii+1; j=jj;}
209 else {i=ii; j=jj+1;}
210 if ((i<=0)||(j<=0)||(i>f.rows)||(j>f.cols)) continue;
211
212 if (f.at<uchar>(i0: i,i1: j)==INSIDE) {
213 dist = min4(a: FastMarching_solve(i1: i-1,j1: j,i2: i,j2: j-1,f,t),
214 b: FastMarching_solve(i1: i+1,j1: j,i2: i,j2: j-1,f,t),
215 c: FastMarching_solve(i1: i-1,j1: j,i2: i,j2: j+1,f,t),
216 d: FastMarching_solve(i1: i+1,j1: j,i2: i,j2: j+1,f,t));
217 t.at<float>(i0: i,i1: j) = dist;
218 f.at<uchar>(i0: i,i1: j) = BAND;
219 Heap->Push(i,j,T: dist);
220 }
221 }
222 }
223
224 if (negate) {
225 for (i=0; i<f.rows; i++) {
226 for(j=0; j<f.cols; j++) {
227 if (f.at<uchar>(i0: i,i1: j) == CHANGE) {
228 f.at<uchar>(i0: i,i1: j) = KNOWN;
229 t.at<float>(i0: i,i1: j) = -t.at<float>(i0: i,i1: j);
230 }
231 }
232 }
233 }
234}
235
236template <typename data_type>
237static void
238icvTeleaInpaintFMM(Mat &f, Mat &t, Mat &out, int range, CvPriorityQueueFloat *Heap ) {
239 int i = 0, j = 0, ii = 0, jj = 0, k, l, q, color = 0;
240 float dist;
241
242 if (out.channels()==3) {
243 typedef Vec<uchar, 3> PixelT;
244
245 while (Heap->Pop(i: &ii,j: &jj)) {
246
247 f.at<uchar>(i0: ii,i1: jj) = KNOWN;
248 for(q=0; q<4; q++) {
249 if (q==0) {i=ii-1; j=jj;}
250 else if(q==1) {i=ii; j=jj-1;}
251 else if(q==2) {i=ii+1; j=jj;}
252 else if(q==3) {i=ii; j=jj+1;}
253 if ((i<=0)||(j<=0)||(i>t.rows-1)||(j>t.cols-1)) continue;
254
255 if (f.at<uchar>(i0: i,i1: j)==INSIDE) {
256 dist = min4(a: FastMarching_solve(i1: i-1,j1: j,i2: i,j2: j-1,f,t),
257 b: FastMarching_solve(i1: i+1,j1: j,i2: i,j2: j-1,f,t),
258 c: FastMarching_solve(i1: i-1,j1: j,i2: i,j2: j+1,f,t),
259 d: FastMarching_solve(i1: i+1,j1: j,i2: i,j2: j+1,f,t));
260 t.at<float>(i0: i,i1: j) = dist;
261
262 cv::Point2f gradT[3];
263 for (color=0; color<=2; color++) {
264 if (f.at<uchar>(i0: i,i1: j+1)!=INSIDE) {
265 if (f.at<uchar>(i0: i,i1: j-1)!=INSIDE) {
266 gradT[color].x=(float)((t.at<float>(i0: i,i1: j+1)-t.at<float>(i0: i,i1: j-1)))*0.5f;
267 } else {
268 gradT[color].x=(float)((t.at<float>(i0: i,i1: j+1)-t.at<float>(i0: i,i1: j)));
269 }
270 } else {
271 if (f.at<uchar>(i0: i,i1: j-1)!=INSIDE) {
272 gradT[color].x=(float)((t.at<float>(i0: i,i1: j)-t.at<float>(i0: i,i1: j-1)));
273 } else {
274 gradT[color].x=0;
275 }
276 }
277 if (f.at<uchar>(i0: i+1,i1: j)!=INSIDE) {
278 if (f.at<uchar>(i0: i-1,i1: j)!=INSIDE) {
279 gradT[color].y=(float)((t.at<float>(i0: i+1,i1: j)-t.at<float>(i0: i-1,i1: j)))*0.5f;
280 } else {
281 gradT[color].y=(float)((t.at<float>(i0: i+1,i1: j)-t.at<float>(i0: i,i1: j)));
282 }
283 } else {
284 if (f.at<uchar>(i0: i-1,i1: j)!=INSIDE) {
285 gradT[color].y=(float)((t.at<float>(i0: i,i1: j)-t.at<float>(i0: i-1,i1: j)));
286 } else {
287 gradT[color].y=0;
288 }
289 }
290 }
291
292 cv::Point2f gradI,r;
293 float Jx[3] = {0,0,0};
294 float Jy[3] = {0,0,0};
295 float Ia[3] = {0,0,0};
296 float s[3] = {1.0e-20f,1.0e-20f,1.0e-20f};
297 float w,dst,lev,dir,sat;
298
299 for (k=i-range; k<=i+range; k++) {
300 int km=k-1+(k==1),kp=k-1-(k==t.rows-2);
301 for (l=j-range; l<=j+range; l++) {
302 int lm=l-1+(l==1),lp=l-1-(l==t.cols-2);
303 if (k>0&&l>0&&k<t.rows-1&&l<t.cols-1) {
304 if ((f.at<uchar>(i0: k,i1: l)!=INSIDE)&&
305 ((l-j)*(l-j)+(k-i)*(k-i)<=range*range)) {
306 for (color=0; color<=2; color++) {
307 r.y = (float)(i-k);
308 r.x = (float)(j-l);
309
310 dst = (float)(1./(VectorLength(v1: r)*sqrt(x: (double)VectorLength(v1: r))));
311 lev = (float)(1./(1+fabs(x: t.at<float>(i0: k,i1: l)-t.at<float>(i0: i,i1: j))));
312
313 dir=VectorScalMult(v1: r,v2: gradT[color]);
314 if (fabs(x: dir)<=0.01) dir=0.000001f;
315 w = (float)fabs(x: dst*lev*dir);
316
317 if (f.at<uchar>(i0: k,i1: l+1)!=INSIDE) {
318 if (f.at<uchar>(i0: k,i1: l-1)!=INSIDE) {
319 gradI.x=(float)((out.at<PixelT>(i0: km,i1: lp+1)[color]-out.at<PixelT>(i0: km,i1: lm-1)[color]))*2.0f;
320 } else {
321 gradI.x=(float)((out.at<PixelT>(i0: km,i1: lp+1)[color]-out.at<PixelT>(i0: km,i1: lm)[color]));
322 }
323 } else {
324 if (f.at<uchar>(i0: k,i1: l-1)!=INSIDE) {
325 gradI.x=(float)((out.at<PixelT>(i0: km,i1: lp)[color]-out.at<PixelT>(i0: km,i1: lm-1)[color]));
326 } else {
327 gradI.x=0;
328 }
329 }
330 if (f.at<uchar>(i0: k+1,i1: l)!=INSIDE) {
331 if (f.at<uchar>(i0: k-1,i1: l)!=INSIDE) {
332 gradI.y=(float)((out.at<PixelT>(i0: kp+1,i1: lm)[color]-out.at<PixelT>(i0: km-1,i1: lm)[color]))*2.0f;
333 } else {
334 gradI.y=(float)((out.at<PixelT>(i0: kp+1,i1: lm)[color]-out.at<PixelT>(i0: km,i1: lm)[color]));
335 }
336 } else {
337 if (f.at<uchar>(i0: k-1,i1: l)!=INSIDE) {
338 gradI.y=(float)((out.at<PixelT>(i0: kp,i1: lm)[color]-out.at<PixelT>(i0: km-1,i1: lm)[color]));
339 } else {
340 gradI.y=0;
341 }
342 }
343 Ia[color] += (float)w * (float)(out.at<PixelT>(i0: k-1,i1: l-1)[color]);
344 Jx[color] -= (float)w * (float)(gradI.x*r.x);
345 Jy[color] -= (float)w * (float)(gradI.y*r.y);
346 s[color] += w;
347 }
348 }
349 }
350 }
351 }
352 for (color=0; color<=2; color++) {
353 sat = (float)(Ia[color]/s[color]+(Jx[color]+Jy[color])/(sqrt(x: Jx[color]*Jx[color]+Jy[color]*Jy[color])+1.0e-20f));
354 out.at<PixelT>(i0: i-1,i1: j-1)[color] = round_cast<uchar>(val: sat);
355 }
356
357 f.at<uchar>(i0: i,i1: j) = BAND;
358 Heap->Push(i,j,T: dist);
359 }
360 }
361 }
362
363 } else if (out.channels()==1) {
364
365 while (Heap->Pop(i: &ii,j: &jj)) {
366
367 f.at<uchar>(i0: ii,i1: jj) = KNOWN;
368 for(q=0; q<4; q++) {
369 if (q==0) {i=ii-1; j=jj;}
370 else if(q==1) {i=ii; j=jj-1;}
371 else if(q==2) {i=ii+1; j=jj;}
372 else if(q==3) {i=ii; j=jj+1;}
373 if ((i<=0)||(j<=0)||(i>t.rows-1)||(j>t.cols-1)) continue;
374
375 if (f.at<uchar>(i0: i,i1: j)==INSIDE) {
376 dist = min4(a: FastMarching_solve(i1: i-1,j1: j,i2: i,j2: j-1,f,t),
377 b: FastMarching_solve(i1: i+1,j1: j,i2: i,j2: j-1,f,t),
378 c: FastMarching_solve(i1: i-1,j1: j,i2: i,j2: j+1,f,t),
379 d: FastMarching_solve(i1: i+1,j1: j,i2: i,j2: j+1,f,t));
380 t.at<float>(i0: i,i1: j) = dist;
381
382 for (color=0; color<=0; color++) {
383 cv::Point2f gradI,gradT,r;
384 float Ia=0,Jx=0,Jy=0,s=1.0e-20f,w,dst,lev,dir,sat;
385
386 if (f.at<uchar>(i0: i,i1: j+1)!=INSIDE) {
387 if (f.at<uchar>(i0: i,i1: j-1)!=INSIDE) {
388 gradT.x=(float)((t.at<float>(i0: i,i1: j+1)-t.at<float>(i0: i,i1: j-1)))*0.5f;
389 } else {
390 gradT.x=(float)((t.at<float>(i0: i,i1: j+1)-t.at<float>(i0: i,i1: j)));
391 }
392 } else {
393 if (f.at<uchar>(i0: i,i1: j-1)!=INSIDE) {
394 gradT.x=(float)((t.at<float>(i0: i,i1: j)-t.at<float>(i0: i,i1: j-1)));
395 } else {
396 gradT.x=0;
397 }
398 }
399 if (f.at<uchar>(i0: i+1,i1: j)!=INSIDE) {
400 if (f.at<uchar>(i0: i-1,i1: j)!=INSIDE) {
401 gradT.y=(float)((t.at<float>(i0: i+1,i1: j)-t.at<float>(i0: i-1,i1: j)))*0.5f;
402 } else {
403 gradT.y=(float)((t.at<float>(i0: i+1,i1: j)-t.at<float>(i0: i,i1: j)));
404 }
405 } else {
406 if (f.at<uchar>(i0: i-1,i1: j)!=INSIDE) {
407 gradT.y=(float)((t.at<float>(i0: i,i1: j)-t.at<float>(i0: i-1,i1: j)));
408 } else {
409 gradT.y=0;
410 }
411 }
412 for (k=i-range; k<=i+range; k++) {
413 int km=k-1+(k==1),kp=k-1-(k==t.rows-2);
414 for (l=j-range; l<=j+range; l++) {
415 int lm=l-1+(l==1),lp=l-1-(l==t.cols-2);
416 if (k>0&&l>0&&k<t.rows-1&&l<t.cols-1) {
417 if ((f.at<uchar>(i0: k,i1: l)!=INSIDE)&&
418 ((l-j)*(l-j)+(k-i)*(k-i)<=range*range)) {
419 r.y = (float)(i-k);
420 r.x = (float)(j-l);
421
422 dst = (float)(1./(VectorLength(v1: r)*sqrt(x: VectorLength(v1: r))));
423 lev = (float)(1./(1+fabs(x: t.at<float>(i0: k,i1: l)-t.at<float>(i0: i,i1: j))));
424
425 dir=VectorScalMult(v1: r,v2: gradT);
426 if (fabs(x: dir)<=0.01) dir=0.000001f;
427 w = (float)fabs(x: dst*lev*dir);
428
429 if (f.at<uchar>(i0: k,i1: l+1)!=INSIDE) {
430 if (f.at<uchar>(i0: k,i1: l-1)!=INSIDE) {
431 gradI.x=(float)((out.at<data_type>(km,lp+1)-out.at<data_type>(km,lm-1)))*2.0f;
432 } else {
433 gradI.x=(float)((out.at<data_type>(km,lp+1)-out.at<data_type>(km,lm)));
434 }
435 } else {
436 if (f.at<uchar>(i0: k,i1: l-1)!=INSIDE) {
437 gradI.x=(float)((out.at<data_type>(km,lp)-out.at<data_type>(km,lm-1)));
438 } else {
439 gradI.x=0;
440 }
441 }
442 if (f.at<uchar>(i0: k+1,i1: l)!=INSIDE) {
443 if (f.at<uchar>(i0: k-1,i1: l)!=INSIDE) {
444 gradI.y=(float)((out.at<data_type>(kp+1,lm)-out.at<data_type>(km-1,lm)))*2.0f;
445 } else {
446 gradI.y=(float)((out.at<data_type>(kp+1,lm)-out.at<data_type>(km,lm)));
447 }
448 } else {
449 if (f.at<uchar>(i0: k-1,i1: l)!=INSIDE) {
450 gradI.y=(float)((out.at<data_type>(kp,lm)-out.at<data_type>(km-1,lm)));
451 } else {
452 gradI.y=0;
453 }
454 }
455 Ia += (float)w * (float)(out.at<data_type>(k-1,l-1));
456 Jx -= (float)w * (float)(gradI.x*r.x);
457 Jy -= (float)w * (float)(gradI.y*r.y);
458 s += w;
459 }
460 }
461 }
462 }
463 sat = (float)(Ia/s+(Jx+Jy)/(sqrt(x: Jx*Jx+Jy*Jy)+1.0e-20f));
464 {
465 out.at<data_type>(i-1,j-1) = round_cast<data_type>(sat);
466 }
467 }
468
469 f.at<uchar>(i0: i,i1: j) = BAND;
470 Heap->Push(i,j,T: dist);
471 }
472 }
473 }
474 }
475}
476
477template <typename data_type>
478static void
479icvNSInpaintFMM(Mat &f, Mat &t, Mat &out, int range, CvPriorityQueueFloat *Heap) {
480 int i = 0, j = 0, ii = 0, jj = 0, k, l, q, color = 0;
481 float dist;
482
483 if (out.channels()==3) {
484 typedef Vec<uchar, 3> PixelT;
485
486 while (Heap->Pop(i: &ii,j: &jj)) {
487
488 f.at<uchar>(i0: ii,i1: jj) = KNOWN;
489 for(q=0; q<4; q++) {
490 if (q==0) {i=ii-1; j=jj;}
491 else if(q==1) {i=ii; j=jj-1;}
492 else if(q==2) {i=ii+1; j=jj;}
493 else if(q==3) {i=ii; j=jj+1;}
494 if ((i<=0)||(j<=0)||(i>t.rows-1)||(j>t.cols-1)) continue;
495
496 if (f.at<uchar>(i0: i,i1: j)==INSIDE) {
497 dist = min4(a: FastMarching_solve(i1: i-1,j1: j,i2: i,j2: j-1,f,t),
498 b: FastMarching_solve(i1: i+1,j1: j,i2: i,j2: j-1,f,t),
499 c: FastMarching_solve(i1: i-1,j1: j,i2: i,j2: j+1,f,t),
500 d: FastMarching_solve(i1: i+1,j1: j,i2: i,j2: j+1,f,t));
501 t.at<float>(i0: i,i1: j) = dist;
502
503 cv::Point2f gradI,r;
504 float Ia[3]={0,0,0};
505 float s[3]={1.0e-20f,1.0e-20f,1.0e-20f};
506 float w,dst,dir;
507
508 for (k=i-range; k<=i+range; k++) {
509 int km=k-1+(k==1),kp=k-1-(k==f.rows-2);
510 for (l=j-range; l<=j+range; l++) {
511 int lm=l-1+(l==1),lp=l-1-(l==f.cols-2);
512 if (k>0&&l>0&&k<f.rows-1&&l<f.cols-1) {
513 if ((f.at<uchar>(i0: k,i1: l)!=INSIDE)&&
514 ((l-j)*(l-j)+(k-i)*(k-i)<=range*range)) {
515 for (color=0; color<=2; color++) {
516 r.y=(float)(k-i);
517 r.x=(float)(l-j);
518
519 dst = 1/(VectorLength(v1: r)*VectorLength(v1: r)+1);
520
521 if (f.at<uchar>(i0: k+1,i1: l)!=INSIDE) {
522 if (f.at<uchar>(i0: k-1,i1: l)!=INSIDE) {
523 gradI.x=(float)(abs(x: out.at<PixelT>(i0: kp+1,i1: lm)[color]-out.at<PixelT>(i0: kp,i1: lm)[color])+
524 abs(x: out.at<PixelT>(i0: kp,i1: lm)[color]-out.at<PixelT>(i0: km-1,i1: lm)[color]));
525 } else {
526 gradI.x=(float)(abs(x: out.at<PixelT>(i0: kp+1,i1: lm)[color]-out.at<PixelT>(i0: kp,i1: lm)[color]))*2.0f;
527 }
528 } else {
529 if (f.at<uchar>(i0: k-1,i1: l)!=INSIDE) {
530 gradI.x=(float)(abs(x: out.at<PixelT>(i0: kp,i1: lm)[color]-out.at<PixelT>(i0: km-1,i1: lm)[color]))*2.0f;
531 } else {
532 gradI.x=0;
533 }
534 }
535 if (f.at<uchar>(i0: k,i1: l+1)!=INSIDE) {
536 if (f.at<uchar>(i0: k,i1: l-1)!=INSIDE) {
537 gradI.y=(float)(abs(x: out.at<PixelT>(i0: km,i1: lp+1)[color]-out.at<PixelT>(i0: km,i1: lm)[color])+
538 abs(x: out.at<PixelT>(i0: km,i1: lm)[color]-out.at<PixelT>(i0: km,i1: lm-1)[color]));
539 } else {
540 gradI.y=(float)(abs(x: out.at<PixelT>(i0: km,i1: lp+1)[color]-out.at<PixelT>(i0: km,i1: lm)[color]))*2.0f;
541 }
542 } else {
543 if (f.at<uchar>(i0: k,i1: l-1)!=INSIDE) {
544 gradI.y=(float)(abs(x: out.at<PixelT>(i0: km,i1: lm)[color]-out.at<PixelT>(i0: km,i1: lm-1)[color]))*2.0f;
545 } else {
546 gradI.y=0;
547 }
548 }
549
550 gradI.x=-gradI.x;
551 dir=VectorScalMult(v1: r,v2: gradI);
552
553 if (fabs(x: dir)<=0.01) {
554 dir=0.000001f;
555 } else {
556 dir = (float)fabs(x: VectorScalMult(v1: r,v2: gradI)/sqrt(x: VectorLength(v1: r)*VectorLength(v1: gradI)));
557 }
558 w = dst*dir;
559 Ia[color] += (float)w * (float)(out.at<PixelT>(i0: k-1,i1: l-1)[color]);
560 s[color] += w;
561 }
562 }
563 }
564 }
565 }
566 for (color=0; color<=2; color++) {
567 out.at<PixelT>(i0: i-1,i1: j-1)[color] = cv::saturate_cast<uchar>(v: (double)Ia[color]/s[color]);
568 }
569
570 f.at<uchar>(i0: i,i1: j) = BAND;
571 Heap->Push(i,j,T: dist);
572 }
573 }
574 }
575
576 } else if (out.channels()==1) {
577
578 while (Heap->Pop(i: &ii,j: &jj)) {
579
580 f.at<uchar>(i0: ii,i1: jj) = KNOWN;
581 for(q=0; q<4; q++) {
582 if (q==0) {i=ii-1; j=jj;}
583 else if(q==1) {i=ii; j=jj-1;}
584 else if(q==2) {i=ii+1; j=jj;}
585 else if(q==3) {i=ii; j=jj+1;}
586 if ((i<=0)||(j<=0)||(i>t.rows-1)||(j>t.cols-1)) continue;
587
588 if (f.at<uchar>(i0: i,i1: j)==INSIDE) {
589 dist = min4(a: FastMarching_solve(i1: i-1,j1: j,i2: i,j2: j-1,f,t),
590 b: FastMarching_solve(i1: i+1,j1: j,i2: i,j2: j-1,f,t),
591 c: FastMarching_solve(i1: i-1,j1: j,i2: i,j2: j+1,f,t),
592 d: FastMarching_solve(i1: i+1,j1: j,i2: i,j2: j+1,f,t));
593 t.at<float>(i0: i,i1: j) = dist;
594
595 {
596 cv::Point2f gradI,r;
597 float Ia=0,s=1.0e-20f,w,dst,dir;
598
599 for (k=i-range; k<=i+range; k++) {
600 int km=k-1+(k==1),kp=k-1-(k==t.rows-2);
601 for (l=j-range; l<=j+range; l++) {
602 int lm=l-1+(l==1),lp=l-1-(l==t.cols-2);
603 if (k>0&&l>0&&k<t.rows-1&&l<t.cols-1) {
604 if ((f.at<uchar>(i0: k,i1: l)!=INSIDE)&&
605 ((l-j)*(l-j)+(k-i)*(k-i)<=range*range)) {
606 r.y=(float)(i-k);
607 r.x=(float)(j-l);
608
609 dst = 1/(VectorLength(v1: r)*VectorLength(v1: r)+1);
610
611 if (f.at<uchar>(i0: k+1,i1: l)!=INSIDE) {
612 if (f.at<uchar>(i0: k-1,i1: l)!=INSIDE) {
613 gradI.x=(float)(std::abs(out.at<data_type>(kp+1,lm)-out.at<data_type>(kp,lm))+
614 std::abs(out.at<data_type>(kp,lm)-out.at<data_type>(km-1,lm)));
615 } else {
616 gradI.x=(float)(std::abs(out.at<data_type>(kp+1,lm)-out.at<data_type>(kp,lm)))*2.0f;
617 }
618 } else {
619 if (f.at<uchar>(i0: k-1,i1: l)!=INSIDE) {
620 gradI.x=(float)(std::abs(out.at<data_type>(kp,lm)-out.at<data_type>(km-1,lm)))*2.0f;
621 } else {
622 gradI.x=0;
623 }
624 }
625 if (f.at<uchar>(i0: k,i1: l+1)!=INSIDE) {
626 if (f.at<uchar>(i0: k,i1: l-1)!=INSIDE) {
627 gradI.y=(float)(std::abs(out.at<data_type>(km,lp+1)-out.at<data_type>(km,lm))+
628 std::abs(out.at<data_type>(km,lm)-out.at<data_type>(km,lm-1)));
629 } else {
630 gradI.y=(float)(std::abs(out.at<data_type>(km,lp+1)-out.at<data_type>(km,lm)))*2.0f;
631 }
632 } else {
633 if (f.at<uchar>(i0: k,i1: l-1)!=INSIDE) {
634 gradI.y=(float)(std::abs(out.at<data_type>(km,lm)-out.at<data_type>(km,lm-1)))*2.0f;
635 } else {
636 gradI.y=0;
637 }
638 }
639
640 gradI.x=-gradI.x;
641 dir=VectorScalMult(v1: r,v2: gradI);
642
643 if (fabs(x: dir)<=0.01) {
644 dir=0.000001f;
645 } else {
646 dir = (float)fabs(x: VectorScalMult(v1: r,v2: gradI)/sqrt(x: VectorLength(v1: r)*VectorLength(v1: gradI)));
647 }
648 w = dst*dir;
649 Ia += (float)w * (float)(out.at<data_type>(k-1,l-1));
650 s += w;
651 }
652 }
653 }
654 }
655 out.at<data_type>(i-1,j-1) = cv::saturate_cast<data_type>((double)Ia/s);
656 }
657
658 f.at<uchar>(i0: i,i1: j) = BAND;
659 Heap->Push(i,j,T: dist);
660 }
661 }
662 }
663
664 }
665}
666
667#define SET_BORDER1_C1(image,type,value) {\
668 int i,j;\
669 for(j=0; j<image.cols; j++) {\
670 image.at<type>(0,j) = value;\
671 }\
672 for (i=1; i<image.rows-1; i++) {\
673 image.at<type>(i,0) = image.at<type>(i,image.cols-1) = value;\
674 }\
675 for(j=0; j<image.cols; j++) {\
676 image.at<type>(erows-1,j) = value;\
677 }\
678 }
679
680#define COPY_MASK_BORDER1_C1(src,dst,type) {\
681 int i,j;\
682 for (i=0; i<src.rows; i++) {\
683 for(j=0; j<src.cols; j++) {\
684 if (src.at<type>(i,j)!=0)\
685 dst.at<type>(i+1,j+1) = INSIDE;\
686 }\
687 }\
688 }
689
690static void
691icvInpaint( const Mat &input_img, const Mat &inpaint_mask, Mat &output_img,
692 double inpaintRange, int flags )
693{
694 cv::Mat mask, band, f, t, out;
695 cv::Ptr<CvPriorityQueueFloat> Heap, Out;
696 cv::Mat el_range, el_cross; // structuring elements for dilate
697
698 int range=cvRound(value: inpaintRange);
699 int erows, ecols;
700
701 if((input_img.size() != output_img.size()) || (input_img.size() != inpaint_mask.size()))
702 CV_Error( cv::Error::StsUnmatchedSizes, "All the input and output images must have the same size" );
703
704 if( (input_img.type() != CV_8U &&
705 input_img.type() != CV_16U &&
706 input_img.type() != CV_32F &&
707 input_img.type() != CV_8UC3) ||
708 (input_img.type() != output_img.type()) )
709 CV_Error( cv::Error::StsUnsupportedFormat,
710 "8-bit, 16-bit unsigned or 32-bit float 1-channel and 8-bit 3-channel input/output images are supported" );
711
712 if( inpaint_mask.type() != CV_8UC1 )
713 CV_Error( cv::Error::StsUnsupportedFormat, "The mask must be 8-bit 1-channel image" );
714
715 range = MAX(range,1);
716 range = MIN(range,100);
717
718 ecols = input_img.cols + 2;
719 erows = input_img.rows + 2;
720
721 f.create(rows: erows, cols: ecols, CV_8UC1);
722 t.create(rows: erows, cols: ecols, CV_32FC1);
723 band.create(rows: erows, cols: ecols, CV_8UC1);
724 mask.create(rows: erows, cols: ecols, CV_8UC1);
725 el_cross = cv::getStructuringElement(shape: cv::MORPH_CROSS, ksize: cv::Size(3, 3), anchor: cv::Point(1, 1));
726
727 input_img.copyTo( m: output_img );
728 mask.setTo(value: Scalar(KNOWN,0,0,0));
729 COPY_MASK_BORDER1_C1(inpaint_mask,mask,uchar);
730 SET_BORDER1_C1(mask,uchar,0);
731 f.setTo(value: Scalar(KNOWN,0,0,0));
732 t.setTo(value: Scalar(1.0e6f,0,0,0));
733 cv::dilate(src: mask, dst: band, kernel: el_cross, anchor: cv::Point(1, 1));
734 Heap=cv::makePtr<CvPriorityQueueFloat>();
735 subtract(src1: band, src2: mask, dst: band);
736 SET_BORDER1_C1(band,uchar,0);
737 if (!Heap->Add(f: band))
738 return;
739
740 f.setTo(value: Scalar(BAND,0,0,0),mask: band);
741 f.setTo(value: Scalar(INSIDE,0,0,0),mask);
742 t.setTo(value: Scalar(0,0,0,0),mask: band);
743
744 if( flags == cv::INPAINT_TELEA )
745 {
746 out.create(rows: erows, cols: ecols, CV_8UC1);
747 el_range = cv::getStructuringElement(shape: cv::MORPH_RECT, ksize: cv::Size(2 * range + 1, 2 * range + 1));
748 cv::dilate(src: mask, dst: out, kernel: el_range);
749 subtract(src1: out, src2: mask, dst: out);
750 Out=cv::makePtr<CvPriorityQueueFloat>();
751 if (!Out->Add(f: band))
752 return;
753 subtract(src1: out, src2: band, dst: out);
754 SET_BORDER1_C1(out,uchar,0);
755 icvCalcFMM(f&: out,t,Heap: Out,negate: true);
756 switch(output_img.depth())
757 {
758 case CV_8U:
759 icvTeleaInpaintFMM<uchar>(f&: mask,t,out&: output_img,range,Heap);
760 break;
761 case CV_16U:
762 icvTeleaInpaintFMM<ushort>(f&: mask,t,out&: output_img,range,Heap);
763 break;
764 case CV_32F:
765 icvTeleaInpaintFMM<float>(f&: mask,t,out&: output_img,range,Heap);
766 break;
767 default:
768 CV_Error( cv::Error::StsBadArg, "Unsupportedformat of the input image" );
769 }
770 }
771 else if (flags == cv::INPAINT_NS) {
772 switch(output_img.depth())
773 {
774 case CV_8U:
775 icvNSInpaintFMM<uchar>(f&: mask,t,out&: output_img,range,Heap);
776 break;
777 case CV_16U:
778 icvNSInpaintFMM<ushort>(f&: mask,t,out&: output_img,range,Heap);
779 break;
780 case CV_32F:
781 icvNSInpaintFMM<float>(f&: mask,t,out&: output_img,range,Heap);
782 break;
783 default:
784 CV_Error( cv::Error::StsBadArg, "Unsupported format of the input image" );
785 }
786 } else {
787 CV_Error( cv::Error::StsBadArg, "The flags argument must be one of INPAINT_TELEA or INPAINT_NS" );
788 }
789}
790
791void cv::inpaint( InputArray _src, InputArray _mask, OutputArray _dst,
792 double inpaintRange, int flags )
793{
794 CV_INSTRUMENT_REGION();
795
796 Mat src = _src.getMat(), mask = _mask.getMat();
797 _dst.create( sz: src.size(), type: src.type() );
798 Mat dst = _dst.getMat();
799 icvInpaint( input_img: src, inpaint_mask: mask, output_img&: dst, inpaintRange, flags );
800}
801

source code of opencv/modules/photo/src/inpaint.cpp