I should read entire posts before sending people down the wrong
pathways. I just happened to be working on a Python equivalent to MATLAB
"triplot" stuff when I read your subject line and made the wrong
assumptions. That program does just plot the edges, as you noted.
I have attached a python program, much of which is a translation of a
program I found at M*B central, contributed from the outside. Given a
triangulation, it allows you to interpolate on a regular rectangular
grid (dx=constant, dy=another constant). In your case, it should allow
you to use your original triangulation, and should avoid the convex hull
artifacts of your original griddata plot. I do not know if this new
program will give you a figure that will look as good as your latest
based on John's suggestion or not.
Just a few warnings:
1. I use an environment that (like matlab but very unpythonic) imports
everything, so the attached may need to import one or two more routines
2. I have only tested with generated data, using triangulations created
with delaunay, but it worked okay as long as my random points did not
lead to bad triangles.
3. The routines use a lot of loops, so they are relatively slow - should
be an ideal program for Cython though.
3. method='linear' uses linear interp,and is very stable. It should not
be influenced in any way by convexity.
method='quadratic' uses a quadratic method, which does well except in
the presence of very thin triangles. I suspect that derivative
estimation at the vertices has problems on the edges and may very well
be influenced by the convexity.
linear on a fine grid will probably look pretty good for you.
4. to get something as nice as the mayavi sample, you will probably need
a fine rectangular mesh (maybe 200x200).
I tested the tinterp program in octave for a number of problems and
found that the quadratic option gave nice results when I did not have
functional discontinuities, and when I had good triangulations. I have
been meaning to translate this to Python, so thanks for the motivation.
At the very least, this program may help you appreciate John's
suggestions more.
Cheers,
Eric
Ondrej Certik wrote:
Ok, I made a progress, it seems it's working. Script and picture
Forgot to attach the script.
Ondrej
------------------------------------------------------------------------
------------------------------------------------------------------------------
Register Now for Creativity and Technology (CaT), June 3rd, NYC. CaT
is a gathering of tech-side developers & brand creativity professionals. Meet
the minds behind Google Creative Lab, Visual Complexity, Processing, &
iPhoneDevCamp asthey present alongside digital heavyweights like Barbarian
Group, R/GA, & Big Spaceship. http://www.creativitycat.com
------------------------------------------------------------------------
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
import numpy as np
from numpy import ones,shape, int32, floor, ceil, array, dot, zeros,reshape
from numpy import logical_not, NaN, any, vstack, sum, mgrid, sin,cos
def mesh_grid_tsearch(X,Y,dx,dy,points,tris):
"""
locate the specific triangle for point in a uniform mesh
INPUT:
X,Y are the x- and y- coodinates from the mesh generated with
meshgrid or mgrid or ogrid (X,Y must have the same shape)
dx is the constant distance between x grid lines and
dy is the constant distance between y grid lines
(dx does not have to equal dy)
points are the vertices of triangulation (shape(points)==[n_points,2])
tris are the triangles, three integer index values pointing to rows from points
shape(tris)==[n_triangles,3]
OUTPUT:
triangle number (row of tris) for each point of the mesh,
or -1 if the mesh point is not located in any triangle
"""
## Compute vectors
tri_loc=int32(-ones(shape(X)))
nr,nc=shape(X)
x_min=np.min(X)
x_max=np.max(X)
y_min=np.min(Y)
y_max=np.max(Y)
n_tris = shape(tris)[0]
for n in range(n_tris):
r=tris[n,:]
A=points[r[0],:]
B=points[r[1],:]
C=points[r[2],:]
tri_x_min=min(A[0],B[0],C[0])
tri_y_min=min(A[1],B[1],C[1])
tri_x_max=max(A[0],B[0],C[0])
tri_y_max=max(A[1],B[1],C[1])
# determination of imin, etc only place where grid regularity required
i_min = int32(max(0,floor((tri_x_min-x_min)/dx)))
j_min = int32(max(0,floor((tri_y_min-y_min)/dy)))
i_max = int32(min(nc,ceil((tri_x_max-x_min)/dx)+1))
j_max = int32(min(nr,ceil((tri_y_max-y_min)/dy)+1))
if not (i_min>=nc or i_max<=0 or j_min>=nc or j_max<=0):
for i in range(i_min,i_max):
for j in range(j_min,j_max):
if tri_loc[i,j]<0:
P=array([X[i,j],Y[i,j]])
v0 = C - A
v1 = B - A
v2 = P - A
## Compute dot products
dot00 = dot(v0, v0)
dot01 = dot(v0, v1)
dot02 = dot(v0, v2)
dot11 = dot(v1, v1)
dot12 = dot(v1, v2)
## Compute barycentric coordinates
invDenom = 1.0 / (dot00 * dot11 - dot01 * dot01)
u = (dot11 * dot02 - dot01 * dot12) * invDenom
v = (dot00 * dot12 - dot01 * dot02) * invDenom
## Check if point is in triangle
if (u >= 0) and (v >= 0) and (u + v <= 1):
tri_loc[i,j]=int32(n)
return tri_loc
def tinterp(p,t,f,xi,yi,method='linear'):
"""
Interpolation of scattered, triangulated data on a uniform grid.
INPUTS:
p = [x1,y1 x2,y2 etc] - an Nx2 array defining the [x,y]
coordinates of the scattered vertices
t = [n11,n12,n13 n21,n22,n23 etc] - an Mx3 array defining the
triangulation of the points in p, as
returned by delaunayn(p)
f = [f1 f2 etc] - an Nx1 vector defining the value
of the function at the scattered
vertices
xi,yi - arrays of identical
size, defining the [x,y] coordinates
where the function is to be
interpolated
method = 'linear','quadratic' - Type of interpolant, method is
optional and if not passed the
'linear' method is used
OUTPUT:
fi - The interpolation evaluated at
all points in xi,yi.
The 'linear' method should return identical results to griddata. The
'quadratic' method is NEW! (I think) and should return more accurate results
when compared with the 'cubic' method in gridata (See tester.m).
Example:
p = rand(100,2)
t = delaunay(p)
f = exp(-20*sum((p-0.5)**2,2))
figure, trimesh(t,p(:,1),p(:,2),f), title('Exact solution on scattered data')
x = linspace(0,1,50)
x,y = meshgrid(x,x)
fL = tinterp(p,t,f,x,y)
fQ = tinterp(p,t,f,x,y,'quadratic')
figure, mesh(x,y,fL), title('Linear interpolation')
figure, mesh(x,y,fQ), title('Quadratic interpolation')
Original MATLAB VERSION:
Any comments? Let me know:
d_engwi...@hotmail.com
Darren Engwirda - 2006
Quick Python Translation from M*B
Eric Carlson - 2009
"""
# Check I/O
# Error checking
sizx = shape(xi)
sizy = shape(yi)
dx=xi[1,0]-xi[0,0]
dy=yi[0,1]-yi[0,0]
# Find enclosing triangle of points in pii
i=mesh_grid_tsearch(xi,yi,dx,dy,p,t)
ig = i.flatten()
pii = vstack((xi.flatten(),yi.flatten())).T
# Allocate output
fi = zeros(shape(pii)[0])
# i = tsearch(p[:,0],p[:,1],t,pii[:,0],pii[:,1])
# Deal with points oustide convex hull
in1 = logical_not(ig<0)
# Keep internal points
pin = pii[in1,:]
tin = t[ig[in1],:]
# Corner nodes
t1 = tin[:,0]
t2 = tin[:,1]
t3 = tin[:,2]
# Linear shape functions
dp1 = pin-p[t1,:]
dp2 = pin-p[t2,:]
dp3 = pin-p[t3,:]
A3 = abs(dp1[:,0]*dp2[:,1]-dp1[:,1]*dp2[:,0])
A2 = abs(dp1[:,0]*dp3[:,1]-dp1[:,1]*dp3[:,0])
A1 = abs(dp3[:,0]*dp2[:,1]-dp3[:,1]*dp2[:,0])
# Interpolant
if method=='linear': # Linear interpolation
# Take weighted average of constant extrapolations
print 'linear interpolation'
fi[in1] = (A1*f[t1]+A2*f[t2]+A3*f[t3]) / (A1+A2+A3)
else: # Quadratic interpolation
print 'quadratic interpolation'
# Form gradients
fx,fy = tgrad(p,t,f)
# Take weighted average of linear extrapolations
fi[in1] = ( A1*(f[t1] + dp1[:,0]*fx[t1] + dp1[:,1]*fy[t1]) +
A2*(f[t2] + dp2[:,0]*fx[t2] + dp2[:,1]*fy[t2]) +
A3*(f[t3] + dp3[:,0]*fx[t3] + dp3[:,1]*fy[t3]) )/ (A1+A2+A3)
# Deal with outside points
if any(logical_not(in1)):
fi[logical_not(in1)] = NaN
fi = np.ma.masked_where(np.isnan(fi),fi)
# Reshape output
fi = reshape(fi,sizx)
return fi
#######################################################################################
def tgrad(p,t,f):
"""
Approximate [df/dx,df/dy] at the vertices of triangulation.
"""
# Nodes
t1 = t[:,0]
t2 = t[:,1]
t3 = t[:,2]
# Evaluate centroidal gradients (piecewise-linear interpolants)
x23 = p[t2,0]-p[t3,0]
y23 = p[t2,1]-p[t3,1]
x21 = p[t2,0]-p[t1,0]
y21 = p[t2,1]-p[t1,1]
# Centroidal values - [df/dx,df/dy]
den = (x23*y21-x21*y23)
fcx = (y23*f[t1]+(y21-y23)*f[t2]-y21*f[t3])/den
fcy = -(x23*f[t1]+(x21-x23)*f[t2]-x21*f[t3])/den
# Calculate simplex quality
q = quality(p,t)
# Form nodal gradients.
# Take the quality weighted average of the neighbouring centroidal values.
# Low quality simplices ("thin" triangles) can contribute inaccurate
# gradient information - quality weighting seems to limit this effect.
fnx = 0.0*f
fny = 0.0*f
den = 0.0*f
for k in range(shape(t)[0]):
# Nodes
n1 = t1[k]
n2 = t2[k]
n3 = t3[k]
# Current values
qk = q[k]
qfx = qk*fcx[k]
qfy = qk*fcy[k]
# Average to n1
fnx[n1] = fnx[n1]+qfx
fny[n1] = fny[n1]+qfy
den[n1] = den[n1]+qk
# Average to n2
fnx[n2] = fnx[n2]+qfx
fny[n2] = fny[n2]+qfy
den[n2] = den[n2]+qk
# Average to n3
fnx[n3] = fnx[n3]+qfx
fny[n3] = fny[n3]+qfy
den[n3] = den[n3]+qk
fnx = fnx/den
fny = fny/den
return fnx,fny
######################################################################################
def quality(p,t):
"""
Approximate simplex quality of triangulation
"""
# Nodes
p1 = p[t[:,0],:]
p2 = p[t[:,1],:]
p3 = p[t[:,2],:]
# Approximate quality
d12 = p2-p1
d13 = p3-p1
q = 3.4641*abs(d12[:,0]*d13[:,1]-d12[:,1]*d13[:,0])/(sum(d12**2,1)+sum(d13**2,1)+sum((p3-p2)**2,1))
return q
import matplotlib.delaunay as delaunay
import matplotlib.pyplot as pp
import scipy as sp
x, y = mgrid[0.0:2.0:101j,0.0:2.0:101j]
dx=x[1,0]-x[0,0]
dy=y[0,1]-y[0,0]
npts=105
#using scipy version >= .7
#scipy version<.7, sp.random(npts)
xpt=sp.random.random_sample(npts)*1.75+.125
ypt=sp.random.random_sample(npts)*1.75+.125
#create triangulation
pts=vstack((xpt,ypt)).T
my_func=lambda x,y:sin(x**2)*cos(y*x)
f = my_func(xpt,ypt)
circumcenters, edges, tri_points, tri_neighbors = delaunay.delaunay(xpt, ypt)
fa1=tinterp(pts,tri_points,f,x,y,method='linear')
pp.contourf(x,y,fa1)
#pp.plot(xpt[edges.T],ypt[edges.T],'k--', linewidth=.5)
pp.figure()
fa=tinterp(pts,tri_points,f,x,y,method='quadratic')
pp.contourf(x,y,fa)
#pp.plot(xpt[edges.T],ypt[edges.T],'k--', linewidth=.5)
pp.figure()
pp.contourf(x,y,my_func(x,y))
pp.figure()
##tri1 = delaunay.Triangulation(xpt,ypt)
### interpolate data
##interp1 = tri1.nn_interpolator(f)
##pp.contourf(x,y,interp1(x,y))
import griddata
pp.contourf(x,y,pp.mlab.griddata(xpt, ypt, f, x, y))
##figure()
##contourf(x,y,fa-fa1)
##pp.plot(xpt[edges.T],ypt[edges.T],'k--', linewidth=.5)
##for n in range(shape(tri_points)[0]):
## r=tri_points[n,:]
## A=pts[r[0],:]
## B=pts[r[1],:]
## C=pts[r[2],:]
## the_center = (A+B+C)/3
## text(the_center[0], the_center[1],str(n))
#plot(x,y,'k.')
------------------------------------------------------------------------------
Register Now for Creativity and Technology (CaT), June 3rd, NYC. CaT
is a gathering of tech-side developers & brand creativity professionals. Meet
the minds behind Google Creative Lab, Visual Complexity, Processing, &
iPhoneDevCamp asthey present alongside digital heavyweights like Barbarian
Group, R/GA, & Big Spaceship. http://www.creativitycat.com
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
