NPL
Neurological Programs and Libraries
accessors.h
Go to the documentation of this file.
1 /******************************************************************************
2  * Copyright 2014 Micah C Chambers (micahc.vt@gmail.com)
3  *
4  * NPL is free software: you can redistribute it and/or modify it under the
5  * terms of the BSD 2-Clause License available in LICENSE or at
6  * http://opensource.org/licenses/BSD-2-Clause
7  *
8  * @file accessors.h Provides accessors to the NDArray data structure. Due to
9  * the fact that dimensionality and pixel type are not carried around in the
10  * container type, it is necessary to provide generic accessors that will cast
11  * input data to the correct type and return it, or take input, cast it to the
12  * underlying type and set the internal pixel. While this might seem
13  * round-about, it allows you to write general purpose algorithms functions
14  * without concerning with dimensionality or pixel type.
15  *
16  *****************************************************************************/
17 
18 #ifndef ACCESSORS_H
19 #define ACCESSORS_H
20 
21 #include <stdexcept>
22 
23 #include "ndarray.h"
24 #include "mrimage.h"
25 #include "basic_functions.h"
26 #include "utility.h"
27 #include "iterators.h"
28 
29 namespace npl {
30 
62 template<typename T>
63 class NDView
64 {
65 public:
66  NDView(std::shared_ptr<NDArray> in) : parent(in)
67  {
68  setArray(in);
69  };
70 
71  NDView() : parent(NULL) {} ;
72 
73  void setArray(ptr<NDArray> in)
74  {
75  parent = in;
76  switch(in->type()) {
77  case UINT8:
78  castget = castgetStatic<uint8_t>;
79  castset = castsetStatic<uint8_t>;
80  break;
81  case INT8:
82  castget = castgetStatic<int8_t>;
83  castset = castsetStatic<int8_t>;
84  break;
85  case UINT16:
86  castget = castgetStatic<uint16_t>;
87  castset = castsetStatic<uint16_t>;
88  break;
89  case INT16:
90  castget = castgetStatic<int16_t>;
91  castset = castsetStatic<int16_t>;
92  break;
93  case UINT32:
94  castget = castgetStatic<uint32_t>;
95  castset = castsetStatic<uint32_t>;
96  break;
97  case INT32:
98  castget = castgetStatic<int32_t>;
99  castset = castsetStatic<int32_t>;
100  break;
101  case UINT64:
102  castget = castgetStatic<uint64_t>;
103  castset = castsetStatic<uint64_t>;
104  break;
105  case INT64:
106  castget = castgetStatic<int64_t>;
107  castset = castsetStatic<int64_t>;
108  break;
109  case FLOAT32:
110  castget = castgetStatic<float>;
111  castset = castsetStatic<float>;
112  break;
113  case FLOAT64:
114  castget = castgetStatic<double>;
115  castset = castsetStatic<double>;
116  break;
117  case FLOAT128:
118  castget = castgetStatic<long double>;
119  castset = castsetStatic<long double>;
120  break;
121  case COMPLEX64:
122  castget = castgetStatic<cfloat_t>;
123  castset = castsetStatic<cfloat_t>;
124  break;
125  case COMPLEX128:
126  castget = castgetStatic<cdouble_t>;
127  castset = castsetStatic<cdouble_t>;
128  break;
129  case COMPLEX256:
130  castget = castgetStatic<cquad_t>;
131  castset = castsetStatic<cquad_t>;
132  break;
133  case RGB24:
134  castget = castgetStatic<rgb_t>;
135  castset = castsetStatic<rgb_t>;
136  break;
137  case RGBA32:
138  castget = castgetStatic<rgba_t>;
139  castset = castsetStatic<rgba_t>;
140  break;
141  default:
142  case UNKNOWN_TYPE:
143  castget = castgetStatic<uint8_t>;
144  castset = castsetStatic<uint8_t>;
145  throw std::invalid_argument("Unknown type to NDView");
146  break;
147  }
148  }
149 
155  T operator[](int64_t index)
156  {
157  auto ptr = this->parent->__getAddr(index);
158  assert(ptr >= this->parent->__getAddr(0) &&
159  ptr < this->parent->__getAddr(this->parent->elements()));
160  return castget(ptr);
161  };
162 
170  T get(const std::vector<int64_t>& index)
171  {
172  auto ptr = this->parent->__getAddr(index);
173  assert(ptr >= this->parent->__getAddr(0) &&
174  ptr < this->parent->__getAddr(this->parent->elements()));
175  return castget(ptr);
176  };
177 
186  T get(size_t len, int64_t* index)
187  {
188  auto ptr = this->parent->__getAddr(len, index);
189  assert(ptr >= this->parent->__getAddr(0) &&
190  ptr < this->parent->__getAddr(this->parent->elements()));
191  return castget(ptr);
192  };
193 
201  T operator[](const std::vector<int64_t>& index)
202  {
203  auto ptr = this->parent->__getAddr(index);
204  assert(ptr >= this->parent->__getAddr(0) &&
205  ptr < this->parent->__getAddr(this->parent->elements()));
206  return castget(ptr);
207  };
208 
218  void set(size_t len, const int64_t* index, T v)
219  {
220  auto ptr = this->parent->__getAddr(len, index);
221  assert(ptr >= this->parent->__getAddr(0) &&
222  ptr < this->parent->__getAddr(this->parent->elements()));
223  return castset(ptr, v);
224  };
225 
234  void set(const std::vector<int64_t>& index, T v)
235  {
236  auto ptr = this->parent->__getAddr(index);
237  assert(ptr >= this->parent->__getAddr(0) &&
238  ptr < this->parent->__getAddr(this->parent->elements()));
239  return castset(ptr, v);
240  };
241 
250  void set(int64_t index, T v)
251  {
252  auto ptr = this->parent->__getAddr(index);
253  assert(ptr >= this->parent->__getAddr(0) &&
254  ptr < this->parent->__getAddr(this->parent->elements()));
255  return castset(ptr, v);
256  };
257 
258  int64_t tlen() { return this->parent->tlen(); };
259 
260 protected:
261 
271  template <typename U>
272  static T castgetStatic(void* ptr)
273  {
274  return (T)(*(static_cast<U*>(ptr)));
275  };
276 
286  template <typename U>
287  static void castsetStatic(void* ptr, const T& val)
288  {
289  (*(static_cast<U*>(ptr))) = (U)val;
290  };
291 
295  std::shared_ptr<NDArray> parent;
296 
303  T (*castget)(void* ptr);
304 
313  void (*castset)(void* ptr, const T& v);
314 
315 };
316 
323 template<typename T>
324 class NDConstView
325 {
326 public:
327  NDConstView(std::shared_ptr<const NDArray> in) : parent(in)
328  {
329  setArray(in);
330  }
331 
332  NDConstView() : parent(NULL) {} ;
333 
334  void setArray(ptr<const NDArray> in)
335  {
336  parent = in;
337  switch(in->type()) {
338  case UINT8:
339  castget = castgetStatic<uint8_t>;
340  break;
341  case INT8:
342  castget = castgetStatic<int8_t>;
343  break;
344  case UINT16:
345  castget = castgetStatic<uint16_t>;
346  break;
347  case INT16:
348  castget = castgetStatic<int16_t>;
349  break;
350  case UINT32:
351  castget = castgetStatic<uint32_t>;
352  break;
353  case INT32:
354  castget = castgetStatic<int32_t>;
355  break;
356  case UINT64:
357  castget = castgetStatic<uint64_t>;
358  break;
359  case INT64:
360  castget = castgetStatic<int64_t>;
361  break;
362  case FLOAT32:
363  castget = castgetStatic<float>;
364  break;
365  case FLOAT64:
366  castget = castgetStatic<double>;
367  break;
368  case FLOAT128:
369  castget = castgetStatic<long double>;
370  break;
371  case COMPLEX64:
372  castget = castgetStatic<cfloat_t>;
373  break;
374  case COMPLEX128:
375  castget = castgetStatic<cdouble_t>;
376  break;
377  case COMPLEX256:
378  castget = castgetStatic<cquad_t>;
379  break;
380  case RGB24:
381  castget = castgetStatic<rgb_t>;
382  break;
383  case RGBA32:
384  castget = castgetStatic<rgba_t>;
385  break;
386  default:
387  case UNKNOWN_TYPE:
388  castget = castgetStatic<uint8_t>;
389  throw std::invalid_argument("Unknown type to NDConstView");
390  break;
391  }
392  };
393 
399  T operator[](int64_t index)
400  {
401  auto ptr = this->parent->__getAddr(index);
402  assert(ptr >= this->parent->__getAddr(0) &&
403  ptr < this->parent->__getAddr(this->parent->elements()));
404  return castget(ptr);
405  };
406 
414  T get(const std::vector<int64_t>& index)
415  {
416  auto ptr = this->parent->__getAddr(index);
417  assert(ptr >= this->parent->__getAddr(0) &&
418  ptr < this->parent->__getAddr(this->parent->elements()));
419  return castget(ptr);
420  };
421 
430  T get(size_t len, int64_t* index)
431  {
432  auto ptr = this->parent->__getAddr(len, index);
433  assert(ptr >= this->parent->__getAddr(0) &&
434  ptr < this->parent->__getAddr(this->parent->elements()));
435  return castget(ptr);
436  };
437 
445  T operator[](const std::vector<int64_t>& index)
446  {
447  auto ptr = this->parent->__getAddr(index);
448  assert(ptr >= this->parent->__getAddr(0) &&
449  ptr < this->parent->__getAddr(this->parent->elements()));
450  return castget(ptr);
451  };
452 
453  int64_t tlen() { return this->parent->tlen(); };
454 
455 protected:
456 
466  template <typename U>
467  static T castgetStatic(void* ptr)
468  {
469  return (T)(*(static_cast<U*>(ptr)));
470  };
471 
475  std::shared_ptr<const NDArray> parent;
476 
483  T (*castget)(void* ptr);
484 
485 };
486 
494 template<typename T>
495 class Pixel3DView : public NDView<T>
496 {
497 public:
498  Pixel3DView(std::shared_ptr<NDArray> in) : NDView<T>(in)
499  { };
500 
501  Pixel3DView();
507  T operator()(int64_t x=0, int64_t y=0, int64_t z=0)
508  {
509  auto ptr = this->parent->__getAddr(x,y,z,0);
510  assert(ptr >= this->parent->__getAddr(0) &&
511  ptr < this->parent->__getAddr(this->parent->elements()));
512  return this->castget(ptr);
513  };
514 
520  T get(int64_t x, int64_t y, int64_t z)
521  {
522  auto ptr = this->parent->__getAddr(x,y,z,0);
523  assert(ptr >= this->parent->__getAddr(0) &&
524  ptr < this->parent->__getAddr(this->parent->elements()));
525  return this->castget(ptr);
526  };
527 
533  void set(int64_t x, int64_t y, int64_t z, T v)
534  {
535  auto ptr = this->parent->__getAddr(x,y,z,0);
536  assert(ptr >= this->parent->__getAddr(0) &&
537  ptr < this->parent->__getAddr(this->parent->elements()));
538  this->castset(ptr, v);
539  };
540 
541 protected:
542 
543  // Remove functions that aren't relevent from NDView
544  T operator[](int64_t i) { (void)(i); return T(); };
545  T get(const std::vector<int64_t>& i) { (void)(i); return T(); };
546  T operator[](const std::vector<int64_t>& i) { (void)(i); return T(); };
547 };
548 
556 template<typename T>
557 class Vector3DConstView : public NDConstView<T>
558 {
559 public:
560  Vector3DConstView(std::shared_ptr<const NDArray> in) : NDConstView<T>(in)
561  { };
562 
563  Vector3DConstView() {};
564 
570  virtual
571  T operator()(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
572  {
573  auto ptr = this->parent->__getAddr(x,y,z,t);
574  assert(ptr >= this->parent->__getAddr(0) &&
575  ptr < this->parent->__getAddr(this->parent->elements()));
576  return this->castget(ptr);
577  };
578 
584  virtual
585  T get(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
586  {
587  auto ptr = this->parent->__getAddr(x,y,z,t);
588  assert(ptr >= this->parent->__getAddr(0) &&
589  ptr < this->parent->__getAddr(this->parent->elements()));
590  return this->castget(ptr);
591  };
592 
598  virtual
599  T operator()(double x=0, double y=0, double z=0, int64_t t=0)
600  {
601  throw INVALID_ARGUMENT("DOUBLE PASSED TO VECTOR3D VIEW!");
602  auto ptr = this->parent->__getAddr(round(x),round(y),round(z),t);
603  assert(ptr >= this->parent->__getAddr(0) &&
604  ptr < this->parent->__getAddr(this->parent->elements()));
605  return this->castget(ptr);
606  };
607 
613  virtual
614  T get(double x=0, double y=0, double z=0, int64_t t=0)
615  {
616  throw INVALID_ARGUMENT("DOUBLE PASSED TO VECTOR3D VIEW!");
617  auto ptr = this->parent->__getAddr(round(x),round(y),round(z),t);
618  assert(ptr >= this->parent->__getAddr(0) &&
619  ptr < this->parent->__getAddr(this->parent->elements()));
620  return this->castget(ptr);
621  };
622 
623 private:
625  // Hide Non-3D Functrions from NDConstView
627 
633  T operator[](int64_t index)
634  {
635  auto ptr = this->parent->__getAddr(index);
636  assert(ptr >= this->parent->__getAddr(0) &&
637  ptr < this->parent->__getAddr(this->parent->elements()));
638  return castget(ptr);
639  };
640 
648  T get(const std::vector<int64_t>& index)
649  {
650  auto ptr = this->parent->__getAddr(index);
651  assert(ptr >= this->parent->__getAddr(0) &&
652  ptr < this->parent->__getAddr(this->parent->elements()));
653  return castget(ptr);
654  };
655 
664  T get(size_t len, int64_t* index)
665  {
666  auto ptr = this->parent->__getAddr(len, index);
667  assert(ptr >= this->parent->__getAddr(0) &&
668  ptr < this->parent->__getAddr(this->parent->elements()));
669  return castget(ptr);
670  };
671 
679  T operator[](const std::vector<int64_t>& index)
680  {
681  auto ptr = this->parent->__getAddr(index);
682  assert(ptr >= this->parent->__getAddr(0) &&
683  ptr < this->parent->__getAddr(this->parent->elements()));
684  return castget(ptr);
685  };
686 
687 };
688 
696 template<typename T>
697 class Vector3DView : public NDView<T>
698 {
699 public:
700  Vector3DView(std::shared_ptr<NDArray> in) : NDView<T>(in)
701  { };
702 
703  Vector3DView() {};
709  T operator()(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
710  {
711  auto ptr = this->parent->__getAddr(x,y,z,t);
712  assert(ptr >= this->parent->__getAddr(0) &&
713  ptr < this->parent->__getAddr(this->parent->elements()));
714  return this->castget(ptr);
715  };
716 
722  T get(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
723  {
724  auto ptr = this->parent->__getAddr(x,y,z,t);
725  assert(ptr >= this->parent->__getAddr(0) &&
726  ptr < this->parent->__getAddr(this->parent->elements()));
727  return this->castget(ptr);
728  };
729 
735  void set(int64_t x, int64_t y, int64_t z, int64_t t, T v)
736  {
737  auto ptr = this->parent->__getAddr(x,y,z,t);
738  assert(ptr >= this->parent->__getAddr(0) &&
739  ptr < this->parent->__getAddr(this->parent->elements()));
740  this->castset(ptr, v);
741  };
742 private:
744  // Hide Non-3D Functrions from NDConstView
746 
752  T operator[](int64_t index)
753  {
754  auto ptr = this->parent->__getAddr(index);
755  assert(ptr >= this->parent->__getAddr(0) &&
756  ptr < this->parent->__getAddr(this->parent->elements()));
757  return castget(ptr);
758  };
759 
767  T get(const std::vector<int64_t>& index)
768  {
769  auto ptr = this->parent->__getAddr(index);
770  assert(ptr >= this->parent->__getAddr(0) &&
771  ptr < this->parent->__getAddr(this->parent->elements()));
772  return castget(ptr);
773  };
774 
783  T get(size_t len, int64_t* index)
784  {
785  auto ptr = this->parent->__getAddr(len, index);
786  assert(ptr >= this->parent->__getAddr(0) &&
787  ptr < this->parent->__getAddr(this->parent->elements()));
788  return castget(ptr);
789  };
790 
798  T operator[](const std::vector<int64_t>& index)
799  {
800  auto ptr = this->parent->__getAddr(index);
801  assert(ptr >= this->parent->__getAddr(0) &&
802  ptr < this->parent->__getAddr(this->parent->elements()));
803  return castget(ptr);
804  };
805 
815  void set(size_t len, const int64_t* index, T v)
816  {
817  auto ptr = this->parent->__getAddr(len, index);
818  assert(ptr >= this->parent->__getAddr(0) &&
819  ptr < this->parent->__getAddr(this->parent->elements()));
820  return castset(ptr, v);
821  };
822 
831  void set(const std::vector<int64_t>& index, T v)
832  {
833  auto ptr = this->parent->__getAddr(index);
834  assert(ptr >= this->parent->__getAddr(0) &&
835  ptr < this->parent->__getAddr(this->parent->elements()));
836  return castset(ptr, v);
837  };
838 
847  void set(int64_t index, T v)
848  {
849  auto ptr = this->parent->__getAddr(index);
850  assert(ptr >= this->parent->__getAddr(0) &&
851  ptr < this->parent->__getAddr(this->parent->elements()));
852  return castset(ptr, v);
853  };
854 
855 };
856 
857 /****************************************
858  * Interpolators
859  ***************************************/
860 
862 //Linear
864 
865 /* Linear Kernel Sampling */
866 inline
867 double linKern(double x)
868 {
869  return fabs(1-fmin(1,fabs(x)));
870 }
871 
878 template<typename T>
879 class LinInterpNDView : public NDConstView<T>
880 {
881 public:
882  LinInterpNDView(std::shared_ptr<const NDArray> in,
883  BoundaryConditionT bound = ZEROFLUX)
884  : NDConstView<T>(in), m_boundmethod(bound), m_ras(false)
885  { };
886 
887  LinInterpNDView() : m_boundmethod(ZEROFLUX), m_ras(false) {} ;
888 
903  T operator()(double x=0, double y=0, double z=0, double t=0, double u = 0,
904  double v = 0, double w = 0, double q = 0)
905  {
906  double tmp[8] = {x,y,z,t,u,v,w,q};
907  return get(8, tmp);
908  };
909 
910 
918  T operator()(const std::vector<float>& index)
919  {
920  assert(index.size() <= 8);
921  double tmp[8];
922  size_t ii=0;
923  for(auto it = index.begin(); it != index.end() && ii<8; ++it, ++ii)
924  tmp[ii] = *it;
925  return get(std::min(8UL, index.size()), tmp);
926  };
927 
935  T operator()(const std::vector<double>& index)
936  {
937  return get(index.size(), index.data());
938  };
939 
947  T operator()(std::initializer_list<double> index)
948  {
949  assert(index.size() <= 8);
950  double tmp[8];
951  size_t ii=0;
952  for(auto it = index.begin(); it != index.end() && ii<8; ++it, ++ii)
953  tmp[ii] = *it;
954  return get(std::min(8UL, index.size()), tmp);
955  };
956 
964  T operator()(std::initializer_list<float> index)
965  {
966  assert(index.size() <= 8);
967  double tmp[8];
968  size_t ii=0;
969  for(auto it = index.begin(); it != index.end() && ii<8; ++it, ++ii)
970  tmp[ii] = *it;
971  return get(std::min(8UL, index.size()), tmp);
972  };
973 
979  T get(const vector<double>& cindex)
980  {
981  return get(cindex.size(), cindex.data());
982  }
983 
991  T get(const std::vector<int64_t>& index)
992  {
993  assert(index.size() <= 8);
994  double tmp[8];
995  size_t ii=0;
996  for(auto it = index.begin(); it != index.end() && ii<8; ++it, ++ii)
997  tmp[ii] = *it;
998  return get(std::min(8UL, index.size()), tmp);
999  };
1000 
1009  T get(size_t len, int64_t* index)
1010  {
1011  assert(len <= 8);
1012  double tmp[8];
1013  for(size_t ii=0; ii<len && ii<8; ii++)
1014  tmp[ii] = index[ii];
1015  return get(std::min(8UL, len), tmp);
1016  };
1017 
1025  T operator[](const std::vector<int64_t>& index)
1026  {
1027  assert(index.size() <= 8);
1028  double tmp[8];
1029  size_t ii=0;
1030  for(auto it = index.begin(); it != index.end() && ii<8; ++it, ++ii)
1031  tmp[ii] = *it;
1032  return get(std::min(8UL, index.size()), tmp);
1033  };
1034 
1040  T get(size_t len, const double* incindex)
1041  {
1042  // Initialize variables
1043  int ndim = this->parent->ndim();
1044  assert(ndim <= MAXDIM);
1045  const size_t* dim = this->parent->dim();
1046  int64_t index[MAXDIM];
1047  double cindex[MAXDIM];
1048 
1049  // convert RAS to index
1050  if(m_ras) {
1051  auto tmp = dPtrCast<const MRImage>(this->parent);
1052  tmp->pointToIndex(len, incindex, cindex);
1053  } else {
1054  for(size_t dd=0; dd<ndim; dd++) {
1055  if(dd < len)
1056  cindex[dd] = incindex[dd];
1057  else
1058  cindex[dd] = 0;
1059  }
1060  }
1061 
1062  Counter<> count(ndim);
1063  for(size_t dd=0; dd<ndim; dd++)
1064  count.sz[dd] = 2;
1065 
1066  // compute weighted pixval by iterating over neighbors
1067  T pixval = 0;
1068  do {
1069  double weight = 1;
1070  bool iioutside = false;
1071 
1072  //set index
1073  for(int dd = 0; dd < ndim; dd++) {
1074  index[dd] = floor(cindex[dd]) + count.pos[dd];
1075  weight *= linKern(index[dd] - cindex[dd]);
1076  iioutside = iioutside || index[dd] < 0 || index[dd] >= dim[dd];
1077  }
1078 
1079  // if the current point maps outside, then we need to deal with it
1080  if(iioutside) {
1081  if(m_boundmethod == ZEROFLUX) {
1082  // clamp
1083  for(size_t dd=0; dd<ndim; dd++)
1084  index[dd] = clamp<int64_t>(0, dim[dd]-1, index[dd]);
1085  } else if(m_boundmethod == WRAP) {
1086  // wrap
1087  for(size_t dd=0; dd<ndim; dd++)
1088  index[dd] = wrap<int64_t>(0, dim[dd]-1, index[dd]);
1089  } else {
1090  // set wieght to zero, then just clamp
1091  weight = 0;
1092  for(size_t dd=0; dd<ndim; dd++)
1093  index[dd] = clamp<int64_t>(0, dim[dd]-1, index[dd]);
1094  }
1095  }
1096 
1097  auto ptr = this->parent->__getAddr(ndim, index);
1098  assert(ptr >= this->parent->__getAddr(0) &&
1099  ptr < this->parent->__getAddr(this->parent->elements()));
1100  T v = this->castget(ptr);
1101  pixval += weight*v;
1102  } while(count.advance());
1103 
1104  return pixval;
1105  }
1106 
1107 
1108  BoundaryConditionT m_boundmethod;
1109 
1114  bool m_ras;
1115 
1116 protected:
1117 
1124  T operator[](int64_t i) { (void)(i); return T(); };
1125 };
1126 
1137 template<typename T>
1138 class LinInterp3DView : public Vector3DConstView<T>
1139 {
1140 public:
1141  LinInterp3DView(std::shared_ptr<const NDArray> in,
1142  BoundaryConditionT bound = ZEROFLUX)
1143  : Vector3DConstView<T>(in), m_boundmethod(bound), m_ras(false)
1144  { };
1145 
1146  LinInterp3DView() : m_boundmethod(ZEROFLUX), m_ras(false) {} ;
1147 
1153  T operator()(double x=0, double y=0, double z=0, int64_t t=0)
1154  {
1155  return get(x,y,z,t);
1156  };
1157 
1163  T get(double x=0, double y=0, double z=0, int64_t t=0)
1164  {
1165  // figure out size of dimensions in parent
1166  size_t dim[4];
1167  dim[0] = this->parent->dim(0);
1168  dim[1] = this->parent->ndim() > 1 ? this->parent->dim(1) : 1;
1169  dim[2] = this->parent->ndim() > 2 ? this->parent->dim(2) : 1;
1170  dim[3] = this->parent->tlen();
1171 
1172  // deal with t being outside bounds
1173  if(t < 0 || t >= dim[3]) {
1174  if(m_boundmethod == ZEROFLUX) {
1175  // clamp
1176  t = clamp<int64_t>(0, dim[3]-1, t);
1177  } else if(m_boundmethod == WRAP) {
1178  // wrap
1179  t = wrap<int64_t>(0, dim[3]-1, t);
1180  } else {
1181  return 0;
1182  }
1183  }
1184 
1185  // initialize variables
1186  double cindex[4] = {x,y,z,(double)t};
1187  int64_t index[4] = {0,0,0,t};
1188 
1189  // convert RAS to cindex
1190  if(m_ras) {
1191  auto tmp = dPtrCast<const MRImage>(this->parent);
1192  tmp->pointToIndex(3, cindex, cindex);
1193  }
1194 
1195  Counter<> count(3);
1196  for(size_t dd=0; dd<3; dd++)
1197  count.sz[dd] = 2;
1198 
1199  T pixval = 0;
1200  bool iioutside = false;
1201  do {
1202  double weight = 1;
1203 
1204  //set index
1205  for(int dd = 0; dd < 3; dd++) {
1206  index[dd] = floor(cindex[dd]) + count.pos[dd];
1207  weight *= linKern(index[dd] - cindex[dd]);
1208  iioutside = iioutside || index[dd] < 0 || index[dd] >= dim[dd];
1209  }
1210 
1211  // if the current point maps outside, then we need to deal with it
1212  if(iioutside) {
1213  if(m_boundmethod == ZEROFLUX) {
1214  // clamp
1215  for(size_t dd=0; dd<3; dd++)
1216  index[dd] = clamp<int64_t>(0, dim[dd]-1, index[dd]);
1217  } else if(m_boundmethod == WRAP) {
1218  // wrap
1219  for(size_t dd=0; dd<3; dd++)
1220  index[dd] = wrap<int64_t>(0, dim[dd]-1, index[dd]);
1221  } else {
1222  // set wieght to zero, then just clamp
1223  weight = 0;
1224  for(size_t dd=0; dd<3; dd++)
1225  index[dd] = clamp<int64_t>(0, dim[dd]-1, index[dd]);
1226  }
1227  }
1228 
1229  auto ptr = this->parent->__getAddr(index[0], index[1],index[2],t);
1230  assert(ptr >= this->parent->__getAddr(0) &&
1231  ptr < this->parent->__getAddr(this->parent->elements()));
1232  T v = this->castget(ptr);
1233  pixval += weight*v;
1234  } while(count.advance());
1235 
1236  return pixval;
1237  }
1238 
1244  T operator()(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
1245  {
1246  return get((double)x,(double)y,(double)z,t);
1247  };
1248 
1254  T get(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
1255  {
1256  return get((double)x,(double)y,(double)z,t);
1257  };
1258 
1259  BoundaryConditionT m_boundmethod;
1260 
1265  bool m_ras;
1266 
1267 };
1268 
1274 template<typename T>
1275 class NNInterpNDView : public NDConstView<T>
1276 {
1277 public:
1278  NNInterpNDView(std::shared_ptr<const NDArray> in,
1279  BoundaryConditionT bound = ZEROFLUX)
1280  : NDConstView<T>(in), m_boundmethod(bound), m_ras(false)
1281  { };
1282 
1283  NNInterpNDView() : m_boundmethod(ZEROFLUX), m_ras(false) {} ;
1284 
1299  T operator()(double x=0, double y=0, double z=0, double t=0, double u = 0,
1300  double v = 0, double w = 0, double q = 0)
1301  {
1302  double tmp[8] = {x,y,z,t,u,v,w,q};
1303  return get(8, tmp);
1304  };
1305 
1306 
1314  T operator()(const std::vector<float>& index)
1315  {
1316  assert(index.size() <= 8);
1317  double tmp[8];
1318  size_t ii=0;
1319  for(auto it = index.begin(); it != index.end() && ii<8; ++it, ++ii)
1320  tmp[ii] = *it;
1321  return get(std::min(8UL, index.size()), tmp);
1322  };
1323 
1331  T get(const std::vector<int64_t>& index)
1332  {
1333  assert(index.size() <= 8);
1334  double tmp[8];
1335  size_t ii=0;
1336  for(auto it = index.begin(); it != index.end() && ii<8; ++it, ++ii)
1337  tmp[ii] = *it;
1338  return get(std::min(8UL, index.size()), tmp);
1339  };
1340 
1349  T get(size_t len, int64_t* index)
1350  {
1351  assert(len <= 8);
1352  double tmp[8];
1353  for(size_t ii=0; ii < len && ii<8; ++ii)
1354  tmp[ii] = index[ii];
1355  return get(std::min(8UL, len), tmp);
1356  };
1357 
1365  T operator[](const std::vector<int64_t>& index)
1366  {
1367  assert(index.size() <= 8);
1368  double tmp[8];
1369  size_t ii=0;
1370  for(auto it = index.begin(); it != index.end() && ii<8; ++it, ++ii)
1371  tmp[ii] = *it;
1372  return get(std::min(8UL, index.size()), tmp);
1373  };
1374 
1382  T operator()(const std::vector<double>& index)
1383  {
1384  return get(index.size(), index.data());
1385  };
1386 
1394  T operator()(std::initializer_list<double> index)
1395  {
1396  assert(index.size() <= 8);
1397  double tmp[8];
1398  size_t ii=0;
1399  for(auto it = index.begin(); it != index.end() && ii<8; ++it, ++ii)
1400  tmp[ii] = *it;
1401  return get(std::min(8UL, index.size()), tmp);
1402  };
1403 
1411  T operator()(std::initializer_list<float> index)
1412  {
1413  assert(index.size() <= 8);
1414  double tmp[8];
1415  size_t ii=0;
1416  for(auto it = index.begin(); it != index.end() && ii<8; ++it, ++ii)
1417  tmp[ii] = *it;
1418  return get(std::min(8UL, index.size()), tmp);
1419  };
1420 
1426  T get(size_t len, const double* incindex)
1427  {
1428  // convert RAS to index
1429  size_t ndim = this->parent->ndim();
1430  assert(ndim <= MAXDIM);
1431  double cindex[MAXDIM];
1432  int64_t index[MAXDIM];
1433 
1434  // convert RAS to index
1435  if(m_ras) {
1436  auto tmp = dPtrCast<const MRImage>(this->parent);
1437  tmp->pointToIndex(len, incindex, cindex);
1438  } else {
1439  for(size_t dd=0; dd<ndim; dd++) {
1440  if(dd < len)
1441  cindex[dd] = incindex[dd];
1442  else
1443  cindex[dd] = 0;
1444  }
1445  }
1446 
1447  // initialize variables
1448  const size_t* dim = this->parent->dim();
1449 
1450  // round values from cindex
1451  if(m_boundmethod == ZEROFLUX) {
1452  // clamp
1453  for(size_t dd=0; dd<ndim; dd++) {
1454  double C = dd < len ? cindex[dd] : 0;
1455  index[dd] = clamp<int64_t>(0, dim[dd]-1, round(C));
1456  }
1457  } else if(m_boundmethod == WRAP) {
1458  // wrap
1459  for(size_t dd=0; dd<ndim; dd++) {
1460  double C = dd < len ? cindex[dd] : 0;
1461  index[dd] = wrap<int64_t>(0, dim[dd]-1, round(C));
1462  }
1463  } else {
1464  for(size_t dd=0; dd<ndim; dd++) {
1465  double C = dd < len ? cindex[dd] : 0;
1466  index[dd] = round(C);
1467  if(index[dd] < 0 || index[dd] > dim[dd]-1)
1468  return (T)0;
1469  }
1470  }
1471 
1472  auto ptr = this->parent->__getAddr(ndim, index);
1473  assert(ptr >= this->parent->__getAddr(0) &&
1474  ptr < this->parent->__getAddr(this->parent->elements()));
1475  return this->castget(ptr);
1476  }
1477 
1483  T get(const vector<double>& cindex)
1484  {
1485  return get(cindex.size(), cindex.data());
1486  }
1487 
1488  BoundaryConditionT m_boundmethod;
1489 
1494  bool m_ras;
1495 
1496 private:
1498  // Hide Unused Functions From Parent
1500 
1506  T operator[](int64_t index)
1507  {
1508  auto ptr = this->parent->__getAddr(index);
1509  assert(ptr >= this->parent->__getAddr(0) &&
1510  ptr < this->parent->__getAddr(this->parent->elements()));
1511  return castget(ptr);
1512  };
1513 };
1514 
1525 template<typename T>
1526 class NNInterp3DView : public Vector3DConstView<T>
1527 {
1528 public:
1529  NNInterp3DView(std::shared_ptr<NDArray> in,
1530  BoundaryConditionT bound = ZEROFLUX)
1531  : Vector3DConstView<T>(in), m_boundmethod(bound), m_ras(false)
1532  { };
1533 
1534  NNInterp3DView() : m_boundmethod(ZEROFLUX), m_ras(false) {} ;
1535 
1541  T operator()(double x=0, double y=0, double z=0, int64_t t=0)
1542  {
1543  return get(x,y,z,t);
1544  };
1545 
1551  T get(double x=0, double y=0, double z=0, int64_t t=0)
1552  {
1553  // convert RAS to cindex
1554  if(m_ras) {
1555  double cindex[3] = {x,y,z};
1556  auto tmp = dPtrCast<const MRImage>(this->parent);
1557  tmp->pointToIndex(3, cindex, cindex);
1558  x = cindex[0];
1559  y = cindex[1];
1560  z = cindex[2];
1561  }
1562 
1563 
1564  // interpolate
1565  int64_t i = round(x);
1566  int64_t j = round(y);
1567  int64_t k = round(z);
1568  size_t xdim = this->parent->dim(0);
1569  size_t ydim = this->parent->ndim() > 1 ? this->parent->dim(1) : 1;
1570  size_t zdim = this->parent->ndim() > 2 ? this->parent->dim(2) : 1;
1571  size_t tdim = this->parent->tlen();
1572 
1573  bool xout = (i < 0 || i >= xdim);
1574  bool yout = (j < 0 || j >= ydim);
1575  bool zout = (k < 0 || k >= zdim);
1576  bool tout = (t < 0 || t >= tdim);
1577 
1578  if(xout || yout || zout || tout) {
1579  // outside = true;
1580  switch(m_boundmethod) {
1581  case ZEROFLUX:
1582  i = clamp<int64_t>(0, xdim-1, i);
1583  j = clamp<int64_t>(0, ydim-1, j);
1584  k = clamp<int64_t>(0, zdim-1, k);
1585  t = clamp<int64_t>(0, tdim-1, t);
1586  break;
1587  case WRAP:
1588  i = wrap<int64_t>(0, xdim-1, i);
1589  j = wrap<int64_t>(0, ydim-1, j);
1590  k = wrap<int64_t>(0, zdim-1, k);
1591  t = wrap<int64_t>(0, tdim-1, t);
1592  break;
1593  case CONSTZERO:
1594  default:
1595  return 0;
1596  break;
1597  }
1598  }
1599 
1600  auto ptr = this->parent->__getAddr(i,j,k,t);
1601  assert(ptr >= this->parent->__getAddr(0) &&
1602  ptr < this->parent->__getAddr(this->parent->elements()));
1603  return this->castget(ptr);
1604  };
1605 
1611  T operator()(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
1612  {
1613  return get((double)x,(double)y,(double)z,t);
1614  };
1615 
1621  T get(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
1622  {
1623  return get((double)x,(double)y,(double)z,t);
1624  };
1625 
1626 
1627  BoundaryConditionT m_boundmethod;
1628 
1633  bool m_ras;
1634 
1635 private:
1636 
1637 };
1638 
1640 // Lanczos
1642 
1653 template<typename T>
1654 class LanczosInterpNDView : public NDConstView<T>
1655 {
1656 public:
1657  LanczosInterpNDView(std::shared_ptr<const NDArray> in,
1658  BoundaryConditionT bound = ZEROFLUX)
1659  : NDConstView<T>(in), m_boundmethod(bound), m_ras(false),
1660  m_radius(2)
1661  { };
1662 
1663  LanczosInterpNDView() : m_boundmethod(ZEROFLUX), m_ras(false), m_radius(2) {} ;
1664 
1665  void setRadius(size_t rad) { m_radius = rad; };
1666  size_t getRadius() { return m_radius; };
1667 
1673  T operator()(double x=0, double y=0, double z=0, double t=0, double u=0,
1674  double v=0, double w=0)
1675  {
1676  double tmp[8] = {x,y,z,t,u,v,w};
1677  return get(8, tmp);
1678  };
1679 
1687  T get(const std::vector<int64_t>& index)
1688  {
1689  assert(index.size() <= 8);
1690  double tmp[8];
1691  size_t ii=0;
1692  for(auto it = index.begin(); it != index.end() && ii<8; ++it, ++ii)
1693  tmp[ii] = *it;
1694  return get(std::min(8UL, index.size()), tmp);
1695  };
1696 
1705  T get(size_t len, int64_t* index)
1706  {
1707  assert(len <= 8);
1708  double tmp[8];
1709  for(size_t ii=0; ii<len && ii<8; ii++)
1710  tmp[ii] = index[ii];
1711  return get(std::min(8UL, len), tmp);
1712  };
1713 
1721  T operator[](const std::vector<int64_t>& index)
1722  {
1723  assert(index.size() <= 8);
1724  double tmp[8];
1725  size_t ii=0;
1726  for(auto it = index.begin(); it != index.end() && ii<8; ++it, ++ii)
1727  tmp[ii] = *it;
1728  return get(std::min(8UL, index.size()), tmp);
1729  };
1730 
1736  T get(const std::vector<double>& incoord)
1737  {
1738  return get(incoord.size(), incoord.data());
1739  }
1740 
1746  T get(size_t len, const double* incoord)
1747  {
1748  // figure out size of dimensions in parent
1749  const size_t* dim = this->parent->dim();
1750  size_t ndim = this->parent->ndim();
1751  assert(ndim < MAXDIM);
1752  int64_t index[MAXDIM];
1753  double cindex[MAXDIM];
1754 
1755  // convert RAS to index
1756  if(m_ras) {
1757  auto tmp = dPtrCast<const MRImage>(this->parent);
1758  tmp->pointToIndex(len, incoord, cindex);
1759  } else {
1760  for(size_t dd=0; dd<ndim; dd++) {
1761  if(dd < len)
1762  cindex[dd] = incoord[dd];
1763  else
1764  cindex[dd] = 0;
1765  }
1766  }
1767 
1768  const int KPOINTS = 1+m_radius*2;
1769 
1770  // 1D version of the weights and indices
1771  double karray[MAXDIM][KPOINTS];
1772  int64_t indarray[MAXDIM][KPOINTS];
1773 
1774  for(int dd = 0; dd < ndim; dd++) {
1775  for(int64_t ii=-m_radius; ii<=m_radius; ii++){
1776  int64_t i = round(cindex[dd])+ii;
1777  indarray[dd][ii+m_radius] = i;
1778  karray[dd][ii+m_radius] = lanczosKern(i-cindex[dd], m_radius);
1779  }
1780  }
1781 
1782  // initialize variables
1783  Counter<int, MAXDIM> count(ndim);
1784  for(size_t dd=0; dd<ndim; dd++)
1785  count.sz[dd] = 1+m_radius*2;
1786 
1787  // compute weighted pixval by iterating over neighbors, which are
1788  // combinations of KPOINTS
1789  T pixval = 0;
1790  do {
1791  double weight = 1;
1792  bool iioutside = false;
1793 
1794  //set index
1795  for(int dd = 0; dd < ndim; dd++) {
1796  index[dd] = indarray[dd][count.pos[dd]];
1797  weight *= karray[dd][count.pos[dd]];
1798  iioutside = iioutside || index[dd] < 0 || index[dd] >= dim[dd];
1799  }
1800 
1801  if(iioutside) {
1802  if(m_boundmethod == ZEROFLUX) {
1803  // clamp
1804  for(size_t dd=0; dd<ndim; dd++)
1805  index[dd] = clamp<int64_t>(0, dim[dd]-1, index[dd]);
1806  } else if(m_boundmethod == WRAP) {
1807  // wrap
1808  for(size_t dd=0; dd<ndim; dd++)
1809  index[dd] = wrap<int64_t>(0, dim[dd]-1, index[dd]);
1810  } else { // zero outside
1811  // set wieght to zero, then just clamp
1812  weight = 0;
1813  for(size_t dd=0; dd<ndim; dd++)
1814  index[dd] = clamp<int64_t>(0, dim[dd]-1, index[dd]);
1815  }
1816  }
1817 
1818  auto ptr = this->parent->__getAddr(ndim, index);
1819  assert(ptr >= this->parent->__getAddr(0) &&
1820  ptr < this->parent->__getAddr(this->parent->elements()));
1821  T v = this->castget(ptr);
1822  pixval += weight*v;
1823  } while(count.advance());
1824 
1825  return pixval;
1826  }
1827 
1828  BoundaryConditionT m_boundmethod;
1829 
1834  bool m_ras;
1835 
1836 protected:
1837 
1843  T operator[](int64_t index)
1844  {
1845  auto ptr = this->parent->__getAddr(index);
1846  assert(ptr >= this->parent->__getAddr(0) &&
1847  ptr < this->parent->__getAddr(this->parent->elements()));
1848  return castget(ptr);
1849  };
1850 
1851  int64_t m_radius;
1852 };
1853 
1854 
1865 template<typename T>
1866 class LanczosInterp3DView : public Vector3DConstView<T>
1867 {
1868 public:
1869  LanczosInterp3DView(std::shared_ptr<const NDArray> in,
1870  BoundaryConditionT bound = ZEROFLUX)
1871  : Vector3DConstView<T>(in), m_boundmethod(bound), m_ras(false),
1872  m_radius(2)
1873  { };
1874 
1875  LanczosInterp3DView() : m_boundmethod(ZEROFLUX), m_ras(false), m_radius(2) {} ;
1876 
1877  void setRadius(size_t rad) { m_radius = rad; };
1878  size_t getRadius() { return m_radius; };
1879 
1885  T operator()(double x=0, double y=0, double z=0, int64_t t=0)
1886  {
1887  return get(x,y,z,t);
1888  };
1889 
1895  T operator()(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
1896  {
1897  return get((double)x,(double)y,(double)z,t);
1898  };
1899 
1905  T get(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
1906  {
1907  return get((double)x,(double)y,(double)z,t);
1908  };
1909 
1910 
1916  T get(double x=0, double y=0, double z=0, int64_t t=0)
1917  {
1918  // figure out size of dimensions in parent
1919  size_t dim[4];
1920  double cindex[3] = {x,y,z};
1921  int64_t index[3];
1922  const int ndim = 3;
1923  dim[0] = this->parent->dim(0);
1924  dim[1] = this->parent->ndim() > 1 ? this->parent->dim(1) : 1;
1925  dim[2] = this->parent->ndim() > 2 ? this->parent->dim(2) : 1;
1926  dim[3] = this->parent->tlen();
1927 
1928  // deal with t being outside bounds
1929  if(t < 0 || t >= dim[3]) {
1930  if(m_boundmethod == ZEROFLUX) {
1931  // clamp
1932  t = clamp<int64_t>(0, dim[3]-1, t);
1933  } else if(m_boundmethod == WRAP) {
1934  // wrap
1935  t = wrap<int64_t>(0, dim[3]-1, t);
1936  } else {
1937  return 0;
1938  }
1939  }
1940 
1941  // convert RAS to cindex
1942  if(m_ras) {
1943  auto tmp = dPtrCast<const MRImage>(this->parent);
1944  tmp->pointToIndex(3, cindex, cindex);
1945  }
1946 
1947  const int KPOINTS = 1+m_radius*2;
1948  const int DIM = 3;
1949 
1950  // 1D version of the weights and indices
1951  double karray[DIM][KPOINTS];
1952  int64_t indarray[DIM][KPOINTS];
1953 
1954  for(int dd = 0; dd < DIM; dd++) {
1955  for(int64_t ii=-m_radius; ii<=m_radius; ii++){
1956  int64_t i = round(cindex[dd])+ii;
1957  indarray[dd][ii+m_radius] = i;
1958  karray[dd][ii+m_radius] = lanczosKern(i-cindex[dd], m_radius);
1959  }
1960  }
1961 
1962  T pixval = 0;
1963  Counter<int, MAXDIM> count(ndim);
1964  for(size_t dd=0; dd<ndim; dd++)
1965  count.sz[dd] = 1+m_radius*2;
1966 
1967  // compute weighted pixval by iterating over neighbors, which are
1968  // combinations of KPOINTS
1969  do {
1970  double weight = 1;
1971  bool iioutside = false;
1972 
1973  //set index
1974  for(int dd = 0; dd < ndim; dd++) {
1975  index[dd] = indarray[dd][count.pos[dd]];
1976  weight *= karray[dd][count.pos[dd]];
1977  iioutside = iioutside || index[dd] < 0 || index[dd] >= dim[dd];
1978  }
1979 
1980  // if the current point maps outside, then we need to deal with it
1981  if(iioutside) {
1982  if(m_boundmethod == ZEROFLUX) {
1983  // clamp
1984  for(size_t dd=0; dd<ndim; dd++)
1985  index[dd] = clamp<int64_t>(0, dim[dd]-1, index[dd]);
1986  } else if(m_boundmethod == WRAP) {
1987  // wrap
1988  for(size_t dd=0; dd<ndim; dd++)
1989  index[dd] = wrap<int64_t>(0, dim[dd]-1, index[dd]);
1990  } else {
1991  // set wieght to zero, then just clamp
1992  weight = 0;
1993  for(size_t dd=0; dd<ndim; dd++)
1994  index[dd] = clamp<int64_t>(0, dim[dd]-1, index[dd]);
1995  }
1996  }
1997 
1998  auto ptr = this->parent->__getAddr(index[0], index[1],index[2],t);
1999  assert(ptr >= this->parent->__getAddr(0) &&
2000  ptr < this->parent->__getAddr(this->parent->elements()));
2001  T v = this->castget(ptr);
2002  pixval += weight*v;
2003  } while(count.advance());
2004 
2005  return pixval;
2006  }
2007 
2008  BoundaryConditionT m_boundmethod;
2009 
2014  bool m_ras;
2015 
2016 protected:
2017 
2023  T operator[](int64_t i) { (void)(i); return T(); };
2024  T get(const std::vector<int64_t>& i) { (void)(i); return T(); };
2025  T operator[](const std::vector<int64_t>& i) { (void)(i); return T(); };
2026  int64_t m_radius;
2027 };
2028 
2029 
2030 
2031 /****************************************************************************
2032  * BSpline
2033  ***************************************************************************/
2034 
2048 template<typename T>
2049 class BSplineView : public NDView<T>
2050 {
2051 public:
2052  BSplineView(std::shared_ptr<MRImage> params,
2053  BoundaryConditionT bound = ZEROFLUX)
2054  : NDView<T>(params), m_boundmethod(bound), m_ras(false)
2055  { };
2056 
2057  BSplineView() : m_boundmethod(ZEROFLUX), m_ras(false){} ;
2058 
2059  void createOverlay(ptr<const MRImage> overlay, double bspace)
2060  {
2061  size_t ndim = overlay->ndim();
2062  assert(ndim <= MAXDIM);
2063  VectorXd spacing(overlay->ndim());
2064  VectorXd origin(overlay->ndim());
2065  size_t osize[MAXDIM];
2066 
2067  // get spacing and size
2068  for(size_t dd=0; dd<ndim; ++dd) {
2069  osize[dd] = 4+ceil(overlay->dim(dd)*overlay->spacing(dd)/bspace);
2070  spacing[dd] = bspace;
2071  }
2072 
2073  auto params = dPtrCast<MRImage>(overlay->createAnother(ndim, osize,
2074  FLOAT64));
2075  this->parent = params;
2076  params->setDirection(overlay->getDirection(), false);
2077  params->setSpacing(spacing, false);
2078 
2079  // compute center of input
2080  VectorXd indc(ndim); // center index
2081  for(size_t dd=0; dd<ndim; dd++)
2082  indc[dd] = (overlay->dim(dd)-1.)/2.;
2083  VectorXd ptc(ndim); // point center
2084  overlay->indexToPoint(ndim, indc.array().data(), ptc.array().data());
2085 
2086  // compute origin from center index (x_c) and center of input (c):
2087  // o = c-R(sx_c)
2088  for(size_t dd=0; dd<ndim; dd++)
2089  indc[dd] = (osize[dd]-1.)/2.;
2090  origin = ptc - overlay->getDirection()*(spacing.asDiagonal()*indc);
2091  params->setOrigin(origin, false);
2092 
2093  this->setArray(params);
2094  };
2095 
2108  bool get(size_t len, const double* point, int dir, double& val, double& dval)
2109  {
2110  assert(this->parent);
2111  // initialize variables
2112  int ndim = this->parent->ndim();
2113  assert(ndim <= MAXDIM);
2114  const size_t* dim = this->parent->dim();
2115  int64_t index[MAXDIM];
2116  double cindex[MAXDIM];
2117 
2118  // convert RAS to index, or just copy into local array
2119  if(m_ras) {
2120  auto tmp = dPtrCast<const MRImage>(this->parent);
2121  tmp->pointToIndex(len, point, cindex);
2122  } else {
2123  for(size_t dd=0; dd<ndim; dd++) {
2124  if(dd < len)
2125  cindex[dd] = point[dd];
2126  else
2127  cindex[dd] = 0;
2128  }
2129  }
2130 
2131  Counter<> count(ndim);
2132  for(size_t dd=0; dd<ndim; dd++)
2133  count.sz[dd] = 5;
2134 
2135  dval = 0;
2136  val = 0;
2137  bool border = false;
2138  do {
2139  double weight = 1;
2140  double dweight = 1;
2141  bool iioutside = false;
2142 
2143  //set index
2144  for(int dd = 0; dd < ndim; dd++) {
2145  index[dd] = round(cindex[dd]) + count.pos[dd] - 2l;
2146  weight *= B3kern(index[dd] - cindex[dd]);
2147  if(dd == dir)
2148  dweight *= -dB3kern(index[dd] - cindex[dd])/
2149  getParams()->spacing(dd);
2150  else
2151  dweight *= B3kern(index[dd] - cindex[dd]);
2152  iioutside = iioutside || index[dd] < 0 || index[dd] >= dim[dd];
2153  }
2154 
2155  // if the current point maps outside, then we need to deal with it
2156  if(iioutside) {
2157  if(m_boundmethod == ZEROFLUX) {
2158  // clamp
2159  for(size_t dd=0; dd<ndim; dd++)
2160  index[dd] = clamp<int64_t>(0, dim[dd]-1, index[dd]);
2161  } else if(m_boundmethod == WRAP) {
2162  // wrap
2163  for(size_t dd=0; dd<ndim; dd++)
2164  index[dd] = wrap<int64_t>(0, dim[dd]-1, index[dd]);
2165  } else {
2166  // set wieght to zero, then just clamp
2167  weight = 0;
2168  for(size_t dd=0; dd<ndim; dd++)
2169  index[dd] = clamp<int64_t>(0, dim[dd]-1, index[dd]);
2170  }
2171  }
2172 
2173  // Update border with outside
2174  border |= iioutside;
2175 
2176  // Compute Values
2177  auto ptr = this->parent->__getAddr(ndim, index);
2178  assert(ptr >= this->parent->__getAddr(0) &&
2179  ptr < this->parent->__getAddr(this->parent->elements()));
2180  T v = this->castget(ptr);
2181  dval += dweight*v;
2182  val += weight*v;
2183  } while(count.advance());
2184 
2185  return border;
2186  };
2187 
2198  double get(size_t len, double* point)
2199  {
2200  assert(this->parent);
2201  // initialize variables
2202  int ndim = this->parent->ndim();
2203  assert(ndim <= MAXDIM);
2204  const size_t* dim = this->parent->dim();
2205  int64_t index[MAXDIM];
2206  double cindex[MAXDIM];
2207 
2208  // convert RAS to index, or just copy into local array
2209  if(m_ras) {
2210  auto tmp = dPtrCast<const MRImage>(this->parent);
2211  tmp->pointToIndex(len, point, cindex);
2212  } else {
2213  for(size_t dd=0; dd<ndim; dd++) {
2214  if(dd < len)
2215  cindex[dd] = point[dd];
2216  else
2217  cindex[dd] = 0;
2218  }
2219  }
2220 
2221  Counter<> count(ndim);
2222  for(size_t dd=0; dd<ndim; dd++)
2223  count.sz[dd] = 5;
2224 
2225  double val = 0;
2226  bool border = false;
2227  do {
2228  double weight = 1;
2229  bool iioutside = false;
2230 
2231  //set index
2232  for(int dd = 0; dd < ndim; dd++) {
2233  index[dd] = round(cindex[dd]) + count.pos[dd] - 2l;
2234  weight *= B3kern(index[dd] - cindex[dd]);
2235  iioutside = iioutside || index[dd] < 0 || index[dd] >= dim[dd];
2236  }
2237 
2238  // if the current point maps outside, then we need to deal with it
2239  if(iioutside) {
2240  if(m_boundmethod == ZEROFLUX) {
2241  // clamp
2242  for(size_t dd=0; dd<ndim; dd++)
2243  index[dd] = clamp<int64_t>(0, dim[dd]-1, index[dd]);
2244  } else if(m_boundmethod == WRAP) {
2245  // wrap
2246  for(size_t dd=0; dd<ndim; dd++)
2247  index[dd] = wrap<int64_t>(0, dim[dd]-1, index[dd]);
2248  } else {
2249  // set wieght to zero, then just clamp
2250  weight = 0;
2251  for(size_t dd=0; dd<ndim; dd++)
2252  index[dd] = clamp<int64_t>(0, dim[dd]-1, index[dd]);
2253  }
2254  }
2255 
2256  border |= iioutside;
2257  auto ptr = this->parent->__getAddr(ndim, index);
2258  assert(ptr >= this->parent->__getAddr(0) &&
2259  ptr < this->parent->__getAddr(this->parent->elements()));
2260  T v = this->castget(ptr);
2261  val += weight*v;
2262  } while(count.advance());
2263 
2264  return val;
2265  };
2266 
2274  ptr<MRImage> reconstruct(ptr<const MRImage> input)
2275  {
2276  assert(this->parent);
2277  auto params = dPtrCast<MRImage>(this->parent);
2278  if(params->getDirection() != input->getDirection())
2279  return reconstructNotAligned(input);
2280  else
2281  return reconstructAligned(input);
2282  };
2283 
2284  ptr<MRImage> reconstructNotAligned(ptr<const MRImage> input)
2285  {
2286  assert(this->parent);
2287  (void)(input);
2288  throw INVALID_ARGUMENT("Not yet implemented");
2289  return NULL;
2290  };
2291 
2292  ptr<MRImage> reconstructAligned(ptr<const MRImage> input)
2293  {
2294  assert(this->parent);
2295  if(input->ndim() != this->parent->ndim()) {
2296  throw INVALID_ARGUMENT("Not sure how to deal with different "
2297  "dimensions for B-Spline right now");
2298  }
2299  auto params = getParams();
2300  auto out = dPtrCast<MRImage>(input->createAnother());
2301 
2302  // for each kernel, iterate over the points in the neighborhood
2303  size_t ndim = input->ndim();
2304  assert(ndim <= MAXDIM);
2305  double cind[MAXDIM];
2306  int64_t ind[MAXDIM];
2307  double scale[MAXDIM];
2308  int64_t center[MAXDIM];
2309  size_t winsize[MAXDIM];
2310  int64_t roistart[MAXDIM];
2311 
2312  for(size_t dd=0; dd<ndim; dd++) {
2313  winsize[dd] = ceil(5*params->spacing(dd)/out->spacing(dd));
2314  scale[dd] = out->spacing(dd)/params->spacing(dd);
2315  }
2316 
2317  // We go through each parameter, and compute the weight of the B-spline
2318  // parameter at each pixel within the range (2 indexes in parameter
2319  // space, 2*S_B/S_I indexs in pixel space)
2320  NDIter<double> oit(out);
2321  for(NDConstIter<double> bit(params); !bit.eof(); ++bit) {
2322 
2323  // get continuous index of pixel
2324  bit.index(ndim, cind);
2325  params->indexToPoint(ndim, cind, cind);
2326  out->pointToIndex(ndim, cind, cind);
2327 
2328  // construct weights / construct ROI
2329  for(size_t dd=0; dd<ndim; dd++) {
2330  center[dd] = round(cind[dd]);
2331  roistart[dd] = center[dd]-((int64_t)winsize[dd]/2);
2332  }
2333 
2334  oit.setROI(ndim, winsize, roistart);
2335  for(oit.goBegin(); !oit.eof(); ++oit) {
2336  double weight = 1;
2337  oit.index(ndim, ind);
2338  for(size_t dd=0; dd<ndim; dd++) {
2339  double dist = (ind[dd] - cind[dd])*scale[dd];
2340  weight *= B3kern(dist);
2341  }
2342  oit.set(oit.get()+(*bit)*weight);
2343  }
2344  }
2345 
2346  return out;
2347  };
2348 
2349  double thinPlateEnergy()
2350  {
2351  using std::pow;
2352 
2353  auto params = getParams();
2354  assert(params->ndim() == 3);
2355 
2356  double dphi2_dx2 = 0;
2357  double dphi2_dy2 = 0;
2358  double dphi2_dz2 = 0;
2359  double dphi2_dxz = 0;
2360  double dphi2_dxy = 0;
2361  double dphi2_dyz = 0;
2362 
2363  // We use NN interpolator because it is fast and handles boundary
2364  // conditions
2365  NNInterpNDView<double> dvw(params);
2366  double ind[3];
2367 
2368  // Create a counter to iterate over [0,4] ie [-2,2] in 6 directions
2369  Counter<int, 3> counter(3);
2370  for(size_t dd=0; dd<3;dd++)
2371  counter.sz[dd] = 9;
2372 
2373  //integrate over all the knots
2374  Counter<int64_t, 3> it(3, (int64_t*)params->dim());
2375  do {
2376  double phi_ijk = dvw.get(3, it.pos);
2377  do {
2378  for(size_t dd=0; dd<3; dd++)
2379  ind[dd] = it.pos[dd] + counter.pos[dd] - 4;
2380  double phi_lmn = dvw.get(3, ind);
2381 
2382  double xv = U200[counter.pos[0]][counter.pos[1]][counter.pos[2]];
2383  double yv = U200[counter.pos[1]][counter.pos[0]][counter.pos[2]];
2384  double zv = U200[counter.pos[2]][counter.pos[1]][counter.pos[0]];
2385  double xyv = U110[counter.pos[0]][counter.pos[1]][counter.pos[2]];
2386  double xzv = U110[counter.pos[0]][counter.pos[2]][counter.pos[1]];
2387  double yzv = U110[counter.pos[1]][counter.pos[2]][counter.pos[0]];
2388 
2389  dphi2_dx2 += phi_ijk*phi_lmn*xv;
2390  dphi2_dy2 += phi_ijk*phi_lmn*yv;
2391  dphi2_dz2 += phi_ijk*phi_lmn*zv;
2392  dphi2_dxy += phi_ijk*phi_lmn*xyv;
2393  dphi2_dyz += phi_ijk*phi_lmn*yzv;
2394  dphi2_dxz += phi_ijk*phi_lmn*xzv;
2395  } while(counter.advance());
2396  } while(it.advance());
2397 
2398  double v = dphi2_dx2/pow(params->spacing(0),4) +
2399  dphi2_dy2/pow(params->spacing(1),4) +
2400  dphi2_dz2/pow(params->spacing(2),4) +
2401  2*dphi2_dxy/(pow(params->spacing(0),2)*pow(params->spacing(1),2)) +
2402  2*dphi2_dxz/(pow(params->spacing(0),2)*pow(params->spacing(2),2)) +
2403  2*dphi2_dyz/(pow(params->spacing(1),2)*pow(params->spacing(2),2));
2404 
2405  return v;
2406  };
2407 
2408  double thinPlateEnergy(size_t len, double* grad)
2409  {
2410  using std::pow;
2411 
2412  auto params = getParams();
2413  if(len != getParams()->elements())
2414  throw INVALID_ARGUMENT("Incorrect length of grad array");
2415 
2416  assert(params->ndim() == 3);
2417 
2418  double dphi2_dx2 = 0;
2419  double dphi2_dy2 = 0;
2420  double dphi2_dz2 = 0;
2421  double dphi2_dxz = 0;
2422  double dphi2_dxy = 0;
2423  double dphi2_dyz = 0;
2424  double reg = 0;
2425 
2426  NNInterpNDView<double> dvw(params);
2427  double ind[3];
2428 
2429  // Create a counter to iterate over [0,4] ie [-2,2] in 6 directions
2430  Counter<int, 3> counter(3);
2431  for(size_t dd=0; dd<3;dd++)
2432  counter.sz[dd] = 9;
2433 
2434  //integrate over all the knots
2435  int ii = 0;
2436  Counter<int64_t, 3> it(3, (int64_t*)params->dim());
2437  ii = 0;
2438  do {
2439  // iterate through abc
2440  double tmp_dphi2_dx2 = 0;
2441  double tmp_dphi2_dy2 = 0;
2442  double tmp_dphi2_dz2 = 0;
2443  double tmp_dphi2_dxz = 0;
2444  double tmp_dphi2_dxy = 0;
2445  double tmp_dphi2_dyz = 0;
2446  double phi_ijk = dvw.get(3, it.pos);
2447 
2448  // for all the neighbors of abc, get
2449  do {
2450  for(size_t dd=0; dd<3; dd++)
2451  ind[dd] = it.pos[dd] + (counter.pos[dd] - 4);
2452  double phi_lmn = dvw.get(3, ind);
2453 
2454  for(size_t dd=0; dd<3; dd++)
2455  ind[dd] = it.pos[dd] - (counter.pos[dd] - 4);
2456  double phi_abc = dvw.get(3, ind);
2457 
2458  double xv = U200[counter.pos[0]][counter.pos[1]][counter.pos[2]];
2459  double yv = U200[counter.pos[1]][counter.pos[0]][counter.pos[2]];
2460  double zv = U200[counter.pos[2]][counter.pos[1]][counter.pos[0]];
2461  double xyv = U110[counter.pos[0]][counter.pos[1]][counter.pos[2]];
2462  double xzv = U110[counter.pos[0]][counter.pos[2]][counter.pos[1]];
2463  double yzv = U110[counter.pos[1]][counter.pos[2]][counter.pos[0]];
2464 
2465  tmp_dphi2_dx2 += phi_abc*xv + phi_lmn*xv;
2466  tmp_dphi2_dy2 += phi_abc*yv + phi_lmn*yv;
2467  tmp_dphi2_dz2 += phi_abc*zv + phi_lmn*zv;
2468  tmp_dphi2_dxy += phi_abc*xyv + phi_lmn*xyv;
2469  tmp_dphi2_dyz += phi_abc*yzv + phi_lmn*yzv;
2470  tmp_dphi2_dxz += phi_abc*xzv + phi_lmn*xzv;
2471 
2472  dphi2_dx2 += phi_ijk*phi_lmn*xv;
2473  dphi2_dy2 += phi_ijk*phi_lmn*yv;
2474  dphi2_dz2 += phi_ijk*phi_lmn*zv;
2475  dphi2_dxy += phi_ijk*phi_lmn*xyv;
2476  dphi2_dyz += phi_ijk*phi_lmn*yzv;
2477  dphi2_dxz += phi_ijk*phi_lmn*xzv;
2478 
2479  } while(counter.advance());
2480 
2481  reg = tmp_dphi2_dx2/pow(params->spacing(0),4) +
2482  tmp_dphi2_dy2/pow(params->spacing(1),4) +
2483  tmp_dphi2_dz2/pow(params->spacing(2),4) +
2484  2*tmp_dphi2_dxy/(pow(params->spacing(0),2)*pow(params->spacing(1),2)) +
2485  2*tmp_dphi2_dxz/(pow(params->spacing(0),2)*pow(params->spacing(2),2)) +
2486  2*tmp_dphi2_dyz/(pow(params->spacing(1),2)*pow(params->spacing(2),2));
2487  grad[ii++] = reg;
2488  } while(it.advance());
2489 
2490  double v = dphi2_dx2/pow(params->spacing(0),4) +
2491  dphi2_dy2/pow(params->spacing(1),4) +
2492  dphi2_dz2/pow(params->spacing(2),4) +
2493  2*dphi2_dxy/(pow(params->spacing(0),2)*pow(params->spacing(1),2)) +
2494  2*dphi2_dxz/(pow(params->spacing(0),2)*pow(params->spacing(2),2)) +
2495  2*dphi2_dyz/(pow(params->spacing(1),2)*pow(params->spacing(2),2));
2496 
2497  return v;
2498  };
2499 
2508  double jacobianDet(int dir)
2509  {
2510  using std::pow;
2511 
2512  auto params = getParams();
2513  assert(params->ndim() == 3);
2514 
2515  // We use NN interpolator because it is fast and handles boundary
2516  // conditions
2517  NNInterpNDView<double> dvw(params);
2518  double ind[3];
2519  double reg = 0;
2520 
2521  // Create a counter to iterate over [0,4] ie [-2,2] in 6 directions
2522  Counter<int, 3> counter(3);
2523  for(size_t dd=0; dd<3;dd++)
2524  counter.sz[dd] = 9;
2525 
2526  //integrate over all the knots
2527  Counter<int64_t, 3> it(3, (int64_t*)params->dim());
2528  do {
2529  double phi_ijk = dvw.get(3, it.pos);
2530  do {
2531  for(size_t dd=0; dd<3; dd++)
2532  ind[dd] = it.pos[dd] + counter.pos[dd] - 4;
2533  double phi_lmn = dvw.get(3, ind);
2534 
2535  switch(dir) {
2536  case 0:
2537  reg += U100[counter.pos[0]][counter.pos[1]][counter.pos[2]]*
2538  phi_ijk*phi_lmn;
2539  break;
2540  case 1:
2541  reg += U100[counter.pos[1]][counter.pos[0]][counter.pos[2]]*
2542  phi_ijk*phi_lmn;
2543  break;
2544  case 2:
2545  reg += U100[counter.pos[2]][counter.pos[1]][counter.pos[0]]*
2546  phi_ijk*phi_lmn;
2547  break;
2548  }
2549  } while(counter.advance());
2550  } while(it.advance());
2551 
2552  return reg/pow(params->spacing(dir),2);
2553  };
2554 
2565  double jacobianDet(int dir, size_t len, double* grad)
2566  {
2567  using std::pow;
2568 
2569  auto params = getParams();
2570  assert(params->ndim() == 3);
2571  if(len != getParams()->elements())
2572  throw INVALID_ARGUMENT("Incorrect length of grad array");
2573 
2574  // We use NN interpolator because it is fast and handles boundary
2575  // conditions
2576  NNInterpNDView<double> dvw(params);
2577  double ind[3];
2578 
2579  // Create a counter to iterate over [0,4] ie [-2,2] in 6 directions
2580  Counter<int, 3> counter(3);
2581  for(size_t dd=0; dd<3;dd++)
2582  counter.sz[dd] = 9;
2583 
2584  //integrate over all the knots
2585  Counter<int64_t, 3> it(3, (int64_t*)params->dim());
2586 
2587  double reg = 0;
2588  size_t ii=0;
2589  do {
2590  double dreg = 0;
2591  double phi_ijk = dvw.get(3, it.pos);
2592  do {
2593  for(size_t dd=0; dd<3; dd++)
2594  ind[dd] = it.pos[dd] + (counter.pos[dd] - 4);
2595  double phi_lmn = dvw.get(3, ind);
2596 
2597  for(size_t dd=0; dd<3; dd++)
2598  ind[dd] = it.pos[dd] - (counter.pos[dd] - 4);
2599  double phi_abc = dvw.get(3, ind);
2600 
2601  double uval = 0;
2602  switch(dir) {
2603  case 0:
2604  uval = U100[counter.pos[0]][counter.pos[1]][counter.pos[2]];
2605  break;
2606  case 1:
2607  uval = U100[counter.pos[1]][counter.pos[0]][counter.pos[2]];
2608  break;
2609  case 2:
2610  uval = U100[counter.pos[2]][counter.pos[1]][counter.pos[0]];
2611  break;
2612  }
2613 
2614  dreg += (phi_abc + phi_lmn)*uval;
2615  reg += uval*phi_ijk*phi_lmn;
2616  } while(counter.advance());
2617 
2618  grad[ii++] = dreg/pow(params->spacing(dir),2);
2619  } while(it.advance());
2620 
2621  double v = reg/pow(params->spacing(dir),2);
2622  return v;
2623  };
2624 
2630  ptr<MRImage> getParams() { return dPtrCast<MRImage>(this->parent); };
2631 
2637  BoundaryConditionT m_boundmethod;
2638 
2643  bool m_ras;
2644 
2645  double vConv[5][5] = {
2646  {0.000031001984126984125,0.0010788690476190477,0.0013950892857142857,0.0000992063492063492,0.},
2647  {0.0010788690476190477,0.058643353174603174,0.11707589285714286,0.02101934523809524,0.0000992063492063492},
2648  {0.0013950892857142857,0.11707589285714286,0.362016369047619,0.11707589285714286,0.0013950892857142857},
2649  {0.0000992063492063492,0.02101934523809524,0.11707589285714286,0.058643353174603174,0.0010788690476190477},
2650  {0.,0.0000992063492063492,0.0013950892857142857,0.0010788690476190477,0.000031001984126984125}};
2651  double dvConv[5][5] = {
2652  {0.0015625,0.013541666666666667,-0.0109375,-0.004166666666666667,0.},
2653  {0.013541666666666667,0.24479166666666666,-0.07604166666666666,-0.178125,-0.004166666666666667},
2654  {-0.0109375,-0.07604166666666666,0.17395833333333333,-0.07604166666666666,-0.0109375},
2655  {-0.004166666666666667,-0.178125,-0.07604166666666666,0.24479166666666666,0.013541666666666667},
2656  {0.,-0.004166666666666667,-0.0109375,0.013541666666666667,0.0015625}};
2657  double ddvConv[5][5] = {
2658  {0.041666666666666664,0.,-0.125,0.08333333333333333,0.},
2659  {0.,0.4166666666666667,-0.75,0.25,0.08333333333333333},
2660  {-0.125,-0.75,1.75,-0.75,-0.125},
2661  {0.08333333333333333,0.25,-0.75,0.4166666666666667,0.},
2662  {0.,0.08333333333333333,-0.125,0.,0.041666666666666664}};
2663 
2664  double U200[9][9][9] =
2665 {{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2666  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2667  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2668  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2669  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2670  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2671  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2672  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2673  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
2674  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2675  {0.0,
2676  6.561266481901402e-09,
2677  7.873519778281682e-07,
2678  7.814468379944569e-06,
2679  1.5852019820273786e-05,
2680  7.814468379944569e-06,
2681  7.873519778281682e-07,
2682  6.561266481901402e-09,
2683  0.0},
2684  {0.0,
2685  7.873519778281682e-07,
2686  9.44822373393802e-05,
2687  0.0009377362055933485,
2688  0.0019022423784328549,
2689  0.0009377362055933485,
2690  9.44822373393802e-05,
2691  7.873519778281682e-07,
2692  0.0},
2693  {0.0,
2694  7.81446837994457e-06,
2695  0.0009377362055933487,
2696  0.00930703184051398,
2697  0.018879755605946087,
2698  0.00930703184051398,
2699  0.0009377362055933487,
2700  7.81446837994457e-06,
2701  0.0},
2702  {0.0,
2703  1.585201982027379e-05,
2704  0.0019022423784328547,
2705  0.018879755605946083,
2706  0.03829847988578147,
2707  0.018879755605946083,
2708  0.0019022423784328547,
2709  1.585201982027379e-05,
2710  0.0},
2711  {0.0,
2712  7.81446837994457e-06,
2713  0.0009377362055933487,
2714  0.00930703184051398,
2715  0.018879755605946087,
2716  0.00930703184051398,
2717  0.0009377362055933487,
2718  7.81446837994457e-06,
2719  0.0},
2720  {0.0,
2721  7.873519778281682e-07,
2722  9.44822373393802e-05,
2723  0.0009377362055933485,
2724  0.0019022423784328549,
2725  0.0009377362055933485,
2726  9.44822373393802e-05,
2727  7.873519778281682e-07,
2728  0.0},
2729  {0.0,
2730  6.561266481901402e-09,
2731  7.873519778281682e-07,
2732  7.814468379944569e-06,
2733  1.5852019820273786e-05,
2734  7.814468379944569e-06,
2735  7.873519778281682e-07,
2736  6.561266481901402e-09,
2737  0.0},
2738  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
2739  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2740  {0.0,
2741  0.0,
2742  8.933550615572699e-23,
2743  6.617444900424221e-22,
2744  -1.4097225589434974e-21,
2745  6.617444900424221e-22,
2746  8.933550615572699e-23,
2747  0.0,
2748  0.0},
2749  {0.0,
2750  -1.3896634290890865e-22,
2751  -2.948733447629033e-20,
2752  3.095640724418451e-19,
2753  -1.2023219095968517e-18,
2754  3.095640724418451e-19,
2755  -2.948733447629033e-20,
2756  -1.3896634290890865e-22,
2757  0.0},
2758  {0.0,
2759  -1.2176098616780567e-21,
2760  -3.814295240604521e-20,
2761  -2.910855193019005e-18,
2762  -5.991280757250155e-18,
2763  -2.910855193019005e-18,
2764  -3.814295240604521e-20,
2765  -1.2176098616780567e-21,
2766  0.0},
2767  {0.0,
2768  -4.5135110123955955e-21,
2769  -1.0007148332605274e-18,
2770  -7.31572264940066e-18,
2771  2.069796097619653e-17,
2772  -7.31572264940066e-18,
2773  -1.0007148332605274e-18,
2774  -4.5135110123955955e-21,
2775  0.0},
2776  {0.0,
2777  -1.2176098616780567e-21,
2778  -3.814295240604521e-20,
2779  -2.910855193019005e-18,
2780  -5.991280757250155e-18,
2781  -2.910855193019005e-18,
2782  -3.814295240604521e-20,
2783  -1.2176098616780567e-21,
2784  0.0},
2785  {0.0,
2786  -1.3896634290890865e-22,
2787  -2.948733447629033e-20,
2788  3.095640724418451e-19,
2789  -1.2023219095968517e-18,
2790  3.095640724418451e-19,
2791  -2.948733447629033e-20,
2792  -1.3896634290890865e-22,
2793  0.0},
2794  {0.0,
2795  0.0,
2796  8.933550615572699e-23,
2797  6.617444900424221e-22,
2798  -1.4097225589434974e-21,
2799  6.617444900424221e-22,
2800  8.933550615572699e-23,
2801  0.0,
2802  0.0},
2803  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
2804  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2805  {0.0,
2806  -5.9051398337112626e-08,
2807  -7.086167800453517e-06,
2808  -7.033021541950112e-05,
2809  -0.00014266817838246408,
2810  -7.033021541950112e-05,
2811  -7.086167800453517e-06,
2812  -5.9051398337112626e-08,
2813  0.0},
2814  {0.0,
2815  -7.086167800453517e-06,
2816  -0.0008503401360544218,
2817  -0.00843962585034014,
2818  -0.017120181405895694,
2819  -0.00843962585034014,
2820  -0.0008503401360544218,
2821  -7.086167800453517e-06,
2822  0.0},
2823  {0.0,
2824  -7.033021541950112e-05,
2825  -0.00843962585034014,
2826  -0.08376328656462587,
2827  -0.16991780045351476,
2828  -0.08376328656462587,
2829  -0.00843962585034014,
2830  -7.033021541950112e-05,
2831  0.0},
2832  {0.0,
2833  -0.00014266817838246406,
2834  -0.017120181405895694,
2835  -0.16991780045351473,
2836  -0.3446863189720333,
2837  -0.16991780045351473,
2838  -0.017120181405895694,
2839  -0.00014266817838246406,
2840  0.0},
2841  {0.0,
2842  -7.033021541950112e-05,
2843  -0.00843962585034014,
2844  -0.08376328656462587,
2845  -0.16991780045351476,
2846  -0.08376328656462587,
2847  -0.00843962585034014,
2848  -7.033021541950112e-05,
2849  0.0},
2850  {0.0,
2851  -7.086167800453517e-06,
2852  -0.0008503401360544218,
2853  -0.00843962585034014,
2854  -0.017120181405895694,
2855  -0.00843962585034014,
2856  -0.0008503401360544218,
2857  -7.086167800453517e-06,
2858  0.0},
2859  {0.0,
2860  -5.9051398337112626e-08,
2861  -7.086167800453517e-06,
2862  -7.033021541950112e-05,
2863  -0.00014266817838246408,
2864  -7.033021541950112e-05,
2865  -7.086167800453517e-06,
2866  -5.9051398337112626e-08,
2867  0.0},
2868  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
2869  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2870  {0.0,
2871  1.0498026371042242e-07,
2872  1.2597631645250694e-05,
2873  0.0001250314940791132,
2874  0.00025363231712438063,
2875  0.0001250314940791132,
2876  1.2597631645250694e-05,
2877  1.0498026371042242e-07,
2878  0.0},
2879  {0.0,
2880  1.259763164525069e-05,
2881  0.001511715797430084,
2882  0.015003779289493578,
2883  0.0304358780549257,
2884  0.015003779289493578,
2885  0.001511715797430084,
2886  1.259763164525069e-05,
2887  0.0},
2888  {0.0,
2889  0.0001250314940791132,
2890  0.015003779289493578,
2891  0.14891250944822357,
2892  0.3020760896951375,
2893  0.14891250944822357,
2894  0.015003779289493578,
2895  0.0001250314940791132,
2896  0.0},
2897  {0.0,
2898  0.0002536323171243806,
2899  0.0304358780549257,
2900  0.3020760896951375,
2901  0.612775678172503,
2902  0.3020760896951375,
2903  0.0304358780549257,
2904  0.0002536323171243806,
2905  0.0},
2906  {0.0,
2907  0.0001250314940791132,
2908  0.015003779289493578,
2909  0.14891250944822357,
2910  0.3020760896951375,
2911  0.14891250944822357,
2912  0.015003779289493578,
2913  0.0001250314940791132,
2914  0.0},
2915  {0.0,
2916  1.259763164525069e-05,
2917  0.001511715797430084,
2918  0.015003779289493578,
2919  0.0304358780549257,
2920  0.015003779289493578,
2921  0.001511715797430084,
2922  1.259763164525069e-05,
2923  0.0},
2924  {0.0,
2925  1.0498026371042242e-07,
2926  1.2597631645250694e-05,
2927  0.0001250314940791132,
2928  0.00025363231712438063,
2929  0.0001250314940791132,
2930  1.2597631645250694e-05,
2931  1.0498026371042242e-07,
2932  0.0},
2933  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
2934  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
2935  {0.0,
2936  -5.9051398337112626e-08,
2937  -7.086167800453517e-06,
2938  -7.033021541950112e-05,
2939  -0.00014266817838246408,
2940  -7.033021541950112e-05,
2941  -7.086167800453517e-06,
2942  -5.9051398337112626e-08,
2943  0.0},
2944  {0.0,
2945  -7.086167800453517e-06,
2946  -0.0008503401360544218,
2947  -0.00843962585034014,
2948  -0.017120181405895694,
2949  -0.00843962585034014,
2950  -0.0008503401360544218,
2951  -7.086167800453517e-06,
2952  0.0},
2953  {0.0,
2954  -7.033021541950112e-05,
2955  -0.00843962585034014,
2956  -0.08376328656462587,
2957  -0.16991780045351476,
2958  -0.08376328656462587,
2959  -0.00843962585034014,
2960  -7.033021541950112e-05,
2961  0.0},
2962  {0.0,
2963  -0.00014266817838246406,
2964  -0.017120181405895694,
2965  -0.16991780045351473,
2966  -0.3446863189720333,
2967  -0.16991780045351473,
2968  -0.017120181405895694,
2969  -0.00014266817838246406,
2970  0.0},
2971  {0.0,
2972  -7.033021541950112e-05,
2973  -0.00843962585034014,
2974  -0.08376328656462587,
2975  -0.16991780045351476,
2976  -0.08376328656462587,
2977  -0.00843962585034014,
2978  -7.033021541950112e-05,
2979  0.0},
2980  {0.0,
2981  -7.086167800453517e-06,
2982  -0.0008503401360544218,
2983  -0.00843962585034014,
2984  -0.017120181405895694,
2985  -0.00843962585034014,
2986  -0.0008503401360544218,
2987  -7.086167800453517e-06,
2988  0.0},
2989  {0.0,
2990  -5.9051398337112626e-08,
2991  -7.086167800453517e-06,
2992  -7.033021541950112e-05,
2993  -0.00014266817838246408,
2994  -7.033021541950112e-05,
2995  -7.086167800453517e-06,
2996  -5.9051398337112626e-08,
2997  0.0},
2998  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
2999  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3000  {0.0,
3001  0.0,
3002  8.933550615572699e-23,
3003  6.617444900424221e-22,
3004  -1.4097225589434974e-21,
3005  6.617444900424221e-22,
3006  8.933550615572699e-23,
3007  0.0,
3008  0.0},
3009  {0.0,
3010  -1.3896634290890865e-22,
3011  -2.948733447629033e-20,
3012  3.095640724418451e-19,
3013  -1.2023219095968517e-18,
3014  3.095640724418451e-19,
3015  -2.948733447629033e-20,
3016  -1.3896634290890865e-22,
3017  0.0},
3018  {0.0,
3019  -1.2176098616780567e-21,
3020  -3.814295240604521e-20,
3021  -2.910855193019005e-18,
3022  -5.991280757250155e-18,
3023  -2.910855193019005e-18,
3024  -3.814295240604521e-20,
3025  -1.2176098616780567e-21,
3026  0.0},
3027  {0.0,
3028  -4.5135110123955955e-21,
3029  -1.0007148332605274e-18,
3030  -7.31572264940066e-18,
3031  2.069796097619653e-17,
3032  -7.31572264940066e-18,
3033  -1.0007148332605274e-18,
3034  -4.5135110123955955e-21,
3035  0.0},
3036  {0.0,
3037  -1.2176098616780567e-21,
3038  -3.814295240604521e-20,
3039  -2.910855193019005e-18,
3040  -5.991280757250155e-18,
3041  -2.910855193019005e-18,
3042  -3.814295240604521e-20,
3043  -1.2176098616780567e-21,
3044  0.0},
3045  {0.0,
3046  -1.3896634290890865e-22,
3047  -2.948733447629033e-20,
3048  3.095640724418451e-19,
3049  -1.2023219095968517e-18,
3050  3.095640724418451e-19,
3051  -2.948733447629033e-20,
3052  -1.3896634290890865e-22,
3053  0.0},
3054  {0.0,
3055  0.0,
3056  8.933550615572699e-23,
3057  6.617444900424221e-22,
3058  -1.4097225589434974e-21,
3059  6.617444900424221e-22,
3060  8.933550615572699e-23,
3061  0.0,
3062  0.0},
3063  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3064  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3065  {0.0,
3066  6.561266481901402e-09,
3067  7.873519778281682e-07,
3068  7.814468379944569e-06,
3069  1.5852019820273786e-05,
3070  7.814468379944569e-06,
3071  7.873519778281682e-07,
3072  6.561266481901402e-09,
3073  0.0},
3074  {0.0,
3075  7.873519778281682e-07,
3076  9.44822373393802e-05,
3077  0.0009377362055933485,
3078  0.0019022423784328549,
3079  0.0009377362055933485,
3080  9.44822373393802e-05,
3081  7.873519778281682e-07,
3082  0.0},
3083  {0.0,
3084  7.81446837994457e-06,
3085  0.0009377362055933487,
3086  0.00930703184051398,
3087  0.018879755605946087,
3088  0.00930703184051398,
3089  0.0009377362055933487,
3090  7.81446837994457e-06,
3091  0.0},
3092  {0.0,
3093  1.585201982027379e-05,
3094  0.0019022423784328547,
3095  0.018879755605946083,
3096  0.03829847988578147,
3097  0.018879755605946083,
3098  0.0019022423784328547,
3099  1.585201982027379e-05,
3100  0.0},
3101  {0.0,
3102  7.81446837994457e-06,
3103  0.0009377362055933487,
3104  0.00930703184051398,
3105  0.018879755605946087,
3106  0.00930703184051398,
3107  0.0009377362055933487,
3108  7.81446837994457e-06,
3109  0.0},
3110  {0.0,
3111  7.873519778281682e-07,
3112  9.44822373393802e-05,
3113  0.0009377362055933485,
3114  0.0019022423784328549,
3115  0.0009377362055933485,
3116  9.44822373393802e-05,
3117  7.873519778281682e-07,
3118  0.0},
3119  {0.0,
3120  6.561266481901402e-09,
3121  7.873519778281682e-07,
3122  7.814468379944569e-06,
3123  1.5852019820273786e-05,
3124  7.814468379944569e-06,
3125  7.873519778281682e-07,
3126  6.561266481901402e-09,
3127  0.0},
3128  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3129  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3130  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3131  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3132  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3133  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3134  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3135  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3136  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3137  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}}};
3138 
3139  double U110[9][9][9] =
3140 {{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3141  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3142  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3143  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3144  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3145  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3146  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3147  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3148  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3149  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3150  {0.0,
3151  1.3778659611992942e-08,
3152  1.6534391534391533e-06,
3153  1.64103835978836e-05,
3154  3.328924162257494e-05,
3155  1.64103835978836e-05,
3156  1.6534391534391533e-06,
3157  1.3778659611992942e-08,
3158  0.0},
3159  {0.0,
3160  3.3068783068783074e-07,
3161  3.9682539682539696e-05,
3162  0.0003938492063492064,
3163  0.0007989417989417992,
3164  0.0003938492063492064,
3165  3.9682539682539696e-05,
3166  3.3068783068783074e-07,
3167  0.0},
3168  {0.0,
3169  2.0667989417989414e-07,
3170  2.4801587301587298e-05,
3171  0.0002461557539682539,
3172  0.0004993386243386242,
3173  0.0002461557539682539,
3174  2.4801587301587298e-05,
3175  2.0667989417989414e-07,
3176  0.0},
3177  {0.0,
3178  -1.1022927689594355e-06,
3179  -0.0001322751322751323,
3180  -0.0013128306878306879,
3181  -0.0026631393298059956,
3182  -0.0013128306878306879,
3183  -0.0001322751322751323,
3184  -1.1022927689594355e-06,
3185  0.0},
3186  {0.0,
3187  2.0667989417989414e-07,
3188  2.4801587301587298e-05,
3189  0.0002461557539682539,
3190  0.0004993386243386242,
3191  0.0002461557539682539,
3192  2.4801587301587298e-05,
3193  2.0667989417989414e-07,
3194  0.0},
3195  {0.0,
3196  3.3068783068783074e-07,
3197  3.9682539682539696e-05,
3198  0.0003938492063492064,
3199  0.0007989417989417992,
3200  0.0003938492063492064,
3201  3.9682539682539696e-05,
3202  3.3068783068783074e-07,
3203  0.0},
3204  {0.0,
3205  1.3778659611992942e-08,
3206  1.6534391534391533e-06,
3207  1.64103835978836e-05,
3208  3.328924162257494e-05,
3209  1.64103835978836e-05,
3210  1.6534391534391533e-06,
3211  1.3778659611992942e-08,
3212  0.0},
3213  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3214  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3215  {0.0,
3216  3.306878306878308e-07,
3217  3.9682539682539696e-05,
3218  0.00039384920634920654,
3219  0.000798941798941799,
3220  0.00039384920634920654,
3221  3.9682539682539696e-05,
3222  3.306878306878308e-07,
3223  0.0},
3224  {0.0,
3225  7.936507936507936e-06,
3226  0.0009523809523809524,
3227  0.009452380952380952,
3228  0.019174603174603164,
3229  0.009452380952380952,
3230  0.0009523809523809524,
3231  7.936507936507936e-06,
3232  0.0},
3233  {0.0,
3234  4.960317460317461e-06,
3235  0.0005952380952380953,
3236  0.005907738095238096,
3237  0.011984126984126979,
3238  0.005907738095238096,
3239  0.0005952380952380953,
3240  4.960317460317461e-06,
3241  0.0},
3242  {0.0,
3243  -2.6455026455026446e-05,
3244  -0.003174603174603174,
3245  -0.03150793650793652,
3246  -0.06391534391534391,
3247  -0.03150793650793652,
3248  -0.003174603174603174,
3249  -2.6455026455026446e-05,
3250  0.0},
3251  {0.0,
3252  4.960317460317461e-06,
3253  0.0005952380952380953,
3254  0.005907738095238096,
3255  0.011984126984126979,
3256  0.005907738095238096,
3257  0.0005952380952380953,
3258  4.960317460317461e-06,
3259  0.0},
3260  {0.0,
3261  7.936507936507936e-06,
3262  0.0009523809523809524,
3263  0.009452380952380952,
3264  0.019174603174603164,
3265  0.009452380952380952,
3266  0.0009523809523809524,
3267  7.936507936507936e-06,
3268  0.0},
3269  {0.0,
3270  3.306878306878308e-07,
3271  3.9682539682539696e-05,
3272  0.00039384920634920654,
3273  0.000798941798941799,
3274  0.00039384920634920654,
3275  3.9682539682539696e-05,
3276  3.306878306878308e-07,
3277  0.0},
3278  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3279  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3280  {0.0,
3281  2.0667989417989414e-07,
3282  2.48015873015873e-05,
3283  0.00024615575396825384,
3284  0.0004993386243386241,
3285  0.00024615575396825384,
3286  2.48015873015873e-05,
3287  2.0667989417989414e-07,
3288  0.0},
3289  {0.0,
3290  4.960317460317462e-06,
3291  0.0005952380952380951,
3292  0.005907738095238095,
3293  0.011984126984126982,
3294  0.005907738095238095,
3295  0.0005952380952380951,
3296  4.960317460317462e-06,
3297  0.0},
3298  {0.0,
3299  3.100198412698412e-06,
3300  0.0003720238095238095,
3301  0.003692336309523805,
3302  0.007490079365079362,
3303  0.003692336309523805,
3304  0.0003720238095238095,
3305  3.100198412698412e-06,
3306  0.0},
3307  {0.0,
3308  -1.6534391534391525e-05,
3309  -0.001984126984126985,
3310  -0.019692460317460303,
3311  -0.039947089947089946,
3312  -0.019692460317460303,
3313  -0.001984126984126985,
3314  -1.6534391534391525e-05,
3315  0.0},
3316  {0.0,
3317  3.100198412698412e-06,
3318  0.0003720238095238095,
3319  0.003692336309523805,
3320  0.007490079365079362,
3321  0.003692336309523805,
3322  0.0003720238095238095,
3323  3.100198412698412e-06,
3324  0.0},
3325  {0.0,
3326  4.960317460317462e-06,
3327  0.0005952380952380951,
3328  0.005907738095238095,
3329  0.011984126984126982,
3330  0.005907738095238095,
3331  0.0005952380952380951,
3332  4.960317460317462e-06,
3333  0.0},
3334  {0.0,
3335  2.0667989417989414e-07,
3336  2.48015873015873e-05,
3337  0.00024615575396825384,
3338  0.0004993386243386241,
3339  0.00024615575396825384,
3340  2.48015873015873e-05,
3341  2.0667989417989414e-07,
3342  0.0},
3343  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3344  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3345  {0.0,
3346  -1.1022927689594353e-06,
3347  -0.00013227513227513228,
3348  -0.001312830687830688,
3349  -0.002663139329805994,
3350  -0.001312830687830688,
3351  -0.00013227513227513228,
3352  -1.1022927689594353e-06,
3353  0.0},
3354  {0.0,
3355  -2.6455026455026446e-05,
3356  -0.003174603174603174,
3357  -0.031507936507936485,
3358  -0.06391534391534393,
3359  -0.031507936507936485,
3360  -0.003174603174603174,
3361  -2.6455026455026446e-05,
3362  0.0},
3363  {0.0,
3364  -1.653439153439153e-05,
3365  -0.001984126984126985,
3366  -0.019692460317460303,
3367  -0.03994708994708995,
3368  -0.019692460317460303,
3369  -0.001984126984126985,
3370  -1.653439153439153e-05,
3371  0.0},
3372  {0.0,
3373  8.818342151675487e-05,
3374  0.010582010582010571,
3375  0.10502645502645507,
3376  0.2130511463844796,
3377  0.10502645502645507,
3378  0.010582010582010571,
3379  8.818342151675487e-05,
3380  0.0},
3381  {0.0,
3382  -1.653439153439153e-05,
3383  -0.001984126984126985,
3384  -0.019692460317460303,
3385  -0.03994708994708995,
3386  -0.019692460317460303,
3387  -0.001984126984126985,
3388  -1.653439153439153e-05,
3389  0.0},
3390  {0.0,
3391  -2.6455026455026446e-05,
3392  -0.003174603174603174,
3393  -0.031507936507936485,
3394  -0.06391534391534393,
3395  -0.031507936507936485,
3396  -0.003174603174603174,
3397  -2.6455026455026446e-05,
3398  0.0},
3399  {0.0,
3400  -1.1022927689594353e-06,
3401  -0.00013227513227513228,
3402  -0.001312830687830688,
3403  -0.002663139329805994,
3404  -0.001312830687830688,
3405  -0.00013227513227513228,
3406  -1.1022927689594353e-06,
3407  0.0},
3408  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3409  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3410  {0.0,
3411  2.0667989417989414e-07,
3412  2.48015873015873e-05,
3413  0.00024615575396825384,
3414  0.0004993386243386241,
3415  0.00024615575396825384,
3416  2.48015873015873e-05,
3417  2.0667989417989414e-07,
3418  0.0},
3419  {0.0,
3420  4.960317460317462e-06,
3421  0.0005952380952380951,
3422  0.005907738095238095,
3423  0.011984126984126982,
3424  0.005907738095238095,
3425  0.0005952380952380951,
3426  4.960317460317462e-06,
3427  0.0},
3428  {0.0,
3429  3.100198412698412e-06,
3430  0.0003720238095238095,
3431  0.003692336309523805,
3432  0.007490079365079362,
3433  0.003692336309523805,
3434  0.0003720238095238095,
3435  3.100198412698412e-06,
3436  0.0},
3437  {0.0,
3438  -1.6534391534391525e-05,
3439  -0.001984126984126985,
3440  -0.019692460317460303,
3441  -0.039947089947089946,
3442  -0.019692460317460303,
3443  -0.001984126984126985,
3444  -1.6534391534391525e-05,
3445  0.0},
3446  {0.0,
3447  3.100198412698412e-06,
3448  0.0003720238095238095,
3449  0.003692336309523805,
3450  0.007490079365079362,
3451  0.003692336309523805,
3452  0.0003720238095238095,
3453  3.100198412698412e-06,
3454  0.0},
3455  {0.0,
3456  4.960317460317462e-06,
3457  0.0005952380952380951,
3458  0.005907738095238095,
3459  0.011984126984126982,
3460  0.005907738095238095,
3461  0.0005952380952380951,
3462  4.960317460317462e-06,
3463  0.0},
3464  {0.0,
3465  2.0667989417989414e-07,
3466  2.48015873015873e-05,
3467  0.00024615575396825384,
3468  0.0004993386243386241,
3469  0.00024615575396825384,
3470  2.48015873015873e-05,
3471  2.0667989417989414e-07,
3472  0.0},
3473  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3474  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3475  {0.0,
3476  3.306878306878308e-07,
3477  3.9682539682539696e-05,
3478  0.00039384920634920654,
3479  0.000798941798941799,
3480  0.00039384920634920654,
3481  3.9682539682539696e-05,
3482  3.306878306878308e-07,
3483  0.0},
3484  {0.0,
3485  7.936507936507936e-06,
3486  0.0009523809523809524,
3487  0.009452380952380952,
3488  0.019174603174603164,
3489  0.009452380952380952,
3490  0.0009523809523809524,
3491  7.936507936507936e-06,
3492  0.0},
3493  {0.0,
3494  4.960317460317461e-06,
3495  0.0005952380952380953,
3496  0.005907738095238096,
3497  0.011984126984126979,
3498  0.005907738095238096,
3499  0.0005952380952380953,
3500  4.960317460317461e-06,
3501  0.0},
3502  {0.0,
3503  -2.6455026455026446e-05,
3504  -0.003174603174603174,
3505  -0.03150793650793652,
3506  -0.06391534391534391,
3507  -0.03150793650793652,
3508  -0.003174603174603174,
3509  -2.6455026455026446e-05,
3510  0.0},
3511  {0.0,
3512  4.960317460317461e-06,
3513  0.0005952380952380953,
3514  0.005907738095238096,
3515  0.011984126984126979,
3516  0.005907738095238096,
3517  0.0005952380952380953,
3518  4.960317460317461e-06,
3519  0.0},
3520  {0.0,
3521  7.936507936507936e-06,
3522  0.0009523809523809524,
3523  0.009452380952380952,
3524  0.019174603174603164,
3525  0.009452380952380952,
3526  0.0009523809523809524,
3527  7.936507936507936e-06,
3528  0.0},
3529  {0.0,
3530  3.306878306878308e-07,
3531  3.9682539682539696e-05,
3532  0.00039384920634920654,
3533  0.000798941798941799,
3534  0.00039384920634920654,
3535  3.9682539682539696e-05,
3536  3.306878306878308e-07,
3537  0.0},
3538  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3539  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3540  {0.0,
3541  1.3778659611992942e-08,
3542  1.6534391534391533e-06,
3543  1.64103835978836e-05,
3544  3.328924162257494e-05,
3545  1.64103835978836e-05,
3546  1.6534391534391533e-06,
3547  1.3778659611992942e-08,
3548  0.0},
3549  {0.0,
3550  3.3068783068783074e-07,
3551  3.9682539682539696e-05,
3552  0.0003938492063492064,
3553  0.0007989417989417992,
3554  0.0003938492063492064,
3555  3.9682539682539696e-05,
3556  3.3068783068783074e-07,
3557  0.0},
3558  {0.0,
3559  2.0667989417989414e-07,
3560  2.4801587301587298e-05,
3561  0.0002461557539682539,
3562  0.0004993386243386242,
3563  0.0002461557539682539,
3564  2.4801587301587298e-05,
3565  2.0667989417989414e-07,
3566  0.0},
3567  {0.0,
3568  -1.1022927689594355e-06,
3569  -0.0001322751322751323,
3570  -0.0013128306878306879,
3571  -0.0026631393298059956,
3572  -0.0013128306878306879,
3573  -0.0001322751322751323,
3574  -1.1022927689594355e-06,
3575  0.0},
3576  {0.0,
3577  2.0667989417989414e-07,
3578  2.4801587301587298e-05,
3579  0.0002461557539682539,
3580  0.0004993386243386242,
3581  0.0002461557539682539,
3582  2.4801587301587298e-05,
3583  2.0667989417989414e-07,
3584  0.0},
3585  {0.0,
3586  3.3068783068783074e-07,
3587  3.9682539682539696e-05,
3588  0.0003938492063492064,
3589  0.0007989417989417992,
3590  0.0003938492063492064,
3591  3.9682539682539696e-05,
3592  3.3068783068783074e-07,
3593  0.0},
3594  {0.0,
3595  1.3778659611992942e-08,
3596  1.6534391534391533e-06,
3597  1.64103835978836e-05,
3598  3.328924162257494e-05,
3599  1.64103835978836e-05,
3600  1.6534391534391533e-06,
3601  1.3778659611992942e-08,
3602  0.0},
3603  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3604  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3605  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3606  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3607  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3608  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3609  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3610  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3611  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3612  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}}};
3613 
3614 double U100[9][9][9] =
3615 {{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3616  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3617  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3618  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3619  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3620  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3621  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3622  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3623  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3624  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3625  {0.0,
3626  -3.280633240950702e-10,
3627  -3.936759889140842e-08,
3628  -3.907234189972286e-07,
3629  -7.926009910136894e-07,
3630  -3.907234189972286e-07,
3631  -3.936759889140842e-08,
3632  -3.280633240950702e-10,
3633  0.0},
3634  {0.0,
3635  -3.936759889140842e-08,
3636  -4.724111866969011e-06,
3637  -4.688681027966742e-05,
3638  -9.511211892164272e-05,
3639  -4.688681027966742e-05,
3640  -4.724111866969011e-06,
3641  -3.936759889140842e-08,
3642  0.0},
3643  {0.0,
3644  -3.9072341899722857e-07,
3645  -4.688681027966742e-05,
3646  -0.0004653515920256992,
3647  -0.000943987780297304,
3648  -0.0004653515920256992,
3649  -4.688681027966742e-05,
3650  -3.9072341899722857e-07,
3651  0.0},
3652  {0.0,
3653  -7.926009910136895e-07,
3654  -9.511211892164274e-05,
3655  -0.000943987780297304,
3656  -0.0019149239942890736,
3657  -0.000943987780297304,
3658  -9.511211892164274e-05,
3659  -7.926009910136895e-07,
3660  0.0},
3661  {0.0,
3662  -3.9072341899722857e-07,
3663  -4.688681027966742e-05,
3664  -0.0004653515920256992,
3665  -0.000943987780297304,
3666  -0.0004653515920256992,
3667  -4.688681027966742e-05,
3668  -3.9072341899722857e-07,
3669  0.0},
3670  {0.0,
3671  -3.936759889140842e-08,
3672  -4.724111866969011e-06,
3673  -4.688681027966742e-05,
3674  -9.511211892164272e-05,
3675  -4.688681027966742e-05,
3676  -4.724111866969011e-06,
3677  -3.936759889140842e-08,
3678  0.0},
3679  {0.0,
3680  -3.280633240950702e-10,
3681  -3.936759889140842e-08,
3682  -3.907234189972286e-07,
3683  -7.926009910136894e-07,
3684  -3.907234189972286e-07,
3685  -3.936759889140842e-08,
3686  -3.280633240950702e-10,
3687  0.0},
3688  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3689  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3690  {0.0,
3691  -7.873519778281682e-09,
3692  -9.44822373393802e-07,
3693  -9.377362055933484e-06,
3694  -1.9022423784328537e-05,
3695  -9.377362055933484e-06,
3696  -9.44822373393802e-07,
3697  -7.873519778281682e-09,
3698  0.0},
3699  {0.0,
3700  -9.448223733938021e-07,
3701  -0.00011337868480725627,
3702  -0.0011252834467120182,
3703  -0.002282690854119424,
3704  -0.0011252834467120182,
3705  -0.00011337868480725627,
3706  -9.448223733938021e-07,
3707  0.0},
3708  {0.0,
3709  -9.377362055933486e-06,
3710  -0.0011252834467120182,
3711  -0.011168438208616769,
3712  -0.022655706727135298,
3713  -0.011168438208616769,
3714  -0.0011252834467120182,
3715  -9.377362055933486e-06,
3716  0.0},
3717  {0.0,
3718  -1.9022423784328537e-05,
3719  -0.002282690854119425,
3720  -0.0226557067271353,
3721  -0.045958175862937774,
3722  -0.0226557067271353,
3723  -0.002282690854119425,
3724  -1.9022423784328537e-05,
3725  0.0},
3726  {0.0,
3727  -9.377362055933486e-06,
3728  -0.0011252834467120182,
3729  -0.011168438208616769,
3730  -0.022655706727135298,
3731  -0.011168438208616769,
3732  -0.0011252834467120182,
3733  -9.377362055933486e-06,
3734  0.0},
3735  {0.0,
3736  -9.448223733938021e-07,
3737  -0.00011337868480725627,
3738  -0.0011252834467120182,
3739  -0.002282690854119424,
3740  -0.0011252834467120182,
3741  -0.00011337868480725627,
3742  -9.448223733938021e-07,
3743  0.0},
3744  {0.0,
3745  -7.873519778281682e-09,
3746  -9.44822373393802e-07,
3747  -9.377362055933484e-06,
3748  -1.9022423784328537e-05,
3749  -9.377362055933484e-06,
3750  -9.44822373393802e-07,
3751  -7.873519778281682e-09,
3752  0.0},
3753  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3754  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3755  {0.0,
3756  -4.92094986142605e-09,
3757  -5.905139833711262e-07,
3758  -5.860851284958427e-06,
3759  -1.1889014865205345e-05,
3760  -5.860851284958427e-06,
3761  -5.905139833711262e-07,
3762  -4.92094986142605e-09,
3763  0.0},
3764  {0.0,
3765  -5.905139833711262e-07,
3766  -7.086167800453512e-05,
3767  -0.0007033021541950114,
3768  -0.0014266817838246414,
3769  -0.0007033021541950114,
3770  -7.086167800453512e-05,
3771  -5.905139833711262e-07,
3772  0.0},
3773  {0.0,
3774  -5.860851284958427e-06,
3775  -0.0007033021541950114,
3776  -0.006980273880385488,
3777  -0.014159816704459559,
3778  -0.006980273880385488,
3779  -0.0007033021541950114,
3780  -5.860851284958427e-06,
3781  0.0},
3782  {0.0,
3783  -1.1889014865205344e-05,
3784  -0.0014266817838246414,
3785  -0.014159816704459562,
3786  -0.028723859914336125,
3787  -0.014159816704459562,
3788  -0.0014266817838246414,
3789  -1.1889014865205344e-05,
3790  0.0},
3791  {0.0,
3792  -5.860851284958427e-06,
3793  -0.0007033021541950114,
3794  -0.006980273880385488,
3795  -0.014159816704459559,
3796  -0.006980273880385488,
3797  -0.0007033021541950114,
3798  -5.860851284958427e-06,
3799  0.0},
3800  {0.0,
3801  -5.905139833711262e-07,
3802  -7.086167800453512e-05,
3803  -0.0007033021541950114,
3804  -0.0014266817838246414,
3805  -0.0007033021541950114,
3806  -7.086167800453512e-05,
3807  -5.905139833711262e-07,
3808  0.0},
3809  {0.0,
3810  -4.92094986142605e-09,
3811  -5.905139833711262e-07,
3812  -5.860851284958427e-06,
3813  -1.1889014865205345e-05,
3814  -5.860851284958427e-06,
3815  -5.905139833711262e-07,
3816  -4.92094986142605e-09,
3817  0.0},
3818  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3819  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3820  {0.0,
3821  2.6245065927605612e-08,
3822  3.1494079113126745e-06,
3823  3.12578735197783e-05,
3824  6.340807928109514e-05,
3825  3.12578735197783e-05,
3826  3.1494079113126745e-06,
3827  2.6245065927605612e-08,
3828  0.0},
3829  {0.0,
3830  3.1494079113126745e-06,
3831  0.0003779289493575207,
3832  0.0037509448223733938,
3833  0.007608969513731416,
3834  0.0037509448223733938,
3835  0.0003779289493575207,
3836  3.1494079113126745e-06,
3837  0.0},
3838  {0.0,
3839  3.1257873519778296e-05,
3840  0.0037509448223733938,
3841  0.03722812736205595,
3842  0.07551902242378432,
3843  0.03722812736205595,
3844  0.0037509448223733938,
3845  3.1257873519778296e-05,
3846  0.0},
3847  {0.0,
3848  6.340807928109514e-05,
3849  0.007608969513731417,
3850  0.07551902242378433,
3851  0.15319391954312578,
3852  0.07551902242378433,
3853  0.007608969513731417,
3854  6.340807928109514e-05,
3855  0.0},
3856  {0.0,
3857  3.1257873519778296e-05,
3858  0.0037509448223733938,
3859  0.03722812736205595,
3860  0.07551902242378432,
3861  0.03722812736205595,
3862  0.0037509448223733938,
3863  3.1257873519778296e-05,
3864  0.0},
3865  {0.0,
3866  3.1494079113126745e-06,
3867  0.0003779289493575207,
3868  0.0037509448223733938,
3869  0.007608969513731416,
3870  0.0037509448223733938,
3871  0.0003779289493575207,
3872  3.1494079113126745e-06,
3873  0.0},
3874  {0.0,
3875  2.6245065927605612e-08,
3876  3.1494079113126745e-06,
3877  3.12578735197783e-05,
3878  6.340807928109514e-05,
3879  3.12578735197783e-05,
3880  3.1494079113126745e-06,
3881  2.6245065927605612e-08,
3882  0.0},
3883  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3884  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3885  {0.0,
3886  -4.92094986142605e-09,
3887  -5.905139833711262e-07,
3888  -5.860851284958427e-06,
3889  -1.1889014865205345e-05,
3890  -5.860851284958427e-06,
3891  -5.905139833711262e-07,
3892  -4.92094986142605e-09,
3893  0.0},
3894  {0.0,
3895  -5.905139833711262e-07,
3896  -7.086167800453512e-05,
3897  -0.0007033021541950114,
3898  -0.0014266817838246414,
3899  -0.0007033021541950114,
3900  -7.086167800453512e-05,
3901  -5.905139833711262e-07,
3902  0.0},
3903  {0.0,
3904  -5.860851284958427e-06,
3905  -0.0007033021541950114,
3906  -0.006980273880385488,
3907  -0.014159816704459559,
3908  -0.006980273880385488,
3909  -0.0007033021541950114,
3910  -5.860851284958427e-06,
3911  0.0},
3912  {0.0,
3913  -1.1889014865205344e-05,
3914  -0.0014266817838246414,
3915  -0.014159816704459562,
3916  -0.028723859914336125,
3917  -0.014159816704459562,
3918  -0.0014266817838246414,
3919  -1.1889014865205344e-05,
3920  0.0},
3921  {0.0,
3922  -5.860851284958427e-06,
3923  -0.0007033021541950114,
3924  -0.006980273880385488,
3925  -0.014159816704459559,
3926  -0.006980273880385488,
3927  -0.0007033021541950114,
3928  -5.860851284958427e-06,
3929  0.0},
3930  {0.0,
3931  -5.905139833711262e-07,
3932  -7.086167800453512e-05,
3933  -0.0007033021541950114,
3934  -0.0014266817838246414,
3935  -0.0007033021541950114,
3936  -7.086167800453512e-05,
3937  -5.905139833711262e-07,
3938  0.0},
3939  {0.0,
3940  -4.92094986142605e-09,
3941  -5.905139833711262e-07,
3942  -5.860851284958427e-06,
3943  -1.1889014865205345e-05,
3944  -5.860851284958427e-06,
3945  -5.905139833711262e-07,
3946  -4.92094986142605e-09,
3947  0.0},
3948  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
3949  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
3950  {0.0,
3951  -7.873519778281682e-09,
3952  -9.44822373393802e-07,
3953  -9.377362055933484e-06,
3954  -1.9022423784328537e-05,
3955  -9.377362055933484e-06,
3956  -9.44822373393802e-07,
3957  -7.873519778281682e-09,
3958  0.0},
3959  {0.0,
3960  -9.448223733938021e-07,
3961  -0.00011337868480725627,
3962  -0.0011252834467120182,
3963  -0.002282690854119424,
3964  -0.0011252834467120182,
3965  -0.00011337868480725627,
3966  -9.448223733938021e-07,
3967  0.0},
3968  {0.0,
3969  -9.377362055933486e-06,
3970  -0.0011252834467120182,
3971  -0.011168438208616769,
3972  -0.022655706727135298,
3973  -0.011168438208616769,
3974  -0.0011252834467120182,
3975  -9.377362055933486e-06,
3976  0.0},
3977  {0.0,
3978  -1.9022423784328537e-05,
3979  -0.002282690854119425,
3980  -0.0226557067271353,
3981  -0.045958175862937774,
3982  -0.0226557067271353,
3983  -0.002282690854119425,
3984  -1.9022423784328537e-05,
3985  0.0},
3986  {0.0,
3987  -9.377362055933486e-06,
3988  -0.0011252834467120182,
3989  -0.011168438208616769,
3990  -0.022655706727135298,
3991  -0.011168438208616769,
3992  -0.0011252834467120182,
3993  -9.377362055933486e-06,
3994  0.0},
3995  {0.0,
3996  -9.448223733938021e-07,
3997  -0.00011337868480725627,
3998  -0.0011252834467120182,
3999  -0.002282690854119424,
4000  -0.0011252834467120182,
4001  -0.00011337868480725627,
4002  -9.448223733938021e-07,
4003  0.0},
4004  {0.0,
4005  -7.873519778281682e-09,
4006  -9.44822373393802e-07,
4007  -9.377362055933484e-06,
4008  -1.9022423784328537e-05,
4009  -9.377362055933484e-06,
4010  -9.44822373393802e-07,
4011  -7.873519778281682e-09,
4012  0.0},
4013  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
4014  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
4015  {0.0,
4016  -3.280633240950702e-10,
4017  -3.936759889140842e-08,
4018  -3.907234189972286e-07,
4019  -7.926009910136894e-07,
4020  -3.907234189972286e-07,
4021  -3.936759889140842e-08,
4022  -3.280633240950702e-10,
4023  0.0},
4024  {0.0,
4025  -3.936759889140842e-08,
4026  -4.724111866969011e-06,
4027  -4.688681027966742e-05,
4028  -9.511211892164272e-05,
4029  -4.688681027966742e-05,
4030  -4.724111866969011e-06,
4031  -3.936759889140842e-08,
4032  0.0},
4033  {0.0,
4034  -3.9072341899722857e-07,
4035  -4.688681027966742e-05,
4036  -0.0004653515920256992,
4037  -0.000943987780297304,
4038  -0.0004653515920256992,
4039  -4.688681027966742e-05,
4040  -3.9072341899722857e-07,
4041  0.0},
4042  {0.0,
4043  -7.926009910136895e-07,
4044  -9.511211892164274e-05,
4045  -0.000943987780297304,
4046  -0.0019149239942890736,
4047  -0.000943987780297304,
4048  -9.511211892164274e-05,
4049  -7.926009910136895e-07,
4050  0.0},
4051  {0.0,
4052  -3.9072341899722857e-07,
4053  -4.688681027966742e-05,
4054  -0.0004653515920256992,
4055  -0.000943987780297304,
4056  -0.0004653515920256992,
4057  -4.688681027966742e-05,
4058  -3.9072341899722857e-07,
4059  0.0},
4060  {0.0,
4061  -3.936759889140842e-08,
4062  -4.724111866969011e-06,
4063  -4.688681027966742e-05,
4064  -9.511211892164272e-05,
4065  -4.688681027966742e-05,
4066  -4.724111866969011e-06,
4067  -3.936759889140842e-08,
4068  0.0},
4069  {0.0,
4070  -3.280633240950702e-10,
4071  -3.936759889140842e-08,
4072  -3.907234189972286e-07,
4073  -7.926009910136894e-07,
4074  -3.907234189972286e-07,
4075  -3.936759889140842e-08,
4076  -3.280633240950702e-10,
4077  0.0},
4078  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
4079  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
4080  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
4081  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
4082  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
4083  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
4084  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
4085  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
4086  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
4087  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}}};
4088 
4089 
4090 };
4091 
4096 } // namespace npl
4097 
4098 #endif //ACCESSORS_H
npl::LanczosInterp3DView::m_ras
bool m_ras
if true, then this assumes the inputs are RAS coordinates rather than indexes. Default is false ...
Definition: accessors.h:2014
INVALID_ARGUMENT
#define INVALID_ARGUMENT(EXP)
Definition: macros.h:18
npl::LinInterp3DView::m_ras
bool m_ras
if true, then this assumes the inputs are RAS coordinates rather than indexes. Default is false ...
Definition: accessors.h:1265
npl::LinInterpNDView::operator()
T operator()(std::initializer_list< double > index)
Gets value at array index and then casts to T.
Definition: accessors.h:947
npl
Definition: accessors.h:29
npl::NDView::set
void set(size_t len, const int64_t *index, T v)
Casts to the appropriate type then sets array at given index.
Definition: accessors.h:218
npl::BSplineView::reconstruct
ptr< MRImage > reconstruct(ptr< const MRImage > input)
Perform full-on reconstruction in the space of the input image.
Definition: accessors.h:2274
npl::Vector3DView::Vector3DView
Vector3DView(std::shared_ptr< NDArray > in)
Definition: accessors.h:700
npl::Pixel3DView::operator[]
T operator[](int64_t i)
Definition: accessors.h:544
npl::NNInterpNDView::operator()
T operator()(double x=0, double y=0, double z=0, double t=0, double u=0, double v=0, double w=0, double q=0)
Gets value at array index and then casts to T.
Definition: accessors.h:1299
npl::Slicer::setROI
void setROI(size_t len, const size_t *roisize, const int64_t *roistart=NULL)
Sets the region of interest, with lower bound of 0. During iteration or any motion the position will ...
npl::NDView::castset
void(* castset)(void *ptr, const T &v)
Function pointer to the correct function for casting to the underlying type. This should be set durin...
Definition: accessors.h:313
npl::LinInterp3DView::operator()
T operator()(double x=0, double y=0, double z=0, int64_t t=0)
Gets value at array index and then casts to T.
Definition: accessors.h:1153
npl::NDView::castgetStatic
static T castgetStatic(void *ptr)
This is a wrapper function that will be called to safely cast from the underlying type...
Definition: accessors.h:272
npl::Slicer::index
void index(size_t len, int64_t *index) const
Places the first len dimension in the given array. If the number of dimensions exceed the len then th...
npl::LanczosInterp3DView::operator[]
T operator[](const std::vector< int64_t > &i)
Definition: accessors.h:2025
npl::LinInterpNDView::operator[]
T operator[](const std::vector< int64_t > &index)
Gets value at array index and then casts to T.
Definition: accessors.h:1025
npl::NNInterpNDView::operator()
T operator()(std::initializer_list< float > index)
Gets value at array index and then casts to T.
Definition: accessors.h:1411
npl::NDIter::set
void set(T v)
Dereference operator.
Definition: iterators.h:1005
npl::NDConstIter::eof
bool eof() const
Are we one past the last element?
Definition: iterators.h:731
npl::Counter::advance
bool advance(bool rorder=false)
Advance through ND-counter. 0,0,0 to 0,0,1 and so on. If roder is true then this will start at 0...
Definition: basic_functions.h:581
npl::NDConstView::setArray
void setArray(ptr< const NDArray > in)
Definition: accessors.h:334
npl::Vector3DView::operator()
T operator()(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
Gets value at array index and then casts to T.
Definition: accessors.h:709
npl::LanczosInterp3DView::setRadius
void setRadius(size_t rad)
Definition: accessors.h:1877
npl::BSplineView::BSplineView
BSplineView(std::shared_ptr< MRImage > params, BoundaryConditionT bound=ZEROFLUX)
Definition: accessors.h:2052
npl::NNInterpNDView::NNInterpNDView
NNInterpNDView(std::shared_ptr< const NDArray > in, BoundaryConditionT bound=ZEROFLUX)
Definition: accessors.h:1278
npl::LinInterpNDView::operator()
T operator()(const std::vector< float > &index)
Gets value at array index and then casts to T.
Definition: accessors.h:918
npl::Counter
Very basic counter that iterates over an ND region.
Definition: basic_functions.h:532
npl::LinInterp3DView::LinInterp3DView
LinInterp3DView(std::shared_ptr< const NDArray > in, BoundaryConditionT bound=ZEROFLUX)
Definition: accessors.h:1141
npl::Vector3DView::set
void set(int64_t x, int64_t y, int64_t z, int64_t t, T v)
Gets value at array index and then casts to T.
Definition: accessors.h:735
npl::BSplineView::jacobianDet
double jacobianDet(int dir)
Computes the regularization term by integrating over the entire space for each knot. Thankfully integrals can be pre-computed (vConv, dvConv, ddvConv). See equations 66-71 in docs/bspline/fmri_dist_correct_2013-12-06.pdf.
Definition: accessors.h:2508
npl::NDIter::get
T get() const
Dereference operator.
Definition: iterators.h:992
npl::Vector3DConstView::operator()
virtual T operator()(double x=0, double y=0, double z=0, int64_t t=0)
Gets value at array index and then casts to T.
Definition: accessors.h:599
npl::NDConstView::operator[]
T operator[](const std::vector< int64_t > &index)
Gets value at array index and then casts to T.
Definition: accessors.h:445
npl::LanczosInterp3DView
The purpose of this class is to view an image as a continuous 3D+vector dimension image rather than a...
Definition: accessors.h:1866
npl::NDConstView::castget
T(* castget)(void *ptr)
Function pointer to the correct function for casting from the underlying type.
Definition: accessors.h:483
npl::NNInterpNDView::operator()
T operator()(const std::vector< float > &index)
Gets value at array index and then casts to T.
Definition: accessors.h:1314
npl::Vector3DConstView
The purpose of this class is to view an image as a 3D+vector dimension image rather than a 4+D image...
Definition: accessors.h:557
npl::LinInterpNDView
The purpose of this class is to view an image as a continuous ND image and to sample at a continuous ...
Definition: accessors.h:879
npl::NNInterpNDView
General purpose Nearest-Neighbor interpolator.
Definition: accessors.h:1275
npl::BSplineView::thinPlateEnergy
double thinPlateEnergy(size_t len, double *grad)
Definition: accessors.h:2408
npl::BSplineView::reconstructNotAligned
ptr< MRImage > reconstructNotAligned(ptr< const MRImage > input)
Definition: accessors.h:2284
npl::NNInterpNDView::m_boundmethod
BoundaryConditionT m_boundmethod
Definition: accessors.h:1488
npl::LinInterp3DView::operator()
T operator()(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
Gets value at array index and then casts to T.
Definition: accessors.h:1244
npl::NDView::castsetStatic
static void castsetStatic(void *ptr, const T &val)
This is a wrapper function that will be called to safely cast to the underlying type.
Definition: accessors.h:287
npl::BoundaryConditionT
BoundaryConditionT
Definition: mrimage.h:37
npl::LanczosInterpNDView::operator()
T operator()(double x=0, double y=0, double z=0, double t=0, double u=0, double v=0, double w=0)
Gets value at array index and then casts to T.
Definition: accessors.h:1673
npl::BSplineView::m_ras
bool m_ras
if true, then this assumes the inputs are RAS coordinates rather than indexes. Default is false ...
Definition: accessors.h:2643
npl::LanczosInterpNDView::m_ras
bool m_ras
if true, then this assumes the inputs are RAS coordinates rather than indexes. Default is false ...
Definition: accessors.h:1834
npl::NNInterpNDView::operator[]
T operator[](const std::vector< int64_t > &index)
Gets value at array index and then casts to T.
Definition: accessors.h:1365
npl::NDIter::goBegin
void goBegin()
Go to beginning of iteration.
Definition: iterators.h:1016
npl::NNInterp3DView::m_ras
bool m_ras
if true, then this assumes the inputs are RAS coordinates rather than indexes. Default is false ...
Definition: accessors.h:1633
npl::dB3kern
double dB3kern(double x)
3rd order B-Spline derivative, radius 2 [-2,2]
Definition: basic_functions.h:294
npl::BSplineView
This is a specialized viewer for computing the value of a cubic B-Spline interpolation from the param...
Definition: accessors.h:2049
npl::NDView::setArray
void setArray(ptr< NDArray > in)
Definition: accessors.h:73
npl::Pixel3DView::operator[]
T operator[](const std::vector< int64_t > &i)
Definition: accessors.h:546
npl::NDView::tlen
int64_t tlen()
Definition: accessors.h:258
npl::LinInterpNDView::operator()
T operator()(std::initializer_list< float > index)
Gets value at array index and then casts to T.
Definition: accessors.h:964
npl::LanczosInterpNDView
The purpose of this class is to view an image as a continuous 3D+vector dimension image rather than a...
Definition: accessors.h:1654
npl::Pixel3DView::set
void set(int64_t x, int64_t y, int64_t z, T v)
Gets value at array index and then casts to T.
Definition: accessors.h:533
npl::NDView::set
void set(int64_t index, T v)
Casts to the appropriate type then sets array at given index.
Definition: accessors.h:250
npl::NNInterpNDView::m_ras
bool m_ras
if true, then this assumes the inputs are RAS coordinates rather than indexes. Default is false ...
Definition: accessors.h:1494
npl::LinInterpNDView::operator()
T operator()(const std::vector< double > &index)
Gets value at array index and then casts to T.
Definition: accessors.h:935
npl::NDConstView::operator[]
T operator[](int64_t index)
Gets value linear position in array, then casts to T.
Definition: accessors.h:399
npl::NNInterpNDView::get
T get(const std::vector< int64_t > &index)
Gets value at array index and then casts to T.
Definition: accessors.h:1331
npl::NDConstView::tlen
int64_t tlen()
Definition: accessors.h:453
npl::LinInterpNDView::operator[]
T operator[](int64_t i)
Gets value at array index and then casts to T doesn't make sense for interpolation.
Definition: accessors.h:1124
npl::NDView::operator[]
T operator[](const std::vector< int64_t > &index)
Gets value at array index and then casts to T.
Definition: accessors.h:201
npl::NDConstView::NDConstView
NDConstView(std::shared_ptr< const NDArray > in)
Definition: accessors.h:327
npl::BSplineView::jacobianDet
double jacobianDet(int dir, size_t len, double *grad)
Computes the gradient of regularization for each of the knots. Thankfully integrals can be pre-comput...
Definition: accessors.h:2565
npl::LinInterpNDView::LinInterpNDView
LinInterpNDView(std::shared_ptr< const NDArray > in, BoundaryConditionT bound=ZEROFLUX)
Definition: accessors.h:882
npl::NNInterpNDView::operator()
T operator()(std::initializer_list< double > index)
Gets value at array index and then casts to T.
Definition: accessors.h:1394
npl::BSplineView::reconstructAligned
ptr< MRImage > reconstructAligned(ptr< const MRImage > input)
Definition: accessors.h:2292
npl::LanczosInterp3DView::operator[]
T operator[](int64_t i)
Gets value at array index and then casts to T.
Definition: accessors.h:2023
npl::NDIter::eof
bool eof() const
Are we one past the last element?
Definition: iterators.h:1031
npl::LanczosInterpNDView::operator[]
T operator[](const std::vector< int64_t > &index)
Gets value at array index and then casts to T.
Definition: accessors.h:1721
npl::linKern
double linKern(double x)
Definition: accessors.h:867
npl::lanczosKern
double lanczosKern(double x, double a)
Derivative of lanczos kernel with respect to x.
Definition: basic_functions.h:189
npl::ptr
std::shared_ptr< T > ptr
Make the shared_ptr name shorter...
Definition: npltypes.h:42
npl::NNInterp3DView
The purpose of this class is to view an image as a continuous 3D+vector dimension image rather than a...
Definition: accessors.h:1526
npl::LinInterpNDView::m_boundmethod
BoundaryConditionT m_boundmethod
Definition: accessors.h:1108
npl::NNInterp3DView::NNInterp3DView
NNInterp3DView(std::shared_ptr< NDArray > in, BoundaryConditionT bound=ZEROFLUX)
Definition: accessors.h:1529
npl::Pixel3DView::Pixel3DView
Pixel3DView(std::shared_ptr< NDArray > in)
Definition: accessors.h:498
npl::NDView::NDView
NDView(std::shared_ptr< NDArray > in)
Definition: accessors.h:66
npl::NNInterp3DView::operator()
T operator()(double x=0, double y=0, double z=0, int64_t t=0)
Gets value at array index and then casts to T.
Definition: accessors.h:1541
npl::BSplineView::createOverlay
void createOverlay(ptr< const MRImage > overlay, double bspace)
Definition: accessors.h:2059
npl::LanczosInterpNDView::operator[]
T operator[](int64_t index)
Gets value linear position in array, then casts to T.
Definition: accessors.h:1843
npl::NDView::parent
std::shared_ptr< NDArray > parent
Where to get the dat a from. Also the shared_ptr prevents dealloc.
Definition: accessors.h:290
npl::Pixel3DView::operator()
T operator()(int64_t x=0, int64_t y=0, int64_t z=0)
Gets value at array index and then casts to T.
Definition: accessors.h:507
npl::LanczosInterp3DView::operator()
T operator()(double x=0, double y=0, double z=0, int64_t t=0)
Gets value at array index and then casts to T.
Definition: accessors.h:1885
npl::NDView
This is a basic accessor class, which allows for accessing array data in the type specified by the te...
Definition: accessors.h:63
npl::NNInterp3DView::operator()
T operator()(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
Gets value at array index and then casts to T.
Definition: accessors.h:1611
npl::MAXDIM
const size_t MAXDIM
Defines the maximum supported dimension by image, used for stack-allocations.
Definition: ndarray.h:46
npl::LinInterpNDView::operator()
T operator()(double x=0, double y=0, double z=0, double t=0, double u=0, double v=0, double w=0, double q=0)
Gets value at array index and then casts to T.
Definition: accessors.h:903
npl::B3kern
double B3kern(double x)
3rd order B-Spline, radius 2 [-2,2]
Definition: basic_functions.h:251
npl::Pixel3DView
The purpose of this class is to view an image as a 3D image rather than a ND image. Therefore all dimensions above the third will be ignored and index 0 will be used.
Definition: accessors.h:495
npl::NDView::set
void set(const std::vector< int64_t > &index, T v)
Casts to the appropriate type then sets array at given index.
Definition: accessors.h:234
npl::NDConstView
This is a basic accessor class, which allows for accessing array data in the type specified by the te...
Definition: accessors.h:324
npl::LinInterp3DView
The purpose of this class is to view an image as a continuous 3D+vector dimension image rather than a...
Definition: accessors.h:1138
npl::NDView::operator[]
T operator[](int64_t index)
Gets value linear position in array, then casts to T.
Definition: accessors.h:155
basic_functions.h
npl::LinInterpNDView::m_ras
bool m_ras
if true, then this assumes the inputs are RAS coordinates rather than indexes. Default is false ...
Definition: accessors.h:1114
npl::LanczosInterp3DView::m_boundmethod
BoundaryConditionT m_boundmethod
Definition: accessors.h:2008
npl::BSplineView::thinPlateEnergy
double thinPlateEnergy()
Definition: accessors.h:2349
npl::Vector3DConstView::Vector3DConstView
Vector3DConstView(std::shared_ptr< const NDArray > in)
Definition: accessors.h:560
npl::Vector3DView
The purpose of this class is to view an image as a 3D+vector dimension image rather than a 4+D image...
Definition: accessors.h:697
npl::NNInterpNDView::operator()
T operator()(const std::vector< double > &index)
Gets value at array index and then casts to T.
Definition: accessors.h:1382
npl::LanczosInterpNDView::m_boundmethod
BoundaryConditionT m_boundmethod
Definition: accessors.h:1828
npl::NDConstView::parent
std::shared_ptr< const NDArray > parent
Where to get the dat a from. Also the shared_ptr prevents dealloc.
Definition: accessors.h:470
npl::LanczosInterp3DView::operator()
T operator()(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
Gets value at array index and then casts to T.
Definition: accessors.h:1895
npl::Vector3DConstView::operator()
virtual T operator()(int64_t x=0, int64_t y=0, int64_t z=0, int64_t t=0)
Gets value at array index and then casts to T.
Definition: accessors.h:571
npl::LanczosInterpNDView::setRadius
void setRadius(size_t rad)
Definition: accessors.h:1665
npl::NDView::castget
T(* castget)(void *ptr)
Function pointer to the correct function for casting from the underlying type.
Definition: accessors.h:303
npl::NDConstView::castgetStatic
static T castgetStatic(void *ptr)
This is a wrapper function that will be called to safely cast from the underlying type...
Definition: accessors.h:467
npl::NDIter
This class is used to iterate through an N-Dimensional array.
Definition: iterators.h:860
npl::LanczosInterpNDView::LanczosInterpNDView
LanczosInterpNDView(std::shared_ptr< const NDArray > in, BoundaryConditionT bound=ZEROFLUX)
Definition: accessors.h:1657
npl::LanczosInterp3DView::LanczosInterp3DView
LanczosInterp3DView(std::shared_ptr< const NDArray > in, BoundaryConditionT bound=ZEROFLUX)
Definition: accessors.h:1869