1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
|
//
// Copyright (C) 2019 Greg Landrum
//
// @@ All Rights Reserved @@
// This file is part of the RDKit.
// The contents are covered by the terms of the BSD license
// which is included in the file license.txt, found at the root
// of the RDKit source tree.
//
#include <vector>
#include <Geometry/point.h>
namespace conrec {
struct ConrecSegment {
RDGeom::Point2D p1;
RDGeom::Point2D p2;
double isoVal;
ConrecSegment(double x1, double y1, double x2, double y2, double isoVal)
: p1(x1, y1), p2(x2, y2), isoVal(isoVal) {}
ConrecSegment(const RDGeom::Point2D &p1, const RDGeom::Point2D &p2,
double isoVal)
: p1(p1), p2(p2), isoVal(isoVal) {}
};
// adapted from conrec.c by Paul Bourke:
// http://paulbourke.net/papers/conrec/conrec.c
/*
Derivation from the fortran version of CONREC by Paul Bourke
d ! matrix of data to contour
ilb,iub,jlb,jub ! index bounds of data matrix
x ! data matrix column coordinates
y ! data matrix row coordinates
nc ! number of contour levels
z ! contour levels in increasing order
*/
inline 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) {
PRECONDITION(d, "no data");
PRECONDITION(x, "no data");
PRECONDITION(y, "no data");
PRECONDITION(z, "no data");
PRECONDITION(nc > 0, "no contours");
PRECONDITION(iub > ilb, "bad bounds");
PRECONDITION(jub > jlb, "bad bounds");
int m1, m2, m3, case_value;
double dmin, dmax, x1 = 0, x2 = 0, y1 = 0, y2 = 0;
int i, j, m;
size_t k;
double h[5];
int sh[5];
double xh[5], yh[5];
int im[4] = {0, 1, 1, 0}, jm[4] = {0, 0, 1, 1};
int castab[3][3][3] = {{{0, 0, 8}, {0, 2, 5}, {7, 6, 9}},
{{0, 3, 4}, {1, 3, 1}, {4, 3, 0}},
{{9, 6, 7}, {5, 2, 0}, {8, 0, 0}}};
double temp1, temp2;
size_t ny = jub - jlb + 1;
auto xsect = [&](int p1, int p2) {
return (h[p2] * xh[p1] - h[p1] * xh[p2]) / (h[p2] - h[p1]);
};
auto ysect = [&](int p1, int p2) {
return (h[p2] * yh[p1] - h[p1] * yh[p2]) / (h[p2] - h[p1]);
};
for (j = (jub - 1); j >= (int)jlb; j--) {
for (i = ilb; i <= (int)iub - 1; i++) {
temp1 = std::min(d[i * ny + j], d[i * ny + j + 1]);
temp2 = std::min(d[(i + 1) * ny + j], d[(i + 1) * ny + j + 1]);
dmin = std::min(temp1, temp2);
temp1 = std::max(d[i * ny + j], d[i * ny + j + 1]);
temp2 = std::max(d[(i + 1) * ny + j], d[(i + 1) * ny + j + 1]);
dmax = std::max(temp1, temp2);
if (dmax < z[0] || dmin > z[nc - 1]) {
continue;
}
for (k = 0; k < nc; k++) {
if (z[k] < dmin || z[k] > dmax) {
continue;
}
for (m = 4; m >= 0; m--) {
if (m > 0) {
h[m] = d[(i + im[m - 1]) * ny + j + jm[m - 1]] - z[k];
xh[m] = x[i + im[m - 1]];
yh[m] = y[j + jm[m - 1]];
} else {
h[0] = 0.25 * (h[1] + h[2] + h[3] + h[4]);
xh[0] = 0.50 * (x[i] + x[i + 1]);
yh[0] = 0.50 * (y[j] + y[j + 1]);
}
if (h[m] > 0.0) {
sh[m] = 1;
} else if (h[m] < 0.0) {
sh[m] = -1;
} else {
sh[m] = 0;
}
}
/*
Note: at this stage the relative heights of the corners and the
centre are in the h array, and the corresponding coordinates are
in the xh and yh arrays. The centre of the box is indexed by 0
and the 4 corners by 1 to 4 as shown below.
Each triangle is then indexed by the parameter m, and the 3
vertices of each triangle are indexed by parameters m1,m2,and m3.
It is assumed that the centre of the box is always vertex 2
though this isimportant only when all 3 vertices lie exactly on
the same contour level, in which case only the side of the box
is drawn.
vertex 4 +-------------------+ vertex 3
| \ / |
| \ m-3 / |
| \ / |
| \ / |
| m=2 X m=2 | the centre is vertex 0
| / \ |
| / \ |
| / m=1 \ |
| / \ |
vertex 1 +-------------------+ vertex 2
*/
/* Scan each triangle in the box */
for (m = 1; m <= 4; m++) {
m1 = m;
m2 = 0;
if (m != 4) {
m3 = m + 1;
} else {
m3 = 1;
}
if ((case_value = castab[sh[m1] + 1][sh[m2] + 1][sh[m3] + 1]) == 0) {
continue;
}
switch (case_value) {
case 1: /* Line between vertices 1 and 2 */
x1 = xh[m1];
y1 = yh[m1];
x2 = xh[m2];
y2 = yh[m2];
break;
case 2: /* Line between vertices 2 and 3 */
x1 = xh[m2];
y1 = yh[m2];
x2 = xh[m3];
y2 = yh[m3];
break;
case 3: /* Line between vertices 3 and 1 */
x1 = xh[m3];
y1 = yh[m3];
x2 = xh[m1];
y2 = yh[m1];
break;
case 4: /* Line between vertex 1 and side 2-3 */
x1 = xh[m1];
y1 = yh[m1];
x2 = xsect(m2, m3);
y2 = ysect(m2, m3);
break;
case 5: /* Line between vertex 2 and side 3-1 */
x1 = xh[m2];
y1 = yh[m2];
x2 = xsect(m3, m1);
y2 = ysect(m3, m1);
break;
case 6: /* Line between vertex 3 and side 1-2 */
x1 = xh[m3];
y1 = yh[m3];
x2 = xsect(m1, m2);
y2 = ysect(m1, m2);
break;
case 7: /* Line between sides 1-2 and 2-3 */
x1 = xsect(m1, m2);
y1 = ysect(m1, m2);
x2 = xsect(m2, m3);
y2 = ysect(m2, m3);
break;
case 8: /* Line between sides 2-3 and 3-1 */
x1 = xsect(m2, m3);
y1 = ysect(m2, m3);
x2 = xsect(m3, m1);
y2 = ysect(m3, m1);
break;
case 9: /* Line between sides 3-1 and 1-2 */
x1 = xsect(m3, m1);
y1 = ysect(m3, m1);
x2 = xsect(m1, m2);
y2 = ysect(m1, m2);
break;
default:
break;
}
/* Finally draw the line */
res.emplace_back(RDGeom::Point2D(x1, y1), RDGeom::Point2D(x2, y2),
z[k]);
} /* m */
} /* k - contour */
} /* i */
} /* j */
}
} // namespace conrec
|
