// --------------------------------------------------------------------
// Ipelet for computing various Voronoi diagrams
// --------------------------------------------------------------------
/*

    This file is part of the extensible drawing editor Ipe.
    Copyright (C) 1993-2004  Otfried Cheong

    Ipe 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 2 of the License, or
    (at your option) any later version.

    As a special exception, you have permission to link Ipe with the
    CGAL library and distribute executables, as long as you follow the
    requirements of the Gnu General Public License in regard to all of
    the software in the executable aside from CGAL.

    Ipe 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 Ipe; if not, you can find it at
    "http://www.gnu.org/copyleft/gpl.html", or write to the Free
    Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

*/

#include <stdio.h>
#include <stdlib.h>

extern "C" {
#include "qhull_a.h"
}

#include "ipelib.h"

#ifdef __BORLANDC__
#pragma warn -8060
#endif

// --------------------------------------------------------------------

struct DelaunayEdge {
public:
  DelaunayEdge(int x, int y) : a(x), b(y) { /* nothing */ }
  DelaunayEdge() : a(-1), b(-1) { /* nothing */ }
public:
  int a,b;
};

inline bool operator<(const DelaunayEdge& x, const DelaunayEdge& y)
{
  return (x.a > y.a || (x.a == y.a && x.b > y.b));
}

inline bool operator!=(const DelaunayEdge& x, const DelaunayEdge& y)
{
  return (x.a != y.a || x.b != y.b);
}

// --------------------------------------------------------------------

class VoronoiIpelet : public Ipelet {
public:
  VoronoiIpelet();
  virtual int IpelibVersion() const { return IPELIB_VERSION; }
  virtual int NumFunctions() const { return 6; }
  virtual const char *Label() const { return "Voronoi diagrams"; }
  virtual const char *SubLabel(int function) const;
  virtual void Run(int function, IpePage *page, IpeletHelper *helper);
private:
  void AddVoronoiEdge(facetT *facet, facetT *neighbor);
  void AddInfiniteEdge(facetT *facet, facetT *neighbor);
  void VoronoiTreatFacet(facetT *facet);
  void AddDelaunayEdge(int from, int to);
  void DelaunayTreatFacet(facetT *facet);
private:
  int iVoronoiSign;
  std::vector<IpeVector> iSites;
  std::vector<IpeSegment> iEdges;
  std::vector<DelaunayEdge> iDelaunay;
  double iInfiniteEdgeLength;
};

VoronoiIpelet::VoronoiIpelet()
{
  iInfiniteEdgeLength = 100.0;
}

// --------------------------------------------------------------------

const char * const sublabel[] = {
  "Delaunay triangulation",
  "Voronoi diagram",
  "Order-2 Voronoi diagram",
  "Order-3 Voronoi diagram",
  "Furthest-point Voronoi diagram",
  "Set length of infinite edges"
};
// "Medial axis of polygon"

const char *VoronoiIpelet::SubLabel(int function) const
{
  return sublabel[function];
}

// --------------------------------------------------------------------

class CollectVisitor : public IpeVisitor {
public:
  CollectVisitor(std::vector<IpeVector> &sites);
  virtual void VisitGroup(const IpeGroup *obj);
  virtual void VisitPath(const IpePath *obj);
  virtual void VisitMark(const IpeMark *obj);
  virtual void VisitReference(const IpeReference *obj);
private:
  std::vector<IpeVector> &iSites;
  std::list<IpeMatrix> iStack;
};

CollectVisitor::CollectVisitor(std::vector<IpeVector> &sites)
  : iSites(sites)
{
  iStack.push_back(IpeMatrix()); // id matrix
}

void CollectVisitor::VisitGroup(const IpeGroup *obj)
{
  iStack.push_back(iStack.back() * obj->Matrix());
  for (IpeGroup::const_iterator it = obj->begin(); it != obj->end(); ++it)
    (*it)->Accept(*this);
  iStack.pop_back();
}

