Line data Source code
1 : /** \file imageUtils.hpp
2 : * \author Jared R. Males
3 : * \brief Header for the image processing utilities
4 : * \ingroup image_processing_files
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2020 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef improc_imageUtils_hpp
28 : #define improc_imageUtils_hpp
29 :
30 : #include <cmath>
31 : #include <limits>
32 : #include <type_traits>
33 :
34 : #include "../math/floatUtils.hpp"
35 :
36 : #include "imageTransforms.hpp"
37 :
38 : namespace mx
39 : {
40 : namespace improc
41 : {
42 :
43 : /// Get the mxlib standard center coordinate of an image
44 : /**
45 : * \returns 0.5*(rows_cols-1)
46 : */
47 : template <typename realT = double>
48 : realT imCen( int rows_cols /**< [in] the number of rows or columns in the image */ )
49 : {
50 : return 0.5 * ( 1.0 * rows_cols - 1 );
51 : }
52 :
53 : /// Get the mxlib standard center x-coordinate of an image
54 : /**
55 : * \returns 0.5*(rows-1)
56 : */
57 : template <typename realT = double, typename imT>
58 : realT imCenX( const imT &im /**< [in] the image to find the center of */ )
59 : {
60 : return imCen<realT>( im.rows() );
61 : }
62 :
63 : /// Get the mxlib standard center y-coordinate of an image
64 : /**
65 : * \returns 0.5*(cols-1)
66 : */
67 : template <typename realT = double, typename imT>
68 : realT imCenY( const imT &im /**< [in] the image to find the center of */ )
69 : {
70 : return imCen<realT>( im.cols() );
71 : }
72 :
73 : /** \ingroup image_utils
74 : *@{
75 : */
76 :
77 : /// Return the configured invalid image-pixel value.
78 : /** The default is a quiet NaN so FITS and image-display tools recognize invalid pixels. Defining
79 : * `MXLIB_INVALID_NUMBER_VALUE` at compile time selects that finite value instead. CMake propagates its
80 : * `MXLIB_INVALID_NUMBER_VALUE` cache setting through the installed pkg-config metadata.
81 : *
82 : * \returns the invalid-pixel value for `realT`
83 : */
84 : template <typename realT>
85 9 : constexpr realT invalidNumber()
86 : {
87 : static_assert( std::is_floating_point_v<realT>, "invalidNumber requires a floating-point type" );
88 :
89 : #ifdef MXLIB_INVALID_NUMBER_VALUE
90 : return static_cast<realT>( MXLIB_INVALID_NUMBER_VALUE );
91 : #else
92 9 : return std::numeric_limits<realT>::quiet_NaN();
93 : #endif
94 : }
95 :
96 : /// Check whether a value represents an invalid image pixel.
97 : /** Detects the configured mxlib invalid-number value as well as NaN and positive or negative infinity. The
98 : * non-finite classification remains reliable under finite-math-only optimization.
99 : *
100 : * \returns true if value is invalid, otherwise false.
101 : */
102 : template <typename realT>
103 198 : bool isInvalidPixel( realT value /**< [in] value to test */ )
104 : {
105 : static_assert( std::is_floating_point_v<realT>, "isInvalidPixel requires a floating-point type" );
106 :
107 : #ifdef MXLIB_INVALID_NUMBER_VALUE
108 : return value == invalidNumber<realT>() || !math::isFinite( value );
109 : #else
110 198 : return !math::isFinite( value );
111 : #endif
112 : }
113 :
114 : /// Reflect pixel coordinates across the given center pixel.
115 : /**
116 : */
117 : template <typename realT>
118 : int reflectImageCoords( int &x1, ///< [out] the reflected x coordinate
119 : int &y1, ///< [out] the reflected y coordinate
120 : int x0, ///< [in] the input x coordinate
121 : int y0, ///< [in] the input y coordinate
122 : realT xc, ///< [in] the center pixel x coordinate
123 : realT yc ///< [in] the center pixel y coordinate
124 : )
125 : {
126 : x1 = xc - ( x0 - xc );
127 : y1 = yc - ( y0 - yc );
128 :
129 : return 0;
130 : }
131 :
132 : /// @}
133 :
134 : /// Zero any NaNs in an image
135 : /**
136 : * \overload
137 : *
138 : * \ingroup eigen_image_processing
139 : */
140 : template <class imageT, typename valueT>
141 2 : void zeroNaNs( imageT &im, ///< [in.out] image which will have any NaN pixels set to zero
142 : valueT val ///< [in] [optional] The value to set NaN pixels too. Default is 0.
143 : )
144 : {
145 20 : for( int c = 0; c < im.cols(); ++c )
146 : {
147 180 : for( int r = 0; r < im.rows(); ++r )
148 : {
149 162 : if( isInvalidPixel( im( r, c ) ) )
150 : {
151 33 : im( r, c ) = val;
152 : }
153 : }
154 : }
155 2 : }
156 :
157 : /// Zero any NaNs in an image
158 : /**
159 : * \ingroup eigen_image_processing
160 : */
161 : template <class imageT>
162 2 : void zeroNaNs( imageT &im /**< [in.out] image which will have any NaN pixels set to zero */ )
163 : {
164 2 : typename imageT::Scalar zero = 0;
165 2 : zeroNaNs<imageT, typename imageT::Scalar>( im, zero );
166 2 : }
167 :
168 : /// Zero any NaNs in an image cube
169 : /** This version fills in a mask with 1s where there were nans, 0s elsewhere.
170 : *
171 : * \ingroup eigen_image_processing
172 : */
173 : template <class cubeT, class maskCubeT>
174 2 : void zeroNaNCube( cubeT &imc, /**< [in.out] cube which will have any NaN pixels set to zero */
175 : maskCubeT *mask /**< [out] a 1/0 mask with 1 indicating which pixels where nan */
176 : )
177 : {
178 2 : if( mask )
179 : {
180 1 : mask->resize( imc.rows(), imc.cols(), imc.planes() );
181 1 : mask->setZero();
182 : }
183 :
184 6 : for( int p = 0; p < imc.planes(); ++p )
185 : {
186 16 : for( int c = 0; c < imc.cols(); ++c )
187 : {
188 36 : for( int r = 0; r < imc.rows(); ++r )
189 : {
190 24 : if( isInvalidPixel( imc.image( p )( r, c ) ) )
191 : {
192 2 : imc.image( p )( r, c ) = 0;
193 2 : if( mask )
194 : {
195 2 : ( *mask ).image( p )( r, c ) = 1;
196 : }
197 : }
198 : }
199 : }
200 : }
201 2 : }
202 :
203 : /// Zero any NaNs in an image cube
204 : /**
205 : * \ingroup eigen_image_processing
206 : */
207 : template <class cubeT>
208 1 : void zeroNaNCube( cubeT &imc /**< [in.out] cube which will have any NaN pixels set to zero */ )
209 : {
210 1 : return zeroNaNCube<cubeT, cubeT>( imc, nullptr );
211 : }
212 :
213 : /// Calculate the mean value of an image
214 : /**
215 : *
216 : * \returns the mean of the input image
217 : *
218 : * \ingroup eigen_image_processing
219 : */
220 : template <class calcT, class imageT>
221 550 : calcT imageMean( imageT &im /**< [in] the image of which to calculate the mean*/ )
222 : {
223 550 : return static_cast<calcT>( im.sum() ) / ( im.rows() * im.cols() );
224 : }
225 :
226 : /// Calculate the mean value of an image over a mask
227 : /**
228 : * \returns the mean of the input image in the masked region
229 : *
230 : * \ingroup eigen_image_processing
231 : */
232 : template <class calcT, class imageT, class maskT>
233 0 : calcT imageMean( imageT &im, /**< [in] the image of which to calculate the mean*/
234 : const maskT &mask /**< [in] a 1/0 mask where 1 defines the good pixels */ )
235 : {
236 0 : return static_cast<calcT>( ( im * mask ).sum() ) / ( mask.sum() );
237 : }
238 :
239 : /// Calculate the variance of an image given its mean
240 : /**
241 : *
242 : * \returns the variance of the input image
243 : *
244 : * \ingroup eigen_image_processing
245 : */
246 : template <typename calcT, class imageT>
247 550 : calcT imageVariance( imageT &im /**< [in] the image of which to calculate the variance*/,
248 : calcT mn /**< [in] the mean value of the image w.r.t. which to calcualate the variance */
249 : )
250 : {
251 550 : return ( im.template cast<calcT>() - mn ).square().sum() / ( im.rows() * im.cols() );
252 : }
253 :
254 : /// Calculate the variance of an image given its mean
255 : /**
256 : *
257 : * \returns the variance of the input image
258 : *
259 : * \ingroup eigen_image_processing
260 : */
261 : template <typename calcT, class imageT, class maskT>
262 0 : calcT imageVariance( imageT &im /**< [in] the image of which to calculate the variance*/,
263 : calcT mn, /**< [in] the mean value of the image w.r.t. which to calcualate the variance */
264 : const maskT &mask /**< [in] a 1/0 mask where 1 defines the good pixels */
265 : )
266 : {
267 0 : return ( im.template cast<calcT>() * mask - mn ).square().sum() / ( mask.sum() );
268 : }
269 :
270 : /// Calculate the median of an Eigen-like array.
271 : /** Calculates the median of the entire array, allowing for some pixels to be ignored using a mask.
272 : * Working memory can be retained between calls.
273 : *
274 : * \tparam imageT is an Eigen-like type
275 : * \tparam maskT is an Eigen-like type
276 : *
277 : * \returns the median of the unmasked pixels of mat, using \ref mx::math::vectorMedianInPlace().
278 : *
279 : * \ingroup eigen_image_processing
280 : *
281 : */
282 : template <typename imageT, typename maskT = imageT>
283 : typename imageT::Scalar
284 8 : imageMedian( const imageT &mat, /**< [in] the image to take the median of*/
285 : const maskT *mask, /**< [in] if non-0, a 1/0 mask where 0 pixels are ignored.*/
286 : std::vector<typename imageT::Scalar> *work = 0 /**< [in] [optional] working memory
287 : can be retained
288 : and re-passed. Is resized.*/
289 : )
290 : {
291 : typename imageT::Scalar med;
292 :
293 8 : bool localWork = false;
294 8 : if( work == 0 )
295 : {
296 3 : work = new std::vector<typename imageT::Scalar>;
297 3 : localWork = true;
298 : }
299 :
300 8 : int sz = mat.size();
301 :
302 8 : if( mask )
303 : {
304 3 : sz = mask->sum();
305 : }
306 :
307 8 : work->resize( sz );
308 :
309 8 : if( mask )
310 : {
311 3 : int ii = 0;
312 11 : for( int i = 0; i < mat.rows(); ++i )
313 : {
314 24 : for( int j = 0; j < mat.cols(); ++j )
315 : {
316 16 : if( ( *mask )( i, j ) == 0 )
317 : {
318 8 : continue;
319 : }
320 8 : ( *work )[ii] = mat( i, j );
321 8 : ++ii;
322 : }
323 : }
324 : }
325 : else
326 : {
327 5 : int ii = 0;
328 21 : for( int i = 0; i < mat.rows(); ++i )
329 : {
330 40 : for( int j = 0; j < mat.cols(); ++j )
331 : {
332 24 : ( *work )[ii] = mat( i, j );
333 24 : ++ii;
334 : }
335 : }
336 : }
337 :
338 8 : med = math::vectorMedianInPlace( *work );
339 :
340 8 : if( localWork )
341 : {
342 3 : delete work;
343 : }
344 :
345 8 : return med;
346 : }
347 :
348 : /// Calculate the median of an Eigen-like array.
349 : /** Calculates the median of the entire array.
350 : * Working memory can be retained between calls.
351 : *
352 : * \tparam imageT is an Eigen-like type
353 : *
354 : * \returns the median of the unmasked pixels of mat, using \ref mx::math::vectorMedianInPlace().
355 : *
356 : * \ingroup eigen_image_processing
357 : *
358 : */
359 : template <typename imageT>
360 5 : typename imageT::Scalar imageMedian( const imageT &mat, /**< [in] the image to take the median of*/
361 : std::vector<typename imageT::Scalar> *work = 0 /**< [in] [optional] working memory
362 : can be retained and re-passed.*/
363 : )
364 : {
365 5 : return imageMedian( mat, static_cast<Eigen::Array<typename imageT::Scalar, -1, -1> *>( nullptr ), work );
366 : }
367 :
368 : /// Calculate the center of light of an image
369 : /** Note that the sum of the image should be > 0.
370 : *
371 : * \ingroup eigen_image_processing
372 : */
373 : template <typename imageT>
374 15 : int imageCenterOfLight( typename imageT::Scalar &x, ///< [out] the x coordinate of the center of light [pixels]
375 : typename imageT::Scalar &y, ///< [out] the y coordinate of hte center of light [pixels]
376 : const imageT &im ///< [in] the image to centroid
377 : )
378 : {
379 15 : x = 0;
380 15 : y = 0;
381 :
382 15 : typename imageT::Scalar sum = im.sum();
383 :
384 15 : if( sum == 0 )
385 : {
386 0 : x = 0.5 * ( im.rows() - 1.0 );
387 0 : y = 0.5 * ( im.cols() - 1.0 );
388 0 : return 0;
389 : }
390 :
391 481 : for( int j = 0; j < im.cols(); ++j )
392 : {
393 20878 : for( int i = 0; i < im.rows(); ++i )
394 : {
395 20412 : x += ( i + 1 ) * im( i, j );
396 20412 : y += ( j + 1 ) * im( i, j );
397 : }
398 : }
399 :
400 15 : x = x / sum - 1;
401 15 : y = y / sum - 1;
402 :
403 15 : return 0;
404 : }
405 :
406 : /// Find the maximum in an image at sub-pixel resolution by interpolation
407 : /** Uses imageMagnify() to expand the image to the desired scale. Because of the
408 : * scaling used in imageMagnify, the desired scale may not be exact. As a result
409 : * the actual scale is returned in scale_x and scale_y.
410 : *
411 : * \ingroup eigen_image_processing
412 : */
413 : template <typename floatT, typename imageT, typename magImageT, typename transformT>
414 : int imageMaxInterp( floatT &x, ///< [out] the x-position of the maximum, in pixels of the input image
415 : floatT &y, ///< [out] the y-position of the maximum, in pixels of the input image
416 : floatT &scale_x, ///< [in.out] the desired scale or resolution, in pixels < 1, in the x direction.
417 : ///< On output contains the actual scale calculated.
418 : floatT &scale_y, ///< [in.out] the desired scale or resolution, in pixels < 1, in the y direction.
419 : ///< On output contains the actual scale calculated.
420 : magImageT &magIm, ///< [in] the magnified image. This is used as working memory, will be resized.
421 : const imageT &im, ///< [in] the image to find the maximum of
422 : transformT trans ///< [in] the transform to use for interpolation
423 : )
424 : {
425 : floatT magSize_x = ceil( ( im.rows() - 1.0 ) / scale_x ) + 1;
426 : floatT magSize_y = ceil( ( im.cols() - 1.0 ) / scale_y ) + 1;
427 :
428 : floatT mag_x = ( (floatT)magSize_x - 1.0 ) / ( (floatT)im.rows() - 1.0 );
429 : floatT mag_y = ( (floatT)magSize_y - 1.0 ) / ( (floatT)im.cols() - 1.0 );
430 :
431 : scale_x = 1.0 / mag_x;
432 : scale_y = 1.0 / mag_y;
433 :
434 : magIm.resize( magSize_x, magSize_y );
435 :
436 : imageMagnify( magIm, im, trans );
437 :
438 : int ix, iy;
439 : magIm.maxCoeff( &ix, &iy );
440 : x = ix * scale_x;
441 : y = iy * scale_y;
442 :
443 : return 0;
444 : }
445 :
446 : /// Find the maximum in an image at sub-pixel resolution by cubic convolution interpolation
447 : /** Uses imageMagnify() to expand the image to the desired scale. Because of the
448 : * scaling used in imageMagnify, the desired scale may not be exact. As a result
449 : * the actual scale is returned in scale_x and scale_y.
450 : *
451 : * \ingroup eigen_image_processing
452 : */
453 : template <typename floatT, typename imageT, typename magImageT>
454 : int imageMaxInterp( floatT &x, ///< [out] the x-position of the maximum, in pixels of the input image
455 : floatT &y, ///< [out] the y-position of the maximum, in pixels of the input image
456 : floatT &scale_x, ///< [in.out] the desired scale or resolution, in pixels < 1, in the x direction.
457 : ///< On output contains the actual scale calculated.
458 : floatT &scale_y, ///< [in.out] the desired scale or resolution, in pixels < 1, in the y direction.
459 : ///< On output contains the actual scale calculated.
460 : magImageT &magIm, ///< [in] the magnified image. This is used as working memory, will be resized.
461 : const imageT &im ///< [in] the image to find the maximum of
462 : )
463 : {
464 : return imageMaxInterp( x, y, scale_x, scale_y, magIm, im, cubicConvolTransform<typename imageT::Scalar>() );
465 : }
466 :
467 : /// Combine two images, each with their own mask defining good pixels.
468 : /** The combined image is made up of the pixels in im1 where mask1 is 1, and the pixels of im2 where mask2 is 1.
469 : * If a pixel in both mask1 and mask2 has a value of 1, that pixel in the combo is the average of im1 and im2.
470 : * All other pixels are set to 0 in the combined image.
471 : *
472 : * Separate template types are used for each argument to allow references, etc.
473 : *
474 : * \tparam imageT the eigen-like array type of the combined image
475 : * \tparam imageT1 the eigen-like array type of image 1
476 : * \tparam imageT2 the eigen-like array type of mask 1
477 : * \tparam imageT3 the eigen-like array type of image 2
478 : * \tparam imageT4 the eigen-like array type of mask 2
479 : *
480 : * \ingroup eigen_image_processing
481 : */
482 : template <typename imageT, typename imageT1, typename imageT2, typename imageT3, typename imageT4>
483 : void combine2ImagesMasked( imageT &combo, ///< [out] the combined image. will be resized.
484 : const imageT1 &im1, ///< [in] the first image
485 : const imageT2 &mask1, ///< [in] the mask for the first image
486 : const imageT3 &im2, ///< [in] the second image
487 : const imageT4 &mask2 ///< [in] the mask for the second image
488 : )
489 : {
490 : combo.resize( im1.rows(), im2.cols() );
491 :
492 : for( int c = 0; c < combo.cols(); ++c )
493 : {
494 : for( int r = 0; r < combo.rows(); ++r )
495 : {
496 : if( mask1( r, c ) == 1 && mask2( r, c ) == 0 )
497 : combo( r, c ) = im1( r, c );
498 : else if( mask2( r, c ) == 1 & mask1( r, c ) == 0 )
499 : combo( r, c ) = im2( r, c );
500 : else if( mask1( r, c ) == 1 && mask2( r, c ) == 1 )
501 : combo( r, c ) = 0.5 * ( im1( r, c ) + im2( r, c ) );
502 : else
503 : combo( r, c ) = 0;
504 : }
505 : }
506 : }
507 :
508 : /// Remove rows and columns
509 : /**
510 : * \ingroup eigen_image_processing
511 : */
512 : template <typename eigenT, typename eigenTin>
513 : void removeRowsAndCols( eigenT &out, const eigenTin &in, int st, int w )
514 : {
515 : out.resize( in.rows() - w, in.cols() - w );
516 :
517 : out.topLeftCorner( st, st ) = in.topLeftCorner( st, st );
518 :
519 : out.bottomLeftCorner( in.rows() - ( st + w ), st ) = in.bottomLeftCorner( in.rows() - ( st + w ), st );
520 :
521 : out.topRightCorner( st, in.cols() - ( st + w ) ) = in.topRightCorner( st, in.cols() - ( st + w ) );
522 :
523 : out.bottomRightCorner( in.rows() - ( st + w ), in.cols() - ( st + w ) ) =
524 : in.bottomRightCorner( in.rows() - ( st + w ), in.cols() - ( st + w ) );
525 : }
526 :
527 : /// Remove rows
528 : /**
529 : * \ingroup eigen_image_processing
530 : */
531 : template <typename eigenT, typename eigenTin>
532 : void removeRows( eigenT &out, const eigenTin &in, int st, int w )
533 : {
534 : out.resize( in.rows() - w, in.cols() );
535 :
536 : out.topLeftCorner( st, in.cols() ) = in.topLeftCorner( st, in.cols() );
537 :
538 : out.bottomLeftCorner( in.rows() - ( st + w ), in.cols() ) =
539 : in.bottomLeftCorner( in.rows() - ( st + w ), in.cols() );
540 : }
541 :
542 : /// Remove columns
543 : /**
544 : * \ingroup eigen_image_processing
545 : */
546 : template <typename eigenT, typename eigenTin>
547 : void removeCols( eigenT &out, const eigenTin &in, int st, int w )
548 : {
549 : out.resize( in.rows(), in.cols() - w );
550 :
551 : out.topLeftCorner( in.rows(), st ) = in.topLeftCorner( in.rows(), st );
552 :
553 : out.topRightCorner( in.rows(), in.cols() - ( st + w ) ) = in.topRightCorner( in.rows(), in.cols() - ( st + w ) );
554 : }
555 :
556 : /** \ingroup image_utils
557 : *@{
558 : */
559 :
560 : /// Copy one image to another, with no transformation
561 : /** This is merely memcpy
562 : *
563 : * \returns dest
564 : */
565 : void *imcpy( void *dest, ///< [out] the address of the first pixel in the destination image
566 : void *src, ///< [in] the address of the first pixel in the source image
567 : size_t width, ///< [in] the width in pixels of size szof
568 : size_t height, ///< [in] the height in pixels of size szof
569 : size_t szof ///< [in] the size in bytes of a one pixel
570 : );
571 :
572 : /// Copy one image to another, flipping up-down
573 : /** This is a reversed row-by-row memcpy
574 : *
575 : * \returns dest
576 : */
577 : void *imcpy_flipUD( void *dest, ///< [out] the address of the first pixel in the destination image
578 : void *src, ///< [in] the address of the first pixel in the source image
579 : size_t width, ///< [in] the width in pixels of size szof
580 : size_t height, ///< [in] the height in pixels of size szof
581 : size_t szof ///< [in] the size in bytes of a one pixel
582 : );
583 :
584 : /// Copy one image to another, flipping left-right
585 : /**
586 : *
587 : * \returns dest
588 : */
589 : void *imcpy_flipLR( void *dest, ///< [out] the address of the first pixel in the destination image
590 : void *src, ///< [in] the address of the first pixel in the source image
591 : size_t width, ///< [in] the width in pixels of size szof
592 : size_t height, ///< [in] the height in pixels of size szof
593 : size_t szof ///< [in] the size in bytes of a one pixel
594 : );
595 :
596 : /// Copy one image to another, flipping up-down and left-right
597 : /**
598 : *
599 : * \returns dest
600 : */
601 : void *imcpy_flipUDLR( void *dest, ///< [out] the address of the first pixel in the destination image
602 : void *src, ///< [in] the address of the first pixel in the source image
603 : size_t width, ///< [in] the width in pixels of size szof
604 : size_t height, ///< [in] the height in pixels of size szof
605 : size_t szof ///< [in] the size in bytes of a one pixel
606 : );
607 :
608 : } // namespace improc
609 : } // namespace mx
610 :
611 : #endif // improc_imageUtils_hpp
|
