Line data Source code
1 : /** \file imageFilters.hpp
2 : * \brief Image filters (smoothing, radial profiles, etc.)
3 : * \ingroup image_processing_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : #ifndef __imageFilters_hpp__
9 : #define __imageFilters_hpp__
10 :
11 : #include <algorithm>
12 : #include <cmath>
13 : #include <cstdlib>
14 : #include <format>
15 : #include <limits>
16 :
17 : #include "../math/floatUtils.hpp"
18 : #include "../math/gslInterpolator.hpp"
19 : #include "../math/vectorUtils.hpp"
20 : #include "../math/geo.hpp"
21 :
22 : #include "imageMasks.hpp"
23 :
24 : namespace mx
25 : {
26 : namespace improc
27 : {
28 :
29 : /** \addtogroup image_filters_kernels
30 : * @{
31 : * The filter function use a kernel that specifies how to filter the image. Filter kernels, usually
32 : * denoted as type kernelT below, must have the following interface:
33 : * \code
34 : * struct filterKernel
35 : * {
36 : * typedef <eigen-like type> arrayT; // arrayT must have an eigen-like interface
37 : *
38 : * typedef <type> arithT; // the type used for arithmetic, normally `typename arrayT::Scalar`
39 : *
40 : * typedef <verbosity-type> verboseT; // the mxlib verbosity for error reports
41 : *
42 : * //The maxWidth function returns the maximum possible full-width (in either direction) of the kernel
43 : * // Called only once for each call to the filter function
44 : * int maxWidth() const
45 : * {
46 : * //returns the maximum half-width given the configuration
47 : * }
48 : *
49 : * //The setKernel function is called for each pixel.
50 : * void setKernel( arithT x, // the pixel x-position relative to the image center (not the pixel index)
51 : * arithT y, // the pixel y-position relative to the image center (not the pixel index)
52 : * arrayT & kernel // the array to be resized and populated
53 : * ) const
54 : * {
55 : * //This must resize and populate the passed in kernel array each time
56 : * //so that it is re-entrant.
57 : *
58 : * //Note: On output kernel array should be normalized so that sum() = 1.0
59 : *
60 : * //Note: the width and height of the kernel array should always be odd
61 : * }
62 : * };
63 : * \endcode
64 : * Additionally `kernelT` must be copyable.
65 : *
66 : */
67 : /// Symetric Gaussian smoothing kernel
68 : /** \ingroup image_filters_kernels
69 : *
70 : */
71 : template <typename _arrayT, size_t _kernW = 4, class _verboseT = verbose::d>
72 : struct gaussKernel
73 : {
74 : typedef _arrayT arrayT;
75 : typedef typename _arrayT::Scalar arithT;
76 : static const int kernW = _kernW;
77 :
78 : typedef _verboseT verboseT;
79 :
80 : arrayT kernel;
81 :
82 : arithT _fwhm;
83 :
84 2 : explicit gaussKernel( arithT fwhm )
85 2 : {
86 2 : _fwhm = fwhm;
87 :
88 2 : int w = kernW * _fwhm;
89 :
90 2 : if( w % 2 == 0 )
91 2 : w++;
92 :
93 2 : kernel.resize( w, w );
94 :
95 2 : arithT kcen = 0.5 * ( w - 1.0 );
96 :
97 2 : arithT sig2 = _fwhm / 2.354820045030949327;
98 2 : sig2 *= sig2;
99 :
100 : arithT r2;
101 32 : for( int i = 0; i < w; ++i )
102 : {
103 552 : for( int j = 0; j < w; ++j )
104 : {
105 522 : r2 = pow( i - kcen, 2 ) + pow( j - kcen, 2 );
106 522 : kernel( i, j ) = exp( -r2 / ( 2.0 * sig2 ) );
107 : }
108 : }
109 :
110 2 : kernel /= kernel.sum();
111 2 : }
112 :
113 2 : int maxWidth() const
114 : {
115 2 : return _kernW * _fwhm;
116 : }
117 :
118 307328 : error_t setKernel( arithT x, /**< [in] x-coordinate relative to image center */
119 : arithT y, /**< [in] x-coordinate relative to image center */
120 : arrayT &kernelArray /**< [in] the array to populate with the kernel, resized */
121 : ) const
122 : {
123 : // Unused parts of interface:
124 : static_cast<void>( x );
125 : static_cast<void>( y );
126 :
127 307328 : kernelArray = kernel;
128 :
129 307328 : return error_t::noerror;
130 : }
131 : };
132 :
133 : /// Azimuthally variable boxcar kernel.
134 : /** Averages the image in a boxcar defined by a radial and azimuthal extent.
135 : *
136 : * \ingroup image_filters_kernels
137 : */
138 : template <typename _arrayT, size_t _kernW = 2, class _verboseT = verbose::d>
139 : struct azBoxKernel
140 : {
141 : typedef _arrayT arrayT;
142 : typedef typename _arrayT::Scalar arithT;
143 :
144 : inline static constexpr int kernW = static_cast<int>( _kernW ); ///< kernel sampling factor.
145 :
146 : typedef _verboseT verboseT;
147 :
148 : arithT m_radWidth{ 0 }; ///< the half-width of the averaging box, in the radial direction, in pixels.
149 : arithT m_azWidth{ 0 }; ///< the half-width of the averaging box, in the azimuthal direction, in pixels.
150 : arithT m_maxAz{ 0 }; ///< maximum azimuthal half-width in radians; 0 means no angular limit.
151 :
152 : int m_maxWidth{ 0 }; ///< maximum kernel half-width needed to keep every generated kernel in bounds.
153 :
154 : /// Construct a kernel without an angular-position limit.
155 20 : azBoxKernel( arithT radWidth, ///< [in] the half-width of the averaging box, in the radial direction, in pixels.
156 : arithT azWidth ///< [in] the half-width of the averaging box, in the azimuthal direction, in pixels.
157 : )
158 20 : : m_radWidth( fabs( radWidth ) ), m_azWidth( fabs( azWidth ) )
159 : {
160 20 : setMaxWidth();
161 20 : }
162 :
163 : /// Construct a kernel with an optional angular-position limit.
164 6 : azBoxKernel( arithT radWidth, ///< [in] the half-width of the averaging box, in the radial direction, in pixels.
165 : arithT azWidth, ///< [in] the half-width of the averaging box, in the azimuthal direction, in pixels.
166 : arithT maxAz ///< [in] the maximum half-width of the averaging box in the azimuthal direction, in
167 : ///< degrees. \>= 0. If 0 or \>= 180, then no maximum is enforced.
168 : )
169 6 : : m_radWidth( fabs( radWidth ) ), m_azWidth( fabs( azWidth ) )
170 : {
171 6 : setMaxWidth();
172 :
173 6 : if( !math::isFinite( maxAz ) )
174 : {
175 1 : m_maxAz = maxAz;
176 1 : return;
177 : }
178 :
179 5 : maxAz = fabs( maxAz );
180 5 : if( maxAz >= 180 )
181 1 : maxAz = 0; // Larger than 180 means no limit.
182 :
183 5 : m_maxAz = math::dtor( maxAz );
184 : }
185 :
186 : /// Sets the max width based on the configured az and rad widths.
187 26 : void setMaxWidth()
188 : {
189 26 : m_maxWidth = 0;
190 :
191 49 : if( kernW <= 0 || !math::isFinite( m_radWidth ) || !math::isFinite( m_azWidth ) ||
192 23 : ( m_radWidth == 0 && m_azWidth == 0 ) )
193 : {
194 4 : return;
195 : }
196 :
197 22 : const arithT maximumDimension =
198 22 : kernW * ( std::floor( std::hypot( m_radWidth, m_azWidth ) ) + static_cast<arithT>( 1 ) );
199 :
200 22 : if( !math::isFinite( maximumDimension ) || maximumDimension > std::numeric_limits<int>::max() )
201 : {
202 1 : return;
203 : }
204 :
205 21 : m_maxWidth = static_cast<int>( maximumDimension / static_cast<arithT>( 2 ) );
206 : }
207 :
208 : /// Get the maximum kernel half-width in either image dimension.
209 12 : int maxWidth() const
210 : {
211 12 : return m_maxWidth;
212 : }
213 :
214 : /// Generate a normalized kernel at the requested image-relative coordinate.
215 312 : error_t setKernel( arithT x, /**< [in] x-coordinate relative to image center */
216 : arithT y, /**< [in] x-coordinate relative to image center */
217 : arrayT &kernel /**< [in] the array to populate with the kernel, resized */
218 : ) const
219 : {
220 312 : kernel.resize( 0, 0 );
221 :
222 620 : if( kernW <= 0 || !math::isFinite( m_radWidth ) || !math::isFinite( m_azWidth ) || !math::isFinite( m_maxAz ) ||
223 308 : m_maxWidth < 0 )
224 : {
225 8 : return internal::mxlib_error_report<verboseT>( error_t::invalidconfig,
226 12 : "kernel widths and maxAz must be finite" );
227 : }
228 :
229 308 : if( !math::isFinite( x ) || !math::isFinite( y ) )
230 : {
231 6 : return internal::mxlib_error_report<verboseT>( error_t::invalidarg, "kernel coordinates must be finite" );
232 : }
233 :
234 306 : if( m_radWidth == 0 && m_azWidth == 0 )
235 : {
236 2 : kernel.resize( 1, 1 );
237 2 : kernel( 0, 0 ) = 1;
238 2 : return error_t::noerror;
239 : }
240 :
241 304 : const arithT rad0 = std::hypot( x, y );
242 :
243 304 : arithT sinq = 0;
244 304 : arithT cosq = 1;
245 304 : if( rad0 > 0 )
246 : {
247 293 : sinq = y / rad0;
248 293 : cosq = x / rad0;
249 : }
250 :
251 : // Only calc q if we're going to use it.
252 304 : arithT q = 0;
253 304 : if( m_maxAz > 0 && rad0 > 0 )
254 : {
255 2 : q = std::atan2( y, x );
256 : }
257 :
258 304 : const arithT width =
259 304 : kernW * ( std::floor( std::fabs( m_azWidth * sinq ) + std::fabs( m_radWidth * cosq ) ) + 1 );
260 304 : const arithT height =
261 304 : kernW * ( std::floor( std::fabs( m_azWidth * cosq ) + std::fabs( m_radWidth * sinq ) ) + 1 );
262 :
263 607 : if( !math::isFinite( width ) || !math::isFinite( height ) || width > std::numeric_limits<int>::max() ||
264 303 : height > std::numeric_limits<int>::max() )
265 : {
266 2 : return internal::mxlib_error_report<verboseT>( error_t::invalidconfig,
267 3 : "kernel dimensions exceed the supported size" );
268 : }
269 :
270 303 : int w = static_cast<int>( width );
271 303 : int h = static_cast<int>( height );
272 :
273 : // Retain the output pixel itself so any valid angular limit has non-empty support. Increasing an even
274 : // dimension by one does not change its integer half-width.
275 303 : if( m_maxAz > 0 )
276 : {
277 5 : if( w % 2 == 0 )
278 : {
279 5 : ++w;
280 : }
281 5 : if( h % 2 == 0 )
282 : {
283 5 : ++h;
284 : }
285 : }
286 :
287 303 : if( w / 2 > m_maxWidth )
288 : {
289 1 : return internal::mxlib_error_report<verboseT>( error_t::sizeerr,
290 1 : std::format( "Width half-width bigger than maxWidth. "
291 : "This is a bug. Details: "
292 : "|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|",
293 : kernW,
294 1 : m_radWidth,
295 1 : m_azWidth,
296 1 : m_maxWidth,
297 : x,
298 : y,
299 : rad0,
300 : sinq,
301 : cosq,
302 : w,
303 1 : h ) );
304 : }
305 :
306 302 : if( h / 2 > m_maxWidth )
307 : {
308 1 : return internal::mxlib_error_report<verboseT>( error_t::sizeerr,
309 1 : std::format( "Height half-width bigger than maxWidth. "
310 : "This is a bug. Details: "
311 : "|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|",
312 : kernW,
313 1 : m_radWidth,
314 1 : m_azWidth,
315 1 : m_maxWidth,
316 : x,
317 : y,
318 : rad0,
319 : sinq,
320 : cosq,
321 : w,
322 1 : h ) );
323 : }
324 :
325 301 : kernel.resize( w, h );
326 :
327 301 : const arithT xcen = 0.5 * ( w - 1.0 );
328 301 : const arithT ycen = 0.5 * ( h - 1.0 );
329 :
330 1614 : for( int j = 0; j < h; ++j )
331 : {
332 7588 : for( int i = 0; i < w; ++i )
333 : {
334 6275 : const arithT dx = i - xcen;
335 6275 : const arithT dy = j - ycen;
336 6275 : const arithT sampleX = x + dx;
337 6275 : const arithT sampleY = y + dy;
338 6275 : const arithT radP = std::hypot( sampleX, sampleY );
339 :
340 6275 : if( std::fabs( radP - rad0 ) > m_radWidth )
341 : {
342 622 : kernel( i, j ) = 0;
343 622 : continue;
344 : }
345 :
346 5653 : const arithT tangentialOffset = dy * cosq - dx * sinq;
347 5653 : if( std::fabs( tangentialOffset ) <= m_azWidth )
348 : {
349 5604 : if( m_maxAz > 0 ) // Only check this if needed.
350 : {
351 104 : arithT q2 = 0;
352 104 : if( radP > 0 )
353 : {
354 100 : q2 = std::atan2( sampleY, sampleX );
355 : }
356 :
357 104 : const arithT dq = math::angleDiff<math::radiansT<arithT>>( q, q2 );
358 104 : if( std::fabs( dq ) > m_maxAz )
359 : {
360 80 : kernel( i, j ) = 0;
361 80 : continue;
362 : }
363 : }
364 5524 : kernel( i, j ) = 1;
365 : }
366 : else
367 49 : kernel( i, j ) = 0;
368 : }
369 : }
370 :
371 301 : const arithT ksum = kernel.sum();
372 301 : if( !math::isFinite( ksum ) || ksum <= 0 )
373 : {
374 1 : kernel.resize( 0, 0 );
375 1 : return internal::mxlib_error_report<verboseT>( error_t::invalidconfig,
376 2 : std::format( "kernel sum 0 at {},{}", x, y ) );
377 : }
378 300 : kernel /= ksum;
379 :
380 300 : return error_t::noerror;
381 : }
382 : };
383 :
384 : /// A kernel that is pre-calculated for the entire image, useful for repeated applications
385 : /** Use this for spatially variable kernels. It is not needed for, e.g., the symmetric Gaussian.
386 : *
387 : * \tparam kernelT is the kernel type. See above for the requirements on kernelT.
388 : */
389 : template <class kernelT>
390 : struct precalcKernel
391 : {
392 : typedef kernelT::arrayT arrayT;
393 : typedef kernelT::arrayT::Scalar arithT;
394 : typedef kernelT::verboseT verboseT;
395 :
396 : kernelT m_kernel; ///< copied production kernel used to populate the cache.
397 :
398 : uint32_t m_rows; ///< number of image rows represented by the cache.
399 : uint32_t m_cols; ///< number of image columns represented by the cache.
400 : arithT m_xcen; ///< pixel x-coordinate of the image center.
401 : arithT m_ycen; ///< pixel y-coordinate of the image center.
402 :
403 : int m_maxWidth; ///< maximum half-width reported by the copied production kernel.
404 :
405 : std::vector<arrayT> m_kernels; ///< generated kernels in column-major image-coordinate order.
406 :
407 : /// Disallow construction without a production kernel and image geometry.
408 : precalcKernel() = delete;
409 :
410 : /// Pre-calculate a production kernel for every coordinate in an image.
411 : /** \throws mx::exception when the production kernel cannot generate any cache entry. */
412 8 : precalcKernel( const kernelT &kernel, /**< [in] A fully initialized kernel. Is copied.*/
413 : uint32_t rows, /**< [in] The rows in the images to be filtered*/
414 : uint32_t cols, /**< [in] The columns in the images to be filtered*/
415 : arithT xcen, /**< [in] The pixel x-coordinate of the center of the image*/
416 : arithT ycen /**< [in] The pixel y-coordinate of the center of the image*/
417 : )
418 8 : : m_kernel( kernel ), m_rows( rows ), m_cols( cols ), m_xcen( xcen ), m_ycen( ycen )
419 : {
420 8 : m_maxWidth = kernel.maxWidth();
421 :
422 8 : m_kernels.resize( static_cast<size_t>( m_rows ) * static_cast<size_t>( m_cols ) );
423 8 : size_t n = 0;
424 51 : for( uint32_t cc = 0; cc < m_cols; ++cc )
425 : {
426 327 : for( uint32_t rr = 0; rr < m_rows; ++rr )
427 : {
428 284 : const error_t result = m_kernel.setKernel( rr - xcen, cc - ycen, m_kernels[n] );
429 284 : if( result != error_t::noerror )
430 : {
431 2 : throw mx::exception<verboseT>(
432 : result,
433 2 : std::format( "failed to pre-calculate kernel at row {}, column {}", rr, cc ) );
434 : }
435 283 : ++n;
436 : }
437 : }
438 8 : }
439 :
440 : /// Get the maximum half-width reported by the cached production kernel.
441 1 : int maxWidth() const
442 : {
443 1 : return m_maxWidth;
444 : }
445 :
446 : /// Retrieve the cached kernel at an integral image-relative coordinate.
447 220 : error_t setKernel( arithT x, /**< [in] x-coordinate relative to image center */
448 : arithT y, /**< [in] x-coordinate relative to image center */
449 : arrayT &kernel /**< [in] the array to populate with the kernel, resized */
450 : ) const
451 : {
452 220 : if( m_kernels.size() == 0 )
453 : {
454 1 : return error_t::invalidconfig;
455 : }
456 :
457 219 : const arithT row = x + m_xcen;
458 219 : const arithT column = y + m_ycen;
459 :
460 437 : if( !math::isFinite( row ) || !math::isFinite( column ) || row < 0 || column < 0 || row >= m_rows ||
461 437 : column >= m_cols || std::floor( row ) != row || std::floor( column ) != column )
462 : {
463 6 : return error_t::invalidarg;
464 : }
465 :
466 213 : const size_t n = static_cast<size_t>( column ) * m_rows + static_cast<size_t>( row );
467 213 : kernel = m_kernels[n];
468 :
469 213 : return error_t::noerror;
470 : }
471 : };
472 :
473 : /// Filter an image with a mean kernel.
474 : /** Applies the kernel to each pixel in the image and sums, storing the filtered result in the output image.
475 : *
476 : * \tparam imageOutT the type of the output image (must have an Eigen-like interface)
477 : * \tparam imageInT the type of the input image (must have an Eigen-like interface)
478 : * \tparam kernelT is the kernel type (see above for requirements)
479 : *
480 : * \ingroup image_filters_kernels
481 : *
482 : */
483 : template <typename imageOutT, typename imageInT, typename kernelT>
484 1 : error_t filterImage( imageOutT &fim, /**< [out] Contains the filtered image, will be resized*/
485 : imageInT im, /**< [in] the image to be filtered*/
486 : const kernelT &kernel, /**< [in] a fully configured kernel object*/
487 : int maxr = 0 /**< [in] [opt] the maximum radius from the image center to
488 : apply the kernel. pixels outside this radius are
489 : set to 0.*/
490 : )
491 : {
492 : typedef typename kernelT::verboseT verboseT;
493 :
494 1 : fim.resize( im.rows(), im.cols() );
495 :
496 1 : float xcen = 0.5 * ( im.rows() - 1 );
497 1 : float ycen = 0.5 * ( im.cols() - 1 );
498 :
499 1 : if( maxr == 0 )
500 : {
501 1 : maxr = 0.5 * std::min( im.rows(), im.cols() ) - kernel.maxWidth();
502 : }
503 :
504 1 : int mini = 0.5 * im.rows() - maxr;
505 1 : int maxi = 0.5 * im.rows() + maxr;
506 1 : int minj = 0.5 * im.cols() - maxr;
507 1 : int maxj = 0.5 * im.cols() + maxr;
508 :
509 1 : typename kernelT::arrayT kernelArray;
510 :
511 1 : error_t top_errc = error_t::noerror;
512 :
513 : // clang-format off
514 1 : #pragma omp parallel private( kernelArray ) // clang-format on
515 : {
516 : int im_i, im_j, im_p, im_q;
517 : int kern_i, kern_j, kern_p, kern_q;
518 : typename imageOutT::Scalar norm;
519 :
520 : error_t errc = error_t::noerror;
521 :
522 : // clang-format off
523 : #pragma omp for // clang-format on
524 : for( int i = 0; i < im.rows(); ++i )
525 : {
526 : if( top_errc != error_t::noerror ) // can't return from omp loop
527 : {
528 : continue;
529 : }
530 :
531 : for( int j = 0; j < im.cols(); ++j )
532 : {
533 : if( errc != error_t::noerror ) // can't return from omp loop
534 : {
535 : continue;
536 : }
537 :
538 : if( ( i >= mini && i < maxi ) && ( j >= minj && j < maxj ) )
539 : {
540 : errc = kernel.setKernel( i - xcen, j - ycen, kernelArray );
541 :
542 : fim( i, j ) = ( im.block( i - 0.5 * ( kernelArray.rows() - 1 ),
543 : j - 0.5 * ( kernelArray.cols() - 1 ),
544 : kernelArray.rows(),
545 : kernelArray.cols() ) *
546 : kernelArray )
547 : .sum();
548 : }
549 : else
550 : {
551 : errc = kernel.setKernel( i - xcen, j - ycen, kernelArray );
552 :
553 : im_i = i - 0.5 * ( kernelArray.rows() - 1 );
554 : if( im_i < 0 )
555 : im_i = 0;
556 :
557 : im_j = j - 0.5 * ( kernelArray.cols() - 1 );
558 : if( im_j < 0 )
559 : im_j = 0;
560 :
561 : im_p = im.rows() - im_i;
562 : if( im_p > kernelArray.rows() )
563 : im_p = kernelArray.rows();
564 :
565 : im_q = im.cols() - im_j;
566 : if( im_q > kernelArray.cols() )
567 : im_q = kernelArray.cols();
568 :
569 : kern_i = 0.5 * ( kernelArray.rows() - 1 ) - i;
570 : if( kern_i < 0 )
571 : kern_i = 0;
572 :
573 : kern_j = 0.5 * ( kernelArray.cols() - 1 ) - j;
574 : if( kern_j < 0 )
575 : kern_j = 0;
576 :
577 : kern_p = kernelArray.rows() - kern_i;
578 : if( kern_p > kernelArray.rows() )
579 : kern_p = kernelArray.rows();
580 :
581 : kern_q = kernelArray.cols() - kern_j;
582 : if( kern_q > kernelArray.cols() )
583 : kern_q = kernelArray.cols();
584 :
585 : // Pick only the smallest widths
586 : if( im_p < kern_p )
587 : kern_p = im_p;
588 : if( im_q < kern_q )
589 : kern_q = im_q;
590 :
591 : norm = kernelArray.block( kern_i, kern_j, kern_p, kern_q ).sum();
592 :
593 : fim( i, j ) =
594 : ( im.block( im_i, im_j, kern_p, kern_q ) * kernelArray.block( kern_i, kern_j, kern_p, kern_q ) )
595 : .sum() /
596 : norm;
597 : if( !std::isfinite( fim( i, j ) ) )
598 : fim( i, j ) = 0.0;
599 : }
600 : } // for rows
601 :
602 : // clang-format off
603 : #pragma omp critical // clang-format on
604 : {
605 : if( errc != error_t::noerror && top_errc == error_t::noerror )
606 : {
607 : top_errc = errc;
608 : }
609 : }
610 :
611 : } // for cols
612 : } // pragma omp parallel
613 :
614 1 : return top_errc;
615 1 : }
616 :
617 : /// Filter an image with a median kernel.
618 : /** Calculates the median of all pixels corresponding to non-zero pixels in the kernel,
619 : * storing the filtered result in the output image.
620 : *
621 : * \tparam imageOutT the type of the output image (must have an Eigen-like interface)
622 : * \tparam imageInT the type of the input image (must have an Eigen-like interface)
623 : * \tparam kernelT is the kernel type (see above for requirements)
624 : *
625 : * \ingroup image_filters_kernels
626 : *
627 : */
628 : template <typename imageOutT, typename imageInT, typename kernelT>
629 1 : void medianFilterImage( imageOutT &fim, /**< [out] Contains the filtered image, will be resized*/
630 : imageInT im, /**< [in] the image to be filtered*/
631 : const kernelT &kernel, /**< [in] a fully configured kernel object*/
632 : int maxr = 0, /**< [in] [opt] the maximum radius from the image center to apply
633 : the kernel. Psixels outside this radius are set to 0.*/
634 : int maxrproc = 1 )
635 : {
636 1 : fim.resize( im.rows(), im.cols() );
637 :
638 1 : float xcen = 0.5 * ( im.rows() - 1 );
639 1 : float ycen = 0.5 * ( im.cols() - 1 );
640 :
641 1 : if( maxr == 0 )
642 1 : maxr = 0.5 * std::min( im.rows(), im.cols() ) - kernel.maxWidth();
643 :
644 1 : int mini = 0.5 * im.rows() - maxr;
645 1 : int maxi = 0.5 * im.rows() + maxr;
646 1 : int minj = 0.5 * im.cols() - maxr;
647 1 : int maxj = 0.5 * im.cols() + maxr;
648 :
649 1 : typename kernelT::arrayT kernelArray;
650 :
651 : // clang-format off
652 1 : #pragma omp parallel private( kernelArray ) // clang-format on
653 : {
654 : int im_i, im_j, im_p, im_q;
655 : int kern_i, kern_j, kern_p, kern_q;
656 : typename imageOutT::Scalar norm;
657 :
658 : std::vector<typename imageOutT::Scalar> pixels;
659 :
660 : // clang-format off
661 : #pragma omp for // clang-format on
662 : for( int i = 0; i < im.rows(); ++i )
663 : {
664 : for( int j = 0; j < im.cols(); ++j )
665 : {
666 : if( ( i >= mini && i < maxi ) && ( j >= minj && j < maxj ) )
667 : {
668 : kernel.setKernel( i - xcen, j - ycen, kernelArray );
669 :
670 : pixels.clear();
671 : for( int cc = 0; cc < kernelArray.cols(); ++cc )
672 : {
673 : for( int rr = 0; rr < kernelArray.rows(); ++rr )
674 : {
675 : if( kernelArray( rr, cc ) != 0 )
676 : pixels.push_back( im.block( i - 0.5 * ( kernelArray.rows() - 1 ),
677 : j - 0.5 * ( kernelArray.cols() - 1 ),
678 : kernelArray.rows(),
679 : kernelArray.cols() )( rr, cc ) );
680 : }
681 : }
682 :
683 : fim( i, j ) = math::vectorMedianInPlace( pixels );
684 : }
685 : else
686 : {
687 : if( maxrproc == 2 )
688 : {
689 : fim( i, j ) = 0;
690 : continue;
691 : }
692 : kernel.setKernel( i - xcen, j - ycen, kernelArray );
693 :
694 : im_i = i - 0.5 * ( kernelArray.rows() - 1 );
695 : if( im_i < 0 )
696 : im_i = 0;
697 :
698 : im_j = j - 0.5 * ( kernelArray.cols() - 1 );
699 : if( im_j < 0 )
700 : im_j = 0;
701 :
702 : im_p = im.rows() - im_i;
703 : if( im_p > kernelArray.rows() )
704 : im_p = kernelArray.rows();
705 :
706 : im_q = im.cols() - im_j;
707 : if( im_q > kernelArray.cols() )
708 : im_q = kernelArray.cols();
709 :
710 : kern_i = 0.5 * ( kernelArray.rows() - 1 ) - i;
711 : if( kern_i < 0 )
712 : kern_i = 0;
713 :
714 : kern_j = 0.5 * ( kernelArray.cols() - 1 ) - j;
715 : if( kern_j < 0 )
716 : kern_j = 0;
717 :
718 : kern_p = kernelArray.rows() - kern_i;
719 : if( kern_p > kernelArray.rows() )
720 : kern_p = kernelArray.rows();
721 :
722 : kern_q = kernelArray.cols() - kern_j;
723 : if( kern_q > kernelArray.cols() )
724 : kern_q = kernelArray.cols();
725 :
726 : // Pick only the smallest widths
727 : if( im_p < kern_p )
728 : kern_p = im_p;
729 : if( im_q < kern_q )
730 : kern_q = im_q;
731 :
732 : norm = kernelArray.block( kern_i, kern_j, kern_p, kern_q ).sum();
733 :
734 : pixels.clear();
735 : for( int cc = 0; cc < kern_q; ++cc )
736 : {
737 : for( int rr = 0; rr < kern_p; ++rr )
738 : {
739 : if( kernelArray.block( kern_i, kern_j, kern_p, kern_q )( rr, cc ) != 0 )
740 : {
741 : pixels.push_back( im.block( im_i, im_j, kern_p, kern_q )( rr, cc ) );
742 : }
743 : }
744 : }
745 :
746 : fim( i, j ) = math::vectorMedianInPlace( pixels );
747 :
748 : // fim(i,j) = ( im.block(im_i, im_j, kern_p, kern_q) * kernelArray.block(kern_i, kern_j, kern_p,
749 : // kern_q )).sum()/norm;
750 : }
751 : }
752 : }
753 : } // pragma omp parallel
754 1 : }
755 :
756 : ///@}
757 :
758 : /// Smooth an image using the mean in a rectangular box, optionally rejecting the highest and lowest values.
759 : /** Calculates the mean value in a rectangular box of imIn, of size meanFullWidth X meanFullWidth and stores it in the
760 : * corresponding center pixel of imOut. For even widths, the window is associated with the higher-index member of the
761 : * central pair: it contains meanFullWidth/2 pixels before the output pixel and one fewer after it. Pixels without a
762 : * complete window are not changed in imOut (i.e. you should initialize them before the call).
763 : *
764 : * imOut is not re-allocated.
765 : *
766 : * If rejectMinMax is `true` then the minimum and maximum value in the box are not included in the mean. Rejection
767 : * makes the algorithm somewhat slower, depending on box width.
768 : *
769 : * \tparam imageTout is an eigen-like image array
770 : * \tparam imageTin is an eigen-like image array
771 : *
772 : * \returns 0 on success
773 : * \returns -1 on error.
774 : *
775 : * \ingroup image_filters_average
776 : */
777 : template <typename imageTout, typename imageTin>
778 3 : int meanSmooth( imageTout &imOut, /**< [out] the smoothed image. Not re-allocated, and the edge
779 : pixels are not modified */
780 : const imageTin &imIn, /**< [in] the image to smooth */
781 : int meanFullWidth, /**< [in] the full width of the smoothing box */
782 : bool rejectMinMax = false /**< [in] whether to reject the minimum and maximum values */
783 : )
784 : {
785 : typedef typename imageTout::Scalar scalarT;
786 :
787 6 : if( meanFullWidth <= 0 || meanFullWidth > imIn.rows() || meanFullWidth > imIn.cols() ||
788 9 : imOut.rows() != imIn.rows() || imOut.cols() != imIn.cols() )
789 : {
790 0 : return -1;
791 : }
792 :
793 3 : const int before = meanFullWidth / 2;
794 3 : const int after = meanFullWidth - before - 1;
795 3 : int nPix = meanFullWidth * meanFullWidth;
796 :
797 3 : if( rejectMinMax && nPix <= 2 )
798 : {
799 0 : return -1;
800 : }
801 :
802 3 : if( rejectMinMax ) // avoid the branch on every pixel
803 : {
804 1 : nPix -= 2;
805 2 : for( int jj = before; jj < imIn.cols() - after; ++jj )
806 : {
807 2 : for( int ii = before; ii < imIn.rows() - after; ++ii )
808 : {
809 1 : scalarT sum = 0;
810 1 : scalarT max = std::numeric_limits<scalarT>::lowest();
811 1 : scalarT min = std::numeric_limits<scalarT>::max();
812 3 : for( int ll = 0; ll < meanFullWidth; ++ll )
813 : {
814 6 : for( int kk = 0; kk < meanFullWidth; ++kk )
815 : {
816 4 : const scalarT value = imIn( ii - before + kk, jj - before + ll );
817 4 : sum += value;
818 4 : if( value > max )
819 3 : max = value;
820 4 : if( value < min )
821 1 : min = value;
822 : }
823 : }
824 1 : imOut( ii, jj ) = ( sum - max - min ) / nPix;
825 : }
826 : }
827 : }
828 : else
829 : {
830 4 : for( int jj = before; jj < imIn.cols() - after; ++jj )
831 : {
832 4 : for( int ii = before; ii < imIn.rows() - after; ++ii )
833 : {
834 2 : scalarT sum = 0;
835 10 : for( int ll = 0; ll < meanFullWidth; ++ll )
836 : {
837 40 : for( int kk = 0; kk < meanFullWidth; ++kk )
838 : {
839 32 : sum += imIn( ii - before + kk, jj - before + ll );
840 : }
841 : }
842 2 : imOut( ii, jj ) = sum / nPix;
843 : }
844 : }
845 : }
846 :
847 3 : return 0;
848 : }
849 :
850 : /** \brief Smooth an image using the mean in a rectangular box, optionally rejecting the highest and lowest values.
851 : * Determines the location and value of the highest pixel.
852 : *
853 : * Calculates the mean value in a rectangular box of imIn, of size meanFullWidth X meanFullWidth and stores it in the
854 : * corresponding center pixel of imOut. For even widths, the window is associated with the higher-index member of the
855 : * central pair: it contains meanFullWidth/2 pixels before the output pixel and one fewer after it. Pixels without a
856 : * complete window are not changed in imOut (i.e. you should initialize them before the call).
857 : *
858 : * imOut is not re-allocated.
859 : *
860 : * If rejectMinMax is `true` then the minimum and maximum value in the box are not included in the mean. Rejection
861 : * makes the somewhat slower, depending on box width.
862 : *
863 : * This version also determines the location and value of the maximum pixel. This adds some overhead, maybe on the
864 : * order of 10% slower than without.
865 : *
866 : * \overload
867 : *
868 : * \tparam imageTout is an eigen-like image array
869 : * \tparam imageTin is an eigen-like image array
870 : *
871 : * \returns 0 on success
872 : * \returns -1 on error.
873 : *
874 : * \ingroup image_filters_average
875 : */
876 : template <typename imageTout, typename imageTin>
877 3 : int meanSmooth(
878 : imageTout &imOut, /**< [out] the smoothed image. Not re-allocated, and the edge pixels are not modified */
879 : int &xMax, /**< [out] the x location of the maximum pixel */
880 : int &yMax, /**< [out] the y location of the maximum pixel */
881 : typename imageTout::Scalar &pMax, /**< [out] the value of the maximum pixel */
882 : const imageTin &imIn, /**< [in] the image to smooth */
883 : int meanFullWidth, /**< [in] the full width of the smoothing box */
884 : bool rejectMinMax = false /**< [in] whether to reject the minimum and maximum values */
885 : )
886 : {
887 : typedef typename imageTout::Scalar scalarT;
888 :
889 3 : xMax = -1;
890 3 : yMax = -1;
891 3 : pMax = std::numeric_limits<scalarT>::lowest();
892 :
893 6 : if( meanFullWidth <= 0 || meanFullWidth > imIn.rows() || meanFullWidth > imIn.cols() ||
894 9 : imOut.rows() != imIn.rows() || imOut.cols() != imIn.cols() )
895 : {
896 0 : return -1;
897 : }
898 :
899 3 : const int before = meanFullWidth / 2;
900 3 : const int after = meanFullWidth - before - 1;
901 3 : int nPix = meanFullWidth * meanFullWidth;
902 :
903 3 : if( rejectMinMax && nPix <= 2 )
904 : {
905 0 : return -1;
906 : }
907 :
908 3 : if( rejectMinMax ) // avoid the branch on every pixel.
909 : {
910 3 : nPix -= 2;
911 6 : for( int jj = before; jj < imIn.cols() - after; ++jj )
912 : {
913 6 : for( int ii = before; ii < imIn.rows() - after; ++ii )
914 : {
915 3 : scalarT sum = 0;
916 3 : scalarT max = std::numeric_limits<scalarT>::lowest();
917 3 : scalarT min = std::numeric_limits<scalarT>::max();
918 13 : for( int ll = 0; ll < meanFullWidth; ++ll )
919 : {
920 46 : for( int kk = 0; kk < meanFullWidth; ++kk )
921 : {
922 36 : const scalarT value = imIn( ii - before + kk, jj - before + ll );
923 36 : sum += value;
924 36 : if( value > max )
925 35 : max = value;
926 36 : if( value < min )
927 3 : min = value;
928 : }
929 : }
930 3 : imOut( ii, jj ) = ( sum - max - min ) / nPix;
931 3 : if( imOut( ii, jj ) > pMax )
932 : {
933 3 : pMax = imOut( ii, jj );
934 3 : xMax = ii;
935 3 : yMax = jj;
936 : }
937 : }
938 : }
939 : }
940 : else
941 : {
942 0 : for( int jj = before; jj < imIn.cols() - after; ++jj )
943 : {
944 0 : for( int ii = before; ii < imIn.rows() - after; ++ii )
945 : {
946 0 : scalarT sum = 0;
947 0 : for( int ll = 0; ll < meanFullWidth; ++ll )
948 : {
949 0 : for( int kk = 0; kk < meanFullWidth; ++kk )
950 : {
951 0 : sum += imIn( ii - before + kk, jj - before + ll );
952 : }
953 : }
954 0 : imOut( ii, jj ) = sum / nPix;
955 0 : if( imOut( ii, jj ) > pMax )
956 : {
957 0 : pMax = imOut( ii, jj );
958 0 : xMax = ii;
959 0 : yMax = jj;
960 : }
961 : }
962 : }
963 : }
964 :
965 3 : return 0;
966 : }
967 :
968 : /// Smooth an image using the median in a rectangular box. Also Determines the location and value of the highest pixel
969 : /// in the smoothed image.
970 : /** Calculates the median value in a rectangular box of imIn, of size medianFullWidth X medianFullWidth and stores it in
971 : * the corresponding center pixel of imOut. For even widths, the window is associated with the higher-index member of
972 : * the central pair: it contains medianFullWidth/2 pixels before the output pixel and one fewer after it. The median of
973 : * an even-sized window is the arithmetic mean of its two central samples. Pixels without a complete window are not
974 : * changed in imOut (i.e. you should initialize them before the call).
975 : *
976 : * imOut is not re-allocated.
977 : *
978 : * Also determines the location and value of the maximum pixel. This is a negligble overhead compared to the median
979 : * operation.
980 : *
981 : *
982 : * \tparam imageTout is an eigen-like image array
983 : * \tparam imageTin is an eigen-like image array
984 : *
985 : * \returns 0 on success
986 : * \returns -1 on error.
987 : *
988 : * \ingroup image_filters_average
989 : */
990 : template <typename imageTout, typename imageTin>
991 7 : int medianSmooth(
992 : imageTout &imOut, /**< [out] the smoothed image. Not re-allocated, and the edge pixels are not modified */
993 : int &xMax, /**< [out] the x location of the maximum pixel */
994 : int &yMax, /**< [out] the y location of the maximum pixel */
995 : typename imageTout::Scalar &pMax, /**< [out] the value of the maximum pixel */
996 : const imageTin &imIn, /**< [in] the image to smooth */
997 : int medianFullWidth /**< [in] the full width of the smoothing box */
998 : )
999 : {
1000 : typedef typename imageTout::Scalar scalarT;
1001 :
1002 7 : xMax = -1;
1003 7 : yMax = -1;
1004 7 : pMax = std::numeric_limits<scalarT>::lowest();
1005 :
1006 9 : if( medianFullWidth <= 0 || medianFullWidth > imIn.rows() || medianFullWidth > imIn.cols() ||
1007 16 : imOut.rows() != imIn.rows() || imOut.cols() != imIn.cols() )
1008 : {
1009 4 : return -1;
1010 : }
1011 :
1012 3 : const int before = medianFullWidth / 2;
1013 3 : const int after = medianFullWidth - before - 1;
1014 3 : const size_t sampleCount = static_cast<size_t>( medianFullWidth ) * static_cast<size_t>( medianFullWidth );
1015 3 : std::vector<scalarT> pixs( sampleCount );
1016 :
1017 8 : for( int jj = before; jj < imIn.cols() - after; ++jj )
1018 : {
1019 13 : for( int ii = before; ii < imIn.rows() - after; ++ii )
1020 : {
1021 8 : size_t n = 0;
1022 21 : for( int ll = 0; ll < medianFullWidth; ++ll )
1023 : {
1024 44 : for( int kk = 0; kk < medianFullWidth; ++kk )
1025 : {
1026 31 : pixs[n] = imIn( ii - before + kk, jj - before + ll );
1027 31 : ++n;
1028 : }
1029 : }
1030 :
1031 8 : imOut( ii, jj ) = math::vectorMedianInPlace( pixs );
1032 8 : if( imOut( ii, jj ) > pMax )
1033 : {
1034 6 : pMax = imOut( ii, jj );
1035 6 : xMax = ii;
1036 6 : yMax = jj;
1037 : }
1038 : }
1039 : }
1040 :
1041 3 : return 0;
1042 3 : }
1043 :
1044 : /// Smooth an image using the median in a rectangular box.
1045 : /** Calculates the median value in a rectangular box of imIn, of size medianFullWidth X medianFullWidth and stores it in
1046 : * the corresponding center pixel of imOut. Even widths use the higher-index member of the central pair as the output
1047 : * pixel. Pixels without a complete window are not changed in imOut (i.e. you should initialize them before the call).
1048 : *
1049 : * imOut is not re-allocated.
1050 : *
1051 : * \overload
1052 : *
1053 : * \tparam imageTout is an eigen-like image array
1054 : * \tparam imageTin is an eigen-like image array
1055 : *
1056 : * \returns 0 on success
1057 : * \returns -1 on error.
1058 : *
1059 : * \ingroup image_filters_average
1060 : */
1061 : template <typename imageTout, typename imageTin>
1062 2 : int medianSmooth(
1063 : imageTout &imOut, /**< [out] the smoothed image. Not re-allocated, and the edge pixels are not modified */
1064 : const imageTin &imIn, /**< [in] the image to smooth */
1065 : int medianFullWidth /**< [in] the full width of the smoothing box */
1066 : )
1067 : {
1068 : int xMax, yMax;
1069 : typename imageTout::Scalar pMax;
1070 4 : return medianSmooth( imOut, xMax, yMax, pMax, imIn, medianFullWidth );
1071 : }
1072 :
1073 : template <typename eigenImT>
1074 : void rowEdgeMedSubtract( eigenImT &im, ///< The image to filter
1075 : int ncols ///< The number of columns on each side of the image to use as the reference
1076 : )
1077 : {
1078 : typedef typename eigenImT::Scalar realT;
1079 :
1080 : std::vector<realT> edge( 2 * ncols );
1081 : for( int rr = 0; rr < im.rows(); ++rr )
1082 : {
1083 : for( int cc = 0; cc < ncols; ++cc )
1084 : {
1085 : edge[cc] = im( rr, cc );
1086 : }
1087 :
1088 : for( int cc = 0; cc < ncols; ++cc )
1089 : {
1090 : edge[ncols + cc] = im( rr, im.cols() - ncols + cc );
1091 : }
1092 :
1093 : realT med = math::vectorMedian( edge );
1094 : im.row( rr ) -= med;
1095 : }
1096 :
1097 : return;
1098 : }
1099 :
1100 : template <typename eigenImT>
1101 : void colEdgeMedSubtract( eigenImT &im, ///< The image to filter
1102 : int nrows ///< The number of rows on each side of the image to use as the reference
1103 : )
1104 : {
1105 : typedef typename eigenImT::Scalar realT;
1106 :
1107 : std::vector<realT> edge( 2 * nrows );
1108 : for( int cc = 0; cc < im.cols(); ++cc )
1109 : {
1110 : for( int rr = 0; rr < nrows; ++rr )
1111 : {
1112 : edge[rr] = im( rr, cc );
1113 : }
1114 :
1115 : for( int rr = 0; rr < nrows; ++rr )
1116 : {
1117 : edge[nrows + rr] = im( im.rows() - nrows + rr, cc );
1118 : }
1119 :
1120 : realT med = math::vectorMedian( edge );
1121 : im.col( cc ) -= med;
1122 : }
1123 :
1124 : return;
1125 : }
1126 :
1127 : //------------ Radial Profile --------------------//
1128 :
1129 : template <typename floatT>
1130 : struct radval
1131 : {
1132 : floatT r;
1133 : floatT v;
1134 : };
1135 :
1136 : template <typename floatT>
1137 : struct radvalRadComp
1138 : {
1139 2307 : bool operator()( radval<floatT> rv1, radval<floatT> rv2 )
1140 : {
1141 2307 : return ( rv1.r < rv2.r );
1142 : }
1143 : };
1144 :
1145 : template <typename floatT>
1146 : struct radvalValComp
1147 : {
1148 439 : bool operator()( radval<floatT> rv1, radval<floatT> rv2 )
1149 : {
1150 439 : return ( rv1.v < rv2.v );
1151 : }
1152 : };
1153 :
1154 : /// Calculate the the radial profile
1155 : /** The median radial profile is calculated by rebinning to a 1 pixel grid.
1156 : *
1157 : *
1158 : * \tparam vecT the std::vector-like type to contain the profile
1159 : * \tparam eigenImT1 the eigen-array-like type of the input image
1160 : * \tparam eigenImT2 the eigen-array-like type of the radius and mask image
1161 : * \tparam eigenImT3 the eigen-array-like type of the mask image
1162 : *
1163 : * \ingroup rad_prof
1164 : */
1165 : template <typename vecT, typename eigenImT1, typename eigenImT2, typename eigenImT3>
1166 2 : void radprof(
1167 : vecT &rad, ///< [out] the radius points for the profile. Should be empty.
1168 : vecT &prof, ///< [out] the median image value at the corresponding radius. Should be empty.
1169 : const eigenImT1 &im, ///< [in] the image of which to calculate the profile
1170 : const eigenImT2 &radim, ///< [in] image of radius values per pixel
1171 : const eigenImT3 *mask, ///< [in] [optional] 1/0 mask, only pixels with a value of 1 are included in the profile. Set
1172 : ///< to 0 to not use.
1173 : bool mean = false, ///< [in] [optional] set to true to use the mean. If false (default) the median is used.
1174 : typename eigenImT1::Scalar minr = 0 )
1175 : {
1176 : typedef typename eigenImT1::Scalar floatT;
1177 :
1178 2 : int dim1 = im.rows();
1179 2 : int dim2 = im.cols();
1180 :
1181 : floatT maxr;
1182 : // floatT minr;
1183 :
1184 : size_t nPix;
1185 2 : if( mask )
1186 : {
1187 1 : nPix = mask->sum();
1188 : }
1189 : else
1190 : {
1191 1 : nPix = dim1 * dim2;
1192 : }
1193 :
1194 : /* A vector of radvals will be sorted, then binned*/
1195 2 : std::vector<radval<floatT>> rv( nPix );
1196 :
1197 2 : size_t i = 0;
1198 :
1199 20 : for( int c = 0; c < im.cols(); ++c )
1200 : {
1201 180 : for( int r = 0; r < im.rows(); ++r )
1202 : {
1203 162 : if( mask )
1204 : {
1205 81 : if( ( *mask )( r, c ) == 0 )
1206 1 : continue;
1207 : }
1208 :
1209 161 : rv[i].r = radim( r, c );
1210 161 : rv[i].v = im( r, c );
1211 161 : ++i;
1212 : }
1213 : }
1214 :
1215 2 : sort( rv.begin(), rv.end(), radvalRadComp<floatT>() );
1216 :
1217 : // for(auto it=rv.begin(); it != rv.end(); ++it)
1218 : // {
1219 : // std::cout << it->r << " " << it->v << "\n";
1220 : // }
1221 : //
1222 : // exit(0);
1223 :
1224 2 : if( minr == 0 )
1225 2 : minr = rv[0].r;
1226 2 : maxr = rv.back().r;
1227 :
1228 : /*Now bin*/
1229 2 : floatT dr = 1;
1230 2 : floatT r0 = minr;
1231 2 : floatT r1 = minr + dr;
1232 2 : int i1 = 0, i2, n;
1233 :
1234 : floatT med;
1235 :
1236 11 : while( r1 < maxr )
1237 : {
1238 9 : while( rv[i1].r < r0 )
1239 0 : ++i1;
1240 9 : i2 = i1;
1241 162 : while( rv[i2].r <= r1 )
1242 153 : ++i2;
1243 :
1244 9 : if( mean )
1245 : {
1246 0 : med = 0;
1247 0 : for( int in = i1; in < i2; ++in )
1248 0 : med += rv[in].v;
1249 0 : med /= ( i2 - i1 );
1250 : }
1251 : else
1252 : {
1253 9 : n = 0.5 * ( i2 - i1 );
1254 :
1255 45 : std::nth_element( rv.begin() + i1, rv.begin() + i1 + n, rv.begin() + i2, radvalValComp<floatT>() );
1256 :
1257 27 : med = ( rv.begin() + i1 + n )->v;
1258 :
1259 : // Average two points if even number of points
1260 9 : if( ( i2 - i1 ) % 2 == 0 )
1261 : {
1262 8 : med =
1263 8 : 0.5 *
1264 40 : ( med + ( *std::max_element( rv.begin() + i1, rv.begin() + i1 + n, radvalValComp<floatT>() ) ).v );
1265 : }
1266 : }
1267 :
1268 9 : rad.push_back( .5 * ( r0 + r1 ) );
1269 9 : prof.push_back( med );
1270 9 : i1 = i2;
1271 9 : r0 += dr;
1272 9 : r1 += dr;
1273 : }
1274 2 : }
1275 :
1276 : /// Calculate the the radial profile
1277 : /** The median radial profile is calculated by rebinning to a 1 pixel grid.
1278 : * This version calculates a centered radius image.
1279 : *
1280 : * \overload
1281 : *
1282 : * \tparam vecT the std::vector-like type to contain the profile
1283 : * \tparam eigenImT1 the eigen-array-like type of the input image
1284 : * \tparam eigenImT2 the eigen-array-like type of the radius and mask image
1285 : * \tparam eigenImT3 the eigen-array-like type of the mask image
1286 : *
1287 : * \ingroup rad_prof
1288 : */
1289 : template <typename vecT, typename eigenImT1, typename eigenImT2>
1290 : void radprof(
1291 : vecT &rad, ///< [out] the radius points for the profile. Should be empty.
1292 : vecT &prof, ///< [out] the median image value at the corresponding radius. Should be empty.
1293 : const eigenImT1 &im, ///< [in] the image of which to calculate the profile
1294 : const eigenImT2 &mask, ///< [in] 1/0 mask, only pixels with a value of 1 are included in the profile
1295 : bool mean = false ///< [in] [optional] set to true to use the mean. If false (default) the median is used.
1296 : )
1297 : {
1298 : eigenImage<typename eigenImT1::Scalar> radim;
1299 : radim.resize( im.cols(), im.rows() );
1300 :
1301 : radiusImage( radim );
1302 :
1303 : radprof( rad, prof, im, radim, &mask, mean );
1304 : }
1305 :
1306 : /// Calculate the the radial profile
1307 : /** The median radial profile is calculated by rebinning to a 1 pixel grid.
1308 : * This version calculates a centered radius image.
1309 : *
1310 : * \overload
1311 : *
1312 : * \tparam vecT the std::vector-like type to contain the profile
1313 : * \tparam eigenImT1 the eigen-array-like type of the input image
1314 : *
1315 : * \ingroup rad_prof
1316 : */
1317 : template <typename vecT, typename eigenImT1>
1318 : void radprof(
1319 : vecT &rad, ///< [out] the radius points for the profile. Should be empty.
1320 : vecT &prof, ///< [out] the median image value at the corresponding radius. Should be empty.
1321 : const eigenImT1 &im, ///< [in] the image of which to calculate the profile
1322 : bool mean = false, ///< [in] [optional] set to true to use the mean. If false (default) the median is used.
1323 : double dr = 1 )
1324 : {
1325 : eigenImage<typename eigenImT1::Scalar> radim;
1326 : radim.resize( im.cols(), im.rows() );
1327 :
1328 : radiusImage( radim );
1329 :
1330 : radprof( rad, prof, im, radim, (eigenImage<typename eigenImT1::Scalar> *)nullptr, mean, dr );
1331 : }
1332 :
1333 : /// Form a radial profile image, and optionally subtract it from the input
1334 : /** The radial profile is calculated using linear interpolation on a 1 pixel grid
1335 : *
1336 : *
1337 : * \tparam radprofT the eigen array type of the output
1338 : * \tparam eigenImT1 the eigen array type of the input image
1339 : * \tparam eigenImT2 the eigen array type of the radius image
1340 : * \tparam eigenImT3 the eigen array type of the mask image
1341 : *
1342 : * \ingroup rad_prof
1343 : */
1344 : template <typename radprofT, typename eigenImT1, typename eigenImT2, typename eigenImT3>
1345 2 : void radprofim( radprofT &radprofIm, ///< [out] the radial profile image. This will be resized.
1346 : eigenImT1 &im, ///< [in the image to form the profile of.
1347 : const eigenImT2 &rad, ///< [in] an array of radius values for each pixel
1348 : const eigenImT3 *mask, /**< [in] [optional 1/0 mask, only pixels with a value of 1 are
1349 : included in the profile. Can be nullptr. */
1350 : bool subtract, /**< [in] if true, then on ouput im will have had its radial
1351 : profile subtracted. */
1352 : bool mean = false /**< [in] [optional] set to true to use the mean.
1353 : If false (default) the median is used. */)
1354 : {
1355 :
1356 2 : std::vector<double> med_r, med_v; // Must be double for GSL interpolator
1357 :
1358 2 : radprof( med_r, med_v, im, rad, mask );
1359 :
1360 : /* And finally, interpolate onto the radius image */
1361 2 : radprofIm.resize( im.rows(), im.cols() );
1362 :
1363 2 : math::gslInterpolator<math::gsl_interp_linear<double>> interp( med_r, med_v );
1364 :
1365 20 : for( int c = 0; c < im.cols(); ++c )
1366 : {
1367 180 : for( int r = 0; r < im.rows(); ++r )
1368 : {
1369 162 : if( mask )
1370 : {
1371 81 : if( ( *mask )( r, c ) == 0 )
1372 : {
1373 1 : radprofIm( r, c ) = 0;
1374 1 : continue;
1375 : }
1376 : }
1377 :
1378 161 : radprofIm( r, c ) = interp( ( (double)rad( r, c ) ) );
1379 161 : if( subtract )
1380 161 : im( r, c ) -= radprofIm( r, c );
1381 : }
1382 : }
1383 2 : }
1384 :
1385 : /// Form a radial profile image, and optionally subtract it from the input
1386 : /** The radial profile is calculated using linear interpolation on a 1 pixel grid.
1387 : * This version calculates a centered radius image.
1388 : *
1389 : * \tparam radprofT the eigen array type of the output
1390 : * \tparam eigenImT the eigen array type of the input
1391 : *
1392 : * \ingroup rad_prof
1393 : */
1394 : template <typename radprofT, typename eigenImT>
1395 1 : void radprofim(
1396 : radprofT &radprof, ///< [out] the radial profile image. This will be resized.
1397 : eigenImT &im, ///< [in] the image to form the profile of.
1398 : bool subtract = false, ///< [in] [optional] if true, then on ouput im will have had its radial profile subtracted.
1399 : bool mean = false ///< [in] [optional] set to true to use the mean. If false (default) the median is used.
1400 : )
1401 : {
1402 1 : eigenImage<typename eigenImT::Scalar> rad;
1403 1 : rad.resize( im.rows(), im.cols() );
1404 :
1405 1 : radiusImage( rad );
1406 :
1407 1 : radprofim( radprof, im, rad, (eigenImage<typename eigenImT::Scalar> *)nullptr, subtract );
1408 1 : }
1409 :
1410 : /** \ingroup std_prof
1411 : * @{
1412 : */
1413 :
1414 : /// Form a standard deviation image, and optionally normalize the input relative to the local mean to form a S/N map.
1415 : /** The standard deviation profile is calculated using linear interpolation on a 1 pixel grid
1416 : *
1417 : * \tparam eigenImT the eigen array type of the output and non-reference images. Each image input can be a different
1418 : * type to allow references, etc.
1419 : *
1420 : */
1421 : template <typename eigenImT, typename eigenImT1, typename eigenImT2, typename eigenImT3>
1422 4 : void stddevImage( eigenImT &stdIm, ///< [out] the standard deviation image. This will be resized.
1423 : const eigenImT1 &im, ///< [in] the image to form the standard deviation profile of, never altered.
1424 : const eigenImT2 &rad, ///< [in] array of radius values
1425 : const eigenImT3 &mask, ///< [in] a 1/0 mask. 0 pixels are excluded from the std-dev calculations.
1426 : typename eigenImT::Scalar minRad, ///< [in] the minimum radius to analyze
1427 : typename eigenImT::Scalar maxRad, ///< [in] the maximum radius to analyze
1428 : bool divide ///< [in] if true, the output is the input image minus the interpolated mean profile,
1429 : ///< divided by the std-dev profile, i.e. a S/N map. default is false.
1430 : )
1431 : {
1432 : typedef typename eigenImT::Scalar floatT;
1433 :
1434 4 : int dim1 = im.cols();
1435 4 : int dim2 = im.rows();
1436 :
1437 4 : floatT mr = rad.maxCoeff();
1438 :
1439 : /* A vector of radvals will be sorted, then binned*/
1440 4 : std::vector<radval<floatT>> rv;
1441 4 : rv.reserve( dim1 * dim2 );
1442 :
1443 260 : for( int i = 0; i < im.size(); ++i )
1444 : {
1445 256 : if( mask( i ) == 0 )
1446 4 : continue;
1447 :
1448 252 : rv.push_back( { rad( i ), im( i ) } );
1449 : }
1450 :
1451 4 : sort( rv.begin(), rv.end(), radvalRadComp<floatT>() );
1452 :
1453 : /*Now bin*/
1454 4 : floatT dr = 1;
1455 4 : floatT r0 = 0;
1456 4 : floatT r1 = dr;
1457 4 : int i1 = 0, i2;
1458 :
1459 4 : std::vector<double> std_r, std_v, mean_v;
1460 :
1461 33 : while( r1 < mr )
1462 : {
1463 31 : while( i1 < rv.size() && rv[i1].r < r0 )
1464 : {
1465 1 : ++i1;
1466 : }
1467 30 : if( i1 == rv.size() )
1468 : {
1469 1 : break;
1470 : }
1471 :
1472 29 : i2 = i1;
1473 256 : while( i2 < rv.size() && rv[i2].r <= r1 )
1474 227 : ++i2;
1475 :
1476 29 : std::vector<double> vals;
1477 :
1478 256 : for( int i = i1; i < i2; ++i )
1479 : {
1480 227 : vals.push_back( rv[i].v );
1481 : }
1482 :
1483 29 : const double mean = math::vectorMean( vals );
1484 29 : std_r.push_back( .5 * ( r0 + r1 ) );
1485 29 : mean_v.push_back( mean );
1486 29 : std_v.push_back( std::sqrt( math::vectorVariance( vals, mean ) ) );
1487 29 : i1 = i2;
1488 29 : r0 += dr;
1489 29 : r1 += dr;
1490 : }
1491 :
1492 : /* And finally, interpolate onto the radius image */
1493 4 : stdIm.resize( dim1, dim2 );
1494 4 : math::gslInterpolator<math::gsl_interp_linear<double>> stdInterp( std_r, std_v );
1495 4 : math::gslInterpolator<math::gsl_interp_linear<double>> meanInterp( std_r, mean_v );
1496 :
1497 36 : for( int i = 0; i < dim1; ++i )
1498 : {
1499 288 : for( int j = 0; j < dim2; ++j )
1500 : {
1501 256 : if( rad( i, j ) < minRad || rad( i, j ) > maxRad )
1502 : {
1503 65 : stdIm( i, j ) = 0;
1504 : }
1505 : else
1506 : {
1507 191 : stdIm( i, j ) = stdInterp( ( (double)rad( i, j ) ) );
1508 191 : if( divide )
1509 191 : stdIm( i, j ) = ( im( i, j ) - meanInterp( ( (double)rad( i, j ) ) ) ) / stdIm( i, j );
1510 : }
1511 : }
1512 : }
1513 4 : }
1514 :
1515 : /// Form a standard deviation image, and optionally normalize the input relative to the local mean to form a S/N map.
1516 : /** The standard deviation profile is calculated using linear interpolation on a 1 pixel grid
1517 : *
1518 : * This version creates a radius map on each call, and calls the above version. This should not
1519 : * be used for repeated alls, rather create a radius map ahead of time.
1520 : *
1521 : * \overload
1522 : *
1523 : * \tparam eigenImT the eigen array type of the output and non-reference images
1524 : *
1525 : */
1526 : template <typename eigenImT, typename eigenImT1, typename eigenImT2>
1527 : void stddevImage(
1528 : eigenImT &stdIm, ///< [out] the standard deviation image. This will be resized.
1529 : const eigenImT1 &im, ///< [in] the image to form the standard deviation profile of, never altered.
1530 : const eigenImT2 &mask, ///< [in] a 1/0 mask. 0 pixels are excluded from the std-dev calculations.
1531 : typename eigenImT::Scalar minRad, ///< [in] the minimum radius to analyze
1532 : typename eigenImT::Scalar maxRad, ///< [in] the maximum radius to analyze
1533 : bool divide = false ///< [in] [optional] if true, the output is the input image minus the interpolated
1534 : ///< mean profile, divided by the std-dev profile, i.e. a S/N map. default is false.
1535 : )
1536 : {
1537 : int dim1 = im.cols();
1538 : int dim2 = im.rows();
1539 :
1540 : eigenImage<typename eigenImT::Scalar> rad;
1541 : rad.resize( dim1, dim2 );
1542 :
1543 : radiusImage( rad );
1544 : stddevImage( stdIm, im, rad, mask, minRad, maxRad, divide );
1545 : }
1546 :
1547 : /// Form a standard deviation image for each image in a cube, and optionally normalize it relative to the local mean
1548 : /// to form a S/N map cube.
1549 : /** The standard deviation profile is calculated using linear interpolation on a 1 pixel grid
1550 : *
1551 : *
1552 : * \tparam eigencubeT is the eigen cube type of the input and output cubes.
1553 : * \tparam eigenImT the eigen array type of the output and non-reference images.
1554 : *
1555 : */
1556 : template <typename eigenCubeT, typename eigenCubeT1, typename eigenCubeT2, typename radT1, typename radT2>
1557 : void stddevImageCube(
1558 : eigenCubeT &stdImc, /**< [out] the standard deviation image cube. This will be resized. */
1559 : const eigenCubeT1 &imc, /**< [in] the image cube to form the standard deviation profile of. */
1560 : const eigenCubeT2 &maskCube, /**< [in] a 1/0 mask. 0 pixels are excluded from the std-dev calculations. */
1561 : radT1 minRad, /**< [in] the minimum radius to analyze */
1562 : radT2 maxRad, /**< [in] the maximum radius to analyze */
1563 : bool divide = false /**< [in] [optional] if true, the output is the input image minus the interpolated
1564 : mean profile, divided by the std-dev profile, i.e. a S/N map. default is false.*/
1565 : )
1566 : {
1567 : int dim1 = imc.cols();
1568 : int dim2 = imc.rows();
1569 :
1570 : typename eigenCubeT::Scalar minRadF = minRad;
1571 : typename eigenCubeT::Scalar maxRadF = maxRad;
1572 :
1573 : eigenImage<typename eigenCubeT::Scalar> rad;
1574 : rad.resize( dim1, dim2 );
1575 :
1576 : radiusImage( rad );
1577 :
1578 : stdImc.resize( imc.rows(), imc.cols(), imc.planes() );
1579 :
1580 : // #pragma omp parallel for
1581 : for( int i = 0; i < imc.planes(); ++i )
1582 : {
1583 : eigenImage<typename eigenCubeT::Scalar> im, stdIm, mask;
1584 :
1585 : im = imc.image( i );
1586 : mask = maskCube.image( i );
1587 :
1588 : stddevImage( stdIm, im, rad, mask, minRadF, maxRadF, divide );
1589 :
1590 : stdImc.image( i ) = stdIm;
1591 : }
1592 : }
1593 :
1594 : ///@}
1595 :
1596 : } // namespace improc
1597 : } // namespace mx
1598 :
1599 : #endif //__imageFilters_hpp__
|
