IsoLine.cpp

Go to the documentation of this file.

/**********************************************************************
zyGrib: meteorological GRIB file viewer
Copyright (C) 2008 - Jacques Zaninetti - http://www.zygrib.org
 
This program 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.
 
This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#include "pi_gl.h"  // Must included before anything using GL stuff
 
#include "wx/wxprec.h"
 
#ifndef WX_PRECOMP
#include "wx/wx.h"
#endif  // precompiled headers
 
// #include "chcanv.h"
// #include "model/georef.h"
#include <wx/graphics.h>
 
#include "IsoLine.h"
#include "GribSettingsDialog.h"
#include "GribOverlayFactory.h"
 
#ifdef __OCPN__ANDROID__
#include "qdebug.h"
#endif
 
// static void GenerateSpline(int n, wxPoint points[]);
// static void ClearSplineList();
wxList ocpn_wx_spline_point_list;
 
#include <wx/listimpl.cpp>
WX_DEFINE_LIST(MySegList);
WX_DEFINE_LIST(MySegListList);
 
#ifndef PI
#define PI 3.14159
#endif
#define CTRUE -1
#define CFALSE 0
 
/* Local variables for cohen_sutherland_line_clip: */
struct LOC_cohen_sutherland_line_clip {
  double xmin, xmax, ymin, ymax;
};
void CompOutCode(double x, double y, outcode *code,
                 struct LOC_cohen_sutherland_line_clip *LINK) {
  /*Compute outcode for the point (x,y) */
  *code = 0;
  if (y > LINK->ymax)
    *code = 1L << ((long)TOP);
  else if (y < LINK->ymin)
    *code = 1L << ((long)BOTTOM);
  if (x > LINK->xmax)
    *code |= 1L << ((long)RIGHT);
  else if (x < LINK->xmin)
    *code |= 1L << ((long)LEFT);
}
 
ClipResult cohen_sutherland_line_clip_d(double *x0, double *y0, double *x1,
                                        double *y1, double xmin_, double xmax_,
                                        double ymin_, double ymax_) {
  /* Cohen-Sutherland clipping algorithm for line P0=(x1,y0) to P1=(x1,y1)
   *    and clip rectangle with diagonal from (xmin,ymin) to (xmax,ymax).*/
  struct LOC_cohen_sutherland_line_clip V;
  int done = CFALSE;
  ClipResult clip = Visible;
  outcode outcode0, outcode1, outcodeOut;
  /*Outcodes for P0,P1, and whichever point lies outside the clip rectangle*/
  double x = 0., y = 0.;
 
  V.xmin = xmin_;
  V.xmax = xmax_;
  V.ymin = ymin_;
  V.ymax = ymax_;
  CompOutCode(*x0, *y0, &outcode0, &V);
  CompOutCode(*x1, *y1, &outcode1, &V);
  do {
    if (outcode0 == 0 && outcode1 == 0) { /*Trivial accept and exit*/
      done = CTRUE;
    } else if ((outcode0 & outcode1) != 0) {
      clip = Invisible;
      done = CTRUE;
    }
    /*Logical intersection is true, so trivial reject and exit.*/
    else {
      clip = Visible;
      /*Failed both tests, so calculate the line segment to clip;
       *            from an outside point to an intersection with clip edge.*/
      /*At least one endpoint is outside the clip rectangle; pick it.*/
      if (outcode0 != 0)
        outcodeOut = outcode0;
      else
        outcodeOut = outcode1;
      /*Now find intersection point;
       *            use formulas y=y0+slope*(x-x0),x=x0+(1/slope)*(y-y0).*/
 
      if (((1L << ((long)TOP)) & outcodeOut) != 0) {
        /*Divide line at top of clip rectangle*/
        x = *x0 + (*x1 - *x0) * (V.ymax - *y0) / (*y1 - *y0);
        y = V.ymax;
      } else if (((1L << ((long)BOTTOM)) & outcodeOut) != 0) {
        /*Divide line at bottom of clip rectangle*/
        x = *x0 + (*x1 - *x0) * (V.ymin - *y0) / (*y1 - *y0);
        y = V.ymin;
      } else if (((1L << ((long)RIGHT)) & outcodeOut) != 0) {
        /*Divide line at right edge of clip rectangle*/
        y = *y0 + (*y1 - *y0) * (V.xmax - *x0) / (*x1 - *x0);
        x = V.xmax;
      } else if (((1L << ((long)LEFT)) & outcodeOut) != 0) {
        /*Divide line at left edge of clip rectangle*/
        y = *y0 + (*y1 - *y0) * (V.xmin - *x0) / (*x1 - *x0);
        x = V.xmin;
      }
      /*Now we move outside point to intersection point to clip,
       *            and get ready for next pass.*/
      if (outcodeOut == outcode0) {
        *x0 = x;
        *y0 = y;
        CompOutCode(*x0, *y0, &outcode0, &V);
      } else {
        *x1 = x;
        *y1 = y;
        CompOutCode(*x1, *y1, &outcode1, &V);
      }
    }
  } while (!done);
  return clip;
}
 
