35 inline void Contour(
const double *d,
size_t ilb,
size_t iub,
size_t jlb,
36 size_t jub,
const double *x,
const double *y,
size_t nc,
37 double *z, std::vector<ConrecSegment> &res) {
46 #define xsect(p1, p2) (h[p2] * xh[p1] - h[p1] * xh[p2]) / (h[p2] - h[p1])
47 #define ysect(p1, p2) (h[p2] * yh[p1] - h[p1] * yh[p2]) / (h[p2] - h[p1])
49 int m1, m2, m3, case_value;
50 double dmin, dmax, x1 = 0, x2 = 0, y1 = 0, y2 = 0;
56 int im[4] = {0, 1, 1, 0}, jm[4] = {0, 0, 1, 1};
57 int castab[3][3][3] = {{{0, 0, 8}, {0, 2, 5}, {7, 6, 9}},
58 {{0, 3, 4}, {1, 3, 1}, {4, 3, 0}},
59 {{9, 6, 7}, {5, 2, 0}, {8, 0, 0}}};
61 size_t ny = jub - jlb + 1;
63 for (j = (jub - 1); j >= (int)jlb; j--) {
64 for (i = ilb; i <= (int)iub - 1; i++) {
65 temp1 = std::min(d[i * ny + j], d[i * ny + j + 1]);
66 temp2 = std::min(d[(i + 1) * ny + j], d[(i + 1) * ny + j + 1]);
67 dmin = std::min(temp1, temp2);
68 temp1 = std::max(d[i * ny + j], d[i * ny + j + 1]);
69 temp2 = std::max(d[(i + 1) * ny + j], d[(i + 1) * ny + j + 1]);
70 dmax = std::max(temp1, temp2);
71 if (dmax < z[0] || dmin > z[nc - 1])
continue;
72 for (k = 0; k < nc; k++) {
73 if (z[k] < dmin || z[k] > dmax)
continue;
74 for (m = 4; m >= 0; m--) {
76 h[m] = d[(i + im[m - 1]) * ny + j + jm[m - 1]] - z[k];
77 xh[m] = x[i + im[m - 1]];
78 yh[m] = y[j + jm[m - 1]];
80 h[0] = 0.25 * (h[1] + h[2] + h[3] + h[4]);
81 xh[0] = 0.50 * (x[i] + x[i + 1]);
82 yh[0] = 0.50 * (y[j] + y[j + 1]);
116 for (m = 1; m <= 4; m++) {
123 if ((case_value = castab[sh[m1] + 1][sh[m2] + 1][sh[m3] + 1]) == 0)
125 switch (case_value) {