Toast++  1.0.2 (r.539)
Forward and inverse modelling in optical tomography
regul.h
1 // -*-C++-*-
2 #ifndef __ST_REGUL_H
3 #define __ST_REGUL_H
4 
5 #include "stoastlib.h"
6 
7 class Raster;
8 class ParamParser;
9 
10 #ifndef SQR
11 #define SQR(_X) ((_X)*(_X))
12 #endif
13 #ifndef CUBE
14 #define CUBE(_X) ((_X)*(_X)*(_X))
15 #endif
16 typedef double (* PSIREG)(const double, const int np, double *p); //pointer to generic function
17 
18 typedef enum {
19  PRIOR_NONE,
20  PRIOR_LAPLACIAN,
21  PRIOR_DIAG,
22  PRIOR_TK1,
23  PRIOR_TV,
24  PRIOR_HUBER,
25  PRIOR_PM,
26  PRIOR_QPM,
27  PRIOR_TK1SIGMA,
28  PRIOR_TVSIGMA,
29  PRIOR_HUBERSIGMA,
30  PRIOR_PMSIGMA,
31  PRIOR_QPMSIGMA,
32  PRIOR_TUKEYSIGMA,
33  PRIOR_GGMRF,
34  PRIOR_GENERIC,
35  PRIOR_GENERICSIGMA,
36  PRIOR_GENERICSCALESPACE,
37  PRIOR_MRF
38 } PRIOR;
39 
40 
41 // ====================== CLASSES =====================================
42 class STOASTLIB Regularisation {
43 public:
44  Regularisation (const Raster *_raster = 0, double _tau = 1,
45  const RVector *_x0 = 0);
46 
47  virtual ~Regularisation ();
48 
49  static Regularisation *Create (ParamParser *pp, const RVector *_x0,
50  const Raster *_raster, const RVector *_xs = 0);
51  // create a regularisation instance from a parameter file
52 
53  virtual void ReadParams (ParamParser *pp);
54  virtual void WriteParams (ParamParser *pp);
55 
56  virtual const char *GetName() const = 0;
57  // return a name to identify the regulariser
58 
59  inline double GetTau() const
60  { return tau; }
61 
62  inline int GetNParam() const
63  { return nset; }
64 
65  virtual double GetValue (const RVector &x) const = 0;
66  // value of prior psi for given coefficient vector x
67 
68  virtual RVector GetGradient (const RVector &x) const = 0;
69  // gradient dpsi/dx_i for given coefficient vector x
70 
71  virtual RVector GetKappa (const RVector &x) const = 0;
72  // diffusivity function for given coefficient vector x
73 
74  virtual void SetHess1 (RCompRowMatrix &Hess1, const RVector &x, const int p) =0;
75  // set the first order Hessian for input vector x
76  virtual void SetHess1FromKappa (RCompRowMatrix &Hess, const RVector &kap) =0;
77  virtual RVector GetHess1f (const RVector &x, const RVector &f) const = 0;
78  // Effect of 1st order Hessian L on vector f at given coefficient vector x
79 
80  virtual void SetFullHess (RCompRowMatrix &Hess, const RVector &x, const int p) =0;
81  // set the total Hessian H for input vector x
82 
83  virtual RVector GetFullHessf (const RVector &x, const RVector &f) const = 0;
84  // Effect of total Hessian H on vector f at given coefficient vector x
85 
86  // virtual void SetHessian (RCompRowMatrix &Hess, const RVector &x) =0;
87 
88  virtual int GetHessianRow (const RVector &x, int i, idxtype *colidx,
89  double *val) const = 0;
90  // returns single row i of the 2nd derivative matrix
91  // d^2 psi / dx_i dx_j in sparse vector format. Return value is
92  // number of nonzeros. column index list is returned in colidx,
93  // matrix values are returned in val. Both colidx and val must
94  // be allocated by the caller with sufficient size
95 
96  virtual RVector GetHessianDiag (const RVector &x) const = 0;
97  // returns the diagonal of the 2nd derivative matrix
98 
99  virtual PRIOR Type() const = 0;
100 
101 protected:
102  void CreateHessStruct1param (RCompRowMatrix &Hess) ;
103  // set up the Hessian structure for a single parameter.
104  // multiple Hessians can be setup by concatenation
105 
106  int nset, pdim;
107  // void CreateNonLinearHessian (RCompRowMatrix &Hess, const RVector &x);
108  double tau; // regularisation hyperparameter
109  PRIOR ptype;
110  const Raster *raster;
111  const RVector *x0;
112  int slen, glen, dim;
113  idxtype *rowptr, *colidx;
114  int nzero;
115  RCompRowMatrix Hs1p; // will hold sparse matrix structure for 1 Hessian
116 };
117 #ifdef REG_1ST_ORDER
118 // =========================================================================
119 class Regularisation_1stOrder : public Regularisation {
120  // base class for Regularisation using gradient information
121  Regularisation_1stOrder (const Raster *_raster = 0,const RVector *_x0 = 0);
122  virtual ~Regularisation_1stOrder () {};
123 
124 protected:
125  void CreateHessStruct1param (RCompRowMatrix &Hess);
126  // set up the Hessian structure for a single parameter.
127  // multiple Hessians can be setup by concatenation
128 
129 };
130 
131 // =========================================================================
132 
133 class STOASTLIB GenericSigma : public Regularisation_1stOrder {
134 #else
135 class STOASTLIB GenericSigma : public Regularisation {
136 #endif
137 public:
138  // GenericSigma (const Raster *_raster, const RVector *_x0, const int _npreg, double * _preg, double _tau, double _sdx = 0, double _sdf = 0, const RVector *_kref=0);
139  GenericSigma (const Raster *_raster, const RVector *_x0, const int _npreg,
140  double * _preg, double _tau = 1e-2, double _sdx = 0, double _sdf = 0,
141  const void *_kref=0, bool _KRefTensFlag = false);
142 
143  GenericSigma (const Raster *_raster, const RVector *_x0, int _npreg,
144  double *_preg, double _tau, const RVector &_refimg, double _sdr,
145  double _sdx = 0, double _sdf = 0, double _fT = 0,
146  bool _KRefTensFlag = false);
147 
148  virtual ~GenericSigma ();
149 
150  const char *GetName() const
151  { return "GenericSigma"; }
152 
153  void ReadParams (ParamParser *pp);
154  void WriteParams (ParamParser *pp);
155 
156  RVector *MakeKref (const RVector &gimgref, double sdr, double fT);
157  // Return scalar field for isotropic diffusivity distribution given
158  // image gimgref. sdr: value for reference sigma; fT: fraction of max to
159  // use for Perona-Malik threshold [0-1]
160 
161  RDenseMatrix *MakeKrefTens (const RVector &gimgref, double sdr, double fT);
162  // Return tensor field for anisotropic diffusivity distribution given
163  // image gimgref. sdr: value for reference sigma; fT: fraction of max to
164  // use for Perona-Malik threshold [0-1]
165 
166  PRIOR Type() const {return PRIOR_GENERICSIGMA; }
167 
168  // void SetHessian (RCompRowMatrix &Hess, const RVector &x) {
169  // CreateNonLinearHessian(raster, x, Hess);
170  // };
171 
172  double GetValue (const RVector &x) const;
173  RVector GetGradient (const RVector &x) const {
174  // MS: modified 100201
175  // the second parameter of GetHess1f doesn't get x0 subtracted
176  // inside GetHess1f, which can lead to a nonzero gradient even
177  // if x = x0
178 
179  //return GetHess1f(x,x);
180  return GetHess1f (x, x - *x0);
181  }
182  RVector GetKappa (const RVector &x) const;
183  void SetHess1 (RCompRowMatrix &Hess1, const RVector &x, const int p);
184  void SetHess1FromKappa (RCompRowMatrix &Hess, const RVector &kap) ;
185  RVector GetHess1f (const RVector &x, const RVector &f) const;
186  void SetFullHess (RCompRowMatrix &Hess, const RVector &x, const int p);
187  RVector GetFullHessf (const RVector &x, const RVector &f) const;
188  // void SetHessian (RCompRowMatrix &Hess, const RVector &x) ;
189 
190  int GetHessianRow (const RVector &x, int i, idxtype *colidx, double *val)const;
191 
192  RVector GetHessianDiag (const RVector &x) const ;
193 
194  inline RDenseMatrix tinterp (const RDenseMatrix *mat, int i, int j) const;
195 
196 protected:
197  double sdr, fT;
198  const double sdx, sdf;
199  virtual double func(const double, const int np, double *p) const = 0;
200  virtual double kfunc(const double, const int np, double *p)const = 0;
201  virtual double d2func(const double, const int np, double *p)const = 0;
202  void CreateDerivativeOperators ();
203  const int npreg;
204  char *krefimg_name;
205  const RVector *kref;
206  const RDenseMatrix *kreftens;
207  const bool KRefTensFlag;
208  bool KRefIsLocal;
209  double *preg;
210  ICompRowMatrix NG; // used only if tensor flag set
211  RCompRowMatrix *Db, *Df; //will hold derivative oeprators
212 };
213 
214 // ===================== Some Utility functions ==========================
215 inline double sinterp(const RVector & v, const int i, const int j)
216 {
217  return 0.5*(v[i] + v[j]);
218 }
219 #ifdef LINTERP
220 inline RDenseMatrix GenericSigma::tinterp (const RDenseMatrix *mat,
221  int i, int j) const
222 {
223  RDenseMatrix tmp(mat[i]);
224  for(int r = 0; r < tmp.nRows(); r++)
225  for(int c = 0; c < tmp.nCols(); c++)
226  tmp(r,c) = 0.5*(mat[i](r,c) + mat[j](r,c));
227  return tmp;
228 }
229 #else
230 #ifdef MINDET
231 //return matrix with minimum determinant!
232 inline RDenseMatrix GenericSigma::tinterp (const RDenseMatrix *mat,
233  int i, int j) const
234 {
235  double d1 = det(mat[i]);
236  double d2 = det(mat[j]);
237  return (d1 < d2 ? mat[i] : mat[j]);
238 }
239 #else
240 inline RDenseMatrix GenericSigma::tinterp (const RDenseMatrix *mat,
241  int i, int j) const
242 // calculate tensor from interpolated gradients
243 {
244 #ifdef TINTERP2D // old version
245  static const int dim = 2; // only 2D for now.
246 
247  double fx = 0.5*( mat[i](0,0) + mat[j](0,0));
248  double fy = 0.5*( mat[i](1,1) + mat[j](1,1));
249  double fxsq = SQR(fx);
250  double fysq = SQR(fy);
251  double gsq = fxsq + fysq;
252  double T = 0.5*( mat[i](0,1) + mat[j](0,1));
253  double Tsq = SQR(T);
254  RDenseMatrix tmp(dim,dim);
255  tmp.Identity(dim);
256  // double gam = 1- exp( -gsq/Tsq);
257  double gam = 1- exp( -sqrt(gsq)/T);
258  tmp(0,0) -= (fx == 0.0 ? 0 : gam*fxsq/gsq); // xx term
259  tmp(1,1) -= (fy == 0.0 ? 0 : gam*fysq/gsq); // yy term
260  tmp(0,1) -= (gsq == 0.0 ? 0 : gam*fx*fy/gsq); // xy term
261  tmp(1,0) = tmp(0,1); // xy term
262  return tmp;
263 #else // new, but untested
264  int m, n;
265  double gsq = 0.0;
266  RVector f(dim);
267  for (m = 0; m < dim; m++) {
268  f[m] = 0.5*(mat[i](m,m) + mat[j](m,m));
269  gsq += SQR(f[m]);
270  }
271  double T = 0.5*(mat[i](0,1) + mat[j](0,1));
272  double Tsq = SQR(T);
273  RDenseMatrix tmp(dim,dim);
274  tmp.Identity (dim);
275  double gam = 1.0-exp(-sqrt(gsq)/T);
276  for (m = 0; m < dim; m++)
277  for (n = 0; n < dim; n++)
278  tmp(m,n) -= (gsq == 0.0 ? 0 : gam*f[m]*f[n]/gsq);
279  return tmp;
280 #endif
281 }
282 #endif
283 #endif
284 
285 // ========================"Classic" TV========================================
286 // Scale Space literature refers to this function as "Charbonnier"
287 //
288 class STOASTLIB TVSigma: public GenericSigma {
289 public:
290  TVSigma (double _tau, double _beta, double _sd, const RVector *_x0,
291  const Raster *_raster, bool _SmoothKappaOnly = true,
292  const void *_kref = 0, bool _KRefTensFlag = false)
293  : GenericSigma (_raster, _x0, 1, (new double [1]), _tau, _sd,
294  (_SmoothKappaOnly ? 0 : _sd), _kref, _KRefTensFlag)
295  {preg[0] = _beta;}
296 
297  TVSigma (double _tau, double _beta, double _sd, const RVector *_x0,
298  const Raster *_raster, const RVector &_refimg, double _sdr,
299  double _fT, bool _SmoothKappaOnly = true, bool _KRefTensFlag = false)
300  : GenericSigma (_raster, _x0, 1, (new double [1]), _tau, _refimg,
301  _sdr, _sd, (_SmoothKappaOnly ? 0 : _sd), _fT, _KRefTensFlag)
302  {preg[0] = _beta;};
303 
304  ~TVSigma ()
305  {}
306 
307  const char *GetName() const
308  { return "TV_SIGMA"; }
309 
310  void ReadParams (ParamParser *pp);
311  void WriteParams (ParamParser *pp);
312 
313 protected:
314  inline double func(const double t, const int np, double *p) const{
315  double betasq = SQR(p[0]);
316  return p[0]*sqrt(SQR(t) + betasq) - betasq ;
317  }
318 
319  inline double kfunc(const double t, const int np, double *p) const{
320  double betasq = SQR(p[0]);
321  return p[0]/sqrt(SQR(t) + betasq) ;
322  }
323 
324  inline double d2func(const double t, const int np, double *p) const{
325  double betasq = SQR(p[0]);
326  double hx =sqrt(SQR(t) + betasq);
327  // return CUBE(p[0]/hx);
328  return p[0]/sqrt(SQR(t) + betasq);
329  // this is "Gauss-Newton" Hessian : convex
330  }
331 };
332 
333 class STOASTLIB TV: public TVSigma {
334 public:
335  TV (double _tau, double _beta, const RVector *_x0, const Raster *_raster,
336  const void *_kref=0, bool _KRefTensFlag = false)
337  : TVSigma ( _tau, _beta, 0, _x0, _raster,false, _kref, _KRefTensFlag)
338  {}
339 
340  TV (double _tau, double _beta, const RVector *_x0, const Raster *_raster,
341  const RVector &_refimg, double _sdr, double _fT,
342  bool _KRefTensFlag = false)
343  : TVSigma (_tau, _beta, 0, _x0, _raster, _refimg, _sdr, _fT, false,
344  _KRefTensFlag)
345  {}
346 
347  ~TV ()
348  {}
349 
350  const char *GetName() const
351  { return "TV"; }
352 };
353 
354 // ======================== Laplacian = 1st order Tikhonov ====================
355 class STOASTLIB TK1Sigma: public GenericSigma {
356 public:
357  TK1Sigma (double _tau, double _sd, const RVector *_x0,
358  const Raster *_raster, bool _SmoothKappaOnly = true,
359  const void *_kref=0, bool _KRefTensFlag =false)
360  : GenericSigma (_raster, _x0, 0, (double*)0, _tau, _sd,
361  (_SmoothKappaOnly ? 0 : _sd), _kref, _KRefTensFlag)
362  {}
363 
364  TK1Sigma (double _tau, double _sd, const RVector *_x0,
365  const Raster *_raster, const RVector &_refimg, double _sdr, double _fT,
366  bool _SmoothKappaOnly = true, bool _KRefTensFlag = false)
367  : GenericSigma (_raster, _x0, 0, (double*)0, _tau, _refimg,
368  _sdr, _sd, (_SmoothKappaOnly ? 0 : _sd), _fT, _KRefTensFlag)
369  {}
370 
371  ~TK1Sigma ()
372  {}
373 
374  const char *GetName() const
375  { return "TK1_SIGMA"; }
376 
377 protected:
378  inline double func(const double t, const int np, double *p) const{
379  return 0.5*SQR(t);
380  }
381 
382  inline double kfunc(const double t, const int np, double *p) const{
383  return 1;
384  }
385 
386  inline double d2func(const double t, const int np, double *p) const{
387  return 1;
388  }
389 };
390 
391 class STOASTLIB TK1: public TK1Sigma {
392 public:
393  TK1 (double _tau, const RVector *_x0, const Raster *_raster,
394  const void *_kref = 0, bool _KRefTensFlag = false)
395  : TK1Sigma ( _tau, 0, _x0, _raster, false, _kref, _KRefTensFlag)
396  {}
397 
398  TK1 (double _tau, const RVector *_x0, const Raster *_raster,
399  const RVector &_refimg, double _sdr, double _fT,
400  bool _KRefTenseFlag = false)
401  : TK1Sigma (_tau, 0, _x0, _raster, _refimg, _sdr, _fT, false,
402  _KRefTenseFlag)
403  {}
404 
405  ~TK1 ()
406  {}
407 
408  const char *GetName() const
409  { return "TK1"; }
410 };
411 
412 // =========================== Huber Function ===========================
413 class STOASTLIB HuberSigma: public GenericSigma {
414 public:
415  HuberSigma (double _tau, double _eps, double _sd, const RVector *_x0,
416  const Raster *_raster, bool _SmoothKappaOnly = true,
417  const void *_kref = 0, bool _KRefTensFlag = false)
418  : GenericSigma (_raster, _x0, 1, (new double [1]), _tau, _sd,
419  (_SmoothKappaOnly ? 0 : _sd), _kref, _KRefTensFlag)
420  {preg[0] = _eps;}
421 
422  HuberSigma (double _tau, double _eps, double _sd, const RVector *_x0,
423  const Raster *_raster, const RVector _refimg, double _sdr, double _fT,
424  bool _SmoothKappaOnly = true, bool _KRefTensFlag = false)
425  : GenericSigma (_raster, _x0, 1, (new double[1]), _tau, _refimg,
426  _sdr, _sd, (_SmoothKappaOnly ? 0 : _sd), _fT, _KRefTensFlag)
427  {preg[0] = _eps;}
428 
429  ~HuberSigma ()
430  {}
431 
432  const char *GetName() const
433  { return "HUBER_SIGMA"; }
434 
435  void ReadParams (ParamParser *pp);
436  void WriteParams (ParamParser *pp);
437 
438 protected:
439  inline double func(const double t, const int np, double *p) const{
440  double eps = p[0];
441  double at = fabs(t);
442  return (at <= eps ? 0.5*SQR(t) : eps*(at - 0.5*eps));
443  }
444 
445  inline double kfunc(const double t, const int np, double *p) const{
446  double eps = p[0];
447  double at = fabs(t);
448  return (at <= eps ? 1 : eps/at);
449  }
450 
451  inline double d2func(const double t, const int np, double *p) const{
452  double eps = p[0];
453  double at = fabs(t);
454  return (at <= eps ? 1 : CUBE(eps/at)); // not standard!
455  }
456 };
457 
458 class STOASTLIB Huber: public HuberSigma {
459 public:
460  Huber (double _tau, double _eps, const RVector *_x0, const Raster *_raster,
461  const void *_kref = 0, bool _KRefTensFlag = false)
462  : HuberSigma (_tau, _eps, 0, _x0, _raster, false, _kref, _KRefTensFlag)
463  {}
464 
465  Huber (double _tau, double _eps, const RVector *_x0, const Raster *_raster,
466  const RVector &_refimg, double _sdr, double _fT,
467  bool _KRefTensFlag = false)
468  : HuberSigma (_tau, _eps, 0, _x0, _raster, _refimg, _sdr, _fT, false,
469  _KRefTensFlag)
470  {}
471 
472  ~Huber ()
473  {}
474 
475  const char *GetName() const
476  { return "HUBER"; }
477 };
478 // =========================================================================
479 class STOASTLIB PMSigma: public GenericSigma {
480 public:
481  PMSigma (double _tau, double _T, double _sd, const RVector *_x0,
482  const Raster *_raster, bool _SmoothKappaOnly = true,
483  const void *_kref = 0, bool _KRefTensFlag = false)
484  : GenericSigma (_raster, _x0, 1, (new double [1]), _tau, _sd,
485  (_SmoothKappaOnly ? 0 : _sd), _kref, _KRefTensFlag)
486  {preg[0] = _T;}
487 
488  PMSigma (double _tau, double _T, double _sd, const RVector *_x0,
489  const Raster *_raster, const RVector &_refimg, double _sdr, double _fT,
490  bool _SmoothKappaOnly = true, bool _KRefTensFlag = false)
491  : GenericSigma (_raster, _x0, 1, (new double[1]), _tau, _refimg,
492  _sdr, _sd, (_SmoothKappaOnly ? 0 : _sd), _fT, _KRefTensFlag)
493  {preg[0] = _T;}
494 
495  ~PMSigma ()
496  {}
497 
498 protected:
499  inline double func(const double t, const int np, double *p) const{
500  double T = p[0];
501  return 0.5*SQR(T)*(1-exp(-SQR(t/T)));
502  }
503 
504  inline double kfunc(const double t, const int np, double *p) const{
505  double T = p[0];
506  return exp(-SQR(t/T));
507  }
508 
509  inline double d2func(const double t, const int np, double *p) const{
510  double T = p[0];
511  // return (1- 2*SQR(t/T))*exp(-SQR(t/T));
512  return exp(-SQR(t/T)); // this is "Gauss-Newton" Hessian : convex
513  }
514 };
515 
516 class STOASTLIB PM: public PMSigma {
517 public:
518  PM (double _tau, double _T, const RVector *_x0, const Raster *_raster,
519  const void *_kref = 0, bool _KRefTensFlag = false)
520  : PMSigma (_tau, _T, 0, _x0, _raster, false, _kref, _KRefTensFlag)
521  {}
522 
523  PM (double _tau, double _T, const RVector *_x0, const Raster *_raster,
524  const RVector &_refimg, double _sdr, double _fT,
525  bool _KRefTensFlag = false)
526  : PMSigma (_tau, _T, 0, _x0, _raster, _refimg, _sdr, _fT, false,
527  _KRefTensFlag)
528  {}
529 
530  ~PM ()
531  {}
532 };
533 // =========================================================================
534 class STOASTLIB QPMSigma: public GenericSigma {
535 public:
536  QPMSigma (double _tau, double _T, double _sd, const RVector *_x0,
537  const Raster *_raster, bool _SmoothKappaOnly = true,
538  const void *_kref = 0, bool _KRefTensFlag = false)
539  : GenericSigma (_raster, _x0, 1, (new double[1]), _tau, _sd,
540  (_SmoothKappaOnly ? 0 : _sd), _kref, _KRefTensFlag)
541  {preg[0] = _T;}
542 
543  QPMSigma (double _tau, double _T, double _sd, const RVector *_x0,
544  const Raster *_raster, const RVector &_refimg, double _sdr, double _fT,
545  bool _SmoothKappaOnly = true, bool _KRefTensFlag = false)
546  : GenericSigma (_raster, _x0, 1, (new double[1]), _tau, _refimg,
547  _sdr, _sd, (_SmoothKappaOnly ? 0 : _sd), _fT, _KRefTensFlag)
548  {preg[0] = _T;}
549 
550  ~QPMSigma ()
551  {}
552 
553 protected:
554  inline double func(const double t, const int np, double *p) const{
555  double T = p[0];
556  return 0.5*SQR(T)*log(SQR(t/T) + 1);
557  }
558 
559  inline double kfunc(const double t, const int np, double *p) const{
560  double T = p[0];
561  return 1/(SQR(t/T) + 1) ;
562  }
563 
564  inline double d2func(const double t, const int np, double *p) const{
565  double T = p[0];
566  double hx = (SQR(t/T) + 1);
567  // return (1-SQR(t/T))/SQR(hx); // this looks a bit dodgy : not sym.pos.def
568  return hx; // this is "Gauss-Newton" Hessian : convex
569  }
570 };
571 
572 class STOASTLIB QPM: public QPMSigma {
573 public:
574  QPM (double _tau, double _T, const RVector *_x0, const Raster *_raster,
575  const void *_kref = 0, bool _KRefTensFlag = false)
576  : QPMSigma (_tau, _T, 0, _x0, _raster, false, _kref, _KRefTensFlag)
577  {}
578 
579  QPM (double _tau, double _T, const RVector *_x0, const Raster *_raster,
580  const RVector &_refimg, double _sdr, double _fT,
581  bool _KRefTensFlag = false)
582  : QPMSigma (_tau, _T, 0, _x0, _raster, _refimg, _sdr, _fT, false,
583  _KRefTensFlag)
584  {}
585 
586  ~QPM ()
587  {}
588 };
589 
590 // ====================== Tukey BiWeight ===============================
591 class STOASTLIB TukeySigma: public GenericSigma {
592 public:
593  TukeySigma (double _tau, double _T, double _sd, const RVector *_x0,
594  const Raster *_raster, bool _SmoothKappaOnly = true,
595  const void *_kref = 0, bool _KRefTensFlag = false)
596  : GenericSigma (_raster, _x0, 1, (new double[1]), _tau, _sd,
597  (_SmoothKappaOnly ? 0 : _sd), _kref, _KRefTensFlag)
598  {preg[0] = _T;}
599 
600  TukeySigma (double _tau, double _T, double _sd, const RVector *_x0,
601  const Raster *_raster, const RVector &_refimg, double _sdr, double _fT,
602  bool _SmoothKappaOnly = true, bool _KRefTensFlag = false)
603  : GenericSigma (_raster, _x0, 1, (new double[1]), _tau, _refimg,
604  _sdr, _sd, (_SmoothKappaOnly ? 0 : _sd), _fT, _KRefTensFlag)
605  {preg[0] = _T;}
606 
607  ~TukeySigma ()
608  {}
609 
610 protected:
611  inline double func(const double t, const int np, double *p) const{
612  double T = p[0];
613  return 0.5*SQR(t)*(1 - SQR(t/T) + SQR(SQR(t/T))/3.0 );
614  }
615 
616  inline double kfunc(const double t, const int np, double *p) const{
617  double T = p[0];
618  double at = fabs(t);
619  return (at <= T ? SQR(1 - (SQR(t/T))) : 0);
620  }
621 
622  inline double d2func(const double t, const int np, double *p) const{
623  double T = p[0];
624  double at = fabs(t);
625  return (at <= T ? SQR(1 - (SQR(t/T))) : 0);
626  }
627 };
628 
629 class STOASTLIB Tukey: public TukeySigma {
630 public:
631  Tukey (double _tau, double _T, const RVector *_x0, const Raster *_raster,
632  const void *_kref = 0, bool _KRefTensFlag = false)
633  : TukeySigma (_tau, _T, 0, _x0, _raster, false, _kref, _KRefTensFlag)
634  {}
635 
636  Tukey (double _tau, double _T, const RVector *_x0, const Raster *_raster,
637  const RVector &_refimg, double _sdr, double _fT,
638  bool _KRefTensFlag = false)
639  : TukeySigma (_tau, _T, 0, _x0, _raster, _refimg, _sdr, _fT, false,
640  _KRefTensFlag)
641  {}
642 
643  ~Tukey ()
644  {}
645 };
646 
647 // =========================================================================
648 
649 class STOASTLIB GenericScaleSpace : public Regularisation {
650 public:
651  GenericScaleSpace (const Raster *_raster, const RVector *_x0, PSIREG _func, PSIREG _kfunc, PSIREG _k1func, const int _npreg, double * _preg, double _tau, double _sdx = 0, double _sdf = 0, const RDenseMatrix *_kreftens=0);
652 
653  ~GenericScaleSpace ();
654 
655  PRIOR Type() const {return PRIOR_GENERICSCALESPACE; }
656 
657  // void SetHessian (RCompRowMatrix &Hess, const RVector &x) {
658  // CreateNonLinearHessian(raster, x, Hess);
659  // };
660 
661  double GetValue (const RVector &x) const;
662  RVector GetGradient (const RVector &x) const {
663  // MS: modified 100201
664  // the second parameter of GetHess1f doesn't get x0 subtracted
665  // inside GetHess1f, which can lead to a nonzero gradient even
666  // if x = x0
667 
668  //return GetHess1f(x,x);
669  return GetHess1f (x, x - *x0);
670  //return GetHess1f(x,x);
671  }
672  RVector GetKappa (const RVector &x) const;
673  void SetHess1 (RCompRowMatrix &Hess1, const RVector &x, const int p);
674  void SetHess1FromKappa (RCompRowMatrix &Hess, const RVector &kap) ;
675  RVector GetHess1f (const RVector &x, const RVector &f) const {
676  if(KappaSet)
677  return GetHess1f_KappaSet(f);
678  else
679  return GetHess1f_KappaNotSet(x,f);
680  }
681  void SetFullHess (RCompRowMatrix &Hess, const RVector &x, const int p);
682  RVector GetFullHessf (const RVector &x, const RVector &f) const;
683  // void SetHessian (RCompRowMatrix &Hess, const RVector &x) ;
684 
685  int GetHessianRow (const RVector &x, int i, idxtype *colidx, double *val) const;
686 
687  RVector GetHessianDiag (const RVector &x) const ;
688 
689 protected:
690  RVector GetHess1f_KappaSet (const RVector &f) const;
691  RVector GetHess1f_KappaNotSet (const RVector &x, const RVector &f) const;
692 
693  const double sdx, sdf ;
694  const PSIREG func, kfunc, k1func;
695  const int npreg;
696  const RDenseMatrix *kreftens;
697  RDenseMatrix **KappaTens;
698  bool KappaSet;
699  double *preg;
700 };
701 // =====================Pretty Much Redundant==================================
702 
703 class STOASTLIB Generic : public GenericSigma {// zero scale version of GenericSigma
704 
705 public:
706  Generic (const Raster *_raster, const RVector *_x0, const int _npreg, double * _preg, double _tau, const RVector *_kap=0) : GenericSigma (_raster, _x0, _npreg, _preg, _tau, 0, 0,_kap) {};
707 };
708 
709 #define OLD_REGUL2
710 #ifdef OLD_REGUL2
711 // =========================================================================
712 
713 class STOASTLIB NullRegularisation: public Regularisation {
714 public:
715  NullRegularisation (): Regularisation() {}
716  PRIOR Type() const { return PRIOR_NONE; }
717  const char *GetName() const { return "NONE"; }
718  double GetValue (const RVector &x) const;
719  RVector GetGradient (const RVector &x) const;
720  RVector GetKappa (const RVector &x) const;
721  RVector GetHess1f (const RVector &x, const RVector &f) const;
722 
723  void SetHessian (RCompRowMatrix &Hess, const RVector &x) {}
724  void SetHessianFromKappa(RCompRowMatrix &Hess, const RVector &kappa) {}
725  int GetHessianRow (const RVector &x, int i, idxtype *colidx, double *val)
726  const;
727  RVector GetHessianDiag (const RVector &x) const;
728 
729  /* undefined functions */
730  void SetHess1 (RCompRowMatrix &Hess1, const RVector &x, const int p) {}
731  void SetHess1FromKappa (RCompRowMatrix &Hess, const RVector &kap) {}
732  void SetFullHess (RCompRowMatrix &Hess, const RVector &x, const int p) {}
733  RVector GetFullHessf (const RVector &x, const RVector &f) const { return RVector(); }
734 
735 };
736 // =========================================================================
737 
738 class STOASTLIB Tikhonov0: public Regularisation {
739 public:
740  Tikhonov0 (double _tau, const RVector *_x0, const RVector *_xs);
741  ~Tikhonov0();
742  PRIOR Type() const {return PRIOR_DIAG; }
743 
744  const char *GetName() const { return "TK0"; }
745  double GetValue (const RVector &x) const;
746  RVector GetGradient (const RVector &x) const;
747  RVector GetKappa (const RVector &x) const;
748  RVector GetHess1f (const RVector &x, const RVector &f) const;
749  void SetHessian (RCompRowMatrix &Hess, const RVector &x);
750  void SetHessianFromKappa(RCompRowMatrix &Hess, const RVector &kappa);
751  int GetHessianRow (const RVector &x, int i, idxtype *colidx, double *val)
752  const;
753  RVector GetHessianDiag (const RVector &x) const;
754  void SetTau (double _tau);
755  void SetTau (const RVector &_tau);
756 
757  /* undefined functions */
758  void SetHess1 (RCompRowMatrix &Hess1, const RVector &x, const int p) {}
759  void SetHess1FromKappa (RCompRowMatrix &Hess, const RVector &kap) {}
760  void SetFullHess (RCompRowMatrix &Hess, const RVector &x, const int p) {}
761  RVector GetFullHessf (const RVector &x, const RVector &f) const { return RVector(); }
762 
763  void ReadParams (ParamParser *pp);
764  void WriteParams (ParamParser *pp);
765 
766 private:
767  RVector tauvec;
768  bool use_tauvec;
769  const RVector *xs;
770  char *kaprefimg_name;
771 };
772 
773 // =========================================================================
774 
775 class STOASTLIB Tikhonov1: public Regularisation {
776 public:
777  Tikhonov1 (double _tau, const RVector *_x0, const Raster *raster, const RVector *_kap=0);
778  ~Tikhonov1 ();
779  PRIOR Type() const { return PRIOR_LAPLACIAN; }
780  const char *GetName() const { return "Tikhonov1"; }
781  double GetValue (const RVector &x) const;
782  RVector GetGradient (const RVector &x) const;
783  RVector GetKappa (const RVector &x) const;
784  RVector GetHess1f (const RVector &x, const RVector &f) const;
785  void SetHessian (RCompRowMatrix &Hess, const RVector &x) {
786  Hess= *Laplacian;
787  };
788  void SetHessianFromKappa(RCompRowMatrix &Hess, const RVector &kappa) {
789  // CreateHessian(raster,kappa,Hess);
790  };
791  int GetHessianRow (const RVector &x, int i, idxtype *colidx, double *val)
792  const;
793  RVector GetHessianDiag (const RVector &x) const;
794  void SetHess1 (RCompRowMatrix &Hess1, const RVector &x, const int p);
795 
796  /* undefined functions */
797  void SetHess1FromKappa (RCompRowMatrix &Hess, const RVector &kap) {}
798  void SetFullHess (RCompRowMatrix &Hess, const RVector &x, const int p) {}
799  RVector GetFullHessf (const RVector &x, const RVector &f) const { return RVector(); }
800 
801 private:
802  void CreateHessian (const Raster *raster, const RVector &kappa,
803  RCompRowMatrix &Hess);
804 
805  const RVector *kappa;
806  RCompRowMatrix *Laplacian;
807  // int nset, pdim;
808 };
809 
810 // =========================================================================
811 
812 class STOASTLIB TVSigmaOLD: public Regularisation {
813 public:
814  TVSigmaOLD (double _tau, double _beta, double _sd, const RVector *_x0, const Raster *_raster, bool _SmoothKappaOnly = true, const RVector *_kap=0);
815  ~TVSigmaOLD ();
816 
817  PRIOR Type() const {return PRIOR_TVSIGMA; }
818  double GetValue (const RVector &x) const;
819  RVector GetGradient (const RVector &x) const;
820  RVector GetKappa (const RVector &x) const;
821  RVector GetHess1f (const RVector &x, const RVector &f) const;
822  void SetHessian (RCompRowMatrix &Hess, const RVector &x) {
823  // CreateNonLinearHessian(raster, x, Hess) ;
824  };
825  void SetHessianFromKappa(RCompRowMatrix &Hess, const RVector &_kap) {
826  // CreateHessian(raster,_kap,Hess);
827  };
828  void OLDSetHessian (const RVector &x) const;
829  int GetHessianRow (const RVector &x, int i, idxtype *colidx, double *val)
830  const;
831  RVector GetHessianDiag (const RVector &x) const;
832  /* undefined functions */
833  void SetHess1 (RCompRowMatrix &Hess1, const RVector &x, const int p) {}
834  void SetHess1FromKappa (RCompRowMatrix &Hess, const RVector &kap) {}
835  void SetFullHess (RCompRowMatrix &Hess, const RVector &x, const int p) {}
836  RVector GetFullHessf (const RVector &x, const RVector &f) const { return RVector(); }
837 
838 protected:
839  double tau, beta, betasq, sd;
840  bool SmoothKappaOnly ;
841  const RVector *kappa;
842  bool delete_kappa; // this isn't at all pretty...
843  RCompRowMatrix *TVhess;
844  // int nset, pdim;
845 };
846 
847 // =========================================================================
848 
849 class STOASTLIB TVOLD: public TVSigmaOLD { // zero scale version of TVSigma
850 public:
851  TVOLD (double _tau, double _beta, const RVector *_x0, const Raster *_raster, const RVector *_kap=0): TVSigmaOLD(_tau,_beta,0,_x0, _raster, false, _kap) {};
852  PRIOR Type() const {return PRIOR_TV; }
853 
854 };
855 
856 // =========================================================================
857 
858 class STOASTLIB MRF: public Regularisation {
859 public:
860  MRF (double _tau, const RVector *_x0, const Raster *raster, const RVector *_kap = 0);
861  ~MRF ();
862  PRIOR Type() const { return PRIOR_MRF; }
863  const char *GetName() const { return "MRF"; }
864 
865  void SetKaprefImg (const RVector *kap = 0);
866 
867  double GetValue (const RVector &x) const;
868  RVector GetGradient (const RVector &x) const;
869 
870  // not implemented (yet)
871  RVector GetKappa (const RVector &x) const;
872  void SetHess1 (RCompRowMatrix &Hess1, const RVector &x, const int p);
873  void SetHess1FromKappa (RCompRowMatrix &Hess, const RVector &kap);
874  RVector GetHess1f (const RVector &x, const RVector &f) const;
875  void SetFullHess (RCompRowMatrix &Hess, const RVector &x, const int p);
876  RVector GetFullHessf (const RVector &x, const RVector &f) const;
877  int GetHessianRow (const RVector &x, int i, idxtype *colidx, double *val)
878  const;
879  RVector GetHessianDiag (const RVector &x) const;
880 
881 private:
882  RCompRowMatrix *MRFa, *MRFb;
883 };
884 
885 #endif // !OLD_REGUL2
886 #endif // !__ST_REGUL_H
Huber
Definition: regul.h:458
TVector
Templated vector class.
Definition: vector.h:39
TVSigma
Definition: regul.h:288
ParamParser
Definition: pparse.h:11
TK1Sigma
Definition: regul.h:355
Tikhonov0
Definition: regul.h:738
QPM
Definition: regul.h:572
GenericScaleSpace
Definition: regul.h:649
Raster
Base class for mapping between mesh and an independent basis representation.
Definition: raster.h:20
Tukey
Definition: regul.h:629
GenericSigma
Definition: regul.h:135
QPMSigma
Definition: regul.h:534
TK1
Definition: regul.h:391
Tikhonov1
Definition: regul.h:775
Generic
Definition: regul.h:703
HuberSigma
Definition: regul.h:413
PM
Definition: regul.h:516
Regularisation
Definition: regul.h:42
TV
Definition: regul.h:333
TVSigmaOLD
Definition: regul.h:812
MRF
Definition: regul.h:858
TDenseMatrix
Dense matrix class.
Definition: crmatrix.h:38
PMSigma
Definition: regul.h:479
TukeySigma
Definition: regul.h:591
TVOLD
Definition: regul.h:849
NullRegularisation
Definition: regul.h:713
Generated on Tue Sep 2 2014 17:14:26 for Toast++ by   doxygen 1.8.6