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