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