ClipResult cohen_sutherland_line_clip_i(int *x0_, int *y0_, int *x1_, int *y1_,
                                        int xmin_, int xmax_, int ymin_,
                                        int ymax_) {
  ClipResult ret;
  double x0, y0, x1, y1;
  x0 = *x0_;
  y0 = *y0_;
  x1 = *x1_;
  y1 = *y1_;
  ret =
      cohen_sutherland_line_clip_d(&x0, &y0, &x1, &y1, (double)xmin_,
                                   (double)xmax_, (double)ymin_, (double)ymax_);
  *x0_ = (int)x0;
  *y0_ = (int)y0;
  *x1_ = (int)x1;
  *y1_ = (int)y1;
  return ret;
}
 
double round_msvc(double x) { return (floor(x + 0.5)); }
 
//---------------------------------------------------------------
IsoLine::IsoLine(double val, double coeff, double offset,
                 const GribRecord *rec_) {
  if (wxGetDisplaySize().x > 0) {
    m_pixelMM = PlugInGetDisplaySizeMM() / wxGetDisplaySize().x;
    m_pixelMM = wxMax(.02, m_pixelMM);  // protect against bad data
  } else
    m_pixelMM = 0.27;  // semi-standard number...
 
  value = val / coeff - offset;
 
  rec = rec_;
  W = rec_->getNi();
  H = rec_->getNj();
 
  //---------------------------------------------------------
  // Génère la liste des segments.
  extractIsoLine(rec_);
 
  value = val;
 
  if (trace.size() == 0) return;
 
  //      Join the isoline segments into a nice list
  //      Which is end-to-end continuous and unidirectional
 
  //      Create a master wxList of the trace list
  std::list<Segment *>::iterator it;
  for (it = trace.begin(); it != trace.end(); it++) {
    Segment *seg = *it;
    seg->bUsed = false;
    m_seglist.Append(*it);
  }
 
  //      Isoline may be discontinuous....
  //      So build a list of continuous segments
  bool bdone = false;
  while (!bdone) {
    MySegList *ps = BuildContinuousSegment();
 
    m_SegListList.Append(ps);
 
    MySegList::Node *node;
    Segment *seg;
 
    // recreate the master list, removing used segs
 
    node = m_seglist.GetFirst();
    while (node) {
      seg = node->GetData();
      if (seg->bUsed) {
        m_seglist.Erase(node);
        node = m_seglist.GetFirst();
      } else
        node = node->GetNext();
    }
 
    if (0 == m_seglist.GetCount()) bdone = true;
  }
 
}
//---------------------------------------------------------------
IsoLine::~IsoLine() {
  // printf("delete Isobar : press=%4.0f long=%d\n", pressure/100,
  // trace.size());
 
  std::list<Segment *>::iterator it;
  for (it = trace.begin(); it != trace.end(); it++) {
    delete *it;
    *it = nullptr;
  }
  trace.clear();
 
  m_SegListList.DeleteContents(true);
  m_SegListList.Clear();
}
 
