Commits

Ruben Martinez-Cantin committed 7264ad9

Working on visualization demo

Comments (0)

Files changed (19)

 add_dependencies(bo_oned bayesopt)
 TARGET_LINK_LIBRARIES(bo_oned bayesopt)
 
+#Display test
+find_package(GLUT)
+find_package(OpenGL)
+INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/matplotpp)
+enable_language(C)
+ADD_EXECUTABLE(bo_display 
+  ./app/bo_display.cpp 
+  ./matplotpp/matplotpp.cc 
+  ./matplotpp/gl2ps.c)
+add_dependencies(bo_display bayesopt)
+TARGET_LINK_LIBRARIES(bo_display bayesopt ${GLUT_LIBRARY} ${OPENGL_LIBRARY})
+
 #Branin
 ADD_EXECUTABLE(bo_branin ./app/bo_branin.cpp )
 add_dependencies(bo_branin bayesopt)

app/bo_display.cpp

+/*
+-------------------------------------------------------------------------
+   This file is part of BayesOpt, an efficient C++ library for 
+   Bayesian optimization.
+
+   Copyright (C) 2011-2013 Ruben Martinez-Cantin <rmcantin@unizar.es>
+ 
+   BayesOpt 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.
+
+   BayesOpt 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 BayesOpt.  If not, see <http://www.gnu.org/licenses/>.
+------------------------------------------------------------------------
+*/
+
+#include <valarray>
+#include "bayesoptcont.hpp"
+
+class ExampleOneD: public bayesopt::ContinuousModel
+{
+public:
+  ExampleOneD(size_t dim, bopt_params par):
+    ContinuousModel(dim,par) {}
+
+  double evaluateSample(const vectord& xin)
+  {
+    if (xin.size() > 1)
+      {
+	std::cout << "WARNING: This only works for 1D inputs." << std::endl
+		  << "WARNING: Using only first component." << std::endl;
+      }
+
+    double x = xin(0);
+    return sqr(x-0.5) + sin(20*x)*0.2;
+  };
+
+  bool checkReachability(const vectord &query)
+  {return true;};
+
+  inline double sqr( double x ){ return x*x; }
+
+  void printOptimal()
+  {
+    std::cout << "Optimal:" << 0.5 << std::endl;
+  }
+
+};
+
+#include "unistd.h"
+using namespace std;
+using namespace bayesopt;
+#include "matplotpp.h"
+
+int is_run=1;
+size_t state_ii = 0;
+BayesOptBase* GLOBAL_MODEL;
+
+class MP :public MatPlot{ 
+void DISPLAY(){
+    // Create test data 
+    int n=40;
+    dvec x,y,z,su,sl,c;
+    x=linspace(0,1,n);
+    y = x; z = x; su = x; sl = x; c= x;
+    vectord q(1);
+    for(int i=0;i<n;++i)
+      {
+	q(0) = x[i];
+	ProbabilityDistribution* pd = GLOBAL_MODEL->getSurrogateModel()->prediction(q);
+	y[i] = pd->getMean();
+	su[i] = y[i] + 2*pd->getStd();
+	sl[i] = y[i] - 2*pd->getStd();
+	c[i] = -GLOBAL_MODEL->evaluateCriteria(q);
+	z[i] = GLOBAL_MODEL->evaluateSample(q);
+      }
+    //plot
+    subplot(2,1,1);
+    title("press r to run and stop");
+    plot(x,z);
+    plot(x,y);
+    plot(x,su);
+    plot(x,sl);
+   
+    //mesh
+    subplot(2,1,2);
+    plot(x,c);
+
+}
+}mp;
+
+void display(){mp.display(); }
+void reshape(int w,int h){ mp.reshape(w,h); }
+void idle( void )
+{
+  if(is_run)
+    {
+      if (state_ii < 300)
+	++state_ii;
+      GLOBAL_MODEL->stepOptimization(state_ii); 
+    }
+  glutPostRedisplay();
+}
+
+void mouse(int button, int state, int x, int y ){ mp.mouse(button,state,x,y); }
+void motion(int x, int y ){mp.motion(x,y); }
+void passive(int x, int y ){mp.passivemotion(x,y); }
+void keyboard(unsigned char key, int x, int y){
+    mp.keyboard(key, x, y); 
+    if(key=='r'){ if(is_run==0){is_run=1;}else{is_run=0;}}
+}
+
+int main(int nargs, char *args[])
+{
+  size_t dim = 1;
+  bopt_params parameters = initialize_parameters_to_default();
+  parameters.n_init_samples = 10;
+  parameters.n_iterations = 300;
+  parameters.kernel.hp_mean[0] = 1.0;
+  parameters.kernel.hp_std[0] = 100.0;
+  parameters.kernel.n_hp = 1;
+  parameters.crit_name = "cHedge(cEI,cLCB,cExpReturn,cOptimisticSampling)";
+  parameters.epsilon = 0.0;
+
+  state_ii = 0;
+
+  ExampleOneD* opt = new ExampleOneD(dim,parameters);
+  vectord result(dim);
+  GLOBAL_MODEL = opt;
+  opt->initializeOptimization();
+  
+  glutInit(&nargs, args);
+  glutCreateWindow(100,100,500,400);
+  glutDisplayFunc( display );
+  glutReshapeFunc( reshape );
+  glutIdleFunc( idle );
+  glutMotionFunc( motion );
+  glutMouseFunc( mouse );
+  glutPassiveMotionFunc(passive);    
+  glutKeyboardFunc( keyboard );        
+  glutMainLoop();    
+
+  delete opt;
+  return 0;
+}

doxygen/reference.dox

 problem. Or, if the knowledge is not available, keep the model as
 general as possible (to avoid bias). 
 
-For that reason, the parameters are bundled in a structure or dictionary, depending on the API that we use. This is a brief explanation of every parameter:
-\li \b n_iterations: Maximum number of iterations of BayesOpt. Each iteration corresponds with a target function evaluation. This is related with the budget of the application 
-\li \b n_init_samples: BayesOpt requires an initial set of samples to learn a preliminar model of the target function. Each sample is also a function evaluation.
-\li \b verbose_level: (integer value) Verbose level 1,2,3 -> stderr, 4,5,6 -> log file
-\li \b log_filename: Filename of the log file (if applicable) 
+For that reason, the parameters are bundled in a structure or dictionary, depending on the API that we use. This is a brief explanation of every parameter
+
+\subsection budgetpar Budget parameters
+\li \b n_iterations: Maximum number of iterations of BayesOpt. Each iteration corresponds with a target function evaluation. This is related with the budget of the application [Default 300]
+\li \b n_inner_iterations: Maximum number of iterations of the inner optimization process. Each iteration corresponds with a criterion evaluation. The inner optimization results in the "most interest point" to run evaluate the target function. This is also used for the kernel hyperparameter computation. In order to scale the process for increasing dimensional spaces, the actual number of iterations is this number times the number of dimensions. [Default 500]
+\li \b n_init_samples: BayesOpt requires an initial set of samples to learn a preliminary model of the target function. Each sample is also a target function evaluation. [Default 30]
+\li \b n_iter_relearn: Although most of the parameters of the model are updated after every iteration, the kernel parameters cannot be updated continuously as it might crash the convergence. This represents the number of iterations between recomputing the kernel parameters. If it is set to 0, they are only learned after the initial set of samples. [Default 0]
+
+\subsection logpar Logging parameters
+\li \b verbose_level: (integer value) Verbose level 0,3 -> warnings, 1,4 -> general information, 2,5 -> debug information, any other value -> only errors. Levels < 3 send the messages to stdout. Levels > 4 send them to a log file. [Default 1].
+\li \b log_filename: Name of the log file (if applicable, verbose > 4)[Default "bayesopt.log"] 
+
+\subsection surrpar Surrogate model parameters
+\li \b surr_name: Name of the hierarchical surrogate function (nonparametric process and the hyperpriors on sigma and w). See section ??? for a detailed description. [Default S_GaussianProcess]
+\li \b sigma_s: Signal variance (if known) [Default 1.0]
+\li \b noise: Observation noise. For computer simulations or deterministic functions, it should be close to 0. However, to avoid numerical instability due to model inaccuracies, do not make it 0. [Default 0.0001]
+\li \b alpha, \b beta: Inverse-Gamma prior hyperparameters (if applicable) [Default 1.0, 1.0] 
+\li \b l_type: Learning method for the kernel hyperparameters. See section ??? for a detailed description [Default L_MAP]
+
+\subsection critpar Exploration parameters
+\li \b epsilon: According to some authors, it is recommendable to include an epsilon-greedy strategy to achieve near optimal convergence rates. Epsilon is the probability of performing a random (blind) evaluation of the target function [Default 0.0 (disabled)]
+\li \b crit_name: Name of the sample selection criterion. It is used to select which points to evaluate for each iteration of the optimization process. See section ??? for the different possibilities. [Default: "cEI"]
+
 \li \b theta: Kernel hyperparameters. Depends on the kernel selected.
 \li \b n_theta: Number of kernel hyperparameters.
 \li \b mu: Mean function hyperparameters. Depends on the mean function selected.
 \li \b n_mu: Number of mean function hyperparameters 
-\li \b alpha, \b beta, \b delta: Inverse-Gamma-Normal prior hyperparameters (if applicable) 
-\li \b noise: Observation noise. For computer simulations or deterministic functions, it should be close to 0. However, to avoid numerical instability due to model inaccuracies, do not make it 0.
-\li \b s_name: Name of the hierarchical surrogate function (gaussian process + hyperpriors)
 \li \b k_name: Name of the kernel function. 
-\li \b c_name: Name of the sample selection criterion. It is used to select which points to evaluate for each iteration of the optimization process.
 \li \b m_name: Name of the mean function (parametric part of the surrogate function).
 
 \section usage Using the library

include/bayesoptbase.hpp

      * @param bestPoint returns point with the optimum value in a ublas::vector.
      * @return 0 if terminate successfully, any other value otherwise
      */
