Application: gvSIG desktop

//

// DelaunayClarkson.java

//

/*

VisAD system for interactive analysis and visualization of numerical

data.  Copyright (C) 1996 - 2002 Bill Hibbard, Curtis Rueden, Tom

Rink, Dave Glowacki, Steve Emmerson, Tom Whittaker, Don Murray, and

Tommy Jasmin.

This library is free software; you can redistribute it and/or

modify it under the terms of 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.

This library 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

Library General Public License for more details.

You should have received a copy of the GNU Library General Public

License along with this library; if not, write to the Free

Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,

MA 02111-1307, USA

*/

package org.gvsig.fmap.algorithm.triangulation.visad;

/* The Delaunay triangulation algorithm in this class

 * is originally from hull by Ken Clarkson:

 *

 * Ken Clarkson wrote this.  Copyright (c) 1995 by AT&T..

 * Permission to use, copy, modify, and distribute this software for any

 * purpose without fee is hereby granted, provided that this entire notice

 * is included in all copies of any software which is or includes a copy

 * or modification of this software and in all copies of the supporting

 * documentation for such software.

 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED

 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY

 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY

 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.

 */

/**

   DelaunayClarkson represents an O(N*logN) method

   with high overhead to find the Delaunay triangulation

   of a set of samples of R^DomainDimension.<P>

*/

public class DelaunayClarkson extends Delaunay {

/* ******* BEGINNING OF CONVERTED HULL CODE ******* */

  // <<<< Constants >>>>

  private static final double DBL_MANT_DIG = 53;

  private static final double FLT_RADIX = 2;

  private static final double DBL_EPSILON = 2.2204460492503131E-16;

  private static final double ln2 = Math.log(2);

  // <<<< Variables >>>>

  /* we need to have two indices for every pointer into basis_s and

     simplex arrays, because they are two-dimensional arrays of

     blocks.  ( _bn = block number ) */

  // for the pseudo-pointers

  private static final int INFINITY = -2;      // replaces infinity

  private static final int NOVAL = -1;         // replaces null

  private double[][] site_blocks;        // copy of samples array

  private int[][]   a3s;                // output array

  private int       a3size;             // output array maximum size

  private int       nts = 0;            // # output objects

  private static final int max_blocks = 10000; // max # basis/simplex blocks

  private static final int Nobj = 10000;

  private static final int MAXDIM = 8;         // max dimension

  private int    dim;

  private int    p;

  private long   pnum;

  private int    rdim,          // region dimension

                 cdim;          // # sites currently specifying region

  private int    exact_bits;

  private double b_err_min,

                 b_err_min_sq;

  private double ldetbound = 0;

  private int    failcount = 0;          // static: reduce_inner

  private int    lscale;                 // static: reduce_inner

  private double max_scale;              // static: reduce_inner

  private float  Sb;                     // static: reduce_inner

  private int    nsb = 0;                // # simplex blocks

  private int    nbb = 0;                // # basis_s blocks

  private int    ss = MAXDIM;            // static: search

  private int    ss2 = 2000;             // static: visit_triang

  private long   vnum = -1;              // static: visit_triang

  private int    p_neigh_vert = NOVAL;   // static: main

  // "void stuff" -- dummy variables to hold unused return information

  private int[] voidp = new int[1];

  private int[] voidp_bn = new int[1];

  // basis_s stuff

  private int[][]      bbt_next = new int[max_blocks][];

  private int[][]      bbt_next_bn = new int[max_blocks][];

  private int[][]      bbt_ref_count = new int[max_blocks][];

  private int[][]      bbt_lscale = new int[max_blocks][];

  private double[][]   bbt_sqa = new double[max_blocks][];

  private double[][]   bbt_sqb = new double[max_blocks][];

  private double[][][] bbt_vecs = new double[max_blocks][][];

  private int ttbp;

  private int ttbp_bn;

  private int ib;

  private int ib_bn;

  private int basis_s_list = NOVAL;

  private int basis_s_list_bn;

  private int pnb = NOVAL;

  private int pnb_bn;

  private int b = NOVAL;              // static: sees

  private int b_bn;

  // simplex stuff

  private int[][]   sbt_next = new int[max_blocks][];

  private int[][]   sbt_next_bn = new int[max_blocks][];

  private long[][]  sbt_visit = new long[max_blocks][];

  private short[][] sbt_mark = new short[max_blocks][];

  private int[][]   sbt_normal = new int[max_blocks][];

  private int[][]   sbt_normal_bn = new int[max_blocks][];

  private int[][]   sbt_peak_vert = new int[max_blocks][];

  private int[][]   sbt_peak_simp = new int[max_blocks][];

  private int[][]   sbt_peak_simp_bn = new int[max_blocks][];

  private int[][]   sbt_peak_basis = new int[max_blocks][];

  private int[][]   sbt_peak_basis_bn = new int[max_blocks][];

  private int[][][] sbt_neigh_vert = new int[max_blocks][][];

  private int[][][] sbt_neigh_simp = new int[max_blocks][][];

  private int[][][] sbt_neigh_simp_bn = new int[max_blocks][][];

  private int[][][] sbt_neigh_basis = new int[max_blocks][][];

  private int[][][] sbt_neigh_basis_bn = new int[max_blocks][][];

  private int   simplex_list = NOVAL;

  private int   simplex_list_bn;

  private int   ch_root;

  private int   ch_root_bn;

  private int   ns;                            // static: make_facets

  private int   ns_bn;

  private int[] st = new int[ss+MAXDIM+1];    // static: search

  private int[] st_bn = new int[ss+MAXDIM+1];

  private int[] st2 = new int[ss2+MAXDIM+1];    // static: visit_triang

  private int[] st2_bn = new int[ss2+MAXDIM+1];

  // <<<< Functions >>>>

  private int new_block_basis_s() {

    bbt_next[nbb] = new int[Nobj];

    bbt_next_bn[nbb] = new int[Nobj];

    bbt_ref_count[nbb] = new int[Nobj];

    bbt_lscale[nbb] = new int[Nobj];

    bbt_sqa[nbb] = new double[Nobj];

    bbt_sqb[nbb] = new double[Nobj];

    bbt_vecs[nbb] = new double[2*rdim][];

    for (int i=0; i<2*rdim; i++) bbt_vecs[nbb][i] = new double[Nobj];

    for (int i=0; i<Nobj; i++) {

      bbt_next[nbb][i] = i+1;

      bbt_next_bn[nbb][i] = nbb;

      bbt_ref_count[nbb][i] = 0;

      bbt_lscale[nbb][i] = 0;

      bbt_sqa[nbb][i] = 0;

      bbt_sqb[nbb][i] = 0;

      for (int j=0; j<2*rdim; j++) bbt_vecs[nbb][j][i] = 0;

    }

    bbt_next[nbb][Nobj-1] = NOVAL;

    basis_s_list = 0;

    basis_s_list_bn = nbb;

    nbb++;

    return basis_s_list;

  }