MySegList *IsoLine::BuildContinuousSegment(void) {
  MySegList::Node *node;
  Segment *seg;
 
  MySegList *ret_list = new MySegList;
 
  //     Build a chain extending from the "2" end of the target segment
  //      The joined list, side 2...
  MySegList segjoin2;
 
  //      Add any first segment to the list
  node = m_seglist.GetFirst();
  Segment *seg0 = node->GetData();
  seg0->bUsed = true;
  segjoin2.Append(seg0);
 
  Segment *tseg = seg0;
 
  while (tseg) {
    bool badded = false;
    Segment *seg;
    node = m_seglist.GetFirst();
    while (node) {
      seg = node->GetData();
 
      if ((!seg->bUsed) && (seg->py1 == tseg->py2) &&
          (seg->px1 == tseg->px2))  // fits without reverse
      {
        seg->bUsed = true;
        segjoin2.Append(seg);
        badded = true;
        break;
      } else if ((!seg->bUsed) && (seg->py2 == tseg->py2) &&
                 (seg->px2 == tseg->px2))  // fits, needs reverse
      {
        seg->bUsed = true;
        double a = seg->px2;
        seg->px2 = seg->px1;
        seg->px1 = a;
        double b = seg->py2;
        seg->py2 = seg->py1;
        seg->py1 = b;
        segjoin2.Append(seg);
        badded = true;
        break;
      }
 
      node = node->GetNext();
    }
    if (badded == true)
      tseg = seg;
    else
      tseg = nullptr;
  }
 
  //     Build a chain extending from the "1" end of the target segment
  //      The joined list, side 1...
  MySegList segjoin1;
 
  //      Add the same first segment to the list
  node = m_seglist.GetFirst();
  seg0 = node->GetData();
  seg0->bUsed = true;
  segjoin1.Append(seg0);
 
  tseg = seg0;
 
  while (tseg) {
    bool badded = false;
    node = m_seglist.GetFirst();
    while (node) {
      seg = node->GetData();
 
      if ((!seg->bUsed) && (seg->py2 == tseg->py1) &&
          (seg->px2 == tseg->px1))  // fits without reverse
      {
        seg->bUsed = true;
        segjoin1.Append(seg);
        badded = true;
        break;
      } else if ((!seg->bUsed) && (seg->py1 == tseg->py1) &&
                 (seg->px1 == tseg->px1))  // fits, needs reverse
      {
        seg->bUsed = true;
        double a = seg->px2;
        seg->px2 = seg->px1;
        seg->px1 = a;
        double b = seg->py2;
        seg->py2 = seg->py1;
        seg->py1 = b;
        segjoin1.Append(seg);
        badded = true;
        break;
      }
 
      node = node->GetNext();
    }
    if (badded == true)
      tseg = seg;
    else
      tseg = nullptr;
  }
 
  //     Now have two lists...
 
  //     Start with "1" side list,
  //    starting from the end, and skipping the first segment
 
  int n1 = segjoin1.GetCount();
  for (int i = n1 - 1; i > 0; i--) {
    node = segjoin1.Item(i);
    seg = node->GetData();
    ret_list->Append(seg);
  }
 
  //     Now add the "2"side list
  int n2 = segjoin2.GetCount();
  for (int i = 0; i < n2; i++) {
    node = segjoin2.Item(i);
    seg = node->GetData();
    ret_list->Append(seg);
  }
 
  //     And there it is
 
  return ret_list;
}
 
