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