  private int reduce_inner(int v, int v_bn, int s, int s_bn, int k) {

    int q, q_bn;

    double dd,

           Sb = 0;

    double scale;

    bbt_sqa[v_bn][v] = 0;

    for (int i=0; i<rdim; i++) {

      bbt_sqa[v_bn][v] += bbt_vecs[v_bn][i][v] * bbt_vecs[v_bn][i][v];

    }

    bbt_sqb[v_bn][v] = bbt_sqa[v_bn][v];

    if (k <= 1) {

      for (int i=0; i<rdim; i++) {

        bbt_vecs[v_bn][i][v] = bbt_vecs[v_bn][rdim+i][v];

      }

      return 1;

    }

    for (int j=0; j<250; j++) {

      int    xx = rdim;

      double labound;

      for (int i=0; i<rdim; i++) {

        bbt_vecs[v_bn][i][v] = bbt_vecs[v_bn][rdim+i][v];

      }

      for (int i=k-1; i>0; i--) {

        q = sbt_neigh_basis[s_bn][i][s];

        q_bn = sbt_neigh_basis_bn[s_bn][i][s];

        dd = 0;

        for (int l=0; l<rdim; l++) {

          dd -= bbt_vecs[q_bn][l][q] * bbt_vecs[v_bn][l][v];

        }

        dd /= bbt_sqb[q_bn][q];

        for (int l=0; l<rdim; l++) {

          bbt_vecs[v_bn][l][v] += dd * bbt_vecs[q_bn][rdim+l][q];

        }

      }

      bbt_sqb[v_bn][v] = 0;

      for (int i=0; i<rdim; i++) {

        bbt_sqb[v_bn][v] += bbt_vecs[v_bn][i][v] * bbt_vecs[v_bn][i][v];

      }

      bbt_sqa[v_bn][v] = 0;

      for (int i=0; i<rdim; i++) {

        bbt_sqa[v_bn][v] += bbt_vecs[v_bn][rdim+i][v]

                          * bbt_vecs[v_bn][rdim+i][v];

      }

      if (2*bbt_sqb[v_bn][v] >= bbt_sqa[v_bn][v]) return 1;

      // scale up vector

      if (j < 10) {

        labound = Math.floor(Math.log(bbt_sqa[v_bn][v])/ln2) / 2;

        max_scale = exact_bits-labound-0.66*(k-2)-1;

        if (max_scale < 1) max_scale = 1;

        if (j == 0) {

          ldetbound = 0;

          Sb = 0;

          for (int l=k-1; l>0; l--) {

            q = sbt_neigh_basis[s_bn][l][s];

            q_bn = sbt_neigh_basis_bn[s_bn][l][s];

            Sb += bbt_sqb[q_bn][q];

            ldetbound += Math.floor(Math.log(bbt_sqb[q_bn][q])/ln2) / 2 + 1;

            ldetbound -= bbt_lscale[q_bn][q];

          }

        }

      }

      if (ldetbound - bbt_lscale[v_bn][v]

        + Math.floor(Math.log(bbt_sqb[v_bn][v])/ln2) / 2 + 1 < 0) {

        scale = 0;

      }

      else {

        lscale = (int) (Math.log(2*Sb/(bbt_sqb[v_bn][v]

                                     + bbt_sqa[v_bn][v]*b_err_min))/ln2) / 2;

        if (lscale > max_scale) lscale = (int) max_scale;

        else if (lscale < 0) lscale = 0;

        bbt_lscale[v_bn][v] += lscale;

        scale = (lscale < 20) ? 1 << lscale : Math.pow(2, lscale);

      }

      while (xx < 2*rdim) bbt_vecs[v_bn][xx++][v] *= scale;

      for (int i=k-1; i>0; i--) {

        q = sbt_neigh_basis[s_bn][i][s];

        q_bn = sbt_neigh_basis_bn[s_bn][i][s];

        dd = 0;

        for (int l=0; l<rdim; l++) {

          dd -= bbt_vecs[q_bn][l][q] * bbt_vecs[v_bn][rdim+l][v];

        }

        dd /= bbt_sqb[q_bn][q];

        dd = Math.floor(dd+0.5);

        for (int l=0; l<rdim; l++) {

          bbt_vecs[v_bn][rdim+l][v] += dd * bbt_vecs[q_bn][rdim+l][q];

        }

      }

    }

    if (failcount++ < 10) System.out.println("reduce_inner failed!");

    return 0;

  }

  private int reduce(int[] v, int[] v_bn, int rp, int s, int s_bn, int k) {

    if (v[0] == NOVAL) {

      v[0] = basis_s_list != NOVAL ? basis_s_list : new_block_basis_s();

      v_bn[0] = basis_s_list_bn;

      basis_s_list = bbt_next[v_bn[0]][v[0]];

      basis_s_list_bn = bbt_next_bn[v_bn[0]][v[0]];

      bbt_ref_count[v_bn[0]][v[0]] = 1;

    }

    else bbt_lscale[v_bn[0]][v[0]] = 0;

    if (rp == INFINITY) {

      bbt_next[v_bn[0]][v[0]] = bbt_next[ib_bn][ib];

      bbt_next_bn[v_bn[0]][v[0]] = bbt_next_bn[ib_bn][ib];

      bbt_ref_count[v_bn[0]][v[0]] = bbt_ref_count[ib_bn][ib];

      bbt_lscale[v_bn[0]][v[0]] = bbt_lscale[ib_bn][ib];

      bbt_sqa[v_bn[0]][v[0]] = bbt_sqa[ib_bn][ib];

      bbt_sqb[v_bn[0]][v[0]] = bbt_sqb[ib_bn][ib];

      for (int i=0; i<2*rdim; i++) {

        bbt_vecs[v_bn[0]][i][v[0]] = bbt_vecs[ib_bn][i][ib];

      }

    }

    else {

      double sum = 0;

      int sbt_nv = sbt_neigh_vert[s_bn][0][s];

      if (sbt_nv == INFINITY) {

        for (int i=0; i<dim; i++) {

          bbt_vecs[v_bn[0]][i+rdim][v[0]] = bbt_vecs[v_bn[0]][i][v[0]]

            = (double) site_blocks[i][rp];

        }

      }

      else {

        for (int i=0; i<dim; i++) {

          bbt_vecs[v_bn[0]][i+rdim][v[0]] = bbt_vecs[v_bn[0]][i][v[0]]

            = (double) (site_blocks[i][rp] - site_blocks[i][sbt_nv]);

        }

      }

      for (int i=0; i<dim; i++) {

        sum += bbt_vecs[v_bn[0]][i][v[0]] * bbt_vecs[v_bn[0]][i][v[0]];

      }

      bbt_vecs[v_bn[0]][2*rdim-1][v[0]] = sum;

      bbt_vecs[v_bn[0]][rdim-1][v[0]] = sum;

    }

    return reduce_inner(v[0], v_bn[0], s, s_bn, k);

  }