//---------------------------------------------------------------
void IsoLine::drawIsoLine(GRIBOverlayFactory *pof, wxDC *dc,
                          PlugIn_ViewPort *vp, bool bHiDef) {
  int nsegs = trace.size();
  if (nsegs < 1) return;
 
  GetGlobalColor(_T ( "UITX1" ), &isoLineColor);
 
#if wxUSE_GRAPHICS_CONTEXT
  wxGraphicsContext *pgc = nullptr;
#endif
 
  if (dc) {
    wxPen ppISO(isoLineColor, 2);
 
#if wxUSE_GRAPHICS_CONTEXT
    wxMemoryDC *pmdc;
    pmdc = dynamic_cast<wxMemoryDC *>(dc);
    pgc = wxGraphicsContext::Create(*pmdc);
    pgc->SetPen(ppISO);
#endif
    dc->SetPen(ppISO);
  } else { /* opengl */
#ifdef ocpnUSE_GL
    // #ifndef USE_ANDROID_GLES2
    //            if(m_pixelMM > 0.2){        // pixel size large enough to
    //            render well
    //            //      Enable anti-aliased lines, at best quality
    //              glEnable( GL_LINE_SMOOTH );
    //              glEnable( GL_BLEND );
    //              glBlendFunc( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA );
    //              glHint( GL_LINE_SMOOTH_HINT, GL_NICEST );
    //              glLineWidth( 2 );
    //            }
    //            else{
    //              glLineWidth( 0.4/m_pixelMM);        //  set a target line
    //              width by MM
    //            }
    // #else
    if (pof->m_oDC) {
      wxPen ppISO(isoLineColor, 2);
      pof->m_oDC->SetPen(ppISO);
    }
#endif
  }
 
  std::list<Segment *>::iterator it;
 
  //---------------------------------------------------------
  // Dessine les segments
  //---------------------------------------------------------
  for (it = trace.begin(); it != trace.end(); it++) {
    Segment *seg = *it;
 
    if (vp->m_projection_type == PI_PROJECTION_MERCATOR ||
        vp->m_projection_type == PI_PROJECTION_EQUIRECTANGULAR) {
      /* skip segments that go the wrong way around the world */
      double sx1 = seg->px1, sx2 = seg->px2;
      if (sx2 - sx1 > 180)
        sx2 -= 360;
      else if (sx1 - sx2 > 180)
        sx1 -= 360;
 
      if ((sx1 + 180 < vp->clon && sx2 + 180 > vp->clon) ||
          (sx1 + 180 > vp->clon && sx2 + 180 < vp->clon) ||
          (sx1 - 180 < vp->clon && sx2 - 180 > vp->clon) ||
          (sx1 - 180 > vp->clon && sx2 - 180 < vp->clon))
        continue;
    }
 
    wxPoint ab;
    GetCanvasPixLL(vp, &ab, seg->py1, seg->px1);
    wxPoint cd;
    GetCanvasPixLL(vp, &cd, seg->py2, seg->px2);
 
    if (dc) {
#if wxUSE_GRAPHICS_CONTEXT
      if (bHiDef && pgc)
        pgc->StrokeLine(ab.x, ab.y, cd.x, cd.y);
      else
#endif
        dc->DrawLine(ab.x, ab.y, cd.x, cd.y);
    } else { /* opengl */
#ifdef ocpnUSE_GL
 
      if (pof->m_oDC) {
        pof->m_oDC->DrawLine(ab.x, ab.y, cd.x, cd.y);
      }
 
#endif
    }
  }
 
#if wxUSE_GRAPHICS_CONTEXT
  delete pgc;
#endif
 
  //      if(!dc) /* opengl */
  //          glEnd();
}
 
//---------------------------------------------------------------
 
void IsoLine::drawIsoLineLabels(GRIBOverlayFactory *pof, wxDC *dc,
                                PlugIn_ViewPort *vp, int density, int first,
                                wxImage &imageLabel)
 
