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

tv_filter.hxx VIGRA

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