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