{
  std::list<Segment *>::iterator it;
  int nb = first;
  wxString label;
 
  //---------------------------------------------------------
  // Ecrit les labels
  //---------------------------------------------------------
  wxRect prev;
  for (it = trace.begin(); it != trace.end(); it++, nb++) {
    if (nb % density == 0) {
      Segment *seg = *it;
 
      //            if(vp->vpBBox.PointInBox((seg->px1 + seg->px2)/2., (seg->py1
      //            + seg->py2)/2., 0.))
      {
        wxPoint ab;
        GetCanvasPixLL(vp, &ab, seg->py1, seg->px1);
        wxPoint cd;
        GetCanvasPixLL(vp, &cd, seg->py1, seg->px1);
 
        int w = imageLabel.GetWidth();
        int h = imageLabel.GetHeight();
 
        int label_offset = 6;
        int xd = (ab.x + cd.x - (w + label_offset * 2)) / 2;
        int yd = (ab.y + cd.y - h) / 2;
 
        int x = xd - label_offset;
        wxRect r(x, yd, w, h);
        r.Inflate(w);
        if (!prev.Intersects(r)) {
          prev = r;
 
          /* don't use alpha for isobars, for some reason draw bitmap ignores
             the 4th argument (true or false has same result) */
          wxImage img(w, h, imageLabel.GetData(), true);
          dc->DrawBitmap(img, xd, yd, false);
        }
      }
    }
  }
}
 
void IsoLine::drawIsoLineLabelsGL(GRIBOverlayFactory *pof, PlugIn_ViewPort *vp,
                                  int density, int first, wxString label,
                                  wxColour &color, TexFont &texfont)
 
{
  std::list<Segment *>::iterator it;
  int nb = first;
 
#ifdef ocpnUSE_GL
  glEnable(GL_BLEND);
  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
 
  //---------------------------------------------------------
  // Ecrit les labels
  //---------------------------------------------------------
  wxRect prev;
  for (it = trace.begin(); it != trace.end(); it++, nb++) {
    if (nb % density == 0) {
      Segment *seg = *it;
 
      //            if(vp->vpBBox.PointInBox((seg->px1 + seg->px2)/2., (seg->py1
      //            + seg->py2)/2., 0.))
      {
        wxPoint ab;
        GetCanvasPixLL(vp, &ab, seg->py1, seg->px1);
        wxPoint cd;
        GetCanvasPixLL(vp, &cd, seg->py1, seg->px1);
 
        int w, h;
        texfont.GetTextExtent(label, &w, &h);
 
        int label_offsetx = 6, label_offsety = 1;
        int xd = (ab.x + cd.x - (w + label_offsetx * 2)) / 2;
        int yd = (ab.y + cd.y - h) / 2;
        int x = xd - label_offsetx, y = yd - label_offsety;
        w += 2 * label_offsetx, h += 2 * label_offsety;
 
        wxRect r(x, y, w, h);
        r.Inflate(w);
        if (!prev.Intersects(r)) {
#if 1
          prev = r;
          if (pof->m_oDC) {
            // pof->m_oDC->SetFont( *mfont );
            pof->m_oDC->SetPen(*wxBLACK_PEN);
            pof->m_oDC->SetBrush(color);
            pof->m_oDC->DrawRectangle(x, y, w, h);
            pof->m_oDC->DrawText(label, xd, yd);
          }
 
#else
          prev = r;
          glColor4ub(color.Red(), color.Green(), color.Blue(), color.Alpha());
 
          /* draw bounding rectangle */
          glBegin(GL_QUADS);
          glVertex2i(x, y);
          glVertex2i(x + w, y);
          glVertex2i(x + w, y + h);
          glVertex2i(x, y + h);
          glEnd();
 
          glColor3ub(0, 0, 0);
 
          glBegin(GL_LINE_LOOP);
          glVertex2i(x, y);
          glVertex2i(x + w, y);
          glVertex2i(x + w, y + h);
          glVertex2i(x, y + h);
          glEnd();
 
          glEnable(GL_TEXTURE_2D);
          texfont.RenderString(label, xd, yd);
          glDisable(GL_TEXTURE_2D);
#endif
        }
      }
    }
  }
  glDisable(GL_BLEND);
#endif
}
 
