Line data Source code
1 : /** \file eigenCube.hpp
2 : * \brief An image cube with an Eigen API
3 : *
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : * \ingroup image_processing_files
7 : *
8 : */
9 :
10 : #ifndef eigenCube_hpp
11 : #define eigenCube_hpp
12 :
13 : #include <vector>
14 :
15 : #pragma GCC system_header
16 : #include <Eigen/Dense>
17 :
18 : #include "../error/exception.hpp"
19 : #include "../math/vectorUtils.hpp"
20 : #include "eigenImage.hpp"
21 : #include "imageUtils.hpp"
22 :
23 : namespace mx
24 : {
25 : namespace improc
26 : {
27 :
28 : /// An image cube with an Eigen-like API
29 : /** \ingroup eigen_image_processing
30 : */
31 : template <typename dataT>
32 : class eigenCube
33 : {
34 : public:
35 : typedef bool is_eigenCube;
36 : typedef dataT Scalar;
37 :
38 : typedef typename Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>::Index Index;
39 :
40 : typedef Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> imageRef;
41 :
42 : typedef Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic> imageT;
43 :
44 : protected:
45 : Index _rows;
46 : Index _cols;
47 : Index _planes;
48 :
49 : dataT *m_data{ nullptr };
50 :
51 : bool _owner;
52 :
53 : /// Validate a mask and good-pixel threshold for a masked combination.
54 : template <typename eigenCubeT>
55 : void
56 : validateMaskedCombination( const eigenCubeT &mask, /**< [in] mask cube whose dimensions must match this cube. */
57 : double minGoodFract /**< [in] inclusive minimum fraction of good pixels in [0, 1]. */
58 : ) const;
59 :
60 : /// Validate the number of weights for a weighted combination.
61 : void
62 : validateWeights( const std::vector<dataT> &weights /**< [in] weights whose count must match the number of planes. */
63 : ) const;
64 :
65 : /// Test whether a nonempty set of good pixels meets the requested fraction.
66 : bool hasEnoughGoodPixels( size_t goodPixels, /**< [in] number of unmasked pixels. */
67 : double minGoodFract /**< [in] inclusive minimum fraction of good pixels. */
68 : ) const;
69 :
70 : public:
71 : eigenCube();
72 :
73 : /// C'tor which will allocate space..
74 : /**
75 : */
76 : eigenCube( Index nrows, ///< [in] Number of rows in the cube
77 : Index ncols, ///< [in] Number of columns the cube
78 : Index nplanes ///< [in] Number of planes in the cube
79 : );
80 :
81 : /// C'tor taking an existing array as an argument.
82 : /** The existing array is used as is, as in a map, and ownership is not taken. You are
83 : * responsible for memory management, e.g. free-ing this array.
84 : */
85 : eigenCube( dataT *ndata, ///< [in] Allocated array with a cube of nrows x ncols x nplanes
86 : size_t nrows, ///< [in] Number of rows in the cube
87 : size_t ncols, ///< [in] Number of columns the cube
88 : size_t nplanes ///< [in] Number of planes in the cube
89 : );
90 :
91 : /// Copy a cube into independently owned storage.
92 : eigenCube( const eigenCube<dataT> &ec /**< [in] cube to copy. */ );
93 :
94 : /// Move a cube, transferring any owned storage.
95 : eigenCube( eigenCube<dataT> &&ec /**< [in,out] cube to move from. */ ) noexcept;
96 :
97 : ~eigenCube();
98 :
99 : void setZero();
100 :
101 : /// Copy a cube into independently owned storage.
102 : eigenCube<dataT> &operator=( const eigenCube<dataT> &ec /**< [in] cube to copy. */ );
103 :
104 : /// Move a cube, transferring any owned storage.
105 : eigenCube<dataT> &operator=( eigenCube<dataT> &&ec /**< [in,out] cube to move from. */ ) noexcept;
106 :
107 : void shallowCopy( eigenCube<dataT> &src, bool takeOwner = false );
108 :
109 : /// De-allocate and set all sizes to 0
110 : void clear();
111 :
112 : void resize( int r, int c, int p );
113 :
114 : void resize( int r, int c );
115 :
116 : dataT *data();
117 :
118 : const dataT *data() const;
119 :
120 : const Index rows() const;
121 :
122 : const Index cols() const;
123 :
124 : const Index planes() const;
125 :
126 : /// Returns a 2D Eigen::Eigen::Map pointed at the entire cube.
127 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> cube();
128 :
129 : /// Returns a 2D Eigen::Eigen::Map pointed at the specified image.
130 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> image( Index n /**< [in] the image number */ );
131 :
132 : /// Returns a 2D Eigen::Eigen::Map pointed at the specified image.
133 : const Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>
134 : image( Index n /**< [in] the image number */ ) const;
135 :
136 : /// Returns an Eigen::Eigen::Map-ed vector of the pixels at the given coordinate
137 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>,
138 : Eigen::Unaligned,
139 : Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>
140 : pixel( Index i, Index j );
141 :
142 : /// Return an Eigen::Eigen::Map of the cube where each image is a vector
143 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> asVectors();
144 :
145 : /// Calculate the covariance matrix of the images in the cube
146 : void Covar( Eigen::Matrix<dataT, Eigen::Dynamic, Eigen::Dynamic> &cv );
147 :
148 : /// Calculate the sum image of the cube
149 : /**
150 : * \tparam eigenT an Eigen-like type.
151 : */
152 : template <typename eigenT>
153 : void sum( eigenT &mim /**< [out] the resultant sum image. Is resized. */ );
154 :
155 : /// Calculate the mean image of the cube
156 : /**
157 : * \tparam eigenT an Eigen-like type.
158 : */
159 : template <typename eigenT>
160 : void mean( eigenT &mim /**< [out] the resultant mean image. Is resized. */ );
161 :
162 : /// Calculate the mean image of the cube with a mask.
163 : /**
164 : * A nonempty pixel sample is accepted when its good-pixel fraction is greater than or equal to
165 : * \p minGoodFract. Rejected pixels are set to invalidNumber<Scalar>().
166 : *
167 : * \tparam eigenT an Eigen-like type.
168 : * \tparam eigenCubeT an eigenCube type.
169 : *
170 : * \throws mx::exception with error_t::sizeerr if the mask dimensions do not match the cube.
171 : * \throws mx::exception with error_t::invalidarg if \p minGoodFract is not finite or is outside [0, 1].
172 : */
173 : template <typename eigenT, typename eigenCubeT>
174 : void mean( eigenT &mim, ///< [out] the resultant mean image. Is resized.
175 : eigenCubeT &mask, /**< [in] a mask cube. Only pixels with value 1 are included in the
176 : mean calculation. */
177 : double minGoodFract = 0.0 /**< [in] inclusive minimum fraction of good pixels in [0, 1]. */
178 : );
179 :
180 : /// Calculate the weighted mean image of the cube
181 : /**
182 : * \tparam eigenT an Eigen-like type.
183 : *
184 : * \throws mx::exception with error_t::sizeerr if the weight count does not match the cube planes.
185 : */
186 : template <typename eigenT>
187 : void mean( eigenT &mim, ///< [out] the resultant mean image. Is resized.
188 : std::vector<dataT> &weights ///< [in] a vector of weights to use for calculating the mean
189 : );
190 :
191 : /// Calculate the weighted mean image of the cube, with a mask cube.
192 : /**
193 : * A nonempty pixel sample is accepted when its good-pixel fraction is greater than or equal to
194 : * \p minGoodFract. Rejected pixels are set to invalidNumber<Scalar>().
195 : *
196 : * \tparam eigenT an Eigen-like type.
197 : * \tparam eigenCubeT an eigenCube type.
198 : *
199 : * \throws mx::exception with error_t::sizeerr if the mask dimensions or weight count do not match the cube.
200 : * \throws mx::exception with error_t::invalidarg if \p minGoodFract is not finite or is outside [0, 1].
201 : */
202 : template <typename eigenT, typename eigenCubeT>
203 : void mean( eigenT &mim, ///< [out] the resultant mean image. Is resized.
204 : std::vector<dataT> &weights, ///< [in] a vector of weights to use for calculating the mean
205 : eigenCubeT &mask, ///< [in] a mask cube. Only pixels with value 1 are included in the mean calculation.
206 : double minGoodFract = 0.0 ///< [in] inclusive minimum fraction of good pixels in [0, 1].
207 : );
208 :
209 : /// Calculate the median image of the cube
210 : /**
211 : * \tparam eigenT an Eigen-like type.
212 : */
213 : template <typename eigenT>
214 : void median( eigenT &mim /**< [out] the resultant median image. Is resized. */ );
215 :
216 : /// Calculate the median image of the cube with a mask.
217 : /**
218 : * A nonempty pixel sample is accepted when its good-pixel fraction is greater than or equal to
219 : * \p minGoodFract. Rejected pixels are set to invalidNumber<Scalar>().
220 : *
221 : * \tparam eigenT an Eigen-like type.
222 : * \tparam eigenCubeT an eigenCube type.
223 : *
224 : * \throws mx::exception with error_t::sizeerr if the mask dimensions do not match the cube.
225 : * \throws mx::exception with error_t::invalidarg if \p minGoodFract is not finite or is outside [0, 1].
226 : */
227 : template <typename eigenT, typename eigenCubeT>
228 : void median( eigenT &mim, ///< [out] the resultant median image. Is resized.
229 : eigenCubeT &mask, ///< [in] a mask cube. Only pixels with value 1 are included.
230 : double minGoodFract = 0.0 /**< [in] the inclusive minimum fraction of good pixels. */
231 : );
232 :
233 : /// Calculate the sigma clipped mean image of the cube
234 : /**
235 : * \tparam eigenT an Eigen-like type.
236 : */
237 : template <typename eigenT>
238 : void sigmaMean( eigenT &mim, ///< [out] the resultant mean image
239 : Scalar sigma ///< [in] the sigma value at which to clip.
240 : );
241 :
242 : /// Calculate the sigma clipped mean image of the cube, with a mask cube
243 : /**
244 : * A nonempty pixel sample is accepted when its good-pixel fraction is greater than or equal to
245 : * \p minGoodFract. Rejected pixels are set to invalidNumber<Scalar>().
246 : *
247 : * \tparam eigenT an Eigen-like type.
248 : * \tparam eigenCubeT an eigenCube type.
249 : *
250 : * \throws mx::exception with error_t::sizeerr if the mask dimensions do not match the cube.
251 : * \throws mx::exception with error_t::invalidarg if \p minGoodFract is not finite or is outside [0, 1].
252 : */
253 : template <typename eigenT, typename eigenCubeT>
254 : void
255 : sigmaMean( eigenT &mim, ///< [out] the resultant mean image. Is resized.
256 : eigenCubeT &mask, ///< [in] a mask cube. Only pixels with value 1 are included in the mean calculation.
257 : Scalar sigma, ///< [in] the sigma value at which to clip.
258 : double minGoodFract = 0.0 ///< [in] inclusive minimum fraction of good pixels in [0, 1].
259 : );
260 :
261 : /// Calculate the sigma clipped weighted mean image of the cube
262 : /**
263 : * \tparam eigenT an Eigen-like type.
264 : *
265 : * \throws mx::exception with error_t::sizeerr if the weight count does not match the cube planes.
266 : */
267 : template <typename eigenT>
268 : void sigmaMean( eigenT &mim, ///< [out] the resultant mean image. Is resized.
269 : std::vector<dataT> &weights, ///< [in] a vector of weights to use for calculating the mean
270 : Scalar sigma ///< [in] the sigma value at which to clip.
271 : );
272 :
273 : /// Calculate the sigma clipped weighted mean image of the cube, with a mask cube.
274 : /**
275 : * A nonempty pixel sample is accepted when its good-pixel fraction is greater than or equal to
276 : * \p minGoodFract. Rejected pixels are set to invalidNumber<Scalar>().
277 : *
278 : * \tparam eigenT an Eigen-like type.
279 : * \tparam eigenCubeT an eigenCube type.
280 : *
281 : * \throws mx::exception with error_t::sizeerr if the mask dimensions or weight count do not match the cube.
282 : * \throws mx::exception with error_t::invalidarg if \p minGoodFract is not finite or is outside [0, 1].
283 : */
284 : template <typename eigenT, typename eigenCubeT>
285 : void
286 : sigmaMean( eigenT &mim, ///< [out] the resultant mean image. Is resized.
287 : std::vector<dataT> &weights, ///< [in] a vector of weights to use for calculating the mean
288 : eigenCubeT &mask, ///< [in] a mask cube. Only pixels with value 1 are included in the mean calculation.
289 : Scalar sigma, ///< [in] the sigma value at which to clip.
290 : double minGoodFract = 0.0 ///< [in] inclusive minimum fraction of good pixels in [0, 1].
291 : );
292 : };
293 :
294 : template <typename dataT>
295 : template <typename eigenCubeT>
296 20 : void eigenCube<dataT>::validateMaskedCombination( const eigenCubeT &mask, double minGoodFract ) const
297 : {
298 20 : if( mask.rows() != _rows || mask.cols() != _cols || mask.planes() != _planes )
299 : {
300 3 : throw mx::exception( error_t::sizeerr, "mask dimensions must match cube dimensions" );
301 : }
302 :
303 19 : if( !math::isFinite( minGoodFract ) || minGoodFract < 0.0 || minGoodFract > 1.0 )
304 : {
305 12 : throw mx::exception( error_t::invalidarg, "minGoodFract must be finite and in [0, 1]" );
306 : }
307 15 : }
308 :
309 : template <typename dataT>
310 6 : void eigenCube<dataT>::validateWeights( const std::vector<dataT> &weights ) const
311 : {
312 6 : if( weights.size() != static_cast<size_t>( _planes ) )
313 : {
314 6 : throw mx::exception( error_t::sizeerr, "weights size must match cube planes" );
315 : }
316 4 : }
317 :
318 : template <typename dataT>
319 15 : bool eigenCube<dataT>::hasEnoughGoodPixels( size_t goodPixels, double minGoodFract ) const
320 : {
321 15 : return goodPixels > 0 && static_cast<double>( goodPixels ) >= minGoodFract * static_cast<double>( _planes );
322 : }
323 :
324 : template <typename dataT>
325 35 : eigenCube<dataT>::eigenCube()
326 : {
327 35 : _rows = 0;
328 35 : _cols = 0;
329 35 : _planes = 0;
330 35 : _owner = false;
331 35 : }
332 :
333 : template <typename dataT>
334 17 : eigenCube<dataT>::eigenCube( Index nrows, Index ncols, Index nplanes )
335 : {
336 17 : _rows = nrows;
337 17 : _cols = ncols;
338 17 : _planes = nplanes;
339 :
340 17 : m_data = new dataT[_rows * _cols * _planes];
341 17 : _owner = true;
342 17 : }
343 :
344 : template <typename dataT>
345 1 : eigenCube<dataT>::eigenCube( dataT *ndata, size_t nrows, size_t ncols, size_t nplanes )
346 : {
347 1 : _rows = nrows;
348 1 : _cols = ncols;
349 1 : _planes = nplanes;
350 :
351 1 : m_data = ndata;
352 1 : _owner = false;
353 1 : }
354 :
355 : template <typename dataT>
356 1 : eigenCube<dataT>::eigenCube( const eigenCube<dataT> &ec ) : eigenCube()
357 : {
358 1 : *this = ec;
359 1 : }
360 :
361 : template <typename dataT>
362 1 : eigenCube<dataT>::eigenCube( eigenCube<dataT> &&ec ) noexcept
363 1 : : _rows( ec._rows ), _cols( ec._cols ), _planes( ec._planes ), m_data( ec.m_data ), _owner( ec._owner )
364 : {
365 1 : ec._rows = 0;
366 1 : ec._cols = 0;
367 1 : ec._planes = 0;
368 1 : ec.m_data = nullptr;
369 1 : ec._owner = false;
370 1 : }
371 :
372 : template <typename dataT>
373 54 : eigenCube<dataT>::~eigenCube()
374 : {
375 54 : if( _owner && m_data )
376 : {
377 44 : delete[] m_data;
378 : }
379 54 : }
380 :
381 : template <typename dataT>
382 7 : void eigenCube<dataT>::setZero()
383 : {
384 7 : int N = _rows * _cols * _planes;
385 :
386 47 : for( int i = 0; i < N; ++i )
387 40 : m_data[i] = ( (dataT)0 );
388 7 : }
389 :
390 : template <typename dataT>
391 4 : eigenCube<dataT> &eigenCube<dataT>::operator=( const eigenCube<dataT> &ec )
392 : {
393 4 : if( this == &ec )
394 : {
395 1 : return *this;
396 : }
397 :
398 3 : resize( ec.rows(), ec.cols(), ec.planes() );
399 :
400 3 : int N = _rows * _cols * _planes;
401 :
402 4194315 : for( int i = 0; i < N; ++i )
403 4194312 : m_data[i] = ec.m_data[i];
404 :
405 3 : return *this;
406 : }
407 :
408 : template <typename dataT>
409 2 : eigenCube<dataT> &eigenCube<dataT>::operator=( eigenCube<dataT> &&ec ) noexcept
410 : {
411 2 : if( this == &ec )
412 : {
413 1 : return *this;
414 : }
415 :
416 1 : clear();
417 :
418 1 : _rows = ec._rows;
419 1 : _cols = ec._cols;
420 1 : _planes = ec._planes;
421 1 : m_data = ec.m_data;
422 1 : _owner = ec._owner;
423 :
424 1 : ec._rows = 0;
425 1 : ec._cols = 0;
426 1 : ec._planes = 0;
427 1 : ec.m_data = nullptr;
428 1 : ec._owner = false;
429 :
430 1 : return *this;
431 : }
432 :
433 : template <typename dataT>
434 3 : void eigenCube<dataT>::shallowCopy( eigenCube<dataT> &src, bool takeOwner )
435 : {
436 3 : if( this == &src )
437 : {
438 1 : return;
439 : }
440 :
441 2 : clear();
442 :
443 2 : _rows = src._rows;
444 2 : _cols = src._cols;
445 2 : _planes = src._planes;
446 2 : m_data = src.m_data;
447 :
448 2 : if( takeOwner == true )
449 : {
450 1 : _owner = src._owner;
451 1 : src._rows = 0;
452 1 : src._cols = 0;
453 1 : src._planes = 0;
454 1 : src.m_data = nullptr;
455 1 : src._owner = false;
456 : }
457 : else
458 : {
459 1 : _owner = false;
460 : }
461 : }
462 :
463 : template <typename dataT>
464 38 : void eigenCube<dataT>::clear()
465 : {
466 38 : if( _owner && m_data )
467 : {
468 7 : delete[] m_data;
469 : }
470 :
471 38 : _rows = 0;
472 38 : _cols = 0;
473 38 : _planes = 0;
474 38 : m_data = nullptr;
475 38 : _owner = false;
476 38 : }
477 :
478 : template <typename dataT>
479 34 : void eigenCube<dataT>::resize( int r, int c, int p )
480 : {
481 34 : clear();
482 :
483 34 : _rows = r;
484 34 : _cols = c;
485 34 : _planes = p;
486 :
487 12714018 : m_data = new dataT[_rows * _cols * _planes];
488 34 : _owner = true;
489 34 : }
490 :
491 : template <typename dataT>
492 1 : void eigenCube<dataT>::resize( int r, int c )
493 : {
494 1 : resize( r, c, 1 );
495 1 : }
496 :
497 : template <typename dataT>
498 814 : dataT *eigenCube<dataT>::data()
499 : {
500 814 : return m_data;
501 : }
502 :
503 : template <typename dataT>
504 3 : const dataT *eigenCube<dataT>::data() const
505 : {
506 3 : return m_data;
507 : }
508 :
509 : template <typename dataT>
510 18303650 : const typename eigenCube<dataT>::Index eigenCube<dataT>::rows() const
511 : {
512 18303650 : return _rows;
513 : }
514 :
515 : template <typename dataT>
516 13976587 : const typename eigenCube<dataT>::Index eigenCube<dataT>::cols() const
517 : {
518 13976587 : return _cols;
519 : }
520 :
521 : template <typename dataT>
522 198454 : const typename eigenCube<dataT>::Index eigenCube<dataT>::planes() const
523 : {
524 198454 : return _planes;
525 : }
526 :
527 : template <typename dataT>
528 3 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> eigenCube<dataT>::cube()
529 : {
530 3 : return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>( m_data, _rows * _cols, _planes );
531 : }
532 :
533 : template <typename dataT>
534 50493946 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> eigenCube<dataT>::image( Index n )
535 : {
536 50493946 : return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>( m_data + n * _rows * _cols, _rows, _cols );
537 : }
538 :
539 : template <typename dataT>
540 4325376 : const Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> eigenCube<dataT>::image( Index n ) const
541 : {
542 4325376 : return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>( m_data + n * _rows * _cols, _rows, _cols );
543 : }
544 :
545 : template <typename dataT>
546 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>,
547 : Eigen::Unaligned,
548 : Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>
549 108 : eigenCube<dataT>::pixel( Index i, Index j )
550 : {
551 : return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>,
552 : Eigen::Unaligned,
553 : Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>(
554 108 : m_data + j * _rows + i,
555 108 : _planes,
556 : 1,
557 216 : Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>( 0, _rows * _cols ) );
558 : }
559 :
560 : template <typename dataT>
561 3 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> eigenCube<dataT>::asVectors()
562 : {
563 3 : return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>( m_data, _rows * _cols, _planes );
564 : }
565 :
566 : template <typename dataT>
567 1 : void eigenCube<dataT>::Covar( Eigen::Matrix<dataT, Eigen::Dynamic, Eigen::Dynamic> &cv )
568 : {
569 1 : cv = asVectors().matrix().transpose() * asVectors().matrix();
570 1 : }
571 :
572 : template <typename dataT>
573 : template <typename eigenT>
574 : void eigenCube<dataT>::sum( eigenT &mim )
575 : {
576 : mim.resize( _rows, _cols );
577 :
578 : #pragma omp parallel for schedule( static, 1 ) num_threads( Eigen::nbThreads() )
579 : for( Index i = 0; i < _rows; ++i )
580 : {
581 : for( Index j = 0; j < _cols; ++j )
582 : {
583 : mim( i, j ) = pixel( i, j ).sum();
584 : }
585 : }
586 : }
587 :
588 : template <typename dataT>
589 : template <typename eigenT>
590 2 : void eigenCube<dataT>::mean( eigenT &mim )
591 : {
592 2 : mim.resize( _rows, _cols );
593 :
594 2 : #pragma omp parallel for schedule( static, 1 ) num_threads( Eigen::nbThreads() )
595 : for( Index i = 0; i < _rows; ++i )
596 : {
597 : for( Index j = 0; j < _cols; ++j )
598 : {
599 : mim( i, j ) = pixel( i, j ).mean();
600 : }
601 : }
602 2 : }
603 :
604 : template <typename dataT>
605 : template <typename eigenT, typename eigenCubeT>
606 8 : void eigenCube<dataT>::mean( eigenT &mim, eigenCubeT &mask, double minGoodFract )
607 : {
608 8 : validateMaskedCombination( mask, minGoodFract );
609 :
610 6 : mim.resize( _rows, _cols );
611 :
612 6 : #pragma omp parallel
613 : {
614 : std::vector<dataT> work;
615 :
616 : #pragma omp for
617 : for( Index i = 0; i < _rows; ++i )
618 : {
619 : for( Index j = 0; j < _cols; ++j )
620 : {
621 :
622 : work.clear();
623 :
624 : for( Index k = 0; k < _planes; ++k )
625 : {
626 : if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
627 : {
628 : work.push_back( ( pixel( i, j ) )( k, 0 ) );
629 : }
630 : }
631 :
632 : if( hasEnoughGoodPixels( work.size(), minGoodFract ) )
633 : {
634 : mim( i, j ) = math::vectorMean( work );
635 : }
636 : else
637 : {
638 : mim( i, j ) = invalidNumber<Scalar>();
639 : }
640 : }
641 : }
642 : }
643 6 : }
644 :
645 : template <typename dataT>
646 : template <typename eigenT>
647 2 : void eigenCube<dataT>::mean( eigenT &mim, std::vector<dataT> &weights )
648 : {
649 2 : validateWeights( weights );
650 :
651 1 : mim.resize( _rows, _cols );
652 :
653 1 : #pragma omp parallel num_threads( Eigen::nbThreads() )
654 : {
655 : std::vector<Scalar> work;
656 :
657 : #pragma omp for schedule( static, 10 )
658 : for( Index i = 0; i < _rows; ++i )
659 : {
660 : for( Index j = 0; j < _cols; ++j )
661 : {
662 : work.resize( _planes ); // work could be smaller after sigmaMean
663 : // int ii=0;
664 : for( int k = 0; k < _planes; ++k )
665 : work[k] = ( pixel( i, j ) )( k, 0 );
666 :
667 : mim( i, j ) = math::vectorMean( work, weights );
668 : }
669 : }
670 : }
671 1 : }
672 :
673 : template <typename dataT>
674 : template <typename eigenT, typename eigenCubeT>
675 1 : void eigenCube<dataT>::mean( eigenT &mim, std::vector<dataT> &weights, eigenCubeT &mask, double minGoodFract )
676 : {
677 1 : validateMaskedCombination( mask, minGoodFract );
678 1 : validateWeights( weights );
679 :
680 1 : mim.resize( _rows, _cols );
681 :
682 1 : #pragma omp parallel num_threads( Eigen::nbThreads() )
683 : {
684 : std::vector<Scalar> work, wwork;
685 :
686 : #pragma omp for schedule( static, 10 )
687 : for( Index i = 0; i < _rows; ++i )
688 : {
689 : for( Index j = 0; j < _cols; ++j )
690 : {
691 : work.clear();
692 : wwork.clear();
693 :
694 : // int ii=0;
695 :
696 : for( Index k = 0; k < _planes; ++k )
697 : {
698 : if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
699 : {
700 : work.push_back( ( pixel( i, j ) )( k, 0 ) );
701 : wwork.push_back( weights[k] );
702 : }
703 : }
704 : if( hasEnoughGoodPixels( work.size(), minGoodFract ) )
705 : {
706 : mim( i, j ) = math::vectorMean( work, wwork );
707 : }
708 : else
709 : mim( i, j ) = invalidNumber<Scalar>();
710 : }
711 : }
712 : }
713 1 : }
714 :
715 : template <typename dataT>
716 : template <typename eigenT>
717 2 : void eigenCube<dataT>::median( eigenT &mim )
718 : {
719 2 : mim.resize( _rows, _cols );
720 :
721 2 : #pragma omp parallel for schedule( static, 10 ) num_threads( Eigen::nbThreads() )
722 : for( Index i = 0; i < _rows; ++i )
723 : {
724 : // Eigen::Eigen::Array<Scalar, Eigen::Eigen::Dynamic, Eigen::Eigen::Dynamic> work;
725 : std::vector<Scalar> work;
726 : for( Index j = 0; j < _cols; ++j )
727 : {
728 : mim( i, j ) = imageMedian( pixel( i, j ), &work );
729 : }
730 : }
731 2 : }
732 :
733 : template <typename dataT>
734 : template <typename eigenT, typename eigenCubeT>
735 6 : void eigenCube<dataT>::median( eigenT &mim, eigenCubeT &mask, double minGoodFract )
736 : {
737 6 : validateMaskedCombination( mask, minGoodFract );
738 :
739 5 : mim.resize( _rows, _cols );
740 :
741 5 : #pragma omp parallel
742 : {
743 : std::vector<Scalar> work;
744 :
745 : #pragma omp for
746 : for( Index i = 0; i < _rows; ++i )
747 : {
748 : for( Index j = 0; j < _cols; ++j )
749 : {
750 : work.clear();
751 :
752 : for( Index k = 0; k < _planes; ++k )
753 : {
754 : if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
755 : {
756 : work.push_back( ( pixel( i, j ) )( k, 0 ) );
757 : }
758 : }
759 :
760 : if( hasEnoughGoodPixels( work.size(), minGoodFract ) )
761 : {
762 : mim( i, j ) = math::vectorMedianInPlace( work );
763 : }
764 : else
765 : {
766 : mim( i, j ) = invalidNumber<Scalar>();
767 : }
768 : }
769 : }
770 : }
771 5 : }
772 :
773 : template <typename dataT>
774 : template <typename eigenT>
775 1 : void eigenCube<dataT>::sigmaMean( eigenT &mim, dataT sigma )
776 : {
777 1 : mim.resize( _rows, _cols );
778 :
779 : // clang-format off
780 1 : #pragma omp parallel //clang-format on
781 : {
782 : std::vector<Scalar> work;
783 :
784 : // clang-format off
785 : #pragma omp for // clang-format off
786 : for( Index i = 0; i < _rows; ++i )
787 : {
788 : for( Index j = 0; j < _cols; ++j )
789 : {
790 : work.resize( _planes );
791 : for( int k = 0; k < _planes; ++k )
792 : {
793 : work[k] = ( pixel( i, j ) )( k, 0 );
794 : }
795 :
796 : mim( i, j ) = math::vectorSigmaMean( work, sigma );
797 : }
798 : }
799 : }
800 1 : }
801 :
802 : template <typename dataT>
803 : template <typename eigenT, typename eigenCubeT>
804 3 : void eigenCube<dataT>::sigmaMean( eigenT &mim, eigenCubeT &mask, dataT sigma, double minGoodFract )
805 : {
806 3 : validateMaskedCombination( mask, minGoodFract );
807 :
808 2 : mim.resize( _rows, _cols );
809 :
810 2 : #pragma omp parallel
811 : {
812 : std::vector<Scalar> work;
813 :
814 : #pragma omp for
815 : for( Index i = 0; i < _rows; ++i )
816 : {
817 : for( Index j = 0; j < _cols; ++j )
818 : {
819 : work.clear();
820 : for( Index k = 0; k < _planes; ++k )
821 : {
822 : if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
823 : {
824 : work.push_back( ( pixel( i, j ) )( k, 0 ) );
825 : }
826 : }
827 :
828 : if( hasEnoughGoodPixels( work.size(), minGoodFract ) )
829 : {
830 : mim( i, j ) = math::vectorSigmaMean( work, sigma );
831 : }
832 : else
833 : {
834 : mim( i, j ) = invalidNumber<Scalar>();
835 : }
836 : }
837 : }
838 : }
839 2 : }
840 :
841 : template <typename dataT>
842 : template <typename eigenT>
843 2 : void eigenCube<dataT>::sigmaMean( eigenT &mim, std::vector<dataT> &weights, dataT sigma )
844 : {
845 2 : validateWeights( weights );
846 :
847 1 : mim.resize( _rows, _cols );
848 :
849 1 : #pragma omp parallel num_threads( Eigen::nbThreads() )
850 : {
851 : std::vector<Scalar> work;
852 :
853 : #pragma omp for schedule( static, 10 )
854 : for( Index i = 0; i < _rows; ++i )
855 : {
856 : for( Index j = 0; j < _cols; ++j )
857 : {
858 : work.resize( _planes ); // work could be smaller after sigmaMean
859 : for( int k = 0; k < _planes; ++k )
860 : work[k] = ( pixel( i, j ) )( k, 0 );
861 :
862 : mim( i, j ) = math::vectorSigmaMean( work, weights, sigma );
863 : }
864 : }
865 : }
866 1 : }
867 :
868 : template <typename dataT>
869 : template <typename eigenT, typename eigenCubeT>
870 2 : void eigenCube<dataT>::sigmaMean(
871 : eigenT &mim, std::vector<dataT> &weights, eigenCubeT &mask, dataT sigma, double minGoodFract )
872 : {
873 2 : validateMaskedCombination( mask, minGoodFract );
874 1 : validateWeights( weights );
875 :
876 1 : mim.resize( _rows, _cols );
877 :
878 1 : #pragma omp parallel num_threads( Eigen::nbThreads() )
879 : {
880 : std::vector<Scalar> work, wwork;
881 :
882 : #pragma omp for schedule( static, 10 )
883 : for( Index i = 0; i < _rows; ++i )
884 : {
885 : for( Index j = 0; j < _cols; ++j )
886 : {
887 : work.clear();
888 : wwork.clear();
889 :
890 : for( Index k = 0; k < _planes; ++k )
891 : {
892 : if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
893 : {
894 : work.push_back( ( pixel( i, j ) )( k, 0 ) );
895 : wwork.push_back( weights[k] );
896 : }
897 : }
898 : if( hasEnoughGoodPixels( work.size(), minGoodFract ) )
899 : {
900 : mim( i, j ) = math::vectorSigmaMean( work, wwork, sigma );
901 : }
902 : else
903 : mim( i, j ) = invalidNumber<Scalar>();
904 : }
905 : }
906 : }
907 1 : }
908 :
909 : } // namespace improc
910 :
911 : } // namespace mx
912 :
913 : #endif // eigenCube_hpp
|
