////////////////////////////////////////////////////////////////////////
//
// Contour lines for function evaluated on a grid.
//
// Copyright (C) 2001-2021 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