  private void get_basis_sede(int s, int s_bn) {

    int   k=1;

    int   q, q_bn;

    int[] curt = new int[1];

    int[] curt_bn = new int[1];

    if (sbt_neigh_vert[s_bn][0][s] == INFINITY && cdim > 1) {

      int t_vert, t_simp, t_simp_bn, t_basis, t_basis_bn;

      t_vert = sbt_neigh_vert[s_bn][0][s];

      t_simp = sbt_neigh_simp[s_bn][0][s];

      t_simp_bn = sbt_neigh_simp_bn[s_bn][0][s];

      t_basis = sbt_neigh_basis[s_bn][0][s];

      t_basis_bn = sbt_neigh_basis_bn[s_bn][0][s];

      sbt_neigh_vert[s_bn][0][s] = sbt_neigh_vert[s_bn][k][s];

      sbt_neigh_simp[s_bn][0][s] = sbt_neigh_simp[s_bn][k][s];

      sbt_neigh_simp_bn[s_bn][0][s] = sbt_neigh_simp_bn[s_bn][k][s];

      sbt_neigh_basis[s_bn][0][s] = sbt_neigh_basis[s_bn][k][s];

      sbt_neigh_basis_bn[s_bn][0][s] = sbt_neigh_basis_bn[s_bn][k][s];

      sbt_neigh_vert[s_bn][k][s] = t_vert;

      sbt_neigh_simp[s_bn][k][s] = t_simp;

      sbt_neigh_simp_bn[s_bn][k][s] = t_simp_bn;

      sbt_neigh_basis[s_bn][k][s] = t_basis;

      sbt_neigh_basis_bn[s_bn][k][s] = t_basis_bn;

      q = sbt_neigh_basis[s_bn][0][s];

      q_bn = sbt_neigh_basis_bn[s_bn][0][s];

      if ((q != NOVAL) && --bbt_ref_count[q_bn][q] == 0) {

        bbt_next[q_bn][q] = basis_s_list;

        bbt_next_bn[q_bn][q] = basis_s_list_bn;

        bbt_ref_count[q_bn][q] = 0;

        bbt_lscale[q_bn][q] = 0;

        bbt_sqa[q_bn][q] = 0;

        bbt_sqb[q_bn][q] = 0;

        for (int j=0; j<2*rdim; j++) bbt_vecs[q_bn][j][q] = 0;

        basis_s_list = q;

        basis_s_list_bn = q_bn;

      }

      sbt_neigh_basis[s_bn][0][s] = ttbp;

      sbt_neigh_basis_bn[s_bn][0][s] = ttbp_bn;

      bbt_ref_count[ttbp_bn][ttbp]++;

    }

    else {

      if (sbt_neigh_basis[s_bn][0][s] == NOVAL) {

        sbt_neigh_basis[s_bn][0][s] = ttbp;

        sbt_neigh_basis_bn[s_bn][0][s] = ttbp_bn;

        bbt_ref_count[ttbp_bn][ttbp]++;

      } else while (k < cdim && sbt_neigh_basis[s_bn][k][s] != NOVAL) k++;

    }

    while (k < cdim) {

      q = sbt_neigh_basis[s_bn][k][s];

      q_bn = sbt_neigh_basis_bn[s_bn][k][s];

      if (q != NOVAL && --bbt_ref_count[q_bn][q] == 0) {

        bbt_next[q_bn][q] = basis_s_list;

        bbt_next_bn[q_bn][q] = basis_s_list_bn;

        bbt_ref_count[q_bn][q] = 0;

        bbt_lscale[q_bn][q] = 0;

        bbt_sqa[q_bn][q] = 0;

        bbt_sqb[q_bn][q] = 0;

        for (int j=0; j<2*rdim; j++) bbt_vecs[q_bn][j][q] = 0;

        basis_s_list = q;

        basis_s_list_bn = q_bn;

      }

      sbt_neigh_basis[s_bn][k][s] = NOVAL;

      curt[0] = sbt_neigh_basis[s_bn][k][s];

      curt_bn[0] = sbt_neigh_basis_bn[s_bn][k][s];

      reduce(curt, curt_bn, sbt_neigh_vert[s_bn][k][s], s, s_bn, k);

      sbt_neigh_basis[s_bn][k][s] = curt[0];

      sbt_neigh_basis_bn[s_bn][k][s] = curt_bn[0];

      k++;

    }

  }

  private int sees(int rp, int s, int s_bn) {

    double  dd, dds;

    int     q, q_bn, q1, q1_bn, q2, q2_bn;

    int[]   curt = new int[1];

    int[]   curt_bn = new int[1];

    if (b == NOVAL) {

      b = (basis_s_list != NOVAL) ? basis_s_list : new_block_basis_s();

      b_bn = basis_s_list_bn;

      basis_s_list = bbt_next[b_bn][b];

      basis_s_list_bn = bbt_next_bn[b_bn][b];

    }

    else bbt_lscale[b_bn][b] = 0;

    if (cdim==0) return 0;

    if (sbt_normal[s_bn][s] == NOVAL) {

      get_basis_sede(s, s_bn);

      if (rdim==3 && cdim==3) {

        sbt_normal[s_bn][s] = basis_s_list != NOVAL ? basis_s_list

                                                    : new_block_basis_s();

        sbt_normal_bn[s_bn][s] = basis_s_list_bn;

        q = sbt_normal[s_bn][s];

        q_bn = sbt_normal_bn[s_bn][s];

        basis_s_list = bbt_next[q_bn][q];

        basis_s_list_bn = bbt_next_bn[q_bn][q];

        q1 = sbt_neigh_basis[s_bn][1][s];

        q1_bn = sbt_neigh_basis_bn[s_bn][1][s];

        q2 = sbt_neigh_basis[s_bn][2][s];

        q2_bn = sbt_neigh_basis_bn[s_bn][2][s];

        bbt_ref_count[q_bn][q] = 1;

        bbt_vecs[q_bn][0][q] = bbt_vecs[q1_bn][1][q1]

                  *bbt_vecs[q2_bn][2][q2]

             - bbt_vecs[q1_bn][2][q1]

                  *bbt_vecs[q2_bn][1][q2];

        bbt_vecs[q_bn][1][q] = bbt_vecs[q1_bn][2][q1]

                  *bbt_vecs[q2_bn][0][q2]

             - bbt_vecs[q1_bn][0][q1]

                  *bbt_vecs[q2_bn][2][q2];

        bbt_vecs[q_bn][2][q] = bbt_vecs[q1_bn][0][q1]

                  *bbt_vecs[q2_bn][1][q2]

             - bbt_vecs[q1_bn][1][q1]

                  *bbt_vecs[q2_bn][0][q2];

        bbt_sqb[q_bn][q] = 0;

        for (int i=0; i<rdim; i++) bbt_sqb[q_bn][q] += bbt_vecs[q_bn][i][q]

                                                     * bbt_vecs[q_bn][i][q];

        for (int i=cdim+1; i>0; i--) {

          int m = (i > 1) ? sbt_neigh_vert[ch_root_bn][i-2][ch_root]

                          : INFINITY;

          int j;

          for (j=0; j<cdim && m != sbt_neigh_vert[s_bn][j][s]; j++);

          if (j < cdim) continue;

          if (m == INFINITY) {

            if (bbt_vecs[q_bn][2][q] > -b_err_min) continue;

          }

          else {

            if (sees(m, s, s_bn) == 0) {

              continue;

            }

          }

          bbt_vecs[q_bn][0][q] = -bbt_vecs[q_bn][0][q];

          bbt_vecs[q_bn][1][q] = -bbt_vecs[q_bn][1][q];

          bbt_vecs[q_bn][2][q] = -bbt_vecs[q_bn][2][q];

          break;

        }

      }

      else {

        for (int i=cdim+1; i>0; i--) {

          int m = (i > 1) ? sbt_neigh_vert[ch_root_bn][i-2][ch_root]

                          : INFINITY;

          int j;

          for (j=0; j<cdim && m != sbt_neigh_vert[s_bn][j][s]; j++);

          if (j < cdim) continue;

          curt[0] = sbt_normal[s_bn][s];

          curt_bn[0] = sbt_normal_bn[s_bn][s];

          reduce(curt, curt_bn, m, s, s_bn, cdim);

          q = sbt_normal[s_bn][s] = curt[0];

          q_bn = sbt_normal_bn[s_bn][s] = curt_bn[0];

          if (bbt_sqb[q_bn][q] != 0) break;

        }

      }

      for (int i=0; i<cdim; i++) {

        q = sbt_neigh_basis[s_bn][i][s];

        q_bn = sbt_neigh_basis_bn[s_bn][i][s];

        if (q != NOVAL && --bbt_ref_count[q_bn][q] == 0) {

          bbt_next[q_bn][q] = basis_s_list;

          bbt_next_bn[q_bn][q] = basis_s_list_bn;

          bbt_ref_count[q_bn][q] = 0;

          bbt_lscale[q_bn][q] = 0;

          bbt_sqa[q_bn][q] = 0;

          bbt_sqb[q_bn][q] = 0;

          for (int l=0; l<2*rdim; l++) bbt_vecs[q_bn][l][q] = 0;

          basis_s_list = q;

          basis_s_list_bn = q_bn;

        }

        sbt_neigh_basis[s_bn][i][s] = NOVAL;

      }

    }

    if (rp == INFINITY) {

      bbt_next[b_bn][b] = bbt_next[ib_bn][ib];

      bbt_next_bn[b_bn][b] = bbt_next_bn[ib_bn][ib];

      bbt_ref_count[b_bn][b] = bbt_ref_count[ib_bn][ib];

      bbt_lscale[b_bn][b] = bbt_lscale[ib_bn][ib];

      bbt_sqa[b_bn][b] = bbt_sqa[ib_bn][ib];

      bbt_sqb[b_bn][b] = bbt_sqb[ib_bn][ib];

      for (int i=0; i<2*rdim; i++) {

        bbt_vecs[b_bn][i][b] = bbt_vecs[ib_bn][i][ib];

      }

    }

    else {

      double sum = 0;

      int sbt_nv = sbt_neigh_vert[s_bn][0][s];

      if (sbt_nv == INFINITY) {

        for (int l=0; l<dim; l++) {

          bbt_vecs[b_bn][l+rdim][b] = bbt_vecs[b_bn][l][b]

            = (double) site_blocks[l][rp];

        }

      }

      else {

        for (int l=0; l<dim; l++) {

          bbt_vecs[b_bn][l+rdim][b] = bbt_vecs[b_bn][l][b]

            = (double) (site_blocks[l][rp] - site_blocks[l][sbt_nv]);

        }

      }

      for (int l=0; l<dim; l++) {

        sum += bbt_vecs[b_bn][l][b] * bbt_vecs[b_bn][l][b];

      }

      bbt_vecs[b_bn][2*rdim-1][b] = bbt_vecs[b_bn][rdim-1][b] = sum;

    }

    q = sbt_normal[s_bn][s];

    q_bn = sbt_normal_bn[s_bn][s];

    for (int i=0; i<3; i++) {

      double sum = 0;

      dd = 0;

      for (int l=0; l<rdim; l++) {

        dd += bbt_vecs[b_bn][l][b] * bbt_vecs[q_bn][l][q];

      }

      if (dd == 0.0) return 0;

      for (int l=0; l<rdim; l++) {

        sum += bbt_vecs[b_bn][l][b] * bbt_vecs[b_bn][l][b];

      }

      dds = dd*dd/bbt_sqb[q_bn][q]/sum;

      if (dds > b_err_min_sq) return (dd < 0 ? 1 : 0);

      get_basis_sede(s, s_bn);

      reduce_inner(b, b_bn, s, s_bn, cdim);

    }

    return 0;

  }

