Logo Search packages:      
Sourcecode: octave-image version File versions  Download package

poly2mask.m

## Copyright (C) 2004 Josep Mones i Teixidor
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {@var{BW} = } poly2mask (@var{x},@var{y},@var{m},@var{n})
## Convert a polygon to a region mask.
##
## BW=poly2mask(x,y,m,n) converts a polygon, specified by a list of
## vertices in @var{x} and @var{y} and returns in a @var{m}-by-@var{n}
## logical mask @var{BW} the filled polygon. Region inside the polygon
## is set to 1, values outside the shape are set to 0.
##
## @var{x} and @var{y} should always represent a closed polygon, first
## and last points should be coincident. If they are not poly2mask will
## close it for you. If @var{x} or @var{y} are fractional they are
## nearest integer.
##
## If all the polygon or part of it falls outside the masking area
## (1:m,1:n), it is discarded or clipped.
##
## This function uses scan-line polygon filling algorithm as described
## in http://www.cs.rit.edu/~icss571/filling/ with some minor
## modifications: capability of clipping and scan order, which can
## affect the results of the algorithm (algorithm is described not to
## reach ymax, xmax border when filling to avoid enlarging shapes). In
## this function we scan the image backwards (we begin at ymax and end
## at ymin), and we don't reach ymin, xmin, which we believe should be
## compatible with MATLAB.
## @end deftypefn


## TODO: check how to create a logical BW without any conversion

## Author:  Josep Mones i Teixidor <jmones@puntbarra.com>

