Line data Source code
1 : /** \file imageTransforms.hpp
2 : * \author Jared R. Males
3 : * \brief Image interpolation and transformation
4 : * \ingroup image_processing_files
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015, 2016, 2017, 2018 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_imageTransforms_hpp
28 : #define improc_imageTransforms_hpp
29 :
30 : #include <cstddef>
31 : #include <cmath>
32 :
33 : #include <iostream>
34 : #include <limits>
35 :
36 : #include "eigenImage.hpp"
37 :
38 : namespace mx
39 : {
40 : namespace improc
41 : {
42 :
43 : /// Transformation by bi-linear interpolation
44 : /** \ingroup image_transforms
45 : */
46 : template <typename _arithT>
47 : struct bilinearTransform
48 : {
49 : typedef _arithT arithT;
50 :
51 : static const size_t width = 2;
52 : static const size_t lbuff = 0;
53 :
54 : template <typename arrT, typename arithT>
55 : void operator()( arrT &kern, arithT x, arithT y )
56 : {
57 : kern.resize( width, width );
58 :
59 : kern( 0, 0 ) = ( 1. - x ) * ( 1. - y );
60 : kern( 0, 1 ) = ( 1. - x ) * y;
61 : kern( 1, 0 ) = x * ( 1. - y );
62 : kern( 1, 1 ) = x * y;
63 : }
64 : };
65 :
66 : /// Typedef for bilinearTransform with single precision
67 : /** \ingroup image_transforms
68 : */
69 : typedef bilinearTransform<float> bilinearTransf;
70 :
71 : /// Typedef for bilinearTransform with double precision
72 : /** \ingroup image_transforms
73 : */
74 : typedef bilinearTransform<double> bilinearTransd;
75 :
76 : /// Transformation by cubic convolution interpolation
77 : /** Uses the cubic convolution interpolation kernel. See <a
78 : * href="https://en.wikipedia.org/wiki/Bicubic_interpolation">https://en.wikipedia.org/wiki/Bicubic_interpolation </a>.
79 : *
80 : * The parameter \ref cubic should be left as the default -0.5 in most cases, which gives the bicubic spline
81 : * interpolator. See <a
82 : * href="https://en.wikipedia.org/wiki/Cubic_Hermite_spline">https://en.wikipedia.org/wiki/Cubic_Hermite_spline</a>.
83 : *
84 : * \tparam _arithT is the type in which to do all calculations. Should be a floating point type.
85 : *
86 : * "[test doc]"
87 : *
88 : * \ingroup image_transforms
89 : */
90 : template <typename _arithT>
91 : struct cubicConvolTransform
92 : {
93 : typedef _arithT arithT; ///< The type in which all calculations are performed.
94 :
95 : static const size_t width = 4;
96 : static const size_t lbuff = 1;
97 :
98 : arithT cubic{ -0.5 }; ///< The kernel parameter. The default value -0.5 gives the bicubic spline interpolator.
99 :
100 : /// Default c'tor.
101 : /**
102 : * This will provide the bicubic spline interpolator
103 : */
104 19 : cubicConvolTransform()
105 19 : {
106 19 : }
107 :
108 : /// Construct setting the kernel parameter.
109 0 : explicit cubicConvolTransform( arithT c /**< [in] [optiona] The kernel parameter. The default value -0.5 gives the bicubic spline interpolator. */)
110 0 : {
111 0 : cubic = c;
112 0 : }
113 :
114 : /// Copy c'tor
115 0 : cubicConvolTransform( const cubicConvolTransform &t )
116 0 : {
117 0 : cubic = t.cubic;
118 0 : }
119 :
120 : /// Calculate the kernel value for a given residual.
121 419329 : arithT cubicConvolKernel( arithT d )
122 : {
123 419329 : if( d <= 1 )
124 : {
125 221228 : return ( cubic + 2. ) * d * d * d - ( cubic + 3. ) * d * d + 1.;
126 : }
127 :
128 198101 : if( d < 2 )
129 : {
130 186537 : return cubic * d * d * d - 5. * cubic * d * d + 8. * cubic * d - 4. * cubic;
131 : }
132 :
133 11564 : return 0;
134 : }
135 :
136 : ///\todo why is this arithT not just the class's arithT?
137 : template <typename arrT, typename arithT>
138 52415 : void operator()( arrT &kern, arithT x, arithT y )
139 : {
140 : arithT km2x, km1x, kp1x, kp2x;
141 : arithT km2y, km1y, kp1y, kp2y;
142 :
143 52415 : km2x = cubicConvolKernel( ( 1. + x ) );
144 52415 : km1x = cubicConvolKernel( x );
145 52415 : kp1x = cubicConvolKernel( 1. - x );
146 52415 : kp2x = cubicConvolKernel( 2. - x );
147 :
148 52415 : km2y = cubicConvolKernel( ( 1. + y ) );
149 52415 : km1y = cubicConvolKernel( y );
150 52415 : kp1y = cubicConvolKernel( 1. - y );
151 52415 : kp2y = cubicConvolKernel( 2. - y );
152 :
153 52415 : kern( 0, 0 ) = km2x * km2y;
154 52415 : kern( 0, 1 ) = km2x * km1y;
155 52415 : kern( 0, 2 ) = km2x * kp1y;
156 52415 : kern( 0, 3 ) = km2x * kp2y;
157 :
158 52415 : kern( 1, 0 ) = km1x * km2y;
159 52415 : kern( 1, 1 ) = km1x * km1y;
160 52415 : kern( 1, 2 ) = km1x * kp1y;
161 52415 : kern( 1, 3 ) = km1x * kp2y;
162 :
163 52415 : kern( 2, 0 ) = kp1x * km2y;
164 52415 : kern( 2, 1 ) = kp1x * km1y;
165 52415 : kern( 2, 2 ) = kp1x * kp1y;
166 52415 : kern( 2, 3 ) = kp1x * kp2y;
167 :
168 52415 : kern( 3, 0 ) = kp2x * km2y;
169 52415 : kern( 3, 1 ) = kp2x * km1y;
170 52415 : kern( 3, 2 ) = kp2x * kp1y;
171 52415 : kern( 3, 3 ) = kp2x * kp2y;
172 52415 : }
173 : };
174 :
175 : /// Typedef for cubicConvolTransform with single precision
176 : /** \ingroup image_transforms
177 : */
178 : typedef cubicConvolTransform<float> cubicConvolTransf;
179 :
180 : /// Typedef for cubicConvolTransform with double precision
181 : /** \ingroup image_transforms
182 : */
183 : typedef cubicConvolTransform<double> cubicConvolTransd;
184 :
185 : /** \ingroup image_transforms
186 : * @{
187 : */
188 :
189 : /// Rotate an image represented as an eigen array
190 : /** Uses the given transformation type to rotate an image.
191 : *
192 : * \tparam transformT specifies the transformation to use [will be resolved by compiler]
193 : * \tparam arrT is the eigen array type of the output [will be resolved by compiler]
194 : * \tparam arrT2 is the eigen array type of the input [will be resolved by compiler]
195 : * \tparam floatT is a floating point type [will be resolved by compiler in most cases]
196 : *
197 : */
198 : template <typename transformT, typename arrT, typename arrT2, typename floatT>
199 3 : void imageRotate( arrT &transim, ///< [out] The rotated image. Must be pre-allocated.
200 : const arrT2 &im, ///< [in] The image to be rotated.
201 : floatT dq, ///< [in] the angle, in radians, by which to rotate in the c.c.w. direction
202 : transformT trans ///< [in] is the transformation to use
203 : )
204 : {
205 : typedef typename transformT::arithT arithT;
206 : arithT cosq, sinq;
207 : arithT x0, y0, x, y;
208 : arithT xcen, ycen;
209 :
210 : int Nrows, Ncols;
211 :
212 : int i0, j0;
213 :
214 3 : const int lbuff = transformT::lbuff;
215 3 : const int width = transformT::width;
216 :
217 3 : cosq = cos( dq );
218 3 : sinq = sin( dq );
219 :
220 3 : Nrows = im.rows();
221 3 : Ncols = im.cols();
222 :
223 3 : transim.resize( Nrows, Ncols );
224 :
225 : // The geometric image center
226 3 : xcen = 0.5 * ( Nrows - 1. );
227 3 : ycen = 0.5 * ( Ncols - 1. );
228 :
229 3 : int xulim = Nrows - width + lbuff; // - 1;
230 3 : int yulim = Ncols - width + lbuff; // - 1;
231 :
232 3 : arithT xc_x_cosq = xcen * cosq;
233 3 : arithT xc_x_sinq = xcen * sinq;
234 3 : arithT yc_x_cosq = ycen * cosq;
235 3 : arithT yc_x_sinq = ycen * sinq;
236 :
237 3 : xc_x_cosq += yc_x_sinq;
238 3 : xc_x_sinq -= yc_x_cosq;
239 :
240 : // clang-format off
241 : #ifdef MXLIB_USE_OMP
242 : #pragma omp parallel private( x0, y0, i0, j0, x, y )
243 : #endif // clang-format on
244 : {
245 : arithT i_x_cosq, i_x_sinq;
246 3 : arrT kern;
247 3 : kern.resize( width, width );
248 :
249 : // clang-format off
250 : #ifdef MXLIB_USE_OMP
251 : #pragma omp for schedule( static, 1 )
252 : #endif // clang-format on
253 54 : for( int i = 0; i < Nrows; ++i )
254 : {
255 51 : i_x_cosq = i * cosq - xc_x_cosq; // + xcen;
256 51 : i_x_sinq = -( i * sinq - xc_x_sinq ); // + ycen;
257 :
258 1302 : for( int j = 0; j < Ncols; ++j )
259 : {
260 : // We are actually doing this rotation matrix:
261 : // x0 = (i-xcen)*cosq + (j-ycen)*sinq;
262 : // y0 = -(i-xcen)*sinq + (j-ycen)*cosq;
263 : // This is the minimum-op representation of the above rotation matrix:
264 1251 : x0 = i_x_cosq + j * sinq;
265 1251 : y0 = i_x_sinq + j * cosq;
266 :
267 : // Get lower left index
268 1251 : i0 = x0 + xcen;
269 1251 : j0 = y0 + ycen;
270 :
271 1251 : if( i0 <= lbuff || i0 >= xulim || j0 <= lbuff || j0 >= yulim )
272 : {
273 450 : transim( i, j ) = 0;
274 450 : continue;
275 : }
276 :
277 : // Get the residual
278 801 : x = x0 + xcen - i0;
279 801 : y = y0 + ycen - j0;
280 :
281 801 : trans( kern, x, y );
282 801 : transim( i, j ) = ( im.block( i0 - lbuff, j0 - lbuff, width, width ) * kern ).sum();
283 : } // for j
284 : } // for i
285 3 : } // #pragma omp parallel
286 :
287 3 : } // void imageRotate(arrT & transim, const arrT2 &im, floatT dq, transformT trans)
288 :
289 : /// Shift an image by whole pixels with (optional) wrapping.
290 : /** The output image can be smaller than the input image, in which case the wrapping (if enabled) still occurs for the
291 : * input image, but only output images worth of pixels are actually shifted. This is useful, for instance, when
292 : * propagating large turbulence phase screens where one only needs a small section at a time.
293 : *
294 : * \tparam outputArrT is the eigen array type of the output [will be resolved by compiler]
295 : * \tparam inputArrT is the eigen array type of the input [will be resolved by compiler]
296 : *
297 : * "[test doc]"
298 : */
299 : template <typename outputArrT, typename inputArrT>
300 13 : void imageShiftWP( outputArrT &out, ///< [out] contains the shifted image. Must be pre-allocated, but can be smaller
301 : ///< than the in array.
302 : inputArrT &in, ///< [in] the image to be shifted.
303 : int dx, ///< [in] the amount to shift in the x direction
304 : int dy, ///< [in] the amount to shift in the y direction
305 : bool wrap = true ///< [in] flag controlling whether or not to wrap around
306 : )
307 : {
308 13 : dx %= in.rows();
309 13 : dy %= in.cols();
310 :
311 13 : int outr = out.rows();
312 13 : int outc = out.cols();
313 13 : int inr = in.rows();
314 13 : int inc = in.cols();
315 :
316 : // clang-format off
317 : #ifdef MXLIB_USE_OMP
318 : #pragma omp parallel
319 : #endif // clang-format on
320 : {
321 : int x, y;
322 :
323 13 : if( wrap )
324 : {
325 : // clang-format off
326 : #ifdef MXLIB_USE_OMP
327 : #pragma omp for
328 : #endif //clang-format on
329 786 : for( int cc = 0; cc < outc; ++cc )
330 : {
331 774 : y = cc - dy;
332 :
333 774 : if( y < 0 )
334 396 : y += inc;
335 378 : else if( y >= inc )
336 0 : y -= inc;
337 :
338 50700 : for( int rr = 0; rr < outr; ++rr )
339 : {
340 49926 : x = rr - dx;
341 :
342 49926 : if( x < 0 )
343 25542 : x += inr;
344 24384 : else if( x >= inr )
345 0 : x -= inr;
346 :
347 49926 : out( rr, cc ) = in( x, y );
348 : }
349 : }
350 : }
351 : else
352 : {
353 : // clang-format off
354 : #ifdef MXLIB_USE_OMP
355 : #pragma omp for
356 : #endif // clang-format on
357 257 : for( int cc = 0; cc < outc; ++cc )
358 : {
359 256 : y = cc - dy;
360 :
361 256 : if( y < 0 || y >= inc )
362 : {
363 257 : for( int rr = 0; rr < outr; ++rr )
364 : {
365 256 : out( rr, cc ) = 0;
366 : }
367 :
368 1 : continue;
369 1 : }
370 :
371 65535 : for( int rr = 0; rr < outr; ++rr )
372 : {
373 65280 : x = rr - dx;
374 :
375 65280 : if( x < 0 || x >= inr )
376 : {
377 255 : out( rr, cc ) = 0;
378 255 : continue;
379 : }
380 :
381 65025 : out( rr, cc ) = in( x, y );
382 : }
383 : }
384 : } // if(wrap)-else
385 : }
386 13 : }
387 :
388 : /// Shift an image by whole pixels, wrapping around, with a scaling image applied to the shifted image.
389 : /** The output image can be smaller than the input image, in which case the wrapping still occurs for the input image,
390 : * but only output images worth of pixels are actually shifted. This is useful, for instance, when propagating large
391 : * turbulence phase screens where one only needs a small section at a time.
392 : *
393 : * The scaling is applied to the output image. The scale image must be the same size as the output image.
394 : *
395 : * \tparam outputArrT is the eigen array type of the output [will be resolved by compiler]
396 : * \tparam inputArrT is the eigen array type of the input [will be resolved by compiler]
397 : * \tparam scaleArrT is the eigen array type of the scale image [will be resolved by compiler]
398 : *
399 : */
400 : template <typename outputArrT, typename inputArrT, typename scaleArrT>
401 : void imageShiftWPScale( outputArrT &out, /**< [out] contains the shifted image. Must be pre-allocated,
402 : but can be smaller than the in array. */
403 : inputArrT &in, ///< [in] the image to be shifted.
404 : scaleArrT &scale, /**< [in] image of scale values applied per-pixel to the output (shifted)
405 : image, same size as out*/
406 : int dx, ///< [in] the amount to shift in the x direction
407 : int dy ///< [in] the amount to shift in the y direction
408 : )
409 : {
410 : dx %= in.rows();
411 : dy %= in.cols();
412 :
413 : int outr = out.rows();
414 : int outc = out.cols();
415 : int inr = in.rows();
416 : int inc = in.cols();
417 :
418 : // clang-format off
419 : #ifdef MXLIB_USE_OMP
420 : // #pragma omp parallel
421 : #endif // clang-format on
422 : {
423 : int x, y;
424 :
425 : // clang-format off
426 : #ifdef MXLIB_USE_OMP
427 : // #pragma omp for
428 : #endif // clang-format on
429 : for( int cc = 0; cc < outc; ++cc )
430 : {
431 : y = cc - dy;
432 :
433 : if( y < 0 )
434 : y += inc;
435 : else if( y >= inc )
436 : y -= inc;
437 :
438 : for( int rr = 0; rr < outr; ++rr )
439 : {
440 : x = rr - dx;
441 :
442 : if( x < 0 )
443 : x += inr;
444 : else if( x >= inr )
445 : x -= inr;
446 :
447 : out( rr, cc ) = in( x, y ) * scale( rr, cc );
448 : }
449 : }
450 : }
451 : }
452 :
453 : /// Shift an image.
454 : /** Uses the given transformation type to shift an image such that objects move by (\p dx,\p dy) pixels.
455 : * The shift is such that an object located at the coordinate
456 : * (\p -dx, \p -dy) from the center of the image will be moved to the center of the image. So to move an object
457 : * located 2 pixels right (dx) and 2 pixels up (dy) from the center to be at the center, use \p dx = -2, \p dy = -2.
458 : *
459 : * Note that this does not treat the edges
460 : * of the image, determined by the buffer width (lbuff) of the kernel and the size of shift. If you wish to
461 : * treat the edges, you must pad the image by at least lbuff+abs(shift) pixels in each direction, and
462 : * implement a strategy (zeros, mirror, wrap) prior to calling this function.
463 : *
464 : * \tparam arrOutT is the Eigen-like array type of the output [will be resolved by compiler]
465 : * \tparam arrInT is the Eigen-like array type of the input [will be resolved by compiler]
466 : * \tparam floatT1 is a floating point type [will be resolved by compiler]
467 : * \tparam floatT2 is a floating point type [will be resolved by compiler]
468 : * \tparam transformT specifies the transformation to use [will be resolved by compiler]
469 : *
470 : * "[test doc]"
471 : */
472 : template <typename arrOutT, typename arrInT, typename floatT1, typename floatT2, typename transformT>
473 6 : void imageShift( arrOutT &transim, ///< [out] Will contain the shifted image. Will be allocated.
474 : const arrInT &im, ///< [in] the image to be shifted.
475 : floatT1 dx, ///< [in] the amount to shift in the x direction
476 : floatT2 dy, ///< [in] the amount to shift in the y direction
477 : transformT trans ///< [in] trans is the transformation to use
478 : )
479 : {
480 : typedef typename transformT::arithT arithT;
481 :
482 : int Nrows, Ncols;
483 :
484 6 : const int lbuff = transformT::lbuff;
485 6 : const int width = transformT::width;
486 :
487 : // If this is a whole pixel, just do that.
488 6 : if( dx == floor( dx ) && dy == floor( dy ) )
489 1 : return imageShiftWP( transim, im, dx, dy, false );
490 :
491 5 : Nrows = im.rows();
492 5 : Ncols = im.cols();
493 :
494 5 : int xulim = Nrows - width + lbuff;
495 5 : int yulim = Ncols - width + lbuff;
496 :
497 5 : transim.resize( Nrows, Ncols );
498 :
499 : #ifdef MXLIB_USE_OMP
500 : #pragma omp parallel
501 : #endif
502 : {
503 : int i0, j0;
504 : // (rx, ry) is fractional residual of shift
505 5 : arithT rx = 1 - ( dx - floor( dx ) );
506 5 : arithT ry = 1 - ( dy - floor( dy ) );
507 :
508 5 : arrOutT kern;
509 5 : kern.resize( width, width );
510 5 : trans( kern, rx, ry );
511 :
512 : #ifdef MXLIB_USE_OMP
513 : #pragma omp for
514 : #endif
515 1285 : for( int i = 0; i < Nrows; ++i )
516 : {
517 : // (i,j) is position in new image
518 : // (i0,j0) is integer position in old image
519 :
520 1280 : i0 = i - dx;
521 :
522 1280 : if( i0 <= lbuff || i0 >= xulim )
523 : {
524 6425 : for( int j = 0; j < Ncols; ++j )
525 : {
526 6400 : transim( i, j ) = 0;
527 : }
528 25 : continue;
529 25 : }
530 :
531 322535 : for( int j = 0; j < Ncols; ++j )
532 : {
533 321280 : j0 = j - dy;
534 :
535 321280 : if( j0 <= lbuff || j0 >= yulim )
536 : {
537 6275 : transim( i, j ) = 0;
538 6275 : continue;
539 : }
540 :
541 315005 : transim( i, j ) = ( im.block( i0 - lbuff, j0 - lbuff, width, width ) * kern ).sum();
542 : } // for j
543 : } // for i
544 5 : } // #pragam omp
545 :
546 : } // imageShift
547 :
548 : /// Magnify an image.
549 : /** Uses the given transformation type to magnify the input image to the size of the output image.
550 : *
551 : * Here we assume that the image center is the mxlib standard:
552 : * \code
553 : x_center = 0.5*(im.rows()-1);
554 : y_center = 0.5*(im.cols()-1);
555 : * \endcode
556 : * Some care is necessary to prevent magnification from shifting the image with respect to this center. The main
557 : result is that the
558 : * magnification factors (which can be different in x and y) are defined thus:
559 : * \code
560 : x_mag = (transim.rows()-1.0) / (im.rows()-1.0);
561 : y_mag = (transim.cols()-1.0) / (im.cols()-1.0);
562 : * \endcode
563 : *
564 : * Example:
565 : * \code
566 : im1.resize(512,512);
567 : //add image to im1
568 : im2.resize(1024,1024);
569 : imageMagnify(im2,im1, cubicConvolTransform<double>());
570 : \endcode
571 : * In this exmple, the image in im1 will be magnified by `1023.0/511.0 = 2.002x` and placed in im2.
572 : *
573 : * This transform function does not handle edges. If treatment of edges is desired, you must pad the input
574 : * image using the desired strategy before calling this function. Note that the padded-size of the input image
575 : * will affect the magnification factor.
576 : *
577 : * \tparam arrOutT is the eigen array type of the output.
578 : * \tparam arrInT is the eigen array type of the input.
579 : * \tparam transformT specifies the transformation to use.
580 : */
581 : template <typename arrOutT, typename arrInT, typename transformT>
582 9 : void imageMagnify( arrOutT &transim, ///< [out] contains the magnified image. Must be pre-allocated.
583 : const arrInT &im, ///< [in] is the image to be magnified.
584 : transformT trans ///< [in] is the transformation to use
585 : )
586 : {
587 : typedef typename transformT::arithT arithT;
588 :
589 : arithT x0, y0, x, y;
590 :
591 : int Nrows, Ncols;
592 :
593 : int i0, j0;
594 :
595 9 : const int lbuff = transformT::lbuff;
596 9 : const int width = transformT::width;
597 :
598 9 : Nrows = transim.rows();
599 9 : Ncols = transim.cols();
600 :
601 9 : int xulim = im.rows() - lbuff - 1;
602 9 : int yulim = im.cols() - lbuff - 1;
603 :
604 9 : arithT x_scale = ( (arithT)im.rows() - 1.0 ) / ( transim.rows() - 1.0 ); // this is 1/x_mag
605 9 : arithT y_scale = ( (arithT)im.cols() - 1.0 ) / ( transim.cols() - 1.0 ); // this is 1/y_mag
606 :
607 9 : arithT xcen = 0.5 * ( (arithT)transim.rows() - 1.0 );
608 9 : arithT ycen = 0.5 * ( (arithT)transim.cols() - 1.0 );
609 :
610 9 : arithT xcen0 = 0.5 * ( (arithT)im.rows() - 1.0 );
611 9 : arithT ycen0 = 0.5 * ( (arithT)im.cols() - 1.0 );
612 :
613 : // #pragma omp parallel private(x0,y0,i0,j0,x,y) num_threads(4)
614 : {
615 9 : arrOutT kern;
616 9 : kern.resize( width, width );
617 :
618 858 : for( int j = 0; j < Ncols; ++j )
619 : {
620 : // (i,j) is position in new image
621 : // (x0,y0) is true position in old image
622 : // (i0,j0) is integer position in old image
623 : // (x, y) is fractional residual of (x0-i0, y0-j0)
624 :
625 849 : y0 = ycen0 + ( j - ycen ) * y_scale;
626 849 : j0 = y0;
627 :
628 849 : if( j0 < lbuff || j0 >= yulim )
629 : {
630 18018 : for( int i = 0; i < Nrows; ++i )
631 : {
632 17829 : transim( i, j ) = 0;
633 : }
634 189 : continue;
635 189 : }
636 :
637 : // #pragma omp for
638 66120 : for( int i = 0; i < Nrows; ++i )
639 : {
640 65460 : x0 = xcen0 + ( i - xcen ) * x_scale;
641 65460 : i0 = x0; // just converting to int
642 :
643 65460 : if( i0 < lbuff || i0 >= xulim )
644 : {
645 13860 : transim( i, j ) = 0;
646 13860 : continue;
647 : }
648 :
649 : // Get the residual
650 51600 : x = x0 - i0;
651 51600 : y = y0 - j0;
652 :
653 51600 : trans( kern, x, y );
654 51600 : transim( i, j ) = ( im.block( i0 - lbuff, j0 - lbuff, width, width ) * kern ).sum();
655 : } // for j
656 : } // for i
657 9 : } // #pragma omp
658 9 : }
659 :
660 : /// Magnify an image with the cubic convolution interpolator.
661 : /** Uses the cubic convolution interpolator to magnify the input image to the size of the output image.
662 : *
663 : * This is a wrapper for imageMagnify with the transform type specified.
664 : *
665 : * \tparam arrOutT is the eigen array type of the output.
666 : * \tparam arrInT is the eigen array type of the input.
667 : *
668 : * \overload
669 : */
670 : template <typename arrOutT, typename arrInT>
671 : void imageMagnify( arrOutT &transim, ///< [out] contains the magnified image. Must be pre-allocated.
672 : const arrInT &im ///< [in] is the image to be magnified.
673 : )
674 : {
675 : return imageMagnify( transim, im, cubicConvolTransform<typename arrInT::Scalar>() );
676 : }
677 :
678 : /// Re-bin an image using the sum, reducing its size while conserving the total flux.
679 : /** Optionally this can be the mean instead of the sum filter, in which case total flux is not conserved.
680 : */
681 : template <typename imageOutT, typename imageInT>
682 : int imageRebinSum( imageOutT &imout, ///< [out] the re-binned image. Must be allocated to size which is an integer
683 : ///< factor smaller than imin.
684 : const imageInT &imin, ///< [in] the image to rebin
685 : bool mean = false ///< [in] if true the output is the mean rather than the sum.
686 : )
687 : {
688 : int rebin = imin.rows() / imout.rows();
689 : if( imin.cols() / imout.cols() != rebin )
690 : return -1;
691 :
692 : int N = 1;
693 : if( mean )
694 : N = rebin * rebin;
695 : for( int i = 0; i < imout.rows(); ++i )
696 : {
697 : for( int j = 0; j < imout.cols(); ++j )
698 : {
699 : imout( i, j ) = imin.block( i * rebin, j * rebin, rebin, rebin ).sum() / N;
700 : }
701 : }
702 :
703 : return 0;
704 : }
705 :
706 : /// Re-bin an image using the sum, reducing its size while conserving the total flux. Records the value and position of
707 : /// the re-binned max pixel.
708 : /** Optionally this can be the mean instead of the sum filter, in which case total flux is not conserved.
709 : *
710 : * \overload
711 : */
712 : template <typename imageOutT, typename imageInT>
713 : int imageRebinSum( imageOutT &imout, ///< [out] the re-binned image. Must be allocated to size which is an integer
714 : ///< factor smaller than imin.
715 : int &xMax, ///< [out] the x-locatioin of the max pixel
716 : int &yMax, ///< [out] the y-locatioin of the max pixel
717 : typename imageOutT::Scalar &pMax, ///< [out] the value of the max pixel
718 : const imageInT &imin, ///< [in] the image to rebin
719 : bool mean = false ///< [in] if true the output is the mean rather than the sum.
720 : )
721 : {
722 : int rebin = imin.rows() / imout.rows();
723 : if( imin.cols() / imout.cols() != rebin )
724 : return -1;
725 :
726 : int N = 1;
727 : if( mean )
728 : N = rebin * rebin;
729 :
730 : xMax = 0;
731 : yMax = 0;
732 : pMax = std::numeric_limits<typename imageOutT::Scalar>::lowest();
733 :
734 : for( int i = 0; i < imout.rows(); ++i )
735 : {
736 : for( int j = 0; j < imout.cols(); ++j )
737 : {
738 : imout( i, j ) = imin.block( i * rebin, j * rebin, rebin, rebin ).sum() / N;
739 : if( imout( i, j ) > pMax )
740 : {
741 : pMax = imout( i, j );
742 : xMax = i;
743 : yMax = j;
744 : }
745 : }
746 : }
747 :
748 : return 0;
749 : }
750 :
751 : /// Re-bin an image using the mean.
752 : /** This is a wrapper for imageRebinSum with `mean=true`.
753 : */
754 : template <typename imageOutT, typename imageInT>
755 : int imageRebinMean( imageOutT &imout, ///< [out] the re-binned image. Must be allocated to size which is an integer
756 : ///< factor smaller than imin.
757 : const imageInT &imin ///< [in] the image to rebin
758 : )
759 : {
760 : return imageRebinSum( imout, imin, true );
761 : }
762 :
763 : /// Re-bin an image using the mean. Records the value and position of the re-binned max pixel.
764 : /** This is a wrapper for imageRebinSum with `mean=true`.
765 : *
766 : * \overload
767 : */
768 : template <typename imageOutT, typename imageInT>
769 : int imageRebinMean( imageOutT &imout, ///< [out] the re-binned image. Must be allocated to size which is an integer
770 : ///< factor smaller than imin.
771 : int &xMax, ///< [out] the x-locatioin of the max pixel
772 : int &yMax, ///< [out] the y-locatioin of the max pixel
773 : typename imageOutT::Scalar &pMax, ///< [out] the value of the max pixel
774 : const imageInT &imin, ///< [in] the image to rebin
775 : bool mean = false ///< [in] if true the output is the mean rather than the sum.
776 : )
777 : {
778 : return imageRebinSum( imout, xMax, yMax, pMax, imin, true );
779 : }
780 :
781 : /// Re-bin an image, takes the mean with a min/max rejection.
782 : /** The mean is calculated after rejecting the minimuma and maximum value.
783 : */
784 : template <typename imageOutT, typename imageInT>
785 : int imageRebinMeanReject( imageOutT &imout, ///< [out] the re-binned image. Must be allocated to size which is an
786 : ///< integer factor smaller than imin.
787 : const imageInT &imin ///< [in] the image to rebin
788 : )
789 : {
790 : int rebin = imin.rows() / imout.rows();
791 : if( imin.cols() / imout.cols() != rebin )
792 : return -1;
793 :
794 : int N = rebin * rebin - 2;
795 :
796 : for( int i = 0; i < imout.rows(); ++i )
797 : {
798 : for( int j = 0; j < imout.cols(); ++j )
799 : {
800 : typename imageOutT::Scalar sum = 0;
801 : typename imageOutT::Scalar max = imin( i * rebin, j * rebin );
802 : typename imageOutT::Scalar min = imin( i * rebin, j * rebin );
803 : for( int k = 0; k < rebin; ++k )
804 : {
805 : for( int l = 0; l < rebin; ++l )
806 : {
807 : sum += imin( i * rebin + k, j * rebin + l );
808 : if( imin( i * rebin + k, j * rebin + l ) > max )
809 : max = imin( i * rebin + k, j * rebin + l );
810 : if( imin( i * rebin + k, j * rebin + l ) < min )
811 : min = imin( i * rebin + k, j * rebin + l );
812 : }
813 : }
814 : imout( i, j ) = ( sum - max - min ) / N; /**/
815 : }
816 : }
817 :
818 : return 0;
819 : }
820 :
821 : /// Re-bin an image, takes the mean with a min/max rejection. Records the value and position of the re-binned max
822 : /// pixel.
823 : /** The mean is calculated after rejecting the minimuma and maximum value.
824 : *
825 : * \overload
826 : */
827 : template <typename imageOutT, typename imageInT>
828 : int imageRebinMeanReject( imageOutT &imout, ///< [out] the re-binned image. Must be allocated to size which is an
829 : ///< integer factor smaller than imin.
830 : int &xMax, ///< [out] the x-locatioin of the max pixel
831 : int &yMax, ///< [out] the y-locatioin of the max pixel
832 : typename imageOutT::Scalar &pMax, ///< [out] the value of the max pixel
833 : const imageInT &imin ///< [in] the image to rebin
834 : )
835 : {
836 : int rebin = imin.rows() / imout.rows();
837 : if( imin.cols() / imout.cols() != rebin )
838 : return -1;
839 :
840 : int N = rebin * rebin - 2;
841 :
842 : xMax = 0;
843 : yMax = 0;
844 : pMax = std::numeric_limits<typename imageOutT::Scalar>::lowest();
845 :
846 : for( int i = 0; i < imout.rows(); ++i )
847 : {
848 : for( int j = 0; j < imout.cols(); ++j )
849 : {
850 : typename imageOutT::Scalar sum = 0;
851 : typename imageOutT::Scalar max = imin( i * rebin, j * rebin );
852 : typename imageOutT::Scalar min = imin( i * rebin, j * rebin );
853 : for( int k = 0; k < rebin; ++k )
854 : {
855 : for( int l = 0; l < rebin; ++l )
856 : {
857 : sum += imin( i * rebin + k, j * rebin + l );
858 : if( imin( i * rebin + k, j * rebin + l ) > max )
859 : max = imin( i * rebin + k, j * rebin + l );
860 : if( imin( i * rebin + k, j * rebin + l ) < min )
861 : min = imin( i * rebin + k, j * rebin + l );
862 : }
863 : }
864 : imout( i, j ) = ( sum - max - min ) / N;
865 :
866 : if( imout( i, j ) > pMax )
867 : {
868 : pMax = imout( i, j );
869 : xMax = i;
870 : yMax = j;
871 : }
872 : }
873 : }
874 :
875 : return 0;
876 : }
877 :
878 : /// Down-sample an image, reducing its size while conserving the total flux.
879 : /** If the old size is an integer multiple of the new size, this is just a re-bin. If not an integer multiple,
880 : * the image is interpolated after performing the closest re-bin, and then re-normalized to conserve flux.
881 : *
882 : * \todo Allow selection of interpolator, providing a default version.
883 : */
884 : template <typename imageOutT, typename imageInT>
885 : void imageDownSample( imageOutT &imout, const imageInT &imin )
886 : {
887 : typedef typename imageOutT::Scalar Scalar;
888 :
889 : // Record this for normalization later
890 : Scalar inputTotal = fabs( imin.sum() );
891 :
892 : // As a first step, rebin to nearest whole pixel factor which is larger than the desired output size
893 : int closestRebin = imin.rows() / imout.rows(); //, imin.cols()/imout.cols() );
894 :
895 : float sample = ( (float)imin.rows() ) / closestRebin;
896 :
897 : while( sample != floor( sample ) )
898 : {
899 : --closestRebin;
900 : if( closestRebin == 1 )
901 : break;
902 : sample = ( (float)imin.rows() ) / closestRebin;
903 : }
904 :
905 : // Eigen::Array<Scalar, Eigen::Dynamic, Eigen::Dynamic> temp;
906 : eigenImage<Scalar> temp;
907 : temp.resize( imin.rows() / closestRebin, imin.cols() / closestRebin );
908 :
909 : for( int i = 0; i < temp.rows(); ++i )
910 : {
911 : for( int j = 0; j < temp.cols(); ++j )
912 : {
913 : temp( i, j ) = imin.block( i * closestRebin, j * closestRebin, closestRebin, closestRebin ).sum();
914 : }
915 : }
916 :
917 : // If the output image is now the requested size return.
918 : if( temp.rows() == imout.rows() && temp.cols() == imout.cols() )
919 : {
920 : imout = temp;
921 : return;
922 : }
923 : // Otherwise, re-sample using bilinear interpolation.
924 : typedef bilinearTransform<Scalar> transformT;
925 :
926 : transformT trans;
927 : // Eigen::Array<Scalar, -1,-1> kern;
928 : eigenImage<Scalar> kern;
929 :
930 : const int lbuff = transformT::lbuff;
931 : const int width = transformT::width;
932 :
933 : for( int i = 0; i < imout.rows(); ++i )
934 : {
935 : for( int j = 0; j < imout.cols(); ++j )
936 : {
937 : double x = ( (double)i / imout.rows() ) * temp.rows();
938 : double y = ( (double)j / imout.cols() ) * temp.cols();
939 :
940 : trans( kern, x - floor( x ), y - floor( y ) );
941 :
942 : imout( i, j ) = ( temp.block( floor( x ) - lbuff, floor( y ) - lbuff, width, width ) * kern ).sum();
943 : }
944 : }
945 :
946 : // Normalize
947 : Scalar outputTotal = fabs( imout.sum() );
948 : imout *= inputTotal / outputTotal;
949 : }
950 :
951 : ///@}
952 :
953 : } // namespace improc
954 : } // namespace mx
955 :
956 : #endif // improc_imageTransforms_hpp
|