  private int new_block_simplex() {

    sbt_next[nsb] = new int[Nobj];

    sbt_next_bn[nsb] = new int[Nobj];

    sbt_visit[nsb] = new long[Nobj];

    sbt_mark[nsb] = new short[Nobj];

    sbt_normal[nsb] = new int[Nobj];

    sbt_normal_bn[nsb] = new int[Nobj];

    sbt_peak_vert[nsb] = new int[Nobj];

    sbt_peak_simp[nsb] = new int[Nobj];

    sbt_peak_simp_bn[nsb] = new int[Nobj];

    sbt_peak_basis[nsb] = new int[Nobj];

    sbt_peak_basis_bn[nsb] = new int[Nobj];

    sbt_neigh_vert[nsb] = new int[rdim][];

    sbt_neigh_simp[nsb] = new int[rdim][];

    sbt_neigh_simp_bn[nsb] = new int[rdim][];

    sbt_neigh_basis[nsb] = new int[rdim][];

    sbt_neigh_basis_bn[nsb] = new int[rdim][];

    for (int i=0; i<rdim; i++) {

      sbt_neigh_vert[nsb][i] = new int[Nobj];

      sbt_neigh_simp[nsb][i] = new int[Nobj];

      sbt_neigh_simp_bn[nsb][i] = new int[Nobj];

      sbt_neigh_basis[nsb][i] = new int[Nobj];

      sbt_neigh_basis_bn[nsb][i] = new int[Nobj];

    }

    for (int i=0; i<Nobj; i++) {

      sbt_next[nsb][i] = i+1;

      sbt_next_bn[nsb][i] = nsb;

      sbt_visit[nsb][i] = 0;

      sbt_mark[nsb][i] = 0;

      sbt_normal[nsb][i] = NOVAL;

      sbt_peak_vert[nsb][i] = NOVAL;

      sbt_peak_simp[nsb][i] = NOVAL;

      sbt_peak_basis[nsb][i] = NOVAL;

      for (int j=0; j<rdim; j++) {

        sbt_neigh_vert[nsb][j][i] = NOVAL;

        sbt_neigh_simp[nsb][j][i] = NOVAL;

        sbt_neigh_basis[nsb][j][i] = NOVAL;

      }

    }

    sbt_next[nsb][Nobj-1] = NOVAL;

    simplex_list = 0;

    simplex_list_bn = nsb;

    nsb++;

    return simplex_list;

  }

  /**

     starting at s, visit simplices t such that test(s,i,0) is true,

     and t is the i'th neighbor of s;

     apply visit function to all visited simplices;

     when visit returns nonnull, exit and return its value.

  */

  private void visit_triang_gen(int s, int s_bn, int whichfunc,

                                int[] ret, int[] ret_bn) {

    int v;

    int v_bn;

    int t;

    int t_bn;

    int tms = 0;

    vnum--;

    if (s != NOVAL) {

      st2[tms] = s;

      st2_bn[tms] = s_bn;

      tms++;

    }

    while (tms != 0) {

      if (tms > ss2) {

        // JAVA: efficiency issue: how much is this stack hammered?

        ss2 += ss2;

        int[] newst2 = new int[ss2+MAXDIM+1];

        int[] newst2_bn = new int[ss2+MAXDIM+1];

        System.arraycopy(st2, 0, newst2, 0, st2.length);

        System.arraycopy(st2_bn, 0, newst2_bn, 0, st2_bn.length);

        st2 = newst2;

        st2_bn = newst2_bn;

      }

      tms--;

      t = st2[tms];

      t_bn = st2_bn[tms];

      if (t == NOVAL || sbt_visit[t_bn][t] == vnum) continue;

      sbt_visit[t_bn][t] = vnum;

      if (whichfunc == 1) {

        if (sbt_peak_vert[t_bn][t] == NOVAL) {

          v = t;

          v_bn = t_bn;

        }

        else {

          v = NOVAL;

          v_bn = NOVAL;

        }

        if (v != NOVAL) {

          ret[0] = v;

          ret_bn[0] = v_bn;

          return;

        }

      }

      else {

        int[] vfp = new int[cdim];

        if (t != NOVAL) {

          for (int j=0; j<cdim; j++) vfp[j] = sbt_neigh_vert[t_bn][j][t];

          for (int j=0; j<cdim; j++) {

            a3s[j][nts] = (vfp[j] == INFINITY) ? -1 : vfp[j];

          }

          nts++;

          if (nts > a3size) {

            // JAVA: efficiency issue, hammering an array

            a3size += a3size;

            int[][] newa3s = new int[rdim][a3size+MAXDIM+1];

            for (int i=0; i<rdim; i++) {

              System.arraycopy(a3s[i], 0, newa3s[i], 0, a3s[i].length);

            }

            a3s = newa3s;

          }

        }

      }

      for (int i=0; i<cdim; i++) {

        int j = sbt_neigh_simp[t_bn][i][t];

        int j_bn = sbt_neigh_simp_bn[t_bn][i][t];

        if ((j != NOVAL) && sbt_visit[j_bn][j] != vnum) {

          st2[tms] = j;

          st2_bn[tms] = j_bn;

          tms++;

        }

      }

    }

    ret[0] = NOVAL;

  }

