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;
141template <
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){
211template <
class str
ide1,
class str
ide2,
class str
ide3>
212void 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){
340template <
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;
468template <
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));
622template <
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
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition separableconvolution.hxx:2253
void setBorderTreatment(BorderTreatmentMode new_mode)
Definition separableconvolution.hxx:2170
Kernel1D & initExplicitly(int left, int right)
Definition separableconvolution.hxx:2100
Base class for, and view to, vigra::MultiArray.
Definition multi_fwd.hxx:127
const difference_type & shape() const
Definition multi_array.hxx:1650
Main MultiArray class containing the memory management.
Definition multi_fwd.hxx:131
Class for fixed size vectors.
Definition tinyvector.hxx:1008
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
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition fixedpoint.hxx:616
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.