# HG changeset patch # User dbateman # Date 1188459755 0 # Node ID 956148c0d3888920358159ab2b67f09881d2b94e # Parent 7911a62a300d8f28570a5c7e68c2a0a8be8723d0 [project @ 2007-08-30 07:39:32 by dbateman] diff --git a/doc/ChangeLog b/doc/ChangeLog --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,3 +1,9 @@ +2007-08-30 David Bateman + + * interpreter/geometryimages.m: Add inpolygon example + * interpreter/Makefile.in (GEOMETRYIMAGES): Add inpolygon example. + * interpret/geometry.txi: Document inpolygon. + 2007-08-27 David Bateman * interpreter/struct.txi: Remove. diff --git a/doc/interpreter/Makefile.in b/doc/interpreter/Makefile.in --- a/doc/interpreter/Makefile.in +++ b/doc/interpreter/Makefile.in @@ -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)) diff --git a/doc/interpreter/geometry.txi b/doc/interpreter/geometry.txi --- a/doc/interpreter/geometry.txi +++ b/doc/interpreter/geometry.txi @@ -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 diff --git a/doc/interpreter/geometryimages.m b/doc/interpreter/geometryimages.m --- a/doc/interpreter/geometryimages.m +++ b/doc/interpreter/geometryimages.m @@ -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 diff --git a/scripts/ChangeLog b/scripts/ChangeLog --- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2007-08-30 David Bateman + + * geometry/inpolygon.m: New file. + * geometry/Makefile.in (SOURCES): Add inpolygon.m. + 2007-08-29 Peter A. Gustafson * plot/__go_draw_axes__.m: Disable linetype in do_linestyle_command. diff --git a/scripts/geometry/Makefile.in b/scripts/geometry/Makefile.in --- a/scripts/geometry/Makefile.in +++ b/scripts/geometry/Makefile.in @@ -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)) diff --git a/scripts/geometry/inpolygon.m b/scripts/geometry/inpolygon.m new file mode 100644 --- /dev/null +++ b/scripts/geometry/inpolygon.m @@ -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 +## Created: 14 November 2006 + +## Vectorized by Søren Hauberg + +## 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."); +