  /**

     make neighbor connections between newly created simplices incident to p.

  */

  private void connect(int s, int s_bn) {

    int xb, xf;

    int sb, sb_bn;

    int sf, sf_bn;

    int tf, tf_bn;

    int ccj, ccj_bn;

    int xfi;

    if (s == NOVAL) return;

    for (int i=0; (sbt_neigh_vert[s_bn][i][s] != p) && (i<cdim); i++);

    if (sbt_visit[s_bn][s] == pnum) return;

    sbt_visit[s_bn][s] = pnum;

    ccj = sbt_peak_simp[s_bn][s];

    ccj_bn = sbt_peak_simp_bn[s_bn][s];

    for (xfi=0; (sbt_neigh_simp[ccj_bn][xfi][ccj] != s

              || sbt_neigh_simp_bn[ccj_bn][xfi][ccj] != s_bn)

                     && (xfi<cdim); xfi++);

    for (int i=0; i<cdim; i++) {

      int l;

      if (p == sbt_neigh_vert[s_bn][i][s]) continue;

      sb = sbt_peak_simp[s_bn][s];

      sb_bn = sbt_peak_simp_bn[s_bn][s];

      sf = sbt_neigh_simp[s_bn][i][s];

      sf_bn = sbt_neigh_simp_bn[s_bn][i][s];

      xf = sbt_neigh_vert[ccj_bn][xfi][ccj];

      if (sbt_peak_vert[sf_bn][sf] == NOVAL) {  // are we done already?

        for (l=0; (sbt_neigh_vert[ccj_bn][l][ccj]

                != sbt_neigh_vert[s_bn][i][s]) && (l<cdim); l++);

        sf = sbt_neigh_simp[ccj_bn][l][ccj];

        sf_bn = sbt_neigh_simp_bn[ccj_bn][l][ccj];

        if (sbt_peak_vert[sf_bn][sf] != NOVAL) continue;

      } else do {

        xb = xf;

        for (l=0; (sbt_neigh_simp[sf_bn][l][sf] != sb

                || sbt_neigh_simp_bn[sf_bn][l][sf] != sb_bn)

                && l<cdim; l++);

        xf = sbt_neigh_vert[sf_bn][l][sf];

        sb = sf;

        sb_bn = sf_bn;

        for (l=0; (sbt_neigh_vert[sb_bn][l][sb] != xb) && (l<cdim); l++);

        tf = sbt_neigh_simp[sf_bn][l][sf];

        tf_bn = sbt_neigh_simp_bn[sf_bn][l][sf];

        sf = tf;

        sf_bn = tf_bn;

      } while (sbt_peak_vert[sf_bn][sf] != NOVAL);

      sbt_neigh_simp[s_bn][i][s] = sf;

      sbt_neigh_simp_bn[s_bn][i][s] = sf_bn;

      for (l=0; (sbt_neigh_vert[sf_bn][l][sf] != xf) && (l<cdim); l++);

      sbt_neigh_simp[sf_bn][l][sf] = s;

      sbt_neigh_simp_bn[sf_bn][l][sf] = s_bn;

      connect(sf, sf_bn);

    }

  }

  /**

     visit simplices s with sees(p,s), and make a facet for every neighbor

     of s not seen by p.

  */

  private void make_facets(int seen, int seen_bn, int[] ret, int[] ret_bn) {

    int n, n_bn;

    int q, q_bn;

    int j;

    if (seen == NOVAL) {

      ret[0] = NOVAL;

      return;

    }

    sbt_peak_vert[seen_bn][seen] = p;

    for (int i=0; i<cdim; i++) {

      n = sbt_neigh_simp[seen_bn][i][seen];

      n_bn = sbt_neigh_simp_bn[seen_bn][i][seen];

      if (pnum != sbt_visit[n_bn][n]) {

        sbt_visit[n_bn][n] = pnum;

        if (sees(p, n, n_bn) != 0) make_facets(n, n_bn, voidp, voidp_bn);

      }

      if (sbt_peak_vert[n_bn][n] != NOVAL) continue;

      ns = (simplex_list != NOVAL) ? simplex_list : new_block_simplex();

      ns_bn = simplex_list_bn;

      simplex_list = sbt_next[ns_bn][ns];

      simplex_list_bn = sbt_next_bn[ns_bn][ns];

      sbt_next[ns_bn][ns] = sbt_next[seen_bn][seen];

      sbt_next_bn[ns_bn][ns] = sbt_next_bn[seen_bn][seen];

      sbt_visit[ns_bn][ns] = sbt_visit[seen_bn][seen];

      sbt_mark[ns_bn][ns] = sbt_mark[seen_bn][seen];

      sbt_normal[ns_bn][ns] = sbt_normal[seen_bn][seen];

      sbt_normal_bn[ns_bn][ns] = sbt_normal_bn[seen_bn][seen];

      sbt_peak_vert[ns_bn][ns] = sbt_peak_vert[seen_bn][seen];

      sbt_peak_simp[ns_bn][ns] = sbt_peak_simp[seen_bn][seen];

      sbt_peak_simp_bn[ns_bn][ns] = sbt_peak_simp_bn[seen_bn][seen];

      sbt_peak_basis[ns_bn][ns] = sbt_peak_basis[seen_bn][seen];

      sbt_peak_basis_bn[ns_bn][ns] = sbt_peak_basis_bn[seen_bn][seen];

      for (j=0; j<rdim; j++) {

        sbt_neigh_vert[ns_bn][j][ns] = sbt_neigh_vert[seen_bn][j][seen];

        sbt_neigh_simp[ns_bn][j][ns] = sbt_neigh_simp[seen_bn][j][seen];

        sbt_neigh_simp_bn[ns_bn][j][ns]

                       = sbt_neigh_simp_bn[seen_bn][j][seen];

        sbt_neigh_basis[ns_bn][j][ns] = sbt_neigh_basis[seen_bn][j][seen];

        sbt_neigh_basis_bn[ns_bn][j][ns]

                        = sbt_neigh_basis_bn[seen_bn][j][seen];

      }

      for (j=0; j<cdim; j++) {

        q = sbt_neigh_basis[seen_bn][j][seen];

        q_bn = sbt_neigh_basis_bn[seen_bn][j][seen];

        if (q != NOVAL) bbt_ref_count[q_bn][q]++;

      }

      sbt_visit[ns_bn][ns] = 0;

      sbt_peak_vert[ns_bn][ns] = NOVAL;

      sbt_normal[ns_bn][ns] = NOVAL;

      sbt_peak_simp[ns_bn][ns] = seen;

      sbt_peak_simp_bn[ns_bn][ns] = seen_bn;

      q = sbt_neigh_basis[ns_bn][i][ns];

      q_bn = sbt_neigh_basis_bn[ns_bn][i][ns];

      if (q != NOVAL && --bbt_ref_count[q_bn][q] == 0) {

        bbt_next[q_bn][q] = basis_s_list;

        bbt_next_bn[q_bn][q] = basis_s_list_bn;

        bbt_ref_count[q_bn][q] = 0;

        bbt_lscale[q_bn][q] = 0;

        bbt_sqa[q_bn][q] = 0;

        bbt_sqb[q_bn][q] = 0;

        for (int l=0; l<2*rdim; l++) bbt_vecs[q_bn][l][q] = 0;

        basis_s_list = q;

        basis_s_list_bn = q_bn;

      }

      sbt_neigh_basis[ns_bn][i][ns] = NOVAL;

      sbt_neigh_vert[ns_bn][i][ns] = p;

      for (j=0; (sbt_neigh_simp[n_bn][j][n] != seen

                  || sbt_neigh_simp_bn[n_bn][j][n] != seen_bn)

                  && j<cdim; j++);

      sbt_neigh_simp[seen_bn][i][seen] = sbt_neigh_simp[n_bn][j][n] = ns;

      sbt_neigh_simp_bn[seen_bn][i][seen] = ns_bn;

      sbt_neigh_simp_bn[n_bn][j][n] = ns_bn;

    }

    ret[0] = ns;

    ret_bn[0] = ns_bn;

  }

