37 #ifndef VIGRA_TV_FILTER_HXX
38 #define VIGRA_TV_FILTER_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"
52 #ifndef VIGRA_MIXED_2ND_DERIVATIVES
53 #define VIGRA_MIXED_2ND_DERIVATIVES 1
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;
141 template <
class str
ide1,
class str
ide2>
144 using namespace multi_math;
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;
160 for (
int i=0;i<steps;i++){
162 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroX(temp1);
164 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroY(temp1);
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));
181 out-=tau*(out-data+alpha*(temp1+temp2));
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));
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;
204 if (f_primal>0 && (f_primal-f_dual)/f_primal<eps){
211 template <
class str
ide1,
class str
ide2,
class str
ide3>
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){
214 using namespace multi_math;
215 int width=data.shape(0),height=data.shape(1);
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;
220 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT);
221 LTx.initExplicitly(0,1)=-1,1;
222 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD);
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;
230 for (
int i=0;i<steps;i++){
231 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroX(temp1);
233 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroY(temp1);
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));
250 out-=tau*(weight*(out-data)+alpha*(temp1+temp2));
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));
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;
273 if (f_primal>0 && (f_primal-f_dual)/f_primal<eps){
340 template <
class str
ide1,
class str
ide2,
class str
ide3,
class str
ide4>
343 double alpha_par,
double beta_par,
double sigma_par,
double rho_par,
double K_par){
345 using namespace multi_math;
370 for (
int y=0;y<data.
shape(1);y++){
371 for (
int x=0;x<data.
shape(0);x++){
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);
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;
468 template <
class str
ide1,
class str
ide2,
class str
ide3,
class str
ide4,
class str
ide5,
class str
ide6>
474 using namespace multi_math;
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));
496 double tau=.9/m/std::sqrt(8.)*0.06;
497 double sigma=.9/m/std::sqrt(8.)/0.06;
500 for (
int i=0;i<steps;i++){
501 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroX(temp1);
503 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroY(temp1);
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;
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);
517 vx(x,y)=skp1*e1-skp2*e2;
518 vy(x,y)=skp1*e2+skp2*e1;
525 out-=tau*(weight*(out-data)+(temp1+temp2));
622 template <
class str
ide1,
class str
ide2,
class str
ide3,
class str
ide4,
class str
ide5,
class str
ide6,
class str
ide7,
class str
ide8,
class str
ide9>
631 using namespace multi_math;
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));
661 for (
int i=0;i<steps;i++){
663 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroX(temp1);
665 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroY(temp1);
670 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroX(temp1);
676 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroY(temp1);
683 #if (VIGRA_MIXED_2ND_DERIVATIVES)
684 separableConvolveY(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroY(temp1);
689 separableConvolveX(srcImageRange(u_bar),destImage(temp1),kernel1d(Lx));setZeroX(temp1);
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;
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;
712 double l=
sqrt(wx(x,y)*wx(x,y)+wy(x,y)*wy(x,y)+wz(x,y)*wz(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;
727 out-=tau*(weight*(out-data)+temp1+temp2);
744 #if (VIGRA_MIXED_2ND_DERIVATIVES)
Generic 1 dimensional convolution kernel.
Definition: separableconvolution.hxx:1367
Kernel1D & initExplicitly(int left, int right)
Definition: separableconvolution.hxx:2096
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2249
void setBorderTreatment(BorderTreatmentMode new_mode)
Definition: separableconvolution.hxx:2166
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:705
const difference_type & shape() const
Definition: multi_array.hxx:1648
image import and export functions
void structureTensor(...)
Calculate the Structure Tensor for each pixel of and image, using Gaussian (derivative) filters.
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
void separableConvolveX(...)
Performs a 1 dimensional convolution in x direction.
double gamma(double x)
The gamma function.
Definition: mathutil.hxx:1587
void totalVariationFilter(...)
Performs standard Total Variation Regularization.
void secondOrderTotalVariationFilter(...)
Performs Anisotropic Total Variation Regularization.
void getAnisotropy(...)
Sets up directional data for anisotropic regularization.
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
void anisotropicTotalVariationFilter(...)
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 separableConvolveY(...)
Performs a 1 dimensional convolution in y direction.