GNU Octave  8.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
__contourc__.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Contour lines for function evaluated on a grid.
4 //
5 // Copyright (C) 2001-2023 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 
50 
51 // FIXME: this looks like trouble...
54 static int elem;
55 
56 // This is the quanta in which we increase this_contour.
57 #define CONTOUR_QUANT 50
58 
59 // Add a coordinate point (x,y) to this_contour.
60 static void
61 add_point (double x, double y)
62 {
63  if (elem % CONTOUR_QUANT == 0)
65 
66  this_contour(0, elem) = x;
67  this_contour(1, elem) = y;
68  elem++;
69 }
70 
71 // Add contents of current contour to contourc.
72 // this_contour.cols () - 1;
73 static void
75 {
76  if (elem > 2)
77  {
78  this_contour(1, 0) = elem - 1;
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
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
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 (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 (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
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
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 
304 DEFUN (__contourc__, args, ,
305  doc: /* -*- texinfo -*-
306 @deftypefn {} {@var{c} =} __contourc__ (@var{x}, @var{y}, @var{z}, @var{levels})
307 Calculate Z-level contours (isolines).
308 @end deftypefn */)
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 */
329 
OCTAVE_END_NAMESPACE(octave)
static void mark_facets(const Matrix &Z, charMatrix &mark, double lvl)
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:98
static void cntr(const RowVector &X, const RowVector &Y, const Matrix &Z, double lvl)
static void add_point(double x, double y)
Definition: __contourc__.cc:61
static Matrix contourc
Definition: __contourc__.cc:53
static Matrix this_contour
Definition: __contourc__.cc:52
static void start_contour(double lvl, double x, double y)
Definition: __contourc__.cc:89
static void end_contour(void)
Definition: __contourc__.cc:74
static int elem
Definition: __contourc__.cc:54
#define CONTOUR_QUANT
Definition: __contourc__.cc:57
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:414
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type rows(void) const
Definition: Array.h:459
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type cols(void) const
Definition: Array.h:469
Definition: dMatrix.h:42
OCTAVE_API Matrix append(const Matrix &a) const
Definition: dMatrix.cc:265
OCTAVE_API Matrix extract_n(octave_idx_type r1, octave_idx_type c1, octave_idx_type nr, octave_idx_type nc) const
Definition: dMatrix.cc:407
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:158
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:56
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
bool isnan(bool)
Definition: lo-mappers.h:178
F77_RET_T const F77_DBLE * x
F77_RET_T const F77_DBLE const F77_DBLE * f
class OCTAVE_API Matrix
Definition: mx-fwd.h:31
T * r
Definition: mx-inlines.cc:773
return octave_value(v1.char_array_value() . concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string()) ? '\'' :'"'))