-    virtual int optimize(vectord &bestPoint) = 0;
-  
+    int optimize(vectord &bestPoint);
+
+    /** 
+     * \brief Execute ONE step the optimization process of the function defined
+     * in evaluateSample.
+     * 
+     * @see evaluateSample
+     * @see checkReachability
+     *
+     * @param ii iteration number.
+     * @return 0 if terminate successfully, any other value otherwise
+     */  
+    int stepOptimization(size_t ii);
+
+    /** 
+     * Initialize the optimization process.
+     * @return error_code
+     */
+    virtual int initializeOptimization() = 0;
+
+    /** 
+     * Once the optimization has been perfomed, return the optimal point.
+     */
+    virtual vectord getFinalResult() = 0;
+
+
     /** 
      * \brief Function that defines the actual function to be optimized.
      * This function must be modified (overriden) according to the
     virtual bool checkReachability( const vectord &query )
     { return true; };
 
-    
-    /** 
-     * Getter to access the underlying surrogate function.
-     * Used mainly for visualization purposed.
-     * 
-     * @return Pointer to the surrogate function object.
-     */  
-    // NonParametricProcess* getSurrogateFunctionPointer()
-    // { return mGP; };
-    
-  protected:
+
     /** 
      * \brief Evaluate the criteria considering if the query is
      * reachable or not.  This is a way to include non-linear
       return (*mCrit)(query);
     };
 
+
+    NonParametricProcess* getSurrogateModel()
+    { return mGP.get(); };
+  protected:
+    /** 
+     * Print data for every step according to the verbose level
+     * 
+     * @param iteration 
+     * @param xNext 
+     * @param yNext 
+     * 
+     * @return error code
+     */
+    virtual int plotStepData(size_t iteration, const vectord& xNext,
+		     double yNext) = 0;
+
+
+    /** 
+     * \brief Wrapper for the target function normalize in the hypercube
+     * [0,1]
+     * @param query point to evaluate in [0,1] hypercube
+     * @return actual return value of the target function
+     */
+    virtual double evaluateSampleInternal( const vectord &query ) = 0;
+
+
     /** 
      * \brief Returns the optimal point acording to certain criteria
      * @see evaluateCriteria
      */
     int nextPoint( vectord &Xnext );  
 
+
+
   protected:
     boost::scoped_ptr<NonParametricProcess> mGP;  ///< Pointer to surrogate model
     boost::scoped_ptr<Criteria> mCrit;                    ///< Metacriteria model

include/bayesoptcont.hpp

     virtual ~ContinuousModel();
   
     /** 
-     * \brief Execute the optimization process of the function defined in
-     * evaluateSample.  If no bounding box is defined, we assume that
-     * the function is defined in the [0,1] hypercube, as a normalized
-     * representation of the bound constrains.
-     * 
-     * @see scaleInput
-     * @see evaluateSample
-     *
-     * @param bestPoint returns the optimum value in a ublas::vector.
-     * @return 0 if terminate successfully, nonzero otherwise
+     * Initialize the optimization process.
+     * @return error_code
      */
-    int optimize(vectord &bestPoint);
-    
+    int initializeOptimization();
 
+    /** 
+     * Once the optimization has been perfomed, return the optimal point.
+     */
+    vectord getFinalResult();
 
     /** 
      * \brief Sets the bounding box. 
      * @param query point to evaluate in [0,1] hypercube
      * @return actual return value of the target function
      */
-    inline double evaluateNormalizedSample( const vectord &query )
+    inline double evaluateSampleInternal( const vectord &query )
     { 
       vectord unnormalizedQuery = mBB->unnormalizeVector(query);
       return evaluateSample(unnormalizedQuery);
-    }; // evaluateNormalizedSample
+    }; // evaluateSampleInternal
 
     /** 
      * \brief Wrapper of the innerOptimization class to find the optimal

include/bayesoptdisc.hpp

     virtual ~DiscreteModel();
 
     /** 
-     * Execute the optimization process of the function defined in 
-     * evaluateSample.
-     * We assume that the function is defined in the [0,1] hypercube, as a 
-     * normalized representation of the bound constrains.
-     * 
-     * @see scaleInput
-     * @see evaluateSample
-     *
-     * @param bestPoint returns the optimum value in a ublas::vector defined in 
-     * the hypercube [0,1], it might also be used as an initial point
-     * 
-     * @return 1 if terminate successfully, 0 otherwise
+     * Initialize the optimization process.
+     * @return error_code
      */
-    int optimize( vectord &bestPoint);
-  
+    int initializeOptimization();
 
+    /** 
+     * Once the optimization has been perfomed, return the optimal point.
+     */
+    vectord getFinalResult();
+
+    
   protected:
     
     
      */
     int sampleInitialPoints();
 
+    /** 
+     * \brief Wrapper for the target function normalize in the hypercube
+     * [0,1]
+     * @param query point to evaluate in [0,1] hypercube
+     * @return actual return value of the target function
+     */
+    inline double evaluateSampleInternal( const vectord &query )
+    { return evaluateSample(query); }; 
+
     int findOptimal(vectord &xOpt);
 
   protected:

include/parameters.h

 
   /** \brief Configuration parameters */
   typedef struct {
-
     size_t n_iterations;         /**< Maximum BayesOpt evaluations (budget) */
     size_t n_inner_iterations;   /**< Maximum inner optimizer evaluations */
     size_t n_init_samples;       /**< Number of samples before optimization */

matlab/bayesoptextras.h

 }
 
 static bopt_params load_parameters(const mxArray* params)