function BW = poly2mask (x, y, m, n)
  if (nargin != 4)
    print_usage ();
  endif

  ## check x and y
  x = round (x (:).');
  y = round (y (:).');
  if (length (x) < 3)
    error ("poly2mask: polygon must have at least 3 vertices.");
  endif
  if (length (x) != length (y))
    error ("poly2mask: length of x doesn't match length of y.");
  endif

  ## create output matrix
  BW = false (m, n);

  ## close polygon if needed
  if ((x (1) != x (length (x))) || (y (1) != y (length (y))))
    x = horzcat (x, x (1));
    y = horzcat (y, y (1));
  endif

  ## build global edge table
  ex = ; ## x values for each edge
  ey = ; ## y values for each edge
  idx = (ey (1, :) != ey (2, :));                 ## eliminate horizontal edges
  ex = ex (:, idx);
  ey = ey (:, idx);
  eminy = min (ey);                               ## minimum y for each edge
  emaxy = max (ey);                               ## maximum y for each edge
  t = (ey == [eminy; eminy]);                     ## values associated to miny
  exminy = ex (:) (t);                            ## x values associated to min y
  exmaxy = ex (:) (!t);                           ## x values associated to max y
  emaxy = emaxy.';                                ## we want them vertical now...
  eminy = eminy.';
  m_inv = (exmaxy - exminy)./(emaxy - eminy);     ## calculate inverse slope
  ge = [emaxy, eminy, exmaxy, m_inv];             ## build global edge table
  ge = sortrows (ge, [1, 3]);                     ## sort on eminy and exminy

  ## we add an extra dummy edge at the end just to avoid checking
  ## while indexing it
  ge = [-Inf, -Inf, -Inf, -Inf; ge];

  ## initial parity is even (0)
  parity = 0;

  ## init scan line set to bottom line
  sl = ge (size (ge, 1), 1);

  ## init active edge table
  ## we use a loop because the table is sorted and edge list could be
  ## huge
  ae = ;
  gei = size (ge, 1);
  while (sl == ge (gei, 1))
    ae = [ge(gei, 2:4); ae];
    gei -= 1;
  endwhile

  ## calc minimum y to draw
  miny = min (y);
  if (miny < 1)
    miny = 1;
  endif

  while (sl >= miny)
    ## check vert clipping
    if (sl <= m)
      ## draw current scan line
      ## we have to round because 1/m is fractional
      ie = round (reshape (ae (:, 2), 2, size (ae, 1)/2));

      ## this discards left border of image (this differs from version at
      ## http://www.cs.rit.edu/~icss571/filling/ which discards right
      ## border) but keeps an exception when the point is a vertex.
      ie (1, :) += (ie (1, :) != ie (2, :));

      ## we'll clip too, just in case m,n is not big enough
      ie (1, (ie (1, :) < 1)) = 1;
      ie (2, (ie (2, :) > n)) = n;

      ## we eliminate segments outside window
      ie = ie (:, (ie (1, :) <= n));
      ie = ie (:, (ie (2, :) >= 1));
      for i = 1:columns (ie)
      BW (sl, ie (1, i):ie (2, i)) = true;
      endfor
    endif

    ## decrement scan line
    sl -= 1;

    ## eliminate edges that eymax==sl
    ## this discards ymin border of image (this differs from version at
    ## http://www.cs.rit.edu/~icss571/filling/ which discards ymax).
    ae = ae ((ae (:, 1) != sl), :);

    ## update x (x1=x0-1/m)
    ae (:, 2) -= ae (:, 3);

    ## update ae with new values
    while (sl == ge (gei, 1))
      ae = vertcat (ae, ge (gei, 2:4));
      gei -= 1;
    endwhile

    ## order the edges in ae by x value
    if (rows (ae) > 0)
      ae = sortrows (ae, 2);
    endif
  endwhile
endfunction

## This should create a filled octagon
%!demo
%! s = [0:pi/4:2*pi];
%! x = cos (s) * 90 + 101;
%! y = sin (s) * 90 + 101;
%! bw = poly2mask(x, y, 200, 200);
%! imshow (bw);

## This should create a 5-vertex star
%!demo
%! s = ;
%! s = s ([1, 3, 5, 2, 4, 6]);
%! x = cos (s) * 90 + 101;
%! y = sin (s) * 90 + 101;
%! bw = poly2mask (x, y, 200, 200);
%! imshow (bw);

%!# Convex polygons

%!shared xs, ys, Rs, xt, yt, Rt
%! xs=;
%! ys=;
%! Rs=zeros(16,14);
%! Rs(5:12,4:10)=1;
%! Rs=logical(Rs);
%! xt=;
%! yt=;
%! Rt=;
%! Rt=logical(Rt);

%!assert(poly2mask(xs,ys,16,14),Rs);          # rectangle
%!assert(poly2mask(xs,ys,8,7),Rs(1:8,1:7));   # clipped
%!assert(poly2mask(xs-7,ys-8,8,7),Rs(9:16,8:14)); # more clipping

%!assert(poly2mask(xt,yt,5,7),Rt);            # triangle
%!assert(poly2mask(xt,yt,3,3),Rt(1:3,1:3));   # clipped


%!# Concave polygons

%!test
%! x=;
%! y=;
%! R=zeros(16,14);
%! R(5:12,4:5)=1;
%! R(5:8,6:8)=1;
%! R(5:11,9:10)=1;
%! R=logical(R);
%! assert(poly2mask(x,y,16,14), R);

%!# Complex polygons
%!test
%! x=;
%! y=;
%! R=;
%! R=logical(R);
%! assert(poly2mask(x,y,5,6), R);





%
% $Log$
% Revision 1.3  2007/03/23 16:14:37  adb014
% Update the FSF address
%
% Revision 1.2  2007/01/04 23:37:54  hauberg
% Minor changes in help text
%
% Revision 1.1  2006/08/20 12:59:35  hauberg
% Changed the structure to match the package system
%
% Revision 1.5  2004/09/07 14:47:50  pkienzle
% Avoid segfaults on pre-2.1.58 octave.  Invisible whitespace changes.
%
% Revision 1.4  2004/09/03 17:12:36  jmones
% Uses uint8 to save some temporal memory (suggested by David Bateman)
%
% Revision 1.3  2004/09/03 13:32:07  jmones
% Work with logical arrays from BW creation
%
% Revision 1.2  2004/08/11 17:39:51  jmones
% Algorithm url in docs corrected.
%
% Revision 1.1  2004/08/11 17:34:11  jmones
% poly2mask added: creates filled polygon bw mask
%
%

Generated by  Doxygen 1.6.0   Back to index