void CollectVisitor::VisitPath(const IpePath *obj)
{
  IpeMatrix m = iStack.back() * obj->Matrix();
  for (int i = 0; i < obj->NumSubPaths(); ++i) {
    const IpeSegmentSubPath *sp = obj->SubPath(i)->AsSegs();
    if (sp) {
      iSites.push_back(m * sp->Segment(0).CP(0));
      for (int j = 0; j < sp->NumSegments(); ++j) {
	IpePathSegment seg = sp->Segment(j);
	iSites.push_back(m * seg.CP(seg.NumCP() - 1));
      }
    }
  }
}

void CollectVisitor::VisitMark(const IpeMark *obj)
{
  iSites.push_back(iStack.back() * obj->Matrix() * obj->Position());
}

void CollectVisitor::VisitReference(const IpeReference *obj)
{
  if (obj->Object()) {
    iStack.push_back(iStack.back() * obj->Matrix());
    obj->Object()->Accept(*this);
    iStack.pop_back();
  }
}

// --------------------------------------------------------------------

//
// readpoints put points into structure for qhull
//
// returns:
//  number of points, array of point coordinates, ismalloc True
//

static coordT *readpoints(const std::vector<IpeVector> &sites, int mode,
			  int *numpoints, boolT *ismalloc)
{
  coordT *points, *coords;

  if (mode < 2 || mode > 3)
    *numpoints= sites.size();
  else if (mode == 2)
    *numpoints = sites.size() * (sites.size() - 1) / 2;
  else if (mode == 3)
    *numpoints = sites.size() * (sites.size() - 1) * (sites.size() - 2) / 6;

  qh normal_size= 3 * sizeof(coordT); /* for tracing with qh_printpoint */
  *ismalloc= True; 		      /* use malloc since memory not setup */
  coords=points=(coordT*)malloc(*numpoints*3*sizeof(coordT));

  if (mode < 2 || mode == 4) {
    for (uint i = 0; i < sites.size(); i++) {
      *(coords++)= sites[i].iX;
      *(coords++)= sites[i].iY;
      *(coords++)= sites[i].iX * sites[i].iX + sites[i].iY * sites[i].iY;
    }
  } else if (mode == 2) {
    for (uint i = 0; i < sites.size() - 1; i++) {
      for (uint j = i + 1; j < sites.size(); j++) {
	*(coords++)= (sites[i].iX + sites[j].iX) / 2.0;
	*(coords++)= (sites[i].iY + sites[j].iY) / 2.0;
	*(coords++)= (sites[i].iX * sites[i].iX + sites[i].iY * sites[i].iY
	  + sites[j].iX * sites[j].iX + sites[j].iY * sites[j].iY) / 2.0;
      }
    }
  } else if (mode == 3) {
    for (uint i = 0; i < sites.size() - 2; i++) {
      for (uint j = i + 1; j < sites.size() - 1; j++) {
	for (uint k = j + 1; k < sites.size(); k++) {
	  *(coords++)= (sites[i].iX + sites[j].iX + sites[k].iX) / 3.0;
	  *(coords++)= (sites[i].iY + sites[j].iY + sites[k].iY) / 3.0;
	  *(coords++)=
	    (sites[i].iX * sites[i].iX + sites[i].iY * sites[i].iY
	     + sites[j].iX * sites[j].iX + sites[j].iY * sites[j].iY
	     + sites[k].iX * sites[k].iX + sites[k].iY * sites[k].iY)
	    / 3.0;
	}
      }
    }
  } else if (mode == 5) {
#if 0
    // for medial axis
    int j = sites.size() - 1;
    pl_unitvec d;
    for (int i = 0; i < sites.size(); i++) {
      d = normalized(sites[i] - sites[j]).normal();
      // facet has equation (d.x, d.y, 1) * (x,y,z) = c = dot(d, sites[i]);
      // dualize it to point (-d.x/2, -d.y/2, -c)
      *(coords++)= -d.iX / 2.0;
      *(coords++)= -d.iY / 2.0;
      *(coords++)= -dot(d, sites[i]);
      j = i;
    }
#endif
  }
  return points;
}

