GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
__contourc__.cc
Go to the documentation of this file.
1 /* Contour lines for function evaluated on a grid.
2 
3 Copyright (C) 2007-2013 Kai Habel
4 Copyright (C) 2004, 2007 Shai Ayal
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
21 Free Software Foundation; either version 3 of the License, or (at your
22 option) any later version.
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
30 along with Octave; see the file COPYING. If not, see
31 <http://www.gnu.org/licenses/>.
32 
33 */
34 
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38 
39 #include <cfloat>
40 
41 #include "quit.h"
42 
43 #include "defun.h"
44 #include "error.h"
45 #include "oct-obj.h"
46 
49 static int elem;
50 
51 // This is the quanta in which we increase this_contour.
52 #define CONTOUR_QUANT 50
53 
54 // Add a coordinate point (x,y) to this_contour.
55 
56 static void
57 add_point (double x, double y)
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 
67 // Add contents of current contour to contourc.
68 // this_contour.cols () - 1;
69 
70 static void
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  }
78 
79  this_contour = Matrix ();
80  elem = 0;
81 }
82 
83 // Start a new contour, and add contents of current one to contourc.
84 
85 static void
86 start_contour (double lvl, double x, double y)
87 {
88  end_contour ();
89  this_contour.resize (2, 0);
90  add_point (lvl, 0);
91  add_point (x, y);
92 }
93 
94 static void
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  unsigned int start_edge, bool first, charMatrix& mark)
98 {
99  double px[4], py[4], pz[4], tmp;
100  unsigned int stop_edge, pt[2];
101 
102  // Continue while next facet is not done yet.
103  while (r >= 0 && c >= 0 && r < mark.rows () && c < mark.cols ()
104  && mark(r, c) > 0)
105  {
106 
107  //get x, y, and z - lvl for current facet
108  px[0] = px[3] = X(c);
109  px[1] = px[2] = X(c+1);
110 
111  py[0] = py[1] = Y(r);
112  py[2] = py[3] = Y(r+1);
113 
114  pz[3] = Z(r+1, c) - lvl;
115  pz[2] = Z(r+1, c + 1) - lvl;
116  pz[1] = Z(r, c+1) - lvl;
117  pz[0] = Z(r, c) - lvl;
118 
119  // Facet edge and point naming assignment.
120  //
121  // 0-----1 .-0-.
122  // | | | |
123  // | | 3 1
124  // | | | |
125  // 3-----2 .-2-.
126 
127  // Get mark value of current facet.
128  char id = static_cast<char> (mark(r, c));
129 
130  // Check startedge s.
131  if (start_edge == 255)
132  {
133  // Find start edge.
134  for (unsigned int k = 0; k < 4; k++)
135  if (static_cast<char> (1 << k) & id)
136  start_edge = k;
137  }
138 
139  if (start_edge == 255)
140  break;
141 
142  // Decrease mark value of current facet for start edge.
143  mark(r, c) -= static_cast<char> (1 << start_edge);
144 
145  // Next point (clockwise).
146  pt[0] = start_edge;
147  pt[1] = (pt[0] + 1) % 4;
148 
149  // Calculate contour segment start if first of contour.
150  if (first)
151  {
152  tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]);
153 
154  if (xisnan (tmp))
155  ct_x = ct_y = 0.5;
156  else
157  {
158  ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp);
159  ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp);
160  }
161 
162  start_contour (lvl, ct_x, ct_y);
163  first = false;
164  }
165 
166  // Find stop edge.
167  // FIXME: perhaps this should use a while loop?
168  for (unsigned int k = 1; k <= 4; k++)
169  {
170  if (start_edge == 0 || start_edge == 2)
171  stop_edge = (start_edge + k) % 4;
172  else
173  stop_edge = (start_edge - k) % 4;
174 
175  if (static_cast<char> (1 << stop_edge) & id)
176  break;
177  }
178 
179  pt[0] = stop_edge;
180  pt[1] = (pt[0] + 1) % 4;
181  tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]);
182 
183  if (xisnan (tmp))
184  ct_x = ct_y = 0.5;
185  else
186  {
187  ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp);
188  ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp);
189  }
190 
191  // Add point to contour.
192  add_point (ct_x, ct_y);
193 
194  // Decrease id value of current facet for start edge.
195  mark(r, c) -= static_cast<char> (1 << stop_edge);
196 
197  // Find next facet.
198  if (stop_edge == 0)
199  r--;
200  else if (stop_edge == 1)
201  c++;
202  else if (stop_edge == 2)
203  r++;
204  else if (stop_edge == 3)
205  c--;
206 
207  // Go to next facet.
208  start_edge = (stop_edge + 2) % 4;
209 
210  }
211 }
212 
213 static void
214 mark_facets (const Matrix& Z, charMatrix& mark, double lvl)
215 {
216  unsigned int nr = mark.rows ();
217  unsigned int nc = mark.cols ();
218 
219  double f[4];
220 
221  for (unsigned int c = 0; c < nc; c++)
222  for (unsigned int r = 0; r < nr; r++)
223  {
224  f[0] = Z(r, c) - lvl;
225  f[1] = Z(r, c+1) - lvl;
226  f[3] = Z(r+1, c) - lvl;
227  f[2] = Z(r+1, c+1) - lvl;
228 
229  for (unsigned int i = 0; i < 4; i++)
230  if (fabs(f[i]) < std::numeric_limits<double>::epsilon ())
231  f[i] = std::numeric_limits<double>::epsilon ();
232 
233  if (f[1] * f[2] < 0)
234  mark(r, c) += 2;
235 
236  if (f[0] * f[3] < 0)
237  mark(r, c) += 8;
238  }
239 
240  for (unsigned int r = 0; r < nr; r++)
241  for (unsigned int c = 0; c < nc; c++)
242  {
243  f[0] = Z(r, c) - lvl;
244  f[1] = Z(r, c+1) - lvl;
245  f[3] = Z(r+1, c) - lvl;
246  f[2] = Z(r+1, c+1) - lvl;
247 
248  for (unsigned int i = 0; i < 4; i++)
249  if (fabs(f[i]) < std::numeric_limits<double>::epsilon ())
250  f[i] = std::numeric_limits<double>::epsilon ();
251 
252  if (f[0] * f[1] < 0)
253  mark(r, c) += 1;
254 
255  if (f[2] * f[3] < 0)
256  mark(r, c) += 4;
257  }
258 }
259 
260 static void
261 cntr (const RowVector& X, const RowVector& Y, const Matrix& Z, double lvl)
262 {
263  unsigned int nr = Z.rows ();
264  unsigned int nc = Z.cols ();
265 
266  charMatrix mark (nr - 1, nc - 1, 0);
267 
268  mark_facets (Z, mark, lvl);
269 
270  // Find contours that start at a domain edge.
271 
272  for (unsigned int c = 0; c < nc - 1; c++)
273  {
274  // Top.
275  if (mark(0, c) & 1)
276  drawcn (X, Y, Z, lvl, 0, c, 0.0, 0.0, 0, true, mark);
277 
278  // Bottom.
279  if (mark(nr - 2, c) & 4)
280  drawcn (X, Y, Z, lvl, nr - 2, c, 0.0, 0.0, 2, true, mark);
281  }
282 
283  for (unsigned int r = 0; r < nr - 1; r++)
284  {
285  // Left.
286  if (mark(r, 0) & 8)
287  drawcn (X, Y, Z, lvl, r, 0, 0.0, 0.0, 3, true, mark);
288 
289  // Right.
290  if (mark(r, nc - 2) & 2)
291  drawcn (X, Y, Z, lvl, r, nc - 2, 0.0, 0.0, 1, true, mark);
292  }
293 
294  for (unsigned int r = 0; r < nr - 1; r++)
295  for (unsigned int c = 0; c < nc - 1; c++)
296  if (mark (r, c) > 0)
297  drawcn (X, Y, Z, lvl, r, c, 0.0, 0.0, 255, true, mark);
298 }
299 
300 DEFUN (__contourc__, args, ,
301  "-*- texinfo -*-\n\
302 @deftypefn {Built-in Function} {} __contourc__ (@var{x}, @var{y}, @var{z}, @var{levels})\n\
303 Undocumented internal function.\n\
304 @end deftypefn")
305 {
306  octave_value retval;
307 
308  if (args.length () == 4)
309  {
310  RowVector X = args (0).row_vector_value ();
311  RowVector Y = args (1).row_vector_value ();
312  Matrix Z = args (2).matrix_value ();
313  RowVector L = args (3).row_vector_value ();
314 
315  if (! error_state)
316  {
317  contourc.resize (2, 0);
318 
319  for (int i = 0; i < L.length (); i++)
320  cntr (X, Y, Z, L (i));
321 
322  end_contour ();
323 
324  retval = contourc;
325  }
326  else
327  error ("__contourc__: invalid argument values");
328  }
329  else
330  print_usage ();
331 
332  return retval;
333 }
334 
335 /*
336 ## No test needed for internal helper function.
337 %!assert (1)
338 */