1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Contour lines for function evaluated on a grid.
4 //
5 // Copyright (C) 2001-2021 The Octave Project Developers
6 //
7 // See the file COPYRIGHT.md in the top-level directory of this
8 // distribution or <https://octave.org/copyright/>.
9 //
10 // Adapted to an oct file from the stand alone contourl by Victro Munoz
11 // Copyright (C) 2004 Victor Munoz
12 //
13 // Based on contour plot routine (plcont.c) in PLPlot package
14 // http://plplot.org/
15 //
16 // Copyright (C) 1995, 2000, 2001 Maurice LeBrun
17 // Copyright (C) 2000, 2002 Joao Cardoso
18 // Copyright (C) 2000, 2001, 2002, 2004 Alan W. Irwin
19 // Copyright (C) 2004 Andrew Ross
20 //
21 //
22 // This file is part of Octave.
23 //
24 // Octave is free software: you can redistribute it and/or modify it
25 // under the terms of the GNU General Public License as published by
26 // the Free Software Foundation, either version 3 of the License, or
27 // (at your option) any later version.
28 //
29 // Octave is distributed in the hope that it will be useful, but
30 // WITHOUT ANY WARRANTY; without even the implied warranty of
31 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 // GNU General Public License for more details.
33 //
34 // You should have received a copy of the GNU General Public License
35 // along with Octave; see the file COPYING. If not, see
36 // <https://www.gnu.org/licenses/>.
37 //
38 ////////////////////////////////////////////////////////////////////////
39
40 #if defined (HAVE_CONFIG_H)
41 # include "config.h"
42 #endif
43
44 #include <limits>
45
46 #include "defun.h"
47 #include "ov.h"
48
49 // FIXME: this looks like trouble...
50 static Matrix this_contour;
51 static Matrix contourc;
52 static int elem;
53
54 // This is the quanta in which we increase this_contour.
55 #define CONTOUR_QUANT 50
56
57 // Add a coordinate point (x,y) to this_contour.
58
59 static void
add_point(double x,double y)60 add_point (double x, double y)
61 {
62 if (elem % CONTOUR_QUANT == 0)
63 this_contour = this_contour.append (Matrix (2, CONTOUR_QUANT, 0));
64
65 this_contour (0, elem) = x;
66 this_contour (1, elem) = y;
67 elem++;
68 }
69
70 // Add contents of current contour to contourc.
71 // this_contour.cols () - 1;
72
73 static void
end_contour(void)74 end_contour (void)
75 {
76 if (elem > 2)
77 {
78 this_contour (1, 0) = elem - 1;
79 contourc = contourc.append (this_contour.extract_n (0, 0, 2, elem));
80 }
81
82 this_contour = Matrix ();
83 elem = 0;
84 }
85
86 // Start a new contour, and add contents of current one to contourc.
87
88 static void
start_contour(double lvl,double x,double y)89 start_contour (double lvl, double x, double y)
90 {
91 end_contour ();
92 this_contour.resize (2, 0);
93 add_point (lvl, 0);
94 add_point (x, y);
95 }
96
97 static void
drawcn(const RowVector & X,const RowVector & Y,const Matrix & Z,double lvl,int r,int c,double ct_x,double ct_y,unsigned int start_edge,bool first,charMatrix & mark)98 drawcn (const RowVector& X, const RowVector& Y, const Matrix& Z,
99 double lvl, int r, int c, double ct_x, double ct_y,
100 unsigned int start_edge, bool first, charMatrix& mark)
101 {
102 double px[4], py[4], pz[4], tmp;
103 unsigned int stop_edge, pt[2];
104
105 // Continue while next facet is not done yet.
106 while (r >= 0 && c >= 0 && r < mark.rows () && c < mark.cols ()
107 && mark(r, c) > 0)
108 {
109
110 //get x, y, and z - lvl for current facet
111 px[0] = px[3] = X(c);
112 px[1] = px[2] = X(c+1);
113
114 py[0] = py[1] = Y(r);
115 py[2] = py[3] = Y(r+1);
116
117 pz[3] = Z(r+1, c) - lvl;
118 pz[2] = Z(r+1, c + 1) - lvl;
119 pz[1] = Z(r, c+1) - lvl;
120 pz[0] = Z(r, c) - lvl;
121
122 // Facet edge and point naming assignment.
123 //
124 // 0-----1 .-0-.
125 // | | | |
126 // | | 3 1
127 // | | | |
128 // 3-----2 .-2-.
129
130 // Get mark value of current facet.
131 char id = static_cast<char> (mark(r, c));
132
133 // Check startedge s.
134 if (start_edge == 255)
135 {
136 // Find start edge.
137 for (unsigned int k = 0; k < 4; k++)
138 if (static_cast<char> (1 << k) & id)
139 start_edge = k;
140 }
141
142 if (start_edge == 255)
143 break;
144
145 // Decrease mark value of current facet for start edge.
146 mark(r, c) -= static_cast<char> (1 << start_edge);
147
148 // Next point (clockwise).
149 pt[0] = start_edge;
150 pt[1] = (pt[0] + 1) % 4;
151
152 // Calculate contour segment start if first of contour.
153 if (first)
154 {
155 tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]);
156
157 if (octave::math::isnan (tmp))
158 ct_x = ct_y = 0.5;
159 else
160 {
161 ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp);
162 ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp);
163 }
164
165 start_contour (lvl, ct_x, ct_y);
166 first = false;
167 }
168
169 // Find stop edge.
170 // FIXME: perhaps this should use a while loop?
171 for (unsigned int k = 1; k <= 4; k++)
172 {
173 if (start_edge == 0 || start_edge == 2)
174 stop_edge = (start_edge + k) % 4;
175 else
176 stop_edge = (start_edge - k) % 4;
177
178 if (static_cast<char> (1 << stop_edge) & id)
179 break;
180 }
181
182 pt[0] = stop_edge;
183 pt[1] = (pt[0] + 1) % 4;
184 tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]);
185
186 if (octave::math::isnan (tmp))
187 ct_x = ct_y = 0.5;
188 else
189 {
190 ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp);
191 ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp);
192 }
193
194 // Add point to contour.
195 add_point (ct_x, ct_y);
196
197 // Decrease id value of current facet for start edge.
198 mark(r, c) -= static_cast<char> (1 << stop_edge);
199
200 // Find next facet.
201 if (stop_edge == 0)
202 r--;
203 else if (stop_edge == 1)
204 c++;
205 else if (stop_edge == 2)
206 r++;
207 else if (stop_edge == 3)
208 c--;
209
210 // Go to next facet.
211 start_edge = (stop_edge + 2) % 4;
212
213 }
214 }
215
216 static void
mark_facets(const Matrix & Z,charMatrix & mark,double lvl)217 mark_facets (const Matrix& Z, charMatrix& mark, double lvl)
218 {
219 unsigned int nr = mark.rows ();
220 unsigned int nc = mark.cols ();
221
222 double f[4];
223
224 for (unsigned int c = 0; c < nc; c++)
225 for (unsigned int r = 0; r < nr; r++)
226 {
227 f[0] = Z(r, c) - lvl;
228 f[1] = Z(r, c+1) - lvl;
229 f[3] = Z(r+1, c) - lvl;
230 f[2] = Z(r+1, c+1) - lvl;
231
232 for (unsigned int i = 0; i < 4; i++)
233 if (fabs(f[i]) < std::numeric_limits<double>::epsilon ())
234 f[i] = std::numeric_limits<double>::epsilon ();
235
236 if (f[1] * f[2] < 0)
237 mark(r, c) += 2;
238
239 if (f[0] * f[3] < 0)
240 mark(r, c) += 8;
241 }
242
243 for (unsigned int r = 0; r < nr; r++)
244 for (unsigned int c = 0; c < nc; c++)
245 {
246 f[0] = Z(r, c) - lvl;
247 f[1] = Z(r, c+1) - lvl;
248 f[3] = Z(r+1, c) - lvl;
249 f[2] = Z(r+1, c+1) - lvl;
250
251 for (unsigned int i = 0; i < 4; i++)
252 if (fabs(f[i]) < std::numeric_limits<double>::epsilon ())
253 f[i] = std::numeric_limits<double>::epsilon ();
254
255 if (f[0] * f[1] < 0)
256 mark(r, c) += 1;
257
258 if (f[2] * f[3] < 0)
259 mark(r, c) += 4;
260 }
261 }
262
263 static void
cntr(const RowVector & X,const RowVector & Y,const Matrix & Z,double lvl)264 cntr (const RowVector& X, const RowVector& Y, const Matrix& Z, double lvl)
265 {
266 unsigned int nr = Z.rows ();
267 unsigned int nc = Z.cols ();
268
269 charMatrix mark (nr - 1, nc - 1, 0);
270
271 mark_facets (Z, mark, lvl);
272
273 // Find contours that start at a domain edge.
274
275 for (unsigned int c = 0; c < nc - 1; c++)
276 {
277 // Top.
278 if (mark(0, c) & 1)
279 drawcn (X, Y, Z, lvl, 0, c, 0.0, 0.0, 0, true, mark);
280
281 // Bottom.
282 if (mark(nr - 2, c) & 4)
283 drawcn (X, Y, Z, lvl, nr - 2, c, 0.0, 0.0, 2, true, mark);
284 }
285
286 for (unsigned int r = 0; r < nr - 1; r++)
287 {
288 // Left.
289 if (mark(r, 0) & 8)
290 drawcn (X, Y, Z, lvl, r, 0, 0.0, 0.0, 3, true, mark);
291
292 // Right.
293 if (mark(r, nc - 2) & 2)
294 drawcn (X, Y, Z, lvl, r, nc - 2, 0.0, 0.0, 1, true, mark);
295 }
296
297 for (unsigned int r = 0; r < nr - 1; r++)
298 for (unsigned int c = 0; c < nc - 1; c++)
299 if (mark (r, c) > 0)
300 drawcn (X, Y, Z, lvl, r, c, 0.0, 0.0, 255, true, mark);
301 }
302
303 DEFUN (__contourc__, args, ,
304 doc: /* -*- texinfo -*-
305 @deftypefn {} {} __contourc__ (@var{x}, @var{y}, @var{z}, @var{levels})
306 Undocumented internal function.
307 @end deftypefn */)
308 {
309 if (args.length () != 4)
310 print_usage ();
311
312 RowVector X = args(0).row_vector_value ();
313 RowVector Y = args(1).row_vector_value ();
314 Matrix Z = args(2).matrix_value ();
315 RowVector L = args(3).row_vector_value ();
316
317 contourc.resize (2, 0);
318
319 for (int i = 0; i < L.numel (); i++)
320 cntr (X, Y, Z, L (i));
321
322 end_contour ();
323
324 return octave_value (contourc);
325 }
326
327 /*
328 ## No test needed for internal helper function.
329 %!assert (1)
330 */
331