Menu
▾
▴
SF.net SVN: matplotlib:[8081] trunk/matplotlib
|
From: <ef...@us...> - 2010-01-16 17:45:38
|
Revision: 8081
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=8081&view=rev
Author: efiring
Date: 2010-01-16 17:45:32 +0000 (Sat, 16 Jan 2010)
Log Message:
-----------
Patch by Ian Thomas fixes 2 major cntr.c problems:
Now line contours coincide with filled contour boundaries,
and interior masked regions are handled correctly by contourf.
Modified Paths:
--------------
trunk/matplotlib/CHANGELOG
trunk/matplotlib/examples/pylab_examples/contourf_demo.py
trunk/matplotlib/src/cntr.c
Modified: trunk/matplotlib/CHANGELOG
===================================================================
--- trunk/matplotlib/CHANGELOG 2010-01-11 19:55:26 UTC (rev 8080)
+++ trunk/matplotlib/CHANGELOG 2010-01-16 17:45:32 UTC (rev 8081)
@@ -1,10 +1,14 @@
-2009-01-11 The color of legend patch follows the rc parameters
+2010-01-16 Applied patch by Ian Thomas to fix two contouring
+ problems: now contourf handles interior masked regions,
+ and the boundaries of line and filled contours coincide. - EF
+
+2009-01-11 The color of legend patch follows the rc parameters
axes.facecolor and axes.edgecolor. -JJL
-2009-01-11 adjustable of Axes can be "box-forced" which allow
+2009-01-11 adjustable of Axes can be "box-forced" which allow
sharing axes. -JJL
-2009-01-11 Add add_click and pop_click methods in
+2009-01-11 Add add_click and pop_click methods in
BlockingContourLabeler. -JJL
2010-01-03 Added rcParams['axes.color_cycle'] - EF
Modified: trunk/matplotlib/examples/pylab_examples/contourf_demo.py
===================================================================
--- trunk/matplotlib/examples/pylab_examples/contourf_demo.py 2010-01-11 19:55:26 UTC (rev 8080)
+++ trunk/matplotlib/examples/pylab_examples/contourf_demo.py 2010-01-16 17:45:32 UTC (rev 8081)
@@ -3,35 +3,14 @@
origin = 'lower'
#origin = 'upper'
-# The following controls only interior masking.
-test_masking = False # There is a bug in filled contour masking with
- # interior masks.
+delta = 0.025
-if test_masking:
- # Use a coarse grid so only a few masked points are needed.
- delta = 0.5
-else:
- delta = 0.025
-
x = y = arange(-3.0, 3.01, delta)
X, Y = meshgrid(x, y)
Z1 = bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)
Z2 = bivariate_normal(X, Y, 1.5, 0.5, 1, 1)
Z = 10 * (Z1 - Z2)
-# interior badmask doesn't work yet for filled contours
-if test_masking:
- badmask = zeros(shape(Z))
-
- badmask[5,5] = 1
- badmask[5,6] = 1
- Z[5,5] = 0
- Z[5,6] = 0
-
- badmask[0,0] = 1
- Z[0,0] = 0
- Z = ma.array(Z, mask=badmask)
-
nr, nc = Z.shape
# put NaNs in one corner:
@@ -43,7 +22,11 @@
# mask another corner:
Z[:nr//6, :nc//6] = ma.masked
+# mask a circle in the middle:
+interior = sqrt((X**2) + (Y**2)) < 0.5
+Z[interior] = ma.masked
+
# We are using automatic selection of contour levels;
# this is usually not such a good idea, because they don't
# occur on nice boundaries, but we do it here for purposes
@@ -63,7 +46,7 @@
origin=origin,
hold='on')
-title('Nonsense (with 2 masked corners)')
+title('Nonsense (3 masked regions)')
xlabel('word length anomaly')
ylabel('sentence length anomaly')
@@ -87,7 +70,7 @@
colors = ('k',),
linewidths = (3,),
origin = origin)
-title('Listed colors (with 2 masked corners)')
+title('Listed colors (3 masked regions)')
clabel(CS4, fmt = '%2.1f', colors = 'w', fontsize=14)
colorbar(CS3)
Modified: trunk/matplotlib/src/cntr.c
===================================================================
--- trunk/matplotlib/src/cntr.c 2010-01-11 19:55:26 UTC (rev 8080)
+++ trunk/matplotlib/src/cntr.c 2010-01-16 17:45:32 UTC (rev 8081)
@@ -50,12 +50,14 @@
* The problem is that two disjoint curves cut through a saddle zone
* (I reject the alternative of connecting the opposite points to make
* a single self-intersecting curve, since those make ugly contour plots
- * -- I've tried it). The real problem with saddle zones is that you
- * need to communicate the connectivity decision you make back to the
- * calling routine, since for the next contour level, we need to tell
- * the contour tracer to make the same decision as on the previous
- * level. The input/output triangulation array is the solution to this
- * nasty problem.
+ * -- I've tried it). The solution is to determine the z value of the
+ * centre of the zone, which is the mean of the z values of the four
+ * corner points. If the centre z is higher than the contour level of
+ * interest and you are moving along the line with higher values on the
+ * left, turn right to leave the saddle zone. If the centre z is lower
+ * than the contour level turn left. Whether the centre z is higher
+ * than the 1 or 2 contour levels is stored in the saddle array so that
+ * it does not need to be recalculated in subsequent passes.
*
* Another complicating factor is that there may be logical holes in
* the mesh -- zones which do not exist. We want our contours to stop
@@ -175,6 +177,11 @@
* or not, z value 0, 1, or 2 -- is kept in a mesh sized data array */
typedef short Cdata;
+/* information to decide on correct contour direction in saddle zones
+ * is stored in a mesh sized array. Only those entries corresponding
+ * to saddle zones have nonzero values in this array. */
+typedef char Saddle;
+
/* here is the minimum structure required to tell where we are in the
* mesh sized data array */
typedef struct Csite Csite;
@@ -189,8 +196,8 @@
long count; /* count of start markers visited */
double zlevel[2]; /* contour levels, zlevel[1]<=zlevel[0]
* signals single level case */
- short *triangle; /* triangulation array for the mesh */
- char *reg; /* region array for the mesh (was int) */
+ Saddle *saddle; /* saddle zone information for the mesh */
+ char *reg; /* region array for the mesh (was int) */
Cdata *data; /* added by EF */
long edge0, left0; /* starting site on this curve for closure */
int level0; /* starting level for closure */
@@ -225,8 +232,6 @@
printf("\n");
}
-/* triangle only takes values of -1, 0, 1, so it could be a signed char. */
-/* most or all of the longs probably could be converted to ints with no loss */
/* the Cdata array consists of the following bits:
* Z_VALUE (2 bits) 0, 1, or 2 function value at point
@@ -243,6 +248,7 @@
* OPEN_END marks an i-edge start point whose other endpoint is
* on a boundary for the single level case
* ALL_DONE marks final start point
+ * SLIT_DN_VISITED this slit downstroke hasn't/has been visited in pass 2
*/
#define Z_VALUE 0x0003
#define ZONE_EX 0x0004
@@ -257,6 +263,7 @@
#define SLIT_DN 0x0800
#define OPEN_END 0x1000
#define ALL_DONE 0x2000
+#define SLIT_DN_VISITED 0x4000
/* some helpful macros to find points relative to a given directed
* edge -- points are designated 0, 1, 2, 3 CCW around zone with 0 and
@@ -272,6 +279,15 @@
enum {kind_zone, kind_edge1, kind_edge2,
kind_slit_up, kind_slit_down, kind_start_slit=16} point_kinds;
+/* Saddle zone array consists of the following bits:
+ * SADDLE_SET whether zone's saddle data has been set.
+ * SADDLE_GT0 whether z of centre of zone is higher than site->level[0].
+ * SADDLE_GT1 whether z of centre of zone is higher than site->level[1].
+ */
+#define SADDLE_SET 0x01
+#define SADDLE_GT0 0x02
+#define SADDLE_GT1 0x04
+
/* ------------------------------------------------------------------------ */
/* these actually mark points */
@@ -313,18 +329,17 @@
long left0 = site->left0;
int level0 = site->level0 == level;
int two_levels = site->zlevel[1] > site->zlevel[0];
- short *triangle = site->triangle;
+ Saddle* saddle = site->saddle;
const double *x = pass2 ? site->x : 0;
const double *y = pass2 ? site->y : 0;
- const double *z = pass2 ? site->z : 0;
- double zlevel = pass2 ? site->zlevel[level] : 0.0;
+ const double *z = site->z;
+ double zlevel = site->zlevel[level];
double *xcp = pass2 ? site->xcp : 0;
double *ycp = pass2 ? site->ycp : 0;
short *kcp = pass2 ? site->kcp : 0;
int z0, z1, z2, z3;
- int keep_left = 0; /* flag to try to minimize curvature in saddles */
int done = 0;
int n_kind;
@@ -402,29 +417,28 @@
{
if (z1 == z3)
{
- /* this is a saddle zone, need triangle to decide
- * -- set triangle if not already decided for this zone */
+ /* this is a saddle zone, determine whether to turn left or
+ * right depending on height of centre of zone relative to
+ * contour level. Set saddle[zone] if not already decided. */
long zone = edge + (left > 0 ? left : 0);
- if (triangle)
+ if (!(saddle[zone] & SADDLE_SET))
{
- if (!triangle[zone])
- {
- if (keep_left)
- triangle[zone] = jedge ? -1 : 1;
- else
- triangle[zone] = jedge ? 1 : -1;
- }
- if (triangle[zone] > 0 ? !jedge : jedge)
- goto bkwd;
+ saddle[zone] = SADDLE_SET;
+ double zcentre = (z[p0] + z[p0+left] + z[p1] + z[p1+left])/4.0;
+ if (zcentre > site->zlevel[0])
+ saddle[zone] |=
+ (two_levels && zcentre > site->zlevel[1])
+ ? SADDLE_GT0 | SADDLE_GT1 : SADDLE_GT0;
}
- else
- {
- if (keep_left)
- goto bkwd;
- }
+
+ int turnRight = level == 2 ? (saddle[zone] & SADDLE_GT1)
+ : (saddle[zone] & SADDLE_GT0);
+ if (z1 ^ (level == 2))
+ turnRight = !turnRight;
+ if (!turnRight)
+ goto bkwd;
}
/* bend forward (right along curve) */
- keep_left = 1;
jedge = !jedge;
edge = p1 + (left > 0 ? left : 0);
{
@@ -437,7 +451,6 @@
{
bkwd:
/* bend backward (left along curve) */
- keep_left = 0;
jedge = !jedge;
edge = p0 + (left > 0 ? left : 0);
{
@@ -590,17 +603,27 @@
if (n_kind) kcp[n_kind] += kind_start_slit;
return slit_cutter (site, 0, pass2);
}
+ if (fwd < 0 && level0 && left < 0)
+ {
+ if (n_kind) kcp[n_kind] += kind_start_slit;
+ return slit_cutter (site, 0, pass2);
+ }
return 3;
}
else if (pass2)
{
if (heads_up || (fwd < 0 && (data[edge] & SLIT_DN)))
{
- site->edge = edge;
- site->left = left;
- site->n = n + marked;
- if (n_kind) kcp[n_kind] += kind_start_slit;
- return slit_cutter (site, heads_up, pass2);
+ if (!heads_up && !(data[edge] & SLIT_DN_VISITED))
+ data[edge] |= SLIT_DN_VISITED;
+ else
+ {
+ site->edge = edge;
+ site->left = left;
+ site->n = n + marked;
+ if (n_kind) kcp[n_kind] += kind_start_slit;
+ return slit_cutter (site, heads_up, pass2);
+ }
}
}
else
@@ -1181,6 +1204,8 @@
/* place immediate stop mark if nothing found */
if (!count)
data[0] |= ALL_DONE;
+ else
+ for (i = 0; i < ijmax; ++i) site->saddle[i] = 0;
/* initialize site */
site->edge0 = site->edge00 = site->edge = 0;
@@ -1252,7 +1277,7 @@
if (site == NULL) return NULL;
site->data = NULL;
site->reg = NULL;
- site->triangle = NULL;
+ site->saddle = NULL;
site->xcp = NULL;
site->ycp = NULL;
site->kcp = NULL;
@@ -1268,7 +1293,6 @@
{
long ijmax = iMax * jMax;
long nreg = iMax * jMax + iMax + 1;
- long i;
site->imax = iMax;
site->jmax = jMax;
@@ -1278,21 +1302,20 @@
PyMem_Free(site);
return -1;
}
- site->triangle = (short *) PyMem_Malloc(sizeof(short) * ijmax);
- if (site->triangle == NULL)
+ site->saddle = (Saddle*) PyMem_Malloc(sizeof(Saddle) * ijmax);
+ if (site->saddle == NULL)
{
PyMem_Free(site->data);
PyMem_Free(site);
return -1;
}
- for (i = 0; i < ijmax; i++) site->triangle[i] = 0;
site->reg = NULL;
if (mask != NULL)
{
site->reg = (char *) PyMem_Malloc(sizeof(char) * nreg);
if (site->reg == NULL)
{
- PyMem_Free(site->triangle);
+ PyMem_Free(site->saddle);
PyMem_Free(site->data);
PyMem_Free(site);
return -1;
@@ -1311,7 +1334,7 @@
void cntr_del(Csite *site)
{
- PyMem_Free(site->triangle);
+ PyMem_Free(site->saddle);
PyMem_Free(site->reg);
PyMem_Free(site->data);
PyMem_Free(site);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
