//
// DelaunayClarkson.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;
/* 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.
*/
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 float[][] 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; i0; 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= 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 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; i0; i--) {
int m = (i > 1) ? sbt_neigh_vert[ch_root_bn][i-2][ch_root]
: INFINITY;
int j;
for (j=0; j -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 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 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 a3size) {
// JAVA: efficiency issue, hammering an array
a3size += a3size;
int[][] newa3s = new int[rdim][a3size+MAXDIM+1];
for (int i=0; i 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 MAXDIM) throw new SetException("DelaunayClarkson: "
+"dimension bound MAXDIM exceeded");
// copy samples
site_blocks = new float[dim][nrs];
for (int j=0; j 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
Ads help maintain this website.