//==================================================================================
// Segment
//==================================================================================
Segment::Segment(int I, int w, int J, char c1, char c2, char c3, char c4,
                 const GribRecord *rec, double pressure) {
  traduitCode(I, w, J, c1, i, j);
  traduitCode(I, w, J, c2, k, l);
  traduitCode(I, w, J, c3, m, n);
  traduitCode(I, w, J, c4, o, p);
 
  intersectionAreteGrille(i, j, k, l, &px1, &py1, rec, pressure);
  intersectionAreteGrille(m, n, o, p, &px2, &py2, rec, pressure);
}
//-----------------------------------------------------------------------
void Segment::intersectionAreteGrille(int i, int j, int k, int l, double *x,
                                      double *y, const GribRecord *rec,
                                      double pressure) {
  double xa, xb, ya, yb, pa, pb, dec;
  pa = rec->getValue(i, j);
  pb = rec->getValue(k, l);
 
  rec->getXY(i, j, &xa, &ya);
  rec->getXY(k, l, &xb, &yb);
 
  // Abscisse
  if (pb != pa)
    dec = (pressure - pa) / (pb - pa);
  else
    dec = 0.5;
  if (fabs(dec) > 1) dec = 0.5;
  double xd = xb - xa;
  if (xd < -180)
    xd += 360;
  else if (xd > 180)
    xd -= 360;
  *x = xa + xd * dec;
 
  // Ordonnée
  if (pb != pa)
    dec = (pressure - pa) / (pb - pa);
  else
    dec = 0.5;
  if (fabs(dec) > 1) dec = 0.5;
  *y = ya + (yb - ya) * dec;
}
//---------------------------------------------------------------
void Segment::traduitCode(int I, int w, int J, char c1, int &i, int &j) {
  int Im1 = I ? I - 1 : w - 1;
  switch (c1) {
    case 'a':
      i = Im1;
      j = J - 1;
      break;
    case 'b':
      i = I;
      j = J - 1;
      break;
    case 'c':
      i = Im1;
      j = J;
      break;
    case 'd':
      i = I;
      j = J;
      break;
    default:
      i = I;
      j = J;
  }
}
 
//-----------------------------------------------------------------------
// Génère la liste des segments.
// Les coordonnées sont les indices dans la grille du GribRecord
//---------------------------------------------------------
void IsoLine::extractIsoLine(const GribRecord *rec) {
  int i, j, W, H;
  double a, b, c, d;
  W = rec->getNi();
  H = rec->getNj();
 
  int We = W;
  if (rec->getLonMax() + rec->getDi() - rec->getLonMin() == 360) We++;
 
  for (j = 1; j < H; j++)  // !!!! 1 to end
  {
    a = rec->getValue(0, j - 1);
    c = rec->getValue(0, j);
    for (i = 1; i < We; i++, a = b, c = d) {
      //            x = rec->getX(i);
      //            y = rec->getY(j);
 
      int ni = i;
      if (i == W) ni = 0;
      b = rec->getValue(ni, j - 1);
      d = rec->getValue(ni, j);
 
      if (a == GRIB_NOTDEF || b == GRIB_NOTDEF || c == GRIB_NOTDEF ||
          d == GRIB_NOTDEF)
        continue;
 
      if ((a < value && b < value && c < value && d < value) ||
          (a > value && b > value && c > value && d > value))
        continue;
 
      // Détermine si 1 ou 2 segments traversent la case ab-cd
      // a  b
      // c  d
      //--------------------------------
      // 1 segment en diagonale
      //--------------------------------
      if ((a <= value && b <= value && c <= value && d > value) ||
          (a > value && b > value && c > value && d <= value))
        trace.push_back(new Segment(ni, W, j, 'c', 'd', 'b', 'd', rec, value));
      else if ((a <= value && c <= value && d <= value && b > value) ||
               (a > value && c > value && d > value && b <= value))
        trace.push_back(new Segment(ni, W, j, 'a', 'b', 'b', 'd', rec, value));
      else if ((c <= value && d <= value && b <= value && a > value) ||
               (c > value && d > value && b > value && a <= value))
        trace.push_back(new Segment(ni, W, j, 'a', 'b', 'a', 'c', rec, value));
      else if ((a <= value && b <= value && d <= value && c > value) ||
               (a > value && b > value && d > value && c <= value))
        trace.push_back(new Segment(ni, W, j, 'a', 'c', 'c', 'd', rec, value));
      //--------------------------------
      // 1 segment H ou V
      //--------------------------------
      else if ((a <= value && b <= value && c > value && d > value) ||
               (a > value && b > value && c <= value && d <= value))
        trace.push_back(new Segment(ni, W, j, 'a', 'c', 'b', 'd', rec, value));
      else if ((a <= value && c <= value && b > value && d > value) ||
               (a > value && c > value && b <= value && d <= value))
        trace.push_back(new Segment(ni, W, j, 'a', 'b', 'c', 'd', rec, value));
      //--------------------------------
      // 2 segments en diagonale
      //--------------------------------
      else if (a <= value && d <= value && c > value && b > value) {
        trace.push_back(new Segment(ni, W, j, 'a', 'b', 'b', 'd', rec, value));
        trace.push_back(new Segment(ni, W, j, 'a', 'c', 'c', 'd', rec, value));
      } else if (a > value && d > value && c <= value && b <= value) {
        trace.push_back(new Segment(ni, W, j, 'a', 'b', 'a', 'c', rec, value));
        trace.push_back(new Segment(ni, W, j, 'b', 'd', 'c', 'd', rec, value));
      }
    }
  }
}
 
