Mercurial > octave-nkf

changeset 6847:956148c0d388

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-08-30 07:39:32 by dbateman]
author dbateman
date Thu, 30 Aug 2007 07:42:35 +0000
parents 7911a62a300d
children 9de60a998cf3
files doc/ChangeLog doc/interpreter/Makefile.in doc/interpreter/geometry.txi doc/interpreter/geometryimages.m scripts/ChangeLog scripts/geometry/Makefile.in scripts/geometry/inpolygon.m
diffstat 7 files changed, 184 insertions(+), 7 deletions(-) [+]
line wrap: on
line diff
--- a/doc/ChangeLog	Thu Aug 30 07:25:57 2007 +0000
+++ b/doc/ChangeLog	Thu Aug 30 07:42:35 2007 +0000
@@ -1,3 +1,9 @@
+2007-08-30  David Bateman  <dbateman@free.fr>
+
+	* interpreter/geometryimages.m: Add inpolygon example
+	* interpreter/Makefile.in (GEOMETRYIMAGES): Add inpolygon example.
+	* interpret/geometry.txi: Document inpolygon.
+
 2007-08-27  David Bateman  <dbateman@free.fr>
 
 	* interpreter/struct.txi: Remove.
--- a/doc/interpreter/Makefile.in	Thu Aug 30 07:25:57 2007 +0000
+++ b/doc/interpreter/Makefile.in	Thu Aug 30 07:42:35 2007 +0000
@@ -43,7 +43,7 @@
 
 EXAMPLE_FILES = $(addprefix $(top_srcdir)/examples/, $(EXAMPLE_FILES_NODIR))
 
-GEOMETRYIMAGES = voronoi triplot griddata convhull delaunay
+GEOMETRYIMAGES = voronoi triplot griddata convhull delaunay inpolygon
 GEOMETRYIMAGES_EPS = $(addsuffix .eps, $(GEOMETRYIMAGES))
 GEOMETRYIMAGES_PDF = $(addsuffix .pdf, $(GEOMETRYIMAGES))
 GEOMETRYIMAGES_PNG = $(addsuffix .png, $(GEOMETRYIMAGES))
--- a/doc/interpreter/geometry.txi	Thu Aug 30 07:25:57 2007 +0000
+++ b/doc/interpreter/geometry.txi	Thu Aug 30 07:42:35 2007 +0000
@@ -298,8 +298,9 @@
 @end float
 @end ifnotinfo
 
-Additional information about the size of the facets of a Voronoi diagram
-can be had with the @code{polyarea} function.
+Additional information about the size of the facets of a Voronoi
+diagram, and which points of a set of points is in a polygon can be had
+with the @code{polyarea} and @code{inpolygon} functions respectively.
 
 @DOCSTRING(polyarea)
 
@@ -320,6 +321,34 @@
 
 Facets of the Voronoi diagram with a vertex at infinity have infinity area.
 
+@DOCSTRING(inpolygon)
+
+An example of the use of @code{inpolygon} might be
+
+@example
+@group
+randn ("state", 2);
+x = randn (100, 1);
+y = randn (100, 1);
+vx = cos (pi * [-1 : 0.1: 1]);
+vy = sin (pi * [-1 : 0.1 : 1]);
+in = inpolygon (x, y, vx, vy);
+plot(vx, vy, x(in), y(in), "r+", x(!in), y(!in), "bo");
+axis ([-2, 2, -2, 2]);
+@end group
+@end example
+
+@ifnotinfo
+@noindent
+The result of which can be seen in @ref{fig:inpolygon}.
+
+@float Figure,fig:inpolygon
+@image{inpolygon,8cm}
+@caption{Demonstration of the @code{inpolygon} function to determine the
+points inside a polygon}
+@end float
+@end ifnotinfo
+
 @node Convex Hull
 @section Convex Hull
 
