Line data Source code
1 : /** \file imagingArray.hpp
2 : * \brief Declares and defines a class for managing images
3 : * \ingroup imaging
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : #ifndef imagingArray_hpp
9 : #define imagingArray_hpp
10 :
11 : #pragma GCC system_header
12 : #include <Eigen/Dense>
13 :
14 : #include "../math/templateBLAS.hpp"
15 : #include "../math/ft/fftT.hpp"
16 :
17 : namespace mx
18 : {
19 : namespace wfp
20 : {
21 :
22 : template <typename Scalar>
23 : struct fftwAllocator
24 : {
25 0 : Scalar *alloc( size_t sz )
26 : {
27 0 : return (Scalar *)::fftw_malloc( sz * sizeof( Scalar ) );
28 : }
29 :
30 0 : void free( Scalar *ptr )
31 : {
32 0 : fftw_free( ptr );
33 0 : }
34 : };
35 :
36 : /// Test whether a type is an imagingArray by testing whether it has a typedef of "is_imagingArray"
37 : /** Used for compile-time determination of type
38 : * Example usage:
39 : * \code
40 : * bool is_eC = is_imagingArray<eigenCube<float> >; //Evaluates to true
41 : * bool is_not_eC = is_imagingArray<float>; //Evaluates to false
42 : * \endcode
43 : */
44 : // This was taken directly from the example at http://en.wikipedia.org/wiki/Substitution_failure_is_not_an_error
45 : template <typename T>
46 : struct is_imagingArray
47 : {
48 : // Types "yes" and "no" are guaranteed to have different sizes,
49 : // specifically sizeof(yes) == 1 and sizeof(no) == 2.
50 : typedef char yes[1];
51 : typedef char no[2];
52 :
53 : template <typename imageT>
54 : static yes &test( typename imageT::imagingArrayT * );
55 :
56 : template <typename>
57 : static no &test( ... );
58 :
59 : // If the "sizeof" of the result of calling test<T>(0) would be equal to sizeof(yes),
60 : // the first overload worked and T has a nested type named "is_mmatrix".
61 : static const bool value = sizeof( test<T>( 0 ) ) == sizeof( yes );
62 : };
63 :
64 : template <typename imageT, typename typeT, bool isScalar = is_imagingArray<typeT>::value>
65 : struct imagingArrayInplaceProduct
66 : {
67 : imagingArrayInplaceProduct()
68 : {
69 : static_assert( is_imagingArray<imageT>::value,
70 : "imagingArrayInplaceProduct template parameter 1 must be an imagingArray type" );
71 : }
72 :
73 : void operator()( imageT &im, typeT &alpha )
74 : {
75 : math::scal<typename imageT::Scalar>( im.cols() * im.rows(), (typeT)alpha, im.data(), 1 );
76 : }
77 : };
78 :
79 : template <typename imageT, typename typeT>
80 : struct imagingArrayInplaceProduct<imageT, typeT, true>
81 : {
82 0 : imagingArrayInplaceProduct()
83 : {
84 : static_assert( is_imagingArray<imageT>::value,
85 : "imagingArrayInplaceProduct template parameter 1 must be an imagingArray type" );
86 0 : }
87 :
88 0 : void operator()( imageT &im, typeT &alpha )
89 : {
90 : typedef typename imageT::Scalar scalarT;
91 : // hadp<scalarT>(im.cols()*im.rows(), alpha.data(), im.data());
92 :
93 0 : Eigen::Map<Eigen::Array<scalarT, -1, -1>> eigY( alpha.data(), alpha.cols(), alpha.rows() );
94 0 : Eigen::Map<Eigen::Array<scalarT, -1, -1>> eigX( im.data(), im.cols(), im.rows() );
95 :
96 0 : eigX *= eigY;
97 0 : }
98 : };
99 :
100 : template <typename imageT, typename typeT, bool isScalar = is_imagingArray<typeT>::value>
101 : struct imagingArrayInplaceDivision
102 : {
103 : imagingArrayInplaceDivision()
104 : {
105 : static_assert( is_imagingArray<imageT>::value,
106 : "imagingArrayInplaceDivision template parameter 1 must be an imagingArray type" );
107 : }
108 :
109 : void operator()( imageT im, typeT alpha )
110 : {
111 : math::scal<typename imageT::Scalar>( im.cols() * im.rows(), ( (typeT)( 1 ) ) / alpha, im.data(), 1 );
112 : }
113 : };
114 :
115 : template <typename imageT, typename typeT>
116 : struct imagingArrayInplaceDivision<imageT, typeT, true>
117 : {
118 : imagingArrayInplaceDivision()
119 : {
120 : static_assert( is_imagingArray<imageT>::value,
121 : "imagingArrayInplaceDivision template parameter 1 must be an imagingArray type" );
122 : }
123 :
124 : void operator()( imageT im, typeT alpha )
125 : {
126 : typedef typename imageT::Scalar scalarT;
127 : math::hadd<scalarT>( im.cols() * im.rows(), (scalarT)1.0, alpha.data(), 1, im.data(), 1 );
128 : }
129 : };
130 :
131 : template <typename imageT, typename typeT, bool isScalar = is_imagingArray<typeT>::value>
132 : struct imagingArrayInplaceAdd
133 : {
134 : imagingArrayInplaceAdd()
135 : {
136 : static_assert( is_imagingArray<imageT>::value,
137 : "imagingArrayInplaceAdd template parameter 1 must be an imagingArray type" );
138 : }
139 :
140 : void operator()( imageT &im, typeT &alpha )
141 : {
142 : typedef typename imageT::Scalar scalarT;
143 : Eigen::Map<Eigen::Array<scalarT, -1, -1>> eigX( im.data(), im.cols(), im.rows() );
144 : eigX += alpha;
145 : }
146 : };
147 :
148 : template <typename imageT, typename typeT>
149 : struct imagingArrayInplaceAdd<imageT, typeT, true>
150 : {
151 : imagingArrayInplaceAdd()
152 : {
153 : static_assert( is_imagingArray<imageT>::value,
154 : "imagingArrayInplaceAdd template parameter 1 must be an imagingArray type" );
155 : }
156 :
157 : void operator()( imageT &im, typeT &alpha )
158 : {
159 : typedef typename imageT::Scalar scalarT;
160 : // hadp<scalarT>(im.cols()*im.rows(), alpha.data(), im.data());
161 :
162 : Eigen::Map<Eigen::Array<scalarT, -1, -1>> eigY( alpha.data(), alpha.cols(), alpha.rows() );
163 : Eigen::Map<Eigen::Array<scalarT, -1, -1>> eigX( im.data(), im.cols(), im.rows() );
164 :
165 : eigX += eigY;
166 : }
167 : };
168 :
169 : template <typename _Scalar, class _allocatorT, int cudaGPU>
170 : class imagingArray;
171 :
172 : template <typename _Scalar, class _allocatorT>
173 : class imagingArray<_Scalar, _allocatorT, 0>
174 : {
175 : public:
176 : typedef bool imagingArrayT;
177 : typedef _Scalar Scalar;
178 : typedef _allocatorT allocatorT;
179 :
180 : protected:
181 : allocatorT allocator;
182 :
183 : Scalar *_data;
184 :
185 : bool _owner;
186 :
187 : int _cols;
188 : int _rows;
189 :
190 : public:
191 0 : imagingArray()
192 : {
193 0 : _data = 0;
194 0 : _owner = false;
195 0 : _cols = 0;
196 0 : _rows = 0;
197 0 : }
198 :
199 0 : imagingArray( int x, int y )
200 : {
201 0 : _data = 0;
202 0 : _owner = false;
203 0 : _cols = 0;
204 0 : _rows = 0;
205 :
206 0 : resize( x, y );
207 0 : }
208 :
209 : template <typename eigenT>
210 : explicit imagingArray( eigenT &im )
211 : {
212 : _data = im.data();
213 : _owner = false;
214 : _cols = im.cols();
215 : _rows = im.rows();
216 : }
217 :
218 0 : ~imagingArray()
219 : {
220 0 : free();
221 0 : }
222 :
223 0 : void resize( int szx, int szy )
224 : {
225 0 : if( szx == _cols && szy == _rows && _owner )
226 0 : return;
227 :
228 0 : free();
229 :
230 0 : _cols = szx;
231 0 : _rows = szy;
232 :
233 0 : _data = allocator.alloc( _cols * _rows );
234 0 : _owner = true;
235 : }
236 :
237 0 : void free()
238 : {
239 0 : if( _data && _owner )
240 : {
241 0 : allocator.free( _data );
242 0 : _data = 0;
243 0 : _cols = 0;
244 0 : _rows = 0;
245 : }
246 0 : }
247 :
248 0 : int cols() const
249 : {
250 0 : return _cols;
251 : }
252 :
253 0 : int rows() const
254 : {
255 0 : return _rows;
256 : }
257 :
258 0 : Scalar *data()
259 : {
260 0 : return _data;
261 : }
262 :
263 0 : Scalar *data() const
264 : {
265 0 : return _data;
266 : }
267 :
268 0 : Scalar operator()( int x, int y ) const
269 : {
270 0 : return _data[x + y * _cols];
271 : }
272 :
273 0 : Scalar &operator()( int x, int y )
274 : {
275 0 : return _data[x + y * _cols];
276 : }
277 :
278 : template <typename argT>
279 0 : void setConstant( argT val )
280 : {
281 0 : int N = _cols * _rows;
282 :
283 0 : Scalar dval = (Scalar)val;
284 0 : for( int i = 0; i < N; ++i )
285 0 : _data[i] = dval;
286 0 : }
287 :
288 0 : void setZero()
289 : {
290 0 : setConstant( (Scalar)0 );
291 0 : }
292 :
293 0 : imagingArray &operator=( const imagingArray &arr )
294 : {
295 0 : resize( arr.rows(), arr.cols() );
296 :
297 0 : for( int i = 0; i < _rows; ++i )
298 : {
299 0 : for( int j = 0; j < _cols; ++j )
300 : {
301 0 : operator()( i, j ) = arr( i, j );
302 : }
303 : }
304 :
305 0 : return ( *this );
306 : }
307 :
308 : /// In-place product. Works for both scalars and imagingArray types.
309 : template <typename typeT>
310 0 : imagingArray &operator*=( typeT &alpha )
311 : {
312 0 : imagingArrayInplaceProduct<imagingArray, typeT> prod;
313 0 : prod( *this, alpha );
314 0 : return *this;
315 : }
316 :
317 : /// In-place division. Works for both scalars and imagingArray types.
318 : template <typename typeT>
319 : imagingArray &operator/=( typeT alpha )
320 : {
321 : imagingArrayInplaceDivision<imagingArray, typeT> div;
322 : div( *this, alpha );
323 : return *this;
324 : }
325 :
326 : /// In-place add. Works for both scalars and imagingArray types.
327 : template <typename typeT>
328 : imagingArray &operator+=( typeT &alpha )
329 : {
330 : imagingArrayInplaceAdd<imagingArray, typeT> prod;
331 : prod( *this, alpha );
332 : return *this;
333 : }
334 :
335 0 : Scalar sum()
336 : {
337 0 : int N = _cols * _rows;
338 0 : Scalar _sum = 0;
339 :
340 0 : for( int i = 0; i < N; ++i )
341 0 : _sum += _data[i];
342 :
343 0 : return _sum;
344 : }
345 :
346 0 : Eigen::Map<Eigen::Array<Scalar, -1, -1>> eigenMap()
347 : {
348 0 : return Eigen::Map<Eigen::Array<Scalar, -1, -1>>( _data, _rows, _cols );
349 : }
350 : };
351 :
352 : template <typename scalarT>
353 : using imagingArrayT = imagingArray<scalarT, fftwAllocator<scalarT>, 0>;
354 :
355 : // typedef imagingArray<float,mx::fftwAllocator<float>,0> imagingArrayf;
356 : // typedef imagingArray<double,mx::fftwAllocator<double>,0> imagingArrayd;
357 : // typedef imagingArray<std::complex<float>,mx::fftwAllocator<std::complex<float> >,0> imagingArraycf;
358 :
359 : } // namespace wfp
360 : } // namespace mx
361 :
362 : #endif // imagingArray_hpp
|
