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