27#ifndef improc_imageMasks_hpp
28#define improc_imageMasks_hpp
30#include "../math/constants.hpp"
31#include "../math/geo.hpp"
47template <
class eigenT>
49 typename eigenT::Scalar xc,
50 typename eigenT::Scalar yc,
51 typename eigenT::Scalar scale = 1
54 typedef typename eigenT::Scalar arithT;
58 size_t dim1 = m.rows();
59 size_t dim2 = m.cols();
61 for(
size_t i = 0; i < dim1; i++ )
63 f_x = ( i - xc ) * ( i - xc );
65 for(
size_t j = 0; j < dim2; j++ )
67 f_y = ( j - yc ) * ( j - yc );
69 m( i, j ) = sqrt( f_x + f_y ) * scale;
81template <
class eigenT>
83 typename eigenT::Scalar scale = 1
86 typedef typename eigenT::Scalar arithT;
90 xc = 0.5 * ( m.rows() - 1 );
91 yc = 0.5 * ( m.cols() - 1 );
102template <
class angleT,
class eigenT>
104 typename angleT::realT xc,
105 typename angleT::realT yc
108 typedef typename angleT::realT realT;
110 for(
size_t j = 0; j < m.cols(); j++ )
112 realT f_y = (
static_cast<realT
>( j ) - yc );
114 for(
size_t i = 0; i < m.rows(); i++ )
116 realT f_x = (
static_cast<realT
>( i ) - xc );
131template <
class angleT,
class eigenT>
134 typedef typename angleT::realT realT;
136 realT xc = 0.5 * ( m.rows() - 1 );
137 realT yc = 0.5 * ( m.cols() - 1 );
139 angleImage<angleT>( m, xc, yc );
150template <
class angleT,
class eigenT1,
class eigenT2>
154 typename angleT::realT xc,
155 typename angleT::realT yc,
156 typename angleT::realT rscale = 1
159 typedef typename angleT::realT realT;
161 qIm.resize( rIm.rows(), rIm.cols() );
163 for(
int cc = 0; cc < rIm.cols(); ++cc )
165 realT f_y = ( (
static_cast<realT
>( cc ) ) - yc );
167 for(
int rr = 0; rr < rIm.rows(); ++rr )
169 realT f_x = ( (
static_cast<realT
>( rr ) ) - xc );
171 rIm( rr, cc ) = std::sqrt( f_x * f_x + f_y * f_y ) * rscale;
188template <
typename angleT,
typename eigenT1,
typename eigenT2,
typename eigenT3 = eigenT1>
193 typename angleT::realT xcen,
194 typename angleT::realT ycen,
195 typename angleT::realT min_r,
196 typename angleT::realT max_r,
197 typename angleT::realT min_q,
198 typename angleT::realT max_q,
205 std::vector<size_t> idx;
212 int x0 = xcen + min_x;
216 int x1 = xcen + max_x + 1;
217 if( x1 > rIm.rows() )
220 int y0 = ycen + min_y;
224 int y1 = ycen + max_y + 1;
225 if( y1 > rIm.cols() )
231 if( fabs( max_q - angleT::full ) > 100 * std::numeric_limits<typename angleT::realT>::epsilon() )
237 typename angleT::realT mid_q;
240 mid_q = 0.5 * ( min_q + max_q );
244 mid_q = 0.5 * ( ( min_q - angleT::full ) + max_q );
247 for(
int j = y0; j < y1; ++j )
249 for(
int i = x0; i < x1; ++i )
251 if( rIm( i, j ) < min_r )
255 if( rIm( i, j ) >= max_r )
269 if( ( *mask )( i, j ) == 0 )
275 idx.push_back( j * rIm.rows() + i );
290template <
typename angleT>
295 typename angleT::realT xcen,
296 typename angleT::realT ycen,
297 typename angleT::realT min_r,
298 typename angleT::realT max_r,
299 typename angleT::realT min_q,
300 typename angleT::realT max_q
303 typedef typename angleT::realT realT;
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 );
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 );
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 } ) );
334template <
typename realT>
342 std::vector<size_t> idxr;
343 idxr.reserve( idxi.size() );
346 for(
size_t n = 0; n < idxi.size(); ++n )
348 int y0 = idxi[n] / w;
349 int x0 = idxi[n] - y0 * w;
353 if( x1 < w && y1 < h )
354 idxr.push_back( y1 * w + x1 );
360template <
typename sizeT>
361void rectangleIndices( std::vector<sizeT> &idx, sizeT rows, sizeT cols, sizeT xmin, sizeT xmax, sizeT ymin, sizeT ymax )
365 if( xmax > rows - 1 )
369 if( ymax > cols - 1 )
372 idx.reserve( ( xmax - xmin + 1 ) * ( ymax - ymin + 1 ) );
374 for( sizeT i = xmin; i <= xmax; ++i )
376 for( sizeT j = ymin; j <= ymax; ++j )
378 idx.push_back( j * rows + i );
383template <
class eigenT>
384void rectangleIndices( std::vector<size_t> &idx, eigenT &mask,
size_t xmin,
size_t xmax,
size_t ymin,
size_t ymax )
386 rectangleIndices<size_t>( idx, (
size_t)mask.rows(), (
size_t)mask.cols(), xmin, xmax, ymin, ymax );
397template <
class eigenT>
399 const std::vector<size_t> &idx,
400 typename eigenT::Scalar maskval
403 for(
size_t i = 0; i < idx.size(); ++i )
405 maskedIm( idx[i] ) = maskval;
417template <
class arrayT>
419 typename arrayT::Scalar xcen,
420 typename arrayT::Scalar ycen,
421 typename arrayT::Scalar rad,
422 typename arrayT::Scalar val,
423 typename arrayT::Scalar pixbuf = 0.5
427 size_t l0 = m.rows();
428 size_t l1 = m.cols();
430 typename arrayT::Scalar rr;
432 for(
size_t c = 0; c < m.cols(); c++ )
434 for(
size_t r = 0; r < m.rows(); r++ )
436 rr = sqrt( pow( r - xcen, 2 ) + pow( c - ycen, 2 ) );
438 if( rr <= rad + pixbuf )
452template <
class arrayT>
455 typename arrayT::Scalar rad,
456 typename arrayT::Scalar val,
457 typename arrayT::Scalar pixbuf = 0.5
460 return maskCircle( m, 0.5 * ( m.rows() - 1.0 ), 0.5 * ( m.cols() - 1.0 ), rad, val, pixbuf );
471template <
class arrayT>
473 typename arrayT::Scalar xcen,
474 typename arrayT::Scalar ycen,
475 typename arrayT::Scalar xrad,
476 typename arrayT::Scalar yrad,
477 typename arrayT::Scalar ang,
478 typename arrayT::Scalar val = 0,
479 typename arrayT::Scalar pixbuf = 0.5
483 typedef typename arrayT::Scalar realT;
485 size_t l0 = m.rows();
486 size_t l1 = m.cols();
494 realT cq = cos( ang );
495 realT sq = sin( ang );
499 for(
size_t i = 0; i < l0; i++ )
502 for(
size_t j = 0; j < l1; j++ )
506 if( x == 0 && y == 0 )
512 xr = x * cq - y * sq;
513 yr = x * sq + y * cq;
516 xe = ( pow( xrad * yrad, 2 ) / ( pow( yrad, 2 ) + pow( xrad * yr / xr, 2 ) ) );
517 ye = ( pow( yrad, 2 ) - xe * pow( yrad / xrad, 2 ) );
519 rad = sqrt( xe + ye );
521 r = sqrt( pow( xr, 2 ) + pow( yr, 2 ) );
523 if( r <= rad + pixbuf )
540template <
class arrayT>
542 typename arrayT::Scalar xcen,
543 typename arrayT::Scalar ycen,
544 typename arrayT::Scalar angCen,
545 typename arrayT::Scalar angHW,
546 typename arrayT::Scalar val = 0
549 size_t l0 = m.rows();
550 size_t l1 = m.cols();
552 typedef typename arrayT::Scalar angleT;
556 for(
size_t c = 0; c < m.cols(); c++ )
558 for(
size_t r = 0; r < m.rows(); r++ )
563 if( dang > -angHW && dang <= angHW )
579template <
typename realT>
588 if( fabs( x1 - x0 ) <= 1 && fabs( y1 - y0 ) <= 1 )
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;
598 realT _x0, _x1, _y0, _y1;
600 realT length = sqrt( pow( x1 - x0, 2 ) + pow( y1 - y0, 2 ) );
603 if( fabs( x1 - x0 ) <= 1 )
620 realT m = ( _x1 - _x0 ) / ( _y1 - _y0 );
621 realT b = _x0 - m * _y0;
623 realT dy = ( _y1 - _y0 ) / ( 2 * length + 1 );
624 for( realT y = _y0; y <= _y1; y += dy )
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;
651 realT m = ( _y1 - _y0 ) / ( _x1 - _x0 );
652 realT b = _y0 - m * _x0;
654 realT dx = ( _x1 - _x0 ) / ( 2 * length + 1 );
655 for( realT x = _x0; x <= _x1; x += dx )
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;
673template <
typename realT>
689 return drawLine( mask, x0, y0, x1, y1, val );
691 realT _x0, _x1, _y0, _y1;
694 if( fabs( x1 - x0 ) <= 1 )
711 realT m = ( _x1 - _x0 ) / ( _y1 - _y0 );
712 realT b = _x0 - m * _y0;
714 realT dy = ( _y1 - _y0 ) / ( mask.cols() + 1 );
715 for( realT y = _y0; y <= _y1; y += dy )
719 realT xs, ys, xe, ye;
722 realT cq = cos( q1 );
723 realT sq = sin( q1 );
725 xs = x + 0.5 * width * cq;
726 ys = y + 0.5 * width * sq;
728 xe = x - 0.5 * width * cq;
729 ye = y - 0.5 * width * sq;
731 drawLine( mask, xs, ys, xe, ye, val );
753 realT m = ( _y1 - _y0 ) / ( _x1 - _x0 );
754 realT b = _y0 - m * _x0;
756 realT dx = ( _x1 - _x0 ) / ( mask.rows() + 1 );
757 for( realT x = _x0; x <= _x1; x += dx )
761 realT xs, ys, xe, ye;
764 realT cq = cos( q1 );
765 realT sq = sin( q1 );
767 xs = x + 0.5 * width * cq;
768 ys = y + 0.5 * width * sq;
770 xe = x - 0.5 * width * cq;
771 ye = y - 0.5 * width * sq;
773 drawLine( mask, xs, ys, xe, ye, val );
786template <
typename imT>
788 typename imT::Scalar x0,
789 typename imT::Scalar y0,
790 typename imT::Scalar rad,
791 typename imT::Scalar height,
792 typename imT::Scalar lwidth,
793 typename imT::Scalar rwidth
796 typename imT::Scalar x, y;
798 for(
int i = 0; i < im.rows(); ++i )
801 for(
int j = 0; j < im.cols(); ++j )
805 if( sqrt( pow( x - x0, 2 ) + pow( y - y0, 2 ) ) <= rad )
809 else if( fabs( y - y0 ) <= height && x >= x0 - lwidth && x <= x0 + rwidth )
831template <
typename imageTout,
typename imageTin>
833 const imageTin &imin,
834 const std::vector<size_t> &idx,
840 imout.resize( idx.size(), 1 );
844 for(
size_t i = 0; i < idx.size(); ++i )
846 imout( i ) = imin( idx[i] );
858template <
typename imageTout,
typename imageTin>
861 const imageTin &imin,
863 const std::vector<size_t> &idx
867 for(
size_t i = 0; i < idx.size(); ++i )
869 imout( idx[i], 0 ) = imin( i, 0 );
876template <
typename imageT,
typename transformT = cubicConvolTransform<
typename imageT::Scalar>>
877void rotateMask( imageT &rotMask, imageT &mask,
typename imageT::Scalar angle )
881 for(
int jj = 0; jj < rotMask.cols(); ++jj )
883 for(
int ii = 0; ii < rotMask.rows(); ++ii )
885 if( rotMask( ii, jj ) < 0.5 )
886 rotMask( ii, jj ) = 0;
888 rotMask( ii, jj ) = 1;
877void rotateMask( imageT &rotMask, imageT &mask,
typename imageT::Scalar angle ) {
…}
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.
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.
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.