| static void Impala::Core::Array::f_iir_tline_filter | ( | DSTTYPE * | src, | |
| DSTTYPE * | dest, | |||
| int | sx, | |||
| int | sy, | |||
| double * | filter, | |||
| double | tanp | |||
| ) | [inline, static] |
Definition at line 216 of file AniGauss.h.
References TriggsM().
Referenced by anigauss().
00218 { 00219 double *p0, *p1, *p2, *p3; 00220 double *buf0, *buf1, *buf2, *buf3; 00221 double *uplusbuf; 00222 int i, j; 00223 double b1, b2, b3; 00224 double sum, sumsq; 00225 double uplus, vplus; 00226 double unp, unp1, unp2; 00227 double M[9]; 00228 double pix, prev, val; 00229 double res, prevres; 00230 double xf; 00231 int x; 00232 double c, d; 00233 double e, f; 00234 int q4 = 0; 00235 00236 /* check filter direction towards first or fourth quadrant */ 00237 if (tanp <= 0.0) { 00238 q4 = 1; 00239 tanp = -tanp; 00240 } 00241 00242 /* alloc buffer for Triggs boundary condition */ 00243 uplusbuf = (double*) malloc(sx*sizeof(*uplusbuf)); 00244 00245 /* alloc recursion line buffers */ 00246 buf0 = (double*) malloc((sx+sy*tanp+2)*sizeof(*buf0)); 00247 buf1 = (double*) malloc((sx+sy*tanp+2)*sizeof(*buf1)); 00248 buf2 = (double*) malloc((sx+sy*tanp+2)*sizeof(*buf2)); 00249 buf3 = (double*) malloc((sx+sy*tanp+2)*sizeof(*buf3)); 00250 00251 if (q4) { 00252 buf0 += (int)(sy*tanp+1); 00253 buf1 += (int)(sy*tanp+1); 00254 buf2 += (int)(sy*tanp+1); 00255 buf3 += (int)(sy*tanp+1); 00256 } 00257 00258 sumsq = filter[3]; 00259 sum = sumsq*sumsq; 00260 00261 /* causal filter */ 00262 b1 = filter[2]; b2 = filter[1]; b3 = filter[0]; 00263 00264 /* first line */ 00265 p1 = buf1; p2 = buf2; p3 = buf3; 00266 for (j = 0; j < sx; j++) { 00267 pix = *src++; 00268 *dest++ = pix; *p1++ = pix; *p2++ = pix; *p3++ = pix; 00269 } 00270 00271 /* calc last line for Triggs boundary condition */ 00272 src += (sy-2)*sx; 00273 for (j = 0; j < sx; j++) 00274 uplusbuf[j] = *src++ * sumsq/(1.0-b1-b2-b3); 00275 src -= (sy-1)*sx; 00276 00277 x = 0; 00278 for (i = 1; i < sy; i++) { 00279 xf = i*tanp; 00280 00281 /* border handling at image corner */ 00282 if (q4) { 00283 p1 = buf1-x; p2 = buf2-x; p3 = buf3-x; 00284 for (j=1; j <= (int)(xf)-x; j++) { 00285 p1[-j] = p1[0]; 00286 p2[-j] = p2[0]; 00287 p3[-j] = p3[0]; 00288 } 00289 } 00290 else { 00291 p1 = buf1+x; p2 = buf2+x; p3 = buf3+x; 00292 for (j=1; j <= (int)(xf)-x; j++) { 00293 p1[sx+j-1] = p1[sx-1]; 00294 p2[sx+j-1] = p2[sx-1]; 00295 p3[sx+j-1] = p3[sx-1]; 00296 } 00297 } 00298 00299 /* calc interpolation coefficients */ 00300 x = (int)xf; 00301 c = xf-(double)x; 00302 d = 1.0-c; 00303 00304 e = c; f = d; 00305 if (!q4) { 00306 res = d; d = c; c = res; 00307 res = f; f = e; e = res; 00308 } 00309 00310 c *= sumsq; d *= sumsq; 00311 00312 /* set buffers at start */ 00313 if (q4) { 00314 p0 = buf0-x; p1 = buf1-x; p2 = buf2-x; p3 = buf3-x; 00315 } 00316 else { 00317 p0 = buf0+x; p1 = buf1+x; p2 = buf2+x; p3 = buf3+x; 00318 } 00319 prev = *src; 00320 prevres = sumsq*prev + b1 * *p1 + b2 * *p2 + b3 * *p3; 00321 00322 /* run the filter */ 00323 for (j = 0; j < sx; j++) { 00324 pix = *src++; 00325 val = c*pix+d*prev; 00326 prev = pix; 00327 00328 res = val + b1 * *p1++ + b2 * *p2++ + b3 * *p3++; 00329 *p0++ = res; 00330 *dest++ = f*res+e*prevres; 00331 prevres = res; 00332 } 00333 00334 /* shift history */ 00335 p0 = buf3; buf3 = buf2; buf2 = buf1; buf1 = buf0; buf0 = p0; 00336 } 00337 00338 /* anti-causal */ 00339 00340 /* apply Triggs border condition */ 00341 b1 = filter[4]; b2 = filter[5]; b3 = filter[6]; 00342 TriggsM(filter, M); 00343 00344 /* first line */ 00345 x = (int)((sy-1)*tanp); 00346 if (q4) { 00347 p1 = buf1+sx-x; p2 = buf2+sx-x; p3 = buf3+sx-x; 00348 } 00349 else { 00350 p1 = buf1+sx+x; p2 = buf2+sx+x; p3 = buf3+sx+x; 00351 } 00352 p0 = uplusbuf+sx; 00353 for (j = 0; j < sx; j++) { 00354 uplus = *(--p0); 00355 vplus = uplus/(1.0-b1-b2-b3); 00356 00357 unp = *(--p1)-uplus; 00358 unp1 = *(--p2)-uplus; 00359 unp2 = *(--p3)-uplus; 00360 pix = M[0]*unp+M[1]*unp1+M[2]*unp2 + vplus; 00361 pix *= sumsq; 00362 *(--dest) = pix; *p1 = pix; 00363 pix = M[3]*unp+M[4]*unp1+M[5]*unp2 + vplus; 00364 *p2 = pix*sumsq; 00365 pix = M[6]*unp+M[7]*unp1+M[8]*unp2 + vplus; 00366 *p3 = pix*sumsq; 00367 } 00368 00369 for (i = sy-2; i >= 0; i--) { 00370 xf = i*tanp; 00371 00372 /* border handling at image corner */ 00373 if (q4) { 00374 p1 = buf1-x; p2 = buf2-x; p3 = buf3-x; 00375 for (j=1; j <= x-(int)(xf); j++) { 00376 p1[sx+j-1] = p1[sx-1]; 00377 p2[sx+j-1] = p2[sx-1]; 00378 p3[sx+j-1] = p3[sx-1]; 00379 } 00380 } 00381 else { 00382 p1 = buf1+x; p2 = buf2+x; p3 = buf3+x; 00383 for (j=1; j <= x-(int)(xf); j++) { 00384 p1[-j] = p1[0]; 00385 p2[-j] = p2[0]; 00386 p3[-j] = p3[0]; 00387 } 00388 } 00389 00390 /* calc interpolation coefficients */ 00391 x = (int)xf; 00392 c = xf-(double)x; 00393 d = 1.0-c; 00394 00395 e = c; f = d; 00396 c *= sumsq; d *= sumsq; 00397 00398 if (!q4) { 00399 res = d; d = c; c = res; 00400 res = f; f = e; e = res; 00401 } 00402 00403 /* set buffers at start */ 00404 if (q4) { 00405 p0 = buf0+sx-x; p1 = buf1+sx-x; p2 = buf2+sx-x; p3 = buf3+sx-x; 00406 } 00407 else { 00408 p0 = buf0+sx+x; p1 = buf1+sx+x; p2 = buf2+sx+x; p3 = buf3+sx+x; 00409 } 00410 prev = *(dest-1); 00411 prevres = sumsq*prev + b1 * *(p1-1) + b2 * *(p2-1) + b3 * *(p3-1); 00412 00413 /* run the filter */ 00414 for (j = 0; j < sx; j++) { 00415 pix = *(--dest); 00416 val = d*pix+c*prev; 00417 prev = pix; 00418 00419 res = val + b1 * *(--p1) + b2 * *(--p2) + b3 * *(--p3); 00420 *(--p0) = res; 00421 *dest = e*res+f*prevres; 00422 prevres = res; 00423 } 00424 00425 /* shift history */ 00426 p0 = buf3; buf3 = buf2; buf2 = buf1; buf1 = buf0; buf0 = p0; 00427 } 00428 00429 if (q4) { 00430 buf0 -= (int)(sy*tanp+1); 00431 buf1 -= (int)(sy*tanp+1); 00432 buf2 -= (int)(sy*tanp+1); 00433 buf3 -= (int)(sy*tanp+1); 00434 } 00435 free(buf0); 00436 free(buf1); 00437 free(buf2); 00438 free(buf3); 00439 free(uplusbuf); 00440 }
Here is the call graph for this function:
