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

tv_filter.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2012 by Frank Lenzen & */
4 /* Ullrich Koethe */
5 /* */
6 /* This file is part of the VIGRA computer vision library. */
7 /* The VIGRA Website is */
8 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
9 /* Please direct questions, bug reports, and contributions to */
10 /* ullrich.koethe@iwr.uni-heidelberg.de or */
11 /* vigra@informatik.uni-hamburg.de */
12 /* */
13 /* Permission is hereby granted, free of charge, to any person */
14 /* obtaining a copy of this software and associated documentation */
15 /* files (the "Software"), to deal in the Software without */
16 /* restriction, including without limitation the rights to use, */
17 /* copy, modify, merge, publish, distribute, sublicense, and/or */
18 /* sell copies of the Software, and to permit persons to whom the */
19 /* Software is furnished to do so, subject to the following */
20 /* conditions: */
21 /* */
22 /* The above copyright notice and this permission notice shall be */
23 /* included in all copies or substantial portions of the */
24 /* Software. */
25 /* */
26 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
27 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
28 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
29 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
30 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
31 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
32 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
33 /* OTHER DEALINGS IN THE SOFTWARE. */
34 /* */
35 /************************************************************************/
36 
37 #ifndef VIGRA_TV_FILTER_HXX
38 #define VIGRA_TV_FILTER_HXX
39 
40 #include <iostream>
41 #include <cmath>
42 #include "config.hxx"
43 #include "impex.hxx"
44 #include "separableconvolution.hxx"
45 #include "multi_array.hxx"
46 #include "multi_math.hxx"
47 #include "eigensystem.hxx"
48 #include "convolution.hxx"
49 #include "fixedpoint.hxx"
50 #include "project2ellipse.hxx"
51 
52 #ifndef VIGRA_MIXED_2ND_DERIVATIVES
53 #define VIGRA_MIXED_2ND_DERIVATIVES 1
54 #endif
55 
56 
57 namespace vigra{
58 
59 
60 
61 /** \addtogroup NonLinearDiffusion
62 */
63 
64 //@{
65 
66 /********************************************************/
67 /* */
68 /* totalVariationFilter */
69 /* */
70 /********************************************************/
71 
72 /** \brief Performs standard Total Variation Regularization
73 
74 The algorithm minimizes
75 
76 \f[
77  \min_u \int_\Omega \frac{1}{2} (u-f)^2\;dx + \alpha TV(u)\qquad\qquad (1)
78 \f]
79 where <em>\f$ f=f(x)\f$</em> are the two dimensional noisy data,
80 <em> \f$ u=u(x)\f$</em> are the smoothed data,<em>\f$ \alpha \ge 0 \f$</em>
81 is the filter parameter and <em>\f$ TV(u)\f$ </em> is the total variation semi-norm.
82 
83 <b> Declarations:</b>
84 
85 \code
86 namespace vigra {
87  template <class stride1,class stride2>
88  void totalVariationFilter(MultiArrayView<2,double,stride1> data,
89  MultiArrayView<2,double,stride2> out,
90  double alpha,
91  int steps,
92  double eps=0);
93  void totalVariationFilter(MultiArrayView<2,double,stride1> data,
94  MultiArrayView<2,double,stride2> weight,
95  MultiArrayView<2,double,stride3> out,
96  double alpha,
97  int steps,
98  double eps=0);
99 }
100 \endcode
101 
102 \ref totalVariationFilter() implements a primal-dual algorithm to solve (1).
103 
104 Input:
105  <table>
106  <tr><td><em>data</em>: </td><td> input data to be smoothed. </td></tr>
107  <tr><td><em>alpha</em>: </td><td> smoothing parameter.</td></tr>
108  <tr><td><em>steps</em>: </td><td> maximal number of iteration steps. </td></tr>
109  <tr><td><em>eps</em>: </td><td> The algorithm stops, if the primal-dual gap is below the threshold <em>eps</em>.
110  </table>
111 
112  Output:
113 
114  <em>out</em> contains the filtered data.
115 
116  In addition, a point-wise weight (\f$ \ge 0 \f$) for the data term can be provided (overloaded function).
117 
118  <b> Usage:</b>
119 
120  <b>\#include</b> <vigra/tv_filter.hxx>
121 
122  \code
123  MultiArray<2,double> data(Shape2(width,height)); // to be initialized
124  MultiArray<2,double> out(Shape2(width,height));
125  MultiArray<2,double> weight(Shape2(width,height))); //optional argument in overloaded function, to be initialized if used
126  int steps; // to be initialized
127  double alpha,eps; // to be initialized
128 
129  totalVariationFilter(data,out,alpha,steps,eps);
130  \endcode
131  or
132  \code
133  totalVariationFilter(data,weight,out,alpha,steps,eps);
134  \endcode
135 
136  */
137 doxygen_overloaded_function(template <...> void totalVariationFilter)
138 
139 template <class stride1,class stride2>
140 void totalVariationFilter(MultiArrayView<2,double,stride1> data,MultiArrayView<2,double,stride2> out, double alpha, int steps, double eps=0){
141 
142  using namespace multi_math;
143 
144  MultiArray<2,double> temp1(data.shape()),temp2(data.shape()),vx(data.shape()),vy(data.shape()),u_bar(data.shape());
145  Kernel1D<double> Lx,LTx;
146  Lx.initExplicitly(0,1)=1,-1; // = Right sided finite differences for d/dx and d/dy
147  Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT); // with hom. Neumann boundary conditions
148  LTx.initExplicitly(-1,0)=-1,1; // = Left sided finite differences for -d/dx and -d/dy
149  LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD); // with hom. Dirichlet b. c.
150 
151  out=data;
152  u_bar=data;
153 
154  double tau=1.0 / std::max(alpha,1.) / std::sqrt(8.0) * 0.06;
155  double sigma=1.0 / std::sqrt(8.0) / 0.06;
156 
157  for (int i=0;i<steps;i++){
158 
159  separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));
160  vx+=(sigma*temp1);
161  separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));
162  vy+=(sigma*temp1);
163 
164  //project to constraint set
165  for (int y=0;y<data.shape(1);y++){
166  for (int x=0;x<data.shape(0);x++){
167  double l=hypot(vx(x,y),vy(x,y));
168  if (l>1){
169  vx(x,y)/=l;
170  vy(x,y)/=l;
171  }
172  }
173  }
174 
175  separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx));
176  separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx));
177  u_bar=out;
178  out-=tau*(out-data+alpha*(temp1+temp2));
179  u_bar=2*out-u_bar; //cf. Chambolle/Pock and Popov's algorithm
180 
181 
182  //stopping criterion
183  if (eps>0){
184  separableConvolveX(srcImageRange(out),destImage(temp1),kernel1d(Lx));
185  separableConvolveY(srcImageRange(out),destImage(temp2),kernel1d(Lx));
186 
187  double f_primal=0,f_dual=0;
188  for (int y=0;y<data.shape(1);y++){
189  for (int x=0;x<data.shape(0);x++){
190  f_primal+=.5*(out(x,y)-data(x,y))*(out(x,y)-data(x,y))+alpha*hypot(temp1(x,y),temp2(x,y));
191  }
192  }
193  separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx));
194  separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx));
195  for (int y=0;y<data.shape(1);y++){
196  for (int x=0;x<data.shape(0);x++){
197  double divv=temp1(x,y)+temp2(x,y);
198  f_dual+=-.5*alpha*alpha*(divv*divv)+alpha*data(x,y)*divv;
199  }
200  }
201  if (f_primal>0 && (f_primal-f_dual)/f_primal<eps){
202  break;
203  }
204  }
205  }
206 }
207 
208 template <class stride1,class stride2, class stride3>
209 void totalVariationFilter(MultiArrayView<2,double,stride1> data,MultiArrayView<2,double,stride2> weight, MultiArrayView<2,double,stride3> out,double alpha, int steps, double eps=0){
210 
211  using namespace multi_math;
212 
213  MultiArray<2,double> temp1(data.shape()),temp2(data.shape()),vx(data.shape()),vy(data.shape()),u_bar(data.shape());
214  Kernel1D<double> Lx,LTx;
215  Lx.initExplicitly(0,1)=1,-1; // = Right sided finite differences for d/dx and d/dy
216  Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT); // with hom. Neumann boundary conditions
217  LTx.initExplicitly(-1,0)=-1,1; // = Left sided finite differences for -d/dx and -d/dy
218  LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD); // with hom. Dirichlet b. c.
219 
220  out=data;
221  u_bar=data;
222 
223  double tau=1.0 / std::max(alpha,1.) / std::sqrt(8.0) * 0.06;
224  double sigma=1.0 / std::sqrt(8.0) / 0.06;
225 
226  for (int i=0;i<steps;i++){
227  separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));
228  vx+=(sigma*temp1);
229  separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));
230  vy+=(sigma*temp1);
231 
232  //project to constraint set
233  for (int y=0;y<data.shape(1);y++){
234  for (int x=0;x<data.shape(0);x++){
235  double l=hypot(vx(x,y),vy(x,y));
236  if (l>1){
237  vx(x,y)/=l;
238  vy(x,y)/=l;
239  }
240  }
241  }
242 
243  separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx));
244  separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx));
245  u_bar=out;
246  out-=tau*(weight*(out-data)+alpha*(temp1+temp2));
247  u_bar=2*out-u_bar;
248 
249 
250  //stopping criterion
251  if (eps>0){
252  separableConvolveX(srcImageRange(out),destImage(temp1),kernel1d(Lx));
253  separableConvolveY(srcImageRange(out),destImage(temp2),kernel1d(Lx));
254 
255  double f_primal=0,f_dual=0;
256  for (int y=0;y<data.shape(1);y++){
257  for (int x=0;x<data.shape(0);x++){
258  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));
259  }
260  }
261  separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx));
262  separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx));
263  for (int y=0;y<data.shape(1);y++){
264  for (int x=0;x<data.shape(0);x++){
265  double divv=temp1(x,y)+temp2(x,y);
266  f_dual+=-.5*alpha*alpha*(weight(x,y)*divv*divv)+alpha*data(x,y)*divv;
267  }
268  }
269  if (f_primal>0 && (f_primal-f_dual)/f_primal<eps){
270  break;
271  }
272  }
273  }
274 }
275 //<!--\f$ \alpha(x)=\beta(x)=\beta_{par}\f$ in homogeneous regions without edges,
276 //and \f$ \alpha(x)=\alpha_{par}\f$ at edges.-->
277 
278 
279 /********************************************************/
280 /* */
281 /* getAnisotropy */
282 /* */
283 /********************************************************/
284 
285 /** \brief Sets up directional data for anisotropic regularization
286 
287 This routine provides a two-dimensional normalized vector field \f$ v \f$, which is normal to edges in the given data,
288 found as the eigenvector of the structure tensor belonging to the largest eigenvalue.
289 \f$ v \f$ is encoded by a scalar field \f$ \varphi \f$ of angles, i.e.
290 \f$ v(x)=(\cos(\varphi(x)),\sin(\varphi(x)))^\top \f$.
291 
292 In addition, two scalar fields \f$ \alpha \f$ and \f$ \beta \f$ are generated from
293 scalar parameters \f$ \alpha_{par}\f$ and \f$ \beta_{par}\f$, such that
294 <center>
295 <table>
296 <tr> <td>\f$ \alpha(x)= \alpha_{par}\f$ at edges,</td></tr>
297 <tr> <td>\f$ \alpha(x)= \beta_{par}\f$ in homogeneous regions,</td></tr>
298 <tr> <td>\f$ \beta(x)=\beta_{par}\f$ .</td></tr>
299 </table>
300 </center>
301 
302 <b> Declarations:</b>
303 
304 \code
305 namespace vigra {
306 void getAnisotropy(MultiArrayView<2,double,stride1> data,
307  MultiArrayView<2,double,stride2> phi,
308  MultiArrayView<2,double,stride3> alpha,
309  MultiArrayView<2,double,stride4> beta,
310  double alpha_par,
311  double beta_par,
312  double sigma_par,
313  double rho_par,
314  double K_par);
315 }
316 \endcode
317 
318 Output:
319 <table>
320  <tr><td>Three scalar fields <em>phi</em>, <em>alpha</em> and <em>beta</em>.</td></tr>
321 </table>
322 
323 Input:
324 <table>
325 <tr><td><em>data</em>:</td><td>two-dimensional scalar field.</td></tr>
326 <tr><td><em>alpha_par,beta_par</em>:</td><td>two positive values for setting up the scalar fields alpha and beta</tr>
327 <tr><td><em>sigma_par</em>:</td><td> non-negative parameter for presmoothing the data.</td></tr>
328 <tr><td> <em>rho_par</em>:</td><td> non-negative parameter for presmoothing the structure tensor.</td></tr>
329 <tr><td><em>K_par</em>:</td><td> positive edge sensitivity parameter.</td></tr>
330  </table>
331 
332 (see \ref anisotropicTotalVariationFilter() and \ref secondOrderTotalVariationFilter() for usage in an application).
333 */
334 doxygen_overloaded_function(template <...> void getAnisotropy)
335 
336 template <class stride1,class stride2,class stride3,class stride4>
337 void getAnisotropy(MultiArrayView<2,double,stride1> data,MultiArrayView<2,double,stride2> phi,
338  MultiArrayView<2,double,stride3> alpha, MultiArrayView<2,double,stride4> beta,
339  double alpha_par, double beta_par, double sigma_par, double rho_par, double K_par){
340 
341  using namespace multi_math;
342 
343  MultiArray<2,double> smooth(data.shape()),tmp(data.shape());
345 
346 
347  gauss.initGaussian(sigma_par);
348  separableConvolveX(srcImageRange(data), destImage(tmp), kernel1d(gauss));
349  separableConvolveY(srcImageRange(tmp), destImage(smooth), kernel1d(gauss));
350 
351  MultiArray<2,double> stxx(data.shape()),stxy(data.shape()),styy(data.shape());
352 
353  // calculate Structure Tensor at inner scale = sigma and outer scale = rho
354  vigra::structureTensor(srcImageRange(smooth),destImage(stxx), destImage(stxy), destImage(styy),1.,1.);
355 
356  gauss.initGaussian(rho_par);
357  separableConvolveX(srcImageRange(stxx), destImage(tmp), kernel1d(gauss));
358  separableConvolveY(srcImageRange(tmp), destImage(stxx), kernel1d(gauss));
359  separableConvolveX(srcImageRange(stxy), destImage(tmp), kernel1d(gauss));
360  separableConvolveY(srcImageRange(tmp), destImage(stxy), kernel1d(gauss));
361  separableConvolveX(srcImageRange(styy), destImage(tmp), kernel1d(gauss));
362  separableConvolveY(srcImageRange(tmp), destImage(styy), kernel1d(gauss));
363 
364  MultiArray<2,double> matrix(Shape2(2,2)),ev(Shape2(2,2)),ew(Shape2(2,1));
365 
366  for (int y=0;y<data.shape(1);y++){
367  for (int x=0;x<data.shape(0);x++){
368 
369  matrix(0,0)=stxx(x,y);
370  matrix(1,1)=styy(x,y);
371  matrix(0,1)=stxy(x,y);
372  matrix(1,0)=stxy(x,y);
373  vigra::linalg::detail::symmetricEigensystem2x2(matrix,ew,ev);
374 
375  phi(x,y)=std::atan2(ev(1,0),ev(0,0));
376  double coherence=ew(0,0)-ew(1,0);
377  double c=std::min(K_par*coherence,1.);
378  alpha(x,y)=alpha_par*c+(1-c)*beta_par;
379  beta(x,y)=beta_par;
380  }
381  }
382 }
383 
384 /********************************************************/
385 /* */
386 /* anisotropicTotalVariationFilter */
387 /* */
388 /********************************************************/
389 
390 /** \brief Performs Anisotropic Total Variation Regularization
391 
392 The algorithm minimizes
393 \f[
394 \min_u \int_\Omega \frac{1}{2} (u-f)^2 + \sqrt{\nabla u^\top A \nabla u}\;dx\qquad\qquad(2)
395 \f]
396 
397 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>
398 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.
399 
400 Matrix <em>\f$ A \f$ </em> is described by providing for each data point a normalized eigenvector (via angle \f$ \phi \f$)
401 and two eigenvalues \f$ \alpha>0 \f$ and \f$ \beta>0 \f$.
402 
403 \ref getAnisotropy() can be use to set up such data \f$ \phi,\alpha,\beta \f$ by providing a vector field normal to edges.
404 
405 <b> Declarations:</b>
406 
407 \code
408 namespace vigra {
409  template <class stride1,class stride2,class stride3,class stride4,class stride5,class stride6>
410  void anisotropicTotalVariationFilter(MultiArrayView<2,double,stride1> data,
411  MultiArrayView<2,double,stride2> weight,
412  MultiArrayView<2,double,stride3> phi,
413  MultiArrayView<2,double,stride4> alpha,
414  MultiArrayView<2,double,stride5> beta,
415  MultiArrayView<2,double,stride6> out,
416  int steps);
417 }
418 \endcode
419 
420 \ref anisotropicTotalVariationFilter() implements a primal-dual algorithm to solve (2).
421 
422 Input:
423 <table>
424 <tr><td><em>data</em>:</td><td>input data to be filtered. </td></tr>
425 <tr><td><em>steps</em>:</td><td>iteration steps.</td></tr>
426 <tr><td><em>weight</em> :</td><td>a point-wise weight (\f$ \ge 0 \f$ ) for the data term.</td></tr>
427 <tr><td><em>phi</em>,<em>alpha</em> and <em>beta</em> :</td><td> describe matrix \f$ A \f$, see above.</td></tr>
428 </table>
429 
430 Output:
431 <table>
432 <tr><td><em>out</em> :</td><td>contains filtered data.</td></tr>
433 </table>
434 
435 <b> Usage:</b>
436 
437 E.g. with a solution-dependent adaptivity cf. [1], by updating the matrix \f$ A=A(u)\f$
438 in an outer loop:
439 
440 <b>\#include</b> <vigra/tv_filter.hxx>
441 
442 \code
443 MultiArray<2,double> data(Shape2(width,height)); //to be initialized
444 MultiArray<2,double> out (Shape2(width,height));
445 MultiArray<2,double> weight(Shape2(width,height)); //to be initialized
446 MultiArray<2,double> phi (Shape2(width,height));
447 MultiArray<2,double> alpha(Shape2(width,height));
448 MultiArray<2,double> beta (Shape2(width,height));
449 double alpha0,beta0,sigma,rho,K; //to be initialized
450 int outer_steps,inner_steps;//to be initialized
451 
452 out=data; // data serves as initial value
453 
454 for (int i=0;i<outer_steps;i++){
455  getAnisotropy(out,phi,alpha,beta,alpha0,beta0,sigma,rho,K); // sets phi, alpha, beta
456  anisotropicTotalVariationFilter(data,weight,phi,alpha,beta,out,inner_steps);
457  }
458 \endcode
459 
460 [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.
461 */
462 doxygen_overloaded_function(template <...> void anisotropicTotalVariationFilter)
463 
464 template <class stride1,class stride2,class stride3,class stride4,class stride5,class stride6>
465 void anisotropicTotalVariationFilter(MultiArrayView<2,double,stride1> data,MultiArrayView<2,double,stride2> weight,
466  MultiArrayView<2,double,stride3> phi,MultiArrayView<2,double,stride4> alpha,
467  MultiArrayView<2,double,stride5> beta,MultiArrayView<2,double,stride6> out,
468  int steps){
469 
470  using namespace multi_math;
471 
472  MultiArray<2,double> temp1(data.shape()),temp2(data.shape()),vx(data.shape()),vy(data.shape()),u_bar(data.shape());
473  MultiArray<2,double> rx(data.shape()),ry(data.shape());
474 
475  Kernel1D<double> Lx,LTx;
476  Lx.initExplicitly(0,1)=1,-1; // = Right sided finite differences for d/dx and d/dy
477  Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT); // with hom. Neumann boundary conditions
478  LTx.initExplicitly(-1,0)=-1,1; // = Left sided finite differences for -d/dx and -d/dy
479  LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD); // with hom. Dirichlet b. c.
480 
481  u_bar=out;
482 
483  double m=0;
484  for (int y=0;y<data.shape(1);y++){
485  for (int x=0;x<data.shape(0);x++){
486  m=std::max(m,alpha(x,y));
487  m=std::max(m,beta (x,y));
488  }
489  }
490  m=std::max(m,1.);
491  double tau=.9/m/std::sqrt(8.)*0.06;
492  double sigma=.9/m/std::sqrt(8.)/0.06;
493 
494 
495  for (int i=0;i<steps;i++){
496  separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));
497  vx+=(sigma*temp1);
498  separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));
499  vy+=(sigma*temp1);
500 
501  //project to constraint set
502  for (int y=0;y<data.shape(1);y++){
503  for (int x=0;x<data.shape(0);x++){
504  double e1,e2,skp1,skp2;
505 
506  e1=std::cos(phi(x,y));
507  e2=std::sin(phi(x,y));
508  skp1=vx(x,y)*e1+vy(x,y)*e2;
509  skp2=vx(x,y)*(-e2)+vy(x,y)*e1;
510  vigra::detail::projectEllipse2D (skp1,skp2,alpha(x,y),beta(x,y),0.001,100);
511 
512  vx(x,y)=skp1*e1-skp2*e2;
513  vy(x,y)=skp1*e2+skp2*e1;
514  }
515  }
516 
517  separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx));
518  separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx));
519  u_bar=out;
520  out-=tau*(weight*(out-data)+(temp1+temp2));
521  u_bar=2*out-u_bar; //cf. Chambolle/Pock and Popov's algorithm
522  }
523 }
524 
525 /********************************************************/
526 /* */
527 /* secondOrderTotalVariationFilter */
528 /* */
529 /********************************************************/
530 
531 /** \brief Performs Anisotropic Total Variation Regularization
532 
533 The algorithm minimizes
534 
535 \f[
536 \min_u \int_\Omega \frac{1}{2} (u-f)^2 + \sqrt{\nabla u^\top A \nabla u} + \gamma |Hu|_F\;dx \qquad\qquad (3)
537 \f]
538 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>
539 is the image gradient in the sense of Total Variation, <em>\f$ A \f$ </em> is a locally varying
540 symmetric, positive-definite 2x2 matrix and <em>\f$ |Hu|_F \f$</em> is the Frobenius norm of the Hessian of \f$ u \f$.
541 
542 Matrix <em>\f$ A \f$ </em> is described by providing for each data point a normalized eigenvector (via angle \f$ \phi \f$)
543 and two eigenvalues \f$ \alpha>0 \f$ and \f$ \beta>0 \f$.
544 \ref getAnisotropy() can be use to set up such data \f$ \phi,\alpha, \beta \f$ by providing a vector field normal to edges.
545 
546 \f$ \gamma>0 \f$ is the locally varying regularization parameter for second order.
547 
548 <b> Declarations:</b>
549 
550 \code
551 namespace vigra {
552  template <class stride1,class stride2,...,class stride9>
553  void secondOrderTotalVariationFilter(MultiArrayView<2,double,stride1> data,
554  MultiArrayView<2,double,stride2> weight,
555  MultiArrayView<2,double,stride3> phi,
556  MultiArrayView<2,double,stride4> alpha,
557  MultiArrayView<2,double,stride5> beta,
558  MultiArrayView<2,double,stride6> gamma,
559  MultiArrayView<2,double,stride7> xedges,
560  MultiArrayView<2,double,stride8> yedges,
561  MultiArrayView<2,double,stride9> out,
562  int steps);
563 }
564 \endcode
565 
566 \ref secondOrderTotalVariationFilter() implements a primal-dual algorithm to solve (3).
567 
568 Input:
569 <table>
570 <tr><td><em>data</em>: </td><td> the input data to be filtered. </td></tr>
571 <tr><td><em>steps</em> : </td><td> number of iteration steps.</td></tr>
572 <tr><td><em>out</em> : </td><td>contains the filtered data.</td></tr>
573 <tr><td><em>weight</em> : </td><td> point-wise weight (\f$ \ge 0\f$ ) for the data term.</td></tr>
574 <tr><td><em>phi</em>,<em>alpha</em>,<em>beta</em>: </td><td> describe matrix \f$ A\f$, see above.</td></tr>
575 <tr><td><em> xedges </em> and <em> yedges </em>: </td><td> binary arrays indicating the
576 presence of horizontal (between (x,y) and (x+1,y)) and vertical edges (between (x,y) and (x,y+1)).
577 These data are considered in the calculation of \f$ Hu\f$, such that
578 finite differences across edges are artificially set to zero to avoid second order smoothing over edges.</td></tr>
579 </table>
580 
581 <b> Usage:</b>
582 
583 E.g. with a solution-dependent adaptivity (cf.[1]), by updating the matrix \f$ A=A(u)\f$
584 in an outer loop:
585 
586 <b>\#include</b> <vigra/tv_filter.hxx>
587 
588 \code
589 MultiArray<2,double> data(Shape2(width,height)); //to be initialized
590 MultiArray<2,double> out(Shape2(width,height));
591 MultiArray<2,double> weight(Shape2(width,height))); //to be initialized
592 MultiArray<2,double> phi(Shape2(width,height);
593 MultiArray<2,double> alpha(Shape2(width,height);
594 MultiArray<2,double> beta(Shape2(width,height));
595 MultiArray<2,double> gamma(Shape2(width,height)); //to be initialized
596 MultiArray<2,double> xedges(Shape2(width,height)); //to be initialized
597 MultiArray<2,double> yedges(Shape2(width,height)); //to be initialized
598 
599 
600 double alpha0,beta0,sigma,rho,K; //to be initialized
601 int outer_steps,inner_steps;//to be initialized
602 
603 out=data; // data serves as initial value
604 
605 for (int i=0;i<outer_steps;i++){
606 
607  getAnisotropy(out,phi,alpha,beta,alpha0,beta0,sigma,rho,K); // sets phi, alpha, beta
608  secondOrderTotalVariationFilter(data,weight,phi,alpha,beta,gamma,xedges,yedges,out,inner_steps);
609 }
610 \endcode
611 
612 
613 [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.
614 */
615 doxygen_overloaded_function(template <...> void secondOrderTotalVariationFilter)
616 
617 template <class stride1,class stride2,class stride3,class stride4,class stride5,class stride6,class stride7,class stride8,class stride9>
618 void secondOrderTotalVariationFilter(MultiArrayView<2,double,stride1> data,
619  MultiArrayView<2,double,stride2> weight,MultiArrayView<2,double,stride3> phi,
620  MultiArrayView<2,double,stride4> alpha,MultiArrayView<2,double,stride5> beta,
621  MultiArrayView<2,double,stride6> gamma,
622  MultiArrayView<2,double,stride7> xedges,MultiArrayView<2,double,stride8> yedges,
623  MultiArrayView<2,double,stride9> out,
624  int steps){
625 
626  using namespace multi_math;
627 
628  MultiArray<2,double> temp1(data.shape()),temp2(data.shape()),vx(data.shape()),vy(data.shape()),u_bar(data.shape());
629  MultiArray<2,double> rx(data.shape()),ry(data.shape());
630  MultiArray<2,double> wx(data.shape()),wy(data.shape()),wz(data.shape());
631 
632 
633  Kernel1D<double> Lx,LTx;
634  Lx.initExplicitly(0,1)=1,-1; // = Right sided finite differences for d/dx and d/dy
635  Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT); // with hom. Neumann boundary conditions
636  LTx.initExplicitly(-1,0)=-1,1; // = Left sided finite differences for -d/dx and -d/dy
637  LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD); // with hom. Dirichlet b. c.
638 
639  u_bar=out;
640 
641  double m=0;
642  for (int y=0;y<data.shape(1);y++){
643  for (int x=0;x<data.shape(0);x++){
644  m=std::max(m,alpha(x,y));
645  m=std::max(m,beta (x,y));
646  }
647  }
648  m=std::max(m,1.);
649  double tau=.1/m;//std::sqrt(8)*0.06;
650  double sigma=.1;//m;/std::sqrt(8)/0.06;
651 
652  for (int i=0;i<steps;i++){
653 
654  separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));
655  vx+=(sigma*temp1);
656  separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));
657  vy+=(sigma*temp1);
658 
659 
660  // update wx
661  separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));
662  temp1*=xedges;
663  separableConvolveX(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
664  wx-=sigma*temp2;//(-Lx'*(xedges.*(Lx*u)));
665 
666  //update wy
667  separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));
668  temp1*=yedges;
669  separableConvolveY(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
670  wy-=sigma*temp2;//(-Ly'*(yedges.*(Ly*u)));
671 
672 
673  //update wz
674  #if (VIGRA_MIXED_2ND_DERIVATIVES)
675  separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));
676  temp1*=yedges;
677  separableConvolveX(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
678  wz-=sigma*temp2;//-Lx'*(yedges.*(Ly*u))
679 
680  separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));
681  temp1*=xedges;
682  separableConvolveY(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
683  wz-=sigma*temp2;//-Ly'*(xedges.*(Lx*u)));
684 
685  #endif
686 
687 
688  //project to constraint sets
689  for (int y=0;y<data.shape(1);y++){
690  for (int x=0;x<data.shape(0);x++){
691  double e1,e2,skp1,skp2;
692 
693  //project v
694  e1=std::cos(phi(x,y));
695  e2=std::sin(phi(x,y));
696  skp1=vx(x,y)*e1+vy(x,y)*e2;
697  skp2=vx(x,y)*(-e2)+vy(x,y)*e1;
698  vigra::detail::projectEllipse2D (skp1,skp2,alpha(x,y),beta(x,y),0.001,100);
699  vx(x,y)=skp1*e1-skp2*e2;
700  vy(x,y)=skp1*e2+skp2*e1;
701 
702  //project w
703  double l=sqrt(wx(x,y)*wx(x,y)+wy(x,y)*wy(x,y)+wz(x,y)*wz(x,y));
704  if (l>gamma(x,y)){
705  wx(x,y)=gamma(x,y)*wx(x,y)/l;
706  wy(x,y)=gamma(x,y)*wy(x,y)/l;
707  #if (VIGRA_MIXED_2ND_DERIVATIVES)
708  wz(x,y)=gamma(x,y)*wz(x,y)/l;
709  #endif
710  }
711  }
712  }
713 
714  separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx));
715  separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx));
716 
717  u_bar=out;
718  out-=tau*(weight*(out-data)+temp1+temp2);
719 
720 
721  // update wx
722  separableConvolveX(srcImageRange(wx),destImage(temp1),kernel1d(Lx));
723  temp1*=xedges;
724  separableConvolveX(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
725  out+=tau*temp2; // (-1)^2
726 
727 
728  //update wy
729  separableConvolveY(srcImageRange(wy),destImage(temp1),kernel1d(Lx));
730  temp1*=yedges;
731  separableConvolveY(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
732  out+=tau*temp2;
733 
734  //update wz
735  #if (VIGRA_MIXED_2ND_DERIVATIVES)
736 
737  separableConvolveY(srcImageRange(wz),destImage(temp1),kernel1d(Lx));
738  temp1*=yedges;
739  separableConvolveX(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
740  out+=tau*temp2;
741 
742  separableConvolveX(srcImageRange(wz),destImage(temp1),kernel1d(Lx));
743  temp1*=xedges;
744  separableConvolveY(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
745  out+=tau*temp2;
746 
747  #endif
748 
749  u_bar=2*out-u_bar; //cf. Chambolle/Pock and Popov's algorithm
750 
751  }
752 }
753 
754 //@}
755 } // closing namespace vigra
756 
757 #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 Oct 22 2013)