1 // --------------------------------------------------------------------------------------------------------------------
2 // <copyright file="Conrec.cs" company="OxyPlot">
3 // The MIT License (MIT)
5 // Copyright (c) 2012 Oystein Bjorke
7 // Permission is hereby granted, free of charge, to any person obtaining a
8 // copy of this software and associated documentation files (the
9 // "Software"), to deal in the Software without restriction, including
10 // without limitation the rights to use, copy, modify, merge, publish,
11 // distribute, sublicense, and/or sell copies of the Software, and to
12 // permit persons to whom the Software is furnished to do so, subject to
13 // the following conditions:
15 // The above copyright notice and this permission notice shall be included
16 // in all copies or substantial portions of the Software.
18 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
19 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
20 // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
21 // IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
22 // CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
23 // TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
24 // SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
27 // Creates contours from a triangular mesh.
29 // --------------------------------------------------------------------------------------------------------------------
35 /// Provides functionality to create contours from a triangular mesh.
39 /// Ported from C / Fortran code by Paul Bourke.
40 /// See <a href="http://paulbourke.net/papers/conrec/">Conrec</a> for
41 /// full description of code and the original source.
44 /// Contouring aids in visualizing three dimensional surfaces on a two dimensional
45 /// medium (on paper or in this case a computer graphics screen). Two most common
46 /// applications are displaying topological features of an area on a map or the air
47 /// pressure on a weather map. In all cases some parameter is plotted as a function
48 /// of two variables, the longitude and latitude or x and y axis. One problem with
49 /// computer contouring is the process is usually CPU intensive and the algorithms
50 /// often use advanced mathematical techniques making them susceptible to error.
53 public static class Conrec
59 /// Start point x-coordinate
62 /// Start point y-coordinate
65 /// End point x-coordinate
68 /// End point y-coordinate
73 public delegate void RendererDelegate(double x1, double y1, double x2, double y2, double z);
76 /// Contour is a contouring subroutine for rectangularily spaced data
77 /// It emits calls to a line drawing subroutine supplied by the user
78 /// which draws a contour map corresponding to data on a randomly
79 /// spaced rectangular grid. The coordinates emitted are in the same
80 /// units given in the x() and y() arrays.
81 /// Any number of contour levels may be specified but they must be
82 /// in order of increasing value.
85 /// Matrix of data to contour.
88 /// Data matrix column coordinates.
91 /// Data matrix row coordinates.
94 /// Contour levels in increasing order.
96 /// <param name="renderer">
99 public static void Contour(double[,] d, double[] x, double[] y, double[] z, RendererDelegate renderer)
106 var h = new double[5];
108 var xh = new double[5];
109 var yh = new double[5];
111 int ilb = d.GetLowerBound(0);
112 int iub = d.GetUpperBound(0);
113 int jlb = d.GetLowerBound(1);
114 int jub = d.GetUpperBound(1);
117 // The indexing of im and jm should be noted as it has to start from zero
118 // unlike the fortran counter part
119 int[] im = { 0, 1, 1, 0 };
120 int[] jm = { 0, 0, 1, 1 };
122 // Note that castab is arranged differently from the FORTRAN code because
123 // Fortran and C/C++ arrays are transposed of each other, in this case
124 // it is more tricky as castab is in 3 dimension
127 { { 0, 0, 8 }, { 0, 2, 5 }, { 7, 6, 9 } }, { { 0, 3, 4 }, { 1, 3, 1 }, { 4, 3, 0 } },
128 { { 9, 6, 7 }, { 5, 2, 0 }, { 8, 0, 0 } }
131 Func<int, int, double> xsect = (p1, p2) => ((h[p2] * xh[p1]) - (h[p1] * xh[p2])) / (h[p2] - h[p1]);
132 Func<int, int, double> ysect = (p1, p2) => ((h[p2] * yh[p1]) - (h[p1] * yh[p2])) / (h[p2] - h[p1]);
134 for (int j = jub - 1; j >= jlb; j--)
137 for (i = ilb; i <= iub - 1; i++)
139 double temp1 = Math.Min(d[i, j], d[i, j + 1]);
140 double temp2 = Math.Min(d[i + 1, j], d[i + 1, j + 1]);
141 double dmin = Math.Min(temp1, temp2);
142 temp1 = Math.Max(d[i, j], d[i, j + 1]);
143 temp2 = Math.Max(d[i + 1, j], d[i + 1, j + 1]);
144 double dmax = Math.Max(temp1, temp2);
146 if (dmax >= z[0] && dmin <= z[nc - 1])
149 for (k = 0; k < nc; k++)
151 if (z[k] >= dmin && z[k] <= dmax)
154 for (m = 4; m >= 0; m--)
158 // The indexing of im and jm should be noted as it has to
160 h[m] = d[i + im[m - 1], j + jm[m - 1]] - z[k];
161 xh[m] = x[i + im[m - 1]];
162 yh[m] = y[j + jm[m - 1]];
166 h[0] = 0.25 * (h[1] + h[2] + h[3] + h[4]);
167 xh[0] = 0.5 * (x[i] + x[i + 1]);
168 yh[0] = 0.5 * (y[j] + y[j + 1]);
185 //// Note: at this stage the relative heights of the corners and the
186 //// centre are in the h array, and the corresponding coordinates are
187 //// in the xh and yh arrays. The centre of the box is indexed by 0
188 //// and the 4 corners by 1 to 4 as shown below.
189 //// Each triangle is then indexed by the parameter m, and the 3
190 //// vertices of each triangle are indexed by parameters m1,m2,and
192 //// It is assumed that the centre of the box is always vertex 2
193 //// though this isimportant only when all 3 vertices lie exactly on
194 //// the same contour level, in which case only the side of the box
196 //// vertex 4 +-------------------+ vertex 3
201 //// | m=2 X m=2 | the centre is vertex 0
206 //// vertex 1 +-------------------+ vertex 2
208 // Scan each triangle in the box
209 for (m = 1; m <= 4; m++)
223 int caseValue = castab[sh[m1] + 1, sh[m2] + 1, sh[m3] + 1];
228 case 1: // Line between vertices 1 and 2
234 case 2: // Line between vertices 2 and 3
240 case 3: // Line between vertices 3 and 1
246 case 4: // Line between vertex 1 and side 2-3
252 case 5: // Line between vertex 2 and side 3-1
258 case 6: // Line between vertex 3 and side 1-2
264 case 7: // Line between sides 1-2 and 2-3
270 case 8: // Line between sides 2-3 and 3-1
276 case 9: // Line between sides 3-1 and 1-2
284 renderer(x1, y1, x2, y2, z[k]);