// ----------------------------------------------------------------------------
// splines code lifted from wxWidgets
// ----------------------------------------------------------------------------
 
// ----------------------------------- spline code
// ----------------------------------------
 
void ocpn_wx_quadratic_spline(double a1, double b1, double a2, double b2,
                              double a3, double b3, double a4, double b4);
void ocpn_wx_clear_stack();
int ocpn_wx_spline_pop(double *x1, double *y1, double *x2, double *y2,
                       double *x3, double *y3, double *x4, double *y4);
void ocpn_wx_spline_push(double x1, double y1, double x2, double y2, double x3,
                         double y3, double x4, double y4);
static bool ocpn_wx_spline_add_point(double x, double y);
 
#define half(z1, z2) ((z1 + z2) / 2.0)
#define THRESHOLD 5
 
/* iterative version */
 
void ocpn_wx_quadratic_spline(double a1, double b1, double a2, double b2,
                              double a3, double b3, double a4, double b4) {
  double xmid, ymid;
  double x1, y1, x2, y2, x3, y3, x4, y4;
 
  ocpn_wx_clear_stack();
  ocpn_wx_spline_push(a1, b1, a2, b2, a3, b3, a4, b4);
 
  while (ocpn_wx_spline_pop(&x1, &y1, &x2, &y2, &x3, &y3, &x4, &y4)) {
    xmid = (double)half(x2, x3);
    ymid = (double)half(y2, y3);
    if (fabs(x1 - xmid) < THRESHOLD && fabs(y1 - ymid) < THRESHOLD &&
        fabs(xmid - x4) < THRESHOLD && fabs(ymid - y4) < THRESHOLD) {
      ocpn_wx_spline_add_point(x1, y1);
      ocpn_wx_spline_add_point(xmid, ymid);
    } else {
      ocpn_wx_spline_push(xmid, ymid, (double)half(xmid, x3),
                          (double)half(ymid, y3), (double)half(x3, x4),
                          (double)half(y3, y4), x4, y4);
      ocpn_wx_spline_push(x1, y1, (double)half(x1, x2), (double)half(y1, y2),
                          (double)half(x2, xmid), (double)half(y2, ymid), xmid,
                          ymid);
    }
  }
}
 
/* utilities used by spline drawing routines */
 
typedef struct ocpn_wx_spline_stack_struct {
  double x1, y1, x2, y2, x3, y3, x4, y4;
} Stack;
 
#define SPLINE_STACK_DEPTH 20
static Stack ocpn_wx_spline_stack[SPLINE_STACK_DEPTH];
static Stack *ocpn_wx_stack_top;
static int ocpn_wx_stack_count;
 
void ocpn_wx_clear_stack() {
  ocpn_wx_stack_top = ocpn_wx_spline_stack;
  ocpn_wx_stack_count = 0;
}
 
