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 : #pragma GCC system_header
14 : #include <Eigen/Dense>
15 :
16 : #include "../math/vectorUtils.hpp"
17 : #include "eigenImage.hpp"
18 : #include "imageUtils.hpp"
19 :
20 : namespace mx
21 : {
22 : namespace improc
23 : {
24 :
25 : /// An image cube with an Eigen-like API
26 : /** \ingroup eigen_image_processing
27 : */
28 : template <typename dataT>
29 : class eigenCube
30 : {
31 : public:
32 : typedef bool is_eigenCube;
33 : typedef dataT Scalar;
34 :
35 : typedef typename Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>::Index Index;
36 :
37 : typedef Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> imageRef;
38 :
39 : typedef Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic> imageT;
40 :
41 : protected:
42 : Index _rows;
43 : Index _cols;
44 : Index _planes;
45 :
46 : dataT *m_data{ nullptr };
47 :
48 : bool _owner;
49 :
50 : public:
51 : eigenCube();
52 :
53 : /// C'tor which will allocate space..
54 : /**
55 : */
56 : eigenCube( Index nrows, ///< [in] Number of rows in the cube
57 : Index ncols, ///< [in] Number of columns the cube
58 : Index nplanes ///< [in] Number of planes in the cube
59 : );
60 :
61 : /// C'tor taking an existing array as an argument.
62 : /** The existing array is used as is, as in a map, and ownership is not taken. You are
63 : * responsible for memory management, e.g. free-ing this array.
64 : */
65 : eigenCube( dataT *ndata, ///< [in] Allocated array with a cube of nrows x ncols x nplanes
66 : size_t nrows, ///< [in] Number of rows in the cube
67 : size_t ncols, ///< [in] Number of columns the cube
68 : size_t nplanes ///< [in] Number of planes in the cube
69 : );
70 :
71 : ~eigenCube();
72 :
73 : void setZero();
74 :
75 : eigenCube<dataT> &operator=( const eigenCube<dataT> &ec );
76 :
77 : void shallowCopy( eigenCube<dataT> &src, bool takeOwner = false );
78 :
79 : /// De-allocate and set all sizes to 0
80 : void clear();
81 :
82 : void resize( int r, int c, int p );
83 :
84 : void resize( int r, int c );
85 :
86 : dataT *data();
87 :
88 : const dataT *data() const;
89 :
90 : const Index rows() const;
91 :
92 : const Index cols() const;
93 :
94 : const Index planes() const;
95 :
96 : /// Returns a 2D Eigen::Eigen::Map pointed at the entire cube.
97 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> cube();
98 :
99 : /// Returns a 2D Eigen::Eigen::Map pointed at the specified image.
100 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> image( Index n /**< [in] the image number */ );
101 :
102 : /// Returns a 2D Eigen::Eigen::Map pointed at the specified image.
103 : const Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>
104 : image( Index n /**< [in] the image number */ ) const;
105 :
106 : /// Returns an Eigen::Eigen::Map-ed vector of the pixels at the given coordinate
107 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>,
108 : Eigen::Unaligned,
109 : Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>
110 : pixel( Index i, Index j );
111 :
112 : /// Return an Eigen::Eigen::Map of the cube where each image is a vector
113 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> asVectors();
114 :
115 : /// Calculate the covariance matrix of the images in the cube
116 : void Covar( Eigen::Matrix<dataT, Eigen::Dynamic, Eigen::Dynamic> &cv );
117 :
118 : /// Calculate the sum image of the cube
119 : /**
120 : * \tparam eigenT an Eigen-like type.
121 : */
122 : template <typename eigenT>
123 : void sum( eigenT &mim /**< [out] the resultant sum image. Is resized. */ );
124 :
125 : /// Calculate the mean image of the cube
126 : /**
127 : * \tparam eigenT an Eigen-like type.
128 : */
129 : template <typename eigenT>
130 : void mean( eigenT &mim /**< [out] the resultant mean image. Is resized. */ );
131 :
132 : /// Calculate the mean image of the cube with a mask.
133 : /**
134 : * \tparam eigenT an Eigen-like type.
135 : * \tparam eigenCubeT an eigenCube type.
136 : */
137 : template <typename eigenT, typename eigenCubeT>
138 : void mean( eigenT &mim, ///< [out] the resultant mean image. Is resized.
139 : eigenCubeT &mask, /**< [in] a mask cube. Only pixels with value 1 are included in the
140 : mean calculation. */
141 : double minGoodFract = 0.0 /**< [in] [optional] the minimum fraction of good pixels, if not met
142 : then the pixel is NaN-ed. */
143 : );
144 :
145 : /// Calculate the weighted mean image of the cube
146 : /**
147 : * \tparam eigenT an Eigen-like type.
148 : */
149 : template <typename eigenT>
150 : void mean( eigenT &mim, ///< [out] the resultant mean image. Is resized.
151 : std::vector<dataT> &weights ///< [in] a vector of weights to use for calculating the mean
152 : );
153 :
154 : /// Calculate the weighted mean image of the cube, with a mask cube.
155 : /**
156 : * \tparam eigenT an Eigen-like type.
157 : * \tparam eigenCubeT an eigenCube type.
158 : */
159 : template <typename eigenT, typename eigenCubeT>
160 : void mean( eigenT &mim, ///< [out] the resultant mean image. Is resized.
161 : std::vector<dataT> &weights, ///< [in] a vector of weights to use for calculating the mean
162 : eigenCubeT &mask, ///< [in] a mask cube. Only pixels with value 1 are included in the mean calculation.
163 : double minGoodFract =
164 : 0.0 ///< [in] [optional] the minimum fraction of good pixels, if not met then the pixel is NaN-ed.
165 : );
166 :
167 : /// Calculate the median image of the cube
168 : /**
169 : * \tparam eigenT an Eigen-like type.
170 : */
171 : template <typename eigenT>
172 : void median( eigenT &mim /**< [out] the resultant median image. Is resized. */ );
173 :
174 : /// Calculate the sigma clipped mean image of the cube
175 : /**
176 : * \tparam eigenT an Eigen-like type.
177 : */
178 : template <typename eigenT>
179 : void sigmaMean( eigenT &mim, ///< [out] the resultant mean image
180 : Scalar sigma ///< [in] the sigma value at which to clip.
181 : );
182 :
183 : /// Calculate the sigma clipped mean image of the cube, with a mask cube
184 : /**
185 : * \tparam eigenT an Eigen-like type.
186 : * \tparam eigenCubeT an eigenCube type.
187 : */
188 : template <typename eigenT, typename eigenCubeT>
189 : void
190 : sigmaMean( eigenT &mim, ///< [out] the resultant mean image. Is resized.
191 : eigenCubeT &mask, ///< [in] a mask cube. Only pixels with value 1 are included in the mean calculation.
192 : Scalar sigma, ///< [in] the sigma value at which to clip.
193 : double minGoodFract =
194 : 0.0 ///< [in] [optional] the minimum fraction of good pixels, if not met then the pixel is NaN-ed.
195 : );
196 :
197 : /// Calculate the sigma clipped weighted mean image of the cube
198 : /**
199 : * \tparam eigenT an Eigen-like type.
200 : */
201 : template <typename eigenT>
202 : void sigmaMean( eigenT &mim, ///< [out] the resultant mean image. Is resized.
203 : std::vector<dataT> &weights, ///< [in] a vector of weights to use for calculating the mean
204 : Scalar sigma ///< [in] the sigma value at which to clip.
205 : );
206 :
207 : /// Calculate the sigma clipped weighted mean image of the cube, with a mask cube.
208 : /**
209 : * \tparam eigenT an Eigen-like type.
210 : * \tparam eigenCubeT an eigenCube type.
211 : */
212 : template <typename eigenT, typename eigenCubeT>
213 : void
214 : sigmaMean( eigenT &mim, ///< [out] the resultant mean image. Is resized.
215 : std::vector<dataT> &weights, ///< [in] a vector of weights to use for calculating the mean
216 : eigenCubeT &mask, ///< [in] a mask cube. Only pixels with value 1 are included in the mean calculation.
217 : Scalar sigma, ///< [in] the sigma value at which to clip.
218 : double minGoodFract =
219 : 0.0 ///< [in] [optional] the minimum fraction of good pixels, if not met then the pixel is NaN-ed.
220 : );
221 : };
222 :
223 : template <typename dataT>
224 15 : eigenCube<dataT>::eigenCube()
225 : {
226 15 : _rows = 0;
227 15 : _cols = 0;
228 15 : _planes = 0;
229 15 : _owner = false;
230 15 : }
231 :
232 : template <typename dataT>
233 2 : eigenCube<dataT>::eigenCube( Index nrows, Index ncols, Index nplanes )
234 : {
235 2 : _rows = nrows;
236 2 : _cols = ncols;
237 2 : _planes = nplanes;
238 :
239 2 : m_data = new dataT[_rows * _cols * _planes];
240 2 : _owner = true;
241 2 : }
242 :
243 : template <typename dataT>
244 : eigenCube<dataT>::eigenCube( dataT *ndata, size_t nrows, size_t ncols, size_t nplanes )
245 : {
246 : _rows = nrows;
247 : _cols = ncols;
248 : _planes = nplanes;
249 :
250 : m_data = ndata;
251 : _owner = false;
252 : }
253 :
254 : template <typename dataT>
255 17 : eigenCube<dataT>::~eigenCube()
256 : {
257 17 : if( _owner && m_data )
258 : {
259 17 : delete[] m_data;
260 : }
261 17 : }
262 :
263 : template <typename dataT>
264 : void eigenCube<dataT>::setZero()
265 : {
266 : int N = _rows * _cols * _planes;
267 :
268 : for( int i = 0; i < N; ++i )
269 : m_data[i] = ( (dataT)0 );
270 : }
271 :
272 : template <typename dataT>
273 1 : eigenCube<dataT> &eigenCube<dataT>::operator=( const eigenCube<dataT> &ec )
274 : {
275 1 : resize( ec.rows(), ec.cols(), ec.planes() );
276 :
277 1 : int N = _rows * _cols * _planes;
278 :
279 4194305 : for( int i = 0; i < N; ++i )
280 4194304 : m_data[i] = ec.m_data[i];
281 :
282 1 : return *this;
283 : }
284 :
285 : template <typename dataT>
286 : void eigenCube<dataT>::shallowCopy( eigenCube<dataT> &src, bool takeOwner )
287 : {
288 : if( _owner && m_data )
289 : {
290 : delete m_data;
291 : _owner = false;
292 : }
293 :
294 : _rows = src._rows;
295 : _cols = src._cols;
296 : _planes = src._planes;
297 : m_data = src.m_data;
298 :
299 : if( takeOwner == true )
300 : {
301 : _owner = true;
302 : src._owner = false;
303 : }
304 : else
305 : {
306 : _owner = false;
307 : }
308 : }
309 :
310 : template <typename dataT>
311 : void eigenCube<dataT>::clear()
312 : {
313 : if( _owner && m_data )
314 : {
315 : delete[] m_data;
316 : }
317 :
318 : _rows = 0;
319 : _cols = 0;
320 : _planes = 0;
321 : }
322 :
323 : template <typename dataT>
324 18 : void eigenCube<dataT>::resize( int r, int c, int p )
325 : {
326 18 : if( _owner && m_data )
327 : {
328 3 : delete[] m_data;
329 : }
330 :
331 18 : _rows = r;
332 18 : _cols = c;
333 18 : _planes = p;
334 :
335 12714002 : m_data = new dataT[_rows * _cols * _planes];
336 18 : _owner = true;
337 18 : }
338 :
339 : template <typename dataT>
340 : void eigenCube<dataT>::resize( int r, int c )
341 : {
342 : resize( r, c, 1 );
343 : }
344 :
345 : template <typename dataT>
346 800 : dataT *eigenCube<dataT>::data()
347 : {
348 800 : return m_data;
349 : }
350 :
351 : template <typename dataT>
352 : const dataT *eigenCube<dataT>::data() const
353 : {
354 : return m_data;
355 : }
356 :
357 : template <typename dataT>
358 18303571 : const typename eigenCube<dataT>::Index eigenCube<dataT>::rows() const
359 : {
360 18303571 : return _rows;
361 : }
362 :
363 : template <typename dataT>
364 13976528 : const typename eigenCube<dataT>::Index eigenCube<dataT>::cols() const
365 : {
366 13976528 : return _cols;
367 : }
368 :
369 : template <typename dataT>
370 198391 : const typename eigenCube<dataT>::Index eigenCube<dataT>::planes() const
371 : {
372 198391 : return _planes;
373 : }
374 :
375 : template <typename dataT>
376 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> eigenCube<dataT>::cube()
377 : {
378 : return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>( m_data, _rows * _cols, _planes );
379 : }
380 :
381 : template <typename dataT>
382 43841925 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> eigenCube<dataT>::image( Index n )
383 : {
384 43841925 : return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>( m_data + n * _rows * _cols, _rows, _cols );
385 : }
386 :
387 : template <typename dataT>
388 4325376 : const Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> eigenCube<dataT>::image( Index n ) const
389 : {
390 4325376 : return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>( m_data + n * _rows * _cols, _rows, _cols );
391 : }
392 :
393 : template <typename dataT>
394 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>,
395 : Eigen::Unaligned,
396 : Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>
397 2 : eigenCube<dataT>::pixel( Index i, Index j )
398 : {
399 : return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>,
400 : Eigen::Unaligned,
401 : Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>(
402 2 : m_data + j * _rows + i,
403 : _planes,
404 : 1,
405 4 : Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>( 0, _rows * _cols ) );
406 : }
407 :
408 : template <typename dataT>
409 : Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>> eigenCube<dataT>::asVectors()
410 : {
411 : return Eigen::Map<Eigen::Array<dataT, Eigen::Dynamic, Eigen::Dynamic>>( m_data, _rows * _cols, _planes );
412 : }
413 :
414 : template <typename dataT>
415 : void eigenCube<dataT>::Covar( Eigen::Matrix<dataT, Eigen::Dynamic, Eigen::Dynamic> &cv )
416 : {
417 : cv = asVectors().matrix().transpose() * asVectors().matrix();
418 : }
419 :
420 : template <typename dataT>
421 : template <typename eigenT>
422 : void eigenCube<dataT>::sum( eigenT &mim )
423 : {
424 : mim.resize( _rows, _cols );
425 :
426 : #pragma omp parallel for schedule( static, 1 ) num_threads( Eigen::nbThreads() )
427 : for( Index i = 0; i < _rows; ++i )
428 : {
429 : for( Index j = 0; j < _cols; ++j )
430 : {
431 : mim( i, j ) = pixel( i, j ).sum();
432 : }
433 : }
434 : }
435 :
436 : template <typename dataT>
437 : template <typename eigenT>
438 : void eigenCube<dataT>::mean( eigenT &mim )
439 : {
440 : mim.resize( _rows, _cols );
441 :
442 : #pragma omp parallel for schedule( static, 1 ) num_threads( Eigen::nbThreads() )
443 : for( Index i = 0; i < _rows; ++i )
444 : {
445 : for( Index j = 0; j < _cols; ++j )
446 : {
447 : mim( i, j ) = pixel( i, j ).mean();
448 : }
449 : }
450 : }
451 :
452 : template <typename dataT>
453 : template <typename eigenT, typename eigenCubeT>
454 : void eigenCube<dataT>::mean( eigenT &mim, eigenCubeT &mask, double minGoodFract )
455 : {
456 : mim.resize( _rows, _cols );
457 :
458 : #pragma omp parallel
459 : {
460 : std::vector<dataT> work;
461 :
462 : #pragma omp for
463 : for( Index i = 0; i < _rows; ++i )
464 : {
465 : for( Index j = 0; j < _cols; ++j )
466 : {
467 :
468 : work.clear();
469 :
470 : for( Index k = 0; k < _planes; ++k )
471 : {
472 : if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
473 : {
474 : work.push_back( ( pixel( i, j ) )( k, 0 ) );
475 : }
476 : }
477 :
478 : if( work.size() > minGoodFract * _planes )
479 : {
480 : mim( i, j ) = math::vectorMean( work );
481 : }
482 : else
483 : {
484 : mim( i, j ) = invalidNumber<float>();
485 : }
486 : }
487 : }
488 : }
489 : }
490 :
491 : template <typename dataT>
492 : template <typename eigenT>
493 : void eigenCube<dataT>::mean( eigenT &mim, std::vector<dataT> &weights )
494 : {
495 : mim.resize( _rows, _cols );
496 :
497 : #pragma omp parallel num_threads( Eigen::nbThreads() )
498 : {
499 : std::vector<Scalar> work;
500 :
501 : #pragma omp for schedule( static, 10 )
502 : for( Index i = 0; i < _rows; ++i )
503 : {
504 : for( Index j = 0; j < _cols; ++j )
505 : {
506 : work.resize( _planes ); // work could be smaller after sigmaMean
507 : // int ii=0;
508 : for( int k = 0; k < _planes; ++k )
509 : work[k] = ( pixel( i, j ) )( k, 0 );
510 :
511 : mim( i, j ) = math::vectorMean( work, weights );
512 : }
513 : }
514 : }
515 : }
516 :
517 : template <typename dataT>
518 : template <typename eigenT, typename eigenCubeT>
519 : void eigenCube<dataT>::mean( eigenT &mim, std::vector<dataT> &weights, eigenCubeT &mask, double minGoodFract )
520 : {
521 : mim.resize( _rows, _cols );
522 :
523 : #pragma omp parallel num_threads( Eigen::nbThreads() )
524 : {
525 : std::vector<Scalar> work, wwork;
526 :
527 : #pragma omp for schedule( static, 10 )
528 : for( Index i = 0; i < _rows; ++i )
529 : {
530 : for( Index j = 0; j < _cols; ++j )
531 : {
532 : work.clear();
533 : wwork.clear();
534 :
535 : // int ii=0;
536 :
537 : for( Index k = 0; k < _planes; ++k )
538 : {
539 : if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
540 : {
541 : work.push_back( ( pixel( i, j ) )( k, 0 ) );
542 : wwork.push_back( weights[k] );
543 : }
544 : }
545 : if( work.size() > minGoodFract * _planes )
546 : {
547 : mim( i, j ) = math::vectorMean( work, wwork );
548 : }
549 : else
550 : mim( i, j ) = invalidNumber<float>();
551 : }
552 : }
553 : }
554 : }
555 :
556 : template <typename dataT>
557 : template <typename eigenT>
558 : void eigenCube<dataT>::median( eigenT &mim )
559 : {
560 : mim.resize( _rows, _cols );
561 :
562 : #pragma omp parallel for schedule( static, 10 ) num_threads( Eigen::nbThreads() )
563 : for( Index i = 0; i < _rows; ++i )
564 : {
565 : // Eigen::Eigen::Array<Scalar, Eigen::Eigen::Dynamic, Eigen::Eigen::Dynamic> work;
566 : std::vector<Scalar> work;
567 : for( Index j = 0; j < _cols; ++j )
568 : {
569 : mim( i, j ) = imageMedian( pixel( i, j ), &work );
570 : }
571 : }
572 : }
573 :
574 : template <typename dataT>
575 : template <typename eigenT>
576 : void eigenCube<dataT>::sigmaMean( eigenT &mim, dataT sigma )
577 : {
578 : mim.resize( _rows, _cols );
579 :
580 : // clang-format off
581 : #pragma omp parallel //clang-format on
582 : {
583 : std::vector<Scalar> work;
584 :
585 : // clang-format off
586 : #pragma omp for // clang-format off
587 : for( Index i = 0; i < _rows; ++i )
588 : {
589 : for( Index j = 0; j < _cols; ++j )
590 : {
591 : work.resize( _planes );
592 : for( int k = 0; k < _planes; ++k )
593 : {
594 : work[k] = ( pixel( i, j ) )( k, 0 );
595 : }
596 :
597 : mim( i, j ) = math::vectorSigmaMean( work, sigma );
598 : }
599 : }
600 : }
601 : }
602 :
603 : template <typename dataT>
604 : template <typename eigenT, typename eigenCubeT>
605 : void eigenCube<dataT>::sigmaMean( eigenT &mim, eigenCubeT &mask, dataT sigma, double minGoodFract )
606 : {
607 : mim.resize( _rows, _cols );
608 :
609 : #pragma omp parallel
610 : {
611 : std::vector<Scalar> work;
612 :
613 : #pragma omp for
614 : for( Index i = 0; i < _rows; ++i )
615 : {
616 : for( Index j = 0; j < _cols; ++j )
617 : {
618 : work.clear();
619 : int ii = 0;
620 :
621 : for( Index k = 0; k < _planes; ++k )
622 : {
623 : if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
624 : {
625 : work.push_back( ( pixel( i, j ) )( k, 0 ) );
626 : }
627 : }
628 :
629 : if( work.size() > minGoodFract * _planes )
630 : {
631 : mim( i, j ) = math::vectorSigmaMean( work, sigma );
632 : }
633 : else
634 : {
635 : mim( i, j ) = invalidNumber<float>();
636 : }
637 : }
638 : }
639 : }
640 : }
641 :
642 : template <typename dataT>
643 : template <typename eigenT>
644 : void eigenCube<dataT>::sigmaMean( eigenT &mim, std::vector<dataT> &weights, dataT sigma )
645 : {
646 : mim.resize( _rows, _cols );
647 :
648 : #pragma omp parallel num_threads( Eigen::nbThreads() )
649 : {
650 : std::vector<Scalar> work;
651 :
652 : #pragma omp for schedule( static, 10 )
653 : for( Index i = 0; i < _rows; ++i )
654 : {
655 : for( Index j = 0; j < _cols; ++j )
656 : {
657 : work.resize( _planes ); // work could be smaller after sigmaMean
658 : int ii = 0;
659 : for( int k = 0; k < _planes; ++k )
660 : work[k] = ( pixel( i, j ) )( k, 0 );
661 :
662 : mim( i, j ) = math::vectorSigmaMean( work, weights, sigma );
663 : }
664 : }
665 : }
666 : }
667 :
668 : template <typename dataT>
669 : template <typename eigenT, typename eigenCubeT>
670 : void eigenCube<dataT>::sigmaMean(
671 : eigenT &mim, std::vector<dataT> &weights, eigenCubeT &mask, dataT sigma, double minGoodFract )
672 : {
673 : mim.resize( _rows, _cols );
674 :
675 : #pragma omp parallel num_threads( Eigen::nbThreads() )
676 : {
677 : std::vector<Scalar> work, wwork;
678 :
679 : #pragma omp for schedule( static, 10 )
680 : for( Index i = 0; i < _rows; ++i )
681 : {
682 : for( Index j = 0; j < _cols; ++j )
683 : {
684 : work.clear();
685 : wwork.clear();
686 :
687 : int ii = 0;
688 :
689 : for( Index k = 0; k < _planes; ++k )
690 : {
691 : if( ( mask.pixel( i, j ) )( k, 0 ) == 1 )
692 : {
693 : work.push_back( ( pixel( i, j ) )( k, 0 ) );
694 : wwork.push_back( weights[k] );
695 : }
696 : }
697 : if( work.size() > minGoodFract * _planes )
698 : {
699 : mim( i, j ) = math::vectorSigmaMean( work, wwork, sigma );
700 : }
701 : else
702 : mim( i, j ) = invalidNumber<float>();
703 : }
704 : }
705 : }
706 : }
707 :
708 : } // namespace improc
709 :
710 : } // namespace mx
711 :
712 : #endif // eigenCube_hpp
|
