Logo Search packages:      
Sourcecode: k3d version File versions

polyhedron.cpp

Go to the documentation of this file.
// K-3D
// Copyright (c) 1995-2004, Timothy M. Shead
//
// Contact: tshead@k-3d.com
//
// 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 2 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, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

/** \file
            \author Romain Behar (romainbehar@yahoo.com)
*/

#include <k3dsdk/classes.h>
#include <k3dsdk/color.h>
#include <k3dsdk/imaterial.h>
#include <k3dsdk/material.h>
#include <k3dsdk/material_collection.h>
#include <k3dsdk/measurement.h>
#include <k3dsdk/mesh_source.h>
#include <k3dsdk/module.h>
#include <k3dsdk/object.h>
#include <k3dsdk/persistence.h>
#include <k3dsdk/plugins.h>
#include <k3dsdk/property.h>
#include <k3dsdk/transform.h>

#include <cctype>
#include <cmath>
#include <cstdlib>
#include <iterator>
#include <sstream>
#include <vector>

#ifndef MAXLONG
#define MAXLONG 0x7FFFFFFF
#endif

#ifndef DBL_EPSILON
#define DBL_EPSILON 2.2204460492503131e-16
#endif

#define BIG_EPSILON 3e-2

#define Err(x) { std::cerr << "Kaleido error : " << x << std::endl; \
 std::cerr << " line : " << __LINE__ << std::endl; return 0; }


#define AZ M_PI/7 // axis azimuth
#define EL M_PI/17 // axis elevation

#define xcompl(x)(frac(x),(double)frax.n/(frax.n-frax.d))

namespace libk3dmesh
{

namespace kaleido
{

/**********************************************************
* kaleido
*
*     Kaleidoscopic construction of uniform polyhedra
*     Copyright  1991-1997 Dr. Zvi Har'El <rl@math.technion.ac.il>
*
*     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 2 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, write to the Free Software
*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
*     Author:
*           Dr. Zvi Har'El,
*           Technion, Israel Institute of Technology,
*           Haifa 32000, Israel.
*           E-Mail: rl@math.technion.ac.il
**********************************************************/

/** \file
            \author Dr. Zvi Har'El <rl@math.technion.ac.il>
*/

class Kvector
{
public:
      double x,y,z;

      Kvector()
      {
            x = y = z = 0;
      }

      Kvector(const double X, const double Y, const double Z)
      {
            x = X;
            y = Y;
            z = Z;
      }

      Kvector& operator = (const Kvector& v)
      {
            x = v.x;
            y = v.y;
            z = v.z;

            return *this;
      }
};

Kvector operator + (const Kvector& a, const Kvector& b);
Kvector operator - (const Kvector& a, const Kvector& b);
Kvector operator * (const Kvector& a, const double d);
Kvector operator * (const double d, const Kvector& a);
Kvector operator ^ (const Kvector& a, const Kvector& b);

inline Kvector operator + (const Kvector& a, const Kvector& b)
{ return Kvector(a.x + b.x, a.y + b.y, a.z + b.z); }

inline Kvector operator - (const Kvector& a, const Kvector& b)
{ return Kvector(a.x - b.x, a.y - b.y, a.z - b.z); }

inline Kvector operator * (const Kvector& a, const double d)
{ return Kvector(d*a.x, d*a.y, d*a.z); }

inline Kvector operator * (const double d, const Kvector& a)
{ return a*d; }

inline Kvector operator * (const Kvector& a, const Kvector& b)
{ return Kvector(a.x * b.x, a.y * b.y, a.z * b.z); }

inline Kvector operator ^ (const Kvector& a, const Kvector& b) {
 return Kvector(a.y*b.z - a.z*b.y,
            a.z*b.x - a.x*b.z,
            a.x*b.y - a.y*b.x);
}


typedef struct { // See uniformnames.cpp for explanation of the fields
      char *Wythoff, *name, *dual;
      short Coxeter, Wenninger;
} Uniform;

typedef struct {long n,d;} Fraction;

Kvector  rotate( Kvector vertex, Kvector axis, double angle );

double  dot( Kvector a, Kvector b );

void  frac( double x );


extern Uniform uniform[];
extern Fraction frax;
extern int last_uniform;

long numerator(double x)
{
      frac(x);
      return frax.n;
}

long denominator(double x)
{
      frac(x);
      return frax.d;
}

int mod(int i, const int j);


class Polyhedron;

Polyhedron* kaleido( unsigned int nsym );


class Polyhedron
{
public:
      Polyhedron()
      {
            index = -1;
            N = M = V = E = F = D = chi = g = 0;
            K = 2;
            hemi = onesided = 0;
            even = -1;
            p[0] = p[1] = p[2] = p[3] = minr = gon = 0.0;
            polyform = "";
            config = "";
            name = "";
            dual_name = "";
      }

      ~Polyhedron()
      {
      }

      /* NOTE: some of the int's can be replaced by short's, char's,
            or even bit fields, at the expense of readability!!!*/
      int index; /* index to the standard list, the array uniform[] */
      int N; /* number of faces types (atmost 5)*/
      int M; /* vertex valency  (may be big for dihedral polyhedra) */
      int V; /* vertex count */
      int E; /* edge count */
      int F; /* face count */
      int D; /* density */
      int chi; /* Euler characteristic */
      int g; /* order of symmetry group */
      int K; /* symmetry type: D=2, T=3, O=4, I=5 */
      int hemi;/* flag hemi polyhedron */
      int onesided;/* flag onesided polyhedron */
      int even; /* removed face in pqr| */
      std::vector<int> Fi; // face counts by type (array N)
      std::vector<int> rot; // vertex configuration (array M of 0..N-1)
      std::vector<int> snub; // snub triangle configuration (array M of 0..1)
      std::vector<int> firstrot; // temporary for vertex generation (array V)
      std::vector<int> ftype; // face types (array F)
      std::vector< std::vector<int> > incid; // vertex-face incidence (matrix M x V of 0..F-1)
      std::vector< std::vector<int> > adj; // vertex-vertex adjacency (matrix M x V of 0..V-1)
      double p[4]; /* p, q and r; |=0 */
      double minr; /* smallest nonzero inradius */
      double gon; /* basis type for dihedral polyhedra */
      std::vector<double> n_sides; // number of side of a face of each type (array N)
      std::vector<double> m_faces; // number of faces at a vertex of each type (array N)
      std::vector<double> gamma; // fundamental angles in radians (array N)
      std::string polyform; // printable Wythoff symbol
      std::string config; // printable vertex configuration
      std::string name; // name, standard or manufactured
      std::string dual_name; // dual name, standard or manufactured
      std::vector<Kvector> v; // vertex coordinates (array V)
      std::vector<Kvector> f; // face coordinates (array F)