--- a/doc/interpreter/geometryimages.m	Thu Aug 30 07:25:57 2007 +0000
+++ b/doc/interpreter/geometryimages.m	Thu Aug 30 07:42:35 2007 +0000
@@ -41,14 +41,24 @@
     print (strcat (nm, ".", typ), strcat ("-d", typ)) 
   elseif (strcmp (nm, "delaunay"))
     rand ("state", 1);
-    x = rand (10, 1);
-    y = rand (10, 1);
+    x = rand (1, 10);
+    y = rand (1, 10);
     T = delaunay (x, y);
     X = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ];
     Y = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ];
     axis ([0, 1, 0, 1]);
     plot(X, Y, "b", x, y, "r*");
     print (strcat (nm, ".", typ), strcat ("-d", typ)) 
+  elseif (strcmp (nm, "inpolygon"))
+    randn ("state", 2);
+    x = randn (100, 1);
+    y = randn (100, 1);
+    vx = cos (pi * [-1 : 0.1: 1]);
+    vy = sin (pi * [-1 : 0.1 : 1]);
+    in = inpolygon (x, y, vx, vy);
+    plot(vx, vy, x(in), y(in), "r+", x(!in), y(!in), "bo");
+    axis ([-2, 2, -2, 2]);
+    print (strcat (nm, ".", typ), strcat ("-d", typ)) 
   else
     error ("unrecognized plot requested");
   endif
--- a/scripts/ChangeLog	Thu Aug 30 07:25:57 2007 +0000
+++ b/scripts/ChangeLog	Thu Aug 30 07:42:35 2007 +0000
@@ -1,3 +1,8 @@
+2007-08-30  David Bateman  <dbateman@free.fr>
+
+	* geometry/inpolygon.m: New file.
+	* geometry/Makefile.in (SOURCES): Add inpolygon.m.
+	
 2007-08-29  Peter A. Gustafson  <petegus@umich.edu>
 
 	* plot/__go_draw_axes__.m: Disable linetype in do_linestyle_command.
--- a/scripts/geometry/Makefile.in	Thu Aug 30 07:25:57 2007 +0000
+++ b/scripts/geometry/Makefile.in	Thu Aug 30 07:42:35 2007 +0000
@@ -21,8 +21,8 @@
 INSTALL_DATA = @INSTALL_DATA@
 
 SOURCES = convhull.m delaunay3.m delaunayn.m delaunay.m dsearch.m dsearchn.m \
