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 int m1, m2, m3, case_value;
47 double dmin, dmax, x1 = 0, x2 = 0, y1 = 0, y2 = 0;
53 int im[4] = {0, 1, 1, 0}, jm[4] = {0, 0, 1, 1};
54 int castab[3][3][3] = {{{0, 0, 8}, {0, 2, 5}, {7, 6, 9}},
55 {{0, 3, 4}, {1, 3, 1}, {4, 3, 0}},
56 {{9, 6, 7}, {5, 2, 0}, {8, 0, 0}}};
58 size_t ny = jub - jlb + 1;
60 auto xsect = [&](
int p1,
int p2) {
61 return (h[p2] * xh[p1] - h[p1] * xh[p2]) / (h[p2] - h[p1]);
63 auto ysect = [&](
int p1,
int p2) {
64 return (h[p2] * yh[p1] - h[p1] * yh[p2]) / (h[p2] - h[p1]);
67 for (j = (jub - 1); j >= (int)jlb; j--) {
68 for (i = ilb; i <= (int)iub - 1; i++) {
69 temp1 = std::min(d[i * ny + j], d[i * ny + j + 1]);
70 temp2 = std::min(d[(i + 1) * ny + j], d[(i + 1) * ny + j + 1]);
71 dmin = std::min(temp1, temp2);
72 temp1 = std::max(d[i * ny + j], d[i * ny + j + 1]);
73 temp2 = std::max(d[(i + 1) * ny + j], d[(i + 1) * ny + j + 1]);
74 dmax = std::max(temp1, temp2);
75 if (dmax < z[0] || dmin > z[nc - 1])
continue;
76 for (k = 0; k < nc; k++) {
77 if (z[k] < dmin || z[k] > dmax)
continue;
78 for (m = 4; m >= 0; m--) {
80 h[m] = d[(i + im[m - 1]) * ny + j + jm[m - 1]] - z[k];
81 xh[m] = x[i + im[m - 1]];
82 yh[m] = y[j + jm[m - 1]];
84 h[0] = 0.25 * (h[1] + h[2] + h[3] + h[4]);
85 xh[0] = 0.50 * (x[i] + x[i + 1]);
86 yh[0] = 0.50 * (y[j] + y[j + 1]);
120 for (m = 1; m <= 4; m++) {
127 if ((case_value = castab[sh[m1] + 1][sh[m2] + 1][sh[m3] + 1]) == 0)
129 switch (case_value) {
#define PRECONDITION(expr, mess)
void Contour(const double *d, size_t ilb, size_t iub, size_t jlb, size_t jub, const double *x, const double *y, size_t nc, double *z, std::vector< ConrecSegment > &res)
ConrecSegment(const RDGeom::Point2D &p1, const RDGeom::Point2D &p2, double isoVal)
ConrecSegment(double x1, double y1, double x2, double y2, double isoVal)