[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/tv_filter.hxx VIGRA

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&ouml;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&ouml;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)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.9.0 (Tue Nov 6 2012)