[ 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#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
59namespace 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
76The 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]
81where <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>
83is the filter parameter and <em>\f$ TV(u)\f$ </em> is the total variation semi-norm.
84
85<b> Declarations:</b>
86
87\code
88namespace 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
106Input:
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 */
139doxygen_overloaded_function(template <...> void totalVariationFilter)
140
141template <class stride1,class stride2>
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());
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
211template <class stride1,class stride2, class stride3>
212void 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
291This routine provides a two-dimensional normalized vector field \f$ v \f$, which is normal to edges in the given data,
292found 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
296In addition, two scalar fields \f$ \alpha \f$ and \f$ \beta \f$ are generated from
297scalar 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
309namespace vigra {
310void 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
322Output:
323<table>
324 <tr><td>Three scalar fields <em>phi</em>, <em>alpha</em> and <em>beta</em>.</td></tr>
325</table>
326
327Input:
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*/
338doxygen_overloaded_function(template <...> void getAnisotropy)
339
340template <class stride1,class stride2,class stride3,class stride4>
343 double alpha_par, double beta_par, double sigma_par, double rho_par, double K_par){
344
345 using namespace multi_math;
346
347 MultiArray<2,double> smooth(data.shape()),tmp(data.shape());
349
350
351 gauss.initGaussian(sigma_par);
352 separableConvolveX(srcImageRange(data), destImage(tmp), kernel1d(gauss));
353 separableConvolveY(srcImageRange(tmp), destImage(smooth), kernel1d(gauss));
354
355 MultiArray<2,double> stxx(data.shape()),stxy(data.shape()),styy(data.shape());
356
357 // calculate Structure Tensor at inner scale = sigma and outer scale = rho
358 vigra::structureTensor(srcImageRange(smooth),destImage(stxx), destImage(stxy), destImage(styy),1.,1.);
359
360 gauss.initGaussian(rho_par);
361 separableConvolveX(srcImageRange(stxx), destImage(tmp), kernel1d(gauss));
362 separableConvolveY(srcImageRange(tmp), destImage(stxx), kernel1d(gauss));
363 separableConvolveX(srcImageRange(stxy), destImage(tmp), kernel1d(gauss));
364 separableConvolveY(srcImageRange(tmp), destImage(stxy), kernel1d(gauss));
365 separableConvolveX(srcImageRange(styy), destImage(tmp), kernel1d(gauss));
366 separableConvolveY(srcImageRange(tmp), destImage(styy), kernel1d(gauss));
367
368 MultiArray<2,double> matrix(Shape2(2,2)),ev(Shape2(2,2)),ew(Shape2(2,1));
369
370 for (int y=0;y<data.shape(1);y++){
371 for (int x=0;x<data.shape(0);x++){
372
373 matrix(0,0)=stxx(x,y);
374 matrix(1,1)=styy(x,y);
375 matrix(0,1)=stxy(x,y);
376 matrix(1,0)=stxy(x,y);
377 vigra::symmetricEigensystemNoniterative(matrix,ew,ev);
378
379 phi(x,y)=std::atan2(ev(1,0),ev(0,0));
380 double coherence=ew(0,0)-ew(1,0);
381 double c=std::min(K_par*coherence,1.);
382 alpha(x,y)=alpha_par*c+(1-c)*beta_par;
383 beta(x,y)=beta_par;
384 }
385 }
386}
387
388/********************************************************/
389/* */
390/* anisotropicTotalVariationFilter */
391/* */
392/********************************************************/
393
394/** \brief Performs Anisotropic Total Variation Regularization
395
396The algorithm minimizes
397\f[
398\min_u \int_\Omega \frac{1}{2} (u-f)^2 + \sqrt{\nabla u^\top A \nabla u}\;dx\qquad\qquad(2)
399\f]
400
401where <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>
402is the image gradient in the sense of Total Variation and <em>\f$ A \f$ </em> is a locally varying symmetric, positive definite 2x2 matrix.
403
404Matrix <em>\f$ A \f$ </em> is described by providing for each data point a normalized eigenvector (via angle \f$ \phi \f$)
405and two eigenvalues \f$ \alpha>0 \f$ and \f$ \beta>0 \f$.
406
407\ref getAnisotropy() can be use to set up such data \f$ \phi,\alpha,\beta \f$ by providing a vector field normal to edges.
408
409<b> Declarations:</b>
410
411\code
412namespace vigra {
413 template <class stride1,class stride2,class stride3,class stride4,class stride5,class stride6>
414 void anisotropicTotalVariationFilter(MultiArrayView<2,double,stride1> data,
415 MultiArrayView<2,double,stride2> weight,
416 MultiArrayView<2,double,stride3> phi,
417 MultiArrayView<2,double,stride4> alpha,
418 MultiArrayView<2,double,stride5> beta,
419 MultiArrayView<2,double,stride6> out,
420 int steps);
421}
422\endcode
423
424\ref anisotropicTotalVariationFilter() implements a primal-dual algorithm to solve (2).
425
426Input:
427<table>
428<tr><td><em>data</em>:</td><td>input data to be filtered. </td></tr>
429<tr><td><em>steps</em>:</td><td>iteration steps.</td></tr>
430<tr><td><em>weight</em> :</td><td>a point-wise weight (\f$ \ge 0 \f$ ) for the data term.</td></tr>
431<tr><td><em>phi</em>,<em>alpha</em> and <em>beta</em> :</td><td> describe matrix \f$ A \f$, see above.</td></tr>
432</table>
433
434Output:
435<table>
436<tr><td><em>out</em> :</td><td>contains filtered data.</td></tr>
437</table>
438
439<b> Usage:</b>
440
441E.g. with a solution-dependent adaptivity cf. [1], by updating the matrix \f$ A=A(u)\f$
442in an outer loop:
443
444<b>\#include</b> <vigra/tv_filter.hxx>
445
446\code
447MultiArray<2,double> data(Shape2(width,height)); //to be initialized
448MultiArray<2,double> out (Shape2(width,height));
449MultiArray<2,double> weight(Shape2(width,height)); //to be initialized
450MultiArray<2,double> phi (Shape2(width,height));
451MultiArray<2,double> alpha(Shape2(width,height));
452MultiArray<2,double> beta (Shape2(width,height));
453double alpha0,beta0,sigma,rho,K; //to be initialized
454int outer_steps,inner_steps;//to be initialized
455
456out=data; // data serves as initial value
457
458for (int i=0;i<outer_steps;i++){
459 getAnisotropy(out,phi,alpha,beta,alpha0,beta0,sigma,rho,K); // sets phi, alpha, beta
460 anisotropicTotalVariationFilter(data,weight,phi,alpha,beta,out,inner_steps);
461 }
462\endcode
463
464[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.
465*/
466doxygen_overloaded_function(template <...> void anisotropicTotalVariationFilter)
467
468template <class stride1,class stride2,class stride3,class stride4,class stride5,class stride6>
472 int steps){
473
474 using namespace multi_math;
475 int width=data.shape(0),height=data.shape(1);
476
477 MultiArray<2,double> temp1(data.shape()),temp2(data.shape()),vx(data.shape()),vy(data.shape()),u_bar(data.shape());
478 MultiArray<2,double> rx(data.shape()),ry(data.shape());
479
481 Lx.initExplicitly(-1,0)=1,-1; // = Right sided finite differences for d/dx and d/dy
482 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT); // with hom. Neumann boundary conditions
483 LTx.initExplicitly(0,1)=-1,1; // = Left sided finite differences for -d/dx and -d/dy
484 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD); // with hom. Dirichlet b. c.
485
486 u_bar=out;
487
488 double m=0;
489 for (int y=0;y<data.shape(1);y++){
490 for (int x=0;x<data.shape(0);x++){
491 m=std::max(m,alpha(x,y));
492 m=std::max(m,beta (x,y));
493 }
494 }
495 m=std::max(m,1.);
496 double tau=.9/m/std::sqrt(8.)*0.06;
497 double sigma=.9/m/std::sqrt(8.)/0.06;
498
499
500 for (int i=0;i<steps;i++){
501 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroX(temp1);
502 vx+=(sigma*temp1);
503 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroY(temp1);
504 vy+=(sigma*temp1);
505
506 //project to constraint set
507 for (int y=0;y<data.shape(1);y++){
508 for (int x=0;x<data.shape(0);x++){
509 double e1,e2,skp1,skp2;
510
511 e1=std::cos(phi(x,y));
512 e2=std::sin(phi(x,y));
513 skp1=vx(x,y)*e1+vy(x,y)*e2;
514 skp2=vx(x,y)*(-e2)+vy(x,y)*e1;
515 vigra::detail::projectEllipse2D (skp1,skp2,alpha(x,y),beta(x,y),0.001,100);
516
517 vx(x,y)=skp1*e1-skp2*e2;
518 vy(x,y)=skp1*e2+skp2*e1;
519 }
520 }
521
522 separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx));
523 separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx));
524 u_bar=out;
525 out-=tau*(weight*(out-data)+(temp1+temp2));
526 u_bar=2*out-u_bar; //cf. Chambolle/Pock and Popov's algorithm
527 }
528}
529
530/********************************************************/
531/* */
532/* secondOrderTotalVariationFilter */
533/* */
534/********************************************************/
535
536/** \brief Performs Anisotropic Total Variation Regularization
537
538The algorithm minimizes
539
540\f[
541\min_u \int_\Omega \frac{1}{2} (u-f)^2 + \sqrt{\nabla u^\top A \nabla u} + \gamma |Hu|_F\;dx \qquad\qquad (3)
542\f]
543where <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>
544is the image gradient in the sense of Total Variation, <em>\f$ A \f$ </em> is a locally varying
545symmetric, positive-definite 2x2 matrix and <em>\f$ |Hu|_F \f$</em> is the Frobenius norm of the Hessian of \f$ u \f$.
546
547Matrix <em>\f$ A \f$ </em> is described by providing for each data point a normalized eigenvector (via angle \f$ \phi \f$)
548and two eigenvalues \f$ \alpha>0 \f$ and \f$ \beta>0 \f$.
549\ref getAnisotropy() can be use to set up such data \f$ \phi,\alpha, \beta \f$ by providing a vector field normal to edges.
550
551\f$ \gamma>0 \f$ is the locally varying regularization parameter for second order.
552
553<b> Declarations:</b>
554
555\code
556namespace vigra {
557 template <class stride1,class stride2,...,class stride9>
558 void secondOrderTotalVariationFilter(MultiArrayView<2,double,stride1> data,
559 MultiArrayView<2,double,stride2> weight,
560 MultiArrayView<2,double,stride3> phi,
561 MultiArrayView<2,double,stride4> alpha,
562 MultiArrayView<2,double,stride5> beta,
563 MultiArrayView<2,double,stride6> gamma,
564 MultiArrayView<2,double,stride7> xedges,
565 MultiArrayView<2,double,stride8> yedges,
566 MultiArrayView<2,double,stride9> out,
567 int steps);
568}
569\endcode
570
571\ref secondOrderTotalVariationFilter() implements a primal-dual algorithm to solve (3).
572
573Input:
574<table>
575<tr><td><em>data</em>: </td><td> the input data to be filtered. </td></tr>
576<tr><td><em>steps</em> : </td><td> number of iteration steps.</td></tr>
577<tr><td><em>out</em> : </td><td>contains the filtered data.</td></tr>
578<tr><td><em>weight</em> : </td><td> point-wise weight (\f$ \ge 0\f$ ) for the data term.</td></tr>
579<tr><td><em>phi</em>,<em>alpha</em>,<em>beta</em>: </td><td> describe matrix \f$ A\f$, see above.</td></tr>
580<tr><td><em> xedges </em> and <em> yedges </em>: </td><td> binary arrays indicating the
581presence of horizontal (between (x,y) and (x+1,y)) and vertical edges (between (x,y) and (x,y+1)).
582These data are considered in the calculation of \f$ Hu\f$, such that
583finite differences across edges are artificially set to zero to avoid second order smoothing over edges.</td></tr>
584</table>
585
586<b> Usage:</b>
587
588E.g. with a solution-dependent adaptivity (cf.[1]), by updating the matrix \f$ A=A(u)\f$
589in an outer loop:
590
591<b>\#include</b> <vigra/tv_filter.hxx>
592
593\code
594MultiArray<2,double> data(Shape2(width,height)); //to be initialized
595MultiArray<2,double> out(Shape2(width,height));
596MultiArray<2,double> weight(Shape2(width,height))); //to be initialized
597MultiArray<2,double> phi(Shape2(width,height);
598MultiArray<2,double> alpha(Shape2(width,height);
599MultiArray<2,double> beta(Shape2(width,height));
600MultiArray<2,double> gamma(Shape2(width,height)); //to be initialized
601MultiArray<2,double> xedges(Shape2(width,height)); //to be initialized
602MultiArray<2,double> yedges(Shape2(width,height)); //to be initialized
603
604
605double alpha0,beta0,sigma,rho,K; //to be initialized
606int outer_steps,inner_steps;//to be initialized
607
608out=data; // data serves as initial value
609
610for (int i=0;i<outer_steps;i++){
611
612 getAnisotropy(out,phi,alpha,beta,alpha0,beta0,sigma,rho,K); // sets phi, alpha, beta
613 secondOrderTotalVariationFilter(data,weight,phi,alpha,beta,gamma,xedges,yedges,out,inner_steps);
614}
615\endcode
616
617
618[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.
619*/
620doxygen_overloaded_function(template <...> void secondOrderTotalVariationFilter)
621
622template <class stride1,class stride2,class stride3,class stride4,class stride5,class stride6,class stride7,class stride8,class stride9>
629 int steps){
630
631 using namespace multi_math;
632 int width=data.shape(0),height=data.shape(1);
633
634 MultiArray<2,double> temp1(data.shape()),temp2(data.shape()),vx(data.shape()),vy(data.shape()),u_bar(data.shape());
635 MultiArray<2,double> rx(data.shape()),ry(data.shape());
636 MultiArray<2,double> wx(data.shape()),wy(data.shape()),wz(data.shape());
637
638
640 Lx.initExplicitly(-1,0)=1,-1; // = Right sided finite differences for d/dx and d/dy
641 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT); // with hom. Neumann boundary conditions
642 LTx.initExplicitly(0,1)=-1,1; // = Left sided finite differences for -d/dx and -d/dy
643 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD); // with hom. Dirichlet b. c.
644
645 u_bar=out;
646
647 double m=0;
648 for (int y=0;y<data.shape(1);y++){
649 for (int x=0;x<data.shape(0);x++){
650 m=std::max(m,alpha(x,y));
651 m=std::max(m,beta (x,y));
652 m=std::max(m,gamma(x,y));
653 }
654 }
655 m=std::max(m,1.);
656 double tau=.1/m;//std::sqrt(8)*0.06;
657 double sigma=.1;//m;/std::sqrt(8)/0.06;
658
659 //std::cout<<"tau= "<<tau<<std::endl;
660
661 for (int i=0;i<steps;i++){
662
663 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroX(temp1);
664 vx+=(sigma*temp1);
665 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroY(temp1);
666 vy+=(sigma*temp1);
667
668
669 // update wx
670 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroX(temp1);
671 temp1*=xedges;
672 separableConvolveX(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
673 wx-=sigma*temp2;//(-Lx'*(xedges.*(Lx*u)));
674
675 //update wy
676 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroY(temp1);
677 temp1*=yedges;
678 separableConvolveY(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
679 wy-=sigma*temp2;//(-Ly'*(yedges.*(Ly*u)));
680
681
682 //update wz
683 #if (VIGRA_MIXED_2ND_DERIVATIVES)
684 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroY(temp1);
685 temp1*=yedges;
686 separableConvolveX(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
687 wz-=sigma*temp2;//-Lx'*(yedges.*(Ly*u))
688
689 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroX(temp1);
690 temp1*=xedges;
691 separableConvolveY(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
692 wz-=sigma*temp2;//-Ly'*(xedges.*(Lx*u)));
693
694 #endif
695
696
697 //project to constraint sets
698 for (int y=0;y<data.shape(1);y++){
699 for (int x=0;x<data.shape(0);x++){
700 double e1,e2,skp1,skp2;
701
702 //project v
703 e1=std::cos(phi(x,y));
704 e2=std::sin(phi(x,y));
705 skp1=vx(x,y)*e1+vy(x,y)*e2;
706 skp2=vx(x,y)*(-e2)+vy(x,y)*e1;
707 vigra::detail::projectEllipse2D (skp1,skp2,alpha(x,y),beta(x,y),0.001,100);
708 vx(x,y)=skp1*e1-skp2*e2;
709 vy(x,y)=skp1*e2+skp2*e1;
710
711 //project w
712 double l=sqrt(wx(x,y)*wx(x,y)+wy(x,y)*wy(x,y)+wz(x,y)*wz(x,y));
713 if (l>gamma(x,y)){
714 wx(x,y)=gamma(x,y)*wx(x,y)/l;
715 wy(x,y)=gamma(x,y)*wy(x,y)/l;
716 #if (VIGRA_MIXED_2ND_DERIVATIVES)
717 wz(x,y)=gamma(x,y)*wz(x,y)/l;
718 #endif
719 }
720 }
721 }
722
723 separableConvolveX(srcImageRange(vx),destImage(temp1),kernel1d(LTx));
724 separableConvolveY(srcImageRange(vy),destImage(temp2),kernel1d(LTx));
725
726 u_bar=out;
727 out-=tau*(weight*(out-data)+temp1+temp2);
728
729
730 // update wx
731 separableConvolveX(srcImageRange(wx),destImage(temp1),kernel1d(Lx));setZeroX(temp1);
732 temp1*=xedges;
733 separableConvolveX(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
734 out+=tau*temp2; // (-1)^2
735
736
737 //update wy
738 separableConvolveY(srcImageRange(wy),destImage(temp1),kernel1d(Lx));setZeroY(temp1);
739 temp1*=yedges;
740 separableConvolveY(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
741 out+=tau*temp2;
742
743 //update wz
744 #if (VIGRA_MIXED_2ND_DERIVATIVES)
745
746 separableConvolveY(srcImageRange(wz),destImage(temp1),kernel1d(Lx));setZeroY(temp1);
747 temp1*=yedges;
748 separableConvolveX(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
749 out+=tau*temp2;
750
751 separableConvolveX(srcImageRange(wz),destImage(temp1),kernel1d(Lx));setZeroX(temp1);
752 temp1*=xedges;
753 separableConvolveY(srcImageRange(temp1),destImage(temp2),kernel1d(LTx));
754 out+=tau*temp2;
755
756 #endif
757
758 u_bar=2*out-u_bar; //cf. Chambolle/Pock and Popov's algorithm
759
760 }
761}
762
763//@}
764} // closing namespace vigra
765
766#endif // VIGRA_TV_FILTER_HXX
Class for a single RGB value.
Definition rgbvalue.hxx:128
image import and export functions
void totalVariationFilter(...)
Performs standard Total Variation Regularization.
void structureTensor(...)
Calculate the Structure Tensor for each pixel of and image, using Gaussian (derivative) filters.
void separableConvolveX(...)
Performs a 1 dimensional convolution in x direction.
double gamma(double x)
The gamma function.
Definition mathutil.hxx:1587
void separableConvolveY(...)
Performs a 1 dimensional convolution in y direction.
void anisotropicTotalVariationFilter(...)
Performs Anisotropic Total Variation Regularization.
void secondOrderTotalVariationFilter(...)
Performs Anisotropic Total Variation Regularization.
FixedPoint16< IntBits, OverflowHandling > hypot(FixedPoint16< IntBits, OverflowHandling > v1, FixedPoint16< IntBits, OverflowHandling > v2)
Length of hypotenuse.
Definition fixedpoint.hxx:1640
void getAnisotropy(...)
Sets up directional data for anisotropic regularization.

© 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.11.2