void ocpn_wx_spline_push(double x1, double y1, double x2, double y2, double x3,
                         double y3, double x4, double y4) {
  ocpn_wx_stack_top->x1 = x1;
  ocpn_wx_stack_top->y1 = y1;
  ocpn_wx_stack_top->x2 = x2;
  ocpn_wx_stack_top->y2 = y2;
  ocpn_wx_stack_top->x3 = x3;
  ocpn_wx_stack_top->y3 = y3;
  ocpn_wx_stack_top->x4 = x4;
  ocpn_wx_stack_top->y4 = y4;
  ocpn_wx_stack_top++;
  ocpn_wx_stack_count++;
}
 
int ocpn_wx_spline_pop(double *x1, double *y1, double *x2, double *y2,
                       double *x3, double *y3, double *x4, double *y4) {
  if (ocpn_wx_stack_count == 0) return (0);
  ocpn_wx_stack_top--;
  ocpn_wx_stack_count--;
  *x1 = ocpn_wx_stack_top->x1;
  *y1 = ocpn_wx_stack_top->y1;
  *x2 = ocpn_wx_stack_top->x2;
  *y2 = ocpn_wx_stack_top->y2;
  *x3 = ocpn_wx_stack_top->x3;
  *y3 = ocpn_wx_stack_top->y3;
  *x4 = ocpn_wx_stack_top->x4;
  *y4 = ocpn_wx_stack_top->y4;
  return (1);
}
 
static bool ocpn_wx_spline_add_point(double x, double y) {
  wxPoint *point = new wxPoint;
  point->x = (int)x;
  point->y = (int)y;
  ocpn_wx_spline_point_list.Append((wxObject *)point);
  return true;
}
 
void GenSpline(wxList *points) {
  wxPoint *p;
  double cx1, cy1, cx2, cy2, cx3, cy3, cx4, cy4;
  double x1, y1, x2, y2;
 
  wxList::compatibility_iterator node = points->GetFirst();
  if (!node)
    // empty list
    return;
 
  p = (wxPoint *)node->GetData();
 
  x1 = p->x;
  y1 = p->y;
 
  node = node->GetNext();
  p = (wxPoint *)node->GetData();
 
  x2 = p->x;
  y2 = p->y;
  cx1 = (double)((x1 + x2) / 2);
  cy1 = (double)((y1 + y2) / 2);
  cx2 = (double)((cx1 + x2) / 2);
  cy2 = (double)((cy1 + y2) / 2);
 
  ocpn_wx_spline_add_point(x1, y1);
 
  while ((node = node->GetNext())
#if !wxUSE_STL
         != nullptr
#endif  // !wxUSE_STL
  ) {
    p = (wxPoint *)node->GetData();
    x1 = x2;
    y1 = y2;
    x2 = p->x;
    y2 = p->y;
    cx4 = (double)(x1 + x2) / 2;
    cy4 = (double)(y1 + y2) / 2;
    cx3 = (double)(x1 + cx4) / 2;
    cy3 = (double)(y1 + cy4) / 2;
 
    ocpn_wx_quadratic_spline(cx1, cy1, cx2, cy2, cx3, cy3, cx4, cy4);
 
    cx1 = cx4;
    cy1 = cy4;
    cx2 = (double)(cx1 + x2) / 2;
    cy2 = (double)(cy1 + y2) / 2;
  }
 
  ocpn_wx_spline_add_point(cx1, cy1);
  ocpn_wx_spline_add_point(x2, y2);
}
 
#if 0
static void GenerateSpline(int n, wxPoint points[])
{
      wxList list;
      for (int i =0; i < n; i++)
      {
            list.Append((wxObject*)&points[i]);
      }
 
      GenSpline(&list);
}
 
static void ClearSplineList()
{
      wxList::compatibility_iterator node = ocpn_wx_spline_point_list.GetFirst();
      while (node)
      {
            wxPoint *point = (wxPoint *)node->GetData();
            delete point;
            ocpn_wx_spline_point_list.Erase(node);
            node = ocpn_wx_spline_point_list.GetFirst();
      }
}
#endif