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 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: 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 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 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 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 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 189 CvPriorityQueueFloat(void) { 190 num=in=0; 191 mem=empty=head=tail=NULL; 192 } 193 194 ~CvPriorityQueueFloat(void) 195 { 196 cvFree( &mem ); 197 } 198 }; 199 200 inline float VectorScalMult(CvPoint2D32f v1,CvPoint2D32f v2) { 201 return v1.x*v2.x+v1.y*v2.y; 202 } 203 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 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 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 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 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 { 721 template<> void cv::DefaultDeleter<IplConvKernel>::operator ()(IplConvKernel* obj) const 722 { 723 cvReleaseStructuringElement(&obj); 724 } 725 } 726 727 void 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 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
