[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/tv_filter.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2012 by Frank Lenzen & */ 00004 /* Ullrich Koethe */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* The VIGRA Website is */ 00008 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00009 /* Please direct questions, bug reports, and contributions to */ 00010 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00011 /* vigra@informatik.uni-hamburg.de */ 00012 /* */ 00013 /* Permission is hereby granted, free of charge, to any person */ 00014 /* obtaining a copy of this software and associated documentation */ 00015 /* files (the "Software"), to deal in the Software without */ 00016 /* restriction, including without limitation the rights to use, */ 00017 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00018 /* sell copies of the Software, and to permit persons to whom the */ 00019 /* Software is furnished to do so, subject to the following */ 00020 /* conditions: */ 00021 /* */ 00022 /* The above copyright notice and this permission notice shall be */ 00023 /* included in all copies or substantial portions of the */ 00024 /* Software. */ 00025 /* */ 00026 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00027 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00028 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00029 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00030 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00031 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00032 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00033 /* OTHER DEALINGS IN THE SOFTWARE. */ 00034 /* */ 00035 /************************************************************************/ 00036 00037 #ifndef VIGRA_TV_FILTER_HXX 00038 #define VIGRA_TV_FILTER_HXX 00039 00040 #include <iostream> 00041 #include <cmath> 00042 #include "config.hxx" 00043 #include "impex.hxx" 00044 #include "separableconvolution.hxx" 00045 #include "multi_array.hxx" 00046 #include "multi_math.hxx" 00047 #include "eigensystem.hxx" 00048 #include "convolution.hxx" 00049 #include "fixedpoint.hxx" 00050 #include "project2ellipse.hxx" 00051 00052 #ifndef VIGRA_MIXED_2ND_DERIVATIVES 00053 #define VIGRA_MIXED_2ND_DERIVATIVES 1 00054 #endif 00055 00056 00057 namespace vigra{ 00058 00059 00060 00061 /** \addtogroup NonLinearDiffusion 00062 */ 00063 00064 //@{ 00065 00066 /********************************************************/ 00067 /* */ 00068 /* totalVariationFilter */ 00069 /* */ 00070 /********************************************************/ 00071 00072 /** \brief Performs standard Total Variation Regularization 00073 00074 The algorithm minimizes 00075 00076 \f[ 00077 \min_u \int_\Omega \frac{1}{2} (u-f)^2\;dx + \alpha TV(u)\qquad\qquad (1) 00078 \f] 00079 where <em>\f$ f=f(x)\f$</em> are the two dimensional noisy data, 00080 <em> \f$ u=u(x)\f$</em> are the smoothed data,<em>\f$ \alpha \ge 0 \f$</em> 00081 is the filter parameter and <em>\f$ TV(u)\f$ </em> is the total variation semi-norm. 00082 00083 <b> Declarations:</b> 00084 00085 \code 00086 namespace vigra { 00087 template <class stride1,class stride2> 00088 void totalVariationFilter(MultiArrayView<2,double,stride1> data, 00089 MultiArrayView<2,double,stride2> out, 00090 double alpha, 00091 int steps, 00092 double eps=0); 00093 void totalVariationFilter(MultiArrayView<2,double,stride1> data, 00094 MultiArrayView<2,double,stride2> weight, 00095 MultiArrayView<2,double,stride3> out, 00096 double alpha, 00097 int steps, 00098 double eps=0); 00099 } 00100 \endcode 00101 00102 \ref totalVariationFilter() implements a primal-dual algorithm to solve (1). 00103 00104 Input: 00105 <table> 00106 <tr><td><em>data</em>: </td><td> input data to be smoothed. </td></tr> 00107 <tr><td><em>alpha</em>: </td><td> smoothing parameter.</td></tr> 00108 <tr><td><em>steps</em>: </td><td> maximal number of iteration steps. </td></tr> 00109 <tr><td><em>eps</em>: </td><td> The algorithm stops, if the primal-dual gap is below the threshold <em>eps</em>. 00110 </table> 00111 00112 Output: 00113 00114 <em>out</em> contains the filtered data. 00115 00116 In addition, a point-wise weight (\f$ \ge 0 \f$) for the data term can be provided (overloaded function). 00117 00118 <b> Usage:</b> 00119 00120 <b>\#include</b> <vigra/tv_filter.hxx> 00121 00122 \code 00123 MultiArray<2,double> data(Shape2(width,height)); // to be initialized 00124 MultiArray<2,double> out(Shape2(width,height)); 00125 MultiArray<2,double> weight(Shape2(width,height))); //optional argument in overloaded function, to be initialized if used 00126 int steps; // to be initialized 00127 double alpha,eps; // to be initialized 00128 00129 totalVariationFilter(data,out,alpha,steps,eps); 00130 \endcode 00131 or 00132 \code 00133 totalVariationFilter(data,weight,out,alpha,steps,eps); 00134 \endcode 00135 00136 */ 00137 doxygen_overloaded_function(template <...> void totalVariationFilter) 00138 00139 template <class stride1,class stride2> 00140 void totalVariationFilter(MultiArrayView<2,double,stride1> data,MultiArrayView<2,double,stride2> out, double alpha, int steps, double eps=0){ 00141 00142 using namespace multi_math; 00143 00144 MultiArray<2,double> temp1(data.shape()),temp2(data.shape()),vx(data.shape()),vy(data.shape()),u_bar(data.shape()); 00145 Kernel1D<double> Lx,LTx; 00146 Lx.initExplicitly(0,1)=1,-1; // = Right sided finite differences for d/dx and d/dy 00147 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT); // with hom. Neumann boundary conditions 00148 LTx.initExplicitly(-1,0)=-1,1; // = Left sided finite differences for -d/dx and -d/dy 00149 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD); // with hom. Dirichlet b. c. 00150 00151 out=data; 00152 u_bar=data; 00153 00154 double tau=1.0 / std::max(alpha,1.) / std::sqrt(8.0) * 0.06; 00155 double sigma=1.0 / std::sqrt(8.0) / 0.06; 00156 00157 for (int i=0;i<steps;i++){ 00158 00159 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx)); 00160 vx+=(sigma*temp1); 00161 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx)); 00162 vy+=(sigma*temp1); 00163 00164 //project to constraint set 00165 for (int y=0;y<data.shape(1);y++){ 00166 for (int x=0;x<data.shape(0);x++){ 00167 double l=hypot(vx(x,y),vy(x,y)); 00168 if (l>1){ 00169 vx(x,y)/=l; 00170 vy(x,y)/=l; 00171 } 00172 } 00173 } 00174 00175 separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx)); 00176 separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx)); 00177 u_bar=out; 00178 out-=tau*(out-data+alpha*(temp1+temp2)); 00179 u_bar=2*out-u_bar; //cf. Chambolle/Pock and Popov's algorithm 00180 00181 00182 //stopping criterion 00183 if (eps>0){ 00184 separableConvolveX(srcImageRange(out),destImage(temp1),kernel1d(Lx)); 00185 separableConvolveY(srcImageRange(out),destImage(temp2),kernel1d(Lx)); 00186 00187 double f_primal=0,f_dual=0; 00188 for (int y=0;y<data.shape(1);y++){ 00189 for (int x=0;x<data.shape(0);x++){ 00190 f_primal+=.5*(out(x,y)-data(x,y))*(out(x,y)-data(x,y))+alpha*hypot(temp1(x,y),temp2(x,y)); 00191 } 00192 } 00193 separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx)); 00194 separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx)); 00195 for (int y=0;y<data.shape(1);y++){ 00196 for (int x=0;x<data.shape(0);x++){ 00197 double divv=temp1(x,y)+temp2(x,y); 00198 f_dual+=-.5*alpha*alpha*(divv*divv)+alpha*data(x,y)*divv; 00199 } 00200 } 00201 if (f_primal>0 && (f_primal-f_dual)/f_primal<eps){ 00202 break; 00203 } 00204 } 00205 } 00206 } 00207 00208 template <class stride1,class stride2, class stride3> 00209 void totalVariationFilter(MultiArrayView<2,double,stride1> data,MultiArrayView<2,double,stride2> weight, MultiArrayView<2,double,stride3> out,double alpha, int steps, double eps=0){ 00210 00211 using namespace multi_math; 00212 00213 MultiArray<2,double> temp1(data.shape()),temp2(data.shape()),vx(data.shape()),vy(data.shape()),u_bar(data.shape()); 00214 Kernel1D<double> Lx,LTx; 00215 Lx.initExplicitly(0,1)=1,-1; // = Right sided finite differences for d/dx and d/dy 00216 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT); // with hom. Neumann boundary conditions 00217 LTx.initExplicitly(-1,0)=-1,1; // = Left sided finite differences for -d/dx and -d/dy 00218 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD); // with hom. Dirichlet b. c. 00219 00220 out=data; 00221 u_bar=data; 00222 00223 double tau=1.0 / std::max(alpha,1.) / std::sqrt(8.0) * 0.06; 00224 double sigma=1.0 / std::sqrt(8.0) / 0.06; 00225 00226 for (int i=0;i<steps;i++){ 00227 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx)); 00228 vx+=(sigma*temp1); 00229 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx)); 00230 vy+=(sigma*temp1); 00231 00232 //project to constraint set 00233 for (int y=0;y<data.shape(1);y++){ 00234 for (int x=0;x<data.shape(0);x++){ 00235 double l=hypot(vx(x,y),vy(x,y)); 00236 if (l>1){ 00237 vx(x,y)/=l; 00238 vy(x,y)/=l; 00239 } 00240 } 00241 } 00242 00243 separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx)); 00244 separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx)); 00245 u_bar=out; 00246 out-=tau*(weight*(out-data)+alpha*(temp1+temp2)); 00247 u_bar=2*out-u_bar; 00248 00249 00250 //stopping criterion 00251 if (eps>0){ 00252 separableConvolveX(srcImageRange(out),destImage(temp1),kernel1d(Lx)); 00253 separableConvolveY(srcImageRange(out),destImage(temp2),kernel1d(Lx)); 00254 00255 double f_primal=0,f_dual=0; 00256 for (int y=0;y<data.shape(1);y++){ 00257 for (int x=0;x<data.shape(0);x++){ 00258 f_primal+=.5*weight(x,y)*(out(x,y)-data(x,y))*(out(x,y)-data(x,y))+alpha*hypot(temp1(x,y),temp2(x,y)); 00259 } 00260 } 00261 separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx)); 00262 separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx)); 00263 for (int y=0;y<data.shape(1);y++){ 00264 for (int x=0;x<data.shape(0);x++){ 00265 double divv=temp1(x,y)+temp2(x,y); 00266 f_dual+=-.5*alpha*alpha*(weight(x,y)*divv*divv)+alpha*data(x,y)*divv; 00267 } 00268 } 00269 if (f_primal>0 && (f_primal-f_dual)/f_primal<eps){ 00270 break; 00271 } 00272 } 00273 } 00274 } 00275 //<!--\f$ \alpha(x)=\beta(x)=\beta_{par}\f$ in homogeneous regions without edges, 00276 //and \f$ \alpha(x)=\alpha_{par}\f$ at edges.--> 00277 00278 00279 /********************************************************/ 00280 /* */ 00281 /* getAnisotropy */ 00282 /* */ 00283 /********************************************************/ 00284 00285 /** \brief Sets up directional data for anisotropic regularization 00286 00287 This routine provides a two-dimensional normalized vector field \f$ v \f$, which is normal to edges in the given data, 00288 found as the eigenvector of the structure tensor belonging to the largest eigenvalue. 00289 \f$ v \f$ is encoded by a scalar field \f$ \varphi \f$ of angles, i.e. 00290 \f$ v(x)=(\cos(\varphi(x)),\sin(\varphi(x)))^\top \f$. 00291 00292 In addition, two scalar fields \f$ \alpha \f$ and \f$ \beta \f$ are generated from 00293 scalar parameters \f$ \alpha_{par}\f$ and \f$ \beta_{par}\f$, such that 00294 <center> 00295 <table> 00296 <tr> <td>\f$ \alpha(x)= \alpha_{par}\f$ at edges,</td></tr> 00297 <tr> <td>\f$ \alpha(x)= \beta_{par}\f$ in homogeneous regions,</td></tr> 00298 <tr> <td>\f$ \beta(x)=\beta_{par}\f$ .</td></tr> 00299 </table> 00300 </center> 00301 00302 <b> Declarations:</b> 00303 00304 \code 00305 namespace vigra { 00306 void getAnisotropy(MultiArrayView<2,double,stride1> data, 00307 MultiArrayView<2,double,stride2> phi, 00308 MultiArrayView<2,double,stride3> alpha, 00309 MultiArrayView<2,double,stride4> beta, 00310 double alpha_par, 00311 double beta_par, 00312 double sigma_par, 00313 double rho_par, 00314 double K_par); 00315 } 00316 \endcode 00317 00318 Output: 00319 <table> 00320 <tr><td>Three scalar fields <em>phi</em>, <em>alpha</em> and <em>beta</em>.</td></tr> 00321 </table> 00322 00323 Input: 00324 <table> 00325 <tr><td><em>data</em>:</td><td>two-dimensional scalar field.</td></tr> 00326 <tr><td><em>alpha_par,beta_par</em>:</td><td>two positive values for setting up the scalar fields alpha and beta</tr> 00327 <tr><td><em>sigma_par</em>:</td><td> non-negative parameter for presmoothing the data.</td></tr> 00328 <tr><td> <em>rho_par</em>:</td><td> non-negative parameter for presmoothing the structure tensor.</td></tr> 00329 <tr><td><em>K_par</em>:</td><td> positive edge sensitivity parameter.</td></tr> 00330 </table> 00331 00332 (see \ref anisotropicTotalVariationFilter() and \ref secondOrderTotalVariationFilter() for usage in an application). 00333 */ 00334 doxygen_overloaded_function(template <...> void getAnisotropy) 00335 00336 template <class stride1,class stride2,class stride3,class stride4> 00337 void getAnisotropy(MultiArrayView<2,double,stride1> data,MultiArrayView<2,double,stride2> phi, 00338 MultiArrayView<2,double,stride3> alpha, MultiArrayView<2,double,stride4> beta, 00339 double alpha_par, double beta_par, double sigma_par, double rho_par, double K_par){ 00340 00341 using namespace multi_math; 00342 00343 MultiArray<2,double> smooth(data.shape()),tmp(data.shape()); 00344 vigra::Kernel1D<double> gauss; 00345 00346 00347 gauss.initGaussian(sigma_par); 00348 separableConvolveX(srcImageRange(data), destImage(tmp), kernel1d(gauss)); 00349 separableConvolveY(srcImageRange(tmp), destImage(smooth), kernel1d(gauss)); 00350 00351 MultiArray<2,double> stxx(data.shape()),stxy(data.shape()),styy(data.shape()); 00352 00353 // calculate Structure Tensor at inner scale = sigma and outer scale = rho 00354 vigra::structureTensor(srcImageRange(smooth),destImage(stxx), destImage(stxy), destImage(styy),1.,1.); 00355 00356 gauss.initGaussian(rho_par); 00357 separableConvolveX(srcImageRange(stxx), destImage(tmp), kernel1d(gauss)); 00358 separableConvolveY(srcImageRange(tmp), destImage(stxx), kernel1d(gauss)); 00359 separableConvolveX(srcImageRange(stxy), destImage(tmp), kernel1d(gauss)); 00360 separableConvolveY(srcImageRange(tmp), destImage(stxy), kernel1d(gauss)); 00361 separableConvolveX(srcImageRange(styy), destImage(tmp), kernel1d(gauss)); 00362 separableConvolveY(srcImageRange(tmp), destImage(styy), kernel1d(gauss)); 00363 00364 MultiArray<2,double> matrix(Shape2(2,2)),ev(Shape2(2,2)),ew(Shape2(2,1)); 00365 00366 for (int y=0;y<data.shape(1);y++){ 00367 for (int x=0;x<data.shape(0);x++){ 00368 00369 matrix(0,0)=stxx(x,y); 00370 matrix(1,1)=styy(x,y); 00371 matrix(0,1)=stxy(x,y); 00372 matrix(1,0)=stxy(x,y); 00373 vigra::linalg::detail::symmetricEigensystem2x2(matrix,ew,ev); 00374 00375 phi(x,y)=std::atan2(ev(1,0),ev(0,0)); 00376 double coherence=ew(0,0)-ew(1,0); 00377 double c=std::min(K_par*coherence,1.); 00378 alpha(x,y)=alpha_par*c+(1-c)*beta_par; 00379 beta(x,y)=beta_par; 00380 } 00381 } 00382 } 00383 00384 /********************************************************/ 00385 /* */ 00386 /* anisotropicTotalVariationFilter */ 00387 /* */ 00388 /********************************************************/ 00389 00390 /** \brief Performs Anisotropic Total Variation Regularization 00391 00392 The algorithm minimizes 00393 \f[ 00394 \min_u \int_\Omega \frac{1}{2} (u-f)^2 + \sqrt{\nabla u^\top A \nabla u}\;dx\qquad\qquad(2) 00395 \f] 00396 00397 where <em> \f$ f=f(x)\f$ </em> are the noisy data, <em> \f$ u=u(x)\f$ </em> are the smoothed data,<em>\f$ \nabla u \f$ </em> 00398 is the image gradient in the sense of Total Variation and <em>\f$ A \f$ </em> is a locally varying symmetric, positive definite 2x2 matrix. 00399 00400 Matrix <em>\f$ A \f$ </em> is described by providing for each data point a normalized eigenvector (via angle \f$ \phi \f$) 00401 and two eigenvalues \f$ \alpha>0 \f$ and \f$ \beta>0 \f$. 00402 00403 \ref getAnisotropy() can be use to set up such data \f$ \phi,\alpha,\beta \f$ by providing a vector field normal to edges. 00404 00405 <b> Declarations:</b> 00406 00407 \code 00408 namespace vigra { 00409 template <class stride1,class stride2,class stride3,class stride4,class stride5,class stride6> 00410 void anisotropicTotalVariationFilter(MultiArrayView<2,double,stride1> data, 00411 MultiArrayView<2,double,stride2> weight, 00412 MultiArrayView<2,double,stride3> phi, 00413 MultiArrayView<2,double,stride4> alpha, 00414 MultiArrayView<2,double,stride5> beta, 00415 MultiArrayView<2,double,stride6> out, 00416 int steps); 00417 } 00418 \endcode 00419 00420 \ref anisotropicTotalVariationFilter() implements a primal-dual algorithm to solve (2). 00421 00422 Input: 00423 <table> 00424 <tr><td><em>data</em>:</td><td>input data to be filtered. </td></tr> 00425 <tr><td><em>steps</em>:</td><td>iteration steps.</td></tr> 00426 <tr><td><em>weight</em> :</td><td>a point-wise weight (\f$ \ge 0 \f$ ) for the data term.</td></tr> 00427 <tr><td><em>phi</em>,<em>alpha</em> and <em>beta</em> :</td><td> describe matrix \f$ A \f$, see above.</td></tr> 00428 </table> 00429 00430 Output: 00431 <table> 00432 <tr><td><em>out</em> :</td><td>contains filtered data.</td></tr> 00433 </table> 00434 00435 <b> Usage:</b> 00436 00437 E.g. with a solution-dependent adaptivity cf. [1], by updating the matrix \f$ A=A(u)\f$ 00438 in an outer loop: 00439 00440 <b>\#include</b> <vigra/tv_filter.hxx> 00441 00442 \code 00443 MultiArray<2,double> data(Shape2(width,height)); //to be initialized 00444 MultiArray<2,double> out (Shape2(width,height)); 00445 MultiArray<2,double> weight(Shape2(width,height)); //to be initialized 00446 MultiArray<2,double> phi (Shape2(width,height)); 00447 MultiArray<2,double> alpha(Shape2(width,height)); 00448 MultiArray<2,double> beta (Shape2(width,height)); 00449 double alpha0,beta0,sigma,rho,K; //to be initialized 00450 int outer_steps,inner_steps;//to be initialized 00451 00452 out=data; // data serves as initial value 00453 00454 for (int i=0;i<outer_steps;i++){ 00455 getAnisotropy(out,phi,alpha,beta,alpha0,beta0,sigma,rho,K); // sets phi, alpha, beta 00456 anisotropicTotalVariationFilter(data,weight,phi,alpha,beta,out,inner_steps); 00457 } 00458 \endcode 00459 00460 [1] Frank Lenzen, Florian Becker, Jan Lellmann, Stefania Petra and Christoph Schnörr, A Class of Quasi-Variational Inequalities for Adaptive Image Denoising and Decomposition, Computational Optimization and Applications, Springer, 2012. 00461 */ 00462 doxygen_overloaded_function(template <...> void anisotropicTotalVariationFilter) 00463 00464 template <class stride1,class stride2,class stride3,class stride4,class stride5,class stride6> 00465 void anisotropicTotalVariationFilter(MultiArrayView<2,double,stride1> data,MultiArrayView<2,double,stride2> weight, 00466 MultiArrayView<2,double,stride3> phi,MultiArrayView<2,double,stride4> alpha, 00467 MultiArrayView<2,double,stride5> beta,MultiArrayView<2,double,stride6> out, 00468 int steps){ 00469 00470 using namespace multi_math; 00471 00472 MultiArray<2,double> temp1(data.shape()),temp2(data.shape()),vx(data.shape()),vy(data.shape()),u_bar(data.shape()); 00473 MultiArray<2,double> rx(data.shape()),ry(data.shape()); 00474 00475 Kernel1D<double> Lx,LTx; 00476 Lx.initExplicitly(0,1)=1,-1; // = Right sided finite differences for d/dx and d/dy 00477 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT); // with hom. Neumann boundary conditions 00478 LTx.initExplicitly(-1,0)=-1,1; // = Left sided finite differences for -d/dx and -d/dy 00479 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD); // with hom. Dirichlet b. c. 00480 00481 u_bar=out; 00482 00483 double m=0; 00484 for (int y=0;y<data.shape(1);y++){ 00485 for (int x=0;x<data.shape(0);x++){ 00486 m=std::max(m,alpha(x,y)); 00487 m=std::max(m,beta (x,y)); 00488 } 00489 } 00490 m=std::max(m,1.); 00491 double tau=.9/m/std::sqrt(8.)*0.06; 00492 double sigma=.9/m/std::sqrt(8.)/0.06; 00493 00494 00495 for (int i=0;i<steps;i++){ 00496 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx)); 00497 vx+=(sigma*temp1); 00498 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx)); 00499 vy+=(sigma*temp1); 00500 00501 //project to constraint set 00502 for (int y=0;y<data.shape(1);y++){ 00503 for (int x=0;x<data.shape(0);x++){ 00504 double e1,e2,skp1,skp2; 00505 00506 e1=std::cos(phi(x,y)); 00507 e2=std::sin(phi(x,y)); 00508 skp1=vx(x,y)*e1+vy(x,y)*e2; 00509 skp2=vx(x,y)*(-e2)+vy(x,y)*e1; 00510 vigra::detail::projectEllipse2D (skp1,skp2,alpha(x,y),beta(x,y),0.001,100); 00511 00512 vx(x,y)=skp1*e1-skp2*e2; 00513 vy(x,y)=skp1*e2+skp2*e1; 00514 } 00515 } 00516 00517 separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx)); 00518 separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx)); 00519 u_bar=out; 00520 out-=tau*(weight*(out-data)+(temp1+temp2)); 00521 u_bar=2*out-u_bar; //cf. Chambolle/Pock and Popov's algorithm 00522 } 00523 } 00524 00525 /********************************************************/ 00526 /* */ 00527 /* secondOrderTotalVariationFilter */ 00528 /* */ 00529 /********************************************************/ 00530 00531 /** \brief Performs Anisotropic Total Variation Regularization 00532 00533 The algorithm minimizes 00534 00535 \f[ 00536 \min_u \int_\Omega \frac{1}{2} (u-f)^2 + \sqrt{\nabla u^\top A \nabla u} + \gamma |Hu|_F\;dx \qquad\qquad (3) 00537 \f] 00538 where <em> \f$ f=f(x)\f$ </em> are the noisy data, <em> \f$ u=u(x)\f$ </em> are the smoothed data,<em>\f$ \nabla u \f$ </em> 00539 is the image gradient in the sense of Total Variation, <em>\f$ A \f$ </em> is a locally varying 00540 symmetric, positive-definite 2x2 matrix and <em>\f$ |Hu|_F \f$</em> is the Frobenius norm of the Hessian of \f$ u \f$. 00541 00542 Matrix <em>\f$ A \f$ </em> is described by providing for each data point a normalized eigenvector (via angle \f$ \phi \f$) 00543 and two eigenvalues \f$ \alpha>0 \f$ and \f$ \beta>0 \f$. 00544 \ref getAnisotropy() can be use to set up such data \f$ \phi,\alpha, \beta \f$ by providing a vector field normal to edges. 00545 00546 \f$ \gamma>0 \f$ is the locally varying regularization parameter for second order. 00547 00548 <b> Declarations:</b> 00549 00550 \code 00551 namespace vigra { 00552 template <class stride1,class stride2,...,class stride9> 00553 void secondOrderTotalVariationFilter(MultiArrayView<2,double,stride1> data, 00554 MultiArrayView<2,double,stride2> weight, 00555 MultiArrayView<2,double,stride3> phi, 00556 MultiArrayView<2,double,stride4> alpha, 00557 MultiArrayView<2,double,stride5> beta, 00558 MultiArrayView<2,double,stride6> gamma, 00559 MultiArrayView<2,double,stride7> xedges, 00560 MultiArrayView<2,double,stride8> yedges, 00561 MultiArrayView<2,double,stride9> out, 00562 int steps); 00563 } 00564 \endcode 00565 00566 \ref secondOrderTotalVariationFilter() implements a primal-dual algorithm to solve (3). 00567 00568 Input: 00569 <table> 00570 <tr><td><em>data</em>: </td><td> the input data to be filtered. </td></tr> 00571 <tr><td><em>steps</em> : </td><td> number of iteration steps.</td></tr> 00572 <tr><td><em>out</em> : </td><td>contains the filtered data.</td></tr> 00573 <tr><td><em>weight</em> : </td><td> point-wise weight (\f$ \ge 0\f$ ) for the data term.</td></tr> 00574 <tr><td><em>phi</em>,<em>alpha</em>,<em>beta</em>: </td><td> describe matrix \f$ A\f$, see above.</td></tr> 00575 <tr><td><em> xedges </em> and <em> yedges </em>: </td><td> binary arrays indicating the 00576 presence of horizontal (between (x,y) and (x+1,y)) and vertical edges (between (x,y) and (x,y+1)). 00577 These data are considered in the calculation of \f$ Hu\f$, such that 00578 finite differences across edges are artificially set to zero to avoid second order smoothing over edges.</td></tr> 00579 </table> 00580 00581 <b> Usage:</b> 00582 00583 E.g. with a solution-dependent adaptivity (cf.[1]), by updating the matrix \f$ A=A(u)\f$ 00584 in an outer loop: 00585 00586 <b>\#include</b> <vigra/tv_filter.hxx> 00587 00588 \code 00589 MultiArray<2,double> data(Shape2(width,height)); //to be initialized 00590 MultiArray<2,double> out(Shape2(width,height)); 00591 MultiArray<2,double> weight(Shape2(width,height))); //to be initialized 00592 MultiArray<2,double> phi(Shape2(width,height); 00593 MultiArray<2,double> alpha(Shape2(width,height); 00594 MultiArray<2,double> beta(Shape2(width,height)); 00595 MultiArray<2,double> gamma(Shape2(width,height)); //to be initialized 00596 MultiArray<2,double> xedges(Shape2(width,height)); //to be initialized 00597 MultiArray<2,double> yedges(Shape2(width,height)); //to be initialized 00598 00599 00600 double alpha0,beta0,sigma,rho,K; //to be initialized 00601 int outer_steps,inner_steps;//to be initialized 00602 00603 out=data; // data serves as initial value 00604 00605 for (int i=0;i<outer_steps;i++){ 00606 00607 getAnisotropy(out,phi,alpha,beta,alpha0,beta0,sigma,rho,K); // sets phi, alpha, beta 00608 secondOrderTotalVariationFilter(data,weight,phi,alpha,beta,gamma,xedges,yedges,out,inner_steps); 00609 } 00610 \endcode 00611 00612 00613 [1] Frank Lenzen, Florian Becker, Jan Lellmann, Stefania Petra and Christoph Schnörr, A Class of Quasi-Variational Inequalities for Adaptive Image Denoising and Decomposition, Computational Optimization and Applications, Springer, 2012. 00614 */ 00615 doxygen_overloaded_function(template <...> void secondOrderTotalVariationFilter) 00616 00617 template <class stride1,class stride2,class stride3,class stride4,class stride5,class stride6,class stride7,class stride8,class stride9> 00618 void secondOrderTotalVariationFilter(MultiArrayView<2,double,stride1> data, 00619 MultiArrayView<2,double,stride2> weight,MultiArrayView<2,double,stride3> phi, 00620 MultiArrayView<2,double,stride4> alpha,MultiArrayView<2,double,stride5> beta, 00621 MultiArrayView<2,double,stride6> gamma, 00622 MultiArrayView<2,double,stride7> xedges,MultiArrayView<2,double,stride8> yedges, 00623 MultiArrayView<2,double,stride9> out, 00624 int steps){ 00625 00626 using namespace multi_math; 00627 00628 MultiArray<2,double> temp1(data.shape()),temp2(data.shape()),vx(data.shape()),vy(data.shape()),u_bar(data.shape()); 00629 MultiArray<2,double> rx(data.shape()),ry(data.shape()); 00630 MultiArray<2,double> wx(data.shape()),wy(data.shape()),wz(data.shape()); 00631 00632 00633 Kernel1D<double> Lx,LTx; 00634 Lx.initExplicitly(0,1)=1,-1; // = Right sided finite differences for d/dx and d/dy 00635 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT); // with hom. Neumann boundary conditions 00636 LTx.initExplicitly(-1,0)=-1,1; // = Left sided finite differences for -d/dx and -d/dy 00637 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD); // with hom. Dirichlet b. c. 00638 00639 u_bar=out; 00640 00641 double m=0; 00642 for (int y=0;y<data.shape(1);y++){ 00643 for (int x=0;x<data.shape(0);x++){ 00644 m=std::max(m,alpha(x,y)); 00645 m=std::max(m,beta (x,y)); 00646 } 00647 } 00648 m=std::max(m,1.); 00649 double tau=.1/m;//std::sqrt(8)*0.06; 00650 double sigma=.1;//m;/std::sqrt(8)/0.06; 00651 00652 for (int i=0;i<steps;i++){ 00653 00654 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx)); 00655 vx+=(sigma*temp1); 00656 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx)); 00657 vy+=(sigma*temp1); 00658 00659 00660 // update wx 00661 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx)); 00662 temp1*=xedges; 00663 separableConvolveX(srcImageRange(temp1),destImage(temp2),kernel1d(LTx)); 00664 wx-=sigma*temp2;//(-Lx'*(xedges.*(Lx*u))); 00665 00666 //update wy 00667 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx)); 00668 temp1*=yedges; 00669 separableConvolveY(srcImageRange(temp1),destImage(temp2),kernel1d(LTx)); 00670 wy-=sigma*temp2;//(-Ly'*(yedges.*(Ly*u))); 00671 00672 00673 //update wz 00674 #if (VIGRA_MIXED_2ND_DERIVATIVES) 00675 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx)); 00676 temp1*=yedges; 00677 separableConvolveX(srcImageRange(temp1),destImage(temp2),kernel1d(LTx)); 00678 wz-=sigma*temp2;//-Lx'*(yedges.*(Ly*u)) 00679 00680 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx)); 00681 temp1*=xedges; 00682 separableConvolveY(srcImageRange(temp1),destImage(temp2),kernel1d(LTx)); 00683 wz-=sigma*temp2;//-Ly'*(xedges.*(Lx*u))); 00684 00685 #endif 00686 00687 00688 //project to constraint sets 00689 for (int y=0;y<data.shape(1);y++){ 00690 for (int x=0;x<data.shape(0);x++){ 00691 double e1,e2,skp1,skp2; 00692 00693 //project v 00694 e1=std::cos(phi(x,y)); 00695 e2=std::sin(phi(x,y)); 00696 skp1=vx(x,y)*e1+vy(x,y)*e2; 00697 skp2=vx(x,y)*(-e2)+vy(x,y)*e1; 00698 vigra::detail::projectEllipse2D (skp1,skp2,alpha(x,y),beta(x,y),0.001,100); 00699 vx(x,y)=skp1*e1-skp2*e2; 00700 vy(x,y)=skp1*e2+skp2*e1; 00701 00702 //project w 00703 double l=sqrt(wx(x,y)*wx(x,y)+wy(x,y)*wy(x,y)+wz(x,y)*wz(x,y)); 00704 if (l>gamma(x,y)){ 00705 wx(x,y)=gamma(x,y)*wx(x,y)/l; 00706 wy(x,y)=gamma(x,y)*wy(x,y)/l; 00707 #if (VIGRA_MIXED_2ND_DERIVATIVES) 00708 wz(x,y)=gamma(x,y)*wz(x,y)/l; 00709 #endif 00710 } 00711 } 00712 } 00713 00714 separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx)); 00715 separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx)); 00716 00717 u_bar=out; 00718 out-=tau*(weight*(out-data)+temp1+temp2); 00719 00720 00721 // update wx 00722 separableConvolveX(srcImageRange(wx),destImage(temp1),kernel1d(Lx)); 00723 temp1*=xedges; 00724 separableConvolveX(srcImageRange(temp1),destImage(temp2),kernel1d(LTx)); 00725 out+=tau*temp2; // (-1)^2 00726 00727 00728 //update wy 00729 separableConvolveY(srcImageRange(wy),destImage(temp1),kernel1d(Lx)); 00730 temp1*=yedges; 00731 separableConvolveY(srcImageRange(temp1),destImage(temp2),kernel1d(LTx)); 00732 out+=tau*temp2; 00733 00734 //update wz 00735 #if (VIGRA_MIXED_2ND_DERIVATIVES) 00736 00737 separableConvolveY(srcImageRange(wz),destImage(temp1),kernel1d(Lx)); 00738 temp1*=yedges; 00739 separableConvolveX(srcImageRange(temp1),destImage(temp2),kernel1d(LTx)); 00740 out+=tau*temp2; 00741 00742 separableConvolveX(srcImageRange(wz),destImage(temp1),kernel1d(Lx)); 00743 temp1*=xedges; 00744 separableConvolveY(srcImageRange(temp1),destImage(temp2),kernel1d(LTx)); 00745 out+=tau*temp2; 00746 00747 #endif 00748 00749 u_bar=2*out-u_bar; //cf. Chambolle/Pock and Popov's algorithm 00750 00751 } 00752 } 00753 00754 //@} 00755 } // closing namespace vigra 00756 00757 #endif // VIGRA_TV_FILTER_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|