Re: [OctDev] signed distance of point from polygon

[c. <car...@gm...>]

On 10 Sep 2012, at 13:20, JuanPi wrote:

> On Mon, Sep 10, 2012 at 1:01 PM, c. <car...@gm...> wrote:
>> Hi,
>> 
>> I need to compute the "signed distance" of a point P from a polygon Q,
>> i.e., the distance of P from the polyline including Q with a "+" sign if
>> the point is inside the polygon and a "-" sign if it is outside.
>> 
>> I see "distancePointPolygon" in the geometry package computes the distance
>> but always returns a positive number. Is there already a function to do
>> what I need?
>> 
>> If not, any idea where I can look-up the algorithm to modify "distancePointPolygon"
>> according to my needs?
>> 
>> Thanks,
>> c.
> 
> Hi,
> 
> I would not modify the function, but use it together with inpolygon
> (in the core geomtry package).
> 
> I guess this is what you want
> http://agora.octave.org/snippet/FRXH/
> 
> Let me know

I found that the function "distancePointPolyline" contains a for loop that can be easily vectorized away.
As I work with a large number of points this has a strong impact on running time for me.
Do you mind if I commit this patch?
c.

Index: inst/polygons2d/distancePointPolyline.m
===================================================================
--- inst/polygons2d/distancePointPolyline.m	(revision 10796)
+++ inst/polygons2d/distancePointPolyline.m	(working copy)
@@ -1,5 +1,6 @@
## Copyright (C) 2003-2011 David Legland <dav...@gr...>
## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <car...@if...>
+## Copyright (C) 2012 Carlo de Falco (Speed up by vectorization)
## All rights reserved.
##
## Redistribution and use in source and binary forms, with or without
@@ -51,17 +52,14 @@
  # allocate memory for result
  minDist = inf * ones(Np, 1);

-  # process each point
-  for p = 1:Np
-      # construct the set of edges
-      edges = [poly(1:end-1, :) poly(2:end, :)];
+  ## construct the set of edges
+  edges = [poly(1:end-1, :) poly(2:end, :)];

-      # compute distance between current each point and all edges
-      dist = distancePointEdge(point(p, :), edges);
+  ## compute distance between current each point and all edges
+  dist = distancePointEdge(point, edges);

-      # update distance if necessary
-      minDist(p) = min(dist);
-  end
+  ## update distance if necessary
+  minDist = min(dist, [], 2);

  # process output arguments
  if nargout<=1