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;
145 int width=data.shape(0),height=data.shape(1);
149 Lx.initExplicitly(-1,0)=1,-1;
150 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT);
151 LTx.initExplicitly(0,1)=-1,1;
152 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD);
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;
168 for (
int y=0;y<data.shape(1);y++){
169 for (
int x=0;x<data.shape(0);x++){
191 for (
int y=0;y<data.shape(1);y++){
192 for (
int x=0;x<data.shape(0);x++){
198 for (
int y=0;y<data.shape(1);y++){
199 for (
int x=0;x<data.shape(0);x++){
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>
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++){
377 vigra::symmetricEigensystemNoniterative(
matrix,
ew,
ev);
379 phi(x,y)=std::atan2(
ev(1,0),
ev(0,0));
468template <
class str
ide1,
class str
ide2,
class str
ide3,
class str
ide4,
class str
ide5,
class str
ide6>
474 using namespace multi_math;
475 int width=data.shape(0),height=data.shape(1);
481 Lx.initExplicitly(-1,0)=1,-1;
482 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT);
483 LTx.initExplicitly(0,1)=-1,1;
484 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD);
489 for (
int y=0;y<data.shape(1);y++){
490 for (
int x=0;x<data.shape(0);x++){
496 double tau=.9/
m/std::sqrt(8.)*0.06;
497 double sigma=.9/
m/std::sqrt(8.)/0.06;
507 for (
int y=0;y<data.shape(1);y++){
508 for (
int x=0;x<data.shape(0);x++){
511 e1=std::cos(
phi(x,y));
512 e2=std::sin(
phi(x,y));
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;
632 int width=data.shape(0),height=data.shape(1);
640 Lx.initExplicitly(-1,0)=1,-1;
641 Lx.setBorderTreatment(BORDER_TREATMENT_REFLECT);
642 LTx.initExplicitly(0,1)=-1,1;
643 LTx.setBorderTreatment(BORDER_TREATMENT_ZEROPAD);
648 for (
int y=0;y<data.shape(1);y++){
649 for (
int x=0;x<data.shape(0);x++){
683 #if (VIGRA_MIXED_2ND_DERIVATIVES)
698 for (
int y=0;y<data.shape(1);y++){
699 for (
int x=0;x<data.shape(0);x++){
703 e1=std::cos(
phi(x,y));
704 e2=std::sin(
phi(x,y));
712 double l=sqrt(
wx(x,y)*
wx(x,y)+
wy(x,y)*
wy(x,y)+
wz(x,y)*
wz(x,y));
716 #if (VIGRA_MIXED_2ND_DERIVATIVES)
744 #if (VIGRA_MIXED_2ND_DERIVATIVES)
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.