-	griddata.m griddata3.m griddatan.m trimesh.m triplot.m tsearchn.m \
-	voronoi.m voronoin.m
+	griddata.m griddata3.m griddatan.m inpolygon.m trimesh.m triplot.m \
+	tsearchn.m voronoi.m voronoin.m
 
 DISTFILES = $(addprefix $(srcdir)/,Makefile.in $(SOURCES))
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/inpolygon.m	Thu Aug 30 07:42:35 2007 +0000
@@ -0,0 +1,127 @@
+## Copyright (C) 2006 Frederick (Rick) A Niles, Søren Hauberg
+##
+## This file is part of Octave.
+##
+## Octave 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, or (at your option)
+## any later version.
+##
+## Octave 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 Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{in}, @var{on}] = } inpolygon (@var{x}, @var{y}, @var{xv}, @var{xy})
+##
+## For a polygon defined by @code{(@var{xv}, @var{yv})} points, determine
+## if the points @code{(@var{x}, @var{y})} are inside or outside the polygon.
+## The variables @var{x}, @var{y}, must have the same dimension. The optional
+## output @var{on} gives the points that are on the polygon.
+##
+## @end deftypefn
+
+## Author: Frederick (Rick) A Niles <niles@rickniles.com>
+## Created: 14 November 2006
+
+## Vectorized by Søren Hauberg <soren@hauberg.org>
+
+## The method for determining if a point is in in a polygon is based on
+## the algorithm shown on
+## http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/ and is
+## credited to Randolph Franklin.
+
+function [IN, ON] = inpolygon (X, Y, xv, yv)
+
+  if ( !(isreal(X) && isreal(Y) && ismatrix(Y) && ...
+	 ismatrix(Y) && size_equal(X,Y)) )
+    error ("inpolygon: first two arguments must be real matrices of same size");
+  elseif ( !(isreal(xv) && isreal(yv) && isvector(xv) && ...
+	     isvector(yv) && size_equal(xv,yv)) )
+    error ("inpolygon: last two arguments must be real vectors of same size");
+  endif
+
+  npol = length(xv);
+  do_boundary = (nargout >= 2);
+  
+  IN = zeros (size(X), "logical");
+  if (do_boundary) 
+    ON = zeros (size(X), "logical"); 
+  endif
+  
+  j = npol;
+  for i = 1 : npol
+    delta_xv = xv(j) - xv(i);
+    delta_yv = yv(j) - yv(i);
+    ## distance = [distance from (X,Y) to edge] * length(edge)
+    distance = delta_xv .* (Y - yv(i)) - (X - xv(i)) .* delta_yv;
+    ##
+    ## is Y between the y-values of edge i,j
+    ##        AND (X,Y) on the left of the edge ?
+    idx1 = ((yv(i) <= Y & Y < yv(j)) | (yv(j) <= Y & Y < yv(i)) ) & ...
+           0 < distance.*delta_yv;
+    IN (idx1) = !IN (idx1);
+
+    ## Check if (X,Y) are actually ON the boundary of the polygon.
+    if (do_boundary)
+       idx2 = ((yv(i) <= Y & Y <= yv(j)) | (yv(j) <= Y & Y <= yv(i))) & ...
+              ((xv(i) <= X & X <= xv(j)) | (xv(j) <= X & X <= xv(i))) & ...
+              (0 == distance | !delta_xv);
+       ON (idx2) = true;
+    endif
+    j = i;
+  endfor
+endfunction
+
+%!demo
+%!  xv=[ 0.05840, 0.48375, 0.69356, 1.47478, 1.32158, \
+%!       1.94545, 2.16477, 1.87639, 1.18218, 0.27615, \
+%!       0.05840 ];
+%!  yv=[ 0.60628, 0.04728, 0.50000, 0.50000, 0.02015, \
+%!       0.18161, 0.78850, 1.13589, 1.33781, 1.04650, \
+%!       0.60628 ];
+%! xa=[0:0.1:2.3];
+%! ya=[0:0.1:1.4];
+%! [x,y]=meshgrid(xa,ya);
+%! [IN,ON]=inpolygon(x,y,xv,yv);
+%! 
+%! inside=IN & !ON;
+%! plot(xv,yv)
+%! hold on
+%! plot(x(inside),y(inside),"@g")
+%! plot(x(~IN),y(~IN),"@m")
+%! plot(x(ON),y(ON),"@b")
+%! hold off
+%! disp("Green points are inside polygon, magenta are outside,");
+%! disp("and blue are on boundary.");
+
+%!demo
+%!  xv=[ 0.05840, 0.48375, 0.69356, 1.47478, 1.32158, \
+%!       1.94545, 2.16477, 1.87639, 1.18218, 0.27615, \
+%!       0.05840, 0.73295, 1.28913, 1.74221, 1.16023, \
+%!       0.73295, 0.05840 ];
+%!  yv=[ 0.60628, 0.04728, 0.50000, 0.50000, 0.02015, \
+%!       0.18161, 0.78850, 1.13589, 1.33781, 1.04650, \
+%!       0.60628, 0.82096, 0.67155, 0.96114, 1.14833, \
+%!       0.82096, 0.60628];
+%! xa=[0:0.1:2.3];
+%! ya=[0:0.1:1.4];
+%! [x,y]=meshgrid(xa,ya);
+%! [IN,ON]=inpolygon(x,y,xv,yv);
+%! 
+%! inside=IN & ~ ON;
+%! plot(xv,yv)
+%! hold on
+%! plot(x(inside),y(inside),"@g")
+%! plot(x(~IN),y(~IN),"@m")
+%! plot(x(ON),y(ON),"@b")
+%! hold off
+%! disp("Green points are inside polygon, magenta are outside,");
+%! disp("and blue are on boundary.");
+