// --------------------------------------------------------------------

// returns Voronoi vertex dual to facet
inline IpeVector voronoi_vertex(facetT *facet)
{
  return IpeVector(-0.5 * facet->normal[0]/facet->normal[2],
		   -0.5 * facet->normal[1]/facet->normal[2]);
}

void VoronoiIpelet::AddVoronoiEdge(facetT *facet, facetT *neighbor)
{
  if (facet->id < neighbor->id) {
    iEdges.push_back(IpeSegment(voronoi_vertex(facet), voronoi_vertex(neighbor)));
  }
}

void VoronoiIpelet::AddInfiniteEdge(facetT *facet, facetT *neighbor)
{
  IpeVector dir;
  IpeVector v = voronoi_vertex(facet);

  if (neighbor->normal[2] == 0.0) {
    // neighboring facet is vertical
    dir = IpeVector(neighbor->normal[0], neighbor->normal[2]);
  } else {
    dir = v - voronoi_vertex(neighbor);
  }
  dir = dir.Normalized();
  iEdges.push_back(IpeSegment(v, v + iInfiniteEdgeLength * dir));
}

void VoronoiIpelet::VoronoiTreatFacet(facetT *facet)
{
  facetT *neighbor, **neighborp;

  if (!facet) return;
  if (qh_skipfacet (facet)) return;
  if (facet == qh_MERGEridge) return;
  if (facet == qh_DUPLICATEridge) return;

  if (iVoronoiSign * facet->normal[2] >= 0.0) return;

  FOREACHneighbor_(facet) {
    if (neighbor != qh_MERGEridge && neighbor != qh_DUPLICATEridge) {
      if (iVoronoiSign * neighbor->normal[2] < 0.0) {
	// make Voronoi edge between the two facets
	AddVoronoiEdge(facet, neighbor);
      } else {
	AddInfiniteEdge(facet, neighbor);
      }
    }
  }
}

// --------------------------------------------------------------------

void VoronoiIpelet::AddDelaunayEdge(int from, int to)
{
  if (from < to)
    iDelaunay.push_back(DelaunayEdge(to, from));
  else
    iDelaunay.push_back(DelaunayEdge(from, to));
}

void VoronoiIpelet::DelaunayTreatFacet(facetT *facet)
{
  setT *vertices;
  vertexT *vertex, **vertexp;

  if (!facet) return;
  if (qh_skipfacet (facet)) return;
  if (facet == qh_MERGEridge) return;
  if (facet == qh_DUPLICATEridge) return;

  if (facet->normal[2] >= 0.0) return;

  vertices= qh_facet3vertex (facet);
  int id, first_id = -1, last_id = -1;
  FOREACHvertex_(vertices) {
    id = qh_pointid(vertex->point);
    if (last_id >= 0) {
      AddDelaunayEdge(last_id, id);
      last_id = id;
    } else {
      last_id = first_id = id;
    }
  }
  AddDelaunayEdge(last_id, first_id);
  qh_settempfree(&vertices);
}

// --------------------------------------------------------------------

#if 0
    // for medial axis: single convex polygon only
    for (IpeObject *ob = ium_input; ob; ob = ob->next) {
      if (ob->type == IPE_LINE) {
	if (sites.size() > 0 || !ob->w.line->closed) {
	  ium_message = "can handle single convex polygon only";
	  ium_end();
	}
	// check whether the polygon is really convex
	pl_polygon pgn(ob->w.line->v);
	if (!pgn.is_convex()) {
	  ium_message = "can handle single convex polygon only";
	  ium_end();
	}
	if (pgn.is_clockwise())
	  pgn.invert_orientation();
	sites = pgn.all_vertices();
      } else if (ob->type != IPE_TEXT) {
	ium_message = "can handle single convex polygon only";
	ium_end();
      }
    }
