octave_packages/image-1.0.15/poly2mask.m

   1 ## Copyright (C) 2004 Josep Mones i Teixidor
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{BW} = } poly2mask (@var{x},@var{y},@var{m},@var{n})
  18 ## Convert a polygon to a region mask.
  19 ##
  20 ## BW=poly2mask(x,y,m,n) converts a polygon, specified by a list of
  21 ## vertices in @var{x} and @var{y} and returns in a @var{m}-by-@var{n}
  22 ## logical mask @var{BW} the filled polygon. Region inside the polygon
  23 ## is set to 1, values outside the shape are set to 0.
  24 ##
  25 ## @var{x} and @var{y} should always represent a closed polygon, first
  26 ## and last points should be coincident. If they are not poly2mask will
  27 ## close it for you. If @var{x} or @var{y} are fractional they are
  28 ## nearest integer.
  29 ##
  30 ## If all the polygon or part of it falls outside the masking area
  31 ## (1:m,1:n), it is discarded or clipped.
  32 ##
  33 ## This function uses scan-line polygon filling algorithm as described
  34 ## in http://www.cs.rit.edu/~icss571/filling/ with some minor
  35 ## modifications: capability of clipping and scan order, which can
  36 ## affect the results of the algorithm (algorithm is described not to
  37 ## reach ymax, xmax border when filling to avoid enlarging shapes). In
  38 ## this function we scan the image backwards (we begin at ymax and end
  39 ## at ymin), and we don't reach ymin, xmin, which we believe should be
  40 ## compatible with MATLAB.
  41 ## @end deftypefn
  42
  43
  44 ## TODO: check how to create a logical BW without any conversion
  45
  46 ## Author:  Josep Mones i Teixidor <jmones@puntbarra.com>
  47
  48 function BW = poly2mask (x, y, m, n)
  49   if (nargin != 4)
  50     print_usage ();
  51   endif
  52
  53   ## check x and y
  54   x = round (x (:).');
  55   y = round (y (:).');
  56   if (length (x) < 3)
  57     error ("poly2mask: polygon must have at least 3 vertices.");
  58   endif
  59   if (length (x) != length (y))
  60     error ("poly2mask: length of x doesn't match length of y.");
  61   endif
  62
  63   ## create output matrix
  64   BW = false (m, n);
  65
  66   ## close polygon if needed
  67   if ((x (1) != x (length (x))) || (y (1) != y (length (y))))
  68     x = horzcat (x, x (1));
  69     y = horzcat (y, y (1));
  70   endif
  71
  72   ## build global edge table
  73   ex = [x(1:length (x) - 1); x(1, 2:length (x))]; ## x values for each edge
  74   ey = [y(1:length (y) - 1); y(1, 2:length (y))]; ## y values for each edge
  75   idx = (ey (1, :) != ey (2, :));                 ## eliminate horizontal edges
  76   ex = ex (:, idx);
  77   ey = ey (:, idx);
  78   eminy = min (ey);                               ## minimum y for each edge
  79   emaxy = max (ey);                               ## maximum y for each edge
  80   t = (ey == [eminy; eminy]);                     ## values associated to miny
  81   exminy = ex (:) (t);                            ## x values associated to min y
  82   exmaxy = ex (:) (!t);                           ## x values associated to max y
  83   emaxy = emaxy.';                                ## we want them vertical now...
  84   eminy = eminy.';
  85   m_inv = (exmaxy - exminy)./(emaxy - eminy);     ## calculate inverse slope
  86   ge = [emaxy, eminy, exmaxy, m_inv];             ## build global edge table
  87   ge = sortrows (ge, [1, 3]);                     ## sort on eminy and exminy
  88
  89   ## we add an extra dummy edge at the end just to avoid checking
  90   ## while indexing it
  91   ge = [-Inf, -Inf, -Inf, -Inf; ge];
  92
  93   ## initial parity is even (0)
  94   parity = 0;
  95
  96   ## init scan line set to bottom line
  97   sl = ge (size (ge, 1), 1);
  98
  99   ## init active edge table
 100   ## we use a loop because the table is sorted and edge list could be
 101   ## huge
 102   ae = [];
 103   gei = size (ge, 1);
 104   while (sl == ge (gei, 1))
 105     ae = [ge(gei, 2:4); ae];
 106     gei -= 1;
 107   endwhile
 108
 109   ## calc minimum y to draw
 110   miny = min (y);
 111   if (miny < 1)
 112     miny = 1;
 113   endif
 114
 115   while (sl >= miny)
 116     ## check vert clipping
 117     if (sl <= m)
 118       ## draw current scan line
 119       ## we have to round because 1/m is fractional
 120       ie = round (reshape (ae (:, 2), 2, size (ae, 1)/2));
 121
 122       ## this discards left border of image (this differs from version at
 123       ## http://www.cs.rit.edu/~icss571/filling/ which discards right
 124       ## border) but keeps an exception when the point is a vertex.
 125       ie (1, :) += (ie (1, :) != ie (2, :));
 126
 127       ## we'll clip too, just in case m,n is not big enough
 128       ie (1, (ie (1, :) < 1)) = 1;
 129       ie (2, (ie (2, :) > n)) = n;
 130
 131       ## we eliminate segments outside window
 132       ie = ie (:, (ie (1, :) <= n));
 133       ie = ie (:, (ie (2, :) >= 1));
 134       for i = 1:columns (ie)
 135         BW (sl, ie (1, i):ie (2, i)) = true;
 136       endfor
 137     endif
 138
 139     ## decrement scan line
 140     sl -= 1;
 141
 142     ## eliminate edges that eymax==sl
 143     ## this discards ymin border of image (this differs from version at
 144     ## http://www.cs.rit.edu/~icss571/filling/ which discards ymax).
 145     ae = ae ((ae (:, 1) != sl), :);
 146
 147     ## update x (x1=x0-1/m)
 148     ae (:, 2) -= ae (:, 3);
 149
 150     ## update ae with new values
 151     while (sl == ge (gei, 1))
 152       ae = vertcat (ae, ge (gei, 2:4));
 153       gei -= 1;
 154     endwhile
 155
 156     ## order the edges in ae by x value
 157     if (rows (ae) > 0)
 158       ae = sortrows (ae, 2);
 159     endif
 160   endwhile
 161 endfunction
 162
 163 ## This should create a filled octagon
 164 %!demo
 165 %! s = [0:pi/4:2*pi];
 166 %! x = cos (s) * 90 + 101;
 167 %! y = sin (s) * 90 + 101;
 168 %! bw = poly2mask(x, y, 200, 200);
 169 %! imshow (bw);
 170
 171 ## This should create a 5-vertex star
 172 %!demo
 173 %! s = [0:2*pi/5:pi*4];
 174 %! s = s ([1, 3, 5, 2, 4, 6]);
 175 %! x = cos (s) * 90 + 101;
 176 %! y = sin (s) * 90 + 101;
 177 %! bw = poly2mask (x, y, 200, 200);
 178 %! imshow (bw);
 179
 180 %!# Convex polygons
 181
 182 %!shared xs, ys, Rs, xt, yt, Rt
 183 %! xs=[3,3,10,10];
 184 %! ys=[4,12,12,4];
 185 %! Rs=zeros(16,14);
 186 %! Rs(5:12,4:10)=1;
 187 %! Rs=logical(Rs);
 188 %! xt=[1,4,7];
 189 %! yt=[1,4,1];
 190 %! Rt=[0,0,0,0,0,0,0;
 191 %!     0,0,1,1,1,1,0;
 192 %!     0,0,0,1,1,0,0;
 193 %!     0,0,0,1,0,0,0;
 194 %!     0,0,0,0,0,0,0];
 195 %! Rt=logical(Rt);
 196
 197 %!assert(poly2mask(xs,ys,16,14),Rs);          # rectangle
 198 %!assert(poly2mask(xs,ys,8,7),Rs(1:8,1:7));   # clipped
 199 %!assert(poly2mask(xs-7,ys-8,8,7),Rs(9:16,8:14)); # more clipping
 200
 201 %!assert(poly2mask(xt,yt,5,7),Rt);            # triangle
 202 %!assert(poly2mask(xt,yt,3,3),Rt(1:3,1:3));   # clipped
 203
 204
 205 %!# Concave polygons
 206
 207 %!test
 208 %! x=[3,3,5,5,8,8,10,10];
 209 %! y=[4,12,12,8,8,11,11,4];
 210 %! R=zeros(16,14);
 211 %! R(5:12,4:5)=1;
 212 %! R(5:8,6:8)=1;
 213 %! R(5:11,9:10)=1;
 214 %! R=logical(R);
 215 %! assert(poly2mask(x,y,16,14), R);
 216
 217 %!# Complex polygons
 218 %!test
 219 %! x=[1,5,1,5];
 220 %! y=[1,1,4,4];
 221 %! R=[0,0,0,0,0,0;
 222 %!    0,0,1,1,0,0;
 223 %!    0,0,1,1,0,0;
 224 %!    0,1,1,1,1,0;
 225 %!    0,0,0,0,0,0];
 226 %! R=logical(R);
 227 %! assert(poly2mask(x,y,5,6), R);
 228
 229
 230
 231
 232
 233 %
 234 % $Log$
 235 % Revision 1.3  2007/03/23 16:14:37  adb014
 236 % Update the FSF address
 237 %
 238 % Revision 1.2  2007/01/04 23:37:54  hauberg
 239 % Minor changes in help text
 240 %
 241 % Revision 1.1  2006/08/20 12:59:35  hauberg
 242 % Changed the structure to match the package system
 243 %
 244 % Revision 1.5  2004/09/07 14:47:50  pkienzle
 245 % Avoid segfaults on pre-2.1.58 octave.  Invisible whitespace changes.
 246 %
 247 % Revision 1.4  2004/09/03 17:12:36  jmones
 248 % Uses uint8 to save some temporal memory (suggested by David Bateman)
 249 %
 250 % Revision 1.3  2004/09/03 13:32:07  jmones
 251 % Work with logical arrays from BW creation
 252 %
 253 % Revision 1.2  2004/08/11 17:39:51  jmones
 254 % Algorithm url in docs corrected.
 255 %
 256 % Revision 1.1  2004/08/11 17:34:11  jmones
 257 % poly2mask added: creates filled polygon bw mask
 258 %
 259 %