      int UnpackSym(unsigned int nsym);
      int Moebius();
      void Decompose();
      void GuessName();
      int Newton();
      void Exceptions();
      void Count();
      void Configuration();
      int Faces();
      int Vertices();

private:
      int Dihedral(const std::string& newname, const std::string& newdual_name);

};


const std::string PrintN(int n)
{
      std::ostringstream stream;
      stream << n;

      return std::string(stream.str());
}

const std::string Sprintfrac(double x)
{
      frac(x);

      if(!frax.d)
            return "infinity";

      if(frax.d == 1)
            return PrintN(frax.n);

      return PrintN(frax.n) + "/" + PrintN(frax.d);
}

double Same(Kvector a, Kvector b, double epsilon)
{
      return fabs(a.x - b.x) < epsilon
            && fabs(a.y - b.y) < epsilon
            && fabs(a.z - b.z) < epsilon;
}

// Compute the polar reciprocal of the plane containing a, b and c:
//
// If this plane does not contain the origin, return p such that
//    dot(p,a) = dot(p,b) = dot(p,b) = r.
//
// Otherwise, return p such that
//    dot(p,a) = dot(p,b) = dot(p,c) = 0
// and
//    dot(p,p) = 1.
Kvector Pole(double r, Kvector a, Kvector b, Kvector c)
{
      Kvector p = (b - a) ^ (c - a);
      double k = dot(p, a);
      if(fabs(k) < 1e-6)
            return (1 / sqrt(dot(p, p)))*p;

      return (r/k)*p;
}

// Compute the mathematical modulus function
int mod(int i, const int j)
{
      return (i %= j) >= 0 ? i : (j < 0 ? i-j : i+j);
}

Polyhedron* kaleido( unsigned int nsym )
{
      // Allocate a Polyhedron structure P
      Polyhedron* P = new Polyhedron();
      if(!P)
            return 0;

      // Unpack input symbol into P
      if(!P->UnpackSym(nsym))
            return 0;

      // Find Mebius triangle, its density and Euler characteristic
      if(!P->Moebius())
            return 0;

      // Decompose Schwarz triangle
      P->Decompose();

      // Find the names of the polyhedron and its dual
      P->GuessName();

      // Solve Fundamental triangles, optionally printing approximations
      if(!P->Newton())
            return 0;

      // Deal with exceptional polyhedra
      P->Exceptions();

      // Count edges and faces, update density and characteristic if needed
      P->Count();

      // Generate printable vertex configuration
      P->Configuration();

      // Compute coordinates
      if(!P->Vertices()) return 0;
      if(!P->Faces()) return 0;

      return P;
}

int Polyhedron::UnpackSym(unsigned int nsym)
/*
* Unpack input symbol: Wythoff symbol from nsym, an index to uniform[].
* The symbol is a three fractions and a bar in some order.
* We allow no bar only if it results from the input symbol #80.
*/
{
      if(nsym < 1 || nsym > 80)
            Err("Bad number");

      char* sym = uniform[index = nsym - 1].Wythoff;

      int i = 0, bars = 0;
      char c;
      for (;;) {
            while((c = *sym++) && std::isspace(c));
            if(!c)
            {
                  if (i == 4 && (bars || index == last_uniform - 1)) return 1;
                  if (!bars) Err("no bars")
                  Err("not enough fractions")
            }
            if (i == 4) Err("data exceeded")
            if (c == '|'){
                  if (++bars > 1) Err("too many bars")
                  p[i++] = 0;
                  continue;
            }
            if (!std::isdigit(c)) Err("not a digit")
            int n = c - '0';
            while ((c = *sym++) && std::isdigit(c)) n = n * 10 + c - '0';
            if (c && std::isspace (c)) while ((c = *sym++) && std::isspace(c))
                  ;
            if (c != '/') {
                  sym--;
                  if ((p[i++] = n) <= 1) Err("fraction<=1")
                  continue;
            }
            while( (c = *sym++) && std::isspace(c));

            if (!c || !std::isdigit (c)) return 0;
            int d = c - '0';
            while ((c = *sym++) && std::isdigit(c)) d = d * 10 + c - '0';
            if (!d) Err ("zero denominator")
            sym--;
            if ((p[i++] = (double) n / d) <= 1) Err("fraction<=1")
      }
}

int Polyhedron::Moebius()
/*
* Using Wythoff symbol (p|qr, pq|r, pqr| or |pqr), find the Moebius
* triangle (2 3 K) (or (2 2 n)) of the Schwarz triangle (pqr),
* the order g of its symmetry group, its Euler characteristic chi,
* and its covering density D.
* g is the number of copies of (2 3 K) covering the sphere, i.e.,
*
*           g * pi * (1/2 + 1/3 + 1/K - 1) = 4 * pi
*
* D is the number of times g copies of (pqr) cover the sphere, i.e.
*
*           D * 4 * pi = g * pi * (1/p + 1/q + 1/r - 1)
*
* chi is V - E + F, where F = g is the number of triangles, E = 3*g/2 is the
* number of triangle edges, and V = Vp+ Vq+ Vr, with Vp = g/(2*np) being the
* number of vertices with angle pi/p (np is the numerator of p).
*/
{
      int twos = 0, j;
/*
* Arrange Wythoff symbol in a presentable form. In the same time check
* the restrictions on the three fractions: They all have to be greater
* then one, and the numerators 4 or 5 cannot occur together.
* We count the ocurrences of 2 in `two', and save the largest numerator
* in `K', since they reflect on the symmetry group.
*/
      K = 2;
      if(index == last_uniform - 1)
            polyform = "|";

      for(j = 0; j < 4; j++) {
            if(p[j]) {
                  std::string s = Sprintfrac(p[j]);
                  if(j && p[j-1])
                        polyform += " ";
                  else
                        polyform += s;

                  if(p[j] != 2) {
                        int k;
                        if ((k = numerator (p[j])) > K) {
                              if (K == 4) break;
                              K = k;
                        } else if (k < K && k == 4) break;
                  } else twos++;
            } else
                  polyform += "|";
      }

/*
* Find the symmetry group K (where 2, 3, 4, 5 represent the dihedral,
* tetrahedral, octahedral and icosahedral groups, respectively),
* and its order g.
*/
      if (twos >= 2) { // dihedral
            g = 4 * K;
            K = 2;
      } else {
            if (K > 5) Err ("numerator too large")
            g = 24 * K / (6 - K);
      }
/*
* Compute the nominal density D and Euler characteristic chi.
* In few exceptional cases, these values will be modified later.
*/
      if (index != last_uniform - 1) {
            int i;
            D = chi = - g;
            for(j = 0; j < 4; j++) if (p[j]) {
                  chi +=  i = g / numerator (p[j]);
                  D += i * denominator (p[j]);
            }
            chi /= 2;
            D /= 4;
            if(D <= 0) Err ("nonpositive density")
      }
      return 1;
}

void Polyhedron::Decompose()
/*
* Decompose Schwarz triangle into N right triangles and compute
* the vertex count V and the vertex valency M.
* V is computed from the number g of Schwarz triangles in the cover, divided
* by the number of triangles which share a vertex. It is halved for one-sided
* polyhedra, because the kaleidoscopic construction really produces a double
* orientable covering of such polyhedra. All q' q|r are of the "hemi" type,
* i.e. have equatorial {2r} faces, and therefore are (except 3/2 3|3 and the
* dihedra 2 2|r) one-sided. A well known example is 3/2 3|4, the "one-sided
* heptahedron". Also, all p q r| with one even denominator have a crossed
* parallelogram as a vertex figure, and thus are one-sided as well.
*/
{
      int j;

      if(!p[1]) { // p|q r
            N = 2;
            M = 2 * numerator (p[0]);
            V = g / M;
            n_sides.resize(N);
            m_faces.resize(N);

            for(j = 0; j < 2; j++) {
                  n_sides[j] = p[j+2];
                  m_faces[j] = p[0];
            }
            for (j = M / 2; j--;) {
                  rot.push_back(0);
                  rot.push_back(1);
            }
      } else if(!p[2]) { // p q|r
            N = 3;
            M = 4;
            V = g / 2;
            n_sides.resize(N);
            m_faces.resize(N);

            n_sides[0] = 2 * p[3];
            m_faces[0] = 2;
            for (j = 1; j < 3; j++) {
                  n_sides[j] = p[j-1];
                  m_faces[j] = 1;
                  rot.push_back(0);
                  rot.push_back(j);
            }
            if (fabs (p[0] - xcompl (p[1])) < DBL_EPSILON) {/* p = q' */
                        /* p[0]==xcompl(p[1]) should work. */
                        /* However, MSDOS yeilds a 7e-17 difference! */
                        /* Reported by Jim Buddenhagen <jb1556@daditz.sbc.com> */
                  hemi = 1;
                  D = 0;
                  if (p[0] != 2 && !(p[3] == 3 && (p[0] == 3 || p[1] == 3))) {
                        onesided = 1;
                        V /= 2;
                        chi /= 2;
                  }
            }
      } else if (!p[3]) { // p q r|
            M = N = 3;
            V = g;
            n_sides.resize(N);
            m_faces.resize(N);

            for (j = 0; j < 3; j++) {
                  if (!(denominator (p[j]) % 2)) {
                        /* what happens if there is more than one even denominator???*/
                        if (p[(j+1)%3] != p[(j+2)%3]) { /*needs postprocessing*/
                              even = j;/*memorize the removed face*/
                              chi -= g / numerator (p[j]) / 2;
                              onesided = 1;
                              D = 0;
                        } else {/*for p = q we get a double 2 2r|p*/
                        /*noted by Roman Maeder <maeder@inf.ethz.ch> for 4 4 3/2|*/
                        /*Euler characteristic is still wrong*/
                              D /= 2;
                        }
                        V /= 2;
                  }
                  n_sides[j] = 2 * p[j];
                  m_faces[j] = 1;
                  rot.push_back(j);
            }
      } else { /* |p q r - snub polyhedron*/
            N = 4;
            M = 6;
            V = g / 2;/* Only "white" triangles carry a vertex */
            n_sides.resize(N);
            m_faces.resize(N);
            snub.resize(M);

            std::vector<int>::iterator t = snub.begin();
                  m_faces[0] = n_sides[0] = 3;
                  for (j = 1; j < 4; j++) {
                        n_sides[j] = p[j];
                        m_faces[j] = 1;
                        rot.push_back(0);
                        rot.push_back(j);
                        *t++ = 1;
                        *t++ = 0;
                  }
            }

      // Sort the fundamental triangles (using bubble sort) according to
      // decreasing n[i], while pushing the trivial triangles (n[i] = 2)
      // to the end.
      int J = N - 1;
      while(J) {
            int last;
            last = J;
            J = 0;
            for (j = 0; j < last; j++) {
                  if ((n_sides[j] < n_sides[j+1] || n_sides[j] == 2) && n_sides[j+1] != 2) {
                        int i;
                        double temp;
                        temp = n_sides[j];
                        n_sides[j] = n_sides[j+1];
                        n_sides[j+1] = temp;
                        temp = m_faces[j];
                        m_faces[j] = m_faces[j+1];
                        m_faces[j+1] = temp;
                        for (i = 0; i < M; i++) {
                              if(rot[i] == j) rot[i] = j+1;
                              else if(rot[i] == j+1) rot[i] = j;
                        }
                        if (even != -1) {
                              if (even == j) even = j+1;
                              else if (even == j+1) even = j;
                        }
                        J = j;
                  }
            }
      }

      // Get rid of repeated triangles
      for (J = 0; J < N && n_sides[J] != 2;J++) {
            int k, i;
            for (j = J+1; j < N && n_sides[j]==n_sides[J]; j++) m_faces[J] += m_faces[j];
            k = j - J - 1;
            if (k) {
                  for (i = j; i < N; i++) {
                        n_sides[i - k] = n_sides[i];
                        m_faces[i - k] = m_faces[i];
                  }
                  N -= k;
                  for (i = 0; i < M; i++) {
                        if(rot[i] >= j) rot[i] -= k;
                        else if(rot[i] > J) rot[i] = J;
                  }
                  if (even >= j) even -= k;
            }
      }

      // Get rid of trivial triangles
      if (!J) J = 1; /*hosohedron*/
      if (J < N) {
            int i;
            N = J;
            for (i = 0; i < M; i++) {
                  if(rot[i] >= N) {
                        for (j = i + 1; j < M; j++) {
                              rot[j-1] = rot[j];
                              if(snub.size()) snub[j-1] = snub[j];
                        }
                        M--;
                  }
            }
      }

      // Truncate arrays
      n_sides.resize(N);
      m_faces.resize(N);
      if(snub.size()) snub.resize(M);
}

void Polyhedron::GuessName()
/*
* Get the polyhedron name, using standard list or guesswork.
* Ideally, we should try to locate the Wythoff symbol in the standard list
* (unless, of course, it is dihedral), after doing few normalizations,
* such as sorting angles and splitting isoceles triangles.
*/
{
      if(index != -1)
      {
            // tabulated
            name = uniform[index].name;
            dual_name = uniform[index].dual;
            return;
      }

      if(K == 2)
      {
            // dihedral nontabulated
            if(!p[0])
            {
                  if(N == 1)
                  {
                        name = "octahedron";
                        dual_name = "cube";
                        return;
                  }

                  gon = n_sides[0] == 3 ? n_sides[1] : n_sides[0];
                  if(gon >= 2)
                        Dihedral("antiprism", "deltohedron");
                  else
                        Dihedral("crossed antiprism", "concave deltohedron");

                  return;
            }

            if(!p[3] || !p[2] && p[3] == 2)
            {
                  if(N == 1) {
                        name = "cube";
                        dual_name = "octahedron";
                        return;
                  }

                  gon = n_sides[0] == 4 ? n_sides[1] : n_sides[0];

                  Dihedral("prism", "dipyramid");
                  return;
            }

            if(!p[1] && p[0] != 2)
            {
                  gon = m_faces[0];
                  Dihedral("hosohedron", "dihedron");
                  return;
            }

            gon = n_sides[0];
            Dihedral("dihedron", "hosohedron");
            return;
      }

      // Other nontabulated
      static std::string pre[] = { "tetr", "oct", "icos" };
      name = pre[K - 3] + "ahedral ";
      if(onesided) name += "one-sided ";
      else if(D == 1) name += "convex ";
      else name += "nonconvex ";

      dual_name = name;
      name += "isogonal polyhedron";
      dual_name += "isohedral polyhedron";
}

int Polyhedron::Dihedral(const std::string& newname, const std::string& newdual_name )
{
      std::string s = Sprintfrac(gon < 2 ? xcompl(gon) : gon);

      name = s + "-gonal " + newname;
      dual_name = s + "-gonal " + newdual_name;

      return 1;
}

int Polyhedron::Newton()
// Solve the fundamental right spherical triangles
{
      // First, we find initial approximations
      int j;

      gamma.resize(N);
      if(N == 1)
      {
            gamma[0] = M_PI / m_faces[0];
            return 1;
      }

      for(j = 0; j < N; j++)
            gamma[j] = M_PI / 2 - M_PI / n_sides[j];

      // Next, iteratively find closer approximations for gamma[0]
      // and compute other gamma[j]'s from Napier's equations.
      for(;;)
      {
            double delta = M_PI;
            for(j = 0; j < N; j++)
                  delta -= m_faces[j] * gamma[j];

            if(fabs(delta) < 11 * DBL_EPSILON)
                  return 1;

            // On a RS/6000, fabs(delta)/DBL_EPSILON may occilate between 8 and 10.
            // Reported by David W. Sanderson <dws@ssec.wisc.edu>
            double sigma = 0;
            for(j = 0; j < N; j++)
                  sigma += m_faces[j] * tan(gamma[j]);

            gamma[0] += delta * tan (gamma[0]) / sigma;
            if(gamma[0] < 0 || gamma[0] > M_PI)
                  Err("gamma out of bounds")

            double cosa = cos(M_PI / n_sides[0]) / sin (gamma[0]);
            for(j = 1; j < N; j++)
                  gamma[j] = asin(cos (M_PI / n_sides[j]) / cosa);
      }
}

void Polyhedron::Exceptions()
{
      // Postprocess pqr| where r has an even denominator (cf. Coxeter &al. Sec.9).
      // Remove the {2r} and add a retrograde {2p} and retrograde {2q}.
      int j;
      if(even != -1)
      {
            M = N = 4;
            n_sides.resize(N);
            m_faces.resize(N);
            gamma.resize(N);
            rot.resize(M);

            for(j = even + 1; j < 3; j++)
            {
                  n_sides[j-1] = n_sides[j];
                  gamma[j-1] = gamma[j];
            }

            n_sides[2] = xcompl (n_sides[1]);
            gamma[2] = - gamma[1];
            n_sides[3] = xcompl (n_sides[0]);
            m_faces[3] = 1;
            gamma[3] = - gamma[0];
            rot[0] = 0;
            rot[1] = 1;
            rot[2] = 3;
            rot[3] = 2;
      }

      // Postprocess the last polyhedron |3/2 5/3 3 5/2 by taking a |5/3 3 5/2,
      // replacing the three snub triangles by four equatorial squares and adding
      // the missing {3/2} (retrograde triangle, cf. Coxeter &al. Sec. 11).
      if(index == last_uniform - 1)
      {
            N = 5;
            M = 8;
            n_sides.resize(N);
            m_faces.resize(N);
            gamma.resize(N);
            rot.resize(M);
            snub.resize(M);

            hemi = 1;
            D = 0;
            for(j = 3; j; j--)
            {
                  m_faces[j] = 1;
                  n_sides[j] = n_sides[j-1];
                  gamma[j] = gamma[j-1];
            }

            m_faces[0] = n_sides[0] = 4;
            gamma[0] = M_PI / 2;
            m_faces[4] = 1;
            n_sides[4] = xcompl (n_sides[1]);
            gamma[4] = - gamma[1];
            for(j = 1; j < 6; j += 2) rot[j]++;
            rot[6] = 0;
            rot[7] = 4;
            snub[6] = 1;
            snub[7] = 0;
      }
}

void Polyhedron::Count()
/*
* Compute edge and face counts, and update D and chi.
* Update D in the few cases the density of the polyhedron is meaningful but
* different than the density of the corresponding Schwarz triangle
* (cf. Coxeter &al., p. 418 and p. 425).
* In these cases, spherical faces of one type are concave (bigger than a
* hemisphere), and the actual density is the number of these faces less the
* computed density.
* Note that if j != 0, the assignment gamma[j] = asin(...)
* implies gamma[j] cannot be obtuse.
* Also, compute chi for the only non-Wythoffian polyhedron.
*/
{
      Fi.resize(N);
      for(int j = 0; j < N; j++)
      {
            int temp = V * numerator(m_faces[j]);
            E += temp;
            F += Fi[j] = temp / numerator(n_sides[j]);
      }

      E /= 2;

      if(D && gamma[0] > M_PI / 2)
            D = Fi[0] - D;

      if(index == last_uniform - 1)
            chi = V - E + F;
}

// Generate a printable vertex configuration symbol
void Polyhedron::Configuration()
{
      int j;
      for(j = 0; j < M; j++)
      {
            std::string s = Sprintfrac(n_sides[rot[j]]);

            if(!j)
                  config = "(";
            else
                  config += ".";

            config += s;
      }

      config += ")";

      if((j = denominator(m_faces[0])) != 1)
            config += "/" + PrintN(j);
}

int Polyhedron::Vertices()
/*
* Compute polyhedron vertices and vertex adjecency lists.
* The vertices adjacent to v[i] are v[adj[0][i], v[adj[1][i], ... v[adj[M-1][i],
* ordered counterclockwise.
* The algorith is a BFS on the vertices, in such a way that the vetices
* adjacent to a givem vertex are obtained from its BFS parent by a cyclic
* sequence of rotations. firstrot[i] points to the first  rotaion in
* the sequence when applied to v[i]. Note that for non-snub polyhedra,
* the rotations at a child are opposite in sense when compared to the
* rotations at the parent. Thus, we fill adj[*][i] from the end to signify
* clockwise rotations. The firstrot[] array is not needed for display
* thus it is freed after being used for face computations below.
*/
{
      int newV = 2;

      v.resize(V);

      adj.resize(M);
      for(int row = 0; row < M; row++)
            (adj[row]).resize(V, 0);

      firstrot.resize(V); /* temporary , put in Polyhedron
                  structure so that may be freed on error */
      double cosa = cos (M_PI / n_sides[0]) / sin (gamma[0]);
      v[0].x = 0;
      v[0].y = 0;
      v[0].z = 1;
      firstrot[0] = 0;
      adj[0][0] = 1;
      v[1].x = 2 * cosa * sqrt (1 - cosa * cosa);
      v[1].y = 0;
      v[1].z = 2 * cosa * cosa - 1;
      if(!snub.size()) {
            firstrot[1] = 0;
            adj[0][1] = -1;/*start the other side*/
            adj[M-1][1] = 0;
      } else {
            firstrot[1] = snub[M-1] ? 0 : M-1;
            adj[0][1] = 0;
      }
      for(int i = 0; i < newV; i++) {
            int j, k;
            int last, one, start, limit;
            if (adj[0][i] == -1) {
                  one = -1; start = M-2; limit = -1;
            } else {
                  one = 1; start = 1; limit = M;
            }
            k = firstrot[i];
            for (j = start; j != limit; j += one) {
                  Kvector temp;
                  int J;
                  temp = rotate (v[adj[j-one][i]], v[i],
                        one * 2 * gamma[rot[k]]);
                  for (J=0; J<newV && !Same(v[J], temp, BIG_EPSILON); J++) /*noop*/;
                  adj[j][i] = J;
                  last = k;
                  if (++k == M) k = 0;
                  if (J == newV) { /*new vertex*/
                        if (newV == V) Err ("too many vertices")
                        v[newV++] = temp;
                        if(!snub.size()) {
                              firstrot[J] = k;
                              if (one > 0) {
                                    adj[0][J] = -1;
                                    adj[M-1][J] = i;
                              } else {
                                    adj[0][J] = i;
                              }
                        } else {
                              firstrot[J] = !snub[last] ? last :
                                    !snub[k] ? (k+1)%M : k ;
                              adj[0][J] = i;
                        }
                  }
            }
      }
      return 1;
}

int Polyhedron::Faces()
/*
* Compute polyhedron faces (dual vertices) and incidence matrices.
* For orientable polyhedra, we can distinguish between the two faces meeting
* at a given directed edge and identify the face on the left and the face on
* the right, as seen from the outside.
* For one-sided polyhedra, the vertex figure is a papillon (in Coxeter &al.
* terminology, a crossed parallelogram) and the two faces meeting at an edge
* can be identified as the side face (n[1] or n[2]) and
* the diagonal face (n[0] or n[3]).
*/
{
      f.resize(F);
      ftype.resize(F);

      // Init incid matrix with -1
      incid.resize(M);
      for(int row = 0; row < M; row++)
            (incid[row]).resize(V, -1);

      minr = 1 / fabs (tan (M_PI / n_sides[hemi]) * tan (gamma[hemi]));

      int newF = 0;
      for(int i = 0; i < V; i++) {
            for(int j = 0; j < M; j++) {
                  int pap = 0; // papillon edge type

                  if (incid[j][i] != -1) continue;
                  incid[j][i] = newF;
                  if (newF == F) Err ("too many faces")
                  f[newF] = Pole(minr, v[i], v[adj[j][i]],
                        v[adj[mod(j + 1, M)][i]]);
                  ftype[newF] = rot [mod(
                        firstrot[i] + (adj[0][i] < adj[M - 1][i] ?
                        j : -j - 2), M)];
                  if (onesided) pap = (firstrot[i] + j) % 2;
                  int i0 = i;
                  int J = j;
                  for (;;) {
                        int k;
                        k = i0;
                        if ((i0 = adj[J][k]) == i) break;
                        for (J = 0; J < M && adj[J][i0] != k; J++) /*noop*/;
                        if (J == M) Err ("too many faces")
                        if (onesided && (J + firstrot[i0]) % 2 == pap) {
                              incid [J][i0] = newF;
                              if (++J >= M) J = 0;
                        } else {
                              if (--J < 0) J = M - 1;
                              incid [J][i0] = newF;
                        }
                  }
            newF++;
            }
      }

      rot.clear();
      snub.clear();

      return 1;
}

Fraction frax;

void  frac( double x )
/*
* Find the numerator and the denominator using the Euclidean algorithm.
*/
{
      static Fraction zero = {0,1}, inf = {1,0};
      Fraction r0, r;
      long fl;
      double s = x;
      r = zero;
      frax = inf;
      for (;;) {
            if (fabs(s) > (double) MAXLONG) return;
            fl = (long) floor (s);
            r0 = r;
            r = frax;
            frax.n = frax.n * fl + r0.n;
            frax.d = frax.d * fl + r0.d;
            if (x == (double)frax.n/(double)frax.d) return;
            s = 1 / (s - fl);
      }
}

double  dot( Kvector a, Kvector b )
{
      return a.x * b.x + a.y * b.y + a.z * b.z;
}

Kvector  rotate( Kvector vertex, Kvector axis, double angle )
{
      Kvector p = dot(axis, vertex) * axis;
      return p + cos(angle)*(vertex - p) +
            sin(angle)*(axis ^ vertex);
}

/****************************************************************************
*     List of Uniform Polyhedra and Their Kaleidoscopic Formulae
*     ==========================================================
*
*     Each entry contains the following items:
*
*     1)    Wythoff symbol.
*     2)    Polyhedron name.
*     3)    Dual name.
*     4)    Coxeter &al. reference figure.
*     5)    Wenninger reference figure.
*
*     Notes:
*
*     (1)   Cundy&Roulette's trapezohedron has been renamed to
*           deltohedron, as its faces are deltoids, not trapezoids.
*     (2)   The names of the non-dihedral polyhedra are those
*           which appear in Wenninger (1984). Some of them are
*           slightly modified versions of those in Wenninger (1971).
*
*     References:
*
*     Coxeter, H.S.M., Longuet-Higgins, M.S. & Miller, J.C.P.,
*           Uniform polyhedra, Phil. Trans. Royal Soc. London, Ser. A,
*           246 (1953), 401-409.
*     Cundy, H.M. & Rollett, A.P.,
*           "Mathematical Models", 3rd Ed., Tarquin, 1981.
*     Har'El, Z.
*           Unifom solution for uniform polyhedra, Geometriae Dedicata,
*           47 (1993), 57-110.
*     Wenninger, M.J.,
*           "Polyhedron Models", Cambridge University Press, 1971.
*           "Dual Models", Cambridge University Press, 1984.
*
****************************************************************************/

Uniform uniform[] = {

/****************************************************************************
*           Dihedral Schwarz Triangles (D5 only)
****************************************************************************/

/* (2 2 5) (D1/5) */

/*1*/ {"2 5|2","pentagonal prism",
                  "pentagonal dipyramid",0,0},
/*2*/ {"|2 2 5","pentagonal antiprism",
                  "pentagonal deltohedron",0,0},

/* (2 2 5/2) (D2/5) */

/*3*/ {"2 5/2|2","pentagrammic prism",
                  "pentagrammic dipyramid",0,0},
/*4*/ {"|2 2 5/2","pentagrammic antiprism",
                  "pentagrammic deltohedron",0,0},

/* (5/3 2 2) (D3/5) */

/*5*/ {"|2 2 5/3","pentagrammic crossed antiprism",
                  "pentagrammic concave deltohedron",0,0},

/****************************************************************************
*           Tetrahedral Schwarz Triangles
****************************************************************************/

/* (2 3 3) (T1) */

/*6*/ {"3|2 3","tetrahedron",
                  "tetrahedron",15,1},
/*7*/ {"2 3|3","truncated tetrahedron",
                  "triakistetrahedron",16,6},

/* (3/2 3 3) (T2) */

/*8*/ {"3/2 3|3","octahemioctahedron",
                  "octahemioctacron",37,68},

/* (3/2 2 3) (T3) */

/*9*/ {"3/2 3|2","tetrahemihexahedron",
                  "tetrahemihexacron",36,67},

/****************************************************************************
*           Octahedral Schwarz Triangles
****************************************************************************/

/* (2 3 4) (O1) */

/*10*/      {"4|2 3","octahedron",
                  "cube",17,2},
/*11*/      {"3|2 4","cube",
                  "octahedron",18,3},
/*12*/      {"2|3 4","cuboctahedron",
                  "rhombic dodecahedron",19,11},
/*13*/      {"2 4|3","truncated octahedron",
                  "tetrakishexahedron",20,7},
/*14*/      {"2 3|4","truncated cube",
                  "triakisoctahedron",21,8},
/*15*/      {"3 4|2","rhombicuboctahedron",
                  "deltoidal icositetrahedron",22,13},
/*16*/      {"2 3 4|","truncated cuboctahedron",
                  "disdyakisdodecahedron",23,15},
/*17*/      {"|2 3 4","snub cube",
                  "pentagonal icositetrahedron",24,17},

/* (3/2 4 4) (O2b) */

/*18*/      {"3/2 4|4","small cubicuboctahedron",
                  "small hexacronic icositetrahedron",38,69},

/* (4/3 3 4) (O4) */

/*19*/      {"3 4|4/3","great cubicuboctahedron",
                  "great hexacronic icositetrahedron",50,77},
/*20*/      {"4/3 4|3","cubohemioctahedron",
                  "hexahemioctacron",51,78},
/*21*/      {"4/3 3 4|","cubitruncated cuboctahedron",
                  "tetradyakishexahedron",52,79},

/* (3/2 2 4) (O5) */

/*22*/      {"3/2 4|2","great rhombicuboctahedron",
                  "great deltoidal icositetrahedron",59,85},
/*23*/      {"3/2 2 4|","small rhombihexahedron",
                  "small rhombihexacron",60,86},

/* (4/3 2 3) (O7) */

/*24*/      {"2 3|4/3","stellated truncated hexahedron",
                  "great triakisoctahedron",66,92},
/*25*/      {"4/3 2 3|","great truncated cuboctahedron",
                  "great disdyakisdodecahedron",67,93},

/* (4/3 3/2 2) (O11) */

/*26*/      {"4/3 3/2 2|","great rhombihexahedron",
                  "great rhombihexacron",82,103},

/****************************************************************************
*           Icosahedral Schwarz Triangles
****************************************************************************/

/* (2 3 5) (I1) */

/*27*/      {"5|2 3","icosahedron",
                  "dodecahedron",25,4},
/*28*/      {"3|2 5","dodecahedron",
                  "icosahedron",26,5},
/*29*/      {"2|3 5","icosidodecahedron",
                  "rhombic triacontahedron",28,12},
/*30*/      {"2 5|3","truncated icosahedron",
                  "pentakisdodecahedron",27,9},
/*31*/      {"2 3|5","truncated dodecahedron",
                  "triakisicosahedron",29,10},
/*32*/      {"3 5|2","rhombicosidodecahedron",
                  "deltoidal hexecontahedron",30,14},
/*33*/      {"2 3 5|","truncated icosidodechedon",
                  "disdyakistriacontahedron",31,16},
/*34*/      {"|2 3 5","snub dodecahedron",
                  "pentagonal hexecontahedron",32,18},

/* (5/2 3 3) (I2a) */

/*35*/      {"3|5/2 3","small ditrigonal icosidodecahedron",
                  "small triambic icosahedron",39,70},
/*36*/      {"5/2 3|3","small icosicosidodecahedron",
                  "small icosacronic hexecontahedron",40,71},
/*37*/      {"|5/2 3 3","small snub icosicosidodecahedron",
                  "small hexagonal hexecontahedron",41,110},

/* (3/2 5 5) (I2b) */

/*38*/      {"3/2 5|5","small dodecicosidodecahedron",
                  "small dodecacronic hexecontahedron",42,72},

/* (2 5/2 5) (I3) */

/*39*/      {"5|2 5/2","small stellated dodecahedron",
                  "great dodecahedron",43,20},
/*40*/      {"5/2|2 5","great dodecahedron",
                  "small stellated dodecahedron",44,21},
/*41*/      {"2|5/2 5","great dodecadodecahedron",
                  "medial rhombic triacontahedron",45,73},
/*42*/      {"2 5/2|5","truncated great dodecahedron",
                  "small stellapentakisdodecahedron",47,75},
/*43*/      {"5/2 5|2","rhombidodecadodecahedron",
                  "medial deltoidal hexecontahedron",48,76},
/*44*/      {"2 5/2 5|","small rhombidodecahedron",
                  "small rhombidodecacron",46,74},
/*45*/      {"|2 5/2 5","snub dodecadodecahedron",
                  "medial pentagonal hexecontahedron",49,111},

/* (5/3 3 5) (I4) */

/*46*/      {"3|5/3 5","ditrigonal dodecadodecahedron",
                  "medial triambic icosahedron",53,80},
/*47*/      {"3 5|5/3","great ditrigonal dodecicosidodecahedron",
                  "great ditrigonal dodecacronic hexecontahedron",54,81},
/*48*/      {"5/3 3|5","small ditrigonal dodecicosidodecahedron",
                  "small ditrigonal dodecacronic hexecontahedron",55,82},
/*49*/      {"5/3 5|3","icosidodecadodecahedron",
                  "medial icosacronic hexecontahedron",56,83},
/*50*/      {"5/3 3 5|","icositruncated dodecadodecahedron",
                  "tridyakisicosahedron",57,84},
/*51*/      {"|5/3 3 5","snub icosidodecadodecahedron",
                  "medial hexagonal hexecontahedron",58,112},

/* (3/2 3 5) (I6b) */

/*52*/      {"3/2|3 5","great ditrigonal icosidodecahedron",
                  "great triambic icosahedron",61,87},
/*53*/      {"3/2 5|3","great icosicosidodecahedron",
                  "great icosacronic hexecontahedron",62,88},
/*54*/      {"3/2 3|5","small icosihemidodecahedron",
                  "small icosihemidodecacron",63,89},
/*55*/      {"3/2 3 5|","small dodecicosahedron",
                  "small dodecicosacron",64,90},

/* (5/4 5 5) (I6c) */

/*56*/      {"5/4 5|5","small dodecahemidodecahedron",
                  "small dodecahemidodecacron",65,91},

/* (2 5/2 3) (I7) */

/*57*/      {"3|2 5/2","great stellated dodecahedron",
                  "great icosahedron",68,22},
/*58*/      {"5/2|2 3","great icosahedron",
                  "great stellated dodecahedron",69,41},
/*59*/      {"2|5/2 3","great icosidodecahedron",
                  "great rhombic triacontahedron",70,94},
/*60*/      {"2 5/2|3","great truncated icosahedron",
                  "great stellapentakisdodecahedron",71,95},
/*61*/      {"2 5/2 3|","rhombicosahedron",
                  "rhombicosacron",72,96},
/*62*/      {"|2 5/2 3","great snub icosidodecahedron",
                  "great pentagonal hexecontahedron",73,113},

/* (5/3 2 5) (I9) */

/*63*/      {"2 5|5/3","small stellated truncated dodecahedron",
                  "great pentakisdodekahedron",74,97},
/*64*/      {"5/3 2 5|","truncated dodecadodecahedron",
                  "medial disdyakistriacontahedron",75,98},
/*65*/      {"|5/3 2 5","inverted snub dodecadodecahedron",
                  "medial inverted pentagonal hexecontahedron",76,114},

/* (5/3 5/2 3) (I10a) */

/*66*/      {"5/2 3|5/3","great dodecicosidodecahedron",
                  "great dodecacronic hexecontahedron",77,99},
/*67*/      {"5/3 5/2|3","small dodecahemicosahedron",
                  "small dodecahemicosacron",78,100},
/*68*/      {"5/3 5/2 3|","great dodecicosahedron",
                  "great dodecicosacron",79,101},
/*69*/      {"|5/3 5/2 3","great snub dodecicosidodecahedron",
                  "great hexagonal hexecontahedron",80,115},

/* (5/4 3 5) (I10b) */

/*70*/      {"5/4 5|3","great dodecahemicosahedron",
                  "great dodecahemicosacron",81,102},

/* (5/3 2 3) (I13) */

/*71*/      {"2 3|5/3","great stellated truncated dodecahedron",
                  "great triakisicosahedron",83,104},
/*72*/      {"5/3 3|2","great rhombicosidodecahedron",
                  "great deltoidal hexecontahedron",84,105},
/*73*/      {"5/3 2 3|","great truncated icosidodecahedron",
                  "great disdyakistriacontahedron",87,108},
/*74*/      {"|5/3 2 3","great inverted snub icosidodecahedron",
                  "great inverted pentagonal hexecontahedron",88,116},

/* (5/3 5/3 5/2) (I18a) */

/*75*/      {"5/3 5/2|5/3","great dodecahemidodecahedron",
                  "great dodecahemidodecacron",86,107},

/* (3/2 5/3 3) (I18b) */

/*76*/      {"3/2 3|5/3","great icosihemidodecahedron",
                  "great icosihemidodecacron",85,106},

/* (3/2 3/2 5/3) (I22) */

/*77*/      {"|3/2 3/2 5/2","small retrosnub icosicosidodecahedron",
                  "small hexagrammic hexecontahedron",91,118},

/* (3/2 5/3 2) (I23) */

/*78*/      {"3/2 5/3 2|","great rhombidodecahedron",
                  "great rhombidodecacron",89,109},
/*79*/      {"|3/2 5/3 2","great retrosnub icosidodecahedron",
                  "great pentagrammic hexecontahedron",90,117},

/****************************************************************************
*           Last But Not Least
****************************************************************************/

/*80*/      {"3/2 5/3 3 5/2","great dirhombicosidodecahedron",
                  "great dirhombicosidodecacron",92,119}
};

int last_uniform = sizeof (uniform) / sizeof (uniform[0]);

} // namespace kaleido

namespace detail
{
      k3d::color colors[13] = {
            k3d::color(0, 0, 0), k3d::color(0, 0, 0), k3d::color(0, 0, 0),
            k3d::color(1, 0, 0), // triangles
            k3d::color(0, 1, 0), // quads
            k3d::color(0, 0, 1), // pentagons
            k3d::color(1, 1, 0), // hexagons
            k3d::color(0, 0, 0),
            k3d::color(1, 0, 1), // octagons
            k3d::color(0, 0, 0), k3d::color(0, 0, 0), k3d::color(0, 0, 0),
            k3d::color(0, 1, 1) // decagons
      };

