[ 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&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)