  /**

     p lies outside flat containing previous sites;

     make p a vertex of every current simplex, and create some new simplices.

  */

  private void extend_simplices(int s, int s_bn, int[] ret, int[] ret_bn) {

    int q, q_bn;

    int ns, ns_bn;

    if (sbt_visit[s_bn][s] == pnum) {

      if (sbt_peak_vert[s_bn][s] != NOVAL) {

        ret[0] = sbt_neigh_simp[s_bn][cdim-1][s];

        ret_bn[0] = sbt_neigh_simp_bn[s_bn][cdim-1][s];

      }

      else {

        ret[0] = s;

        ret_bn[0] = s_bn;

      }

      return;

    }

    sbt_visit[s_bn][s] = pnum;

    sbt_neigh_vert[s_bn][cdim-1][s] = p;

    q = sbt_normal[s_bn][s];

    q_bn = sbt_normal_bn[s_bn][s];

    if (q != NOVAL && --bbt_ref_count[q_bn][q] == 0) {

      bbt_next[q_bn][q] = basis_s_list;

      bbt_next_bn[q_bn][q] = basis_s_list_bn;

      bbt_ref_count[q_bn][q] = 0;

      bbt_lscale[q_bn][q] = 0;

      bbt_sqa[q_bn][q] = 0;

      bbt_sqb[q_bn][q] = 0;

      for (int j=0; j<2*rdim; j++) bbt_vecs[q_bn][j][q] = 0;

      basis_s_list = q;

      basis_s_list_bn = q_bn;

    }

    sbt_normal[s_bn][s] = NOVAL;

    q = sbt_neigh_basis[s_bn][0][s];

    q_bn = sbt_neigh_basis_bn[s_bn][0][s];

    if (q != NOVAL && --bbt_ref_count[q_bn][q] == 0) {

      bbt_next[q_bn][q] = basis_s_list;

      bbt_ref_count[q_bn][q] = 0;

      bbt_lscale[q_bn][q] = 0;

      bbt_sqa[q_bn][q] = 0;

      bbt_sqb[q_bn][q] = 0;

      for (int j=0; j<2*rdim; j++) bbt_vecs[q_bn][j][q] = 0;

      basis_s_list = q;

      basis_s_list_bn = q_bn;

    }

    sbt_neigh_basis[s_bn][0][s] = NOVAL;

    if (sbt_peak_vert[s_bn][s] == NOVAL) {

      int[] esretp = new int[1];

      int[] esretp_bn = new int[1];

      extend_simplices(sbt_peak_simp[s_bn][s],

                       sbt_peak_simp_bn[s_bn][s], esretp, esretp_bn);

      sbt_neigh_simp[s_bn][cdim-1][s] = esretp[0];

      sbt_neigh_simp_bn[s_bn][cdim-1][s] = esretp_bn[0];

      ret[0] = s;

      ret_bn[0] = s_bn;

      return;

    }

    else {

      ns = (simplex_list != NOVAL) ? simplex_list : new_block_simplex();

      ns_bn = simplex_list_bn;

      simplex_list = sbt_next[ns_bn][ns];

      simplex_list_bn = sbt_next_bn[ns_bn][ns];

      sbt_next[ns_bn][ns] = sbt_next[s_bn][s];

      sbt_next_bn[ns_bn][ns] = sbt_next_bn[s_bn][s];

      sbt_visit[ns_bn][ns] = sbt_visit[s_bn][s];

      sbt_mark[ns_bn][ns] = sbt_mark[s_bn][s];

      sbt_normal[ns_bn][ns] = sbt_normal[s_bn][s];

      sbt_normal_bn[ns_bn][ns] = sbt_normal_bn[s_bn][s];

      sbt_peak_vert[ns_bn][ns] = sbt_peak_vert[s_bn][s];

      sbt_peak_simp[ns_bn][ns] = sbt_peak_simp[s_bn][s];

      sbt_peak_simp_bn[ns_bn][ns] = sbt_peak_simp_bn[s_bn][s];

      sbt_peak_basis[ns_bn][ns] = sbt_peak_basis[s_bn][s];

      sbt_peak_basis_bn[ns_bn][ns] = sbt_peak_basis_bn[s_bn][s];

      for (int j=0; j<rdim; j++) {

        sbt_neigh_vert[ns_bn][j][ns] = sbt_neigh_vert[s_bn][j][s];

        sbt_neigh_simp[ns_bn][j][ns] = sbt_neigh_simp[s_bn][j][s];

        sbt_neigh_simp_bn[ns_bn][j][ns] = sbt_neigh_simp_bn[s_bn][j][s];

        sbt_neigh_basis[ns_bn][j][ns] = sbt_neigh_basis[s_bn][j][s];

        sbt_neigh_basis_bn[ns_bn][j][ns] = sbt_neigh_basis_bn[s_bn][j][s];

      }

      for (int j=0; j<cdim; j++) {

        q = sbt_neigh_basis[s_bn][j][s];

        q_bn = sbt_neigh_basis_bn[s_bn][j][s];

        if (q != NOVAL) bbt_ref_count[q_bn][q]++;

      }

      sbt_neigh_simp[s_bn][cdim-1][s] = ns;

      sbt_neigh_simp_bn[s_bn][cdim-1][s] = ns_bn;

      sbt_peak_vert[ns_bn][ns] = NOVAL;

      sbt_peak_simp[ns_bn][ns] = s;

      sbt_peak_simp_bn[ns_bn][ns] = s_bn;

      sbt_neigh_vert[ns_bn][cdim-1][ns] = sbt_peak_vert[s_bn][s];

      sbt_neigh_simp[ns_bn][cdim-1][ns] = sbt_peak_simp[s_bn][s];

      sbt_neigh_simp_bn[ns_bn][cdim-1][ns] = sbt_peak_simp_bn[s_bn][s];

      sbt_neigh_basis[ns_bn][cdim-1][ns] = sbt_peak_basis[s_bn][s];

      sbt_neigh_basis_bn[ns_bn][cdim-1][ns] = sbt_peak_basis_bn[s_bn][s];

      q = sbt_peak_basis[s_bn][s];

      q_bn = sbt_peak_basis_bn[s_bn][s];

      if (q != NOVAL) bbt_ref_count[q_bn][q]++;

      for (int i=0; i<cdim; i++) {

        int[] esretp = new int[1];

        int[] esretp_bn = new int[1];

        extend_simplices(sbt_neigh_simp[ns_bn][i][ns],

                         sbt_neigh_simp_bn[ns_bn][i][ns], esretp, esretp_bn);

        sbt_neigh_simp[ns_bn][i][ns] = esretp[0];

        sbt_neigh_simp_bn[ns_bn][i][ns] = esretp_bn[0];

      }

    }

    ret[0] = ns;

    ret_bn[0] = ns_bn;

    return;

  }

