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
139 template <
class str
ide1,
class str
ide2>
140 void totalVariationFilter(MultiArrayView<2,double,stride1> data,MultiArrayView<2,double,stride2> out,
double alpha,
int steps,
double eps=0){
142 using namespace multi_math;
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;
147 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT);
148 LTx.initExplicitly(-1,0)=-1,1;
149 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD);
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;
157 for (
int i=0;i<steps;i++){
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));
178 out-=tau*(out-data+alpha*(temp1+temp2));
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));
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;
201 if (f_primal>0 && (f_primal-f_dual)/f_primal<eps){
208 template <
class str
ide1,
class str
ide2,
class str
ide3>
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){
211 using namespace multi_math;
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;
216 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT);
217 LTx.initExplicitly(-1,0)=-1,1;
218 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD);
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;
226 for (
int i=0;i<steps;i++){
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));
246 out-=tau*(weight*(out-data)+alpha*(temp1+temp2));
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));
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;
269 if (f_primal>0 && (f_primal-f_dual)/f_primal<eps){
334 doxygen_overloaded_function(template <...>
void getAnisotropy)
336 template <
class str
ide1,
class str
ide2,
class str
ide3,
class str
ide4>
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){
341 using namespace multi_math;
343 MultiArray<2,double> smooth(data.shape()),tmp(data.shape());
351 MultiArray<2,double> stxx(data.shape()),stxy(data.shape()),styy(data.shape());
356 gauss.initGaussian(rho_par);
366 for (
int y=0;y<data.shape(1);y++){
367 for (
int x=0;x<data.shape(0);x++){
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);
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;
464 template <
class str
ide1,
class str
ide2,
class str
ide3,
class str
ide4,
class str
ide5,
class str
ide6>
466 MultiArrayView<2,double,stride3> phi,MultiArrayView<2,double,stride4> alpha,
467 MultiArrayView<2,double,stride5> beta,MultiArrayView<2,double,stride6> out,
470 using namespace multi_math;
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());
475 Kernel1D<double> Lx,LTx;
476 Lx.initExplicitly(0,1)=1,-1;
477 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT);
478 LTx.initExplicitly(-1,0)=-1,1;
479 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD);
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));
495 for (
int i=0;i<steps;i++){
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;
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);
512 vx(x,y)=skp1*e1-skp2*e2;
513 vy(x,y)=skp1*e2+skp2*e1;
520 out-=tau*(weight*(out-data)+(temp1+temp2));
617 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>
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,
626 using namespace multi_math;
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());
633 Kernel1D<double> Lx,LTx;
634 Lx.initExplicitly(0,1)=1,-1;
635 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT);
636 LTx.initExplicitly(-1,0)=-1,1;
637 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD);
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));
652 for (
int i=0;i<steps;i++){
674 #if (VIGRA_MIXED_2ND_DERIVATIVES)
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;
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;
703 double l=
sqrt(wx(x,y)*wx(x,y)+wy(x,y)*wy(x,y)+wz(x,y)*wz(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;
718 out-=tau*(weight*(out-data)+temp1+temp2);
735 #if (VIGRA_MIXED_2ND_DERIVATIVES)
757 #endif // VIGRA_TV_FILTER_HXX