Milena (Olena)
User documentation 2.0a Id
|
00001 // Copyright (C) 2009 EPITA Research and Development Laboratory (LRDE) 00002 // 00003 // This file is part of Olena. 00004 // 00005 // Olena is free software: you can redistribute it and/or modify it under 00006 // the terms of the GNU General Public License as published by the Free 00007 // Software Foundation, version 2 of the License. 00008 // 00009 // Olena is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // General Public License for more details. 00013 // 00014 // You should have received a copy of the GNU General Public License 00015 // along with Olena. If not, see <http://www.gnu.org/licenses/>. 00016 // 00017 // As a special exception, you may use this file as part of a free 00018 // software project without restriction. Specifically, if other files 00019 // instantiate templates or use macros or inline functions from this 00020 // file, or you compile this file and link it with other files to produce 00021 // an executable, this file does not by itself cause the resulting 00022 // executable to be covered by the GNU General Public License. This 00023 // exception does not however invalidate any other reasons why the 00024 // executable file might be covered by the GNU General Public License. 00025 00026 #ifndef MLN_LINEAR_GAUSSIAN_DIRECTIONAL_2D_HH 00027 # define MLN_LINEAR_GAUSSIAN_DIRECTIONAL_2D_HH 00028 00034 00035 #include <mln/core/image/image2d.hh> 00036 #include <mln/extension/adjust_fill.hh> 00037 #include <mln/extension/adjust_duplicate.hh> 00038 00039 00040 00041 namespace mln 00042 { 00043 00044 namespace linear 00045 { 00046 00047 00048 template <typename I> 00049 mln_concrete(I) 00050 gaussian_directional_2d(const Image<I>& input, 00051 unsigned dir, double sigma, 00052 const mln_value(I)& bdr); 00053 00054 00055 00056 # ifndef MLN_INCLUDE_ONLY 00057 00058 namespace my 00059 { 00060 00061 struct recursivefilter_coef_ 00062 { 00063 enum FilterType { DericheGaussian, 00064 DericheGaussianFirstDerivative, 00065 DericheGaussianSecondDerivative }; 00066 00067 std::vector<double> n, d, nm, dm; 00068 double sumA, sumC; 00069 00070 recursivefilter_coef_(double a0, double a1, 00071 double b0, double b1, 00072 double c0, double c1, 00073 double w0, double w1, 00074 double s, FilterType filter_type) 00075 { 00076 n.reserve(5); 00077 d.reserve(5); 00078 nm.reserve(5); 00079 dm.reserve(5); 00080 00081 b0 /= s; 00082 b1 /= s; 00083 w0 /= s; 00084 w1 /= s; 00085 00086 double sin0 = sin(w0); 00087 double sin1 = sin(w1); 00088 double cos0 = cos(w0); 00089 double cos1 = cos(w1); 00090 00091 switch (filter_type) { 00092 00093 case DericheGaussian : 00094 { 00095 sumA = 00096 (2.0 * a1 * exp( b0 ) * cos0 * cos0 - a0 * sin0 * exp( 2.0 * b0 ) 00097 + a0 * sin0 - 2.0 * a1 * exp( b0 )) / 00098 (( 2.0 * cos0 * exp( b0 ) - exp( 2.0 * b0 ) - 1 ) * sin0); 00099 00100 sumC = 00101 (2.0 * c1 * exp( b1 ) * cos1 * cos1 - c0 * sin1 * exp( 2.0 * b1 ) 00102 + c0 * sin1 - 2.0 * c1 * exp( b1 )) 00103 / (( 2.0 * cos1 * exp( b1 ) - exp( 2.0 * b1 ) - 1 ) * sin1); 00104 break; 00105 } 00106 00107 case DericheGaussianFirstDerivative : 00108 { 00109 sumA = -2.f * 00110 (a0 * cos0 - a1 * sin0 + a1 * sin0 * exp( 2.0 * b0 ) 00111 + a0 * cos0 * exp( 2.0 * b0 ) - 2.0 * a0 * exp( b0 )) 00112 * exp( b0 ) 00113 / (exp( 4.0 * b0 ) - 4.0 * cos0 * exp( 3.0 * b0 ) 00114 + 2.0 * exp( 2.0 * b0 ) + 4.0 * cos0 * cos0 * exp( 2.0 * b0 ) 00115 + 1.0 - 4.0 * cos0 * exp( b0 )); 00116 sumC = - 2.f * 00117 (c0 * cos1 - c1 * sin1 + c1 * sin1 * exp( 2.0 * b1 ) 00118 + c0 * cos1 * exp( 2.0 * b1 ) - 2.0 * c0 * exp( b1 )) 00119 * exp( b1 ) / 00120 (exp( 4.0 * b1 ) - 4.0 * cos1 * exp( 3.0 * b1 ) 00121 + 2.0 * exp( 2.0 * b1 ) + 4.0 * cos1 * cos1 * exp( 2.0 * b1 ) 00122 + 1.0 - 4.0 * cos1 * exp( b1 )); 00123 break; 00124 } 00125 00126 case DericheGaussianSecondDerivative : 00127 { 00128 double aux; 00129 aux = 00130 12.0 * cos0 * exp( 3.0 * b0 ) - 3.0 * exp( 2.0 * b0 ) 00131 + 8.0 * cos0 * cos0 * cos0 * exp( 3.0 * b0 ) - 12.0 * cos0 * cos0 * 00132 exp( 4.0 * b0 ) 00133 - (3.0 * exp( 4.0 * b0 )) 00134 + 6.0 * cos0 * exp( 5.0 * b0 ) - exp( 6.0 * b0 ) + 6.0 * cos0 * exp 00135 ( b0 ) 00136 - ( 1.0 + 12.0 * cos0 * cos0 * exp( 2.0 * b0 ) ); 00137 sumA = 00138 4.0 * a0 * sin0 * exp( 3.0 * b0 ) + a1 * cos0 * cos0 * exp( 4.0 * b0 ) 00139 - ( 4.0 * a0 * sin0 * exp( b0 ) + 6.0 * a1 * cos0 * cos0 * exp( 2.0 * b0 ) ) 00140 + 2.0 * a1 * cos0 * cos0 * cos0 * exp( b0 ) - 2.0 * a1 * cos0 * exp(b0) 00141 + 2.0 * a1 * cos0 * cos0 * cos0 * exp( 3.0 * b0 ) - 2.0 * a1 * cos0 00142 * exp( 3.0 * b0 ) 00143 + a1 * cos0 * cos0 - a1 * exp( 4.0 * b0 ) 00144 + 2.0 * a0 * sin0 * cos0 * cos0 * exp( b0 ) - 2.0 * a0 * sin0 * cos0 00145 * cos0 * exp( 3.0 * b0 ) 00146 - ( a0 * sin0 * cos0 * exp( 4.0 * b0 ) + a1 ) 00147 + 6.0 * a1 * exp( 2.0 * b0 ) + a0 * cos0 * sin0 00148 * 2.0 * exp( b0 ) / ( aux * sin0 ); 00149 aux = 00150 12.0 * cos1 * exp( 3.0 * b1 ) - 3.0 * exp( 2.0 * b1 ) 00151 + 8.0 * cos1 * cos1 * cos1 * exp( 3.0 * b1 ) - 12.0 * cos1 * cos1 * 00152 exp( 4.0 * b1 ) 00153 - 3.0 * exp( 4.0 * b1 ) 00154 + 6.0 * cos1 * exp( 5.0 * b1 ) - exp( 6.0 * b1 ) + 6.0 * cos1 * exp 00155 ( b1 ) 00156 - ( 1.0 + 12.0 * cos1 * cos1 * exp( 2.0 * b1 ) ); 00157 sumC = 4.0 * c0 * sin1 * exp( 3.0 * b1 ) + c1 * cos1 * cos1 * exp( 4.0 * b1 ) 00158 - ( 4.0 * c0 * sin1 * exp( b1 ) + 6.0 * c1 * cos1 * cos1 * exp( 2.0 * b1 ) ) 00159 + 2.0 * c1 * cos1 * cos1 * cos1 * exp( b1 ) - 2.0 * c1 * cos1 * exp( b1 ) 00160 + 2.0 * c1 * cos1 * cos1 * cos1 * exp( 3.0 * b1 ) - 2.0 * c1 * cos1 00161 * exp( 3.0 * b1 ) 00162 + c1 * cos1 * cos1 - c1 * exp( 4.0 * b1 ) 00163 + 2.0 * c0 * sin1 * cos1 * cos1 * exp( b1 ) - 2.0 * c0 * sin1 * cos1 00164 * cos1 * exp( 3.0 * b1 ) 00165 - ( c0 * sin1 * cos1 * exp( 4.0 * b1 ) + c1 ) 00166 + 6.0 * c1 * exp( 2.0 * b1 ) + c0 * cos1 * sin1 00167 * 2.0 * exp( b1 ) / ( aux * sin1 ); 00168 sumA /= 2; 00169 sumC /= 2; 00170 break; 00171 } 00172 } 00173 00174 a0 /= (sumA + sumC); 00175 a1 /= (sumA + sumC); 00176 c0 /= (sumA + sumC); 00177 c1 /= (sumA + sumC); 00178 00179 n[3] = 00180 exp( -b1 - 2*b0 ) * (c1 * sin1 - cos1 * c0) + 00181 exp( -b0 - 2*b1 ) * (a1 * sin0 - cos0 * a0); 00182 n[2] = 00183 2 * exp(-b0 - b1) * ((a0 + c0) * cos1 * cos0 - 00184 cos1 * a1 * sin0 - 00185 cos0 * c1 * sin1) + 00186 c0 * exp(-2 * b0) + a0 * exp(-2 * b1); 00187 n[1] = 00188 exp(-b1) * (c1 * sin1 - (c0 + 2*a0) * cos1) + 00189 exp(-b0) * (a1 * sin0 - (2*c0 + a0) * cos0); 00190 n[0] = 00191 a0 + c0; 00192 00193 d[4] = exp(-2 * b0 - 2 * b1); 00194 d[3] = 00195 -2 * cos0 * exp(-b0 - 2*b1) - 00196 2 * cos1 * exp(-b1 - 2*b0); 00197 d[2] = 00198 4 * cos1 * cos0 * exp(-b0 - b1) + 00199 exp(-2*b1) + exp(-2*b0); 00200 d[1] = 00201 -2*exp(-b1) * cos1 - 2 * exp(-b0) * cos0; 00202 00203 switch (filter_type) { 00204 case DericheGaussian : 00205 case DericheGaussianSecondDerivative : 00206 { 00207 for (unsigned i = 1; i <= 3; ++i) 00208 { 00209 dm[i] = d[i]; 00210 nm[i] = n[i] - d[i] * n[0]; 00211 } 00212 dm[4] = d[4]; 00213 nm[4] = -d[4] * n[0]; 00214 break; 00215 } 00216 case DericheGaussianFirstDerivative : 00217 { 00218 for (unsigned i = 1; i <= 3; ++i) 00219 { 00220 dm[i] = d[i]; 00221 nm[i] = - (n[i] - d[i] * n[0]); 00222 } 00223 dm[4] = d[4]; 00224 nm[4] = d[4] * n[0]; 00225 break; 00226 } 00227 } 00228 } 00229 }; 00230 00231 00232 } // end of namespace mln::linear::my 00233 00234 00235 00236 00237 // FIXME: in the "generic" code below there is no test that we 00238 // actually stay in the image domain... 00239 00240 00241 template <typename I, typename C> 00242 inline 00243 void 00244 recursivefilter_directional_generic(I& ima, 00245 const C& c, 00246 const mln_psite(I)& start, 00247 const mln_psite(I)& finish, 00248 int len, 00249 const mln_deduce(I, psite, delta)& d) 00250 { 00251 std::vector<double> tmp1(len); 00252 std::vector<double> tmp2(len); 00253 00254 tmp1[0] = 00255 c.n[0]*ima(start); 00256 00257 tmp1[1] = 00258 c.n[0]*ima(start + d) 00259 + c.n[1]*ima(start) 00260 - c.d[1]*tmp1[0]; 00261 00262 tmp1[2] = 00263 c.n[0]*ima(start + d + d) 00264 + c.n[1]*ima(start + d) 00265 + c.n[2]*ima(start) 00266 - c.d[1]*tmp1[1] 00267 - c.d[2]*tmp1[0]; 00268 00269 tmp1[3] = 00270 c.n[0]*ima(start + d + d + d) 00271 + c.n[1]*ima(start + d + d) 00272 + c.n[2]*ima(start + d) 00273 + c.n[3]*ima(start) 00274 - c.d[1]*tmp1[2] - c.d[2]*tmp1[1] 00275 - c.d[3]*tmp1[0]; 00276 00277 mln_site(I) current = start + d + d + d + d; 00278 for (int i = 4; i < len; ++i) 00279 { 00280 tmp1[i] = 00281 c.n[0]*ima(current) 00282 + c.n[1]*ima(current - d) 00283 + c.n[2]*ima(current - d - d) 00284 + c.n[3]*ima(current - d - d - d) 00285 - c.d[1]*tmp1[i - 1] - c.d[2]*tmp1[i - 2] 00286 - c.d[3]*tmp1[i - 3] - c.d[4]*tmp1[i - 4]; 00287 current += d; 00288 } 00289 00290 // Non causal part 00291 00292 tmp2[len - 1] = 0; 00293 00294 tmp2[len - 2] = 00295 c.nm[1]*ima(finish); 00296 00297 tmp2[len - 3] = 00298 c.nm[1]*ima(finish - d) 00299 + c.nm[2]*ima(finish) 00300 - c.dm[1]*tmp2[len-2]; 00301 00302 tmp2[len - 4] = 00303 c.nm[1]*ima(finish - d - d) 00304 + c.nm[2]*ima(finish - d) 00305 + c.nm[3]*ima(finish) 00306 - c.dm[1]*tmp2[len-3] 00307 - c.dm[2]*tmp2[len-2]; 00308 00309 current = finish - d - d - d ; 00310 00311 for (int i = len - 5; i >= 0; --i) 00312 { 00313 tmp2[i] = 00314 c.nm[1]*ima(current) 00315 + c.nm[2]*ima(current + d) 00316 + c.nm[3]*ima(current + d + d) 00317 + c.nm[4]*ima(current + d + d + d) 00318 - c.dm[1]*tmp2[i+1] - c.dm[2]*tmp2[i+2] 00319 - c.dm[3]*tmp2[i+3] - c.dm[4]*tmp2[i+4]; 00320 current -= d; 00321 } 00322 00323 // Combine results from causal and non-causal parts. 00324 00325 current = start; 00326 for (int i = 0; i < len; ++i) 00327 { 00328 ima(current) = (tmp1[i] + tmp2[i]); 00329 current += d; 00330 } 00331 } 00332 00333 00334 00335 00336 template <typename I, typename C> 00337 inline 00338 void 00339 recursivefilter_directional_fastest(I& ima, 00340 const C& c, 00341 const mln_psite(I)& start, 00342 const mln_psite(I)& finish, 00343 int len, 00344 const mln_deduce(I, psite, delta)& d, 00345 const mln_value(I)& bdr) 00346 { 00347 // extension::adjust_fill(ima, 5 * int(151 + .50001) + 1, bdr); 00348 // extension::fill(ima, bdr); 00349 00350 std::vector<double> tmp1(len); 00351 std::vector<double> tmp2(len); 00352 00353 unsigned delta_offset = ima.delta_index(d); 00354 unsigned 00355 o_start = ima.index_of_point(start), 00356 o_start_d = o_start + delta_offset, 00357 o_start_dd = o_start + 2 * delta_offset, 00358 o_start_ddd = o_start + 3 * delta_offset; 00359 00360 tmp1[0] = 00361 c.n[0] * ima.element(o_start); 00362 00363 tmp1[1] = 0 00364 + c.n[0] * ima.element(o_start_d) 00365 + c.n[1] * ima.element(o_start) 00366 - c.d[1] * tmp1[0]; 00367 00368 tmp1[2] = 0 00369 + c.n[0] * ima.element(o_start_dd) 00370 + c.n[1] * ima.element(o_start_d) 00371 + c.n[2] * ima.element(o_start) 00372 - c.d[1] * tmp1[1] 00373 - c.d[2] * tmp1[0]; 00374 00375 tmp1[3] = 0 00376 + c.n[0] * ima.element(o_start_ddd) 00377 + c.n[1] * ima.element(o_start_dd) 00378 + c.n[2] * ima.element(o_start_d) 00379 + c.n[3] * ima.element(o_start) 00380 - c.d[1] * tmp1[2] - c.d[2] * tmp1[1] 00381 - c.d[3] * tmp1[0]; 00382 00383 unsigned 00384 o_current = o_start + 4 * delta_offset, 00385 o_current_d = o_current - delta_offset, 00386 o_current_dd = o_current - 2 * delta_offset, 00387 o_current_ddd = o_current - 3 * delta_offset; 00388 00389 for (int i = 4; i < len; ++i) 00390 { 00391 tmp1[i] = 0 00392 + c.n[0] * ima.element(o_current) 00393 + c.n[1] * ima.element(o_current_d) 00394 + c.n[2] * ima.element(o_current_dd) 00395 + c.n[3] * ima.element(o_current_ddd) 00396 - c.d[1] * tmp1[i - 1] - c.d[2] * tmp1[i - 2] 00397 - c.d[3] * tmp1[i - 3] - c.d[4] * tmp1[i - 4]; 00398 o_current += delta_offset; 00399 o_current_d += delta_offset; 00400 o_current_dd += delta_offset; 00401 o_current_ddd += delta_offset; 00402 } 00403 00404 // Non causal part 00405 00406 (void) bdr; 00407 // extension::fill(ima, bdr); 00408 00409 unsigned 00410 o_finish = ima.index_of_point(finish), 00411 o_finish_d = o_finish - delta_offset, 00412 o_finish_dd = o_finish - 2 * delta_offset; 00413 00414 tmp2[len - 1] = 0; 00415 00416 tmp2[len - 2] = 00417 c.nm[1] * ima.element(o_finish); 00418 00419 tmp2[len - 3] = 0 00420 + c.nm[1] * ima.element(o_finish_d) 00421 + c.nm[2] * ima.element(o_finish) 00422 - c.dm[1] * tmp2[len-2]; 00423 00424 tmp2[len - 4] = 0 00425 + c.nm[1] * ima.element(o_finish_dd) 00426 + c.nm[2] * ima.element(o_finish_d) 00427 + c.nm[3] * ima.element(o_finish) 00428 - c.dm[1] * tmp2[len-3] 00429 - c.dm[2] * tmp2[len-2]; 00430 00431 o_current = o_finish - 3 * delta_offset; 00432 o_current_d = o_current + delta_offset; 00433 o_current_dd = o_current + 2 * delta_offset; 00434 o_current_ddd = o_current + 3 * delta_offset; 00435 00436 for (int i = len - 5; i >= 0; --i) 00437 { 00438 tmp2[i] = 0 00439 + c.nm[1] * ima.element(o_current) 00440 + c.nm[2] * ima.element(o_current_d) 00441 + c.nm[3] * ima.element(o_current_dd) 00442 + c.nm[4] * ima.element(o_current_ddd) 00443 - c.dm[1] * tmp2[i+1] - c.dm[2] * tmp2[i+2] 00444 - c.dm[3] * tmp2[i+3] - c.dm[4] * tmp2[i+4]; 00445 o_current -= delta_offset; 00446 o_current_d -= delta_offset; 00447 o_current_dd -= delta_offset; 00448 o_current_ddd -= delta_offset; 00449 } 00450 00451 // Combine results from causal and non-causal parts. 00452 00453 o_current = o_start; 00454 for (int i = 0; i < len; ++i) 00455 { 00456 ima.element(o_current) = static_cast<mln_value(I)>(tmp1[i] + tmp2[i]); 00457 o_current += delta_offset; 00458 } 00459 } 00460 00461 00462 template <typename I> 00463 inline 00464 mln_concrete(I) 00465 gaussian_directional_2d(const Image<I>& input_, 00466 unsigned dir, double sigma, 00467 const mln_value(I)& bdr) 00468 { 00469 trace::entering("linear::gaussian_directional_2d"); 00470 00471 typedef mln_site(I) P; 00472 mlc_bool(P::dim == 2)::check(); 00473 00474 const I& input = exact(input_); 00475 00476 mln_precondition(dir == 0 || dir == 1); 00477 mln_precondition(input.is_valid()); 00478 00479 my::recursivefilter_coef_ coef(1.68f, 3.735f, 00480 1.783f, 1.723f, 00481 -0.6803f, -0.2598f, 00482 0.6318f, 1.997f, 00483 sigma, 00484 my::recursivefilter_coef_::DericheGaussian); 00485 00486 extension::adjust_fill(input, 5 * int(sigma + .50001) + 1, bdr); 00487 mln_concrete(I) output = duplicate(input); 00488 00489 if (sigma < 0.006) 00490 return output; 00491 00492 int 00493 nrows = geom::nrows(input), 00494 ncols = geom::ncols(input), 00495 b = input.border(); 00496 00497 if (dir == 0) 00498 { 00499 for (int j = 0; j < ncols; ++j) 00500 recursivefilter_directional_fastest(output, coef, 00501 point2d(- b, j), 00502 point2d(nrows - 1 + b, j), 00503 nrows + 2 * b, 00504 dpoint2d(1, 0), 00505 bdr); 00506 } 00507 00508 if (dir == 1) 00509 { 00510 for (int i = 0; i < nrows; ++i) 00511 recursivefilter_directional_fastest(output, coef, 00512 point2d(i, - b), 00513 point2d(i, ncols - 1 + b), 00514 ncols + 2 * b, 00515 dpoint2d(0, 1), 00516 bdr); 00517 } 00518 00519 trace::exiting("linear::gaussian_directional_2d"); 00520 return output; 00521 } 00522 00523 # endif // ! MLN_INCLUDE_ONLY 00524 00525 } // end of namespace mln::linear 00526 00527 } // end of namespace mln 00528 00529 00530 #endif // ! MLN_LINEAR_GAUSSIAN_DIRECTIONAL_2D_HH