  /**

     return a simplex s that corresponds to a facet of the

     current hull, and sees(p, s).

  */

  private void search(int root, int root_bn, int[] ret, int[] ret_bn) {

    int s, s_bn;

    int tms = 0;

    st[tms] = sbt_peak_simp[root_bn][root];

    st_bn[tms] = sbt_peak_simp_bn[root_bn][root];

    tms++;

    sbt_visit[root_bn][root] = pnum;

    if (sees(p, root, root_bn) == 0) {

      for (int i=0; i<cdim; i++) {

        st[tms] = sbt_neigh_simp[root_bn][i][root];

        st_bn[tms] = sbt_neigh_simp_bn[root_bn][i][root];

        tms++;

      }

    }

    while (tms != 0) {

      if (tms > ss) {

        // JAVA: efficiency issue: how much is this stack hammered?

        ss += ss;

        int[] newst = new int[ss+MAXDIM+1];

        int[] newst_bn = new int[ss+MAXDIM+1];

        System.arraycopy(st, 0, newst, 0, st.length);

        System.arraycopy(st_bn, 0, newst_bn, 0, st_bn.length);

        st = newst;

        st_bn = newst_bn;

      }

      tms--;

      s = st[tms];

      s_bn = st_bn[tms];

      if (sbt_visit[s_bn][s] == pnum) continue;

      sbt_visit[s_bn][s] = pnum;

      if (sees(p, s, s_bn) == 0) continue;

      if (sbt_peak_vert[s_bn][s] == NOVAL) {

        ret[0] = s;

        ret_bn[0] = s_bn;

        return;

      }

      for (int i=0; i<cdim; i++) {

        st[tms] = sbt_neigh_simp[s_bn][i][s];

        st_bn[tms] = sbt_neigh_simp_bn[s_bn][i][s];

        tms++;

      }

    }

    ret[0] = NOVAL;

    return;

  }

  // <<<< Constructor >>>>

