GNU Octave  4.2.1
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-2017 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 #if defined (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 "ovl.h"
46 
47 // FIXME: this looks like trouble...
50 static int elem;
51 
52 // This is the quanta in which we increase this_contour.
53 #define CONTOUR_QUANT 50
54 
55 // Add a coordinate point (x,y) to this_contour.
56 
57 static void
58 add_point (double x, double y)
59 {
60  if (elem % CONTOUR_QUANT == 0)
61  this_contour = this_contour.append (Matrix (2, CONTOUR_QUANT, 0));
62 
63  this_contour (0, elem) = x;
64  this_contour (1, elem) = y;
65  elem++;
66 }
67 
68 // Add contents of current contour to contourc.
69 // this_contour.cols () - 1;
70 
71 static void
73 {
74  if (elem > 2)
75  {
76  this_contour (1, 0) = elem - 1;
77  contourc = contourc.append (this_contour.extract_n (0, 0, 2, elem));
78  }
79 
80  this_contour = Matrix ();
81  elem = 0;
82 }
83 
84 // Start a new contour, and add contents of current one to contourc.
85 
86 static void
87 start_contour (double lvl, double x, double y)
88 {
89  end_contour ();
90  this_contour.resize (2, 0);
91  add_point (lvl, 0);
92  add_point (x, y);
93 }
94 
95 static void
96 drawcn (const RowVector& X, const RowVector& Y, const Matrix& Z,
97  double lvl, int r, int c, double ct_x, double ct_y,
98  unsigned int start_edge, bool first, charMatrix& mark)
99 {
100  double px[4], py[4], pz[4], tmp;
101  unsigned int stop_edge, pt[2];
102 
103  // Continue while next facet is not done yet.
104  while (r >= 0 && c >= 0 && r < mark.rows () && c < mark.cols ()
105  && mark(r, c) > 0)
106  {
107 
108  //get x, y, and z - lvl for current facet
109  px[0] = px[3] = X(c);
110  px[1] = px[2] = X(c+1);
111 
112  py[0] = py[1] = Y(r);
113  py[2] = py[3] = Y(r+1);
114 
115  pz[3] = Z(r+1, c) - lvl;
116  pz[2] = Z(r+1, c + 1) - lvl;
117  pz[1] = Z(r, c+1) - lvl;
118  pz[0] = Z(r, c) - lvl;
119 
120  // Facet edge and point naming assignment.
121  //
122  // 0-----1 .-0-.
123  // | | | |
124  // | | 3 1
125  // | | | |
126  // 3-----2 .-2-.
127 
128  // Get mark value of current facet.
129  char id = static_cast<char> (mark(r, c));
130 
131  // Check startedge s.
132  if (start_edge == 255)
133  {
134  // Find start edge.
135  for (unsigned int k = 0; k < 4; k++)
136  if (static_cast<char> (1 << k) & id)
137  start_edge = k;
138  }
139 
140  if (start_edge == 255)
141  break;
142 
143  // Decrease mark value of current facet for start edge.
144  mark(r, c) -= static_cast<char> (1 << start_edge);
145 
146  // Next point (clockwise).
147  pt[0] = start_edge;
148  pt[1] = (pt[0] + 1) % 4;
149 
150  // Calculate contour segment start if first of contour.
151  if (first)
152  {
153  tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]);
154 
155  if (octave::math::isnan (tmp))
156  ct_x = ct_y = 0.5;
157  else
158  {
159  ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp);
160  ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp);
161  }
162 
163  start_contour (lvl, ct_x, ct_y);
164  first = false;
165  }
166 
167  // Find stop edge.
168  // FIXME: perhaps this should use a while loop?
169  for (unsigned int k = 1; k <= 4; k++)
170  {
171  if (start_edge == 0 || start_edge == 2)
172  stop_edge = (start_edge + k) % 4;
173  else
174  stop_edge = (start_edge - k) % 4;
175 
176  if (static_cast<char> (1 << stop_edge) & id)
177  break;
178  }
179 
180  pt[0] = stop_edge;
181  pt[1] = (pt[0] + 1) % 4;
182  tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]);
183 
184  if (octave::math::isnan (tmp))
185  ct_x = ct_y = 0.5;
186  else
187  {
188  ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp);
189  ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp);
190  }
191 
192  // Add point to contour.
193  add_point (ct_x, ct_y);
194 
195  // Decrease id value of current facet for start edge.
196  mark(r, c) -= static_cast<char> (1 << stop_edge);
197 
198  // Find next facet.
199  if (stop_edge == 0)
200  r--;
201  else if (stop_edge == 1)
202  c++;
203  else if (stop_edge == 2)
204  r++;
205  else if (stop_edge == 3)
206  c--;
207 
208  // Go to next facet.
209  start_edge = (stop_edge + 2) % 4;
210 
211  }
212 }
213 
214 static void
215 mark_facets (const Matrix& Z, charMatrix& mark, double lvl)
216 {
217  unsigned int nr = mark.rows ();
218  unsigned int nc = mark.cols ();
219 
220  double f[4];
221 
222  for (unsigned int c = 0; c < nc; c++)
223  for (unsigned int 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 (unsigned int i = 0; i < 4; i++)
231  if (fabs(f[i]) < std::numeric_limits<double>::epsilon ())
232  f[i] = std::numeric_limits<double>::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  }
240 
241  for (unsigned int r = 0; r < nr; r++)
242  for (unsigned int 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 (unsigned int i = 0; i < 4; i++)
250  if (fabs(f[i]) < std::numeric_limits<double>::epsilon ())
251  f[i] = std::numeric_limits<double>::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  unsigned int nr = Z.rows ();
265  unsigned int 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 (unsigned int c = 0; c < nc - 1; c++)
274  {
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);
282  }
283 
284  for (unsigned int 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 (unsigned int r = 0; r < nr - 1; r++)
296  for (unsigned int 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);
299 }
300 
301 DEFUN (__contourc__, args, ,
302  doc: /* -*- texinfo -*-
303 @deftypefn {} {} __contourc__ (@var{x}, @var{y}, @var{z}, @var{levels})
304 Undocumented internal function.
305 @end deftypefn */)
306 {
307  if (args.length () != 4)
308  print_usage ();
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  contourc.resize (2, 0);
316 
317  for (int i = 0; i < L.numel (); i++)
318  cntr (X, Y, Z, L (i));
319 
320  end_contour ();
321 
322  return octave_value (contourc);
323 }
324 
325 /*
326 ## No test needed for internal helper function.
327 %!assert (1)
328 */
uint32_t id
Definition: graphics.cc:11587
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:145
static Matrix this_contour
Definition: __contourc__.cc:48
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:363
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE * f
bool isnan(double x)
Definition: lo-mappers.cc:347
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup so it is best to provide a consistent set and leave this option set to zero tem it may help to set this parameter to a nonzero value it is probably best to try leaving this option set to zero first
Definition: DASSL-opts.cc:419
for large enough k
Definition: lu.cc:606
Matrix extract_n(octave_idx_type r1, octave_idx_type c1, octave_idx_type nr, octave_idx_type nc) const
Definition: dMatrix.cc:408
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
octave_idx_type rows(void) const
Definition: Array.h:401
static void start_contour(double lvl, double x, double y)
Definition: __contourc__.cc:87
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup calculate Y_a and Y _d item Given Y
Definition: DASPK-opts.cc:739
JNIEnv void * args
Definition: ov-java.cc:67
#define CONTOUR_QUANT
Definition: __contourc__.cc:53
static int elem
Definition: __contourc__.cc:50
double tmp
Definition: data.cc:6300
static void cntr(const RowVector &X, const RowVector &Y, const Matrix &Z, double lvl)
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Z
Definition: dMatrix.h:37
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
Definition: lu.cc:138
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)
Definition: __contourc__.cc:96
static void end_contour(void)
Definition: __contourc__.cc:72
static Matrix contourc
Definition: __contourc__.cc:49
static void mark_facets(const Matrix &Z, charMatrix &mark, double lvl)
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
the element is set to zero In other the statement xample y
Definition: data.cc:5342
static void add_point(double x, double y)
Definition: __contourc__.cc:58
octave_idx_type cols(void) const
Definition: Array.h:409
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))
Matrix append(const Matrix &a) const
Definition: dMatrix.cc:266