      std::string names[13] = {
            "", "", "",
            "Triangles material",
            "Quadrilaterals material",
            "Pentagons material",
            "Hexagons material",
            "",
            "Octagons material",
            "", "", "",
            "Decagons material"
      };

} // namespace detail

/////////////////////////////////////////////////////////////////////////////
// polyhedron_implementation

class polyhedron_implementation :
      public k3d::material_collection<k3d::mesh_source<k3d::persistent<k3d::object> > >
{
      typedef k3d::material_collection<k3d::mesh_source<k3d::persistent<k3d::object> > > base;

public:
      polyhedron_implementation(k3d::idocument& Document) :
            base(Document),
            m_document(Document),
            m_number(k3d::init_name("number") + k3d::init_description("Polyhedron Number [number]") + k3d::init_value(1) + k3d::init_constraint(k3d::constraint::minimum(1UL, k3d::constraint::maximum(80UL))) + k3d::init_document(Document) + k3d::init_precision(0) + k3d::init_step_increment(1) + k3d::init_units(typeid(k3d::measurement::scalar))),
            m_dual(k3d::init_name("dual") + k3d::init_description("") + k3d::init_value(false) + k3d::init_document(Document)),
            m_size(k3d::init_name("size") + k3d::init_description("Size [number]") + k3d::init_document(Document) + k3d::init_value(5.0) + k3d::init_precision(2) + k3d::init_step_increment(0.1) + k3d::init_units(typeid(k3d::measurement::scalar)))//,
            //m_previous_materials(k3d::init_name("previousmaterials") + k3d::init_description("Use previous materials [boolean]") + k3d::init_value(true) + k3d::init_document(Document))
      {
            enable_serialization(k3d::persistence::proxy(m_number));
            enable_serialization(k3d::persistence::proxy(m_dual));
            enable_serialization(k3d::persistence::proxy(m_size));
            //enable_serialization(k3d::persistence::proxy(m_previous_materials));

            register_property(m_number);
            register_property(m_dual);
            register_property(m_size);
            //register_property(m_previous_materials);

            m_material.changed_signal().connect(SigC::slot(*this, &polyhedron_implementation::on_reset_geometry));

            m_number.changed_signal().connect(SigC::slot(*this, &polyhedron_implementation::on_reset_geometry));
            m_dual.changed_signal().connect(SigC::slot(*this, &polyhedron_implementation::on_reset_geometry));
            m_size.changed_signal().connect(SigC::slot(*this, &polyhedron_implementation::on_reset_geometry));
            //m_previous_materials.changed_signal().connect(SigC::slot(*this, &polyhedron_implementation::on_reset_geometry));

            m_output_mesh.need_data_signal().connect(SigC::slot(*this, &polyhedron_implementation::on_create_geometry));
      }

      void on_reset_geometry()
      {
            m_output_mesh.reset();
      }

      k3d::mesh* on_create_geometry()
      {
            std::auto_ptr<k3d::mesh> mesh(new k3d::mesh());

            const unsigned long number = m_number.property_value();
            const bool dual = m_dual.property_value();
            const double size= m_size.property_value();

            const std::string name = dual ? kaleido::uniform[number - 1].dual : kaleido::uniform[number - 1].name;

            LoadPolyhedron(number, dual);

            // Add geometry to the mesh ...
            for(unsigned long i = 0; i != m_PolyhedronPoints.size(); ++i)
                  mesh->points.push_back(new k3d::point(m_PolyhedronPoints[i] * m_PolyhedronScaling * size));

            mesh->polyhedra.push_back(new k3d::polyhedron());
            k3d::polyhedron& polyhedron = *mesh->polyhedra.back();
            polyhedron.material = m_material.interface();

            unsigned long polygon_index = 0;
            for(unsigned long i = 0; i != m_PolyhedronPaths.size(); ++i)
                  {
                        if(!m_PolyhedronPaths[i].size())
                              continue;

                        k3d::polyhedron::edges_t edges;
                        for(unsigned long j = 0; j != m_PolyhedronPaths[i].size(); ++j)
                              edges.push_back(new k3d::split_edge(mesh->points[m_PolyhedronPaths[i][j]]));

                        k3d::loop_edges(edges.begin(), edges.end());

                        // Set material ...
                        unsigned long color_index = *m_polygon_colors.begin();
                        // Trick to handle special case #39 and #57
                        if(polygon_index < m_polygon_colors.size())
                              color_index = m_polygon_colors[polygon_index];

                        polyhedron.edges.insert(polyhedron.edges.end(), edges.begin(), edges.end());

                        k3d::face* new_face = new k3d::face(edges.front());
                        return_val_if_fail(new_face, 0);
                        new_face->uniform_data["Cs"] = detail::colors[color_index];
                        polyhedron.faces.push_back(new_face);

                        polygon_index++;
                  }

            return_val_if_fail(is_valid(polyhedron), 0);
            //assert_warning(is_solid(polyhedron));

            return mesh.release();
      }

      k3d::iplugin_factory& factory()
      {
            return get_factory();
      }

      static k3d::iplugin_factory& get_factory()
      {
            static k3d::plugin_factory<k3d::document_plugin<polyhedron_implementation>, k3d::interface_list<k3d::imesh_source > > factory(
                  k3d::uuid(0x916c98dd, 0x74f24362, 0x98364c94, 0x2a9f9ae7),
                  "Polyhedron",
                  "Generates a uniform polyhedron",
                  "Objects",
                  k3d::iplugin_factory::STABLE);

            return factory;
      }

private:
      k3d::idocument& m_document;
      k3d_measurement_property(unsigned long, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::with_constraint) m_number;
      k3d_data_property(bool, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::no_constraint) m_dual;

      k3d_measurement_property(double, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::no_constraint) m_size;

      //k3d_data_property(bool, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::no_constraint) m_previous_materials;

      /// Special variables and methods to generate polyhedra ...
      typedef std::vector<k3d::vector3> PolyhedronPoints;
      PolyhedronPoints m_PolyhedronPoints;

      std::string m_PolyhedronProperties;
      double m_PolyhedronScaling;

      typedef std::vector<int> PolyhedronPath;
      typedef std::vector<PolyhedronPath> PolyhedronPaths;
      PolyhedronPaths m_PolyhedronPaths;

      // Triangulated self-intersecting polygons have to
      // use their real color: reate a triangulation list
      std::vector<unsigned long> m_polygon_colors;
      typedef std::map<unsigned long, k3d::imaterial*> material_map_t;
      material_map_t m_materials;

      void LoadPolyhedron(int PolyhedronIndex, bool Dual)
      {
            // Clear out old data ...
            m_PolyhedronPoints.clear();
            m_PolyhedronPaths.clear();
            m_PolyhedronScaling = 0.5;

            // Solve polyhedron
            kaleido::Polyhedron* polyhedron = kaleido::kaleido(PolyhedronIndex);
            return_if_fail(polyhedron);

            double freeze = 0;

            if(!Dual)
                  vrmodel(polyhedron, polyhedron->v, polyhedron->V, polyhedron->f, polyhedron->F, false, AZ, EL, freeze);
            else
                  vrmodel(polyhedron, polyhedron->f, polyhedron->F, polyhedron->v, polyhedron->V, true, AZ, EL, freeze);

            // Write the polyhedron geometric properties string
            m_PolyhedronProperties = "";
            if(Dual && (polyhedron->hemi || polyhedron->index == kaleido::last_uniform - 1))
                  {
                        m_PolyhedronProperties += "infinite dual";

                        // Infinite duals have to scaled to a unit sphere
                        double maxlength = 0.0;
                        for(unsigned long i = 0; i < m_PolyhedronPoints.size(); i++)
                              {
                                    double length = m_PolyhedronPoints[i].Length();
                                    if(length > maxlength)
                                          maxlength = length;
                              }

                        m_PolyhedronScaling /= maxlength;
                  }

            // Free the polyhedron data ...
            delete polyhedron;
      }

      void PolyhedronPoint(kaleido::Kvector v)
      {
            m_PolyhedronPoints.push_back(k3d::vector3(v.x, v.y, v.z));
      }

      void PolyhedronFace()
      {
            m_PolyhedronPaths.push_back(PolyhedronPath());
      }

      void PolyhedronPathPoint(int PointIndex)
      {
            m_PolyhedronPaths.back().push_back(PointIndex);
      }

      void PolyhedronTriangle(int PointIndex1, int PointIndex2, int PointIndex3)
      {
            m_PolyhedronPaths.push_back(PolyhedronPath());

            m_PolyhedronPaths.back().push_back(PointIndex1);
            m_PolyhedronPaths.back().push_back(PointIndex2);
            m_PolyhedronPaths.back().push_back(PointIndex3);
      }

      void PolyhedronQuad(int PointIndex1, int PointIndex2, int PointIndex3, int PointIndex4)
      {
            m_PolyhedronPaths.push_back(PolyhedronPath());

            m_PolyhedronPaths.back().push_back(PointIndex1);
            m_PolyhedronPaths.back().push_back(PointIndex2);
            m_PolyhedronPaths.back().push_back(PointIndex3);
            m_PolyhedronPaths.back().push_back(PointIndex4);
      }

      void vrmodel(kaleido::Polyhedron *P, const std::vector<kaleido::Kvector>& v, int V, const std::vector<kaleido::Kvector>& f, int F, const bool dual, double azimuth, double elevation, double angle)
      {
            // Rotate frame
            static kaleido::Kvector X(1,0,0), Y(0,1,0), Z(0,0,1);
            kaleido::Kvector axis = kaleido::rotate (kaleido::rotate (X, Y, elevation), Z, azimuth);
            kaleido::Kvector rx = kaleido::rotate(X, axis, angle);
            kaleido::Kvector ry = kaleido::rotate(Y, axis, angle);
            kaleido::Kvector rz = kaleido::rotate(Z, axis, angle);

            // Rotate polyhedron
            std::vector<kaleido::Kvector> vb;
            for(int r = 0; r < V; r++)
            {
                  kaleido::Kvector rotv = v[r].x * rx + v[r].y * ry + v[r].z * rz;
                  vb.push_back(rotv);
                  PolyhedronPoint(rotv);
            }

            std::vector<kaleido::Kvector> fb;
            for(int r = 0; r < F; r++)
                  fb.push_back(f[r].x * rx + f[r].y * ry + f[r].z * rz);

      // Auxiliary vertices (needed because current VRML browsers cannot handle
      // non-simple paths, i.e., paths with self intersections):
      // Each non-simple face is assigned an auxiliary vertex. By connecting it to
      // the rest of the vertices the face is triangulated. The circum-center is used
      // for the regular star faces of uniform polyhedra. The in-center is used for
      // the pentagram (#79) and hexagram (#77) of the high-density snub duals, and
      // for the pentagrams (#40, #58) and hexagram (#52) of the stellated duals with
      // configuration (....)/2. Finally, the self-intersection of the crossed
      // parallelogram is used for duals with form p q r| with an even denominator.

      // This method do not work for the hemi-duals, whose faces are not star-shaped
      // and have two self-intersections each.
      // Thus, for each face we need six auxiliary vertices: The self intersections
      // and the terminal points of the truncations of the infinite edges. The ideal
      // vertices are listed, but are not used by the face-list.
      // Note that the face of the last dual (#80) is octagonal, and consists of two
      // quadrilaterals of the infinite type.

            int i, j, l, ll, ii;
            std::vector<int> hit;

            if(dual && P->even != -1)
                  hit.resize(F);

            for(i = 0; i < F; i++)
                  if(!dual && (kaleido::frac (P->n_sides[P->ftype[i]]),
                        kaleido::frax.d != 1 && kaleido::frax.d != kaleido::frax.n - 1) ||
                        dual && (P->K == 5 && P->D > 30 ||
                        kaleido::denominator (P->m_faces[0]) != 1)) {
                        // find the center of the face
                        double h;
                        if (!dual && P->hemi && !P->ftype[i])
                              h = 0;
                        else
                              h = P->minr / dot (fb[i],fb[i]);
                        PolyhedronPoint(h * fb[i]);
                  } else if (dual && P->even != -1) {
                        // find the self-intersection of a crossed parallelogram. hit is set if v0v1 intersects v2v3
                        kaleido::Kvector v0, v1, v2, v3, c0, c1, p;
                        double d0, d1;
                        v0 = vb[P->incid[0][i]];
                        v1 = vb[P->incid[1][i]];
                        v2 = vb[P->incid[2][i]];
                        v3 = vb[P->incid[3][i]];
                        d0 = sqrt (dot (v0 - v2, v0 - v2));
                        d1 = sqrt (dot (v1 - v3, v1 - v3));
                        c0 = d1 *(v0 + v2);
                        c1 = d0 *(v1 + v3);
                        p = (0.5 / (d0 + d1))*(c0 + c1);
                        PolyhedronPoint(p);
                        p = (p - v2)^(p - v3);
                        hit[i] = (dot (p, p) < 1e-6);
                  } else if (dual && P->hemi && P->index != kaleido::last_uniform - 1) {
                        // find the terminal points of the truncation and the
                        // self-intersections.
                        //  v23       v0       v21
                        //  |  \     /  \     /  |
                        //  |   v0123    v0321   |
                        //  |  /     \  /     \  |
                        //  v01       v2       v03

                        kaleido::Kvector v0, v1, v2, v3, v01, v03, v21, v23, v0123, v0321 ;
                        kaleido::Kvector u;
                        double t = 1.5; //truncation adjustment factor
                        j = !P->ftype[P->incid[0][i]];
                        v0 = vb[P->incid[j][i]];//real vertex
                        v1 = vb[P->incid[j+1][i]];//ideal vertex (unit kaleido::Kvector)
                        v2 = vb[P->incid[j+2][i]];//real
                        v3 = vb[P->incid[(j+3)%4][i]];//ideal
                        //compute intersections
                        //    this uses the following linear algebra:
                        //    v0123 = v0 + a v1 = v2 + b v3
                        //    v0 x v3 + a (v1 x v3) = v2 x v3
                        //    a (v1 x v3) = (v2 - v0) x v3
                        //    a (v1 x v3) . (v1 x v3) =
                        //          (v2 - v0) x v3 . (v1 x v3)
                        u = v1 ^ v3;
                        v0123 = v0 + (dot((v2 - v0)^v3,
                              u) / dot(u,u))*v1;
                        v0321 = v0 + (dot((v0 - v2)^v1,
                              u) / dot (u,u))*v3;
                        //compute truncations
                        v01 = v0 + t*(v0123 - v0);
                        v23 = v2 + t*(v0123 - v2);
                        v03 = v0 + t*(v0321 - v0);
                        v21 = v2 + t*(v0321 - v2);
                        PolyhedronPoint(v01);
                        PolyhedronPoint(v23);
                        PolyhedronPoint(v0123);
                        PolyhedronPoint(v03);
                        PolyhedronPoint(v21);
                        PolyhedronPoint(v0321);
                  } else if (dual && P->index == kaleido::last_uniform - 1) {
                        // find the terminal points of the truncation and the
                        // self-intersections.
                        //  v23       v0       v21
                        //  |  \     /  \     /  |
                        //  |   v0123    v0721   |
                        //  |  /     \  /     \  |
                        //  v01       v2       v07
                        //
                        //  v65       v4       v67
                        //  |  \     /  \     /  |
                        //  |   v4365    v4567   |
                        //  |  /     \  /     \  |
                        //  v43       v6       v45

                        kaleido::Kvector v0, v1, v2, v3, v4, v5, v6, v7, v01, v07, v21, v23;
                        kaleido::Kvector v43, v45, v65, v67, v0123, v0721, v4365, v4567;
                        double t = 1.5;//truncation adjustment factor
                        kaleido::Kvector u;
                        for (j = 0; j < 8; j++)
                              if (P->ftype[P->incid[j][i]] == 3)
                                    break;
                        v0 = vb[P->incid[j][i]];//real {5/3}
                        v1 = vb[P->incid[(j+1)%8][i]];//ideal
                        v2 = vb[P->incid[(j+2)%8][i]];//real {3}
                        v3 = vb[P->incid[(j+3)%8][i]];//ideal
                        v4 = vb[P->incid[(j+4)%8][i]];//real {5/2}
                        v5 = vb[P->incid[(j+5)%8][i]];//ideal
                        v6 = vb[P->incid[(j+6)%8][i]];//real {3/2}
                        v7 = vb[P->incid[(j+7)%8][i]];//ideal
                        //compute intersections
                        u = v1 ^ v3;
                        v0123 = v0 + (dot ((v2 - v0)^v3,
                              u) / dot (u,u))*v1;
                        u = v7 ^ v1;
                        v0721 = v0 + (dot ((v2 - v0)^v1,
                              u) / dot (u,u))*v7;
                        u = v5 ^ v7;
                        v4567 = v4 + (dot ((v6 - v4)^v7,
                              u) / dot (u,u))*v5;
                        u = v3 ^ v5;
                        v4365 = v4 + (dot ((v6 - v4)^v5,
                              u) / dot (u,u))*v3;
                        //compute truncations
                        v01 = v0 + t*(v0123 - v0);
                        v23 = v2 + t*(v0123 - v2);
                        v07 = v0 + t*(v0721 - v0);
                        v21 = v2 + t*(v0721 - v2);
                        v45 = v4 + t*(v4567 - v4);
                        v67 = v6 + t*(v4567 - v6);
                        v43 = v4 + t*(v4365 - v4);
                        v65 = v6 + t*(v4365 - v6);
                        PolyhedronPoint(v01);
                        PolyhedronPoint(v23);
                        PolyhedronPoint(v0123);
                        PolyhedronPoint(v07);
                        PolyhedronPoint(v21);
                        PolyhedronPoint(v0721);
                        PolyhedronPoint(v45);
                        PolyhedronPoint(v67);
                        PolyhedronPoint(v4567);
                        PolyhedronPoint(v43);
                        PolyhedronPoint(v65);
                        PolyhedronPoint(v4365);
                  }


            // Face list:
            // Each face is printed in a separate line, by listing the indices of its
            // vertices. In the non-simple case, the path is represented by the
            // triangulation, each triangle consists of two polyhedron vertices and one
            // auxiliary vertex.

            ii = V;
            for (i = 0; i < F; i++) {

                  PolyhedronFace();

                  if (dual) {
                        if(P->K == 5 && P->D > 30 || kaleido::denominator (P->m_faces[0]) != 1) {
                              for (j = 0; j < P->M - 1; j++)
                                    PolyhedronTriangle(P->incid[j][i], P->incid[j+1][i], ii);

                              PolyhedronTriangle(P->incid[j][i], P->incid[0][i], ii++);

                        } else if (P->even != -1) {
                              if(hit[i])
                              {
                                    PolyhedronTriangle(P->incid[3][i], P->incid[0][i], ii);
                                    PolyhedronTriangle(P->incid[1][i], P->incid[2][i], ii);
                              }
                              else
                              {
                                    PolyhedronTriangle(P->incid[0][i], P->incid[1][i], ii);
                                    PolyhedronTriangle(P->incid[2][i], P->incid[3][i], ii);
                              }
                              ii++;
                        } else if (P->hemi && P->index != kaleido::last_uniform - 1) {
                              j = !P->ftype[P->incid[0][i]];
                              PolyhedronTriangle(ii, ii+1, ii+2);
                              PolyhedronQuad(P->incid[j][i], ii+2, P->incid[j+2][i], ii+5);
                              PolyhedronTriangle(ii+3, ii+4, ii+5);
                              ii += 6;
                        } else if (P->index == kaleido::last_uniform - 1) {
                              for (j = 0; j < 8; j++)
                                    if (P->ftype[P->incid[j][i]] == 3)
                                          break;
                              PolyhedronTriangle(ii, ii+1, ii+2);
                              PolyhedronQuad(P->incid[j][i], ii+2, P->incid[(j+2)%8][i], ii+5);
                              PolyhedronTriangle(ii+3, ii+4, ii+5);

                              PolyhedronTriangle(ii+6, ii+7, ii+8);
                              PolyhedronQuad(P->incid[(j+4)%8][i], ii+8, P->incid[(j+6)%8][i], ii+11);
                              PolyhedronTriangle(ii+9, ii+10, ii+11);
                              ii += 12;
                        } else {
                              for (j = 0; j < P->M; j++)
                                    PolyhedronPathPoint(P->incid[j][i]);
                        }
                  } else {
                        int split = (kaleido::frac (P->n_sides[P->ftype[i]]),
                               kaleido::frax.d != 1 && kaleido::frax.d != kaleido::frax.n - 1);
                        int k = 0;
                        for (j = 0; j < V; j++) {
                              for (k = 0; k < P->M; k++)
                                    if (P->incid[k][j] == i)
                                          break;
                              if (k != P->M)
                                    break;
                        }
                        if (!split)
                              PolyhedronPathPoint(j);
                        ll = j;
                        for (l = P->adj[k][j]; l != j; l = P->adj[k][l]) {
                              for (k = 0; k < P->M; k++)
                                    if (P->incid[k][l] == i)
                                          break;
                              if (P->adj[k][l] == ll) k = kaleido::mod (k + 1 , P->M);
                              if (!split)
                                    PolyhedronPathPoint(l);
                              else
                                    PolyhedronTriangle(ll, l, ii);
                              ll = l;
                        }

                        if (split)
                              PolyhedronTriangle(ll, j, ii++);
                  }
            }

            // Face color indices - for polyhedra with multiple face types
            // For non-simple faces, the index is repeated as many times
            // as needed by the triangulation
            m_polygon_colors.clear();

            if(dual)
            {
                  m_polygon_colors.push_back(P->M);
            }
            else if(P->N == 1)
            {
                  for(i = 0; i < F; i++)
                        m_polygon_colors.push_back(kaleido::numerator(P->n_sides.front()));
            }
            else
            {
                  for(i = 0; i < F; i++)
                        if(kaleido::frac(P->n_sides[P->ftype[i]]), kaleido::frax.d == 1 || kaleido::frax.d == kaleido::frax.n - 1)
                              m_polygon_colors.push_back(kaleido::numerator(P->n_sides[P->ftype[i]]));
                        else
                              for(j = 0; j < kaleido::frax.n; j++)
                                    m_polygon_colors.push_back(kaleido::numerator(P->n_sides[P->ftype[i]]));
            }
      }
};

/////////////////////////////////////////////////////////////////////////////
// polyhedron_factory

k3d::iplugin_factory& polyhedron_factory()
{
      return polyhedron_implementation::get_factory();
}

} // namespace libk3dmesh



Generated by  Doxygen 1.6.0   Back to index