Total Variation Deconvolution using Split Bregman
tvregopt.h
Go to the documentation of this file.
1 
14 #ifndef _TVREGOPT_H_
15 #define _TVREGOPT_H_
16 
17 #if defined(TVREG_NONGAUSSIAN) || (defined(TVREG_INPAINT) && defined(TVREG_DECONV))
18 #define TVREG_USEZ
19 #endif
20 
21 #ifdef TVREG_DECONV
22 #include <fftw3.h>
23 #endif
24 #include "tvreg.h"
25 
27 #define ALGSTRING_SIZE 128
28 
41 #define _TVREG_CONCAT(A,B) A ## B
42 
43 #ifdef NUM_SINGLE
44 #define FFT(S) _TVREG_CONCAT(fftwf_,S)
45 #else
46 #define FFT(S) _TVREG_CONCAT(fftw_,S)
47 #endif
48 
49 
50 /* Internal type definitions */
51 
53 typedef struct
54 {
55  num x;
56  num y;
57 } numvec2;
58 
60 typedef num numcomplex[2];
61 
63 typedef enum {
64  NOISEMODEL_L2,
65  NOISEMODEL_L1,
66  NOISEMODEL_POISSON
67 } noisemodel;
68 
71 {
72  num Lambda;
73  const num *VaryingLambda;
74  int LambdaWidth;
75  int LambdaHeight;
76  const num *Kernel;
77  int KernelWidth;
78  int KernelHeight;
79  num Tol;
80  num Gamma1;
81  num Gamma2;
82  int MaxIter;
83  noisemodel NoiseModel;
84  int (*PlotFun)(int, int, num, const num*, int, int, int, void*);
85  void *PlotParam;
86  char *AlgString;
87 };
88 
95 typedef struct tag_tvregsolver
96 {
97  num *u;
98  const num *f;
101  num *Ku;
103  num fNorm;
104  num Alpha;
105  int Width;
106  int Height;
107  int PadWidth;
108  int PadHeight;
111  int UseZ;
113 #ifdef TVREG_USEZ
114  num *z;
115  num *ztilde;
116 #endif
117 
118 #ifdef TVREG_DECONV
119  num *A, *B;
120  num *ATrans, *BTrans;
121  num *DenomTrans;
122  num *KernelTrans;
123  FFT(plan) TransformA;
124  FFT(plan) TransformB;
125  FFT(plan) InvTransformA;
126  FFT(plan) InvTransformB;
127 #endif
128 } tvregsolver;
129 
130 typedef num (*usolver)(tvregsolver*);
131 typedef void (*zsolver)(tvregsolver*);
132 
137  TvRestoreSimplePlot, NULL, NULL};
138 
139 static int TvRestoreChooseAlgorithm(int *UseZ, int *DeconvFlag, int *DctFlag,
140  usolver *USolveFun, zsolver *ZSolveFun, const tvregopt *Opt);
141 
142 
143 /* If GNU C language extensions are available, apply the "unused" attribute
144  to avoid warnings. TvRestoreSimplePlot is a plotting callback function
145  for TvRestore, so the unused arguments are indeed required. */
146 int TvRestoreSimplePlot(int State, int Iter, num Delta,
147  ATTRIBUTE_UNUSED const num *u,
148  ATTRIBUTE_UNUSED int Width,
149  ATTRIBUTE_UNUSED int Height,
150  ATTRIBUTE_UNUSED int NumChannels,
151  ATTRIBUTE_UNUSED void *Param)
152 {
153  switch(State)
154  {
155  case 0: /* TvRestore is running */
156  /* We print to stderr so that messages are displayed on the console
157  immediately, during the TvRestore computation. If we use stdout,
158  messages might be buffered and not displayed until after TvRestore
159  completes, which would defeat the point of having this real-time
160  plot callback. */
161  fprintf(stderr, " Iteration %4d Delta %7.4f\r", Iter, Delta);
162  break;
163  case 1: /* Converged successfully */
164  fprintf(stderr, "Converged in %d iterations. \n", Iter);
165  break;
166  case 2: /* Maximum iterations exceeded */
167  fprintf(stderr, "Maximum number of iterations exceeded.\n");
168  break;
169  }
170  return 1;
171 }
172 
173 
183 {
184  tvregopt *Opt;
185 
186  if((Opt = (tvregopt *)Malloc(sizeof(tvregopt))))
187  *Opt = TvRegDefaultOpt;
188 
189  if(!(Opt->AlgString = (char *)Malloc(sizeof(char)*ALGSTRING_SIZE)))
190  {
191  Free(Opt);
192  return NULL;
193  }
194 
195  return Opt;
196 }
197 
198 
204 {
205  if(Opt)
206  {
207  if(Opt->AlgString)
208  Free(Opt->AlgString);
209  Free(Opt);
210  }
211 }
212 
213 
219 void TvRegSetLambda(tvregopt *Opt, num Lambda)
220 {
221  if(Opt)
222  Opt->Lambda = Lambda;
223 }
224 
225 
248  const num *VaryingLambda, int LambdaWidth, int LambdaHeight)
249 {
250  if(Opt)
251  {
252  Opt->VaryingLambda = VaryingLambda;
253  Opt->LambdaWidth = LambdaWidth;
254  Opt->LambdaHeight = LambdaHeight;
255  }
256 }
257 
258 
271  const num *Kernel, int KernelWidth, int KernelHeight)
272 {
273  if(Opt)
274  {
275  Opt->Kernel = Kernel;
276  Opt->KernelWidth = KernelWidth;
277  Opt->KernelHeight = KernelHeight;
278  }
279 }
280 
281 
287 void TvRegSetTol(tvregopt *Opt, num Tol)
288 {
289  if(Opt)
290  Opt->Tol = Tol;
291 }
292 
293 
299 void TvRegSetGamma1(tvregopt *Opt, num Gamma1)
300 {
301  if(Opt)
302  Opt->Gamma1 = Gamma1;
303 }
304 
305 
311 void TvRegSetGamma2(tvregopt *Opt, num Gamma2)
312 {
313  if(Opt)
314  Opt->Gamma2 = Gamma2;
315 }
316 
317 
323 void TvRegSetMaxIter(tvregopt *Opt, int MaxIter)
324 {
325  if(Opt)
326  Opt->MaxIter = MaxIter;
327 }
328 
329 
346 int TvRegSetNoiseModel(tvregopt *Opt, const char *NoiseModel)
347 {
348  if(!Opt)
349  return 0;
350 
351  if(!NoiseModel || !strcmp(NoiseModel, "L2") || !strcmp(NoiseModel, "l2")
352  || !strcmp(NoiseModel, "Gaussian") || !strcmp(NoiseModel, "gaussian"))
353  Opt->NoiseModel = NOISEMODEL_L2;
354  else if(!strcmp(NoiseModel, "L1") || !strcmp(NoiseModel, "l1")
355  || !strcmp(NoiseModel, "Laplace") || !strcmp(NoiseModel, "laplace")
356  || !strcmp(NoiseModel, "Laplacian") || !strcmp(NoiseModel, "laplacian"))
357  Opt->NoiseModel = NOISEMODEL_L1;
358  else if(!strcmp(NoiseModel, "Poisson") || !strcmp(NoiseModel, "poisson"))
359  Opt->NoiseModel = NOISEMODEL_POISSON;
360  else
361  return 0;
362 
363  return 1;
364 }
365 
366 
405  int (*PlotFun)(int, int, num, const num*, int, int, int, void*),
406  void *PlotParam)
407 {
408  if(Opt)
409  {
410  Opt->PlotFun = PlotFun;
411  Opt->PlotParam = PlotParam;
412  }
413 }
414 
415 
420 void TvRegPrintOpt(const tvregopt *Opt)
421 {
422  if(!Opt)
423  Opt = &TvRegDefaultOpt;
424 
425  printf("lambda : ");
426 
427  if(!Opt->VaryingLambda)
428  printf("%g\n", Opt->Lambda);
429  else
430  printf("[%d x %d]\n",
431  Opt->LambdaWidth, Opt->LambdaHeight);
432 
433  printf("K : ");
434 
435  if(!Opt->Kernel)
436  printf("(identity)\n");
437  else
438  printf("[%d x %d]\n", Opt->KernelWidth, Opt->KernelHeight);
439 
440  printf("tol : %g\n", (double)Opt->Tol);
441  printf("max iter : %d\n", Opt->MaxIter);
442  printf("gamma1 : %g\n", (double)Opt->Gamma1);
443  printf("gamma2 : %g\n", (double)Opt->Gamma2);
444  printf("noise : ");
445 
446  switch(Opt->NoiseModel)
447  {
448  case NOISEMODEL_L2:
449  printf("L2\n");
450  break;
451  case NOISEMODEL_L1:
452  printf("L1\n");
453  break;
454  case NOISEMODEL_POISSON:
455  printf("Poisson\n");
456  break;
457  default:
458  printf("(invalid)\n");
459  break;
460  }
461 
462  printf("plotting : ");
463 
464  if(Opt->PlotFun == TvRestoreSimplePlot)
465  printf("default\n");
466  else if(!Opt->PlotFun)
467  printf("none\n");
468  else
469  printf("custom\n");
470 
471  printf("algorithm : %s\n", TvRegGetAlgorithm(Opt));
472 }
473 
474 
484 const char *TvRegGetAlgorithm(const tvregopt *Opt)
485 {
486  static const char *DefaultAlgorithm =
487  (char *)"split Bregman (d = grad u) Gauss-Seidel u-solver";
488  static const char *Invalid = (char *)"(invalid)";
489  usolver USolveFun;
490  zsolver ZSolveFun;
491  int UseZ, DeconvFlag, DctFlag;
492 
493  if(!Opt)
494  return DefaultAlgorithm;
495 
496  if(!TvRestoreChooseAlgorithm(&UseZ, &DeconvFlag,
497  &DctFlag, &USolveFun, &ZSolveFun, Opt))
498  return Invalid;
499 
500  sprintf(Opt->AlgString, "split Bregman (%s) %s u-solver",
501  (UseZ) ?
502  "d = grad u, z = Ku" :
503  "d = grad u",
504  (!DeconvFlag) ?
505  "Gauss-Seidel" :
506  ((DctFlag) ?
507  "DCT" :
508  "Fourier"));
509  return Opt->AlgString;
510 }
511 
512 #endif /* _TVREGOPT_H_ */