Loading [MathJax]/extensions/tex2jax.js
mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Modules Pages
imageMasks.hpp
Go to the documentation of this file.
1/** \file imageMasks.hpp
2 * \brief Declares and defines functions to work with image masks
3 * \ingroup image_processing_files
4 * \author Jared R. Males (jaredmales@gmail.com)
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_imageMasks_hpp
28#define improc_imageMasks_hpp
29
30#include "../math/constants.hpp"
31#include "../math/geo.hpp"
32
33#include "eigenImage.hpp"
34#include "imageTransforms.hpp"
35
36namespace mx
37{
38
39namespace improc
40{
41
42/// Fills in the cells of an Eigen 2D Array with their radius from the center
43/** \ingroup image_masks
44 *
45 * \tparam eigenT is an Eigen-like 2D array type
46 */
47template <class eigenT>
48void radiusImage( eigenT &m, ///< [out] the allocated radius array, will be filled in with radius values.
49 typename eigenT::Scalar xc, ///< [in] the x center
50 typename eigenT::Scalar yc, ///< [in] the y center
51 typename eigenT::Scalar scale = 1 ///< [in] [optional] a scaling to apply to each value (default = 1)
52)
53{
54 typedef typename eigenT::Scalar arithT;
55
56 arithT f_x, f_y;
57
58 size_t dim1 = m.rows();
59 size_t dim2 = m.cols();
60
61 for( size_t i = 0; i < dim1; i++ )
62 {
63 f_x = ( i - xc ) * ( i - xc );
64
65 for( size_t j = 0; j < dim2; j++ )
66 {
67 f_y = ( j - yc ) * ( j - yc );
68
69 m( i, j ) = sqrt( f_x + f_y ) * scale;
70 }
71 }
72}
73
74/// Fills in the cells of Eigen-like 2D Array with their radius from the canonical center
75/** \ingroup image_masks
76 *
77 * The center is @f$ (x_c, y_c) = (0.5*(dim_1-1), 0.5*(dim_2 -1)) @f$.
78 *
79 * \tparam eigenT is an Eigen-like 2D array type
80 */
81template <class eigenT>
82void radiusImage( eigenT &m, ///< [out] the allocated radius array, will be filled in with radius values.
83 typename eigenT::Scalar scale = 1 ///< [in] [optional] a scaling to apply to each value (default = 1)
84)
85{
86 typedef typename eigenT::Scalar arithT;
87
88 arithT xc, yc;
89
90 xc = 0.5 * ( m.rows() - 1 );
91 yc = 0.5 * ( m.cols() - 1 );
92
93 radiusImage( m, xc, yc, scale );
94}
95
96/// Fills in the cells of an Eigen-like 2D Array with their angle relative to the center
97/** \ingroup image_masks
98 *
99 * \tparam angleT is the angle type, either radiansT<realT> or degreesT<realT>. Note that realT sets the type for all
100 * arithmetic. \tparam eigenT is an Eigen-like 2D array type
101 */
102template <class angleT, class eigenT>
103void angleImage( eigenT &m, ///< [out] the allocated angle array. Will be filled in with angle values.
104 typename angleT::realT xc, ///< [in] the x center
105 typename angleT::realT yc ///< [in] the y center
106)
107{
108 typedef typename angleT::realT realT;
109
110 for( size_t j = 0; j < m.cols(); j++ )
111 {
112 realT f_y = ( static_cast<realT>( j ) - yc );
113
114 for( size_t i = 0; i < m.rows(); i++ )
115 {
116 realT f_x = ( static_cast<realT>( i ) - xc );
117
118 m( i, j ) = fmod( atan2( f_y, f_x ) + math::two_pi<realT>(), math::two_pi<realT>() ) * angleT::scale;
119 }
120 }
121}
122
123/// Fills in the cells of Eigen 2D Array with their angle relative the canonical center
124/** \ingroup image_masks
125 *
126 * The center is @f$ (x_c, y_c) = (0.5*(dim_1-1), 0.5*(dim_2 -1)) @f$.
127 *
128 * \tparam angleT is the angle type, either radiansT<realT> or degreesT<realT>. Note that realT sets the type for all
129 * arithmetic. \tparam eigenT is an Eigen-like 2D array type
130 */
131template <class angleT, class eigenT>
132void angleImage( eigenT &m /** < [out] the allocated angle array. Will be filled in with angle values. */ )
133{
134 typedef typename angleT::realT realT;
135
136 realT xc = 0.5 * ( m.rows() - 1 );
137 realT yc = 0.5 * ( m.cols() - 1 );
138
139 angleImage<angleT>( m, xc, yc );
140}
141
142/// Fills in the cells of Eigen-like arrays with their radius amd angle relative to the center
143/** \ingroup image_masks
144 *
145 * \tparam angleT is the angle type, either radiansT<realT> or degreesT<realT>. Note that realT sets the type for all
146 * arithmetic. \tparam eigenT1 is an Eigen-like 2D array type. Should be resolved by compiler. \tparam eigenT2 is an
147 * Eigen-like 2D array type. Should be resolved by compiler.
148 *
149 */
150template <class angleT, class eigenT1, class eigenT2>
152 eigenT1 &rIm, ///< [out] the allocated radius array, will be filled in with radius values.
153 eigenT2 &qIm, ///< [out] the angle array, will be re-sized to match rIm. Will be filled in with angle values.
154 typename angleT::realT xc, ///< [in] the x center
155 typename angleT::realT yc, ///< [in] the y center
156 typename angleT::realT rscale = 1 ///< [in] [optional] a scaling to apply to each radius value. Default is 1.0.
157)
158{
159 typedef typename angleT::realT realT;
160
161 qIm.resize( rIm.rows(), rIm.cols() );
162
163 for( int cc = 0; cc < rIm.cols(); ++cc )
164 {
165 realT f_y = ( ( static_cast<realT>( cc ) ) - yc );
166
167 for( int rr = 0; rr < rIm.rows(); ++rr )
168 {
169 realT f_x = ( ( static_cast<realT>( rr ) ) - xc );
170
171 rIm( rr, cc ) = std::sqrt( f_x * f_x + f_y * f_y ) * rscale;
172 qIm( rr, cc ) = fmod( atan2( f_y, f_x ) + math::two_pi<realT>(), math::two_pi<realT>() ) * angleT::scale;
173 }
174 }
175}
176
177/// Get the vector indices of an annular region in an image
178/**
179 * \ingroup image_masks
180 *
181 * \tparam angleT is the angle type, either radiansT<realT> or degreesT<realT>. Note that realT sets the type for all
182 * arithmetic. \tparam eigenT1 is an Eigen-like 2D array type. Should be resolved by compiler. \tparam eigenT2 is an
183 * Eigen-like 2D array type. Should be resolved by compiler. \tparam eigenT3 is an Eigen-like 2D array type. Should be
184 * resolved by compiler.
185 *
186 * \returns a vector containing the 1D indices of the region defined by the input parameters
187 */
188template <typename angleT, typename eigenT1, typename eigenT2, typename eigenT3 = eigenT1>
189std::vector<size_t> annulusIndices( const eigenT1 &rIm, /**< [in] a radius image of the type produced by
190 \ref radiusImage */
191 const eigenT2 &qIm, /**< [in] an angle image of the type produce by
192 \ref angleImage */
193 typename angleT::realT xcen, /**< [in] the x center of the image */
194 typename angleT::realT ycen, /**< [in] the y center of the image */
195 typename angleT::realT min_r, /**< [in] the minimum radius of the region */
196 typename angleT::realT max_r, /**< [in] the maximum radius of the region */
197 typename angleT::realT min_q, /**< [in] the minimum angle of the region. */
198 typename angleT::realT max_q, /**< [in] the maximum angle of the region. */
199 eigenT3 *mask = 0 /**< [in] [optional] pointer to a mask image, only
200 pixels of value 1 are included
201 in the indices. */
202)
203{
204
205 std::vector<size_t> idx;
206
207 int min_x = -max_r;
208 int max_x = max_r;
209 int min_y = -max_r;
210 int max_y = max_r;
211
212 int x0 = xcen + min_x;
213 if( x0 < 0 )
214 x0 = 0;
215
216 int x1 = xcen + max_x + 1;
217 if( x1 > rIm.rows() )
218 x1 = rIm.rows();
219
220 int y0 = ycen + min_y;
221 if( y0 < 0 )
222 y0 = 0;
223
224 int y1 = ycen + max_y + 1;
225 if( y1 > rIm.cols() )
226 y1 = rIm.cols();
227
228 // Normalize the angles by mod
229 min_q = math::angleMod<angleT>( min_q );
230 // If max_q is exactly 360/2pi (to within a tol of 100*eps) we don't mod, other mod.
231 if( fabs( max_q - angleT::full ) > 100 * std::numeric_limits<typename angleT::realT>::epsilon() )
232 {
233 max_q = math::angleMod<angleT>( max_q );
234 }
235
236 // Find the mid-point to enable sectors wider than 180.
237 typename angleT::realT mid_q;
238 if( min_q <= max_q )
239 {
240 mid_q = 0.5 * ( min_q + max_q );
241 }
242 else
243 {
244 mid_q = 0.5 * ( ( min_q - angleT::full ) + max_q );
245 }
246
247 for( int j = y0; j < y1; ++j )
248 {
249 for( int i = x0; i < x1; ++i )
250 {
251 if( rIm( i, j ) < min_r )
252 {
253 continue;
254 }
255 if( rIm( i, j ) >= max_r )
256 {
257 continue;
258 }
259 if( !( ( math::angleDiff<angleT>( min_q, qIm( i, j ) ) >= 0 &&
260 math::angleDiff<angleT>( qIm( i, j ), mid_q ) > 0 ) ||
261 ( math::angleDiff<angleT>( mid_q, qIm( i, j ) ) >= 0 &&
262 math::angleDiff<angleT>( qIm( i, j ), max_q ) > 0 ) ) )
263 {
264 continue;
265 }
266
267 if( mask )
268 {
269 if( ( *mask )( i, j ) == 0 )
270 {
271 continue;
272 }
273 }
274
275 idx.push_back( j * rIm.rows() + i );
276 }
277 }
278
279 return idx;
280}
281
282/// Get the coordinates of the bounding rectangle of an annulus.
283/**
284 * \ingroup image_masks
285 *
286 * \tparam angleT is the angle type, either radiansT<realT> or degreesT<realT>. Note that realT sets the type for all
287 * arithmetic.
288 *
289 */
290template <typename angleT>
291void annulusBoundingRect( int &x0, ///< [out] The lower left x-coordinate of the bounding rect.
292 int &y0, ///< [out] The lower left y-coordinate of the bounding rect.
293 int &x1, ///< [out] The upper right x-coordinate of the bounding rect.
294 int &y1, ///< [out] The upper right y-coordinate of the bounding rect.
295 typename angleT::realT xcen, ///< [in] the x center of the image
296 typename angleT::realT ycen, ///< [in] the y center of the image
297 typename angleT::realT min_r, ///< [in] the minimum radius of the region
298 typename angleT::realT max_r, ///< [in] the maximum radius of the region
299 typename angleT::realT min_q, ///< [in] the minimum angle of the region.
300 typename angleT::realT max_q ///< [in] the maximum angle of the region.
301)
302{
303 typedef typename angleT::realT realT;
304
305 // Get the corners
306 realT x00 = xcen + min_r * cos( min_q * angleT::radians );
307 realT y00 = ycen + min_r * sin( min_q * angleT::radians );
308 realT x01 = xcen + max_r * cos( min_q * angleT::radians );
309 realT y01 = ycen + max_r * sin( min_q * angleT::radians );
310
311 realT x10 = xcen + min_r * cos( max_q * angleT::radians );
312 realT y10 = ycen + min_r * sin( max_q * angleT::radians );
313 realT x11 = xcen + max_r * cos( max_q * angleT::radians );
314 realT y11 = ycen + max_r * sin( max_q * angleT::radians );
315
316 // vertex of min_r is probably not necessary, but might as well
317 realT x20 = xcen + min_r * cos( math::angleMean<angleT>( { min_q, max_q } ) * angleT::radians );
318 realT y20 = ycen + min_r * sin( math::angleMean<angleT>( { min_q, max_q } ) * angleT::radians );
319
320 // vertex of max_r
321 realT x21 = xcen + max_r * cos( math::angleMean<angleT>( { min_q, max_q } ) * angleT::radians );
322 realT y21 = ycen + max_r * sin( math::angleMean<angleT>( { min_q, max_q } ) * angleT::radians );
323
324 x0 = std::ceil( std::min( { x00, x01, x10, x11, x20, x21 } ) );
325 y0 = std::ceil( std::min( { y00, y01, y10, y11, y20, y21 } ) );
326 x1 = std::ceil( std::max( { x00, x01, x10, x11, x20, x21 } ) );
327 y1 = std::ceil( std::max( { y00, y01, y10, y11, y20, y21 } ) );
328}
329
330/// Reflect vector indices across the given center pixel.
331/** This assumes column major order.
332 * \returns the vector of reflected indices
333 */
334template <typename realT>
335std::vector<size_t> reflectImageIndices( const std::vector<size_t> &idxi, ///< [in] the vector indices to reflect
336 int w, ///< [in] the image width
337 int h, ///< [in] the image height
338 realT xc, ///< [in] the image center x coordinate
339 realT yc ///< [in] the image center y coordinate
340)
341{
342 std::vector<size_t> idxr;
343 idxr.reserve( idxi.size() );
344
345 int x0, y0, x1, y1;
346 for( size_t n = 0; n < idxi.size(); ++n )
347 {
348 int y0 = idxi[n] / w;
349 int x0 = idxi[n] - y0 * w;
350
351 reflectImageCoords( x1, y1, x0, y0, xc, yc );
352
353 if( x1 < w && y1 < h )
354 idxr.push_back( y1 * w + x1 );
355 }
356
357 return idxr;
358}
359
360template <typename sizeT>
361void rectangleIndices( std::vector<sizeT> &idx, sizeT rows, sizeT cols, sizeT xmin, sizeT xmax, sizeT ymin, sizeT ymax )
362{
363
364 // if(xmin < 0) xmin = 0;
365 if( xmax > rows - 1 )
366 xmax = rows - 1;
367
368 // if(ymin < 0) ymin = 0;
369 if( ymax > cols - 1 )
370 ymax = cols - 1;
371
372 idx.reserve( ( xmax - xmin + 1 ) * ( ymax - ymin + 1 ) );
373
374 for( sizeT i = xmin; i <= xmax; ++i )
375 {
376 for( sizeT j = ymin; j <= ymax; ++j )
377 {
378 idx.push_back( j * rows + i );
379 }
380 }
381}
382
383template <class eigenT>
384void rectangleIndices( std::vector<size_t> &idx, eigenT &mask, size_t xmin, size_t xmax, size_t ymin, size_t ymax )
385{
386 rectangleIndices<size_t>( idx, (size_t)mask.rows(), (size_t)mask.cols(), xmin, xmax, ymin, ymax );
387}
388
389/// Apply a mask to an image
390/** The pixels indicated by the indices are set to a value.
391 *
392 * \ingroup image_masks
393 *
394 * \tparam eigenT is an Eigen-like 2D array type
395 *
396 */
397template <class eigenT>
398void applyMask( eigenT &maskedIm, ///< [out] the image to mask (will be modified)
399 const std::vector<size_t> &idx, ///< [in] the indices of the pixels to mask
400 typename eigenT::Scalar maskval ///< [in] the mask value.
401)
402{
403 for( size_t i = 0; i < idx.size(); ++i )
404 {
405 maskedIm( idx[i] ) = maskval;
406 }
407}
408
409/// Mask a circle in an image.
410/** The circle is describe by its center coordinates and radius. Any value can be set for the mask.
411 * Pixels outside the masked circle are not altered.
412 *
413 * \tparam arrayT is an Eigen-like type.
414 *
415 * \ingroup image_masks
416 */
417template <class arrayT>
418void maskCircle( arrayT &m, /**< [in.out] the image to be masked, is modified. */
419 typename arrayT::Scalar xcen, /**< [in] the x coordinate of the center of the circle */
420 typename arrayT::Scalar ycen, /**< [in] the y coordinate of the center of the circle */
421 typename arrayT::Scalar rad, /**< [in] the radius of the circle */
422 typename arrayT::Scalar val, /**< [in] the mask value. */
423 typename arrayT::Scalar pixbuf = 0.5 /**< [in] [optional] buffer for radius comparison.
424 Default is 0.5 pixels.*/
425)
426{
427 size_t l0 = m.rows();
428 size_t l1 = m.cols();
429
430 typename arrayT::Scalar rr;
431
432 for( size_t c = 0; c < m.cols(); c++ )
433 {
434 for( size_t r = 0; r < m.rows(); r++ )
435 {
436 rr = sqrt( pow( r - xcen, 2 ) + pow( c - ycen, 2 ) );
437
438 if( rr <= rad + pixbuf )
439 m( r, c ) = val;
440 }
441 }
442}
443
444/// Mask a circle in an image at the standard center.
445/** The circle is centered at 0.5*(rows-1) and 0.5*(cols-1), and described by its radius. Any value can be set for the
446 * mask. Pixels outside the masked circle are not altered.
447 *
448 * \tparam arrayT is an Eigen-like type.
449 *
450 * \ingroup image_masks
451 */
452template <class arrayT>
454 arrayT &m, ///< [in.out] the image to be masked, is modified.
455 typename arrayT::Scalar rad, ///< [in] the radius of the circle
456 typename arrayT::Scalar val, ///< [in] the mask value.
457 typename arrayT::Scalar pixbuf = 0.5 ///< [in] [optional] buffer for radius comparison. Default is 0.5 pixels.
458)
459{
460 return maskCircle( m, 0.5 * ( m.rows() - 1.0 ), 0.5 * ( m.cols() - 1.0 ), rad, val, pixbuf );
461}
462
463/// Mask an ellipse in an image.
464/** The ellipse is describe by its center coordinates and x and y direction radii (the semi-major and -minor axes,
465 * in either order). Any value can be set for the mask,cwith 0 being the default.
466 *
467 * \tparam arrayT is an Eigen-like type.
468 *
469 * \ingroup image_masks
470 */
471template <class arrayT>
472void maskEllipse( arrayT &m, ///< [in.out] the image to be masked, is modified.
473 typename arrayT::Scalar xcen, ///< [in] the x coordinate of the center of the ellipse
474 typename arrayT::Scalar ycen, ///< [in] the y coordinate of the center of the ellipse
475 typename arrayT::Scalar xrad, ///< [in] the x radius of the ellipse
476 typename arrayT::Scalar yrad, ///< [in] the y radius of the ellipse
477 typename arrayT::Scalar ang, ///< [in] the c.c.w. angle to rotate the ellipse by
478 typename arrayT::Scalar val = 0, ///< [in] [optional] the mask value. Default is 0.
479 typename arrayT::Scalar pixbuf = 0.5 /**< [in] [optional] buffer for radius comparison.
480 Default is 0.5 pixels.*/
481)
482{
483 typedef typename arrayT::Scalar realT;
484
485 size_t l0 = m.rows();
486 size_t l1 = m.cols();
487
488 realT r;
489 realT x;
490 realT y;
491 realT xe, ye;
492 realT rad;
493
494 realT cq = cos( ang );
495 realT sq = sin( ang );
496
497 realT xr, yr;
498
499 for( size_t i = 0; i < l0; i++ )
500 {
501 x = i - xcen;
502 for( size_t j = 0; j < l1; j++ )
503 {
504 y = j - ycen;
505
506 if( x == 0 && y == 0 )
507 {
508 m( i, j ) = val;
509 }
510 else
511 {
512 xr = x * cq - y * sq;
513 yr = x * sq + y * cq;
514
515 // Coordinate on the ellipse where the line to the point intersects
516 xe = ( pow( xrad * yrad, 2 ) / ( pow( yrad, 2 ) + pow( xrad * yr / xr, 2 ) ) );
517 ye = ( pow( yrad, 2 ) - xe * pow( yrad / xrad, 2 ) );
518
519 rad = sqrt( xe + ye );
520
521 r = sqrt( pow( xr, 2 ) + pow( yr, 2 ) );
522
523 if( r <= rad + pixbuf )
524 {
525 m( i, j ) = val;
526 }
527 }
528 }
529 }
530}
531
532/// Mask a wedge in an image.
533/** The wedge is describe by its center coordinate, central angle, and angular half-width. Any value can be set for the
534 * mask. Pixels outside the masked circle are not altered.
535 *
536 * \tparam arrayT is an Eigen-like type.
537 *
538 * \ingroup image_masks
539 */
540template <class arrayT>
541void maskWedge( arrayT &m, ///< [in.out] the image to be masked, is modified.
542 typename arrayT::Scalar xcen, ///< [in] the x coordinate of the center of the circle
543 typename arrayT::Scalar ycen, ///< [in] the y coordinate of the center of the circle
544 typename arrayT::Scalar angCen, ///< [in] the central angle of the wedge, in degrees
545 typename arrayT::Scalar angHW, ///< [in] the angular half-wdith of the wedge, in degrees
546 typename arrayT::Scalar val = 0 ///< [in] [optional] the mask value. Default is 0.
547)
548{
549 size_t l0 = m.rows();
550 size_t l1 = m.cols();
551
552 typedef typename arrayT::Scalar angleT;
553
554 angleT dang;
555
556 for( size_t c = 0; c < m.cols(); c++ )
557 {
558 for( size_t r = 0; r < m.rows(); r++ )
559 {
560
561 dang = math::angleDiff<math::degreesT<angleT>>( math::rtod( atan2( c - ycen, r - xcen ) ), angCen );
562
563 if( dang > -angHW && dang <= angHW )
564 {
565 m( r, c ) = val;
566 }
567 }
568 }
569}
570
571/// Draw a thin (1-pixel) line from one point to another.
572/** Calculates the line connecting the two points and sets the pixels along that line to the
573 * supplied value.
574 *
575 * \tparam realT a real floating point type
576 *
577 * \ingroup image_masks
578 */
579template <typename realT>
580int drawLine( eigenImage<realT> &mask, ///< [in/out] [pre-allocated] The array in which to draw the line.
581 realT x0, ///< [in] The x coordinate of the first point
582 realT y0, ///< [in] The y coordinate of the first point
583 realT x1, ///< [in] The x coordinate of the second point
584 realT y1, ///< [in] The y coordinate of the second point
585 realT val ///< [in] The value to set on the pixels in the line
586)
587{
588 if( fabs( x1 - x0 ) <= 1 && fabs( y1 - y0 ) <= 1 )
589 {
590 // If it is a single pixel, plot at the rounded average point.
591 realT xa = 0.5 * ( x0 + x1 ) + 0.5;
592 realT ya = 0.5 * ( y0 + y1 ) + 0.5;
593 if( xa >= 0 && xa < mask.rows() && ya >= 0 && ya < mask.cols() )
594 mask( (int)xa, (int)ya ) = val;
595 return 0;
596 }
597
598 realT _x0, _x1, _y0, _y1;
599
600 realT length = sqrt( pow( x1 - x0, 2 ) + pow( y1 - y0, 2 ) );
601
602 // If X doesn't change enough, we try to do it as a function of y.
603 if( fabs( x1 - x0 ) <= 1 )
604 {
605 if( y0 < y1 )
606 {
607 _x0 = x0;
608 _x1 = x1;
609 _y0 = y0;
610 _y1 = y1;
611 }
612 else
613 {
614 _x0 = x1;
615 _x1 = x0;
616 _y0 = y1;
617 _y1 = y0;
618 }
619
620 realT m = ( _x1 - _x0 ) / ( _y1 - _y0 );
621 realT b = _x0 - m * _y0;
622
623 realT dy = ( _y1 - _y0 ) / ( 2 * length + 1 );
624 for( realT y = _y0; y <= _y1; y += dy )
625 {
626 realT x = m * y + b;
627
628 if( x + 0.5 >= 0 && x + 0.5 < mask.rows() && y + 0.5 >= 0 && y + 0.5 < mask.cols() )
629 mask( (int)( x + 0.5 ), (int)( y + 0.5 ) ) = val;
630 }
631
632 return 0;
633 }
634
635 // Ordinarily draw in x.
636 if( x0 < x1 )
637 {
638 _x0 = x0;
639 _x1 = x1;
640 _y0 = y0;
641 _y1 = y1;
642 }
643 else
644 {
645 _x0 = x1;
646 _x1 = x0;
647 _y0 = y1;
648 _y1 = y0;
649 }
650
651 realT m = ( _y1 - _y0 ) / ( _x1 - _x0 );
652 realT b = _y0 - m * _x0;
653
654 realT dx = ( _x1 - _x0 ) / ( 2 * length + 1 );
655 for( realT x = _x0; x <= _x1; x += dx )
656 {
657 realT y = m * x + b;
658
659 if( x + 0.5 >= 0 && x + 0.5 < mask.rows() && y + 0.5 >= 0 && y + 0.5 < mask.cols() )
660 mask( (int)( x + 0.5 ), (int)( y + 0.5 ) ) = val;
661 }
662}
663
664/// Draw a thick line from one point to another.
665/** Calculates the line connecting the two points and sets the pixels along that line to the
666 * supplied value. Makes the line thick by drawing lines perpindicular to each point with
667 * length equal to the specified width.
668 *
669 * \tparam realT a real floating point type
670 *
671 * \ingroup image_masks
672 */
673template <typename realT>
674int drawLine( eigenImage<realT> &mask, ///< [in.out] [pre-allocated] The array in which to draw the line.
675 realT x0, ///< [in] The x coordinate of the first point
676 realT y0, ///< [in] The y coordinate of the first point
677 realT x1, ///< [in] The x coordinate of the second point
678 realT y1, ///< [in] The y coordinate of the second point
679 realT width, ///< [in] The desired full-width of the line
680 realT val ///< [in] The value to set on the pixels in the line
681)
682{
683 // You can't do this.
684 if( width <= 0 )
685 return -1;
686
687 // If no reason to try this hard, don't.
688 if( width < 1 )
689 return drawLine( mask, x0, y0, x1, y1, val );
690
691 realT _x0, _x1, _y0, _y1;
692
693 // If X doesn't change enough, we try to do it as a function of y.
694 if( fabs( x1 - x0 ) <= 1 )
695 {
696 if( y0 < y1 )
697 {
698 _x0 = x0;
699 _x1 = x1;
700 _y0 = y0;
701 _y1 = y1;
702 }
703 else
704 {
705 _x0 = x1;
706 _x1 = x0;
707 _y0 = y1;
708 _y1 = y0;
709 }
710
711 realT m = ( _x1 - _x0 ) / ( _y1 - _y0 );
712 realT b = _x0 - m * _y0;
713
714 realT dy = ( _y1 - _y0 ) / ( mask.cols() + 1 );
715 for( realT y = _y0; y <= _y1; y += dy )
716 {
717 realT x = m * y + b;
718
719 realT xs, ys, xe, ye; // The start and end point of the perpindicular line.
720
721 realT q1 = atan( m ) + 0.5 * math::pi<realT>();
722 realT cq = cos( q1 );
723 realT sq = sin( q1 );
724
725 xs = x + 0.5 * width * cq;
726 ys = y + 0.5 * width * sq;
727
728 xe = x - 0.5 * width * cq;
729 ye = y - 0.5 * width * sq;
730
731 drawLine( mask, xs, ys, xe, ye, val );
732 }
733
734 return 0;
735 }
736
737 // Ordinarily draw in x.
738 if( x0 < x1 )
739 {
740 _x0 = x0;
741 _x1 = x1;
742 _y0 = y0;
743 _y1 = y1;
744 }
745 else
746 {
747 _x0 = x1;
748 _x1 = x0;
749 _y0 = y1;
750 _y1 = y0;
751 }
752
753 realT m = ( _y1 - _y0 ) / ( _x1 - _x0 );
754 realT b = _y0 - m * _x0;
755
756 realT dx = ( _x1 - _x0 ) / ( mask.rows() + 1 );
757 for( realT x = _x0; x <= _x1; x += dx )
758 {
759 realT y = m * x + b;
760
761 realT xs, ys, xe, ye; // The start and end point of the perpindicular line.
762
763 realT q1 = atan( m ) + 0.5 * math::pi<realT>();
764 realT cq = cos( q1 );
765 realT sq = sin( q1 );
766
767 xs = x + 0.5 * width * cq;
768 ys = y + 0.5 * width * sq;
769
770 xe = x - 0.5 * width * cq;
771 ye = y - 0.5 * width * sq;
772
773 drawLine( mask, xs, ys, xe, ye, val );
774 }
775}
776
777/// Populate a mask based on a typical CCD bleeding pattern.
778/** Masks a circle for saturated pixels, as well as a horizontal bar for bleeding (in both directions if needed).
779 *
780 * \returns 0 on success, -1 otherwise.
781 *
782 * \tparam imT is an Eigen-like 2D array type.
783 *
784 * \ingroup image_masks
785 */
786template <typename imT>
787int ccdBleedMask( imT &im, ///< [out] the mask, on output will be 1/0. Must be allocated prior to call.
788 typename imT::Scalar x0, ///< [in] the x-coordinate of the center of the mask
789 typename imT::Scalar y0, ///< [in] the y-coordinate of the center of the mask
790 typename imT::Scalar rad, ///< [in] the radius of the central circle.
791 typename imT::Scalar height, ///< [in] the half-height of the saturation bar.
792 typename imT::Scalar lwidth, ///< [in] the length of the saturation bar to the left.
793 typename imT::Scalar rwidth ///< [in] the length of the saturation bar to the right.
794)
795{
796 typename imT::Scalar x, y;
797
798 for( int i = 0; i < im.rows(); ++i )
799 {
800 x = i;
801 for( int j = 0; j < im.cols(); ++j )
802 {
803 y = j;
804
805 if( sqrt( pow( x - x0, 2 ) + pow( y - y0, 2 ) ) <= rad )
806 {
807 im( i, j ) = 0;
808 }
809 else if( fabs( y - y0 ) <= height && x >= x0 - lwidth && x <= x0 + rwidth )
810 {
811 im( i, j ) = 0;
812 }
813 else
814 {
815 im( i, j ) = 1;
816 }
817 }
818 }
819
820 return 0;
821}
822
823/// Cut out a region of an image specified by an index-mask.
824/** The output will be a row-image containing the pixel values.
825 *
826 * \tparam imageTout is the output image Eigen-like type.
827 * \tparam imageTin is the input image Eigen-like type.
828 *
829 * \ingroup image_masks
830 */
831template <typename imageTout, typename imageTin>
832void cutImageRegion( imageTout &imout, ///< [out] a row-image containing the pixel values specified by the indices.
833 const imageTin &imin, ///< [in] a 2-D image
834 const std::vector<size_t> &idx, ///< [in] the linear indices of the pixel values
835 bool resize = true ///< [in] [optional] flag controls whether imout is resized. It is if true.
836)
837{
838 if( resize )
839 {
840 imout.resize( idx.size(), 1 );
841 }
842
843 // #pragma omp parallel for schedule(static, 1)
844 for( size_t i = 0; i < idx.size(); ++i )
845 {
846 imout( i ) = imin( idx[i] );
847 }
848}
849
850/// Insert a region of an image specified by an index-mask.
851/** Both the input and the output are row-images containing the pixel values.
852 *
853 * \tparam imageTout is the output image Eigen-like type.
854 * \tparam imageTin is the input image Eigen-like type.
855 *
856 * \ingroup image_masks
857 */
858template <typename imageTout, typename imageTin>
859void insertImageRegion( imageTout imout, /**< [out] a row-image into which the pixel values
860 specified by the indices are inserted. */
861 const imageTin &imin, /**< [in] a row-image containing the pixel values,
862 same size as idx*/
863 const std::vector<size_t> &idx /**< [in] the linear indices of the pixel values */
864)
865{
866 // #pragma omp parallel for schedule(static, 1)
867 for( size_t i = 0; i < idx.size(); ++i )
868 {
869 imout( idx[i], 0 ) = imin( i, 0 );
870 }
871}
872
873/// Rotate a binary mask
874/** Sets edge pixels to 0 or 1 depending on the interpolation, being below/above 0.5.
875 */
876template <typename imageT, typename transformT = cubicConvolTransform<typename imageT::Scalar>>
877void rotateMask( imageT &rotMask, imageT &mask, typename imageT::Scalar angle )
878{
879 imageRotate( rotMask, mask, angle, transformT() );
880
881 for( int jj = 0; jj < rotMask.cols(); ++jj )
882 {
883 for( int ii = 0; ii < rotMask.rows(); ++ii )
884 {
885 if( rotMask( ii, jj ) < 0.5 )
886 rotMask( ii, jj ) = 0;
887 else
888 rotMask( ii, jj ) = 1;
889 }
890 }
891}
892
893} // namespace improc
894} // namespace mx
895
896#endif // improc_imageMasks_hpp
Tools for using the eigen library for image processing.
Eigen::Array< scalarT, -1, -1 > eigenImage
Definition of the eigenImage type, which is an alias for Eigen::Array.
constexpr floatT six_fifths()
Return 6/5 in the specified precision.
realT rtod(realT q)
Convert from radians to degrees.
Definition geo.hpp:153
void radiusImage(eigenT &m, typename eigenT::Scalar xc, typename eigenT::Scalar yc, typename eigenT::Scalar scale=1)
Fills in the cells of an Eigen 2D Array with their radius from the center.
void maskCircle(arrayT &m, typename arrayT::Scalar xcen, typename arrayT::Scalar ycen, typename arrayT::Scalar rad, typename arrayT::Scalar val, typename arrayT::Scalar pixbuf=0.5)
Mask a circle in an image.
void cutImageRegion(imageTout &imout, const imageTin &imin, const std::vector< size_t > &idx, bool resize=true)
Cut out a region of an image specified by an index-mask.
void maskEllipse(arrayT &m, typename arrayT::Scalar xcen, typename arrayT::Scalar ycen, typename arrayT::Scalar xrad, typename arrayT::Scalar yrad, typename arrayT::Scalar ang, typename arrayT::Scalar val=0, typename arrayT::Scalar pixbuf=0.5)
Mask an ellipse in an image.
int drawLine(eigenImage< realT > &mask, realT x0, realT y0, realT x1, realT y1, realT val)
Draw a thin (1-pixel) line from one point to another.
void applyMask(eigenT &maskedIm, const std::vector< size_t > &idx, typename eigenT::Scalar maskval)
Apply a mask to an image.
int ccdBleedMask(imT &im, typename imT::Scalar x0, typename imT::Scalar y0, typename imT::Scalar rad, typename imT::Scalar height, typename imT::Scalar lwidth, typename imT::Scalar rwidth)
Populate a mask based on a typical CCD bleeding pattern.
void insertImageRegion(imageTout imout, const imageTin &imin, const std::vector< size_t > &idx)
Insert a region of an image specified by an index-mask.
void angleImage(eigenT &m, typename angleT::realT xc, typename angleT::realT yc)
Fills in the cells of an Eigen-like 2D Array with their angle relative to the center.
void maskWedge(arrayT &m, typename arrayT::Scalar xcen, typename arrayT::Scalar ycen, typename arrayT::Scalar angCen, typename arrayT::Scalar angHW, typename arrayT::Scalar val=0)
Mask a wedge in an image.
std::vector< size_t > annulusIndices(const eigenT1 &rIm, const eigenT2 &qIm, typename angleT::realT xcen, typename angleT::realT ycen, typename angleT::realT min_r, typename angleT::realT max_r, typename angleT::realT min_q, typename angleT::realT max_q, eigenT3 *mask=0)
Get the vector indices of an annular region in an image.
void annulusBoundingRect(int &x0, int &y0, int &x1, int &y1, typename angleT::realT xcen, typename angleT::realT ycen, typename angleT::realT min_r, typename angleT::realT max_r, typename angleT::realT min_q, typename angleT::realT max_q)
Get the coordinates of the bounding rectangle of an annulus.
void radAngImage(eigenT1 &rIm, eigenT2 &qIm, typename angleT::realT xc, typename angleT::realT yc, typename angleT::realT rscale=1)
Fills in the cells of Eigen-like arrays with their radius amd angle relative to the center.
void imageRotate(arrT &transim, const arrT2 &im, floatT dq, transformT trans)
Rotate an image represented as an eigen array.
int reflectImageCoords(int &x1, int &y1, int x0, int y0, realT xc, realT yc)
Reflect pixel coordinates across the given center pixel.
std::vector< size_t > reflectImageIndices(const std::vector< size_t > &idxi, int w, int h, realT xc, realT yc)
Reflect vector indices across the given center pixel.
void rotateMask(imageT &rotMask, imageT &mask, typename imageT::Scalar angle)
Rotate a binary mask.
Image interpolation and transformation.
The mxlib c++ namespace.
Definition mxError.hpp:106