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