Ok, forwarding it to the matplotlib-devel list.
Best wishes,
Konrad (on behalf of our workgroup)
-------- Original Message --------
Subject: Source of inaccuracies in the matplotlib library
Date: Fri, 8 Apr 2011 18:12:47 +0200
From: Bartkowski, Konrad
<k.bartkow...@fz-juelich.de><mailto:k.bartkow...@fz-juelich.de>
To: d...@cornell.edu<mailto:d...@cornell.edu> <d...@cornell.edu><mailto:d...@cornell.edu>, md...@stsci.edu<mailto:md...@stsci.edu>
<md...@stsci.edu><mailto:md...@stsci.edu>, efir...@hawaii.edu<mailto:efir...@hawaii.edu> <efir...@hawaii.edu><mailto:efir...@hawaii.edu>,
jdhun...@ace.bsd.uchicago.edu<mailto:jdhun...@ace.bsd.uchicago.edu> <jdhun...@ace.bsd.uchicago.edu><mailto:jdhun...@ace.bsd.uchicago.edu>,
jdh2...@gmail.com<mailto:jdh2...@gmail.com> <jdh2...@gmail.com><mailto:jdh2...@gmail.com>
CC: Bartkowski, Konrad <k.bartkow...@fz-juelich.de><mailto:k.bartkow...@fz-juelich.de>,
el...@interia.pl<mailto:el...@interia.pl> <el...@interia.pl><mailto:el...@interia.pl>, Matthias Bolten
<bol...@math.uni-wuppertal.de><mailto:bol...@math.uni-wuppertal.de>, Grotendorst, Johannes
<j.grotendo...@fz-juelich.de><mailto:j.grotendo...@fz-juelich.de>, Steffen, Bernhard
<b.stef...@fz-juelich.de><mailto:b.stef...@fz-juelich.de>
Dear Matplotlib developers,
I am writing about the matplotlib library with the mpl_toolkits. First
of all let me emphasize how great software it is. Recently, in one of
our projects we were rendering big surfaces and encountered the
following problem:
http://www.mail-archive.com/matplotlib-devel@lists.sourceforge.net/msg06869.html
It's not a bug (which all in all is a natural and unavoidable ingredient
of the software, and especially in such a big and complex system like
matplotlib would be fully natural), since the software does exactly the
projection mathematics that it is expected to do, but a source of the
inaccuracies, which is especially visible in the critical examples. For
the profit of the Python community we are sending You a proposition of a
modification of the surface plotting rendering system, in case You find
it interesting enough to include in the consecutive version of the
library. In the source code from the attachment we redesigned a little
bit the computation process – since the computations are especially
sensible to numerical errors, that are for example amplified while
norming or processing the quaterions in the various stages (for example
division over coordinate in the perspective projection). Therefore the
computational focus can be shifted from the Polygon collection to the
polygons itself. In the example from the above forum or the slightly
modified one, one can observe a big difference in the numerical
precision while the speed of the computations does not decrease (at
least visibly). While instead of the surfaces from the forum, the
following surfaces are rendered:
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = 10 * np.outer(np.cos(u), np.sin(v))
y = 10 * np.outer(np.sin(u), np.sin(v))
z = 10 * np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, rstride=8, cstride=8, color='y', alpha=0.5)
shiftX=28
shiftY=28
X,Y=np.meshgrid(range(-20+shiftX,20+shiftX),range(-20+shiftY,20+shiftY))
Z=np.ones((X.shape[0], Y.shape[1]))
ax.plot_surface(X, Y, Z, color='r', rstride=10, cstride=10, alpha=1.0)
the issue is visible for example at the azimuth=40 , elevation=70 – with
those parameters the mentioned case is visible on the red surface, while
with elevation=68 not. Moreover, now also the stride is big (in the new
approach the influence of increasing stride on the numerical precision
grows).
So again let me use this opportunity to thank You for empowering the
Python community worldwide in a great, powerful scientific visualization
tool.
Best wishes,
Konrad Bartkowski
------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------
Forschungszentrum Juelich GmbH
52425 Juelich
Sitz der Gesellschaft: Juelich
Eingetragen im Handelsregister des Amtsgerichts Dueren Nr. HR B 3498
Vorsitzender des Aufsichtsrats: MinDirig Dr. Karl Eugen Huthmacher
Geschaeftsfuehrung: Prof. Dr. Achim Bachem (Vorsitzender),
Dr. Ulrich Krafft (stellv. Vorsitzender), Prof. Dr.-Ing. Harald Bolt,
Prof. Dr. Sebastian M. Schmidt
------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------
Besuchen Sie uns auf unserem neuen Webauftritt unter www.fz-juelich.de
def plot_surface(self, X, Y, Z, *args, **kwargs):
'''
Create a surface plot.
By default it will be colored in shades of a solid color,
but it also supports color mapping by supplying the *cmap*
argument.
============= ================================================
Argument Description
============= ================================================
*X*, *Y*, *Z* Data values as numpy.arrays
*rstride* Array row stride (step size)
*cstride* Array column stride (step size)
*color* Color of the surface patches
*cmap* A colormap for the surface patches.
*facecolors* Face colors for the individual patches
*norm* An instance of Normalize to map values to colors
*vmin* Minimum value to map
*vmax* Maximum value to map
*shade* Whether to shade the facecolors
============= ================================================
Other arguments are passed on to
:func:`~mpl_toolkits.mplot3d.art3d.Poly3DCollection.__init__`
'''
had_data = self.has_data()
rows, cols = Z.shape
tX, tY, tZ = np.transpose(X), np.transpose(Y), np.transpose(Z)
rstride = kwargs.pop('rstride', 10)
cstride = kwargs.pop('cstride', 10)
if 'facecolors' in kwargs:
fcolors = kwargs.pop('facecolors')
else:
color = np.array(colorConverter.to_rgba(kwargs.pop('color', 'b')))
fcolors = None
cmap = kwargs.get('cmap', None)
norm = kwargs.pop('norm', None)
vmin = kwargs.pop('vmin', None)
vmax = kwargs.pop('vmax', None)
linewidth = kwargs.get('linewidth', None)
shade = kwargs.pop('shade', cmap is None)
lightsource = kwargs.pop('lightsource', None)
# Shade the data
if shade and cmap is not None and fcolors is not None:
fcolors = self._shade_colors_lightsource(Z, cmap, lightsource)
POLYS = []
POLYC = []
COLSET=[]
CMAPCOLSET=[]
NORMALS=[]
CMAPNORMALS=[]
#colset contains the data for coloring: either average z or the facecolor
rsCounter=0
for rs in np.arange(0, rows-1, rstride):
POLYS.append([])
POLYC.append([])
COLSET.append([])
NORMALS.append([])
for cs in np.arange(0, cols-1, cstride):
ps = []
corners = []
for a, ta in [(X, tX), (Y, tY), (Z, tZ)]:
ztop = a[rs][cs:min(cols, cs+cstride+1)]
zleft = ta[min(cols-1, cs+cstride)][rs:min(rows, rs+rstride+1)]
zbase = a[min(rows-1, rs+rstride)][cs:min(cols, cs+cstride+1):]
zbase = zbase[::-1]
zright = ta[cs][rs:min(rows, rs+rstride+1):]
zright = zright[::-1]
corners.append([ztop[0], ztop[-1], zbase[0], zbase[-1]])
z = np.concatenate((ztop, zleft, zbase, zright))
ps.append(z)
# The construction leaves the array with duplicate points, which
# are removed here.
ps = zip(*ps)
lastp = np.array([])
ps2 = []
avgzsum = 0.0
for p in ps:
if p != lastp:
ps2.append(p)
lastp = p
avgzsum += p[2]
POLYS[rsCounter].append(ps2)
if fcolors is not None:
COLSET[rsCounter].append(fcolors[rs][cs])
else:
COLSET[rsCounter].append(avgzsum / len(ps2))
CMAPCOLSET.append(avgzsum / len(ps2))
# Only need vectors to shade if no cmap
if cmap is None and shade:
v1 = np.array(ps2[0]) - np.array(ps2[1])
v2 = np.array(ps2[2]) - np.array(ps2[0])
NORMALS[rsCounter].append(np.cross(v1, v2))
rsCounter+=1
rsCounter=0
for rs in np.arange(0, rows-1, rstride):
csCounter=0
for cs in np.arange(0, cols-1, cstride):
POLYC[rsCounter].append(art3d.Poly3DCollection([POLYS[rsCounter][csCounter]], *args, **kwargs))
csCounter+=1
rsCounter+=1
rsCounter=0
cmapcounter=0
if 0!=len(CMAPCOLSET):
minCMAPCOLSET=min(CMAPCOLSET)
maxCMAPCOLSET=max(CMAPCOLSET)
for rs in np.arange(0, rows-1, rstride):
csCounter=0
for cs in np.arange(0, cols-1, cstride):
if fcolors is not None:
if shade:
COLSET[rsCounter][csCounter] = self._shade_colors(np.array([COLSET[rsCounter][csCounter]]), [NORMALS[rsCounter][csCounter]])
element=COLSET[rsCounter][csCounter]
if len(element)==1:
element=element[0]
POLYC[rsCounter][csCounter].set_facecolors(element)
POLYC[rsCounter][csCounter].set_edgecolors(element)
elif cmap:
CMAPCOLSET2=[0.,]*3
CMAPCOLSET2[0]=CMAPCOLSET[cmapcounter]
CMAPCOLSET2[1]=minCMAPCOLSET
CMAPCOLSET2[2]=maxCMAPCOLSET
CMAPCOLSET2= np.array(CMAPCOLSET2)
POLYC[rsCounter][csCounter].set_array(CMAPCOLSET2)
if vmin is not None or vmax is not None:
POLYC[rsCounter][csCounter].set_clim(vmin, vmax)
if norm is not None:
POLYC[rsCounter][csCounter].set_norm(norm)
else:
if shade:
COLSET[rsCounter][csCounter]= self._shade_colors(color, [NORMALS[rsCounter][csCounter]])
else:
COLSET[rsCounter][csCounter]= color
POLYC[rsCounter][csCounter].set_facecolors(COLSET[rsCounter][csCounter])
cmapcounter+=1
csCounter+=1
rsCounter+=1
rsCounter=0
for rs in np.arange(0, rows-1, rstride):
csCounter=0
for cs in np.arange(0, cols-1, cstride):
self.add_collection(POLYC[rsCounter][csCounter])
csCounter+=1
rsCounter+=1
self.auto_scale_xyz(X, Y, Z, had_data)
return POLYC
def plot_surface(self, X, Y, Z, *args, **kwargs):
'''
Create a surface plot.
By default it will be colored in shades of a solid color,
but it also supports color mapping by supplying the *cmap*
argument.
============= ================================================
Argument Description
============= ================================================
*X*, *Y*, *Z* Data values as numpy.arrays
*rstride* Array row stride (step size)
*cstride* Array column stride (step size)
*color* Color of the surface patches
*cmap* A colormap for the surface patches.
*facecolors* Face colors for the individual patches
*norm* An instance of Normalize to map values to colors
*vmin* Minimum value to map
*vmax* Maximum value to map
*shade* Whether to shade the facecolors
============= ================================================
Other arguments are passed on to
:func:`~mpl_toolkits.mplot3d.art3d.Poly3DCollection.__init__`
'''
had_data = self.has_data()
rows, cols = Z.shape
tX, tY, tZ = np.transpose(X), np.transpose(Y), np.transpose(Z)
rstride = kwargs.pop('rstride', 10)
cstride = kwargs.pop('cstride', 10)
if 'facecolors' in kwargs:
fcolors = kwargs.pop('facecolors')
else:
color = np.array(colorConverter.to_rgba(kwargs.pop('color', 'b')))
fcolors = None
cmap = kwargs.get('cmap', None)
norm = kwargs.pop('norm', None)
vmin = kwargs.pop('vmin', None)
vmax = kwargs.pop('vmax', None)
linewidth = kwargs.get('linewidth', None)
shade = kwargs.pop('shade', cmap is None)
lightsource = kwargs.pop('lightsource', None)
# Shade the data
if shade and cmap is not None and fcolors is not None:
fcolors = self._shade_colors_lightsource(Z, cmap, lightsource)
POLYS = []
POLYC = []
COLSET=[]
CMAPCOLSET=[]
NORMALS=[]
CMAPNORMALS=[]
#colset contains the data for coloring: either average z or the facecolor
rsCounter=0
for rs in np.arange(0, rows-1, rstride):
POLYS.append([])
POLYC.append([])
COLSET.append([])
NORMALS.append([])
for cs in np.arange(0, cols-1, cstride):
ps = []
corners = []
for a, ta in [(X, tX), (Y, tY), (Z, tZ)]:
ztop = a[rs][cs:min(cols, cs+cstride+1)]
zleft = ta[min(cols-1, cs+cstride)][rs:min(rows, rs+rstride+1)]
zbase = a[min(rows-1, rs+rstride)][cs:min(cols, cs+cstride+1):]
zbase = zbase[::-1]
zright = ta[cs][rs:min(rows, rs+rstride+1):]
zright = zright[::-1]
corners.append([ztop[0], ztop[-1], zbase[0], zbase[-1]])
z = np.concatenate((ztop, zleft, zbase, zright))
ps.append(z)
# The construction leaves the array with duplicate points, which
# are removed here.
ps = zip(*ps)
lastp = np.array([])
ps2 = []
avgzsum = 0.0
for p in ps:
if p != lastp:
ps2.append(p)
lastp = p
avgzsum += p[2]
POLYS[rsCounter].append(ps2)
if fcolors is not None:
COLSET[rsCounter].append(fcolors[rs][cs])
else:
COLSET[rsCounter].append(avgzsum / len(ps2))
CMAPCOLSET.append(avgzsum / len(ps2))
# Only need vectors to shade if no cmap
if cmap is None and shade:
v1 = np.array(ps2[0]) - np.array(ps2[1])
v2 = np.array(ps2[2]) - np.array(ps2[0])
NORMALS[rsCounter].append(np.cross(v1, v2))
rsCounter+=1
rsCounter=0
for rs in np.arange(0, rows-1, rstride):
csCounter=0
for cs in np.arange(0, cols-1, cstride):
POLYC[rsCounter].append(art3d.Poly3DCollection([POLYS[rsCounter][csCounter]], *args, **kwargs))
csCounter+=1
rsCounter+=1
rsCounter=0
cmapcounter=0
if 0!=len(CMAPCOLSET):
minCMAPCOLSET=min(CMAPCOLSET)
maxCMAPCOLSET=max(CMAPCOLSET)
for rs in np.arange(0, rows-1, rstride):
csCounter=0
for cs in np.arange(0, cols-1, cstride):
if fcolors is not None:
if shade:
COLSET[rsCounter][csCounter] = self._shade_colors(np.array([COLSET[rsCounter][csCounter]]), [NORMALS[rsCounter][csCounter]])
element=COLSET[rsCounter][csCounter]
if len(element)==1:
element=element[0]
POLYC[rsCounter][csCounter].set_facecolors(element)
POLYC[rsCounter][csCounter].set_edgecolors(element)
elif cmap:
CMAPCOLSET2=[0.,]*3
CMAPCOLSET2[0]=CMAPCOLSET[cmapcounter]
CMAPCOLSET2[1]=minCMAPCOLSET
CMAPCOLSET2[2]=maxCMAPCOLSET
CMAPCOLSET2= np.array(CMAPCOLSET2)
POLYC[rsCounter][csCounter].set_array(CMAPCOLSET2)
if vmin is not None or vmax is not None:
POLYC[rsCounter][csCounter].set_clim(vmin, vmax)
if norm is not None:
POLYC[rsCounter][csCounter].set_norm(norm)
else:
if shade:
COLSET[rsCounter][csCounter]= self._shade_colors(color, [NORMALS[rsCounter][csCounter]])
else:
COLSET[rsCounter][csCounter]= color
POLYC[rsCounter][csCounter].set_facecolors(COLSET[rsCounter][csCounter])
cmapcounter+=1
csCounter+=1
rsCounter+=1
rsCounter=0
for rs in np.arange(0, rows-1, rstride):
csCounter=0
for cs in np.arange(0, cols-1, cstride):
self.add_collection(POLYC[rsCounter][csCounter])
csCounter+=1
rsCounter+=1
self.auto_scale_xyz(X, Y, Z, had_data)
return POLYC
------------------------------------------------------------------------------
Xperia(TM) PLAY
It's a major breakthrough. An authentic gaming
smartphone on the nation's most reliable network.
And it wants your games.
http://p.sf.net/sfu/verizon-sfdev
_______________________________________________
Matplotlib-devel mailing list
Matplotlib-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-devel
