utils.h
1 /***************************************************************************
2  mri_utils.h - description
3  -------------------
4  begin : Wed Feb 9 2005
5  copyright : (C) 2000-2015 by Thies H. Jochimsen
6  email :
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 #ifndef UTILS_H
18 #define UTILS_H
19 
20 #include <odindata/data.h>
21 
32 Data<float,1> unwrap_phase(const Data<float,1>& phase, int startindex=0);
33 
34 
36 
37 
41 template<typename T>
42 Array<T,1> matrix_product(const Array<T,2>& matrix, const Array<T,1>& vector) {
43  Log<OdinData> odinlog("","matrix_product");
44  int nrows=matrix.extent()(0);
45  int ncols=matrix.extent()(1);
46 
47  Array<T,1> result(nrows);
48  result=0;
49 
50  int vector_extent=vector.extent(0);
51  if(vector.extent(0)!=ncols) {
52  ODINLOG(odinlog,errorLog) << "size mismatch (vector_extent=" << vector_extent << ") != (ncols=" << ncols << ")" << STD_endl;
53  return result;
54  }
55 
56 
57  for(int icol=0; icol<ncols; icol++) {
58  for(int irow=0; irow<nrows; irow++) {
59  result(irow)+=matrix(irow,icol)*vector(icol);
60  }
61  }
62 
63  return result;
64 }
65 
67 
71 template<typename T>
72 Array<T,2> matrix_product(const Array<T,2>& matrix1, const Array<T,2>& matrix2) {
73  Log<OdinData> odinlog("","matrix_product");
74  int nrows=matrix1.extent()(0);
75  int ncols=matrix2.extent()(1);
76  ODINLOG(odinlog,normalDebug) << "nrows/ncols=" << nrows << "/" << ncols << STD_endl;
77 
78  Array<T,2> result(nrows,ncols);
79  result=0;
80 
81  if(matrix1.extent(1)!=matrix2.extent(0)) {
82  ODINLOG(odinlog,errorLog) << "size mismatch (matrix1=" << matrix1.shape() << ", matrix2=" << matrix2.shape() << ")" << STD_endl;
83  return result;
84  }
85 
86  int nprod=matrix1.extent(1);
87  ODINLOG(odinlog,normalDebug) << "nprod=" << nprod << STD_endl;
88  for(int irow=0; irow<nrows; irow++) {
89  for(int icol=0; icol<ncols; icol++) {
90 
91  T scalprod(0);
92  for(int iprod=0; iprod<nprod; iprod++) {
93  scalprod+= matrix1(irow,iprod)*matrix2(iprod,icol);
94  }
95  result(irow,icol)=scalprod;
96  }
97  }
98 
99  return result;
100 }
101 
103 
104 
108 template<typename T>
109 Array<T,1> vector_product(const Array<T,1>& u, const Array<T,1>& v) {
110  Log<OdinData> odinlog("","vector_product");
111  Array<T,1> result(3);
112  if(u.extent(0)!=3 || v.extent(0)!=3) {
113  ODINLOG(odinlog,errorLog) << "input size != 3" << STD_endl;
114  return result;
115  }
116  result(0)=u(1)*v(2)-u(2)*v(1);
117  result(1)=u(2)*v(0)-u(0)*v(2);
118  result(2)=u(0)*v(1)-u(1)*v(0);
119  return result;
120 }
121 
122 
124 
128 template<typename T, int N_rank>
129 bool operator == (const TinyVector<T,N_rank>& t1, const TinyVector<T,N_rank>& t2) {
130  return sum(abs(t1-t2))==0;
131 }
132 
134 
138 template<typename T, int N_rank>
139 bool operator != (const TinyVector<T,N_rank>& t1, const TinyVector<T,N_rank>& t2) {
140  return !(t1==t2);
141 }
142 
144 
148 template<typename T, int N_rows, int N_columns>
149 TinyVector<T,N_rows> operator * (const TinyMatrix<T, N_rows, N_columns>& matrix, const TinyVector<T, N_columns>& vector) {
150  TinyVector<T,N_rows> result;
151  result=0;
152 
153  for(int icol=0; icol<N_columns; icol++) {
154  for(int irow=0; irow<N_rows; irow++) {
155  result(irow)+=matrix(irow,icol)*vector(icol);
156  }
157  }
158 
159  return result;
160 }
161 
163 
167 template<typename T, int N_rank>
168 bool same_shape(const Array<T,N_rank>& a1, const Array<T,N_rank>& a2, const TinyVector<int,N_rank>& dimmask=1) {
169  for(int i=0; i<N_rank; i++) {
170  if(dimmask(i) && (a1.extent(i)!=a2.extent(i)) ) return false;
171  }
172  return true;
173 }
174 
175 
176 
178 
182 template <typename T>
183 bool check_range(T& val, T min, T max) {
184  bool in_range=true;
185  if(val<min) {val=min; in_range=false;}
186  if(val>max) {val=max; in_range=false;}
187  return in_range;
188 }
189 
190 
191 
192 
197 #endif
bool operator==(const TinyVector< T, N_rank > &t1, const TinyVector< T, N_rank > &t2)
Definition: utils.h:129
bool operator!=(const TinyVector< T, N_rank > &t1, const TinyVector< T, N_rank > &t2)
Definition: utils.h:139
Definition: tjlog.h:218
bool same_shape(const Array< T, N_rank > &a1, const Array< T, N_rank > &a2, const TinyVector< int, N_rank > &dimmask=1)
Definition: utils.h:168
Array< T, 1 > matrix_product(const Array< T, 2 > &matrix, const Array< T, 1 > &vector)
Definition: utils.h:42
Array< T, 1 > vector_product(const Array< T, 1 > &u, const Array< T, 1 > &v)
Definition: utils.h:109
TinyVector< T, N_rows > operator*(const TinyMatrix< T, N_rows, N_columns > &matrix, const TinyVector< T, N_columns > &vector)
Definition: utils.h:149
fvector phase(const cvector &cv)
Data< float, 1 > unwrap_phase(const Data< float, 1 > &phase, int startindex=0)
bool check_range(T &val, T min, T max)
Definition: utils.h:183