-{
+{  
+  
+  /* See parameters.h for the available options */
+  
   char log_str[100], k_s_str[100];
   char c_str[100], s_str[100], k_str[100], m_str[100], l_str[100];
-  size_t n_theta, n_mu;
+  size_t n_hp_test, n_coef_test;
 
   bopt_params parameters = initialize_parameters_to_default();
-  n_theta = parameters.kernel.n_theta;
-  n_mu = parameters.mean.n_mu;
+
+  n_hp_test = parameters.kernel.n_hp;
+  n_coef_test = parameters.mean.n_coef;
 
   struct_size(params,"n_iterations", &parameters.n_iterations);
   struct_size(params,"n_inner_iterations", &parameters.n_inner_iterations);
-  struct_size(params, "n_init_iterations", &parameters.n_init_samples);
+  struct_size(params, "n_init_samples", &parameters.n_init_samples);
+  struct_size(params, "n_iter_relearn", &parameters.n_iter_relearn);
+  
   struct_size(params, "verbose_level", &parameters.verbose_level);
-
-  struct_value(params, "alpha", &parameters.alpha);
-  struct_value(params, "beta",  &parameters.beta);
-  struct_value(params, "noise", &parameters.noise);
-
-  struct_array(params, "theta", &parameters.kernel.n_theta, 
-	       &parameters.kernel.theta[0]);
-
-  struct_array(params, "s_theta", &n_theta, 
-	       &parameters.kernel.s_theta[0]);
-
-  CHECK0(parameters.kernel.n_theta == n_theta, 
-	 "Error processing kernel parameters");
-
-  struct_array(params, "mu", &parameters.mean.n_mu, 
-	       &parameters.mean.mu[0]);
-
-  struct_array(params, "s_mu", &n_mu, 
-	       &parameters.mean.s_mu[0]);
-
-  CHECK0(parameters.mean.n_mu == n_mu, 
-	 "Error processing mean parameters");
-
-  /* Extra configuration
-  See parameters.h for the available options */
-
   struct_string(params, "log_filename", parameters.log_filename);
-  struct_string(params, "kernel_name", parameters.kernel.name);
-  struct_string(params, "mean_name", parameters.mean.name);
-  struct_string(params, "crit_name", parameters.crit_name);
 
   strcpy( s_str, surrogate2str(parameters.surr_name));
   struct_string(params, "surr_name", s_str);
   parameters.surr_name = str2surrogate(s_str);
 
+  struct_value(params, "sigma_s", &parameters.sigma_s);
+  struct_value(params, "noise", &parameters.noise);
+  struct_value(params, "alpha", &parameters.alpha);
+  struct_value(params, "beta",  &parameters.beta);
+  
+
   strcpy( l_str, learn2str(parameters.l_type));
   struct_string(params, "l_type", l_str);
   parameters.l_type = str2learn(l_str);
 
+  struct_value(params, "epsilon",  &parameters.epsilon);
+
+  struct_string(params, "crit_name", parameters.crit_name);
+  struct_array(params, "crit_params", &parameters.n_crit_params, 
+	       &parameters.crit_params[0]);
+
+  /* Kernel parameters */
+  struct_string(params, "kernel_name", parameters.kernel.name);
+  struct_array(params, "kernel_hp_mean", &parameters.kernel.n_hp, 
+	       &parameters.kernel.hp_mean[0]);
+  struct_array(params, "kernel_hp_std", &n_hp_test, 
+	       &parameters.kernel.hp_std[0]);
+
+  CHECK0(parameters.kernel.n_hp == n_hp_test, 
+	 "Error processing kernel parameters");
+
+  /* Mean function parameters */
+  struct_string(params, "mean_name", parameters.mean.name);
+  struct_array(params, "mean_coef_mean", &parameters.mean.n_coef, 
+	       &parameters.mean.coef_mean[0]);
+
+  struct_array(params, "mean_coef_std", &n_coef_test, 
+	       &parameters.mean.coef_std[0]);
+
+  CHECK0(parameters.mean.n_coef == n_coef_test, 
+	 "Error processing mean parameters");
+
+
   return parameters;
 }
 
 params.surr_name = 'GAUSSIAN_PROCESS';
 params.noise = 0.005;
 params.kernel_name = 'kMaternISO3';
-params.theta = [0.5];
-params.s_theta = [100];
+params.kernel_hp_mean = [0.5];
+params.kernel_hp_std = [10];
 params.verbose_level = 1;
 params.log_filename = 'matbopt.log';
 
 
 lb = ones(n,1)*pi/2;
 ub = ones(n,1)*pi;
+fun = 'michalewicz';
+
+disp('Continuous optimization');pause;
+tic;
+bayesopt(fun,n,params,lb,ub)
+toc;
+
+disp('Discrete optimization');pause;
+% The set of points must be nDim x nPoints.
+xset = repmat((ub-lb),1,100) .* rand(n,100) - repmat(lb,1,100);
 
 tic;
-bayesopt('michalewicz',n,params,lb,ub)
-toc;
-
-% The set of points must be nDim x nPoints.
-xset = rand(n,100);
-
-tic;
-bayesoptdisc('quadratic',xset, params);
+bayesoptdisc(fun,xset, params);
 toc;
+MATPLOT++ (MATLAB-like plotting tool in C++)
+Copyright (c) 2011 Yuichi Katori All Rights Reserved
+Author: Yuichi Katori (yuichi.katori@gmail.com)
+
+Getting Started
+------------------------------------------
+
+To install a required library on ubuntu
+
+$ sudo apt-get install freeglut3-dev
+
+To compile MATPLOT++ library, type as following:
+
+$ tar xvf matplotpp-0.3.x.tar
+$ cd matplotpp
+$ make
+
+To compile examples:
+
+$ cd examples
+$ make

matplotpp/gl2ps.c

+/*
+ * GL2PS, an OpenGL to PostScript Printing Library
+ * Copyright (C) 1999-2009 C. Geuzaine
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of either:
+ *
+ * a) the GNU Library General Public License as published by the Free
+ * Software Foundation, either version 2 of the License, or (at your
+ * option) any later version; or
+ *
+ * b) the GL2PS License as published by Christophe Geuzaine, 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 either
+ * the GNU Library General Public License or the GL2PS License for
+ * more details.
+ *
+ * You should have received a copy of the GNU Library General Public
+ * License along with this library in the file named "COPYING.LGPL";
+ * if not, write to the Free Software Foundation, Inc., 675 Mass Ave,
+ * Cambridge, MA 02139, USA.
+ *
+ * You should have received a copy of the GL2PS License with this
+ * library in the file named "COPYING.GL2PS"; if not, I will be glad
+ * to provide one.
+ *
+ * For the latest info about gl2ps and a full list of contributors,
+ * see http://www.geuz.org/gl2ps/.
+ *
+ * Please report all bugs and problems to <gl2ps@geuz.org>.
+ */
+
+#include "gl2ps.h"
+
+#include <math.h>
+#include <string.h>
+#include <sys/types.h>
+#include <stdarg.h>
+#include <time.h>
+#include <float.h>
+
+#if defined(GL2PS_HAVE_ZLIB)
+#include <zlib.h>
+#endif
+
+#if defined(GL2PS_HAVE_LIBPNG)
+#include <png.h>
+#endif
+
+/********************************************************************* 
+ *
+ * Private definitions, data structures and prototypes
+ *
+ *********************************************************************/
+
+/* Magic numbers (assuming that the order of magnitude of window
+   coordinates is 10^3) */
+
+#define GL2PS_EPSILON       5.0e-3F
+#define GL2PS_ZSCALE        1000.0F
+#define GL2PS_ZOFFSET       5.0e-2F
+#define GL2PS_ZOFFSET_LARGE 20.0F
+#define GL2PS_ZERO(arg)     (fabs(arg) < 1.e-20)
+
+/* Primitive types */
+
+#define GL2PS_NO_TYPE          -1
+#define GL2PS_TEXT             1
+#define GL2PS_POINT            2
+#define GL2PS_LINE             3
+#define GL2PS_QUADRANGLE       4
+#define GL2PS_TRIANGLE         5
+#define GL2PS_PIXMAP           6
+#define GL2PS_IMAGEMAP         7
+#define GL2PS_IMAGEMAP_WRITTEN 8
+#define GL2PS_IMAGEMAP_VISIBLE 9
+#define GL2PS_SPECIAL          10
+
+/* BSP tree primitive comparison */
+
+#define GL2PS_COINCIDENT  1
+#define GL2PS_IN_FRONT_OF 2
+#define GL2PS_IN_BACK_OF  3
+#define GL2PS_SPANNING    4
+
+/* 2D BSP tree primitive comparison */
+
+#define GL2PS_POINT_COINCIDENT 0
+#define GL2PS_POINT_INFRONT    1
+#define GL2PS_POINT_BACK       2
+
+/* Internal feedback buffer pass-through tokens */
+
+#define GL2PS_BEGIN_OFFSET_TOKEN   1
+#define GL2PS_END_OFFSET_TOKEN     2
+#define GL2PS_BEGIN_BOUNDARY_TOKEN 3
+#define GL2PS_END_BOUNDARY_TOKEN   4
+#define GL2PS_BEGIN_STIPPLE_TOKEN  5
+#define GL2PS_END_STIPPLE_TOKEN    6
+#define GL2PS_POINT_SIZE_TOKEN     7
+#define GL2PS_LINE_WIDTH_TOKEN     8
+#define GL2PS_BEGIN_BLEND_TOKEN    9
+#define GL2PS_END_BLEND_TOKEN      10
+#define GL2PS_SRC_BLEND_TOKEN      11
+#define GL2PS_DST_BLEND_TOKEN      12
+#define GL2PS_IMAGEMAP_TOKEN       13
+#define GL2PS_DRAW_PIXELS_TOKEN    14
+#define GL2PS_TEXT_TOKEN           15
+
+typedef enum {
+  T_UNDEFINED    = -1,
+  T_CONST_COLOR  = 1,
+  T_VAR_COLOR    = 1<<1,
+  T_ALPHA_1      = 1<<2,
+  T_ALPHA_LESS_1 = 1<<3,
+  T_VAR_ALPHA    = 1<<4
+} GL2PS_TRIANGLE_PROPERTY;
+
+typedef GLfloat GL2PSxyz[3];
+typedef GLfloat GL2PSplane[4];
+
+typedef struct _GL2PSbsptree2d GL2PSbsptree2d;
+
+struct _GL2PSbsptree2d {
+  GL2PSplane plane;
+  GL2PSbsptree2d *front, *back;
+};
+
+typedef struct {
+  GLint nmax, size, incr, n;
+  char *array;
+} GL2PSlist;
+
+typedef struct _GL2PSbsptree GL2PSbsptree;
+
+struct _GL2PSbsptree {
+  GL2PSplane plane;
+  GL2PSlist *primitives;
+  GL2PSbsptree *front, *back;
+};
+
+typedef struct {
+  GL2PSxyz xyz;
+  GL2PSrgba rgba;
+} GL2PSvertex;
+
+typedef struct {
+  GL2PSvertex vertex[3];
+  int prop;
+} GL2PStriangle;
+
+typedef struct {
+  GLshort fontsize;
+  char *str, *fontname;
+  /* Note: for a 'special' string, 'alignment' holds the format
+     (PostScript, PDF, etc.) of the special string */
+  GLint alignment;
+  GLfloat angle;
+} GL2PSstring;
+
+typedef struct {
+  GLsizei width, height;
+  /* Note: for an imagemap, 'type' indicates if it has already been
+     written to the file or not, and 'format' indicates if it is
+     visible or not */
+  GLenum format, type;
+  GLfloat *pixels;
+} GL2PSimage;
+
+typedef struct _GL2PSimagemap GL2PSimagemap;
+
+struct _GL2PSimagemap {
+  GL2PSimage *image;
+  GL2PSimagemap *next;
+};
+
+typedef struct {
+  GLshort type, numverts;
+  GLushort pattern;
+  char boundary, offset, culled;
+  GLint factor;
+  GLfloat width;
+  GL2PSvertex *verts;
+  union {
+    GL2PSstring *text;
+    GL2PSimage *image;
+  } data;
+} GL2PSprimitive;
+
+typedef struct {
+#if defined(GL2PS_HAVE_ZLIB)
+  Bytef *dest, *src, *start;
+  uLongf destLen, srcLen;
+#else
+  int dummy;
+#endif
+} GL2PScompress;
+
+typedef struct{
+  GL2PSlist* ptrlist;
+  int gsno, fontno, imno, shno, maskshno, trgroupno;
+  int gsobjno, fontobjno, imobjno, shobjno, maskshobjno, trgroupobjno;
+} GL2PSpdfgroup;
+
+typedef struct {
+  /* General */
+  GLint format, sort, options, colorsize, colormode, buffersize;
+  char *title, *producer, *filename;
+  GLboolean boundary, blending;
+  GLfloat *feedback, offset[2], lastlinewidth;
+  GLint viewport[4], blendfunc[2], lastfactor;
+  GL2PSrgba *colormap, lastrgba, threshold, bgcolor;
+  GLushort lastpattern;
+  GL2PSvertex lastvertex;
+  GL2PSlist *primitives, *auxprimitives;
+  FILE *stream;
+  GL2PScompress *compress;
+  GLboolean header;
+
+  /* BSP-specific */
+  GLint maxbestroot;
+
+  /* Occlusion culling-specific */
+  GLboolean zerosurfacearea;
+  GL2PSbsptree2d *imagetree;
+  GL2PSprimitive *primitivetoadd;
+  
+  /* PDF-specific */
+  int streamlength;
+  GL2PSlist *pdfprimlist, *pdfgrouplist;
+  int *xreflist;
+  int objects_stack; /* available objects */
+  int extgs_stack; /* graphics state object number */
+  int font_stack; /* font object number */
+  int im_stack; /* image object number */
+  int trgroupobjects_stack; /* xobject numbers */
+  int shader_stack; /* shader object numbers */
+  int mshader_stack; /* mask shader object numbers */
+
+  /* for image map list */
+  GL2PSimagemap *imagemap_head;
+  GL2PSimagemap *imagemap_tail;
+} GL2PScontext;
+
+typedef struct {
+  void  (*printHeader)(void);
+  void  (*printFooter)(void);
+  void  (*beginViewport)(GLint viewport[4]);
+  GLint (*endViewport)(void);
+  void  (*printPrimitive)(void *data);
+  void  (*printFinalPrimitive)(void);
+  const char *file_extension;
+  const char *description;
+} GL2PSbackend;
+
+/* The gl2ps context. gl2ps is not thread safe (we should create a
+   local GL2PScontext during gl2psBeginPage) */
+
+static GL2PScontext *gl2ps = NULL;
+
+/* Need to forward-declare this one */
+
+static GLint gl2psPrintPrimitives(void);
+
+/********************************************************************* 
+ *
+ * Utility routines
+ *
+ *********************************************************************/
+
+static void gl2psMsg(GLint level, const char *fmt, ...)
+{
+  va_list args;
+
+  if(!(gl2ps->options & GL2PS_SILENT)){
+    switch(level){
+    case GL2PS_INFO : fprintf(stderr, "GL2PS info: "); break;
+    case GL2PS_WARNING : fprintf(stderr, "GL2PS warning: "); break;
+    case GL2PS_ERROR : fprintf(stderr, "GL2PS error: "); break;
+    }
+    va_start(args, fmt);
+    vfprintf(stderr, fmt, args); 
+    va_end(args);
+    fprintf(stderr, "\n");
+  }
+  /* if(level == GL2PS_ERROR) exit(1); */
+}
+
+static void *gl2psMalloc(size_t size)
+{
+  void *ptr;
+
+  if(!size) return NULL;
+  ptr = malloc(size);
+  if(!ptr){
+    gl2psMsg(GL2PS_ERROR, "Couldn't allocate requested memory");
+    return NULL;
+  }
+  return ptr;
+}
+
+static void *gl2psRealloc(void *ptr, size_t size)
+{
+  if(!size) return NULL;
+  ptr = realloc(ptr, size);
+  if(!ptr){
+    gl2psMsg(GL2PS_ERROR, "Couldn't reallocate requested memory");
+    return NULL;
+  }
+  return ptr;
+}
+
+static void gl2psFree(void *ptr)
+{
+  if(!ptr) return;
+  free(ptr);
+}
+
+static size_t gl2psWriteBigEndian(unsigned long data, size_t bytes)
+{
+  size_t i;
+  size_t size = sizeof(unsigned long);
+  for(i = 1; i <= bytes; ++i){
+    fputc(0xff & (data >> (size-i) * 8), gl2ps->stream);
+  }
+  return bytes;
+}
+
+/* zlib compression helper routines */
+
+#if defined(GL2PS_HAVE_ZLIB)
+
+static void gl2psSetupCompress(void)
+{
+  gl2ps->compress = (GL2PScompress*)gl2psMalloc(sizeof(GL2PScompress));
+  gl2ps->compress->src = NULL;
+  gl2ps->compress->start = NULL;
+  gl2ps->compress->dest = NULL;
+  gl2ps->compress->srcLen = 0;
+  gl2ps->compress->destLen = 0;
+}
+
+static void gl2psFreeCompress(void)
+{
+  if(!gl2ps->compress)
+    return;
+  gl2psFree(gl2ps->compress->start);
+  gl2psFree(gl2ps->compress->dest);
+  gl2ps->compress->src = NULL;
+  gl2ps->compress->start = NULL;
+  gl2ps->compress->dest = NULL;
+  gl2ps->compress->srcLen = 0;
+  gl2ps->compress->destLen = 0;
+}
+
+static int gl2psAllocCompress(unsigned int srcsize)
+{
+  gl2psFreeCompress();
+  
+  if(!gl2ps->compress || !srcsize)
+    return GL2PS_ERROR;
+  
+  gl2ps->compress->srcLen = srcsize;
+  gl2ps->compress->destLen = (int)ceil(1.001 * gl2ps->compress->srcLen + 12);
+  gl2ps->compress->src = (Bytef*)gl2psMalloc(gl2ps->compress->srcLen);
+  gl2ps->compress->start = gl2ps->compress->src;
+  gl2ps->compress->dest = (Bytef*)gl2psMalloc(gl2ps->compress->destLen);
+  
+  return GL2PS_SUCCESS;
+}
+
+static void *gl2psReallocCompress(unsigned int srcsize)
+{
+  if(!gl2ps->compress || !srcsize)
+    return NULL;
+  
+  if(srcsize < gl2ps->compress->srcLen)
+    return gl2ps->compress->start;
+  
+  gl2ps->compress->srcLen = srcsize;
+  gl2ps->compress->destLen = (int)ceil(1.001 * gl2ps->compress->srcLen + 12);
+  gl2ps->compress->src = (Bytef*)gl2psRealloc(gl2ps->compress->src, 
+                                              gl2ps->compress->srcLen);
+  gl2ps->compress->start = gl2ps->compress->src;
+  gl2ps->compress->dest = (Bytef*)gl2psRealloc(gl2ps->compress->dest, 
+                                               gl2ps->compress->destLen);
+  
+  return gl2ps->compress->start;
+}
+
+static size_t gl2psWriteBigEndianCompress(unsigned long data, size_t bytes)
+{
+  size_t i;
+  size_t size = sizeof(unsigned long);
+  for(i = 1; i <= bytes; ++i){
+    *gl2ps->compress->src = (Bytef)(0xff & (data >> (size-i) * 8));
+    ++gl2ps->compress->src;
+  }
+  return bytes;
+}
+
+static int gl2psDeflate(void)
+{
+  /* For compatibility with older zlib versions, we use compress(...)
+     instead of compress2(..., Z_BEST_COMPRESSION) */
+  return compress(gl2ps->compress->dest, &gl2ps->compress->destLen, 
+                  gl2ps->compress->start, gl2ps->compress->srcLen);  
+}
+
+#endif
+
+static int gl2psPrintf(const char* fmt, ...)
+{
+  int ret;
+  va_list args;
+
+#if defined(GL2PS_HAVE_ZLIB)
+  unsigned int oldsize = 0;
+  static char buf[1000];
+  if(gl2ps->options & GL2PS_COMPRESS){
+    va_start(args, fmt);
+    ret = vsprintf(buf, fmt, args);
+    va_end(args);
+    oldsize = gl2ps->compress->srcLen;
+    gl2ps->compress->start = (Bytef*)gl2psReallocCompress(oldsize + ret);
+    memcpy(gl2ps->compress->start+oldsize, buf, ret);
+    ret = 0;
+  }
+  else{
+#endif
+    va_start(args, fmt);
+    ret = vfprintf(gl2ps->stream, fmt, args);
+    va_end(args);
+#if defined(GL2PS_HAVE_ZLIB)
+  }
+#endif
+  return ret;
+}
+
+static void gl2psPrintGzipHeader()
+{
+#if defined(GL2PS_HAVE_ZLIB)
+  char tmp[10] = {'\x1f', '\x8b', /* magic numbers: 0x1f, 0x8b */
+                  8, /* compression method: Z_DEFLATED */
+                  0, /* flags */
+                  0, 0, 0, 0, /* time */
+                  2, /* extra flags: max compression */
+                  '\x03'}; /* OS code: 0x03 (Unix) */
+
+  if(gl2ps->options & GL2PS_COMPRESS){
+    gl2psSetupCompress();
+    /* add the gzip file header */
+    fwrite(tmp, 10, 1, gl2ps->stream);
+  }
+#endif  
+}
+
+static void gl2psPrintGzipFooter()
+{
+#if defined(GL2PS_HAVE_ZLIB)
+  int n;
+  uLong crc, len;
+  char tmp[8];
+
+  if(gl2ps->options & GL2PS_COMPRESS){
+    if(Z_OK != gl2psDeflate()){
+      gl2psMsg(GL2PS_ERROR, "Zlib deflate error");
+    }
+    else{
+      /* determine the length of the header in the zlib stream */
+      n = 2; /* CMF+FLG */
+      if(gl2ps->compress->dest[1] & (1<<5)){
+        n += 4; /* DICTID */
+      }
+      /* write the data, without the zlib header and footer */
+      fwrite(gl2ps->compress->dest+n, gl2ps->compress->destLen-(n+4), 
+             1, gl2ps->stream);
+      /* add the gzip file footer */
+      crc = crc32(0L, gl2ps->compress->start, gl2ps->compress->srcLen);
+      for(n = 0; n < 4; ++n){
+        tmp[n] = (char)(crc & 0xff);
+        crc >>= 8;
+      }
+      len = gl2ps->compress->srcLen;
+      for(n = 4; n < 8; ++n){
+        tmp[n] = (char)(len & 0xff);
+        len >>= 8;
+      }
+      fwrite(tmp, 8, 1, gl2ps->stream);
+    }
+    gl2psFreeCompress();
+    gl2psFree(gl2ps->compress);
+    gl2ps->compress = NULL;
+  }
+#endif 
+}
+
+/* The list handling routines */
+
+static void gl2psListRealloc(GL2PSlist *list, GLint n)
+{
+  if(!list){
+    gl2psMsg(GL2PS_ERROR, "Cannot reallocate NULL list");
+    return;
+  }
+  if(n <= 0) return;
+  if(!list->array){
+    list->nmax = n;
+    list->array = (char*)gl2psMalloc(list->nmax * list->size);
+  }
+  else{
+    if(n > list->nmax){
+      list->nmax = ((n - 1) / list->incr + 1) * list->incr;
+      list->array = (char*)gl2psRealloc(list->array,
+                                        list->nmax * list->size);
+    }
+  }
+}
+
+static GL2PSlist *gl2psListCreate(GLint n, GLint incr, GLint size)
+{
+  GL2PSlist *list;
+
+  if(n < 0) n = 0;
+  if(incr <= 0) incr = 1;
+  list = (GL2PSlist*)gl2psMalloc(sizeof(GL2PSlist));
+  list->nmax = 0;
+  list->incr = incr;
+  list->size = size;
+  list->n = 0;
+  list->array = NULL;
+  gl2psListRealloc(list, n);
+  return list;
+}
+
+static void gl2psListReset(GL2PSlist *list)
+{
+  if(!list) return;
+  list->n = 0;
+}
+
+static void gl2psListDelete(GL2PSlist *list)
+{
+  if(!list) return;  
+  gl2psFree(list->array);
+  gl2psFree(list);
+}
+
+static void gl2psListAdd(GL2PSlist *list, void *data)
+{
+  if(!list){
+    gl2psMsg(GL2PS_ERROR, "Cannot add into unallocated list");
+    return;
+  }
+  list->n++;
+  gl2psListRealloc(list, list->n);
+  memcpy(&list->array[(list->n - 1) * list->size], data, list->size);
+}
+
+static int gl2psListNbr(GL2PSlist *list)
+{
+  if(!list)
+    return 0;
+  return list->n;
+}
+
+static void *gl2psListPointer(GL2PSlist *list, GLint index)
+{
+  if(!list){
+    gl2psMsg(GL2PS_ERROR, "Cannot point into unallocated list");
+    return NULL;
+  }
+  if((index < 0) || (index >= list->n)){
+    gl2psMsg(GL2PS_ERROR, "Wrong list index in gl2psListPointer");
+    return NULL;
+  }
+  return &list->array[index * list->size];
+}
+
+static void gl2psListSort(GL2PSlist *list,
+                          int (*fcmp)(const void *a, const void *b))
+{
+  if(!list)
+    return;
+  qsort(list->array, list->n, list->size, fcmp);
+}
+
+static void gl2psListAction(GL2PSlist *list, void (*action)(void *data))
+{
+  GLint i;
+
+  for(i = 0; i < gl2psListNbr(list); i++){
+    (*action)(gl2psListPointer(list, i));
+  }
+}
+
+static void gl2psListActionInverse(GL2PSlist *list, void (*action)(void *data))
+{
+  GLint i;
+
+  for(i = gl2psListNbr(list); i > 0; i--){
+    (*action)(gl2psListPointer(list, i-1));
+  }
+}
+
+#if defined(GL2PS_HAVE_LIBPNG)
+
+static void gl2psListRead(GL2PSlist *list, int index, void *data)
+{
+  if((index < 0) || (index >= list->n))
+    gl2psMsg(GL2PS_ERROR, "Wrong list index in gl2psListRead");
+  memcpy(data, &list->array[index * list->size], list->size);
+}
+
+static void gl2psEncodeBase64Block(unsigned char in[3], unsigned char out[4], int len)
+{
+  static const char cb64[] = 
+    "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";
+
+  out[0] = cb64[ in[0] >> 2 ];
+  out[1] = cb64[ ((in[0] & 0x03) << 4) | ((in[1] & 0xf0) >> 4) ];
+  out[2] = (len > 1) ? cb64[ ((in[1] & 0x0f) << 2) | ((in[2] & 0xc0) >> 6) ] : '=';
+  out[3] = (len > 2) ? cb64[ in[2] & 0x3f ] : '=';
+}
+
+static void gl2psListEncodeBase64(GL2PSlist *list)
+{
+  unsigned char *buffer, in[3], out[4];
+  int i, n, index, len;
+
+  n = list->n * list->size;
+  buffer = (unsigned char*)gl2psMalloc(n * sizeof(unsigned char));
+  memcpy(buffer, list->array, n * sizeof(unsigned char));
+  gl2psListReset(list);
+
+  index = 0;
+  while(index < n) {
+    len = 0;
+    for(i = 0; i < 3; i++) {
+      if(index < n){
+        in[i] = buffer[index];
+        len++;
+      }
+      else{
+        in[i] = 0;
+      }
+      index++;
+    }
+    if(len) {
+      gl2psEncodeBase64Block(in, out, len);
+      for(i = 0; i < 4; i++)
+        gl2psListAdd(list, &out[i]);
+    }
+  }
+  gl2psFree(buffer);
+}
+
+#endif
+
+/* Helpers for rgba colors */
+
+static GLboolean gl2psSameColor(GL2PSrgba rgba1, GL2PSrgba rgba2)
+{
+  if(!GL2PS_ZERO(rgba1[0] - rgba2[0]) ||
+     !GL2PS_ZERO(rgba1[1] - rgba2[1]) ||
+     !GL2PS_ZERO(rgba1[2] - rgba2[2]))
+    return GL_FALSE;
+  return GL_TRUE;
+}
+  
+static GLboolean gl2psVertsSameColor(const GL2PSprimitive *prim)
+{
+  int i;
+
+  for(i = 1; i < prim->numverts; i++){
+    if(!gl2psSameColor(prim->verts[0].rgba, prim->verts[i].rgba)){
+      return GL_FALSE;
+    }
+  }
+  return GL_TRUE;
+}
+
+static GLboolean gl2psSameColorThreshold(int n, GL2PSrgba rgba[],
+                                         GL2PSrgba threshold)
+{
+  int i;
+
+  if(n < 2) return GL_TRUE;
+  
+  for(i = 1; i < n; i++){
+    if(fabs(rgba[0][0] - rgba[i][0]) > threshold[0] ||
+       fabs(rgba[0][1] - rgba[i][1]) > threshold[1] ||
+       fabs(rgba[0][2] - rgba[i][2]) > threshold[2])
+      return GL_FALSE;
+  }
+  
+  return GL_TRUE;
+}
+
+static void gl2psSetLastColor(GL2PSrgba rgba)
+{
+  int i;        
+  for(i = 0; i < 3; ++i){
+    gl2ps->lastrgba[i] = rgba[i];
+  }
+}
+
+static GLfloat gl2psGetRGB(GL2PSimage *im, GLuint x, GLuint y,
+                           GLfloat *red, GLfloat *green, GLfloat *blue)
+{
+  
+  GLsizei width = im->width;
+  GLsizei height = im->height;
+  GLfloat *pixels = im->pixels;
+  GLfloat *pimag;
+
+  /* OpenGL image is from down to up, PS image is up to down */  
+  switch(im->format){
+  case GL_RGBA:
+    pimag = pixels + 4 * (width * (height - 1 - y) + x);
+    break;
+  case GL_RGB:
+  default:
+    pimag = pixels + 3 * (width * (height - 1 - y) + x);
+    break;
+  }
+  *red = *pimag; pimag++;
+  *green = *pimag; pimag++;
+  *blue = *pimag; pimag++;
+
+  return (im->format == GL_RGBA) ? *pimag : 1.0F;
+}
+
+/* Helper routines for pixmaps */
+
+static GL2PSimage *gl2psCopyPixmap(GL2PSimage *im)
+{
+  int size;
+  GL2PSimage *image = (GL2PSimage*)gl2psMalloc(sizeof(GL2PSimage));
+  
+  image->width = im->width;
+  image->height = im->height;
+  image->format = im->format;
+  image->type = im->type;
+
+  switch(image->format){
+  case GL_RGBA:
+    size = image->height * image->width * 4 * sizeof(GLfloat);
+    break;
+  case GL_RGB:
+  default:
+    size = image->height * image->width * 3 * sizeof(GLfloat);
+    break;
+  }
+
+  image->pixels = (GLfloat*)gl2psMalloc(size);
+  memcpy(image->pixels, im->pixels, size);
+  
+  return image;
+}
+
+static void gl2psFreePixmap(GL2PSimage *im)
+{
+  if(!im)
+    return;
+  gl2psFree(im->pixels);
+  gl2psFree(im);
+}
+
+#if defined(GL2PS_HAVE_LIBPNG)
+
+#if !defined(png_jmpbuf)
+#  define png_jmpbuf(png_ptr) ((png_ptr)->jmpbuf)
+#endif
+
+static void gl2psUserWritePNG(png_structp png_ptr, png_bytep data, png_size_t length)
+{
+  unsigned int i;
+  GL2PSlist *png = (GL2PSlist*)png_get_io_ptr(png_ptr);
+  for(i = 0; i < length; i++) 
+    gl2psListAdd(png, &data[i]);
+}
+
+static void gl2psUserFlushPNG(png_structp png_ptr)
+{
+}
+
+static void gl2psConvertPixmapToPNG(GL2PSimage *pixmap, GL2PSlist *png)
+{
+  png_structp png_ptr;
+  png_infop info_ptr;
+  unsigned char *row_data;
+  GLfloat dr, dg, db;
+  int row, col;
+
+  if(!(png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL)))
+    return;
+  
+  if(!(info_ptr = png_create_info_struct(png_ptr))){
+    png_destroy_write_struct(&png_ptr, NULL);
+    return;
+  }
+  
+  if(setjmp(png_jmpbuf(png_ptr))) {
+    png_destroy_write_struct(&png_ptr, &info_ptr);
+    return;
+  }
+  
+  png_set_write_fn(png_ptr, (void *)png, gl2psUserWritePNG, gl2psUserFlushPNG);
+  png_set_compression_level(png_ptr, Z_DEFAULT_COMPRESSION);
+  png_set_IHDR(png_ptr, info_ptr, pixmap->width, pixmap->height, 8, 
+               PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, 
+               PNG_FILTER_TYPE_BASE);
+  png_write_info(png_ptr, info_ptr);
+
+  row_data = (unsigned char*)gl2psMalloc(3 * pixmap->width * sizeof(unsigned char));
+  for(row = 0; row < pixmap->height; row++){
+    for(col = 0; col < pixmap->width; col++){
+      gl2psGetRGB(pixmap, col, row, &dr, &dg, &db);
+      row_data[3*col] = (unsigned char)(255. * dr);
+      row_data[3*col+1] = (unsigned char)(255. * dg);
+      row_data[3*col+2] = (unsigned char)(255. * db);
+    }
+    png_write_row(png_ptr, (png_bytep)row_data);
+  }
+  gl2psFree(row_data);
+
+  png_write_end(png_ptr, info_ptr);
+  png_destroy_write_struct(&png_ptr, &info_ptr);
+}
+
+#endif
+
+/* Helper routines for text strings */
+
+static GLint gl2psAddText(GLint type, const char *str, const char *fontname, 
+                          GLshort fontsize, GLint alignment, GLfloat angle)
+{
+  GLfloat pos[4];
+  GL2PSprimitive *prim;
+  GLboolean valid;
+
+  if(!gl2ps || !str || !fontname) return GL2PS_UNINITIALIZED;
+
+  if(gl2ps->options & GL2PS_NO_TEXT) return GL2PS_SUCCESS;
+
+  glGetBooleanv(GL_CURRENT_RASTER_POSITION_VALID, &valid);
+  if(GL_FALSE == valid) return GL2PS_SUCCESS; /* the primitive is culled */
+
+  glGetFloatv(GL_CURRENT_RASTER_POSITION, pos);
+
+  prim = (GL2PSprimitive*)gl2psMalloc(sizeof(GL2PSprimitive));
+  prim->type = type;
+  prim->boundary = 0;
+  prim->numverts = 1;
+  prim->verts = (GL2PSvertex*)gl2psMalloc(sizeof(GL2PSvertex));
+  prim->verts[0].xyz[0] = pos[0];
+  prim->verts[0].xyz[1] = pos[1];
+  prim->verts[0].xyz[2] = pos[2];
+  prim->culled = 0;
+  prim->offset = 0;
+  prim->pattern = 0;
+  prim->factor = 0;
+  prim->width = 1;
+  glGetFloatv(GL_CURRENT_RASTER_COLOR, prim->verts[0].rgba);
+  prim->data.text = (GL2PSstring*)gl2psMalloc(sizeof(GL2PSstring));
+  prim->data.text->str = (char*)gl2psMalloc((strlen(str)+1)*sizeof(char));
+  strcpy(prim->data.text->str, str); 
+  prim->data.text->fontname = (char*)gl2psMalloc((strlen(fontname)+1)*sizeof(char));
+  strcpy(prim->data.text->fontname, fontname);
+  prim->data.text->fontsize = fontsize;
+  prim->data.text->alignment = alignment;
+  prim->data.text->angle = angle;
+
+  gl2psListAdd(gl2ps->auxprimitives, &prim);
+  glPassThrough(GL2PS_TEXT_TOKEN);
+    
+  return GL2PS_SUCCESS;
+}
+
+static GL2PSstring *gl2psCopyText(GL2PSstring *t)
+{
+  GL2PSstring *text = (GL2PSstring*)gl2psMalloc(sizeof(GL2PSstring));
+  text->str = (char*)gl2psMalloc((strlen(t->str)+1)*sizeof(char));
+  strcpy(text->str, t->str); 
+  text->fontname = (char*)gl2psMalloc((strlen(t->fontname)+1)*sizeof(char));
+  strcpy(text->fontname, t->fontname);
+  text->fontsize = t->fontsize;
+  text->alignment = t->alignment;
+  text->angle = t->angle;
+  
+  return text;
+}
+
+static void gl2psFreeText(GL2PSstring *text)
+{
+  if(!text)
+    return;
+  gl2psFree(text->str);
+  gl2psFree(text->fontname);
+  gl2psFree(text);
+}
+
+/* Helpers for blending modes */
+
+static GLboolean gl2psSupportedBlendMode(GLenum sfactor, GLenum dfactor)
+{
+  /* returns TRUE if gl2ps supports the argument combination: only two
+     blending modes have been implemented so far */
+
+  if( (sfactor == GL_SRC_ALPHA && dfactor == GL_ONE_MINUS_SRC_ALPHA) || 
+      (sfactor == GL_ONE && dfactor == GL_ZERO) )
+    return GL_TRUE;
+  return GL_FALSE;
+}
+
+static void gl2psAdaptVertexForBlending(GL2PSvertex *v)
+{
+  /* Transforms vertex depending on the actual blending function -
+     currently the vertex v is considered as source vertex and his
+     alpha value is changed to 1.0 if source blending GL_ONE is
+     active. This might be extended in the future */
+
+  if(!v || !gl2ps)
+    return;
+
+  if(gl2ps->options & GL2PS_NO_BLENDING || !gl2ps->blending){
+    v->rgba[3] = 1.0F;
+    return;
+  }
+  
+  switch(gl2ps->blendfunc[0]){
+  case GL_ONE:
+    v->rgba[3] = 1.0F;
+    break;
+  default:
+    break;
+  }
+}
+
+static void gl2psAssignTriangleProperties(GL2PStriangle *t)
+{
+  /* int i; */
+
+  t->prop = T_VAR_COLOR;
+
+  /* Uncommenting the following lines activates an even more fine
+     grained distinction between triangle types - please don't delete,
+     a remarkable amount of PDF handling code inside this file depends
+     on it if activated */
+  /*
+  t->prop = T_CONST_COLOR;    
+  for(i = 0; i < 3; ++i){
+    if(!GL2PS_ZERO(t->vertex[0].rgba[i] - t->vertex[1].rgba[i]) || 
+       !GL2PS_ZERO(t->vertex[1].rgba[i] - t->vertex[2].rgba[i])){
+      t->prop = T_VAR_COLOR;
+      break;
+    }
+  }
+  */
+
+  if(!GL2PS_ZERO(t->vertex[0].rgba[3] - t->vertex[1].rgba[3]) || 
+     !GL2PS_ZERO(t->vertex[1].rgba[3] - t->vertex[2].rgba[3])){
+    t->prop |= T_VAR_ALPHA;
+  }
+  else{
+    if(t->vertex[0].rgba[3] < 1)
+      t->prop |= T_ALPHA_LESS_1;
+    else
+      t->prop |= T_ALPHA_1;
+  }
+}
+
+static void gl2psFillTriangleFromPrimitive(GL2PStriangle *t, GL2PSprimitive *p,
+                                           GLboolean assignprops)
+{
+  t->vertex[0] = p->verts[0];
+  t->vertex[1] = p->verts[1];
+  t->vertex[2] = p->verts[2];
+  if(GL_TRUE == assignprops)
+    gl2psAssignTriangleProperties(t);
+}
+
+static void gl2psInitTriangle(GL2PStriangle *t)
+{
+  int i;
+  GL2PSvertex vertex = { {-1.0F, -1.0F, -1.0F}, {-1.0F, -1.0F, -1.0F, -1.0F} };
+  for(i = 0; i < 3; i++)
+    t->vertex[i] = vertex;
+  t->prop = T_UNDEFINED;
+}
+
+/* Miscellaneous helper routines */
+
+static GL2PSprimitive *gl2psCopyPrimitive(GL2PSprimitive *p)
+{
+  GL2PSprimitive *prim;
+
+  if(!p){
+    gl2psMsg(GL2PS_ERROR, "Trying to copy an empty primitive");
+    return NULL;
+  }
+
+  prim = (GL2PSprimitive*)gl2psMalloc(sizeof(GL2PSprimitive));
+  
+  prim->type = p->type;
+  prim->numverts = p->numverts;
+  prim->boundary = p->boundary;
+  prim->offset = p->offset;
+  prim->pattern = p->pattern;
+  prim->factor = p->factor;
+  prim->culled = p->culled;
+  prim->width = p->width;
+  prim->verts = (GL2PSvertex*)gl2psMalloc(p->numverts*sizeof(GL2PSvertex));
+  memcpy(prim->verts, p->verts, p->numverts * sizeof(GL2PSvertex));
+
+  switch(prim->type){
+  case GL2PS_PIXMAP :
+    prim->data.image = gl2psCopyPixmap(p->data.image);
+    break;
+  case GL2PS_TEXT :
+  case GL2PS_SPECIAL :
+    prim->data.text = gl2psCopyText(p->data.text);
+    break;
+  default:
+    break;
+  }
+
+  return prim;
+}
+
+static GLboolean gl2psSamePosition(GL2PSxyz p1, GL2PSxyz p2)
+{
+  if(!GL2PS_ZERO(p1[0] - p2[0]) ||
+     !GL2PS_ZERO(p1[1] - p2[1]) ||
+     !GL2PS_ZERO(p1[2] - p2[2]))
+    return GL_FALSE;
+  return GL_TRUE;
+}
+
+/********************************************************************* 
+ *
+ * 3D sorting routines 
+ *
+ *********************************************************************/
+
+static GLfloat gl2psComparePointPlane(GL2PSxyz point, GL2PSplane plane)
+{
+  return (plane[0] * point[0] + 
+          plane[1] * point[1] + 
+          plane[2] * point[2] + 
+          plane[3]);
+}
+
+static GLfloat gl2psPsca(GLfloat *a, GLfloat *b)
+{
+  return (a[0]*b[0] + a[1]*b[1] + a[2]*b[2]);
+}
+
+static void gl2psPvec(GLfloat *a, GLfloat *b, GLfloat *c)
+{
+  c[0] = a[1]*b[2] - a[2]*b[1];
+  c[1] = a[2]*b[0] - a[0]*b[2];
+  c[2] = a[0]*b[1] - a[1]*b[0];
+}
+
+static GLfloat gl2psNorm(GLfloat *a)
+{
+  return (GLfloat)sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
+}
+
+static void gl2psGetNormal(GLfloat *a, GLfloat *b, GLfloat *c)
+{
+  GLfloat norm;
+
+  gl2psPvec(a, b, c);
+  if(!GL2PS_ZERO(norm = gl2psNorm(c))){
+    c[0] = c[0] / norm;
+    c[1] = c[1] / norm;
+    c[2] = c[2] / norm;
+  }
+  else{
+    /* The plane is still wrong despite our tests in gl2psGetPlane.
+       Let's return a dummy value for now (this is a hack: we should
+       do more intelligent tests in GetPlane) */
+    c[0] = c[1] = 0.0F;
+    c[2] = 1.0F;
+  }
+}
+
+static void gl2psGetPlane(GL2PSprimitive *prim, GL2PSplane plane)
+{
+  GL2PSxyz v = {0.0F, 0.0F, 0.0F}, w = {0.0F, 0.0F, 0.0F};
+
+  switch(prim->type){
+  case GL2PS_TRIANGLE :
+  case GL2PS_QUADRANGLE :
+    v[0] = prim->verts[1].xyz[0] - prim->verts[0].xyz[0]; 
+    v[1] = prim->verts[1].xyz[1] - prim->verts[0].xyz[1]; 
+    v[2] = prim->verts[1].xyz[2] - prim->verts[0].xyz[2]; 
+    w[0] = prim->verts[2].xyz[0] - prim->verts[0].xyz[0]; 
+    w[1] = prim->verts[2].xyz[1] - prim->verts[0].xyz[1]; 
+    w[2] = prim->verts[2].xyz[2] - prim->verts[0].xyz[2]; 
+    if((GL2PS_ZERO(v[0]) && GL2PS_ZERO(v[1]) && GL2PS_ZERO(v[2])) || 
+       (GL2PS_ZERO(w[0]) && GL2PS_ZERO(w[1]) && GL2PS_ZERO(w[2]))){
+      plane[0] = plane[1] = 0.0F;
+      plane[2] = 1.0F;
+      plane[3] = -prim->verts[0].xyz[2];
+    }
+    else{
+      gl2psGetNormal(v, w, plane);
+      plane[3] = 
+        - plane[0] * prim->verts[0].xyz[0] 
+        - plane[1] * prim->verts[0].xyz[1] 
+        - plane[2] * prim->verts[0].xyz[2];
+    }
+    break;
+  case GL2PS_LINE :
+    v[0] = prim->verts[1].xyz[0] - prim->verts[0].xyz[0]; 
+    v[1] = prim->verts[1].xyz[1] - prim->verts[0].xyz[1]; 
+    v[2] = prim->verts[1].xyz[2] - prim->verts[0].xyz[2]; 
+    if(GL2PS_ZERO(v[0]) && GL2PS_ZERO(v[1]) && GL2PS_ZERO(v[2])){
+      plane[0] = plane[1] = 0.0F;
+      plane[2] = 1.0F;
+      plane[3] = -prim->verts[0].xyz[2];
+    }
+    else{
+      if(GL2PS_ZERO(v[0]))      w[0] = 1.0F;
+      else if(GL2PS_ZERO(v[1])) w[1] = 1.0F;
+      else                      w[2] = 1.0F;
+      gl2psGetNormal(v, w, plane);
+      plane[3] = 
+        - plane[0] * prim->verts[0].xyz[0] 
+        - plane[1] * prim->verts[0].xyz[1] 
+        - plane[2] * prim->verts[0].xyz[2];
+    }
+    break;
+  case GL2PS_POINT :
+  case GL2PS_PIXMAP :
+  case GL2PS_TEXT :
+  case GL2PS_SPECIAL :
+  case GL2PS_IMAGEMAP:
+    plane[0] = plane[1] = 0.0F;
+    plane[2] = 1.0F;
+    plane[3] = -prim->verts[0].xyz[2];
+    break;
+  default :
+    gl2psMsg(GL2PS_ERROR, "Unknown primitive type in BSP tree");
+    plane[0] = plane[1] = plane[3] = 0.0F;
+    plane[2] = 1.0F;
+    break;
+  }
+}
+
+static void gl2psCutEdge(GL2PSvertex *a, GL2PSvertex *b, GL2PSplane plane,
+                         GL2PSvertex *c)
+{
+  GL2PSxyz v;
+  GLfloat sect, psca;
+
+  v[0] = b->xyz[0] - a->xyz[0];
+  v[1] = b->xyz[1] - a->xyz[1];
+  v[2] = b->xyz[2] - a->xyz[2];
+
+  if(!GL2PS_ZERO(psca = gl2psPsca(plane, v)))
+    sect = -gl2psComparePointPlane(a->xyz, plane) / psca;
+  else
+    sect = 0.0F;
+  
+  c->xyz[0] = a->xyz[0] + v[0] * sect;
+  c->xyz[1] = a->xyz[1] + v[1] * sect;
+  c->xyz[2] = a->xyz[2] + v[2] * sect;
+  
+  c->rgba[0] = (1 - sect) * a->rgba[0] + sect * b->rgba[0];
+  c->rgba[1] = (1 - sect) * a->rgba[1] + sect * b->rgba[1];
+  c->rgba[2] = (1 - sect) * a->rgba[2] + sect * b->rgba[2];
+  c->rgba[3] = (1 - sect) * a->rgba[3] + sect * b->rgba[3];
+}
+
+static void gl2psCreateSplitPrimitive(GL2PSprimitive *parent, GL2PSplane plane,
+                                      GL2PSprimitive *child, GLshort numverts,
+                                      GLshort *index0, GLshort *index1)
+{
+  GLshort i;
+
+  if(parent->type == GL2PS_IMAGEMAP){
+    child->type = GL2PS_IMAGEMAP;
+    child->data.image = parent->data.image;
+  }
+  else{
+    if(numverts > 4){
+      gl2psMsg(GL2PS_WARNING, "%d vertices in polygon", numverts);
+      numverts = 4;
+    }
+    switch(numverts){
+    case 1 : child->type = GL2PS_POINT; break; 
+    case 2 : child->type = GL2PS_LINE; break; 
+    case 3 : child->type = GL2PS_TRIANGLE; break; 
+    case 4 : child->type = GL2PS_QUADRANGLE; break;    
+    default: child->type = GL2PS_NO_TYPE; break;
+    }
+  }
+
+  child->boundary = 0; /* FIXME: not done! */
+  child->culled = parent->culled;
+  child->offset = parent->offset;
+  child->pattern = parent->pattern;
+  child->factor = parent->factor;
+  child->width = parent->width;
+  child->numverts = numverts;
+  child->verts = (GL2PSvertex*)gl2psMalloc(numverts * sizeof(GL2PSvertex));
+
+  for(i = 0; i < numverts; i++){
+    if(index1[i] < 0){
+      child->verts[i] = parent->verts[index0[i]];
+    }
+    else{
+      gl2psCutEdge(&parent->verts[index0[i]], &parent->verts[index1[i]], 
+                   plane, &child->verts[i]);
+    }
+  }
+}
+
+static void gl2psAddIndex(GLshort *index0, GLshort *index1, GLshort *nb, 
+                          GLshort i, GLshort j)
+{
+  GLint k;
+
+  for(k = 0; k < *nb; k++){
+    if((index0[k] == i && index1[k] == j) ||
+       (index1[k] == i && index0[k] == j)) return;
+  }
+  index0[*nb] = i;
+  index1[*nb] = j;
+  (*nb)++;
+}
+
+static GLshort gl2psGetIndex(GLshort i, GLshort num)
+{
+  return (i < num - 1) ? i + 1 : 0;
+}
+
+static GLint gl2psTestSplitPrimitive(GL2PSprimitive *prim, GL2PSplane plane)
+{
+  GLint type = GL2PS_COINCIDENT;
+  GLshort i, j;
+  GLfloat d[5]; 
+
+  for(i = 0; i < prim->numverts; i++){  
+    d[i] = gl2psComparePointPlane(prim->verts[i].xyz, plane);
+  }
+
+  if(prim->numverts < 2){
+    return 0;
+  }
+  else{
+    for(i = 0; i < prim->numverts; i++){
+      j = gl2psGetIndex(i, prim->numverts);
+      if(d[j] > GL2PS_EPSILON){
+        if(type == GL2PS_COINCIDENT)      type = GL2PS_IN_BACK_OF;
+        else if(type != GL2PS_IN_BACK_OF) return 1; 
+        if(d[i] < -GL2PS_EPSILON)         return 1;
+      }
+      else if(d[j] < -GL2PS_EPSILON){
+        if(type == GL2PS_COINCIDENT)       type = GL2PS_IN_FRONT_OF;   
+        else if(type != GL2PS_IN_FRONT_OF) return 1;
+        if(d[i] > GL2PS_EPSILON)           return 1;
+      }
+    }
+  }
+  return 0;
+}
+
+static GLint gl2psSplitPrimitive(GL2PSprimitive *prim, GL2PSplane plane, 
+                                 GL2PSprimitive **front, GL2PSprimitive **back)
+{
+  GLshort i, j, in = 0, out = 0, in0[5], in1[5], out0[5], out1[5];
+  GLint type;
+  GLfloat d[5]; 
+
+  type = GL2PS_COINCIDENT;
+
+  for(i = 0; i < prim->numverts; i++){  
+    d[i] = gl2psComparePointPlane(prim->verts[i].xyz, plane);
+  }
+
+  switch(prim->type){
+  case GL2PS_POINT :
+    if(d[0] > GL2PS_EPSILON)       type = GL2PS_IN_BACK_OF;
+    else if(d[0] < -GL2PS_EPSILON) type = GL2PS_IN_FRONT_OF;
+    else                           type = GL2PS_COINCIDENT;
+    break;
+  default :
+    for(i = 0; i < prim->numverts; i++){
+      j = gl2psGetIndex(i, prim->numverts);
+      if(d[j] > GL2PS_EPSILON){
+        if(type == GL2PS_COINCIDENT)      type = GL2PS_IN_BACK_OF;
+        else if(type != GL2PS_IN_BACK_OF) type = GL2PS_SPANNING; 
+        if(d[i] < -GL2PS_EPSILON){
+          gl2psAddIndex(in0, in1, &in, i, j);
+          gl2psAddIndex(out0, out1, &out, i, j);
+          type = GL2PS_SPANNING;
+        }
+        gl2psAddIndex(out0, out1, &out, j, -1);
+      }
+      else if(d[j] < -GL2PS_EPSILON){
+        if(type == GL2PS_COINCIDENT)       type = GL2PS_IN_FRONT_OF;   
+        else if(type != GL2PS_IN_FRONT_OF) type = GL2PS_SPANNING;
+        if(d[i] > GL2PS_EPSILON){
+          gl2psAddIndex(in0, in1, &in, i, j);
+          gl2psAddIndex(out0, out1, &out, i, j);
+          type = GL2PS_SPANNING;
+        }
+        gl2psAddIndex(in0, in1, &in, j, -1);
+      }
+      else{
+        gl2psAddIndex(in0, in1, &in, j, -1);
+        gl2psAddIndex(out0, out1, &out, j, -1);
+      }
+    }
+    break;
+  }
+
+  if(type == GL2PS_SPANNING){
+    *back = (GL2PSprimitive*)gl2psMalloc(sizeof(GL2PSprimitive));
+    *front = (GL2PSprimitive*)gl2psMalloc(sizeof(GL2PSprimitive));
+    gl2psCreateSplitPrimitive(prim, plane, *back, out, out0, out1);
+    gl2psCreateSplitPrimitive(prim, plane, *front, in, in0, in1);
+  }
+
+  return type;
+}
+
+static void gl2psDivideQuad(GL2PSprimitive *quad, 
+                            GL2PSprimitive **t1, GL2PSprimitive **t2)
+{
+  *t1 = (GL2PSprimitive*)gl2psMalloc(sizeof(GL2PSprimitive));
+  *t2 = (GL2PSprimitive*)gl2psMalloc(sizeof(GL2PSprimitive));
+  (*t1)->type = (*t2)->type = GL2PS_TRIANGLE;
+  (*t1)->numverts = (*t2)->numverts = 3;
+  (*t1)->culled = (*t2)->culled = quad->culled;
+  (*t1)->offset = (*t2)->offset = quad->offset;
+  (*t1)->pattern = (*t2)->pattern = quad->pattern;
+  (*t1)->factor = (*t2)->factor = quad->factor;
+  (*t1)->width = (*t2)->width = quad->width;
+  (*t1)->verts = (GL2PSvertex*)gl2psMalloc(3 * sizeof(GL2PSvertex));
+  (*t2)->verts = (GL2PSvertex*)gl2psMalloc(3 * sizeof(GL2PSvertex));
+  (*t1)->verts[0] = quad->verts[0];
+  (*t1)->verts[1] = quad->verts[1];
+  (*t1)->verts[2] = quad->verts[2];
+  (*t1)->boundary = ((quad->boundary & 1) ? 1 : 0) | ((quad->boundary & 2) ? 2 : 0);
+  (*t2)->verts[0] = quad->verts[0];
+  (*t2)->verts[1] = quad->verts[2];
+  (*t2)->verts[2] = quad->verts[3];
+  (*t2)->boundary = ((quad->boundary & 4) ? 2 : 0) | ((quad->boundary & 4) ? 2 : 0);
+}
+
+static int gl2psCompareDepth(const void *a, const void *b)
+{
+  GL2PSprimitive *q, *w;
+  GLfloat dq = 0.0F, dw = 0.0F, diff;
+  int i;
+  
+  q = *(GL2PSprimitive**)a;
+  w = *(GL2PSprimitive**)b;
+
+  for(i = 0; i < q->numverts; i++){
+    dq += q->verts[i].xyz[2]; 
+  }
+  dq /= (GLfloat)q->numverts;
+
+  for(i = 0; i < w->numverts; i++){
+    dw += w->verts[i].xyz[2]; 
+  }
+  dw /= (GLfloat)w->numverts;
+
+  diff = dq - dw;
+  if(diff > 0.){
+    return -1;
+  }
+  else if(diff < 0.){
+    return 1;
+  }
+  else{
+    return 0;
+  }
+}
+
+static int gl2psTrianglesFirst(const void *a, const void *b)
+{
+  GL2PSprimitive *q, *w;
+
+  q = *(GL2PSprimitive**)a;
+  w = *(GL2PSprimitive**)b;
+  return (q->type < w->type ? 1 : -1);
+}
+
+static GLint gl2psFindRoot(GL2PSlist *primitives, GL2PSprimitive **root)
+{
+  GLint i, j, count, best = 1000000, index = 0;
+  GL2PSprimitive *prim1, *prim2;
+  GL2PSplane plane;
+  GLint maxp;
+
+  if(!gl2psListNbr(primitives)){
+    gl2psMsg(GL2PS_ERROR, "Cannot fint root in empty primitive list");
+    return 0;
+  }
+
+  *root = *(GL2PSprimitive**)gl2psListPointer(primitives, 0);
+
+  if(gl2ps->options & GL2PS_BEST_ROOT){
+    maxp = gl2psListNbr(primitives);
+    if(maxp > gl2ps->maxbestroot){
+      maxp = gl2ps->maxbestroot;
+    }
+    for(i = 0; i < maxp; i++){
+      prim1 = *(GL2PSprimitive**)gl2psListPointer(primitives, i);
+      gl2psGetPlane(prim1, plane);
+      count = 0;
+      for(j = 0; j < gl2psListNbr(primitives); j++){
+        if(j != i){
+          prim2 = *(GL2PSprimitive**)gl2psListPointer(primitives, j);
+          count += gl2psTestSplitPrimitive(prim2, plane); 
+        }
+        if(count > best) break;
+      }
+      if(count < best){
+        best = count;
+        index = i;
+        *root = prim1;
+        if(!count) return index;
+      }
+    }
+    /* if(index) gl2psMsg(GL2PS_INFO, "GL2PS_BEST_ROOT was worth it: %d", index); */
+    return index;
+  }
+  else{
+    return 0;
+  }
+}
+
+static void gl2psFreeImagemap(GL2PSimagemap *list)
+{
+  GL2PSimagemap *next;
+  while(list != NULL){
+    next = list->next;
+    gl2psFree(list->image->pixels);
+    gl2psFree(list->image);
+    gl2psFree(list);
+    list = next;
+  }
+}
+
+static void gl2psFreePrimitive(void *data)
+{
+  GL2PSprimitive *q;
+  
+  q = *(GL2PSprimitive**)data;
+  gl2psFree(q->verts);
+  if(q->type == GL2PS_TEXT || q->type == GL2PS_SPECIAL){
+    gl2psFreeText(q->data.text);
+  }
+  else if(q->type == GL2PS_PIXMAP){
+    gl2psFreePixmap(q->data.image);
+  }
+  gl2psFree(q);
+}
+
+static void gl2psAddPrimitiveInList(GL2PSprimitive *prim, GL2PSlist *list)
+{
+  GL2PSprimitive *t1, *t2;
+
+  if(prim->type != GL2PS_QUADRANGLE){
+    gl2psListAdd(list, &prim);
+  }
+  else{
+    gl2psDivideQuad(prim, &t1, &t2);
+    gl2psListAdd(list, &t1);
+    gl2psListAdd(list, &t2);
+    gl2psFreePrimitive(&prim);
+  }
+  
+}
+
+static void gl2psFreeBspTree(GL2PSbsptree **tree)
+{
+  if(*tree){
+    if((*tree)->back) gl2psFreeBspTree(&(*tree)->back);
+    if((*tree)->primitives){
+      gl2psListAction((*tree)->primitives, gl2psFreePrimitive);
+      gl2psListDelete((*tree)->primitives);
+    }
+    if((*tree)->front) gl2psFreeBspTree(&(*tree)->front);
+    gl2psFree(*tree);
+    *tree = NULL;
+  }
+}
+
+static GLboolean gl2psGreater(GLfloat f1, GLfloat f2)
+{
+  if(f1 > f2) return GL_TRUE;
+  else return GL_FALSE;
+}
+
+static GLboolean gl2psLess(GLfloat f1, GLfloat f2)
+{
+  if(f1 < f2) return GL_TRUE;
+  else return GL_FALSE;
+}
+
+static void gl2psBuildBspTree(GL2PSbsptree *tree, GL2PSlist *primitives)
+{
+  GL2PSprimitive *prim, *frontprim = NULL, *backprim = NULL;
+  GL2PSlist *frontlist, *backlist;
+  GLint i, index;
+
+  tree->front = NULL;
+  tree->back = NULL;
+  tree->primitives = gl2psListCreate(1, 2, sizeof(GL2PSprimitive*));
+  index = gl2psFindRoot(primitives, &prim);
+  gl2psGetPlane(prim, tree->plane);
+  gl2psAddPrimitiveInList(prim, tree->primitives);
+
+  frontlist = gl2psListCreate(1, 2, sizeof(GL2PSprimitive*));
+  backlist = gl2psListCreate(1, 2, sizeof(GL2PSprimitive*));
+
+  for(i = 0; i < gl2psListNbr(primitives); i++){
+    if(i != index){
+      prim = *(GL2PSprimitive**)gl2psListPointer(primitives,i);
+      switch(gl2psSplitPrimitive(prim, tree->plane, &frontprim, &backprim)){
+      case GL2PS_COINCIDENT:
+        gl2psAddPrimitiveInList(prim, tree->primitives);
+        break;
+      case GL2PS_IN_BACK_OF:
+        gl2psAddPrimitiveInList(prim, backlist);
+        break;
+      case GL2PS_IN_FRONT_OF:
+        gl2psAddPrimitiveInList(prim, frontlist);
+        break;
+      case GL2PS_SPANNING:
+        gl2psAddPrimitiveInList(backprim, backlist);
+        gl2psAddPrimitiveInList(frontprim, frontlist);
+        gl2psFreePrimitive(&prim);
+        break;
+      }
+    }
+  }
+
+  if(gl2psListNbr(tree->primitives)){
+    gl2psListSort(tree->primitives, gl2psTrianglesFirst);
+  }
+
+  if(gl2psListNbr(frontlist)){
+    gl2psListSort(frontlist, gl2psTrianglesFirst);
+    tree->front = (GL2PSbsptree*)gl2psMalloc(sizeof(GL2PSbsptree));
+    gl2psBuildBspTree(tree->front, frontlist);
+  }
+  else{
+    gl2psListDelete(frontlist);
+  }
+
+  if(gl2psListNbr(backlist)){
+    gl2psListSort(backlist, gl2psTrianglesFirst);
+    tree->back = (GL2PSbsptree*)gl2psMalloc(sizeof(GL2PSbsptree));
+    gl2psBuildBspTree(tree->back, backlist);
+  }
+  else{
+    gl2psListDelete(backlist);
+  }
+