1 /*
2 cntr.c
3 General purpose contour tracer for quadrilateral meshes.
4 Handles single level contours, or region between a pair of levels.
5
6 The routines that do all the work, as well as the explanatory
7 comments, came from gcntr.c, part of the GIST package. The
8 original mpl interface was also based on GIST. The present
9 interface uses parts of the original, but places them in
10 the entirely different framework of a Python type. It
11 was written by following the Python "Extending and Embedding"
12 tutorial.
13 */
14
15 /*
16 LICENSE AGREEMENT FOR MATPLOTLIB 0.84
17 --------------------------------------
18
19 1. This LICENSE AGREEMENT is between John D. Hunter ("JDH"), and the
20 Individual or Organization ("Licensee") accessing and otherwise using
21 matplotlib software in source or binary form and its associated
22 documentation.
23
24 2. Subject to the terms and conditions of this License Agreement, JDH
25 hereby grants Licensee a nonexclusive, royalty-free, world-wide license
26 to reproduce, analyze, test, perform and/or display publicly, prepare
27 derivative works, distribute, and otherwise use matplotlib 0.84
28 alone or in any derivative version, provided, however, that JDH's
29 License Agreement and JDH's notice of copyright, i.e., "Copyright (c)
30 2002-2005 John D. Hunter; All Rights Reserved" are retained in
31 matplotlib 0.84 alone or in any derivative version prepared by
32 Licensee.
33
34 3. In the event Licensee prepares a derivative work that is based on or
35 incorporates matplotlib 0.84 or any part thereof, and wants to
36 make the derivative work available to others as provided herein, then
37 Licensee hereby agrees to include in any such work a brief summary of
38 the changes made to matplotlib 0.84.
39
40 4. JDH is making matplotlib 0.84 available to Licensee on an "AS
41 IS" basis. JDH MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR
42 IMPLIED. BY WAY OF EXAMPLE, BUT NOT LIMITATION, JDH MAKES NO AND
43 DISCLAIMS ANY REPRESENTATION OR WARRANTY OF MERCHANTABILITY OR FITNESS
44 FOR ANY PARTICULAR PURPOSE OR THAT THE USE OF MATPLOTLIB 0.84
45 WILL NOT INFRINGE ANY THIRD PARTY RIGHTS.
46
47 5. JDH SHALL NOT BE LIABLE TO LICENSEE OR ANY OTHER USERS OF MATPLOTLIB
48 0.84 FOR ANY INCIDENTAL, SPECIAL, OR CONSEQUENTIAL DAMAGES OR
49 LOSS AS A RESULT OF MODIFYING, DISTRIBUTING, OR OTHERWISE USING
50 MATPLOTLIB 0.84, OR ANY DERIVATIVE THEREOF, EVEN IF ADVISED OF
51 THE POSSIBILITY THEREOF.
52
53 6. This License Agreement will automatically terminate upon a material
54 breach of its terms and conditions.
55
56 7. Nothing in this License Agreement shall be deemed to create any
57 relationship of agency, partnership, or joint venture between JDH and
58 Licensee. This License Agreement does not grant permission to use JDH
59 trademarks or trade name in a trademark sense to endorse or promote
60 products or services of Licensee, or any third party.
61
62 8. By copying, installing or otherwise using matplotlib 0.84,
63 Licensee agrees to be bound by the terms and conditions of this License
64 Agreement.
65 */
66
67 #include <Python.h>
68 #include "structmember.h"
69 #include <stdlib.h>
70 #include <stdio.h>
71
72 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
73
74 #include "numpy/arrayobject.h"
75
76 /* Note that all arrays in these routines are Fortran-style,
77 in the sense that the "i" index varies fastest; the dimensions
78 of the corresponding C array are z[jmax][imax] in the notation
79 used here. We can identify i and j with the x and y dimensions,
80 respectively.
81
82 */
83
84 /* What is a contour?
85 *
86 * Given a quadrilateral mesh (x,y), and values of a z at the points
87 * of that mesh, we seek a set of polylines connecting points at a
88 * particular value of z. Each point on such a contour curve lies
89 * on an edge of the mesh, at a point linearly interpolated to the
90 * contour level z0 between the given values of z at the endpoints
91 * of the edge.
92 *
93 * Identifying these points is easy. Figuring out how to connect them
94 * into a curve -- or possibly a set of disjoint curves -- is difficult.
95 * Each disjoint curve may be either a closed circuit, or it may begin
96 * and end on a mesh boundary.
97 *
98 * One of the problems with a quadrilateral mesh is that when the z
99 * values at one pair of diagonally opposite points lie below z0, and
100 * the values at the other diagonal pair of the same zone lie above z0,
101 * all four edges of the zone are cut, and there is an ambiguity in
102 * how we should connect the points. I call this a saddle zone.
103 * The problem is that two disjoint curves cut through a saddle zone
104 * (I reject the alternative of connecting the opposite points to make
105 * a single self-intersecting curve, since those make ugly contour plots
106 * -- I've tried it). The real problem with saddle zones is that you
107 * need to communicate the connectivity decision you make back to the
108 * calling routine, since for the next contour level, we need to tell
109 * the contour tracer to make the same decision as on the previous
110 * level. The input/output triangulation array is the solution to this
111 * nasty problem.
112 *
113 * Another complicating factor is that there may be logical holes in
114 * the mesh -- zones which do not exist. We want our contours to stop
115 * if they hit the edge of such a zone, just as if they'd hit the edge
116 * of the whole mesh. The input region array addresses this issue.
117 *
118 * Yet another complication: We may want a list of closed polygons which
119 * outline the region between two contour levels z0 and z1. These may
120 * include sections of the mesh boundary (including edges of logical
121 * holes defined by the region array), in addition to sections of the
122 * contour curves at one or both levels. This introduces a huge
123 * topological problem -- if one of the closed contours (possibly
124 * including an interior logical hole in the mesh, but not any part of
125 * the boundary of the whole mesh) encloses a region which is not
126 * between z0 and z1, that curve must be connected by a slit (or "branch
127 * cut") to the enclosing curve, so that the list of disjoint polygons
128 * we return is each simply connected.
129 *
130 * Okay, one final stunning difficulty: For the two level case, no
131 * individual polygon should have more than a few thousand sides, since
132 * huge filled polygons place an inordinate load on rendering software,
133 * which needs an amount of scratch space proportional to the number
134 * of sides it needs to fill. So in the two level case, we want to
135 * chunk the mesh into rectangular pieces of no more than, say, 30x30
136 * zones, which keeps each returned polygon to less than a few thousand
137 * sides (the worst case is very very bad -- you can easily write down
138 * a function and two level values which produce a polygon that cuts
139 * every edge of the mesh twice).
140 */
141
142 /*
143 * Here is the numbering scheme for points, edges, and zones in
144 * the mesh -- note that each ij corresponds to one point, one zone,
145 * one i-edge (i=constant edge) and one j-edge (j=constant edge):
146 *
147 * (ij-1)-------(ij)-------(ij)
148 * | |
149 * | |
150 * | |
151 * (ij-1) (ij) (ij)
152 * | |
153 * | |
154 * | |
155 * (ij-iX-1)----(ij-iX)----(ij-iX)
156 *
157 * At each point, the function value is either 0, 1, or 2, depending
158 * on whether it is below z0, between z0 and z1, or above z1.
159 * Each zone either exists (1) or not (0).
160 * From these three bits of data, all of the curve connectivity follows.
161 *
162 * The tracing algorithm is naturally edge-based: Either you are at a
163 * point where a level cuts an edge, ready to step across a zone to
164 * another edge, or you are drawing the edge itself, if it happens to
165 * be a boundary with at least one section between z0 and z1.
166 *
167 * In either case, the edge is a directed edge -- either the zone
168 * you are advancing into is to its left or right, or you are actually
169 * drawing it. I always trace curves keeping the region between z0 and
170 * z1 to the left of the curve. If I'm tracing a boundary, I'm always
171 * moving CCW (counter clockwise) around the zone that exists. And if
172 * I'm about to cross a zone, I'll make the direction of the edge I'm
173 * sitting on be such that the zone I'm crossing is to its left.
174 *
175 * I start tracing each curve near its lower left corner (mesh oriented
176 * as above), which is the first point I encounter scanning through the
177 * mesh in order. When I figure the 012 z values and zonal existence,
178 * I also mark the potential starting points: Each edge may harbor a
179 * potential starting point corresponding to either direction, so there
180 * are four start possibilities at each ij point. Only the following
181 * possibilities need to be marked as potential starting edges:
182 *
183 * +-+-+-+
184 * | | | |
185 * A-0-C-+ One or both levels cut E and have z=1 above them, and
186 * | EZ| | 0A is cut and either 0C is cut or CD is cut.
187 * +-B-D-+ Or, one or both levels cut E and E is a boundary edge.
188 * | | | | (and Z exists)
189 * +-+-+-+
190 *
191 * +-+-+-+
192 * | | | |
193 * +-A-0-C One or both levels cut E and have z=1 below them, and
194 * | |ZE | 0A is cut and either 0C is cut or CD is cut.
195 * +-+-B-D Or, one or both levels cut E and E is a boundary edge.
196 * | | | | (and Z exists)
197 * +-+-+-+
198 *
199 * +-+-+-+
200 * | | | |
201 * +-+-+-+ E is a boundary edge, Z exists, at some point on E
202 * | |Z| | lies between the levels.
203 * +-+E+-+
204 * | | | |
205 * +-+-+-+
206 *
207 * +-+-+-+
208 * | | | |
209 * +-+E+-+ E is a boundary edge, Z exists, at some point on E
210 * | |Z| | lies between the levels.
211 * +-+-+-+
212 * | | | |
213 * +-+-+-+
214 *
215 * During the first tracing pass, the start mark is erased whenever
216 * any non-starting edge is encountered, reducing the number of points
217 * that need to be considered for the second pass. The first pass
218 * makes the basic connectivity decisions. It figures out how many
219 * disjoint curves there will be, and identifies slits for the two level
220 * case or open contours for the single level case, and removes all but
221 * the actual start markers. A second tracing pass can perform the
222 * actual final trace.
223 */
224
225 /* ------------------------------------------------------------------------ */
226
227 /* the data about edges, zones, and points -- boundary or not, exists
228 * or not, z value 0, 1, or 2 -- is kept in a mesh sized data array */
229 typedef short Cdata;
230
231 /* here is the minimum structure required to tell where we are in the
232 * mesh sized data array */
233 typedef struct Csite Csite;
234 struct Csite
235 {
236 long edge; /* ij of current edge */
237 long left; /* +-1 or +-imax as the zone is to right, left, below,
238 * or above the edge */
239 long imax; /* imax for the mesh */
240 long jmax; /* jmax for the mesh */
241 long n; /* number of points marked on this curve so far */
242 long count; /* count of start markers visited */
243 double zlevel[2]; /* contour levels, zlevel[1]<=zlevel[0]
244 * signals single level case */
245 short *triangle; /* triangulation array for the mesh */
246 char *reg; /* region array for the mesh (was int) */
247 Cdata *data; /* added by EF */
248 long edge0, left0; /* starting site on this curve for closure */
249 int level0; /* starting level for closure */
250 long edge00; /* site needing START_ROW mark */
251
252 /* making the actual marks requires a bunch of other stuff */
253 const double *x, *y, *z; /* mesh coordinates and function values */
254 double *xcp, *ycp; /* output contour points */
255 };
256
257 #if 0
258 static void print_Csite(Csite *Csite)
259 {
260 Cdata *data = Csite->data;
261 int i, j, ij;
262 int nd = Csite->imax * (Csite->jmax + 1) + 1;
263 printf("zlevels: %8.2lg %8.2lg\n", Csite->zlevel[0], Csite->zlevel[1]);
264 printf("edge %ld, left %ld, n %ld, count %ld, edge0 %ld, left0 %ld\n",
265 Csite->edge, Csite->left, Csite->n, Csite->count,
266 Csite->edge0, Csite->left0);
267 printf(" level0 %d, edge00 %ld\n", Csite->level0, Csite->edge00);
268 printf("%04x\n", data[nd-1]);
269 for (j = Csite->jmax; j >= 0; j--)
270 {
271 for (i=0; i < Csite->imax; i++)
272 {
273 ij = i + j * Csite->imax;
274 printf("%04x ", data[ij]);
275 }
276 printf("\n");
277 }
278 printf("\n");
279 }
280 #endif
281
282 /* triangle only takes values of -1, 0, 1, so it could be a signed char. */
283 /* most or all of the longs probably could be converted to ints with no loss */
284
285 /* the Cdata array consists of the following bits:
286 * Z_VALUE (2 bits) 0, 1, or 2 function value at point
287 * ZONE_EX 1 zone exists, 0 zone doesn't exist
288 * I_BNDY this i-edge (i=constant edge) is a mesh boundary
289 * J_BNDY this j-edge (i=constant edge) is a mesh boundary
290 * I0_START this i-edge is a start point into zone to left
291 * I1_START this i-edge is a start point into zone to right
292 * J0_START this j-edge is a start point into zone below
293 * J1_START this j-edge is a start point into zone above
294 * START_ROW next start point is in current row (accelerates 2nd pass)
295 * SLIT_UP marks this i-edge as the beginning of a slit upstroke
296 * SLIT_DN marks this i-edge as the beginning of a slit downstroke
297 * OPEN_END marks an i-edge start point whose other endpoint is
298 * on a boundary for the single level case
299 * ALL_DONE marks final start point
300 */
301 #define Z_VALUE 0x0003
302 #define ZONE_EX 0x0004
303 #define I_BNDY 0x0008
304 #define J_BNDY 0x0010
305 #define I0_START 0x0020
306 #define I1_START 0x0040
307 #define J0_START 0x0080
308 #define J1_START 0x0100
309 #define START_ROW 0x0200
310 #define SLIT_UP 0x0400
311 #define SLIT_DN 0x0800
312 #define OPEN_END 0x1000
313 #define ALL_DONE 0x2000
314
315 /* some helpful macros to find points relative to a given directed
316 * edge -- points are designated 0, 1, 2, 3 CCW around zone with 0 and
317 * 1 the endpoints of the current edge */
318 #define FORWARD(left,ix) ((left)>0?((left)>1?1:-(ix)):((left)<-1?-1:(ix)))
319 #define POINT0(edge,fwd) ((edge)-((fwd)>0?fwd:0))
320 #define POINT1(edge,fwd) ((edge)+((fwd)<0?fwd:0))
321 #define IS_JEDGE(edge,left) ((left)>0?((left)>1?1:0):((left)<-1?1:0))
322 #define ANY_START (I0_START|I1_START|J0_START|J1_START)
323 #define START_MARK(left) \
324 ((left)>0?((left)>1?J1_START:I1_START):((left)<-1?J0_START:I0_START))
325
326 /* ------------------------------------------------------------------------ */
327
328 /* these actually mark points */
329 static int zone_crosser (Csite * site, int level, int pass2);
330 static int edge_walker (Csite * site, int pass2);
331 static int slit_cutter (Csite * site, int up, int pass2);
332
333 /* this calls the first three to trace the next disjoint curve
334 * -- return value is number of points on this curve, or
335 * 0 if there are no more curves this pass
336 * -(number of points) on first pass if:
337 * this is two level case, and the curve closed on a hole
338 * this is single level case, curve is open, and will start from
339 * a different point on the second pass
340 * -- in both cases, this curve will be combined with another
341 * on the second pass */
342 static long curve_tracer (Csite * site, int pass2);
343
344 /* this initializes the data array for curve_tracer */
345 static void data_init (Csite * site, int region, long nchunk);
346
347 /* ------------------------------------------------------------------------ */
348
349 /* zone_crosser assumes you are sitting at a cut edge about to cross
350 * the current zone. It always marks the initial point, crosses at
351 * least one zone, and marks the final point. On non-boundary i-edges,
352 * it is responsible for removing start markers on the first pass. */
353 static int
zone_crosser(Csite * site,int level,int pass2)354 zone_crosser (Csite * site, int level, int pass2)
355 {
356 Cdata * data = site->data;
357 long edge = site->edge;
358 long left = site->left;
359 long n = site->n;
360 long fwd = FORWARD (left, site->imax);
361 long p0, p1;
362 int jedge = IS_JEDGE (edge, left);
363 long edge0 = site->edge0;
364 long left0 = site->left0;
365 int level0 = site->level0 == level;
366 int two_levels = site->zlevel[1] > site->zlevel[0];
367 short *triangle = site->triangle;
368
369 const double *x = pass2 ? site->x : 0;
370 const double *y = pass2 ? site->y : 0;
371 const double *z = pass2 ? site->z : 0;
372 double zlevel = pass2 ? site->zlevel[level] : 0.0;
373 double *xcp = pass2 ? site->xcp : 0;
374 double *ycp = pass2 ? site->ycp : 0;
375
376 int z0, z1, z2, z3;
377 int keep_left = 0; /* flag to try to minimize curvature in saddles */
378 int done = 0;
379
380 if (level)
381 level = 2;
382
383 for (;;)
384 {
385 /* set edge endpoints */
386 p0 = POINT0 (edge, fwd);
387 p1 = POINT1 (edge, fwd);
388
389 /* always mark cut on current edge */
390 if (pass2)
391 {
392 /* second pass actually computes and stores the point */
393 double zcp = (zlevel - z[p0]) / (z[p1] - z[p0]);
394 xcp[n] = zcp * (x[p1] - x[p0]) + x[p0];
395 ycp[n] = zcp * (y[p1] - y[p0]) + y[p0];
396 }
397 if (!done && !jedge)
398 {
399 if (n)
400 {
401 /* if this is not the first point on the curve, and we're
402 * not done, and this is an i-edge, check several things */
403 if (!two_levels && !pass2 && (data[edge] & OPEN_END))
404 {
405 /* reached an OPEN_END mark, skip the n++ */
406 done = 4; /* same return value 4 used below */
407 break;
408 }
409
410 /* check for curve closure -- if not, erase any start mark */
411 if (edge == edge0 && left == left0)
412 {
413 /* may signal closure on a downstroke */
414 if (level0)
415 done = (!pass2 && two_levels && left < 0) ? 5 : 3;
416 }
417 else if (!pass2)
418 {
419 Cdata start =
420 data[edge] & (fwd > 0 ? I0_START : I1_START);
421 if (start)
422 {
423 data[edge] &= ~start;
424 site->count--;
425 }
426 if (!two_levels)
427 {
428 start = data[edge] & (fwd > 0 ? I1_START : I0_START);
429 if (start)
430 {
431 data[edge] &= ~start;
432 site->count--;
433 }
434 }
435 }
436 }
437 }
438 n++;
439 if (done)
440 break;
441
442 /* cross current zone to another cut edge */
443 z0 = (data[p0] & Z_VALUE) != level; /* 1 if fill toward p0 */
444 z1 = !z0; /* know level cuts edge */
445 z2 = (data[p1 + left] & Z_VALUE) != level;
446 z3 = (data[p0 + left] & Z_VALUE) != level;
447 if (z0 == z2)
448 {
449 if (z1 == z3)
450 {
451 /* this is a saddle zone, need triangle to decide
452 * -- set triangle if not already decided for this zone */
453 long zone = edge + (left > 0 ? left : 0);
454 if (triangle)
455 {
456 if (!triangle[zone])
457 {
458 if (keep_left)
459 triangle[zone] = jedge ? -1 : 1;
460 else
461 triangle[zone] = jedge ? 1 : -1;
462 }
463 if (triangle[zone] > 0 ? !jedge : jedge)
464 goto bkwd;
465 }
466 else
467 {
468 if (keep_left)
469 goto bkwd;
470 }
471 }
472 /* bend forward (right along curve) */
473 keep_left = 1;
474 jedge = !jedge;
475 edge = p1 + (left > 0 ? left : 0);
476 {
477 long tmp = fwd;
478 fwd = -left;
479 left = tmp;
480 }
481 }
482 else if (z1 == z3)
483 {
484 bkwd:
485 /* bend backward (left along curve) */
486 keep_left = 0;
487 jedge = !jedge;
488 edge = p0 + (left > 0 ? left : 0);
489 {
490 long tmp = fwd;
491 fwd = left;
492 left = -tmp;
493 }
494 }
495 else
496 {
497 /* straight across to opposite edge */
498 edge += left;
499 }
500 /* after crossing zone, edge/left/fwd is oriented CCW relative to
501 * the next zone, assuming we will step there */
502
503 /* now that we've taken a step, check for the downstroke
504 * of a slit on the second pass (upstroke checked above)
505 * -- taking step first avoids a race condition */
506 if (pass2 && two_levels && !jedge)
507 {
508 if (left > 0)
509 {
510 if (data[edge] & SLIT_UP)
511 done = 6;
512 }
513 else
514 {
515 if (data[edge] & SLIT_DN)
516 done = 5;
517 }
518 }
519
520 if (!done)
521 {
522 /* finally, check if we are on a boundary */
523 if (data[edge] & (jedge ? J_BNDY : I_BNDY))
524 {
525 done = two_levels ? 2 : 4;
526 /* flip back into the zone that exists */
527 left = -left;
528 fwd = -fwd;
529 if (!pass2 && (edge != edge0 || left != left0))
530 {
531 Cdata start = data[edge] & START_MARK (left);
532 if (start)
533 {
534 data[edge] &= ~start;
535 site->count--;
536 }
537 }
538 }
539 }
540 }
541
542 site->edge = edge;
543 site->n = n;
544 site->left = left;
545 return done > 4 ? slit_cutter (site, done - 5, pass2) : done;
546 }
547
548 /* edge_walker assumes that the current edge is being drawn CCW
549 * around the current zone. Since only boundary edges are drawn
550 * and we always walk around with the filled region to the left,
551 * no edge is ever drawn CW. We attempt to advance to the next
552 * edge on this boundary, but if current second endpoint is not
553 * between the two contour levels, we exit back to zone_crosser.
554 * Note that we may wind up marking no points.
555 * -- edge_walker is never called for single level case */
556 static int
edge_walker(Csite * site,int pass2)557 edge_walker (Csite * site, int pass2)
558 {
559 Cdata * data = site->data;
560 long edge = site->edge;
561 long left = site->left;
562 long n = site->n;
563 long fwd = FORWARD (left, site->imax);
564 long p0 = POINT0 (edge, fwd);
565 long p1 = POINT1 (edge, fwd);
566 int jedge = IS_JEDGE (edge, left);
567 long edge0 = site->edge0;
568 long left0 = site->left0;
569 int level0 = site->level0 == 2;
570 int marked;
571
572 const double *x = pass2 ? site->x : 0;
573 const double *y = pass2 ? site->y : 0;
574 double *xcp = pass2 ? site->xcp : 0;
575 double *ycp = pass2 ? site->ycp : 0;
576
577 int z0, z1, heads_up = 0;
578
579 for (;;)
580 {
581 /* mark endpoint 0 only if value is 1 there, and this is a
582 * two level task */
583 z0 = data[p0] & Z_VALUE;
584 z1 = data[p1] & Z_VALUE;
585 marked = 0;
586 if (z0 == 1)
587 {
588 /* mark current boundary point */
589 if (pass2)
590 {
591 xcp[n] = x[p0];
592 ycp[n] = y[p0];
593 }
594 marked = 1;
595 }
596 else if (!n)
597 {
598 /* if this is the first point is not between the levels
599 * must do the job of the zone_crosser and mark the first cut here,
600 * so that it will be marked again by zone_crosser as it closes */
601 if (pass2)
602 {
603 double zcp = site->zlevel[(z0 != 0)];
604 zcp = (zcp - site->z[p0]) / (site->z[p1] - site->z[p0]);
605 xcp[n] = zcp * (x[p1] - x[p0]) + x[p0];
606 ycp[n] = zcp * (y[p1] - y[p0]) + y[p0];
607 }
608 marked = 1;
609 }
610 if (n)
611 {
612 /* check for closure */
613 if (level0 && edge == edge0 && left == left0)
614 {
615 site->edge = edge;
616 site->left = left;
617 site->n = n + marked;
618 /* if the curve is closing on a hole, need to make a downslit */
619 if (fwd < 0 && !(data[edge] & (jedge ? J_BNDY : I_BNDY)))
620 return slit_cutter (site, 0, pass2);
621 return 3;
622 }
623 else if (pass2)
624 {
625 if (heads_up || (fwd < 0 && (data[edge] & SLIT_DN)))
626 {
627 site->edge = edge;
628 site->left = left;
629 site->n = n + marked;
630 return slit_cutter (site, heads_up, pass2);
631 }
632 }
633 else
634 {
635 /* if this is not first point, clear start mark for this edge */
636 Cdata start = data[edge] & START_MARK (left);
637 if (start)
638 {
639 data[edge] &= ~start;
640 site->count--;
641 }
642 }
643 }
644 if (marked)
645 n++;
646
647 /* if next endpoint not between levels, need to exit to zone_crosser */
648 if (z1 != 1)
649 {
650 site->edge = edge;
651 site->left = left;
652 site->n = n;
653 return (z1 != 0); /* return level closest to p1 */
654 }
655
656 /* step to p1 and find next edge
657 * -- turn left if possible, else straight, else right
658 * -- check for upward slit beginning at same time */
659 edge = p1 + (left > 0 ? left : 0);
660 if (pass2 && jedge && fwd > 0 && (data[edge] & SLIT_UP))
661 {
662 jedge = !jedge;
663 heads_up = 1;
664 }
665 else if (data[edge] & (jedge ? I_BNDY : J_BNDY))
666 {
667 long tmp = fwd;
668 fwd = left;
669 left = -tmp;
670 jedge = !jedge;
671 }
672 else
673 {
674 edge = p1 + (fwd > 0 ? fwd : 0);
675 if (pass2 && !jedge && fwd > 0 && (data[edge] & SLIT_UP))
676 {
677 heads_up = 1;
678 }
679 else if (!(data[edge] & (jedge ? J_BNDY : I_BNDY)))
680 {
681 edge = p1 - (left < 0 ? left : 0);
682 jedge = !jedge;
683 {
684 long tmp = fwd;
685 fwd = -left;
686 left = tmp;
687 }
688 }
689 }
690 p0 = p1;
691 p1 = POINT1 (edge, fwd);
692 }
693 }
694
695 /* -- slit_cutter is never called for single level case */
696 static int
slit_cutter(Csite * site,int up,int pass2)697 slit_cutter (Csite * site, int up, int pass2)
698 {
699 Cdata * data = site->data;
700 long imax = site->imax;
701 long n = site->n;
702
703 const double *x = pass2 ? site->x : 0;
704 const double *y = pass2 ? site->y : 0;
705 double *xcp = pass2 ? site->xcp : 0;
706 double *ycp = pass2 ? site->ycp : 0;
707
708 if (up)
709 {
710 /* upward stroke of slit proceeds up left side of slit until
711 * it hits a boundary or a point not between the contour levels
712 * -- this never happens on the first pass */
713 long p1 = site->edge;
714 int z1;
715 for (;;)
716 {
717 z1 = data[p1] & Z_VALUE;
718 if (z1 != 1)
719 {
720 site->edge = p1;
721 site->left = -1;
722 site->n = n;
723 return (z1 != 0);
724 }
725 else if (data[p1] & J_BNDY)
726 {
727 /* this is very unusual case of closing on a mesh hole */
728 site->edge = p1;
729 site->left = -imax;
730 site->n = n;
731 return 2;
732 }
733 xcp[n] = x[p1];
734 ycp[n] = y[p1];
735 n++;
736 p1 += imax;
737 }
738
739 }
740 else
741 {
742 /* downward stroke proceeds down right side of slit until it
743 * hits a boundary or point not between the contour levels */
744 long p0 = site->edge;
745 int z0;
746 /* at beginning of first pass, mark first i-edge with SLIT_DN */
747 data[p0] |= SLIT_DN;
748 p0 -= imax;
749 for (;;)
750 {
751 z0 = data[p0] & Z_VALUE;
752 if (!pass2)
753 {
754 if (z0 != 1 || (data[p0] & I_BNDY) || (data[p0 + 1] & J_BNDY))
755 {
756 /* at end of first pass, mark final i-edge with SLIT_UP */
757 data[p0 + imax] |= SLIT_UP;
758 /* one extra count for splicing at outer curve */
759 site->n = n + 1;
760 return 4; /* return same special value as for OPEN_END */
761 }
762 }
763 else
764 {
765 if (z0 != 1)
766 {
767 site->edge = p0 + imax;
768 site->left = 1;
769 site->n = n;
770 return (z0 != 0);
771 }
772 else if (data[p0 + 1] & J_BNDY)
773 {
774 site->edge = p0 + 1;
775 site->left = imax;
776 site->n = n;
777 return 2;
778 }
779 else if (data[p0] & I_BNDY)
780 {
781 site->edge = p0;
782 site->left = 1;
783 site->n = n;
784 return 2;
785 }
786 }
787 if (pass2)
788 {
789 xcp[n] = x[p0];
790 ycp[n] = y[p0];
791 n++;
792 }
793 else
794 {
795 /* on first pass need to count for upstroke as well */
796 n += 2;
797 }
798 p0 -= imax;
799 }
800 }
801 }
802
803 /* ------------------------------------------------------------------------ */
804
805 /* curve_tracer finds the next starting point, then traces the curve,
806 * returning the number of points on this curve
807 * -- in a two level trace, the return value is negative on the
808 * first pass if the curve closed on a hole
809 * -- in a single level trace, the return value is negative on the
810 * first pass if the curve is an incomplete open curve
811 * -- a return value of 0 indicates no more curves */
812 static long
curve_tracer(Csite * site,int pass2)813 curve_tracer (Csite * site, int pass2)
814 {
815 Cdata * data = site->data;
816 long imax = site->imax;
817 long edge0 = site->edge0;
818 long left0 = site->left0;
819 long edge00 = site->edge00;
820 int two_levels = site->zlevel[1] > site->zlevel[0];
821 int level, level0, mark_row;
822 long n;
823
824 /* it is possible for a single i-edge to serve as two actual start
825 * points, one to the right and one to the left
826 * -- for the two level case, this happens on the first pass for
827 * a doubly cut edge, or on a chunking boundary
828 * -- for single level case, this is impossible, but a similar
829 * situation involving open curves is handled below
830 * a second two start possibility is when the edge0 zone does not
831 * exist and both the i-edge and j-edge boundaries are cut
832 * yet another possibility is three start points at a junction
833 * of chunk cuts
834 * -- sigh, several other rare possibilities,
835 * allow for general case, just go in order i1, i0, j1, j0 */
836 int two_starts;
837 /* printf("curve_tracer pass %d\n", pass2); */
838 /* print_Csite(site); */
839 if (left0 == 1)
840 two_starts = data[edge0] & (I0_START | J1_START | J0_START);
841 else if (left0 == -1)
842 two_starts = data[edge0] & (J1_START | J0_START);
843 else if (left0 == imax)
844 two_starts = data[edge0] & J0_START;
845 else
846 two_starts = 0;
847
848 if (pass2 || edge0 == 0)
849 {
850 /* zip up to row marked on first pass (or by data_init if edge0==0)
851 * -- but not for double start case */
852 if (!two_starts)
853 {
854 /* final start point marked by ALL_DONE marker */
855 int first = (edge0 == 0 && !pass2);
856 long e0 = edge0;
857 if (data[edge0] & ALL_DONE)
858 return 0;
859 while (!(data[edge0] & START_ROW))
860 edge0 += imax;
861 if (e0 == edge0)
862 edge0++; /* two starts handled specially */
863 if (first)
864 /* if this is the very first start point, we want to remove
865 * the START_ROW marker placed by data_init */
866 data[edge0 - edge0 % imax] &= ~START_ROW;
867 }
868
869 }
870 else
871 {
872 /* first pass ends when all potential start points visited */
873 if (site->count <= 0)
874 {
875 /* place ALL_DONE marker for second pass */
876 data[edge00] |= ALL_DONE;
877 /* reset initial site for second pass */
878 site->edge0 = site->edge00 = site->left0 = 0;
879 return 0;
880 }
881 if (!two_starts)
882 edge0++;
883 }
884
885 if (two_starts)
886 {
887 /* trace second curve with this start immediately */
888 if (left0 == 1 && (data[edge0] & I0_START))
889 {
890 left0 = -1;
891 level = (data[edge0] & I_BNDY) ? 2 : 0;
892 }
893 else if ((left0 == 1 || left0 == -1) && (data[edge0] & J1_START))
894 {
895 left0 = imax;
896 level = 2;
897 }
898 else
899 {
900 left0 = -imax;
901 level = 2;
902 }
903
904 }
905 else
906 {
907 /* usual case is to scan for next start marker
908 * -- on second pass, this is at most one row of mesh, but first
909 * pass hits nearly every point of the mesh, since it can't
910 * know in advance which potential start marks removed */
911 while (!(data[edge0] & ANY_START))
912 edge0++;
913
914 if (data[edge0] & I1_START)
915 left0 = 1;
916 else if (data[edge0] & I0_START)
917 left0 = -1;
918 else if (data[edge0] & J1_START)
919 left0 = imax;
920 else /*data[edge0]&J0_START */
921 left0 = -imax;
922
923 if (data[edge0] & (I1_START | I0_START))
924 level = (data[edge0] & I_BNDY) ? 2 : 0;
925 else
926 level = 2;
927 }
928
929 /* this start marker will not be unmarked, but it has been visited */
930 if (!pass2)
931 site->count--;
932
933 /* if this curve starts on a non-boundary i-edge, we need to
934 * determine the level */
935 if (!level && two_levels)
936 level = left0 > 0 ?
937 ((data[edge0 - imax] & Z_VALUE) !=
938 0) : ((data[edge0] & Z_VALUE) != 0);
939
940 /* initialize site for this curve */
941 site->edge = site->edge0 = edge0;
942 site->left = site->left0 = left0;
943 site->level0 = level0 = level; /* for open curve detection only */
944
945 /* single level case just uses zone_crosser */
946 if (!two_levels)
947 level = 0;
948
949 /* to generate the curve, alternate between zone_crosser and
950 * edge_walker until closure or first call to edge_walker in
951 * single level case */
952 site->n = 0;
953 for (;;)
954 {
955 if (level < 2)
956 level = zone_crosser (site, level, pass2);
957 else if (level < 3)
958 level = edge_walker (site, pass2);
959 else
960 break;
961 }
962 n = site->n;
963
964 /* single level case may have ended at a boundary rather than closing
965 * -- need to recognize this case here in order to place the
966 * OPEN_END mark for zone_crosser, remove this start marker,
967 * and be sure not to make a START_ROW mark for this case
968 * two level case may close with slit_cutter, in which case start
969 * must also be removed and no START_ROW mark made
970 * -- change sign of return n to inform caller */
971 if (!pass2 && level > 3 && (two_levels || level0 == 0))
972 {
973 if (!two_levels)
974 data[edge0] |= OPEN_END;
975 data[edge0] &= ~(left0 > 0 ? I1_START : I0_START);
976 mark_row = 0; /* do not mark START_ROW */
977 n = -n;
978 }
979 else
980 {
981 if (two_levels)
982 mark_row = !two_starts;
983 else
984 mark_row = 1;
985 }
986
987 /* on first pass, must apply START_ROW mark in column above previous
988 * start marker
989 * -- but skip if we just did second of two start case */
990 if (!pass2 && mark_row)
991 {
992 data[edge0 - (edge0 - edge00) % imax] |= START_ROW;
993 site->edge00 = edge0;
994 }
995
996 return n;
997 }
998
999 /* ------------------------------------------------------------------------ */
1000
1001 /* The sole function of the "region" argument is to specify the
1002 value in Csite.reg that denotes a missing zone. We always
1003 use zero.
1004 */
1005
1006 static void
data_init(Csite * site,int region,long nchunk)1007 data_init (Csite * site, int region, long nchunk)
1008 {
1009 Cdata * data = site->data;
1010 long imax = site->imax;
1011 long jmax = site->jmax;
1012 long ijmax = imax * jmax;
1013 const double *z = site->z;
1014 double zlev0 = site->zlevel[0];
1015 double zlev1 = site->zlevel[1];
1016 int two_levels = zlev1 > zlev0;
1017 char *reg = site->reg;
1018 long count = 0;
1019 int started = 0;
1020 int ibndy, jbndy, i_was_chunk;
1021
1022 long icsize = imax - 1;
1023 long jcsize = jmax - 1;
1024 long ichunk, jchunk, irem, jrem, i, j, ij;
1025
1026 if (nchunk && two_levels)
1027 {
1028 /* figure out chunk sizes
1029 * -- input nchunk is square root of maximum allowed zones per chunk
1030 * -- start points for single level case are wrong, so don't try it */
1031 long inum = (nchunk * nchunk) / (jmax - 1);
1032 long jnum = (nchunk * nchunk) / (imax - 1);
1033 if (inum < nchunk)
1034 inum = nchunk;
1035 if (jnum < nchunk)
1036 jnum = nchunk;
1037 /* ijnum= actual number of chunks,
1038 * ijrem= number of those chunks needing one more zone (ijcsize+1) */
1039 inum = (imax - 2) / inum + 1;
1040 icsize = (imax - 1) / inum;
1041 irem = (imax - 1) % inum;
1042 jnum = (jmax - 2) / jnum + 1;
1043 jcsize = (jmax - 1) / jnum;
1044 jrem = (jmax - 1) % jnum;
1045 /* convert ijrem into value of i or j at which to begin adding an
1046 * extra zone */
1047 irem = (inum - irem) * icsize;
1048 jrem = (jnum - jrem) * jcsize;
1049 }
1050 else
1051 {
1052 irem = imax;
1053 jrem = jmax;
1054 }
1055
1056 /* do everything in a single pass through the data array to
1057 * minimize cache faulting (z, reg, and data are potentially
1058 * very large arrays)
1059 * access to the z and reg arrays is strictly sequential,
1060 * but we need two rows (+-imax) of the data array at a time */
1061 if (z[0] > zlev0)
1062 data[0] = (two_levels && z[0] > zlev1) ? 2 : 1;
1063 else
1064 data[0] = 0;
1065 jchunk = 0;
1066 for (j = ij = 0; j < jmax; j++)
1067 {
1068 ichunk = i_was_chunk = 0;
1069 for (i = 0; i < imax; i++, ij++)
1070 {
1071 /* transfer zonal existence from reg to data array
1072 * -- get these for next row so we can figure existence of
1073 * points and j-edges for this row */
1074 data[ij + imax + 1] = 0;
1075 if (reg)
1076 {
1077 if (region ? (reg[ij + imax + 1] == region)
1078 : (reg[ij + imax + 1] != 0))
1079 data[ij + imax + 1] = ZONE_EX;
1080 }
1081 else
1082 {
1083 if (i < imax - 1 && j < jmax - 1)
1084 data[ij + imax + 1] = ZONE_EX;
1085 }
1086
1087 /* translate z values to 0, 1, 2 flags */
1088 if (ij < imax)
1089 data[ij + 1] = 0;
1090 if (ij < ijmax - 1 && z[ij + 1] > zlev0)
1091 data[ij + 1] |= (two_levels && z[ij + 1] > zlev1) ? 2 : 1;
1092
1093 /* apply edge boundary marks */
1094 ibndy = i == ichunk
1095 || (data[ij] & ZONE_EX) != (data[ij + 1] & ZONE_EX);
1096 jbndy = j == jchunk
1097 || (data[ij] & ZONE_EX) != (data[ij + imax] & ZONE_EX);
1098 if (ibndy)
1099 data[ij] |= I_BNDY;
1100 if (jbndy)
1101 data[ij] |= J_BNDY;
1102
1103 /* apply i-edge start marks
1104 * -- i-edges are only marked when actually cut
1105 * -- no mark is necessary if one of the j-edges which share
1106 * the lower endpoint is also cut
1107 * -- no I0 mark necessary unless filled region below some cut,
1108 * no I1 mark necessary unless filled region above some cut */
1109 if (j)
1110 {
1111 int v0 = (data[ij] & Z_VALUE);
1112 int vb = (data[ij - imax] & Z_VALUE);
1113 if (v0 != vb)
1114 { /* i-edge is cut */
1115 if (ibndy)
1116 {
1117 if (data[ij] & ZONE_EX)
1118 {
1119 data[ij] |= I0_START;
1120 count++;
1121 }
1122 if (data[ij + 1] & ZONE_EX)
1123 {
1124 data[ij] |= I1_START;
1125 count++;
1126 }
1127 }
1128 else
1129 {
1130 int va = (data[ij - 1] & Z_VALUE);
1131 int vc = (data[ij + 1] & Z_VALUE);
1132 int vd = (data[ij - imax + 1] & Z_VALUE);
1133 if (v0 != 1 && va != v0
1134 && (vc != v0 || vd != v0) && (data[ij] & ZONE_EX))
1135 {
1136 data[ij] |= I0_START;
1137 count++;
1138 }
1139 if (vb != 1 && va == vb
1140 && (vc == vb || vd == vb)
1141 && (data[ij + 1] & ZONE_EX))
1142 {
1143 data[ij] |= I1_START;
1144 count++;
1145 }
1146 }
1147 }
1148 }
1149
1150 /* apply j-edge start marks
1151 * -- j-edges are only marked when they are boundaries
1152 * -- all cut boundary edges marked
1153 * -- for two level case, a few uncut edges must be marked
1154 */
1155 if (i && jbndy)
1156 {
1157 int v0 = (data[ij] & Z_VALUE);
1158 int vb = (data[ij - 1] & Z_VALUE);
1159 if (v0 != vb)
1160 {
1161 if (data[ij] & ZONE_EX)
1162 {
1163 data[ij] |= J0_START;
1164 count++;
1165 }
1166 if (data[ij + imax] & ZONE_EX)
1167 {
1168 data[ij] |= J1_START;
1169 count++;
1170 }
1171 }
1172 else if (two_levels && v0 == 1)
1173 {
1174 if (data[ij + imax] & ZONE_EX)
1175 {
1176 if (i_was_chunk || !(data[ij + imax - 1] & ZONE_EX))
1177 {
1178 /* lower left is a drawn part of boundary */
1179 data[ij] |= J1_START;
1180 count++;
1181 }
1182 }
1183 else if (data[ij] & ZONE_EX)
1184 {
1185 if (data[ij + imax - 1] & ZONE_EX)
1186 {
1187 /* weird case of open hole at lower left */
1188 data[ij] |= J0_START;
1189 count++;
1190 }
1191 }
1192 }
1193 }
1194
1195 i_was_chunk = (i == ichunk);
1196 if (i_was_chunk)
1197 ichunk += icsize + (ichunk >= irem);
1198 }
1199
1200 if (j == jchunk)
1201 jchunk += jcsize + (jchunk >= jrem);
1202
1203 /* place first START_ROW marker */
1204 if (count && !started)
1205 {
1206 data[ij - imax] |= START_ROW;
1207 started = 1;
1208 }
1209 }
1210
1211 /* place immediate stop mark if nothing found */
1212 if (!count)
1213 data[0] |= ALL_DONE;
1214
1215 /* initialize site */
1216 site->edge0 = site->edge00 = site->edge = 0;
1217 site->left0 = site->left = 0;
1218 site->n = 0;
1219 site->count = count;
1220 }
1221
1222 /* ------------------------------------------------------------------------
1223 Original (slightly modified) core contour generation routines are above;
1224 below are new routines for interfacing to mpl.
1225
1226 ------------------------------------------------------------------------ */
1227
1228 /* Note: index order gets switched in the Python interface;
1229 python Z[i,j] -> C z[j,i]
1230 so if the array has shape Mi, Nj in python,
1231 we have iMax = Nj, jMax = Mi in gcntr.c.
1232 On the Python side: Ny, Nx = shape(z),
1233 so in C, the x-dimension is the first index, the y-dimension
1234 the second.
1235 */
1236
1237 /* reg should have the same dimensions as data, which
1238 has an extra iMax + 1 points relative to Z.
1239 It differs from mask in being the opposite (True
1240 where a region exists, versus the mask, which is True
1241 where a data point is bad), and in that it marks
1242 zones, not points. All four zones sharing a bad
1243 point must be marked as not existing.
1244 */
1245 static void
mask_zones(long iMax,long jMax,char * mask,char * reg)1246 mask_zones (long iMax, long jMax, char *mask, char *reg)
1247 {
1248 long i, j, ij;
1249 long nreg = iMax * jMax + iMax + 1;
1250
1251 for (ij = iMax+1; ij < iMax*jMax; ij++)
1252 {
1253 reg[ij] = 1;
1254 }
1255
1256 ij = 0;
1257 for (j = 0; j < jMax; j++)
1258 {
1259 for (i = 0; i < iMax; i++, ij++)
1260 {
1261 if (i == 0 || j == 0) reg[ij] = 0;
1262 if (mask[ij] != 0)
1263 {
1264 reg[ij] = 0;
1265 reg[ij + 1] = 0;
1266 reg[ij + iMax] = 0;
1267 reg[ij + iMax + 1] = 0;
1268 }
1269 }
1270 }
1271 for (; ij < nreg; ij++)
1272 {
1273 reg[ij] = 0;
1274 }
1275 }
1276
1277 static Csite *
cntr_new(void)1278 cntr_new(void)
1279 {
1280 Csite *site;
1281 site = (Csite *) PyMem_Malloc(sizeof(Csite));
1282 if (site == NULL) return NULL;
1283 site->data = NULL;
1284 site->reg = NULL;
1285 site->triangle = NULL;
1286 site->xcp = NULL;
1287 site->ycp = NULL;
1288 site->x = NULL;
1289 site->y = NULL;
1290 site->z = NULL;
1291 return site;
1292 }
1293
1294 static int
cntr_init(Csite * site,long iMax,long jMax,double * x,double * y,double * z,char * mask)1295 cntr_init(Csite *site, long iMax, long jMax, double *x, double *y,
1296 double *z, char *mask)
1297 {
1298 long ijmax = iMax * jMax;
1299 long nreg = iMax * jMax + iMax + 1;
1300 long i;
1301
1302 site->imax = iMax;
1303 site->jmax = jMax;
1304 site->data = (Cdata *) PyMem_Malloc(sizeof(Cdata) * nreg);
1305 if (site->data == NULL)
1306 {
1307 PyMem_Free(site);
1308 return -1;
1309 }
1310 site->triangle = (short *) PyMem_Malloc(sizeof(short) * ijmax);
1311 if (site->triangle == NULL)
1312 {
1313 PyMem_Free(site->data);
1314 PyMem_Free(site);
1315 return -1;
1316 }
1317 for (i = 0; i < ijmax; i++) site->triangle[i] = 0;
1318 site->reg = NULL;
1319 if (mask != NULL)
1320 {
1321 site->reg = (char *) PyMem_Malloc(sizeof(char) * nreg);
1322 if (site->reg == NULL)
1323 {
1324 PyMem_Free(site->triangle);
1325 PyMem_Free(site->data);
1326 PyMem_Free(site);
1327 return -1;
1328 }
1329 mask_zones(iMax, jMax, mask, site->reg);
1330 }
1331 /* I don't think we need to initialize site->data. */
1332 site->x = x;
1333 site->y = y;
1334 site->z = z;
1335 site->xcp = NULL;
1336 site->ycp = NULL;
1337 return 0;
1338 }
1339
cntr_del(Csite * site)1340 static void cntr_del(Csite *site)
1341 {
1342 PyMem_Free(site->triangle);
1343 PyMem_Free(site->reg);
1344 PyMem_Free(site->data);
1345 PyMem_Free(site);
1346 site = NULL;
1347 }
1348
1349 /* Build a list of lists of points, where each point is an (x,y)
1350 tuple.
1351 */
1352 static PyObject *
build_cntr_list_p(long * np,double * xp,double * yp,int nparts,long ntotal)1353 build_cntr_list_p(long *np, double *xp, double *yp, int nparts, long ntotal)
1354 {
1355 PyObject *point, *contourList, *all_contours;
1356 int start = 0, end = 0;
1357 int i, j, k;
1358
1359 all_contours = PyList_New(nparts);
1360 if (all_contours == NULL) return NULL;
1361
1362 for (i = 0; i < nparts; i++)
1363 {
1364 start = end;
1365 end += np[i];
1366 contourList = PyList_New(np[i]);
1367 if (contourList == NULL) goto error;
1368 for (k = 0, j = start; j < end; j++, k++)
1369 {
1370 point = Py_BuildValue("(dd)", xp[j], yp[j]);
1371 if (PyList_SetItem(contourList, k, point)) goto error;
1372 }
1373 if (PyList_SetItem(all_contours, i, contourList)) goto error;
1374 }
1375 return all_contours;
1376
1377 error:
1378 Py_XDECREF(all_contours);
1379 return NULL;
1380 }
1381
1382 /* Build a list of tuples (X, Y), where X and Y are 1-D arrays. */
1383 #if 0
1384 static PyObject *
1385 build_cntr_list_v(long *np, double *xp, double *yp, int nparts, long ntotal)
1386 {
1387 PyObject *point, *all_contours;
1388 PyArrayObject *xv, *yv;
1389 npy_intp dims[1];
1390 int i;
1391 long j, k;
1392
1393 all_contours = PyList_New(nparts);
1394
1395 k = 0;
1396 for (i = 0; i < nparts; i++)
1397 {
1398 dims[0] = np[i];
1399 xv = (PyArrayObject *) PyArray_SimpleNew(1, dims, NPY_DOUBLE);
1400 yv = (PyArrayObject *) PyArray_SimpleNew(1, dims, NPY_DOUBLE);
1401 if (xv == NULL || yv == NULL) goto error;
1402 for (j = 0; j < dims[0]; j++)
1403 {
1404 ((double *)xv->data)[j] = xp[k];
1405 ((double *)yv->data)[j] = yp[k];
1406 k++;
1407 }
1408 point = Py_BuildValue("(NN)", xv, yv);
1409 /* "O" increments ref count; "N" does not. */
1410 if (PyList_SetItem(all_contours, i, point)) goto error;
1411 }
1412 return all_contours;
1413
1414 error:
1415 Py_XDECREF(all_contours);
1416 return NULL;
1417 }
1418 #endif
1419
1420 /* Build a list of XY 2-D arrays, shape (N,2) */
1421 static PyObject *
build_cntr_list_v2(long * np,double * xp,double * yp,int nparts,long ntotal)1422 build_cntr_list_v2(long *np, double *xp, double *yp, int nparts, long ntotal)
1423 {
1424 PyObject *all_contours;
1425 PyArrayObject *xyv;
1426 npy_intp dims[2];
1427 int i;
1428 long j, k;
1429
1430 all_contours = PyList_New(nparts);
1431 if (all_contours == NULL) return NULL;
1432
1433 k = 0;
1434 for (i = 0; i < nparts; i++)
1435 {
1436 double *data;
1437 dims[0] = np[i];
1438 dims[1] = 2;
1439 xyv = (PyArrayObject *) PyArray_SimpleNew(2, dims, NPY_DOUBLE);
1440 if (xyv == NULL) goto error;
1441 data = (double*)PyArray_DATA(xyv);
1442 for (j = 0; j < dims[0]; j++)
1443 {
1444 data[2*j] = xp[k];
1445 data[2*j+1] = yp[k];
1446 k++;
1447 }
1448 if (PyList_SetItem(all_contours, i, (PyObject *)xyv)) goto error;
1449 }
1450 return all_contours;
1451
1452 error:
1453 Py_XDECREF(all_contours);
1454 return NULL;
1455 }
1456
1457
1458
1459 /* cntr_trace is called once per contour level or level pair.
1460 If nlevels is 1, a set of contour lines will be returned; if nlevels
1461 is 2, the set of polygons bounded by the levels will be returned.
1462 If points is True, the lines will be returned as a list of list
1463 of points; otherwise, as a list of tuples of vectors.
1464 */
1465
1466 static PyObject *
cntr_trace(Csite * site,double levels[],int nlevels,int points,long nchunk)1467 cntr_trace(Csite *site, double levels[], int nlevels, int points, long nchunk)
1468 {
1469 PyObject *c_list;
1470 double *xp0;
1471 double *yp0;
1472 long *nseg0;
1473 int iseg;
1474
1475 /* long nchunk = 30; was hardwired */
1476 long n;
1477 long nparts = 0;
1478 long ntotal = 0;
1479 long nparts2 = 0;
1480 long ntotal2 = 0;
1481
1482 site->zlevel[0] = levels[0];
1483 site->zlevel[1] = levels[0];
1484 if (nlevels == 2)
1485 {
1486 site->zlevel[1] = levels[1];
1487 }
1488 site->n = site->count = 0;
1489 data_init (site, 0, nchunk);
1490
1491 /* make first pass to compute required sizes for second pass */
1492 for (;;)
1493 {
1494 n = curve_tracer (site, 0);
1495
1496 if (!n)
1497 break;
1498 if (n > 0)
1499 {
1500 nparts++;
1501 ntotal += n;
1502 }
1503 else
1504 {
1505 ntotal -= n;
1506 }
1507 }
1508 xp0 = (double *) PyMem_Malloc(ntotal * sizeof(double));
1509 yp0 = (double *) PyMem_Malloc(ntotal * sizeof(double));
1510 nseg0 = (long *) PyMem_Malloc(nparts * sizeof(long));
1511 if (xp0 == NULL || yp0 == NULL || nseg0 == NULL)
1512 {
1513 PyErr_NoMemory();
1514 goto error;
1515 }
1516
1517 /* second pass */
1518 site->xcp = xp0;
1519 site->ycp = yp0;
1520 iseg = 0;
1521 for (;;iseg++)
1522 {
1523 n = curve_tracer (site, 1);
1524 if (ntotal2 + n > ntotal)
1525 {
1526 PyErr_SetString(PyExc_RuntimeError,
1527 "curve_tracer: ntotal2, pass 2 exceeds ntotal, pass 1");
1528 goto error;
1529 }
1530 if (n == 0)
1531 break;
1532 if (n > 0)
1533 {
1534 /* could add array bounds checking */
1535 nseg0[iseg] = n;
1536 site->xcp += n;
1537 site->ycp += n;
1538 ntotal2 += n;
1539 nparts2++;
1540 }
1541 else
1542 {
1543 PyErr_SetString(PyExc_RuntimeError,
1544 "Negative n from curve_tracer in pass 2");
1545 goto error;
1546 }
1547 }
1548
1549
1550 if (points)
1551 {
1552 c_list = build_cntr_list_p(nseg0, xp0, yp0, nparts, ntotal);
1553 }
1554 else
1555 {
1556 c_list = build_cntr_list_v2(nseg0, xp0, yp0, nparts, ntotal);
1557 }
1558 PyMem_Free(xp0); PyMem_Free(yp0); PyMem_Free(nseg0);
1559 site->xcp = NULL; site->ycp = NULL;
1560 return c_list;
1561
1562 error:
1563 PyMem_Free(xp0); PyMem_Free(yp0); PyMem_Free(nseg0);
1564 site->xcp = NULL; site->ycp = NULL;
1565 return NULL;
1566 }
1567
1568 /******* Make an extension type. Based on the tutorial.************/
1569
1570 /* site points to the data arrays in the arrays pointed to
1571 by xpa, ypa, zpa, and mpa, so we include them in the
1572 structure so we can ensure they are not deleted until
1573 we have finished using them.
1574 */
1575 typedef struct {
1576 PyObject_HEAD
1577 PyArrayObject *xpa, *ypa, *zpa, *mpa;
1578 Csite *site;
1579 } Cntr;
1580
1581
1582 static int
Cntr_clear(Cntr * self)1583 Cntr_clear(Cntr* self)
1584 {
1585 PyArrayObject *tmp;
1586
1587 cntr_del(self->site);
1588
1589 tmp = self->xpa;
1590 self->xpa = NULL;
1591 Py_XDECREF(tmp);
1592
1593 tmp = self->ypa;
1594 self->ypa = NULL;
1595 Py_XDECREF(tmp);
1596
1597 tmp = self->zpa;
1598 self->zpa = NULL;
1599 Py_XDECREF(tmp);
1600
1601 tmp = self->mpa;
1602 self->mpa = NULL;
1603 Py_XDECREF(tmp);
1604 return 0;
1605 }
1606
1607 static void
Cntr_dealloc(Cntr * self)1608 Cntr_dealloc(Cntr* self)
1609 {
1610 Cntr_clear(self);
1611 Py_TYPE(self)->tp_free((PyObject*)self);
1612 }
1613
1614 static PyObject *
Cntr_new(PyTypeObject * type,PyObject * args,PyObject * kwds)1615 Cntr_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
1616 {
1617 Cntr *self;
1618
1619 self = (Cntr *)type->tp_alloc(type, 0);
1620 if (self != NULL)
1621 {
1622 self->site = cntr_new();
1623 if (self->site == NULL)
1624 {
1625 PyErr_SetString(PyExc_MemoryError,
1626 "Memory allocation failed in cntr_new.");
1627 Py_XDECREF(self);
1628 return NULL;
1629 }
1630 self->xpa = NULL;
1631 self->ypa = NULL;
1632 self->zpa = NULL;
1633 self->mpa = NULL;
1634 }
1635
1636 return (PyObject *)self;
1637 }
1638
1639 static int
Cntr_init(Cntr * self,PyObject * args,PyObject * kwds)1640 Cntr_init(Cntr *self, PyObject *args, PyObject *kwds)
1641 {
1642 static char *kwlist[] = {"x", "y", "z", "mask", NULL};
1643 PyObject *xarg, *yarg, *zarg, *marg;
1644 PyArrayObject *xpa, *ypa, *zpa, *mpa;
1645 long iMax, jMax;
1646 char *mask;
1647
1648 marg = NULL;
1649
1650 if (! PyArg_ParseTupleAndKeywords(args, kwds, "OOO|O", kwlist,
1651 &xarg, &yarg, &zarg, &marg))
1652 return -1;
1653 if (marg == Py_None)
1654 marg = NULL;
1655
1656 if (!PyArray_Check(xarg) || !PyArray_Check(yarg) ||
1657 !PyArray_Check(zarg) || (marg && !PyArray_Check(marg)))
1658 {
1659 PyErr_SetString(PyExc_TypeError,
1660 "Arguments x, y, z, (optional) mask must be arrays.");
1661 return -1;
1662 }
1663
1664 xpa = (PyArrayObject *) PyArray_ContiguousFromObject(xarg,
1665 NPY_DOUBLE, 2, 2);
1666 ypa = (PyArrayObject *) PyArray_ContiguousFromObject(yarg,
1667 NPY_DOUBLE,
1668 2, 2);
1669 zpa = (PyArrayObject *) PyArray_ContiguousFromObject(zarg, NPY_DOUBLE,
1670 2, 2);
1671 if (marg)
1672 mpa = (PyArrayObject *) PyArray_ContiguousFromObject(marg,
1673 NPY_BYTE,
1674 2, 2);
1675 else
1676 mpa = NULL;
1677
1678 if (xpa == NULL || ypa == NULL || zpa == NULL || (marg && mpa == NULL))
1679 {
1680 PyErr_SetString(PyExc_ValueError,
1681 "Arguments x, y, z, mask (if present) must be 2D arrays.");
1682 goto error;
1683 }
1684 iMax = PyArray_DIMS(zpa)[1];
1685 jMax = PyArray_DIMS(zpa)[0];
1686 if (PyArray_DIMS(xpa)[0] != jMax || PyArray_DIMS(xpa)[1] != iMax ||
1687 PyArray_DIMS(ypa)[0] != jMax || PyArray_DIMS(ypa)[1] != iMax ||
1688 (mpa && (PyArray_DIMS(mpa)[0] != jMax || PyArray_DIMS(mpa)[1] != iMax)))
1689 {
1690 PyErr_SetString(PyExc_ValueError,
1691 "Arguments x, y, z, mask (if present)"
1692 " must have the same dimensions.");
1693 goto error;
1694 }
1695 if (mpa) mask = PyArray_DATA(mpa);
1696 else mask = NULL;
1697 if ( cntr_init(self->site, iMax, jMax, (double *)PyArray_DATA(xpa),
1698 (double *)PyArray_DATA(ypa),
1699 (double *)PyArray_DATA(zpa), mask))
1700 {
1701 PyErr_SetString(PyExc_MemoryError,
1702 "Memory allocation failure in cntr_init");
1703 goto error;
1704 }
1705 self->xpa = xpa;
1706 self->ypa = ypa;
1707 self->zpa = zpa;
1708 self->mpa = mpa;
1709 return 0;
1710
1711 error:
1712 Py_XDECREF(xpa);
1713 Py_XDECREF(ypa);
1714 Py_XDECREF(zpa);
1715 Py_XDECREF(mpa);
1716 return -1;
1717 }
1718
1719 static PyObject *
Cntr_trace(Cntr * self,PyObject * args,PyObject * kwds)1720 Cntr_trace(Cntr *self, PyObject *args, PyObject *kwds)
1721 {
1722 double levels[2] = {0.0, -1e100};
1723 int nlevels = 2;
1724 int points = 0;
1725 long nchunk = 0L;
1726 static char *kwlist[] = {"level0", "level1", "points", "nchunk", NULL};
1727
1728 if (! PyArg_ParseTupleAndKeywords(args, kwds, "d|dil", kwlist,
1729 levels, levels+1, &points, &nchunk))
1730 {
1731 return NULL;
1732 }
1733 if (levels[1] == -1e100 || levels[1] <= levels[0])
1734 nlevels = 1;
1735 return cntr_trace(self->site, levels, nlevels, points, nchunk);
1736 }
1737
1738 static PyMethodDef Cntr_methods[] = {
1739 {"trace", (PyCFunction)Cntr_trace, METH_VARARGS | METH_KEYWORDS,
1740 "Return a list of contour line segments or polygons.\n\n"
1741 " Required argument: level0, a contour level\n"
1742 " Optional argument: level1; if given, and if level1 > level0,\n"
1743 " then the contours will be polygons surrounding areas between\n"
1744 " the levels.\n"
1745 " Optional argument: points; if 0 (default), return a list of\n"
1746 " vector pairs; otherwise, return a list of lists of points.\n"
1747 " Optional argument: nchunk; approximate number of grid points\n"
1748 " per chunk. 0 (default) for no chunking.\n"
1749 },
1750 {NULL} /* Sentinel */
1751 };
1752
1753 static PyTypeObject CntrType = {
1754 PyVarObject_HEAD_INIT(NULL, 0)
1755 "_nc_cntr.Cntr", /*tp_name*/
1756 sizeof(Cntr), /*tp_basicsize*/
1757 0, /*tp_itemsize*/
1758 (destructor)Cntr_dealloc, /*tp_dealloc*/
1759 0, /*tp_print*/
1760 0, /*tp_getattr*/
1761 0, /*tp_setattr*/
1762 0, /*tp_compare*/
1763 0, /*tp_repr*/
1764 0, /*tp_as_number*/
1765 0, /*tp_as_sequence*/
1766 0, /*tp_as_mapping*/
1767 0, /*tp_hash */
1768 0, /*tp_call*/
1769 0, /*tp_str*/
1770 0, /*tp_getattro*/
1771 0, /*tp_setattro*/
1772 0, /*tp_as_buffer*/
1773 Py_TPFLAGS_DEFAULT, /*tp_flags*/
1774 "Contour engine", /* tp_doc */
1775 0, /* tp_traverse */
1776 (inquiry)Cntr_clear, /* tp_clear */
1777 0, /* tp_richcompare */
1778 0, /* tp_weaklistoffset */
1779 0, /* tp_iter */
1780 0, /* tp_iternext */
1781 Cntr_methods, /* tp_methods */
1782 0, /* tp_members */
1783 0, /* tp_getset */
1784 0, /* tp_base */
1785 0, /* tp_dict */
1786 0, /* tp_descr_get */
1787 0, /* tp_descr_set */
1788 0, /* tp_dictoffset */
1789 (initproc)Cntr_init, /* tp_init */
1790 0, /* tp_alloc */
1791 Cntr_new, /* tp_new */
1792 };
1793
1794 static PyMethodDef module_methods[] = {
1795 {NULL} /* Sentinel */
1796 };
1797
1798 /* see http://python3porting.com/cextensions.html */
1799 #if PY_MAJOR_VERSION >= 3
1800 # define MOD_ERROR_VAL NULL
1801 # define MOD_SUCCESS_VAL(val) val
1802 # define MOD_INIT(name) PyMODINIT_FUNC PyInit_##name(void)
1803 # define MOD_DEF(ob, name, doc, methods) \
1804 { \
1805 static struct PyModuleDef moduledef = { \
1806 PyModuleDef_HEAD_INIT, name, doc, -1, methods, }; \
1807 ob = PyModule_Create(&moduledef); \
1808 }
1809 #else
1810 # define MOD_ERROR_VAL
1811 # define MOD_SUCCESS_VAL(val)
1812 # define MOD_INIT(name) void init##name(void)
1813 # define MOD_DEF(ob, name, doc, methods) \
1814 ob = Py_InitModule3(name, methods, doc)
1815 #endif
1816
1817 /* returns different type depending on python version */
1818 #if PY_MAJOR_VERSION >= 3
donumpyinit(void)1819 static void* donumpyinit(void) { import_array(); return NULL; }
1820 #else
donumpyinit(void)1821 static void donumpyinit(void) { import_array(); }
1822 #endif
1823
MOD_INIT(_nc_cntr)1824 MOD_INIT(_nc_cntr)
1825 {
1826 PyObject *m;
1827
1828 if( PyType_Ready(&CntrType) < 0 )
1829 return MOD_ERROR_VAL;
1830
1831 MOD_DEF(m, "_nc_cntr", "Contour 2D data", module_methods);
1832
1833 if( m == NULL )
1834 return MOD_ERROR_VAL;
1835
1836 donumpyinit();
1837
1838 Py_INCREF(&CntrType);
1839 PyModule_AddObject(m, "Cntr", (PyObject *)&CntrType);
1840
1841 return MOD_SUCCESS_VAL(m);
1842 }
1843