Line data Source code
1 : #include "DataHandlerCPU.h"
2 : #include "DataHandlerMapCPU.h"
3 : #include "Error.h"
4 : #include "TypeCasts.hpp"
5 :
6 : #include "DataDescriptor.h"
7 :
8 : #include <iostream>
9 :
10 : #if WITH_FFTW
11 : #define EIGEN_FFTW_DEFAULT
12 : #endif
13 : #include <unsupported/Eigen/FFT>
14 :
15 : namespace elsa
16 : {
17 :
18 : template <typename data_t>
19 0 : DataHandlerCPU<data_t>::DataHandlerCPU(index_t size)
20 0 : : _data(std::make_shared<DataVector_t>(size))
21 : {
22 0 : }
23 :
24 : template <typename data_t>
25 0 : DataHandlerCPU<data_t>::DataHandlerCPU(DataVector_t const& vector)
26 0 : : _data{std::make_shared<DataVector_t>(vector)}
27 : {
28 0 : }
29 :
30 : template <typename data_t>
31 0 : DataHandlerCPU<data_t>::DataHandlerCPU(const DataHandlerCPU<data_t>& other)
32 0 : : _data{other._data}, _associatedMaps{}
33 : {
34 0 : }
35 :
36 : template <typename data_t>
37 0 : DataHandlerCPU<data_t>::DataHandlerCPU(DataHandlerCPU<data_t>&& other)
38 0 : : _data{std::move(other._data)}, _associatedMaps{std::move(other._associatedMaps)}
39 : {
40 0 : for (auto& map : _associatedMaps)
41 0 : map->_dataOwner = this;
42 0 : }
43 :
44 : template <typename data_t>
45 0 : DataHandlerCPU<data_t>::~DataHandlerCPU()
46 : {
47 0 : for (auto& map : _associatedMaps)
48 0 : map->_dataOwner = nullptr;
49 0 : }
50 :
51 : template <typename data_t>
52 0 : index_t DataHandlerCPU<data_t>::getSize() const
53 : {
54 0 : return static_cast<index_t>(_data->size());
55 : }
56 :
57 : template <typename data_t>
58 0 : data_t& DataHandlerCPU<data_t>::operator[](index_t index)
59 : {
60 0 : detach();
61 0 : return (*_data)[index];
62 : }
63 :
64 : template <typename data_t>
65 0 : const data_t& DataHandlerCPU<data_t>::operator[](index_t index) const
66 : {
67 0 : return (*_data)[index];
68 : }
69 :
70 : template <typename data_t>
71 0 : data_t DataHandlerCPU<data_t>::dot(const DataHandler<data_t>& v) const
72 : {
73 0 : if (v.getSize() != getSize())
74 0 : throw InvalidArgumentError("DataHandlerCPU: dot product argument has wrong size");
75 :
76 : // use Eigen if the other handler is CPU or Map, otherwise use the slow fallback version
77 0 : if (auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&v)) {
78 0 : return _data->dot(*otherHandler->_data);
79 0 : } else if (auto otherHandler = downcast_safe<DataHandlerMapCPU<data_t>>(&v)) {
80 0 : return _data->dot(otherHandler->_map);
81 : } else {
82 0 : return this->slowDotProduct(v);
83 : }
84 : }
85 :
86 : template <typename data_t>
87 0 : GetFloatingPointType_t<data_t> DataHandlerCPU<data_t>::squaredL2Norm() const
88 : {
89 0 : return _data->squaredNorm();
90 : }
91 :
92 : template <typename data_t>
93 0 : GetFloatingPointType_t<data_t> DataHandlerCPU<data_t>::l2Norm() const
94 : {
95 0 : return _data->norm();
96 : }
97 :
98 : template <typename data_t>
99 0 : index_t DataHandlerCPU<data_t>::l0PseudoNorm() const
100 : {
101 : using FloatType = GetFloatingPointType_t<data_t>;
102 0 : return (_data->array().cwiseAbs() >= std::numeric_limits<FloatType>::epsilon()).count();
103 : }
104 :
105 : template <typename data_t>
106 0 : GetFloatingPointType_t<data_t> DataHandlerCPU<data_t>::l1Norm() const
107 : {
108 0 : return _data->array().abs().sum();
109 : }
110 :
111 : template <typename data_t>
112 0 : GetFloatingPointType_t<data_t> DataHandlerCPU<data_t>::lInfNorm() const
113 : {
114 0 : return _data->array().abs().maxCoeff();
115 : }
116 :
117 : template <typename data_t>
118 0 : data_t DataHandlerCPU<data_t>::sum() const
119 : {
120 0 : return _data->sum();
121 : }
122 :
123 : template <typename data_t>
124 0 : data_t DataHandlerCPU<data_t>::minElement() const
125 : {
126 : if constexpr (isComplex<data_t>) {
127 0 : throw LogicError("DataHandlerCPU: minElement of complex type not supported");
128 : } else {
129 0 : return _data->minCoeff();
130 : }
131 : }
132 :
133 : template <typename data_t>
134 0 : data_t DataHandlerCPU<data_t>::maxElement() const
135 : {
136 : if constexpr (isComplex<data_t>) {
137 0 : throw LogicError("DataHandlerCPU: maxElement of complex type not supported");
138 : } else {
139 0 : return _data->maxCoeff();
140 : }
141 : }
142 :
143 : template <typename data_t>
144 0 : DataHandler<data_t>& DataHandlerCPU<data_t>::fft(const DataDescriptor& source_desc,
145 : FFTNorm norm)
146 : {
147 0 : this->base_fft<true>(source_desc, norm);
148 0 : return *this;
149 : }
150 :
151 : template <typename data_t>
152 0 : DataHandler<data_t>& DataHandlerCPU<data_t>::ifft(const DataDescriptor& source_desc,
153 : FFTNorm norm)
154 : {
155 0 : this->base_fft<false>(source_desc, norm);
156 0 : return *this;
157 : }
158 :
159 : template <typename data_t>
160 0 : DataHandler<data_t>& DataHandlerCPU<data_t>::operator+=(const DataHandler<data_t>& v)
161 : {
162 0 : if (v.getSize() != getSize())
163 0 : throw InvalidArgumentError("DataHandler: addition argument has wrong size");
164 :
165 0 : detach();
166 :
167 : // use Eigen if the other handler is CPU or Map, otherwise use the slow fallback version
168 0 : if (auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&v)) {
169 0 : *_data += *otherHandler->_data;
170 0 : } else if (auto otherHandler = downcast_safe<DataHandlerMapCPU<data_t>>(&v)) {
171 0 : *_data += otherHandler->_map;
172 : } else {
173 0 : this->slowAddition(v);
174 : }
175 :
176 0 : return *this;
177 : }
178 :
179 : template <typename data_t>
180 0 : DataHandler<data_t>& DataHandlerCPU<data_t>::operator-=(const DataHandler<data_t>& v)
181 : {
182 0 : if (v.getSize() != getSize())
183 0 : throw InvalidArgumentError("DataHandler: subtraction argument has wrong size");
184 :
185 0 : detach();
186 :
187 : // use Eigen if the other handler is CPU or Map, otherwise use the slow fallback version
188 0 : if (auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&v)) {
189 0 : *_data -= *otherHandler->_data;
190 0 : } else if (auto otherHandler = downcast_safe<DataHandlerMapCPU<data_t>>(&v)) {
191 0 : *_data -= otherHandler->_map;
192 : } else {
193 0 : this->slowSubtraction(v);
194 : }
195 :
196 0 : return *this;
197 : }
198 :
199 : template <typename data_t>
200 0 : DataHandler<data_t>& DataHandlerCPU<data_t>::operator*=(const DataHandler<data_t>& v)
201 : {
202 0 : if (v.getSize() != getSize())
203 0 : throw InvalidArgumentError("DataHandler: multiplication argument has wrong size");
204 :
205 0 : detach();
206 :
207 : // use Eigen if the other handler is CPU or Map, otherwise use the slow fallback version
208 0 : if (auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&v)) {
209 0 : _data->array() *= otherHandler->_data->array();
210 0 : } else if (auto otherHandler = downcast_safe<DataHandlerMapCPU<data_t>>(&v)) {
211 0 : _data->array() *= otherHandler->_map.array();
212 : } else {
213 0 : this->slowMultiplication(v);
214 : }
215 :
216 0 : return *this;
217 : }
218 :
219 : template <typename data_t>
220 0 : DataHandler<data_t>& DataHandlerCPU<data_t>::operator/=(const DataHandler<data_t>& v)
221 : {
222 0 : if (v.getSize() != getSize())
223 0 : throw InvalidArgumentError("DataHandler: division argument has wrong size");
224 :
225 0 : detach();
226 :
227 : // use Eigen if the other handler is CPU or Map, otherwise use the slow fallback version
228 0 : if (auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&v)) {
229 0 : _data->array() /= otherHandler->_data->array();
230 0 : } else if (auto otherHandler = downcast_safe<DataHandlerMapCPU<data_t>>(&v)) {
231 0 : _data->array() /= otherHandler->_map.array();
232 : } else {
233 0 : this->slowDivision(v);
234 : }
235 :
236 0 : return *this;
237 : }
238 :
239 : template <typename data_t>
240 0 : DataHandlerCPU<data_t>& DataHandlerCPU<data_t>::operator=(const DataHandlerCPU<data_t>& v)
241 : {
242 0 : if (v.getSize() != getSize())
243 0 : throw InvalidArgumentError("DataHandler: assignment argument has wrong size");
244 :
245 0 : attach(v._data);
246 0 : return *this;
247 : }
248 :
249 : template <typename data_t>
250 0 : DataHandlerCPU<data_t>& DataHandlerCPU<data_t>::operator=(DataHandlerCPU<data_t>&& v)
251 : {
252 0 : if (v.getSize() != getSize())
253 0 : throw InvalidArgumentError("DataHandler: assignment argument has wrong size");
254 :
255 0 : attach(std::move(v._data));
256 :
257 0 : for (auto& map : v._associatedMaps)
258 0 : map->_dataOwner = this;
259 :
260 0 : _associatedMaps.splice(_associatedMaps.end(), std::move(v._associatedMaps));
261 :
262 : // make sure v no longer owns the object
263 0 : v._data.reset();
264 0 : return *this;
265 : }
266 :
267 : template <typename data_t>
268 0 : DataHandler<data_t>& DataHandlerCPU<data_t>::operator+=(data_t scalar)
269 : {
270 0 : detach();
271 0 : _data->array() += scalar;
272 0 : return *this;
273 : }
274 :
275 : template <typename data_t>
276 0 : DataHandler<data_t>& DataHandlerCPU<data_t>::operator-=(data_t scalar)
277 : {
278 0 : detach();
279 0 : _data->array() -= scalar;
280 0 : return *this;
281 : }
282 :
283 : template <typename data_t>
284 0 : DataHandler<data_t>& DataHandlerCPU<data_t>::operator*=(data_t scalar)
285 : {
286 0 : detach();
287 0 : _data->array() *= scalar;
288 0 : return *this;
289 : }
290 :
291 : template <typename data_t>
292 0 : DataHandler<data_t>& DataHandlerCPU<data_t>::operator/=(data_t scalar)
293 : {
294 0 : detach();
295 0 : _data->array() /= scalar;
296 0 : return *this;
297 : }
298 :
299 : template <typename data_t>
300 0 : DataHandler<data_t>& DataHandlerCPU<data_t>::operator=(data_t scalar)
301 : {
302 0 : detach();
303 0 : _data->setConstant(scalar);
304 0 : return *this;
305 : }
306 :
307 : template <typename data_t>
308 0 : std::unique_ptr<DataHandler<data_t>> DataHandlerCPU<data_t>::getBlock(index_t startIndex,
309 : index_t numberOfElements)
310 : {
311 0 : if (startIndex >= getSize() || numberOfElements > getSize() - startIndex)
312 0 : throw InvalidArgumentError("DataHandler: requested block out of bounds");
313 :
314 : return std::unique_ptr<DataHandlerMapCPU<data_t>>{
315 0 : new DataHandlerMapCPU{this, _data->data() + startIndex, numberOfElements}};
316 : }
317 :
318 : template <typename data_t>
319 : std::unique_ptr<const DataHandler<data_t>>
320 0 : DataHandlerCPU<data_t>::getBlock(index_t startIndex, index_t numberOfElements) const
321 : {
322 0 : if (startIndex >= getSize() || numberOfElements > getSize() - startIndex)
323 0 : throw InvalidArgumentError("DataHandler: requested block out of bounds");
324 :
325 : // using a const_cast here is fine as long as the DataHandlers never expose the internal
326 : // Eigen objects
327 0 : auto mutableThis = const_cast<DataHandlerCPU<data_t>*>(this);
328 0 : auto mutableData = const_cast<data_t*>(_data->data() + startIndex);
329 : return std::unique_ptr<const DataHandlerMapCPU<data_t>>{
330 0 : new DataHandlerMapCPU{mutableThis, mutableData, numberOfElements}};
331 : }
332 :
333 : template <typename data_t>
334 0 : DataHandlerCPU<data_t>* DataHandlerCPU<data_t>::cloneImpl() const
335 : {
336 0 : return new DataHandlerCPU<data_t>(*this);
337 : }
338 :
339 : template <typename data_t>
340 0 : bool DataHandlerCPU<data_t>::isEqual(const DataHandler<data_t>& other) const
341 : {
342 0 : if (const auto otherHandler = downcast_safe<DataHandlerMapCPU<data_t>>(&other)) {
343 :
344 0 : if (_data->size() != otherHandler->_map.size())
345 0 : return false;
346 :
347 0 : if (_data->data() != otherHandler->_map.data() && *_data != otherHandler->_map)
348 0 : return false;
349 :
350 0 : return true;
351 0 : } else if (const auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&other)) {
352 :
353 0 : if (_data->size() != otherHandler->_data->size())
354 0 : return false;
355 :
356 0 : if (_data->data() != otherHandler->_data->data() && *_data != *otherHandler->_data)
357 0 : return false;
358 :
359 0 : return true;
360 : } else
361 0 : return false;
362 : }
363 :
364 : template <typename data_t>
365 0 : void DataHandlerCPU<data_t>::assign(const DataHandler<data_t>& other)
366 : {
367 :
368 0 : if (const auto otherHandler = downcast_safe<DataHandlerMapCPU<data_t>>(&other)) {
369 0 : if (getSize() == otherHandler->_dataOwner->getSize()) {
370 0 : attach(otherHandler->_dataOwner->_data);
371 : } else {
372 0 : detachWithUninitializedBlock(0, getSize());
373 0 : *_data = otherHandler->_map;
374 : }
375 0 : } else if (const auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&other)) {
376 0 : attach(otherHandler->_data);
377 : } else
378 0 : this->slowAssign(other);
379 0 : }
380 :
381 : template <typename data_t>
382 0 : void DataHandlerCPU<data_t>::assign(DataHandler<data_t>&& other)
383 : {
384 0 : if (const auto otherHandler = downcast_safe<DataHandlerMapCPU<data_t>>(&other)) {
385 0 : if (getSize() == otherHandler->_dataOwner->getSize()) {
386 0 : attach(otherHandler->_dataOwner->_data);
387 : } else {
388 0 : detachWithUninitializedBlock(0, getSize());
389 0 : *_data = otherHandler->_map;
390 : }
391 0 : } else if (const auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&other)) {
392 0 : attach(std::move(otherHandler->_data));
393 :
394 0 : for (auto& map : otherHandler->_associatedMaps)
395 0 : map->_dataOwner = this;
396 :
397 0 : _associatedMaps.splice(_associatedMaps.end(), std::move(otherHandler->_associatedMaps));
398 :
399 : // make sure v no longer owns the object
400 0 : otherHandler->_data.reset();
401 : } else
402 0 : this->slowAssign(other);
403 0 : }
404 :
405 : template <typename data_t>
406 0 : void DataHandlerCPU<data_t>::detach()
407 : {
408 0 : if (_data.use_count() != 1) {
409 0 : #pragma omp barrier
410 : #pragma omp single
411 : {
412 0 : data_t* oldData = _data->data();
413 :
414 : // create deep copy of vector
415 0 : _data = std::make_shared<DataVector_t>(*_data);
416 :
417 : // modify all associated maps
418 0 : for (auto map : _associatedMaps)
419 0 : new (&map->_map) Eigen::Map<DataVector_t>(
420 0 : _data->data() + (map->_map.data() - oldData), map->getSize());
421 : }
422 : }
423 0 : }
424 :
425 : template <typename data_t>
426 0 : void DataHandlerCPU<data_t>::detachWithUninitializedBlock(index_t startIndex,
427 : index_t numberOfElements)
428 : {
429 0 : if (_data.use_count() != 1) {
430 : // allocate new vector
431 0 : auto newData = std::make_shared<DataVector_t>(getSize());
432 :
433 : // copy elements before start of block
434 0 : newData->head(startIndex) = _data->head(startIndex);
435 0 : newData->tail(getSize() - startIndex - numberOfElements) =
436 0 : _data->tail(getSize() - startIndex - numberOfElements);
437 :
438 : // modify all associated maps
439 0 : for (auto map : _associatedMaps)
440 0 : new (&map->_map) Eigen::Map<DataVector_t>(
441 0 : newData->data() + (map->_map.data() - _data->data()), map->getSize());
442 :
443 0 : _data = newData;
444 0 : }
445 0 : }
446 :
447 : template <typename data_t>
448 0 : void DataHandlerCPU<data_t>::attach(const std::shared_ptr<DataVector_t>& data)
449 : {
450 0 : data_t* oldData = _data->data();
451 :
452 : // shallow copy
453 0 : _data = data;
454 :
455 : // modify all associated maps
456 0 : for (auto& map : _associatedMaps)
457 0 : new (&map->_map) Eigen::Map<DataVector_t>(_data->data() + (map->_map.data() - oldData),
458 0 : map->getSize());
459 0 : }
460 :
461 : template <typename data_t>
462 0 : void DataHandlerCPU<data_t>::attach(std::shared_ptr<DataVector_t>&& data)
463 : {
464 0 : data_t* oldData = _data->data();
465 :
466 : // shallow copy
467 0 : _data = std::move(data);
468 :
469 : // modify all associated maps
470 0 : for (auto& map : _associatedMaps)
471 0 : new (&map->_map) Eigen::Map<DataVector_t>(_data->data() + (map->_map.data() - oldData),
472 0 : map->getSize());
473 0 : }
474 :
475 : template <typename data_t>
476 0 : typename DataHandlerCPU<data_t>::DataMap_t DataHandlerCPU<data_t>::accessData() const
477 : {
478 0 : return DataMap_t(&(_data->operator[](0)), getSize());
479 : }
480 :
481 : template <typename data_t>
482 0 : typename DataHandlerCPU<data_t>::DataMap_t DataHandlerCPU<data_t>::accessData()
483 : {
484 0 : detach();
485 0 : return DataMap_t(&(_data->operator[](0)), getSize());
486 : }
487 :
488 : template <typename data_t>
489 : template <bool is_forward>
490 0 : void DataHandlerCPU<data_t>::base_fft(const DataDescriptor& source_desc, FFTNorm norm)
491 : {
492 : if constexpr (isComplex<data_t>) {
493 :
494 : // now that we change this datahandler,
495 : // copy-on-write it.
496 0 : this->detach();
497 :
498 0 : const auto& src_shape = source_desc.getNumberOfCoefficientsPerDimension();
499 0 : const auto& src_dims = source_desc.getNumberOfDimensions();
500 :
501 0 : typename DataVector_t::Scalar* this_data = this->_data->data();
502 :
503 : // TODO: fftw variant
504 :
505 : // generalization of an 1D-FFT
506 : // walk over each dimenson and 1d-fft one 'line' of data
507 0 : for (index_t dim_idx = 0; dim_idx < src_dims; ++dim_idx) {
508 : // jumps in the data for the current dimension's data
509 : // dim_size[0] * dim_size[1] * ...
510 : // 1 for dim_idx == 0.
511 0 : const index_t stride = src_shape.head(dim_idx).prod();
512 :
513 : // number of coefficients for the current dimension
514 0 : const index_t dim_size = src_shape(dim_idx);
515 :
516 : // number of coefficients for the other dimensions
517 : // this is the number of 1d-ffts we'll do
518 : // e.g. shape=[2, 3, 4] and we do dim_idx=2 (=shape 4)
519 : // -> src_shape.prod() == 24 / 4 = 6 == 2*3
520 0 : const index_t other_dims_size = src_shape.prod() / dim_size;
521 :
522 : #ifndef EIGEN_FFTW_DEFAULT
523 : // when using eigen+fftw, this corrupts the memory, so don't parallelize.
524 : // error messages may include:
525 : // * double free or corruption (fasttop)
526 : // * malloc_consolidate(): unaligned fastbin chunk detected
527 0 : #pragma omp parallel for
528 : #endif
529 : // do all the 1d-ffts along the current dimensions axis
530 : for (index_t i = 0; i < other_dims_size; ++i) {
531 :
532 : index_t ray_start = i;
533 : // each time i is a multiple of stride,
534 : // jump forward the current+previous dimensions' shape product
535 : // (draw an indexed 3d cube to visualize this)
536 : ray_start += (stride * (dim_size - 1)) * ((i - (i % stride)) / stride);
537 :
538 : // this is one "ray" through the volume
539 : Eigen::Map<DataVector_t, Eigen::AlignmentType::Unaligned, Eigen::InnerStride<>>
540 : input_map(this_data + ray_start, dim_size, Eigen::InnerStride<>(stride));
541 :
542 : using inner_t = GetFloatingPointType_t<typename DataVector_t::Scalar>;
543 :
544 : Eigen::FFT<inner_t> fft_op;
545 :
546 : // disable any scaling in eigen - normally it does 1/n for ifft
547 : fft_op.SetFlag(Eigen::FFT<inner_t>::Flag::Unscaled);
548 :
549 : Eigen::Matrix<std::complex<inner_t>, Eigen::Dynamic, 1> fft_in{dim_size};
550 : Eigen::Matrix<std::complex<inner_t>, Eigen::Dynamic, 1> fft_out{dim_size};
551 :
552 : // eigen internally copies the fwd input matrix anyway if
553 : // it doesn't have stride == 1
554 : fft_in =
555 : input_map.block(0, 0, dim_size, 1).template cast<std::complex<inner_t>>();
556 :
557 : if (unlikely(dim_size == 1)) {
558 : // eigen kiss-fft crashes for size=1...
559 : fft_out = fft_in;
560 : } else {
561 : // arguments for in and out _must not_ be the same matrix!
562 : // they will corrupt wildly otherwise.
563 : if constexpr (is_forward) {
564 : fft_op.fwd(fft_out, fft_in);
565 : if (norm == FFTNorm::FORWARD) {
566 : fft_out /= dim_size;
567 : } else if (norm == FFTNorm::ORTHO) {
568 : fft_out /= std::sqrt(dim_size);
569 : }
570 : } else {
571 : fft_op.inv(fft_out, fft_in);
572 : if (norm == FFTNorm::BACKWARD) {
573 : fft_out /= dim_size;
574 : } else if (norm == FFTNorm::ORTHO) {
575 : fft_out /= std::sqrt(dim_size);
576 : }
577 : }
578 : }
579 :
580 : // we can't directly use the map as fft output,
581 : // since Eigen internally just uses the pointer to
582 : // the map's first element, and doesn't respect stride at all..
583 : input_map.block(0, 0, dim_size, 1) = fft_out.template cast<data_t>();
584 : }
585 : }
586 0 : } else {
587 0 : throw Error{"fft with non-complex input container not supported"};
588 : }
589 0 : }
590 :
591 : // ------------------------------------------
592 : // explicit template instantiation
593 : template class DataHandlerCPU<float>;
594 : template class DataHandlerCPU<complex<float>>;
595 : template class DataHandlerCPU<double>;
596 : template class DataHandlerCPU<complex<double>>;
597 : template class DataHandlerCPU<index_t>;
598 : } // namespace elsa
|
