[Numpy-discussion] Creating parallel curves
Andrea Gavana
my apologies for my deep ignorance about math stuff; I guess I
should be able to find this out but I keep getting impossible results.
Basically I have a set of x, y data (around 1,000 elements each) and I
want to create 2 parallel "curves" (offset curves) to the original
one; "parallel" means curves which are displaced from the base curve
by a constant offset, either positive or negative, in the direction of
the curve's normal. Something like this:
http://pyx.sourceforge.net/examples/drawing2/parallel.html
The point in the x variable are monotonically increasing.
For every point in the curve, I also have the angle this point creates
to the vertical (y) direction, so I naively thought of doing this:
angles = numpy.deg2rad(angles)
x_low = x - DISTANCE*numpy.cos(angles)
y_low = y + DISTANCE*numpy.sin(angles)
x_high = x + DISTANCE*numpy.cos(angles)
y_high = y - DISTANCE*numpy.sin(angles)
But by plotting these thing out with matplotlib it seems to me they
don't really look very parallel nor very constant-distance. I admit
the matplotlib axis scales can play a significant role here, but I was
wondering if some of you could spot my dumb mistake (and maybe provide
a way to nicely take into account the different scales of the
matplotlib axes so that the curves really look parallel...).
Thank you in advance for any suggestion.
Andrea.
"Imagination Is The Only Weapon In The War Against Reality."
http://xoomer.alice.it/infinity77/
_______________________________________________
NumPy-Discussion mailing list
NumPy-Di...@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Chris Barker
> Basically I have a set of x, y data (around 1,000 elements each) and I
> want to create 2 parallel "curves" (offset curves) to the original
> one; "parallel" means curves which are displaced from the base curve
> by a constant offset, either positive or negative, in the direction of
> the curve's normal. Something like this:
>
> http://pyx.sourceforge.net/examples/drawing2/parallel.html
THis is called "buffering" in GIS parlance -- there are functions
available to do it in GIS an computational geometry libraries: you
might look in the shapely package:
https://github.com/sgillies/shapely
or CGAL
If the overhead of these packages is too much, and you still want to
write your own code, try googling:
"buffering a line GIS algorithm" or something like that, and you'll
find pointers.
> But by plotting these thing out with matplotlib it seems to me they
> don't really look very parallel nor very constant-distance.
as we say on the wxPython list -- post a fully functional example, so
we can check it out.
-Chris
--
Christopher Barker, Ph.D.
Oceanographer
Emergency Response Division
NOAA/NOS/OR&R (206) 526-6959 voice
7600 Sand Point Way NE (206) 526-6329 fax
Seattle, WA 98115 (206) 526-6317 main reception
Andrea Gavana
On 10 February 2012 17:53, Chris Barker wrote:
> Andrea,
>
>> Basically I have a set of x, y data (around 1,000 elements each) and I
>> want to create 2 parallel "curves" (offset curves) to the original
>> one; "parallel" means curves which are displaced from the base curve
>> by a constant offset, either positive or negative, in the direction of
>> the curve's normal. Something like this:
>>
>> http://pyx.sourceforge.net/examples/drawing2/parallel.html
>
> THis is called "buffering" in GIS parlance -- there are functions
> available to do it in GIS an computational geometry libraries: you
> might look in the shapely package:
>
> https://github.com/sgillies/shapely
Thanks, I hoped this one would prove itself useful, but unfortunately
to me it looks like one of the most impenetrable Python code I have
ever seen. Or maybe my Python is just too weak. The end result is the
same.
I have surfed quite a lot and found many reference to Bezier curves,
with a lot of mathematical acrobatics but little of practical value
(i.e., code :-) ). In the end I stumbled upon a nice function in the
matplotlib library (obviously) which can give me the normal line to
every point in the curve, and I almost got there.
I still have 2 problems (attached sample), picture at
http://img689.imageshack.us/img689/7246/exampleplot.png
1) If you run the attached sample, you'll see a plot of a cubic
polynomial with two "almost" parallel lines to it on the first
subplot. I say "almost" because at the inflection points of the cubic
something funny happens (see suplot number 2);
2) You can see in subplot 3 that the 3 lines are nicely shown as
"almost" parallel (except the problem at point 1), as the X and Y axes
have the same scale. Unfortunately I can't keep the same scales in my
real plot and I can't use axis=square in matplotlib either. How do I
"normalize" the get_normal_points method to get the same visual
appearance of parallelism on X and Y having different X and Y scales
(i.e., the same X and Y scales as in subplot 1)?
>> But by plotting these thing out with matplotlib it seems to me they
>> don't really look very parallel nor very constant-distance.
>
> as we say on the wxPython list -- post a fully functional example, so
> we can check it out.
Attached as promised.
Jonathan Hilmer
Here is how to do it with splines. I would be more standard to return
an array of normals, rather than two arrays of x and y components, but
it actually requires less housekeeping this way. As an aside, I would
prefer to work with rotations via matrices, but it looks like there's
no support for that built in to Numpy or Scipy?
def normal_vectors(x, y, scalar=1.0):
tck = scipy.interpolate.splrep(x, y)
y_deriv = scipy.interpolate.splev(x, tck, der=1)
normals_rad = np.arctan(y_deriv)+np.pi/2.
return np.cos(normals_rad)*scalar, np.sin(normals_rad)*scalar
#
Jonathan
Charles R Harris
Hi All,
my apologies for my deep ignorance about math stuff; I guess I
should be able to find this out but I keep getting impossible results.
Basically I have a set of x, y data (around 1,000 elements each) and I
want to create 2 parallel "curves" (offset curves) to the original
one; "parallel" means curves which are displaced from the base curve
by a constant offset, either positive or negative, in the direction of
the curve's normal. Something like this:
Note that curves produced in this way aren't actually 'parallel' and can even cross themselves.
<snip>
Chuck
Andrea Gavana
On 12 February 2012 20:53, Jonathan Hilmer wrote:
> Andrea,
>
> Here is how to do it with splines. I would be more standard to return
> an array of normals, rather than two arrays of x and y components, but
> it actually requires less housekeeping this way. As an aside, I would
> prefer to work with rotations via matrices, but it looks like there's
> no support for that built in to Numpy or Scipy?
>
> def normal_vectors(x, y, scalar=1.0):
> tck = scipy.interpolate.splrep(x, y)
>
> y_deriv = scipy.interpolate.splev(x, tck, der=1)
>
> normals_rad = np.arctan(y_deriv)+np.pi/2.
>
> return np.cos(normals_rad)*scalar, np.sin(normals_rad)*scalar
Thank you for this, I'll give it a go in a few minutes (hopefully I
will also be able to correctly understand what you did). One thing
though, at first glance, it appears to me that your approach is very
similar to mine (meaning it will give "parallel" curves that cross
themselves as in the example I posted). But maybe I am wrong, my
apologies if I missed something.
Thank you so much for your answer.
Andrea.
Andrea Gavana
I know, my definition of "parallel" was probably not orthodox enough.
What I am looking for is to generate 2 curves that look "graphically
parallel enough" to the original one, and not "parallel" in the true
mathematical sense.
Charles R Harris
Jonathan,
Thank you for this, I'll give it a go in a few minutes (hopefully I
On 12 February 2012 20:53, Jonathan Hilmer wrote:
> Andrea,
>
> Here is how to do it with splines. I would be more standard to return
> an array of normals, rather than two arrays of x and y components, but
> it actually requires less housekeeping this way. As an aside, I would
> prefer to work with rotations via matrices, but it looks like there's
> no support for that built in to Numpy or Scipy?
>
> def normal_vectors(x, y, scalar=1.0):
> tck = scipy.interpolate.splrep(x, y)
>
> y_deriv = scipy.interpolate.splev(x, tck, der=1)
>
> normals_rad = np.arctan(y_deriv)+np.pi/2.
>
> return np.cos(normals_rad)*scalar, np.sin(normals_rad)*scalar
will also be able to correctly understand what you did). One thing
though, at first glance, it appears to me that your approach is very
similar to mine (meaning it will give "parallel" curves that cross
themselves as in the example I posted). But maybe I am wrong, my
apologies if I missed something.
Thank you so much for your answer.
Crossing curves is the correct behavior for this method, think propagating wavefronts and the appearance of light on the bottom of a swimming pool when there are ripples on the surface. I think you need to define what you really want from 'parallel'.
Chuck
Charles R Harris
Charles,
I know, my definition of "parallel" was probably not orthodox enough.
On 12 February 2012 21:00, Charles R Harris wrote:
>
>
> On Fri, Feb 10, 2012 at 9:38 AM, Andrea Gavana <andrea...@gmail.com>
> wrote:
>>
>> Hi All,
>>
>> my apologies for my deep ignorance about math stuff; I guess I
>> should be able to find this out but I keep getting impossible results.
>>
>> Basically I have a set of x, y data (around 1,000 elements each) and I
>> want to create 2 parallel "curves" (offset curves) to the original
>> one; "parallel" means curves which are displaced from the base curve
>> by a constant offset, either positive or negative, in the direction of
>> the curve's normal. Something like this:
>>
>
> Note that curves produced in this way aren't actually 'parallel' and can
> even cross themselves.
What I am looking for is to generate 2 curves that look "graphically
parallel enough" to the original one, and not "parallel" in the true
mathematical sense.
You could try setting a point and 'contracting' the curve towards the point. A point a infinity would give the usual parallel curves. There are probably a lot of perspective like transformations that would do something similar.
Chuck
Jonathan Hilmer
I realized that my answer wouldn't be complete, but as people have
pointed out that's a substantially more difficult question, so I
wanted to give you a complete answer to just a subset of your problem.
I'm currently writing a variant that avoids the overlapping normal
vectors by interatively 1.) expanding along normals then 2.)
condensing points, for every iteration. However, I'm doing it mostly
for my own interest since I'm pretty sure it will not be functional
when complete: your problem is that calculation of derivatives/normals
is going to become unstable in acute convex regions, and the overlap
issue there will become more severe. I would strongly recommend
adapting some existing library for this problem.
Jonathan
Andrea Gavana
On 12 February 2012 21:59, Jonathan Hilmer wrote:
> Andrea,
>
> I realized that my answer wouldn't be complete, but as people have
> pointed out that's a substantially more difficult question, so I
> wanted to give you a complete answer to just a subset of your problem.
>
> I'm currently writing a variant that avoids the overlapping normal
> vectors by interatively 1.) expanding along normals then 2.)
> condensing points, for every iteration. However, I'm doing it mostly
> for my own interest since I'm pretty sure it will not be functional
> when complete: your problem is that calculation of derivatives/normals
> is going to become unstable in acute convex regions, and the overlap
> issue there will become more severe. I would strongly recommend
> adapting some existing library for this problem.
Thank you for your clear explanation.
I feared something like that, but then again I was almost expecting
it... I do wonder, however, if the problem I am having could be
somehow simplified/modified by looking at another situation, in 3D
this time: streamtubes and 3D representation of "tubes"/"pipes" are
relatively common (i.e., VTKTubeFilter from the VTK library,
streamtubes in Mayavi). They all build some kind of "cylindrical"
shape around the main path (curve in 3D) to give the visual effect of
some "tube" surrounding the main path.
What I am trying to do is basically the same thing but in 2D: now, I
have no idea if this can be done, and not even how it could be done.
My math is relatively weak, I can't fathom how to build a 3D tube let
alone splatter it back on a 2D map to give that effect.
I am not sure my original problem and this one are related or not, but
again every suggestion is most welcome.
Thank you.
Andrea.
--
Andrea.
"Imagination Is The Only Weapon In The War Against Reality."
http://xoomer.alice.it/infinity77/
>>> import PyQt4.QtGui
Traceback (most recent call last):
File "<interactive input>", line 1, in <module>
ImportError: No module named PyQt4.QtGui
>>>
>>> import pygtk
Traceback (most recent call last):
File "<interactive input>", line 1, in <module>
ImportError: No module named pygtk
>>>
>>> import wx
Robert Kern
> I know, my definition of "parallel" was probably not orthodox enough.
> What I am looking for is to generate 2 curves that look "graphically
> parallel enough" to the original one, and not "parallel" in the true
> mathematical sense.
There is a rigorous way to define the curve that you are looking for,
and fortunately it gives some hints for implementation. For each point
(x,y) in space, associate with it the nearest distance D from that
point to the reference curve. The "parallel" curves are just two sides
of the level set where D(x,y) is equal to the specified distance
(possibly removing the circular caps that surround the ends of the
reference curve).
If performance is not a constraint, then you could just evaluate that
D(x,y) function on a fine-enough grid and do marching squares to find
the level set. matplotlib's contour plotting routines can help here.
There is a hint in the PyX page that you linked to that you should
consider. Angles in the reference curve become circular arcs in the
"parallel" curves. So if your reference curve is just a bunch of line
segments, then what you can do is take each line segment, and make
parallel copies the same length to either side. Now you just need to
connect up these parallel segments with each other. You do this by
using circular arcs centered on the vertices of the reference curve.
Do this on both sides. On the "outer" side, the arcs will go "forward"
while on the "inner" side, the arcs will go "backwards" just like the
cusps that you saw in your attempt. Now let's take care of that. You
will have two self-intersecting curves consisting of alternating line
segments and circular arcs. Parts of these curves will be too close to
the reference curve. You will have to go through these curves to find
the locations of self-intersection and remove the parts of the
segments and arcs that are too close to the reference curve. This is
tricky to do, but the formulae for segment-segment, segment-arc, and
arc-arc intersection can be found online.
--
Robert Kern
"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
-- Umberto Eco
Niki Spahiev
make cython interface to it.
HTH
Niki
Chris Barker
> You can get polygon buffer from http://angusj.com/delphi/clipper.php and
> make cython interface to it.
This should be built into GEOS as well, and the shapely package
provides a python wrapper already.
-Chris
> HTH
>
> Niki
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Di...@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
--
Christopher Barker, Ph.D.
Oceanographer
Emergency Response Division
NOAA/NOS/OR&R (206) 526-6959 voice
7600 Sand Point Way NE (206) 526-6329 fax
Seattle, WA 98115 (206) 526-6317 main reception
Andrea Gavana
From: "Andrea Gavana" <andrea...@gmail.com>
Date: Feb 13, 2012 11:31 PM
Subject: Re: [Numpy-discussion] Creating parallel curves
To: "Jonathan Hilmer" <jkhi...@gmail.com>
Thank you Jonathan for this, it's exactly what I was looking for. I' ll try it tomorrow on the 768 well trajectories I have and I'll let you know if I stumble upon any issue.
If someone could shed some light on my problem number 2 (how to adjust the scaling/distance) so that the curves look parallel on a matplotlib graph even though the axes scales are different, I'd be more than grateful.
Thank you in advance.
Andrea.
Andrea,
This is playing some tricks with 2D array expansion to make a tradeoff
in memory for speed. Given two sets of matching vectors (one
reference, given first, and a newly-expanded one, given second), it
removes all points from the expanded vectors that aren't needed to
describe the new contour.
def filter_expansion(x, y, x_expan, y_expan, distance_target, tol=1e-6):
target_xx, expansion_xx = scipy.meshgrid(x, x_expan)
target_yy, expansion_yy = scipy.meshgrid(y, y_expan)
distance = scipy.sqrt((expansion_yy - target_yy)**2 + (expansion_xx -
target_xx)**2)
valid = distance.min(axis=1) > distance_target*(1.-tol)
return x_expan.compress(valid), y_expan.compress(valid)
#
Jonathan
Kiko
---------- Forwarded message ----------
From: "Andrea Gavana" <andrea...@gmail.com>
Date: Feb 13, 2012 11:31 PM
Subject: Re: [Numpy-discussion] Creating parallel curves
To: "Jonathan Hilmer" <jkhi...@gmail.com>Thank you Jonathan for this, it's exactly what I was looking for. I' ll try it tomorrow on the 768 well trajectories I have and I'll let you know if I stumble upon any issue.
If someone could shed some light on my problem number 2 (how to adjust the scaling/distance) so that the curves look parallel on a matplotlib graph even though the axes scales are different, I'd be more than grateful.
Thank you in advance.
Hi.
Maybe this could help you as a starting point.
from Shapely.geometry import LineString
from matplotlib import pyplot
myline = LineString(...)
x, y = myline.xy
xx, yy = myline.buffer(distancefrommyline).exterior.xy # coordinates around myline
pyplot.plot(x, y)
pyplot.plot(xx,yy)
pyplot.show()
Best.