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 
13 namespace conrec {
14 struct ConrecSegment {
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 */
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) {
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 #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])
48 
49  int m1, m2, m3, case_value;
50  double dmin, dmax, x1 = 0, x2 = 0, y1 = 0, y2 = 0;
51  int i, j, m;
52  size_t k;
53  double h[5];
54  int sh[5];
55  double xh[5], yh[5];
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}}};
60  double temp1, temp2;
61  size_t ny = jub - jlb + 1;
62 
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--) {
75  if (m > 0) {
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]];
79  } else {
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]);
83  }
84  if (h[m] > 0.0)
85  sh[m] = 1;
86  else if (h[m] < 0.0)
87  sh[m] = -1;
88  else
89  sh[m] = 0;
90  }
91 
92  /*
93  Note: at this stage the relative heights of the corners and the
94  centre are in the h array, and the corresponding coordinates are
95  in the xh and yh arrays. The centre of the box is indexed by 0
96  and the 4 corners by 1 to 4 as shown below.
97  Each triangle is then indexed by the parameter m, and the 3
98  vertices of each triangle are indexed by parameters m1,m2,and m3.
99  It is assumed that the centre of the box is always vertex 2
100  though this isimportant only when all 3 vertices lie exactly on
101  the same contour level, in which case only the side of the box
102  is drawn.
103  vertex 4 +-------------------+ vertex 3
104  | \ / |
105  | \ m-3 / |
106  | \ / |
107  | \ / |
108  | m=2 X m=2 | the centre is vertex 0
109  | / \ |
110  | / \ |
111  | / m=1 \ |
112  | / \ |
113  vertex 1 +-------------------+ vertex 2
114  */
115  /* Scan each triangle in the box */
116  for (m = 1; m <= 4; m++) {
117  m1 = m;
118  m2 = 0;
119  if (m != 4)
120  m3 = m + 1;
121  else
122  m3 = 1;
123  if ((case_value = castab[sh[m1] + 1][sh[m2] + 1][sh[m3] + 1]) == 0)
124  continue;
125  switch (case_value) {
126  case 1: /* Line between vertices 1 and 2 */
127  x1 = xh[m1];
128  y1 = yh[m1];
129  x2 = xh[m2];
130  y2 = yh[m2];
131  break;
132  case 2: /* Line between vertices 2 and 3 */
133  x1 = xh[m2];
134  y1 = yh[m2];
135  x2 = xh[m3];
136  y2 = yh[m3];
137  break;
138  case 3: /* Line between vertices 3 and 1 */
139  x1 = xh[m3];
140  y1 = yh[m3];
141  x2 = xh[m1];
142  y2 = yh[m1];
143  break;
144  case 4: /* Line between vertex 1 and side 2-3 */
145  x1 = xh[m1];
146  y1 = yh[m1];
147  x2 = xsect(m2, m3);
148  y2 = ysect(m2, m3);
149  break;
150  case 5: /* Line between vertex 2 and side 3-1 */
151  x1 = xh[m2];
152  y1 = yh[m2];
153  x2 = xsect(m3, m1);
154  y2 = ysect(m3, m1);
155  break;
156  case 6: /* Line between vertex 3 and side 1-2 */
157  x1 = xh[m3];
158  y1 = yh[m3];
159  x2 = xsect(m1, m2);
160  y2 = ysect(m1, m2);
161  break;
162  case 7: /* Line between sides 1-2 and 2-3 */
163  x1 = xsect(m1, m2);
164  y1 = ysect(m1, m2);
165  x2 = xsect(m2, m3);
166  y2 = ysect(m2, m3);
167  break;
168  case 8: /* Line between sides 2-3 and 3-1 */
169  x1 = xsect(m2, m3);
170  y1 = ysect(m2, m3);
171  x2 = xsect(m3, m1);
172  y2 = ysect(m3, m1);
173  break;
174  case 9: /* Line between sides 3-1 and 1-2 */
175  x1 = xsect(m3, m1);
176  y1 = ysect(m3, m1);
177  x2 = xsect(m1, m2);
178  y2 = ysect(m1, m2);
179  break;
180  default:
181  break;
182  }
183 
184  /* Finally draw the line */
185  res.push_back(ConrecSegment(RDGeom::Point2D(x1, y1),
186  RDGeom::Point2D(x2, y2), z[k]));
187  } /* m */
188  } /* k - contour */
189  } /* i */
190  } /* j */
191 }
192 } // namespace conrec
conrec
Definition: Conrec.h:13
conrec::ConrecSegment::p1
RDGeom::Point2D p1
Definition: Conrec.h:15
point.h
conrec::ConrecSegment::p2
RDGeom::Point2D p2
Definition: Conrec.h:16
ysect
#define ysect(p1, p2)
conrec::Contour
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
conrec::ConrecSegment
Definition: Conrec.h:14
conrec::ConrecSegment::isoVal
double isoVal
Definition: Conrec.h:17
RDGeom::Point2D
Definition: point.h:258
PRECONDITION
#define PRECONDITION(expr, mess)
Definition: Invariant.h:109
xsect
#define xsect(p1, p2)
conrec::ConrecSegment::ConrecSegment
ConrecSegment(const RDGeom::Point2D &p1, const RDGeom::Point2D &p2, double isoVal)
Definition: Conrec.h:20
conrec::ConrecSegment::ConrecSegment
ConrecSegment(double x1, double y1, double x2, double y2, double isoVal)
Definition: Conrec.h:18