Mercurial > octave
view libinterp/corefcn/__contourc__.cc @ 30564:796f54d4ddbf stable
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2021.
In all .txi and .texi files except gpl.txi and gpl.texi in the
doc/liboctave and doc/interpreter directories, change the copyright
to "Octave Project Developers", the same as used for other source
files. Update copyright notices for 2022 (not done since 2019). For
gpl.txi and gpl.texi, change the copyright notice to be "Free Software
Foundation, Inc." and leave the date at 2007 only because this file
only contains the text of the GPL, not anything created by the Octave
Project Developers.
Add Paul Thomas to contributors.in.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Tue, 28 Dec 2021 18:22:40 -0500 |
parents | 7d6709900da7 |
children | e88a07dec498 |
line wrap: on
line source
//////////////////////////////////////////////////////////////////////// // // Contour lines for function evaluated on a grid. // // Copyright (C) 2001-2022 The Octave Project Developers // // See the file COPYRIGHT.md in the top-level directory of this // distribution or <https://octave.org/copyright/>. // // Adapted to an oct file from the stand alone contourl by Victro Munoz // Copyright (C) 2004 Victor Munoz // // Based on contour plot routine (plcont.c) in PLPlot package // http://plplot.org/ // // Copyright (C) 1995, 2000, 2001 Maurice LeBrun // Copyright (C) 2000, 2002 Joao Cardoso // Copyright (C) 2000, 2001, 2002, 2004 Alan W. Irwin // Copyright (C) 2004 Andrew Ross // // // This file is part of Octave. // // Octave is free software: you can redistribute it and/or modify it // under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // Octave is distributed in the hope that it will be useful, but // WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with Octave; see the file COPYING. If not, see // <https://www.gnu.org/licenses/>. // //////////////////////////////////////////////////////////////////////// #if defined (HAVE_CONFIG_H) # include "config.h" #endif #include <limits> #include "defun.h" #include "ov.h" OCTAVE_NAMESPACE_BEGIN // FIXME: this looks like trouble... static Matrix this_contour; static Matrix contourc; static int elem; // This is the quanta in which we increase this_contour. #define CONTOUR_QUANT 50 // Add a coordinate point (x,y) to this_contour. static void add_point (double x, double y) { if (elem % CONTOUR_QUANT == 0) this_contour = this_contour.append (Matrix (2, CONTOUR_QUANT, 0)); this_contour(0, elem) = x; this_contour(1, elem) = y; elem++; } // Add contents of current contour to contourc. // this_contour.cols () - 1; static void end_contour (void) { if (elem > 2) { this_contour(1, 0) = elem - 1; contourc = contourc.append (this_contour.extract_n (0, 0, 2, elem)); } this_contour = Matrix (); elem = 0; } // Start a new contour, and add contents of current one to contourc. static void start_contour (double lvl, double x, double y) { end_contour (); this_contour.resize (2, 0); add_point (lvl, 0); add_point (x, y); } 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) { double px[4], py[4], pz[4], tmp; unsigned int stop_edge, pt[2]; // Continue while next facet is not done yet. while (r >= 0 && c >= 0 && r < mark.rows () && c < mark.cols () && mark(r, c) > 0) { //get x, y, and z - lvl for current facet px[0] = px[3] = X(c); px[1] = px[2] = X(c+1); py[0] = py[1] = Y(r); py[2] = py[3] = Y(r+1); pz[3] = Z(r+1, c) - lvl; pz[2] = Z(r+1, c + 1) - lvl; pz[1] = Z(r, c+1) - lvl; pz[0] = Z(r, c) - lvl; // Facet edge and point naming assignment. // // 0-----1 .-0-. // | | | | // | | 3 1 // | | | | // 3-----2 .-2-. // Get mark value of current facet. char id = static_cast<char> (mark(r, c)); // Check startedge s. if (start_edge == 255) { // Find start edge. for (unsigned int k = 0; k < 4; k++) if (static_cast<char> (1 << k) & id) start_edge = k; } if (start_edge == 255) break; // Decrease mark value of current facet for start edge. mark(r, c) -= static_cast<char> (1 << start_edge); // Next point (clockwise). pt[0] = start_edge; pt[1] = (pt[0] + 1) % 4; // Calculate contour segment start if first of contour. if (first) { tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]); if (math::isnan (tmp)) ct_x = ct_y = 0.5; else { ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp); ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp); } start_contour (lvl, ct_x, ct_y); first = false; } // Find stop edge. // FIXME: perhaps this should use a while loop? for (unsigned int k = 1; k <= 4; k++) { if (start_edge == 0 || start_edge == 2) stop_edge = (start_edge + k) % 4; else stop_edge = (start_edge - k) % 4; if (static_cast<char> (1 << stop_edge) & id) break; } pt[0] = stop_edge; pt[1] = (pt[0] + 1) % 4; tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]); if (math::isnan (tmp)) ct_x = ct_y = 0.5; else { ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp); ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp); } // Add point to contour. add_point (ct_x, ct_y); // Decrease id value of current facet for start edge. mark(r, c) -= static_cast<char> (1 << stop_edge); // Find next facet. if (stop_edge == 0) r--; else if (stop_edge == 1) c++; else if (stop_edge == 2) r++; else if (stop_edge == 3) c--; // Go to next facet. start_edge = (stop_edge + 2) % 4; } } static void mark_facets (const Matrix& Z, charMatrix& mark, double lvl) { unsigned int nr = mark.rows (); unsigned int nc = mark.cols (); double f[4]; for (unsigned int c = 0; c < nc; c++) for (unsigned int r = 0; r < nr; r++) { f[0] = Z(r, c) - lvl; f[1] = Z(r, c+1) - lvl; f[3] = Z(r+1, c) - lvl; f[2] = Z(r+1, c+1) - lvl; for (unsigned int i = 0; i < 4; i++) if (fabs(f[i]) < std::numeric_limits<double>::epsilon ()) f[i] = std::numeric_limits<double>::epsilon (); if (f[1] * f[2] < 0) mark(r, c) += 2; if (f[0] * f[3] < 0) mark(r, c) += 8; } for (unsigned int r = 0; r < nr; r++) for (unsigned int c = 0; c < nc; c++) { f[0] = Z(r, c) - lvl; f[1] = Z(r, c+1) - lvl; f[3] = Z(r+1, c) - lvl; f[2] = Z(r+1, c+1) - lvl; for (unsigned int i = 0; i < 4; i++) if (fabs(f[i]) < std::numeric_limits<double>::epsilon ()) f[i] = std::numeric_limits<double>::epsilon (); if (f[0] * f[1] < 0) mark(r, c) += 1; if (f[2] * f[3] < 0) mark(r, c) += 4; } } static void cntr (const RowVector& X, const RowVector& Y, const Matrix& Z, double lvl) { unsigned int nr = Z.rows (); unsigned int nc = Z.cols (); charMatrix mark (nr - 1, nc - 1, 0); mark_facets (Z, mark, lvl); // Find contours that start at a domain edge. for (unsigned int c = 0; c < nc - 1; c++) { // Top. if (mark(0, c) & 1) drawcn (X, Y, Z, lvl, 0, c, 0.0, 0.0, 0, true, mark); // Bottom. if (mark(nr - 2, c) & 4) drawcn (X, Y, Z, lvl, nr - 2, c, 0.0, 0.0, 2, true, mark); } for (unsigned int r = 0; r < nr - 1; r++) { // Left. if (mark(r, 0) & 8) drawcn (X, Y, Z, lvl, r, 0, 0.0, 0.0, 3, true, mark); // Right. if (mark(r, nc - 2) & 2) drawcn (X, Y, Z, lvl, r, nc - 2, 0.0, 0.0, 1, true, mark); } for (unsigned int r = 0; r < nr - 1; r++) for (unsigned int c = 0; c < nc - 1; c++) if (mark (r, c) > 0) drawcn (X, Y, Z, lvl, r, c, 0.0, 0.0, 255, true, mark); } DEFUN (__contourc__, args, , doc: /* -*- texinfo -*- @deftypefn {} {@var{c} =} __contourc__ (@var{x}, @var{y}, @var{z}, @var{levels}) Calculate Z-level contours (isolines). @end deftypefn */) { RowVector X = args(0).row_vector_value (); RowVector Y = args(1).row_vector_value (); Matrix Z = args(2).matrix_value (); RowVector L = args(3).row_vector_value (); contourc.resize (2, 0); for (int i = 0; i < L.numel (); i++) cntr (X, Y, Z, L(i)); end_contour (); return octave_value (contourc); } /* ## No test needed for internal helper function. %!assert (1) */ OCTAVE_NAMESPACE_END