  public DelaunayClarkson(double[][] samples) throws VisADException {

    int s, s_bn, q, q_bn;

    int root, root_bn;

    int k=0;

    int[] retp = new int[1];

    int[] retp_bn = new int[1];

    int[] ret2p = new int[1];

    int[] ret2p_bn = new int[1];

    int[] curt = new int[1];

    int[] curt_bn = new int[1];

    int s_num = 0;

    int nrs;

    // Start of main hull triangulation algorithm

    dim = samples.length;

    nrs = samples[0].length;

    for (int i=1; i<dim; i++) nrs = Math.min(nrs, samples[i].length);

    if (nrs <= dim) throw new SetException("DelaunayClarkson: "

                                          +"not enough samples");

    if (dim > MAXDIM) throw new SetException("DelaunayClarkson: "

                               +"dimension bound MAXDIM exceeded");

    // copy samples

    site_blocks = new double[dim][nrs];

    for (int j=0; j<dim; j++) {

      System.arraycopy(samples[j], 0, site_blocks[j], 0, nrs);

    }

/* WLH 29 Jan 98 - scale samples values as discussed in Delaunay.factory

    for (int j=0; j<dim; j++) {

      for (int kk=0; kk<nrs; kk++) {

        site_blocks[j][kk] = 100.0f * samples[j][kk];

      }

    }

*/

    exact_bits = (int) (DBL_MANT_DIG*Math.log(FLT_RADIX)/ln2);

    b_err_min = DBL_EPSILON*MAXDIM*(1<<MAXDIM)*MAXDIM*3.01;

    b_err_min_sq = b_err_min * b_err_min;

    cdim = 0;

    rdim = dim+1;

    if (rdim > MAXDIM) throw new SetException(

              "dimension bound MAXDIM exceeded; rdim="+rdim+"; dim="+dim);

    pnb = basis_s_list != NOVAL ? basis_s_list : new_block_basis_s();

    pnb_bn = basis_s_list_bn;

    basis_s_list = bbt_next[pnb_bn][pnb];

    basis_s_list_bn = bbt_next_bn[pnb_bn][pnb];

    bbt_next[pnb_bn][pnb] = NOVAL;

    ttbp = basis_s_list != NOVAL ? basis_s_list : new_block_basis_s();

    ttbp_bn = basis_s_list_bn;

    basis_s_list = bbt_next[ttbp_bn][ttbp];

    basis_s_list_bn = bbt_next_bn[ttbp_bn][ttbp];

    bbt_next[ttbp_bn][ttbp] = NOVAL;

    bbt_ref_count[ttbp_bn][ttbp] = 1;

    bbt_lscale[ttbp_bn][ttbp] = -1;

    bbt_sqa[ttbp_bn][ttbp] = 0;

    bbt_sqb[ttbp_bn][ttbp] = 0;

    for (int j=0; j<2*rdim; j++) bbt_vecs[ttbp_bn][j][ttbp] = 0;

    root = NOVAL;

    p = INFINITY;

    ib = (basis_s_list != NOVAL) ? basis_s_list : new_block_basis_s();

    ib_bn = basis_s_list_bn;

    basis_s_list = bbt_next[ib_bn][ib];

    basis_s_list_bn = bbt_next_bn[ib_bn][ib];

    bbt_ref_count[ib_bn][ib] = 1;

    bbt_vecs[ib_bn][2*rdim-1][ib] = bbt_vecs[ib_bn][rdim-1][ib] = 1;

    bbt_sqa[ib_bn][ib] = bbt_sqb[ib_bn][ib] = 1;

    root = (simplex_list != NOVAL) ? simplex_list : new_block_simplex();

    root_bn = simplex_list_bn;

    simplex_list = sbt_next[root_bn][root];

    simplex_list_bn = sbt_next_bn[root_bn][root];

    ch_root = root;

    ch_root_bn = root_bn;

    s = (simplex_list != NOVAL) ? simplex_list : new_block_simplex();

    s_bn = simplex_list_bn;

    simplex_list = sbt_next[s_bn][s];

    simplex_list_bn = sbt_next_bn[s_bn][s];

    sbt_next[s_bn][s] = sbt_next[root_bn][root];

    sbt_next_bn[s_bn][s] = sbt_next_bn[root_bn][root];

    sbt_visit[s_bn][s] = sbt_visit[root_bn][root];

    sbt_mark[s_bn][s] = sbt_mark[root_bn][root];

    sbt_normal[s_bn][s] = sbt_normal[root_bn][root];

    sbt_normal_bn[s_bn][s] = sbt_normal_bn[root_bn][root];

    sbt_peak_vert[s_bn][s] = sbt_peak_vert[root_bn][root];

    sbt_peak_simp[s_bn][s] = sbt_peak_simp[root_bn][root];

    sbt_peak_simp_bn[s_bn][s] = sbt_peak_simp_bn[root_bn][root];

    sbt_peak_basis[s_bn][s] = sbt_peak_basis[root_bn][root];

    sbt_peak_basis_bn[s_bn][s] = sbt_peak_basis_bn[root_bn][root];

    for (int i=0; i<rdim; i++) {

      sbt_neigh_vert[s_bn][i][s] = sbt_neigh_vert[root_bn][i][root];

      sbt_neigh_simp[s_bn][i][s] = sbt_neigh_simp[root_bn][i][root];

      sbt_neigh_simp_bn[s_bn][i][s] = sbt_neigh_simp_bn[root_bn][i][root];

      sbt_neigh_basis[s_bn][i][s] = sbt_neigh_basis[root_bn][i][root];

      sbt_neigh_basis_bn[s_bn][i][s] = sbt_neigh_basis_bn[root_bn][i][root];

    }

    for (int i=0;i<cdim;i++) {

      q = sbt_neigh_basis[root_bn][i][root];

      q_bn = sbt_neigh_basis_bn[root_bn][i][root];

      if (q != NOVAL) bbt_ref_count[q_bn][q]++;

    }

    sbt_peak_vert[root_bn][root] = p;

    sbt_peak_simp[root_bn][root] = s;

    sbt_peak_simp_bn[root_bn][root] = s_bn;

    sbt_peak_simp[s_bn][s] = root;

    sbt_peak_simp_bn[s_bn][s] = root_bn;

    while (cdim < rdim) {

      int oof = 0;

      if (s_num == 0) p = 0;

      else p++;

      for (int i=0; i<dim; i++) {

        site_blocks[i][p] = (double) Math.floor(site_blocks[i][p]+0.5);

      }

      s_num++;

      pnum = (s_num*dim-1)/dim + 2;

      cdim++;

      sbt_neigh_vert[root_bn][cdim-1][root] = sbt_peak_vert[root_bn][root];

      q = sbt_neigh_basis[root_bn][cdim-1][root];

      q_bn = sbt_neigh_basis_bn[root_bn][cdim-1][root];

      if (q != NOVAL && --bbt_ref_count[q_bn][q] == 0) {

        bbt_next[q_bn][q] = basis_s_list;

        bbt_next_bn[q_bn][q] = basis_s_list_bn;

        bbt_ref_count[q_bn][q] = 0;

        bbt_lscale[q_bn][q] = 0;

        bbt_sqa[q_bn][q] = 0;

        bbt_sqb[q_bn][q] = 0;

        for (int l=0; l<2*rdim; l++) bbt_vecs[q_bn][l][q] = 0;

        basis_s_list = q;

        basis_s_list_bn = q_bn;

      }

      sbt_neigh_basis[root_bn][cdim-1][root] = NOVAL;

      get_basis_sede(root, root_bn);

      if (sbt_neigh_vert[root_bn][0][root] == INFINITY) oof = 1;

      else {

        curt[0] = pnb;

        curt_bn[0] = pnb_bn;

        reduce(curt, curt_bn, p, root, root_bn, cdim);

        pnb = curt[0];

        pnb_bn = curt_bn[0];

        if (bbt_sqa[pnb_bn][pnb] != 0) oof = 1;

        else cdim--;

      }

      if (oof != 0) extend_simplices(root, root_bn, voidp, voidp_bn);

      else {

        search(root, root_bn, retp, retp_bn);

        make_facets(retp[0], retp_bn[0], ret2p, ret2p_bn);

        connect(ret2p[0], ret2p_bn[0]);

      }

    }

    for (int i=s_num; i<nrs; i++) {

      p++;

      s_num++;

      for (int j=0; j<dim; j++) {

        site_blocks[j][p] = (double) Math.floor(site_blocks[j][p]+0.5);

      }

      pnum = (s_num*dim-1)/dim + 2;

      search(root, root_bn, retp, retp_bn);

      make_facets(retp[0], retp_bn[0], ret2p, ret2p_bn);

      connect(ret2p[0], ret2p_bn[0]);

    }

    a3size = rdim*nrs;

    a3s = new int[rdim][a3size+MAXDIM+1];

    visit_triang_gen(root, root_bn, 1, retp, retp_bn);

    visit_triang_gen(retp[0], retp_bn[0], 0, voidp, voidp_bn);

    // deallocate memory

    /* NOTE: If this deallocation is not performed, more points

       could theoretically be added to the triangulation later */

    site_blocks = null;

    st = null;

    st_bn = null;

    st2 = null;

    st2_bn = null;

    sbt_next = null;

    sbt_next_bn = null;

    sbt_visit = null;

    sbt_mark = null;

    sbt_normal = null;

    sbt_normal_bn = null;

    sbt_peak_vert = null;

    sbt_peak_simp = null;

    sbt_peak_simp_bn = null;

    sbt_peak_basis = null;

    sbt_peak_basis_bn = null;

    sbt_neigh_vert = null;

    sbt_neigh_simp = null;

    sbt_neigh_simp_bn = null;

    sbt_neigh_basis = null;

    sbt_neigh_basis_bn = null;

    bbt_next = null;

    bbt_next_bn = null;

    bbt_ref_count = null;

    bbt_lscale = null;

    bbt_sqa = null;

    bbt_sqb = null;

    bbt_vecs = null;

/* ********** END OF CONVERTED HULL CODE ********** */

/*          (but still inside constructor)          */

    // compute number of triangles or tetrahedra

    int[] nverts = new int[nrs];

    for (int i=0; i<nrs; i++) nverts[i] = 0;

    int ntris = 0;

    boolean positive;

    for (int i=0; i<nts; i++) {

      positive = true;

      for (int j=0; j<rdim; j++) {

        if (a3s[j][i] < 0) positive = false;

      }

      if (positive) {

        ntris++;

        for (int j=0; j<rdim; j++) nverts[a3s[j][i]]++;

      }

    }

    Vertices = new int[nrs][];

    for (int i=0; i<nrs; i++) Vertices[i] = new int[nverts[i]];

    for (int i=0; i<nrs; i++) nverts[i] = 0;

    // build Tri & Vertices components

    Tri = new int[ntris][rdim];

    int a, b, c, d;

    int itri = 0;

    for (int i=0; i<nts; i++) {

      positive = true;

      for (int j=0; j<rdim; j++) {

        if (a3s[j][i] < 0) positive = false;

      }

      if (positive) {

        for (int j=0; j<rdim; j++) {

          Vertices[a3s[j][i]][nverts[a3s[j][i]]++] = itri;

          Tri[itri][j] = a3s[j][i];

        }

        itri++;

      }

    }

    // Deallocate remaining helper information

    a3s = null;

    // call more generic method for constructing Walk and Edges arrays

    finish_triang(samples);

  }

}