#endif

void VoronoiIpelet::Run(int function, IpePage *page, IpeletHelper *helper)
{
  if (function == 5) {
    char buf[32];
    std::sprintf(buf, "%g", iInfiniteEdgeLength);
    IpeString el(buf);
    if (helper->GetString("Length of infinite edges (in points):", el))
      iInfiniteEdgeLength = std::strtod(el.CString(), 0);
    return;
  }
  iVoronoiSign = (function == 4) ? -1 : 1;

  iSites.clear();
  CollectVisitor vis(iSites);

  for (IpePage::const_iterator it = page->begin(); it != page->end(); ++it) {
    if (it->Select())
      vis(*it);
  }
  if (iSites.size() < 4) {
    helper->MessageBox("You need to select at least four sites.",
		       "Dismiss", 0, 0);
    return;
  }

  int numpoints;
  coordT *points;
  boolT ismalloc;

  qh_meminit(stderr);
  qh_initqhull_start(stdin, stdout, stderr);

  if (!setjmp(qh errexit)) {
    // the command(s) for qhull

    // strcpy (qh qhull_command, qhull_flags);
    // fprintf(stderr, "Qhull flags = %s\n", qh qhull_command);
    // qh_initflags (qh qhull_command);

    points = readpoints(iSites, function, &numpoints, &ismalloc);
    qh_initqhull_globals (points, numpoints, 3, ismalloc);
    qh_initqhull_mem();
    /* mem.c and set.c are initialized */
    qh_initqhull_buffers();
    qh_initthresholds (qh qhull_command);
    if (qh SCALEinput)
      qh_scaleinput();
    if (qh ROTATErandom >= 0) {
      qh_randommatrix(qh gm_matrix, qh hull_dim, qh gm_row);
      qh_gram_schmidt(qh hull_dim, qh gm_row);
      qh_rotateinput(qh gm_row);
    }
    qh_qhull();
    qh_check_output();
    qh_produce_output();
    if (qh VERIFYoutput && !qh FORCEoutput && !qh STOPpoint && !qh STOPcone)
      qh_check_points();

    // now create segments for Ipe

    facetT *facet, *facetlist = qh facet_list;

    FORALLfacet_(facetlist) {
      if (!function)
	DelaunayTreatFacet(facet);
      else
	VoronoiTreatFacet(facet);
    }

    IpeGroup *group = new IpeGroup;
    IpeAllAttributes attr; // all null attributes

    if (!function) {
      std::sort(iDelaunay.begin(), iDelaunay.end());
      for (uint j = 0; j < iDelaunay.size(); j++) {
	if (!j || iDelaunay[j] != iDelaunay[j-1]) {
 	  IpeVector a(points[3*iDelaunay[j].a], points[3*iDelaunay[j].a + 1]);
 	  IpeVector b(points[3*iDelaunay[j].b], points[3*iDelaunay[j].b + 1]);
 	  IpeSegment seg(a,b);
	  group->push_back(new IpePath(attr, seg));
	}
      }
    } else {
      for (uint i = 0; i < iEdges.size(); ++i)
	group->push_back(new IpePath(attr, iEdges[i]));
    }
    group->SetStroke(helper->Attributes().iStroke);
    group->SetDashStyle(helper->Attributes().iDashStyle);
    group->SetLineWidth(helper->Attributes().iLineWidth);
    page->push_back(IpePgObject(IpePgObject::ESecondary,
				helper->CurrentLayer(), group));
  }
  qh NOerrexit= True;  /* no more setjmp */
  iEdges.clear();
  iSites.clear();
  iDelaunay.clear();
}

// --------------------------------------------------------------------

IPELET_DECLARE Ipelet *NewIpelet()
{
  return new VoronoiIpelet;
}

// --------------------------------------------------------------------