//
// DelaunayWatson.java
//

/*
VisAD system for interactive analysis and visualization of numerical
data.  Copyright (C) 1996 - 2017 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 visad;

import java.util.*;

/* The Delaunay triangulation/tetrahedralization algorithm in this class is
 * originally from nnsort.c by David F. Watson:
 *
 * nnsort() finds the Delaunay triangulation of the two- or three-component
 * vectors in 'data_list' and returns a list of simplex vertices in
 * 'vertices' with the corresponding circumcentre and squared radius in the
 * rows of 'circentres'.  nnsort() also can be used to find the ordered
 * convex hull of the two- or three-component vectors in 'data_list' and
 * returns a list of (d-1)-facet vertices in 'vertices' (dummy filename for
 * 'circentres' must be used).
 * nnsort() was written by Dave Watson and uses the algorithm described in -
 *    Watson, D.F., 1981, Computing the n-dimensional Delaunay tessellation
 *          with application to Voronoi polytopes:
 *                                      The Computer J., 24(2), p. 167-172.
 *
 * additional information about this algorithm can be found in -
 *    CONTOURING: A guide to the analysis and display of spatial data,
 *    by David F. Watson, Pergamon Press, 1992, ISBN 0 08 040286 0
 *                                                                        */

/**
   DelaunayWatson represents an O(N^2) method with low overhead
   to find the Delaunay triangulation or tetrahedralization of
   a set of samples of R^2 or R^3.

*/
public class DelaunayWatson extends Delaunay {

  private static final float BIGNUM = (float) 1E37;
  private static final float EPSILON = 0.00001f;
  // temporary storage size factor
  private static final int TSIZE = 75;
  // factor (>=1) for radius of control points
  private static final float RANGE = 10.0f;

  /**
   * construct a Delaunay triangulation of the points in the
   * samples array using Watson's algorithm
   * @param samples locations of points for topology - dimensioned
   *                float[dimension][number_of_points]
   * @throws VisADException a VisAD error occurred
   */
  public DelaunayWatson(float[][] samples) throws VisADException {
    int dim = samples.length;
    int nrs = samples[0].length;

    float xx, yy, bgs;
    int i0, i1, i2, i3, i4, i5, i6, i7, i8, i9, i11;
    int[] ii = new int[3];
    int dm, dim1, nts, tsz;

    float[][] mxy = new float[2][dim];
    for (i0=0; i0 pts[i0][i1]) mxy[1][i1] = pts[i0][i1]; // min
      }
    }

    for (bgs=0, i0=0; i0dm) {
            i8 = i7;
Loop1:
            for (i3=0; i3<=i8; i3++) {
              for (i5=0; i5 tsz) {
            int newtsz = 2 * tsz;
            int[][] newtmp = new int[newtsz + 1][dim];
            System.arraycopy(tmp, 0, newtmp, 0, tsz);
            tsz = newtsz;
            tmp = newtmp;
            // WLH 23 july 97
            // throw new VisADException(
            //                "DelaunayWatson: Temporary storage exceeded");
          }
          for (i3=0; i3 EPSILON) {
          if (xx < 0)
            fprintf(afile,"%d %d %d\n",a3s[i0][0],a3s[i0][2],a3s[i0][1]);
          else fprintf(afile,"%d %d %d\n",a3s[i0][0],a3s[i0][1],a3s[i0][2]);
          fprintf(bfile,"%e %e %e\n",ccr[i0][0],ccr[i0][1],ccr[i0][2]);
        }
      }
      else {
        xx = ((wrk[0][0] * (wrk[1][1] * wrk[2][2] - wrk[2][1] * wrk[1][2]))
           -  (wrk[0][1] * (wrk[1][0] * wrk[2][2] - wrk[2][0] * wrk[1][2]))
           +  (wrk[0][2] * (wrk[1][0] * wrk[2][1] - wrk[2][0] * wrk[1][1])));
        if (fabs(xx) > EPSILON) {
          if (xx < 0)
            fprintf(afile,"%d %d %d %d\n",
                    a3s[i0][0],a3s[i0][1],a3s[i0][3],a3s[i0][2]);
          else fprintf(afile,"%d %d %d %d\n",
                       a3s[i0][0],a3s[i0][1],a3s[i0][2],a3s[i0][3]);
          fprintf(bfile,"%e %e %e %e\n",
                  ccr[i0][0],ccr[i0][1],ccr[i0][2],ccr[i0][3]);
        }
      }
    }
  }
*/

}