//
// Irregular3DSet.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;
/**
Irregular3DSet represents a finite set of samples of R^3.
No Irregular3DSet with ManifoldDimension = 1. Use
Gridded3DSet with ManifoldDimension = 1 instead.
*/
public class Irregular3DSet extends IrregularSet {
private float LowX, HiX, LowY, HiY, LowZ, HiZ;
/** a 3-D irregular set with null errors, CoordinateSystem
and Units are defaults from type; topology is computed
by the constructor */
public Irregular3DSet(MathType type, float[][] samples)
throws VisADException {
this(type, samples, null, null, null, null, true);
}
/** a 3-D irregular set; samples array is organized
float[3][number_of_samples]; no geometric constraint on
samples; if delan is non-null it defines the topology of
samples (which may have manifold dimension 2 or 3), else
the constructor computes a topology with manifold dimension
3; note that Gridded3DSet can be used for an irregular set
with domain dimension 3 and manifold dimension 1;
coordinate_system and units must be compatible with
defaults for type, or may be null; errors may be null */
public Irregular3DSet(MathType type, float[][] samples,
CoordinateSystem coord_sys, Unit[] units,
ErrorEstimate[] errors, Delaunay delan)
throws VisADException {
this(type, samples, coord_sys, units, errors, delan, true);
}
public Irregular3DSet(MathType type, float[][] samples,
CoordinateSystem coord_sys, Unit[] units,
ErrorEstimate[] errors, Delaunay delan,
boolean copy) throws VisADException {
/* ManifoldDimension might not be equal to samples.length
if a 2D triangulation has been specified */
super(type, samples, (delan == null) ? samples.length
: delan.Tri[0].length-1,
coord_sys, units, errors, delan, copy);
LowX = Low[0];
HiX = Hi[0];
LowY = Low[1];
HiY = Hi[1];
LowZ = Low[2];
HiZ = Hi[2];
oldToNew = null;
newToOld = null;
}
/** construct Irregular3DSet using sort from existing
Irregular1DSet */
public Irregular3DSet(MathType type, float[][] samples,
int[] new2old, int[] old2new) throws VisADException {
this(type, samples, new2old, old2new, null, null, null, true);
}
/** construct Irregular3DSet using sort from existing
Irregular1DSet */
public Irregular3DSet(MathType type, float[][] samples,
int[] new2old, int[] old2new,
CoordinateSystem coord_sys, Unit[] units,
ErrorEstimate[] errors) throws VisADException {
this(type, samples, new2old, old2new, coord_sys, units, errors, true);
}
public Irregular3DSet(MathType type, float[][] samples,
int[] new2old, int[] old2new,
CoordinateSystem coord_sys, Unit[] units,
ErrorEstimate[] errors, boolean copy)
throws VisADException {
super(type, samples, 1, coord_sys, units, errors, null, copy);
if (Length != new2old.length || Length != old2new.length) {
throw new SetException("Irregular3DSet: sort lengths do not match");
}
newToOld = new int[Length];
oldToNew = new int[Length];
System.arraycopy(new2old, 0, newToOld, 0, Length);
System.arraycopy(old2new, 0, oldToNew, 0, Length);
LowX = Low[0];
HiX = Hi[0];
LowY = Low[1];
HiY = Hi[1];
LowZ = Low[2];
HiZ = Hi[2];
Delan = null;
}
public Set makeSpatial(SetType type, float[][] samples) throws VisADException {
if (samples.length == 3) {
if (ManifoldDimension == 1) {
return new Irregular3DSet(type, samples, newToOld, oldToNew,
null, null, null, false);
}
else {
if (Delan.Tri == null || Delan.Tri.length == 0) return null;
return new Irregular3DSet(type, samples, null, null, null,
Delan, false);
}
}
else {
throw new SetException("Irregular3DSet.makeSpatial: bad samples length");
}
}
/** convert an array of 1-D indices to an array of values in
R^DomainDimension */
public float[][] indexToValue(int[] index) throws VisADException {
float[][] value = new float[3][index.length];
float[][]mySamples = getMySamples();
for (int i=0; i= 0) && (index[i] < Length) ) {
value[0][i] = mySamples[0][index[i]];
value[1][i] = mySamples[1][index[i]];
value[2][i] = mySamples[2][index[i]];
}
else {
value[0][i] = value[1][i] = value[2][i] = Float.NaN;
}
}
return value;
}
/** valueToTri returns an array of containing triangles given
an array of points in R^DomainDimension */
public int[] valueToTri(float[][] value) throws VisADException {
if (ManifoldDimension != 3) {
throw new SetException("Irregular3DSet.valueToTri: " +
"ManifoldDimension must be 3, not " +
ManifoldDimension);
}
int length = value[0].length;
if (length != value[1].length || length != value[2].length) {
throw new SetException("Irregular3DSet.valueToTri: lengths " +
"don't match");
}
boolean nonConvex = Delan.getNonConvex();
float[] PA = new float[3];
float[] PB = new float[3];
float[] PC = new float[3];
float[] PD = new float[3];
float[] BAxCB = new float[3];
float[] CBxDC = new float[3];
float[] DCxAD = new float[3];
float[] ADxBA = new float[3];
float sum_BAxCB;
float sum_CBxDC;
float sum_DCxAD;
float sum_ADxBA;
boolean[] fail_tri = new boolean[Delan.Tri.length];
for ( int kk = 0; kk < fail_tri.length; kk++ ) {
fail_tri[kk] = false;
}
int [] fail_list = null;
int fail_length = 0;
int[] tri = new int[length];
int curtri = 0;
// System.out.println("length = " + length + " Delan.Tri.length = " +
// Delan.Tri.length);
float[][]mySamples = getMySamples();
for (int i=0; i 0) == (
BAxCB[0]*(Dx-Ax)
+ BAxCB[1]*(Dy-Ay)
+ BAxCB[2]*(Dz-Az) > 0) )) && (sum_BAxCB != 0);
boolean test2 = ((tval2 == 0.0f) || ( (tval2 > 0) == (
CBxDC[0]*(Ax-Bx)
+ CBxDC[1]*(Ay-By)
+ CBxDC[2]*(Az-Bz) > 0) )) && (sum_CBxDC != 0);
boolean test3 = ((tval3 == 0.0f) || ( (tval3 > 0) == (
DCxAD[0]*(Bx-Cx)
+ DCxAD[1]*(By-Cy)
+ DCxAD[2]*(Bz-Cz) > 0) )) && (sum_DCxAD != 0);
boolean test4 = ((tval4 == 0.0f) || ( (tval4 > 0) == (
ADxBA[0]*(Cx-Dx)
+ ADxBA[1]*(Cy-Dy)
+ ADxBA[2]*(Cz-Dz) > 0) )) && (sum_ADxBA != 0);
// System.out.println("i: "+i+" curtri: "+curtri+" test1: "+test1+" test2: "+test2+" test3: "+test3+" test4: "+test4);
// figure out which triangle to go to next
if (!test1 || !test2 || !test3 || !test4) {
// record failed tri
fail_tri[curtri] = true;
// add to list of failed tris for efficient reset
if (fail_list == null) {
fail_list = new int[4];
fail_length = 0;
}
else if (fail_length >= fail_list.length) {
int[] new_fail_list = new int[2 * fail_list.length];
System.arraycopy(fail_list, 0, new_fail_list, 0, fail_list.length);
fail_list = new_fail_list;
}
fail_list[fail_length] = curtri;
fail_length++;
int t = -1;
boolean fail = true;
if (!test1 && fail) {
t = Delan.Walk[curtri][0];
if (t != -1) fail = fail_tri[t];
}
if (!test2 && fail) {
t = Delan.Walk[curtri][1];
if (t != -1) fail = fail_tri[t];
}
if (!test3 && fail) {
t = Delan.Walk[curtri][2];
if (t != -1) fail = fail_tri[t];
}
if (!test4 && fail) {
t = Delan.Walk[curtri][3];
if (t != -1) fail = fail_tri[t];
}
if (!fail || t == -1) {
curtri = t;
}
if (fail) curtri = -1;
if (nonConvex) {
// to deal with non-convex Set, but very slow
if (curtri == -1) {
for (int jj=0; jj 0) {
if (color_length > 3) color_length = 3; // no alpha for lines
color_levels = new byte[color_length][maxv];
}
float[] vx = new float[maxv];
float[] vy = new float[maxv];
float[] vz = new float[maxv];
int numv = 0;
for (int jj=0; jj 0) {
auxa = new byte[color_length];
auxb = new byte[color_length];
auxc = new byte[color_length];
for (int i=0; i gb ? ga : gb;
gx = gc > gx ? gc : gx;
for (int il=0; il= maxv) {
// allocate more space
maxv = 2 * maxv;
byte[][] t = color_levels;
color_levels = new byte[color_length][maxv];
for (int i=0; i gx) break;
if (gg < lowlimit) continue;
if (gg > highlimit) break;
// compute orientation of lines inside box
ii = 0;
if (gg > ga) ii = 1;
if (gg > gb) ii += 2;
if (gg > gc) ii += 4;
if (ii > 3) ii = 7 - ii;
if (ii <= 0) continue;
switch (ii) {
case 1:
gba = gb-ga;
gca = gc-ga;
ratioba = (gg-ga)/gba;
ratioca = (gg-ga)/gca;
if (color_length > 0) {
for (int i=0; i 255) ? -1 : ((t < 128) ? t : t - 256) ) );
t = (int) ( (1.0f - ratioca) * ((auxa[i] < 0) ?
((float) auxa[i]) + 256.0f : ((float) auxa[i]) ) +
ratioca * ((auxc[i] < 0) ?
((float) auxc[i]) + 256.0f : ((float) auxc[i]) ) );
color_levels[i][numv+1] = (byte)
( (t < 0) ? 0 : ((t > 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
color_levels[i][numv] = auxa[i] + (auxb[i]-auxa[i]) * ratioba;
color_levels[i][numv+1] = auxa[i] + (auxc[i]-auxa[i]) * ratioca;
*/
}
}
vx[numv] = samples[0][va] + (samples[0][vb]-samples[0][va]) * ratioba;
vy[numv] = samples[1][va] + (samples[1][vb]-samples[1][va]) * ratioba;
vz[numv] = samples[2][va] + (samples[2][vb]-samples[2][va]) * ratioba;
numv++;
vx[numv] = samples[0][va] + (samples[0][vc]-samples[0][va]) * ratioca;
vy[numv] = samples[1][va] + (samples[1][vc]-samples[1][va]) * ratioca;
vz[numv] = samples[2][va] + (samples[2][vc]-samples[2][va]) * ratioca;
numv++;
break;
case 2:
gba = gb-ga;
gcb = gc-gb;
ratioba = (gg-ga)/gba;
ratiocb = (gg-gb)/gcb;
if (color_length > 0) {
for (int i=0; i 255) ? -1 : ((t < 128) ? t : t - 256) ) );
t = (int) ( (1.0f - ratiocb) * ((auxb[i] < 0) ?
((float) auxb[i]) + 256.0f : ((float) auxb[i]) ) +
ratiocb * ((auxc[i] < 0) ?
((float) auxc[i]) + 256.0f : ((float) auxc[i]) ) );
color_levels[i][numv+1] = (byte)
( (t < 0) ? 0 : ((t > 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
color_levels[i][numv] = auxa[i] + (auxb[i]-auxa[i]) * ratioba;
color_levels[i][numv+1] = auxb[i] + (auxc[i]-auxb[i]) * ratiocb;
*/
}
}
vx[numv] = samples[0][va] + (samples[0][vb]-samples[0][va]) * ratioba;
vy[numv] = samples[1][va] + (samples[1][vb]-samples[1][va]) * ratioba;
vz[numv] = samples[2][va] + (samples[2][vb]-samples[2][va]) * ratioba;
numv++;
vx[numv] = samples[0][vb] + (samples[0][vc]-samples[0][vb]) * ratiocb;
vy[numv] = samples[1][vb] + (samples[1][vc]-samples[1][vb]) * ratiocb;
vz[numv] = samples[2][vb] + (samples[2][vc]-samples[2][vb]) * ratiocb;
numv++;
break;
case 3:
gca = gc-ga;
gcb = gc-gb;
ratioca = (gg-ga)/gca;
ratiocb = (gg-gb)/gcb;
if (color_length > 0) {
for (int i=0; i 255) ? -1 : ((t < 128) ? t : t - 256) ) );
t = (int) ( (1.0f - ratiocb) * ((auxb[i] < 0) ?
((float) auxb[i]) + 256.0f : ((float) auxb[i]) ) +
ratiocb * ((auxc[i] < 0) ?
((float) auxc[i]) + 256.0f : ((float) auxc[i]) ) );
color_levels[i][numv+1] = (byte)
( (t < 0) ? 0 : ((t > 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
color_levels[i][numv] = auxa[i] + (auxc[i]-auxa[i]) * ratioca;
color_levels[i][numv+1] = auxb[i] + (auxc[i]-auxb[i]) * ratiocb;
*/
}
}
vx[numv] = samples[0][va] + (samples[0][vc]-samples[0][va]) * ratioca;
vy[numv] = samples[1][va] + (samples[1][vc]-samples[1][va]) * ratioca;
vz[numv] = samples[2][va] + (samples[2][vc]-samples[2][va]) * ratioca;
numv++;
vx[numv] = samples[0][vb] + (samples[0][vc]-samples[0][vb]) * ratiocb;
vy[numv] = samples[1][vb] + (samples[1][vc]-samples[1][vb]) * ratiocb;
vz[numv] = samples[2][vb] + (samples[2][vc]-samples[2][vb]) * ratiocb;
numv++;
break;
} // end switch (ii)
} // end for (int il=0; il 0) {
colors = new byte[3][numv];
System.arraycopy(color_levels[0], 0, colors[0], 0, numv);
System.arraycopy(color_levels[1], 0, colors[1], 0, numv);
System.arraycopy(color_levels[2], 0, colors[2], 0, numv);
/* MEM_WLH
colors = new byte[3][numv];
for (int i=0; i<3; i++) {
for (int j=0; j 255) ? 255 : k;
colors[i][j] = (byte) ((k < 128) ? k : k - 256);
}
}
*/
color_levels = null;
}
setGeometryArray(lineArray, coordinates, 4, colors);
return new VisADLineArray[][] {{lineArray}, null, null};
}
public VisADGeometryArray makeIsoSurface(float isolevel,
float[] fieldValues, byte[][] color_values, boolean indexed)
throws VisADException {
if (ManifoldDimension != 3) {
throw new DisplayException("Irregular3DSet.main_isosurf: " +
"ManifoldDimension must be 3, not " +
ManifoldDimension);
}
float[][] fieldVertices = new float[3][];
byte[][] color_levels = null;
if (color_values != null) {
color_levels = new byte[color_values.length][];
/* MEM_WLH
cfloat = new float[color_values.length][];
for (int i = 0; i< color_levels.length; i++) {
if (color_values[i] != null) {
cfloat[i] = new float[color_values[i].length];
for (int j=0; j 255) ? 255 : k;
c[i][j] = (byte) ((k < 128) ? k : k - 256);
}
}
}
// FREE
color_levels = null;
}
*/
int nvertex = vertToPoly[0].length;
int npolygons = polyToVert[0].length;
float[] NX = new float[nvertex];
float[] NY = new float[nvertex];
float[] NZ = new float[nvertex];
if (nvertex == 0 || npolygons == 0) return null;
// with make_normals
float[] NxA = new float[npolygons];
float[] NxB = new float[npolygons];
float[] NyA = new float[npolygons];
float[] NyB = new float[npolygons];
float[] NzA = new float[npolygons];
float[] NzB = new float[npolygons];
float[] Pnx = new float[npolygons];
float[] Pny = new float[npolygons];
float[] Pnz = new float[npolygons];
make_normals(fieldVertices[0], fieldVertices[1], fieldVertices[2],
NX, NY, NZ, nvertex, npolygons, Pnx, Pny, Pnz,
NxA, NxB, NyA, NyB, NzA, NzB, vertToPoly[0], polyToVert[0]);
// take the garbage out
NxA = NxB = NyA = NyB = NzA = NzB = Pnx = Pny = Pnz = null;
// with poly_triangle_stripe
float[] normals = new float[3 * nvertex];
int j = 0;
for (int i=0; i 0) {
if (auxLevels == null || auxLevels.length != naux) {
throw new DisplayException("Irregular3DSet.makeIsosurface: "
+"auxLevels length doesn't match");
}
for (int i=0; i " + fieldValues[i]);
}
System.out.println(Delan.sampleString(getMySamples()));
}
int trilength = Delan.Tri.length;
// temporary storage of polyToVert structure
int[][] polys = new int[trilength][4];
// pointers from global edge number to nvertex
int[] globalToVertex = new int[Delan.NumEdges];
for (int i=0; i 0) ? new byte[naux][Delan.NumEdges] : null;
int t;
int nvertex = 0;
int npolygons = 0;
float[][]mySamples = getMySamples();
for (int i=0; i isolevel) ? 1 : 0)
+ ((f1 > isolevel) ? 2 : 0)
+ ((f2 > isolevel) ? 4 : 0)
+ ((f3 > isolevel) ? 8 : 0);
// apply signature to index
if (sig < 0.0f) index = 15 - index;
switch (index) {
case 0:
case 15: // plane does not intersect this tetrahedron
break;
case 1:
case 14: // plane slices a triangle
// define edge values needed
e0 = Delan.Edges[i][0];
e1 = Delan.Edges[i][1];
e2 = Delan.Edges[i][2];
// fill in edge interpolation values if they haven't been found
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e0])) {
*/
// test for missing
if (edgeInterp[0][e0] != edgeInterp[0][e0]) {
float a = (isolevel - f1)/(f0 - f1);
if (a < 0) a = -a;
edgeInterp[0][e0] = (float) a*mySamples[0][v0] + (1-a)*mySamples[0][v1];
edgeInterp[1][e0] = (float) a*mySamples[1][v0] + (1-a)*mySamples[1][v1];
edgeInterp[2][e0] = (float) a*mySamples[2][v0] + (1-a)*mySamples[2][v1];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e0] = (float) a*auxValues[j][v0] + (1-a)*auxValues[j][v1];
*/
}
globalToVertex[e0] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e1])) {
*/
// test for missing
if (edgeInterp[0][e1] != edgeInterp[0][e1]) {
float a = (isolevel - f2)/(f0 - f2);
if (a < 0) a = -a;
edgeInterp[0][e1] = (float) a*mySamples[0][v0] + (1-a)*mySamples[0][v2];
edgeInterp[1][e1] = (float) a*mySamples[1][v0] + (1-a)*mySamples[1][v2];
edgeInterp[2][e1] = (float) a*mySamples[2][v0] + (1-a)*mySamples[2][v2];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e1] = (float) a*auxValues[j][v0] + (1-a)*auxValues[j][v2];
*/
}
globalToVertex[e1] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e2])) {
*/
// test for missing
if (edgeInterp[0][e2] != edgeInterp[0][e2]) {
float a = (isolevel - f3)/(f0 - f3);
if (a < 0) a = -a;
edgeInterp[0][e2] = (float) a*mySamples[0][v0] + (1-a)*mySamples[0][v3];
edgeInterp[1][e2] = (float) a*mySamples[1][v0] + (1-a)*mySamples[1][v3];
edgeInterp[2][e2] = (float) a*mySamples[2][v0] + (1-a)*mySamples[2][v3];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e2] = (float) a*auxValues[j][v0] + (1-a)*auxValues[j][v3];
*/
}
globalToVertex[e2] = nvertex;
nvertex++;
}
// fill in the polys and vertToPoly arrays
polys[npolygons][0] = e0;
if (index == 1) {
polys[npolygons][1] = e1;
polys[npolygons][2] = e2;
}
else { // index == 14
polys[npolygons][1] = e2;
polys[npolygons][2] = e1;
}
polys[npolygons][3] = -1;
// on to the next tetrahedron
npolygons++;
break;
case 2:
case 13: // plane slices a triangle
// define edge values needed
e0 = Delan.Edges[i][0];
e3 = Delan.Edges[i][3];
e4 = Delan.Edges[i][4];
// fill in edge interpolation values if they haven't been found
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e0])) {
*/
// test for missing
if (edgeInterp[0][e0] != edgeInterp[0][e0]) {
float a = (isolevel - f1)/(f0 - f1);
if (a < 0) a = -a;
edgeInterp[0][e0] = (float) a*mySamples[0][v0] + (1-a)*mySamples[0][v1];
edgeInterp[1][e0] = (float) a*mySamples[1][v0] + (1-a)*mySamples[1][v1];
edgeInterp[2][e0] = (float) a*mySamples[2][v0] + (1-a)*mySamples[2][v1];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e0] = (float) a*auxValues[j][v0] + (1-a)*auxValues[j][v1];
*/
}
globalToVertex[e0] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e3])) {
*/
// test for missing
if (edgeInterp[0][e3] != edgeInterp[0][e3]) {
float a = (isolevel - f2)/(f1 - f2);
if (a < 0) a = -a;
edgeInterp[0][e3] = (float) a*mySamples[0][v1] + (1-a)*mySamples[0][v2];
edgeInterp[1][e3] = (float) a*mySamples[1][v1] + (1-a)*mySamples[1][v2];
edgeInterp[2][e3] = (float) a*mySamples[2][v1] + (1-a)*mySamples[2][v2];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e3] = (float) a*auxValues[j][v1] + (1-a)*auxValues[j][v2];
*/
}
globalToVertex[e3] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e4])) {
*/
// test for missing
if (edgeInterp[0][e4] != edgeInterp[0][e4]) {
float a = (isolevel - f3)/(f1 - f3);
if (a < 0) a = -a;
edgeInterp[0][e4] = (float) a*mySamples[0][v1] + (1-a)*mySamples[0][v3];
edgeInterp[1][e4] = (float) a*mySamples[1][v1] + (1-a)*mySamples[1][v3];
edgeInterp[2][e4] = (float) a*mySamples[2][v1] + (1-a)*mySamples[2][v3];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e4] = (float) a*auxValues[j][v1] + (1-a)*auxValues[j][v3];
*/
}
globalToVertex[e4] = nvertex;
nvertex++;
}
// fill in the polys array
polys[npolygons][0] = e0;
if (index == 2) {
polys[npolygons][1] = e4;
polys[npolygons][2] = e3;
}
else { // index == 13
polys[npolygons][1] = e3;
polys[npolygons][2] = e4;
}
polys[npolygons][3] = -1;
// on to the next tetrahedron
npolygons++;
break;
case 3:
case 12: // plane slices a quadrilateral
// define edge values needed
e1 = Delan.Edges[i][1];
e2 = Delan.Edges[i][2];
e3 = Delan.Edges[i][3];
e4 = Delan.Edges[i][4];
// fill in edge interpolation values if they haven't been found
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e1])) {
*/
// test for missing
if (edgeInterp[0][e1] != edgeInterp[0][e1]) {
float a = (isolevel - f2)/(f0 - f2);
if (a < 0) a = -a;
edgeInterp[0][e1] = (float) a*mySamples[0][v0] + (1-a)*mySamples[0][v2];
edgeInterp[1][e1] = (float) a*mySamples[1][v0] + (1-a)*mySamples[1][v2];
edgeInterp[2][e1] = (float) a*mySamples[2][v0] + (1-a)*mySamples[2][v2];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e1] = (float) a*auxValues[j][v0] + (1-a)*auxValues[j][v2];
*/
}
globalToVertex[e1] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e2])) {
*/
// test for missing
if (edgeInterp[0][e2] != edgeInterp[0][e2]) {
float a = (isolevel - f3)/(f0 - f3);
if (a < 0) a = -a;
edgeInterp[0][e2] = (float) a*mySamples[0][v0] + (1-a)*mySamples[0][v3];
edgeInterp[1][e2] = (float) a*mySamples[1][v0] + (1-a)*mySamples[1][v3];
edgeInterp[2][e2] = (float) a*mySamples[2][v0] + (1-a)*mySamples[2][v3];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e2] = (float) a*auxValues[j][v0] + (1-a)*auxValues[j][v3];
*/
}
globalToVertex[e2] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e4])) {
*/
// test for missing
if (edgeInterp[0][e4] != edgeInterp[0][e4]) {
float a = (isolevel - f3)/(f1 - f3);
if (a < 0) a = -a;
edgeInterp[0][e4] = (float) a*mySamples[0][v1] + (1-a)*mySamples[0][v3];
edgeInterp[1][e4] = (float) a*mySamples[1][v1] + (1-a)*mySamples[1][v3];
edgeInterp[2][e4] = (float) a*mySamples[2][v1] + (1-a)*mySamples[2][v3];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e4] = (float) a*auxValues[j][v1] + (1-a)*auxValues[j][v3];
*/
}
globalToVertex[e4] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e3])) {
*/
// test for missing
if (edgeInterp[0][e3] != edgeInterp[0][e3]) {
float a = (isolevel - f2)/(f1 - f2);
if (a < 0) a = -a;
edgeInterp[0][e3] = (float) a*mySamples[0][v1] + (1-a)*mySamples[0][v2];
edgeInterp[1][e3] = (float) a*mySamples[1][v1] + (1-a)*mySamples[1][v2];
edgeInterp[2][e3] = (float) a*mySamples[2][v1] + (1-a)*mySamples[2][v2];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e3] = (float) a*auxValues[j][v1] + (1-a)*auxValues[j][v2];
*/
}
globalToVertex[e3] = nvertex;
nvertex++;
}
// fill in the polys array
polys[npolygons][0] = e1;
if (index == 3) {
polys[npolygons][1] = e2;
polys[npolygons][2] = e4;
polys[npolygons][3] = e3;
}
else { // index == 12
polys[npolygons][1] = e3;
polys[npolygons][2] = e4;
polys[npolygons][3] = e2;
}
// on to the next tetrahedron
npolygons++;
break;
case 4:
case 11: // plane slices a triangle
// define edge values needed
e1 = Delan.Edges[i][1];
e3 = Delan.Edges[i][3];
e5 = Delan.Edges[i][5];
// fill in edge interpolation values if they haven't been found
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e1])) {
*/
// test for missing
if (edgeInterp[0][e1] != edgeInterp[0][e1]) {
float a = (isolevel - f2)/(f0 - f2);
if (a < 0) a = -a;
edgeInterp[0][e1] = (float) a*mySamples[0][v0] + (1-a)*mySamples[0][v2];
edgeInterp[1][e1] = (float) a*mySamples[1][v0] + (1-a)*mySamples[1][v2];
edgeInterp[2][e1] = (float) a*mySamples[2][v0] + (1-a)*mySamples[2][v2];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e1] = (float) a*auxValues[j][v0] + (1-a)*auxValues[j][v2];
*/
}
globalToVertex[e1] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e3])) {
*/
// test for missing
if (edgeInterp[0][e3] != edgeInterp[0][e3]) {
float a = (isolevel - f2)/(f1 - f2);
if (a < 0) a = -a;
edgeInterp[0][e3] = (float) a*mySamples[0][v1] + (1-a)*mySamples[0][v2];
edgeInterp[1][e3] = (float) a*mySamples[1][v1] + (1-a)*mySamples[1][v2];
edgeInterp[2][e3] = (float) a*mySamples[2][v1] + (1-a)*mySamples[2][v2];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e3] = (float) a*auxValues[j][v1] + (1-a)*auxValues[j][v2];
*/
}
globalToVertex[e3] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e5])) {
*/
// test for missing
if (edgeInterp[0][e5] != edgeInterp[0][e5]) {
float a = (isolevel - f3)/(f2 - f3);
if (a < 0) a = -a;
edgeInterp[0][e5] = (float) a*mySamples[0][v2] + (1-a)*mySamples[0][v3];
edgeInterp[1][e5] = (float) a*mySamples[1][v2] + (1-a)*mySamples[1][v3];
edgeInterp[2][e5] = (float) a*mySamples[2][v2] + (1-a)*mySamples[2][v3];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e5] = (float) a*auxValues[j][v2] + (1-a)*auxValues[j][v3];
*/
}
globalToVertex[e5] = nvertex;
nvertex++;
}
// fill in the polys array
polys[npolygons][0] = e1;
if (index == 4) {
polys[npolygons][1] = e3;
polys[npolygons][2] = e5;
}
else { // index == 11
polys[npolygons][1] = e5;
polys[npolygons][2] = e3;
}
polys[npolygons][3] = -1;
// on to the next tetrahedron
npolygons++;
break;
case 5:
case 10: // plane slices a quadrilateral
// define edge values needed
e0 = Delan.Edges[i][0];
e2 = Delan.Edges[i][2];
e3 = Delan.Edges[i][3];
e5 = Delan.Edges[i][5];
// fill in edge interpolation values if they haven't been found
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e0])) {
*/
// test for missing
if (edgeInterp[0][e0] != edgeInterp[0][e0]) {
float a = (isolevel - f1)/(f0 - f1);
if (a < 0) a = -a;
edgeInterp[0][e0] = (float) a*mySamples[0][v0] + (1-a)*mySamples[0][v1];
edgeInterp[1][e0] = (float) a*mySamples[1][v0] + (1-a)*mySamples[1][v1];
edgeInterp[2][e0] = (float) a*mySamples[2][v0] + (1-a)*mySamples[2][v1];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e0] = (float) a*auxValues[j][v0] + (1-a)*auxValues[j][v1];
*/
}
globalToVertex[e0] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e2])) {
*/
// test for missing
if (edgeInterp[0][e2] != edgeInterp[0][e2]) {
float a = (isolevel - f3)/(f0 - f3);
if (a < 0) a = -a;
edgeInterp[0][e2] = (float) a*mySamples[0][v0] + (1-a)*mySamples[0][v3];
edgeInterp[1][e2] = (float) a*mySamples[1][v0] + (1-a)*mySamples[1][v3];
edgeInterp[2][e2] = (float) a*mySamples[2][v0] + (1-a)*mySamples[2][v3];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e2] = (float) a*auxValues[j][v0] + (1-a)*auxValues[j][v3];
*/
}
globalToVertex[e2] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e5])) {
*/
// test for missing
if (edgeInterp[0][e5] != edgeInterp[0][e5]) {
float a = (isolevel - f3)/(f2 - f3);
if (a < 0) a = -a;
edgeInterp[0][e5] = (float) a*mySamples[0][v2] + (1-a)*mySamples[0][v3];
edgeInterp[1][e5] = (float) a*mySamples[1][v2] + (1-a)*mySamples[1][v3];
edgeInterp[2][e5] = (float) a*mySamples[2][v2] + (1-a)*mySamples[2][v3];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e5] = (float) a*auxValues[j][v2] + (1-a)*auxValues[j][v3];
*/
}
globalToVertex[e5] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e3])) {
*/
// test for missing
if (edgeInterp[0][e3] != edgeInterp[0][e3]) {
float a = (isolevel - f2)/(f1 - f2);
if (a < 0) a = -a;
edgeInterp[0][e3] = (float) a*mySamples[0][v1] + (1-a)*mySamples[0][v2];
edgeInterp[1][e3] = (float) a*mySamples[1][v1] + (1-a)*mySamples[1][v2];
edgeInterp[2][e3] = (float) a*mySamples[2][v1] + (1-a)*mySamples[2][v2];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e3] = (float) a*auxValues[j][v1] + (1-a)*auxValues[j][v2];
*/
}
globalToVertex[e3] = nvertex;
nvertex++;
}
// fill in the polys array
polys[npolygons][0] = e0;
if (index == 5) {
polys[npolygons][1] = e3;
polys[npolygons][2] = e5;
polys[npolygons][3] = e2;
}
else { // index == 10
polys[npolygons][1] = e2;
polys[npolygons][2] = e5;
polys[npolygons][3] = e3;
}
// on to the next tetrahedron
npolygons++;
break;
case 6:
case 9: // plane slices a quadrilateral
// define edge values needed
e0 = Delan.Edges[i][0];
e1 = Delan.Edges[i][1];
e4 = Delan.Edges[i][4];
e5 = Delan.Edges[i][5];
// fill in edge interpolation values if they haven't been found
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e0])) {
*/
// test for missing
if (edgeInterp[0][e0] != edgeInterp[0][e0]) {
float a = (isolevel - f1)/(f0 - f1);
if (a < 0) a = -a;
edgeInterp[0][e0] = (float) a*mySamples[0][v0] + (1-a)*mySamples[0][v1];
edgeInterp[1][e0] = (float) a*mySamples[1][v0] + (1-a)*mySamples[1][v1];
edgeInterp[2][e0] = (float) a*mySamples[2][v0] + (1-a)*mySamples[2][v1];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e0] = (float) a*auxValues[j][v0] + (1-a)*auxValues[j][v1];
*/
}
globalToVertex[e0] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e1])) {
*/
// test for missing
if (edgeInterp[0][e1] != edgeInterp[0][e1]) {
float a = (isolevel - f2)/(f0 - f2);
if (a < 0) a = -a;
edgeInterp[0][e1] = (float) a*mySamples[0][v0] + (1-a)*mySamples[0][v2];
edgeInterp[1][e1] = (float) a*mySamples[1][v0] + (1-a)*mySamples[1][v2];
edgeInterp[2][e1] = (float) a*mySamples[2][v0] + (1-a)*mySamples[2][v2];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e1] = (float) a*auxValues[j][v0] + (1-a)*auxValues[j][v2];
*/
}
globalToVertex[e1] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e5])) {
*/
// test for missing
if (edgeInterp[0][e5] != edgeInterp[0][e5]) {
float a = (isolevel - f3)/(f2 - f3);
if (a < 0) a = -a;
edgeInterp[0][e5] = (float) a*mySamples[0][v2] + (1-a)*mySamples[0][v3];
edgeInterp[1][e5] = (float) a*mySamples[1][v2] + (1-a)*mySamples[1][v3];
edgeInterp[2][e5] = (float) a*mySamples[2][v2] + (1-a)*mySamples[2][v3];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e5] = (float) a*auxValues[j][v2] + (1-a)*auxValues[j][v3];
*/
}
globalToVertex[e5] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e4])) {
*/
// test for missing
if (edgeInterp[0][e4] != edgeInterp[0][e4]) {
float a = (isolevel - f3)/(f1 - f3);
if (a < 0) a = -a;
edgeInterp[0][e4] = (float) a*mySamples[0][v1] + (1-a)*mySamples[0][v3];
edgeInterp[1][e4] = (float) a*mySamples[1][v1] + (1-a)*mySamples[1][v3];
edgeInterp[2][e4] = (float) a*mySamples[2][v1] + (1-a)*mySamples[2][v3];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e4] = (float) a*auxValues[j][v1] + (1-a)*auxValues[j][v3];
*/
}
globalToVertex[e4] = nvertex;
nvertex++;
}
// fill in the polys array
polys[npolygons][0] = e0;
if (index == 6) {
polys[npolygons][1] = e4;
polys[npolygons][2] = e5;
polys[npolygons][3] = e1;
}
else { // index == 9
polys[npolygons][1] = e1;
polys[npolygons][2] = e5;
polys[npolygons][3] = e4;
}
// on to the next tetrahedron
npolygons++;
break;
case 7:
case 8: // plane slices a triangle
// interpolate between 3:0, 3:1, 3:2 for tri edges, same as case 8
// define edge values needed
e2 = Delan.Edges[i][2];
e4 = Delan.Edges[i][4];
e5 = Delan.Edges[i][5];
// fill in edge interpolation values if they haven't been found
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e2])) {
*/
// test for missing
if (edgeInterp[0][e2] != edgeInterp[0][e2]) {
float a = (isolevel - f3)/(f0 - f3);
if (a < 0) a = -a;
edgeInterp[0][e2] = (float) a*mySamples[0][v0] + (1-a)*mySamples[0][v3];
edgeInterp[1][e2] = (float) a*mySamples[1][v0] + (1-a)*mySamples[1][v3];
edgeInterp[2][e2] = (float) a*mySamples[2][v0] + (1-a)*mySamples[2][v3];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e2] = (float) a*auxValues[j][v0] + (1-a)*auxValues[j][v3];
*/
}
globalToVertex[e2] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e4])) {
*/
// test for missing
if (edgeInterp[0][e4] != edgeInterp[0][e4]) {
float a = (isolevel - f3)/(f1 - f3);
if (a < 0) a = -a;
edgeInterp[0][e4] = (float) a*mySamples[0][v1] + (1-a)*mySamples[0][v3];
edgeInterp[1][e4] = (float) a*mySamples[1][v1] + (1-a)*mySamples[1][v3];
edgeInterp[2][e4] = (float) a*mySamples[2][v1] + (1-a)*mySamples[2][v3];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e4] = (float) a*auxValues[j][v1] + (1-a)*auxValues[j][v3];
*/
}
globalToVertex[e4] = nvertex;
nvertex++;
}
/* WLH 25 Oct 97
if (Double.isNaN(edgeInterp[0][e5])) {
*/
// test for missing
if (edgeInterp[0][e5] != edgeInterp[0][e5]) {
float a = (isolevel - f3)/(f2 - f3);
if (a < 0) a = -a;
edgeInterp[0][e5] = (float) a*mySamples[0][v2] + (1-a)*mySamples[0][v3];
edgeInterp[1][e5] = (float) a*mySamples[1][v2] + (1-a)*mySamples[1][v3];
edgeInterp[2][e5] = (float) a*mySamples[2][v2] + (1-a)*mySamples[2][v3];
for (int j=0; j 255) ? -1 : ((t < 128) ? t : t - 256) ) );
/* MEM_WLH
auxInterp[j][e5] = (float) a*auxValues[j][v2] + (1-a)*auxValues[j][v3];
*/
}
globalToVertex[e5] = nvertex;
nvertex++;
}
// fill in the polys array
polys[npolygons][0] = e2;
if (index == 7) {
polys[npolygons][1] = e4;
polys[npolygons][2] = e5;
}
else { // index == 8
polys[npolygons][1] = e5;
polys[npolygons][2] = e4;
}
polys[npolygons][3] = -1;
// on to the next tetrahedron
npolygons++;
break;
} // end switch (index)
} // end for (int i=0; i global edges) " + npolygons + "\n");
for (int i=0; i ";
for (int j=0; j<4; j++) {
s = s + " " + polys[i][j];
}
System.out.println(s + "\n");
}
}
// transform polys array into polyToVert array
polyToVert[0] = new int[npolygons][];
for (int i=0; i vertices) " + npolygons + "\n");
for (int i=0; i ";
for (int j=0; j 3) nverts[polyToVert[0][i][3]]++;
}
// initialize vertToPoly array
vertToPoly[0] = new int[nvertex][];
for (int i=0; i 3) {
int d = polyToVert[0][i][3];
if (d != -1) vertToPoly[0][d][nverts[d]++] = i;
}
}
if (DEBUG) {
System.out.println("\nvertToPoly (vertices -> polys) " + nvertex + "\n");
for (int i=0; i ";
for (int j=0; j= 0) {
fieldVertices[0][k] = edgeInterp[0][i];
fieldVertices[1][k] = edgeInterp[1][i];
fieldVertices[2][k] = edgeInterp[2][i];
for (int j=0; j ";
for (int j=0; j<3; j++) {
s = s + " " + fieldVertices[j][i];
}
System.out.println(s + "\n");
}
}
}
/*
make_normals and poly_triangle_stripe altered according to:
Pol_f_Vert[9*i + 8] --> vertToPoly[i].length
Pol_f_Vert[9*i + off] --> vertToPoly[i][off]
Vert_f_Pol[7*i + 6] --> polyToVert[i].length
Vert_f_Pol[7*i + off] --> polyToVert[i][off]
*/
/* from Contour3D.java, used by make_normals */
static final float EPS_0 = (float) 1.0e-5;
/* copied from Contour3D.java */
private static void make_normals(float[] VX, float[] VY, float[] VZ,
float[] NX, float[] NY, float[] NZ, int nvertex,
int npolygons, float[] Pnx, float[] Pny, float[] Pnz,
float[] NxA, float[] NxB, float[] NyA, float[] NyB,
float[] NzA, float[] NzB,
int[][] vertToPoly, int[][] polyToVert)
/* WLH 25 Oct 97
int[] Pol_f_Vert, int[] Vert_f_Pol )
*/
throws VisADException {
int i, k, n;
int i1, i2, ix, iy, iz, ixb, iyb, izb;
int max_vert_per_pol, swap_flag;
float x, y, z, a, minimum_area, len;
int iv[] = new int[3];
if (nvertex <= 0) return;
for (i = 0; i < nvertex; i++) {
NX[i] = 0;
NY[i] = 0;
NZ[i] = 0;
}
// WLH 12 Nov 2001
// minimum_area = (float) ((1.e-4 > EPS_0) ? 1.e-4 : EPS_0);
minimum_area = Float.MIN_VALUE;
/* Calculate maximum number of vertices per polygon */
/* WLH 25 Oct 97
k = 6;
n = 7*npolygons;
*/
k = 0;
while ( true ) {
/* WLH 25 Oct 97
for (i=k+7; i Vert_f_Pol[k]) break;
*/
for (i=k+1; i polyToVert[k].length) break;
/* WLH 25 Oct 97
if (i >= n) break;
*/
if (i >= npolygons) break;
k = i;
}
/* WLH 25 Oct 97
max_vert_per_pol = Vert_f_Pol[k];
*/
max_vert_per_pol = polyToVert[k].length;
/* Calculate the Normals vector components for each Polygon */
for ( i=0; i0) {
NxA[i] = VX[Vert_f_Pol[1+i*7]] - VX[Vert_f_Pol[0+i*7]];
NyA[i] = VY[Vert_f_Pol[1+i*7]] - VY[Vert_f_Pol[0+i*7]];
NzA[i] = VZ[Vert_f_Pol[1+i*7]] - VZ[Vert_f_Pol[0+i*7]];
}
*/
if (polyToVert[i].length>0) {
int j1 = polyToVert[i][1];
int j0 = polyToVert[i][0];
NxA[i] = VX[j1] - VX[j0];
NyA[i] = VY[j1] - VY[j0];
NzA[i] = VZ[j1] - VZ[j0];
}
}
swap_flag = 0;
for ( k = 2; k < max_vert_per_pol; k++ ) {
if (swap_flag==0) {
/*$dir no_recurrence */ /* Vectorized */
for ( i=0; i= 0 ) {
NxB[i] = VX[Vert_f_Pol[k+i*7]] - VX[Vert_f_Pol[0+i*7]];
NyB[i] = VY[Vert_f_Pol[k+i*7]] - VY[Vert_f_Pol[0+i*7]];
NzB[i] = VZ[Vert_f_Pol[k+i*7]] - VZ[Vert_f_Pol[0+i*7]];
*/
if ( k < polyToVert[i].length ) {
int jk = polyToVert[i][k];
int j0 = polyToVert[i][0];
NxB[i] = VX[jk] - VX[j0];
NyB[i] = VY[jk] - VY[j0];
NzB[i] = VZ[jk] - VZ[j0];
Pnx[i] = NyA[i]*NzB[i] - NzA[i]*NyB[i];
Pny[i] = NzA[i]*NxB[i] - NxA[i]*NzB[i];
Pnz[i] = NxA[i]*NyB[i] - NyA[i]*NxB[i];
NxA[i] = Pnx[i]*Pnx[i] + Pny[i]*Pny[i] + Pnz[i]*Pnz[i];
if (NxA[i] > minimum_area) {
Pnx[i] /= NxA[i];
Pny[i] /= NxA[i];
Pnz[i] /= NxA[i];
}
}
}
}
else { /* swap_flag!=0 */
/*$dir no_recurrence */ /* Vectorized */
for (i=0; i= 0 ) {
NxA[i] = VX[Vert_f_Pol[k+i*7]] - VX[Vert_f_Pol[0+i*7]];
NyA[i] = VY[Vert_f_Pol[k+i*7]] - VY[Vert_f_Pol[0+i*7]];
NzA[i] = VZ[Vert_f_Pol[k+i*7]] - VZ[Vert_f_Pol[0+i*7]];
*/
if ( k < polyToVert[i].length ) {
int jk = polyToVert[i][k];
int j0 = polyToVert[i][0];
NxA[i] = VX[jk] - VX[j0];
NyA[i] = VY[jk] - VY[j0];
NzA[i] = VZ[jk] - VZ[j0];
Pnx[i] = NyB[i]*NzA[i] - NzB[i]*NyA[i];
Pny[i] = NzB[i]*NxA[i] - NxB[i]*NzA[i];
Pnz[i] = NxB[i]*NyA[i] - NyB[i]*NxA[i];
NxB[i] = Pnx[i]*Pnx[i] + Pny[i]*Pny[i] + Pnz[i]*Pnz[i];
if (NxB[i] > minimum_area) {
Pnx[i] /= NxB[i];
Pny[i] /= NxB[i];
Pnz[i] /= NxB[i];
}
}
} // end for (i=0; i be Vectorized */
for ( i=0; i= 0) {
iv[0] = Vert_f_Pol[0+i*7];
iv[1] = Vert_f_Pol[(k-1)+i*7];
iv[2] = Vert_f_Pol[k+i*7];
*/
if (k < polyToVert[i].length) {
iv[0] = polyToVert[i][0];
iv[1] = polyToVert[i][k-1];
iv[2] = polyToVert[i][k];
x = Pnx[i]; y = Pny[i]; z = Pnz[i];
/*
System.out.println("vertices: " + iv[0] + " " + iv[1] + " " + iv[2]);
System.out.println(" normal: " + x + " " + y + " " + z + "\n");
*/
// Update the origin vertex
NX[iv[0]] += x; NY[iv[0]] += y; NZ[iv[0]] += z;
// Update the vertex that defines the first vector
NX[iv[1]] += x; NY[iv[1]] += y; NZ[iv[1]] += z;
// Update the vertex that defines the second vector
NX[iv[2]] += x; NY[iv[2]] += y; NZ[iv[2]] += z;
}
} // end for ( i=0; i EPS_0) {
NX[i] /= len;
NY[i] /= len;
NZ[i] /= len;
}
}
}
/* from Contour3D.java, used by poly_triangle_stripe */
static final int NTAB[] =
{ 0,1,2, 1,2,0, 2,0,1,
0,1,3,2, 1,2,0,3, 2,3,1,0, 3,0,2,1,
0,1,4,2,3, 1,2,0,3,4, 2,3,1,4,0, 3,4,2,0,1, 4,0,3,1,2,
0,1,5,2,4,3, 1,2,0,3,5,4, 2,3,1,4,0,5, 3,4,2,5,1,0, 4,5,3,0,2,1,
5,0,4,1,3,2
};
/* from Contour3D.java, used by poly_triangle_stripe */
static final int ITAB[] =
{ 0,2,1, 1,0,2, 2,1,0,
0,3,1,2, 1,0,2,3, 2,1,3,0, 3,2,0,1,
0,4,1,3,2, 1,0,2,4,3, 2,1,3,0,4, 3,2,4,1,0, 4,3,0,2,1,
0,5,1,4,2,3, 1,0,2,5,3,4, 2,1,3,0,4,5, 3,2,4,1,5,0, 4,3,5,2,0,1,
5,4,0,3,1,2
};
/* from Contour3D.java, used by poly_triangle_stripe */
static final int STAB[] = { 0, 9, 25, 50 };
/* copied from Contour3D.java */
static int poly_triangle_stripe( int[] Tri_Stripe,
int nvertex, int npolygons,
int[][] vertToPoly, int[][] polyToVert)
/* WLH 25 Oct 97
int[] Pol_f_Vert, int[] Vert_f_Pol )
*/
throws VisADException {
int i, j, k, m, ii, npol, cpol, idx, off, Nvt,
vA, vB, ivA, ivB, iST, last_pol;
boolean f_line_conection = false;
boolean[] vet_pol = new boolean[npolygons];
last_pol = 0;
npol = 0;
iST = 0;
ivB = 0;
for (i=0; i=0 && Vert_f_Pol[ivB+off]>=0) {
i=Vert_f_Pol[ivA+off]*9;
k=i+Pol_f_Vert [i+8];
j=Vert_f_Pol[ivB+off]*9;
m=j+Pol_f_Vert [j+8];
while (i>0 && j>0 && i= 0) break;
}
/* insert polygon alone */
if (npol < 0)
{ /*ptT = NTAB + STAB[Nvt-3];*/
idx = STAB[Nvt-3];
if (iST > 0)
{ Tri_Stripe[iST] = Tri_Stripe[iST-1]; iST++;
/* WLH 25 Oct 97
Tri_Stripe[iST++] = Vert_f_Pol[NTAB[idx]+off];
*/
// System.out.println("1 Tri_Stripe[" + iST + "] from poly " + cpol);
Tri_Stripe[iST++] = polyToVert[cpol][NTAB[idx]];
}
else f_line_conection = true; /* WLH 3-9-95 added */
for (ii=0; ii< ((Nvt < 6) ? Nvt:6); ii++) {
/* WLH 25 Oct 97
Tri_Stripe[iST++] = Vert_f_Pol[NTAB[idx++]+off];
*/
// System.out.println("2 Tri_Stripe[" + iST + "] from poly " + cpol);
Tri_Stripe[iST++] = polyToVert[cpol][NTAB[idx++]];
}
continue;
}
if (( (ivB != 0) && ivA==(ivB-1)) || ( !(ivB != 0) && ivA==Nvt-1)) {
/* ptT = ITAB + STAB[Nvt-3] + (ivB+1)*Nvt; */
idx = STAB[Nvt-3] + (ivB+1)*Nvt;
if (f_line_conection)
{ Tri_Stripe[iST] = Tri_Stripe[iST-1]; iST++;
/* WLH 25 Oct 97
Tri_Stripe[iST++] = Vert_f_Pol[ITAB[idx-1]+off];
*/
// System.out.println("3 Tri_Stripe[" + iST + "] from poly " + cpol);
Tri_Stripe[iST++] = polyToVert[cpol][ITAB[idx-1]];
f_line_conection = false;
}
for (ii=0; ii<((Nvt < 6) ? Nvt:6); ii++) {
/* WLH 25 Oct 97
Tri_Stripe[iST++] = Vert_f_Pol[ITAB[--idx]+off];
*/
// System.out.println("4 Tri_Stripe[" + iST + "] from poly " + cpol);
Tri_Stripe[iST++] = polyToVert[cpol][ITAB[--idx]];
}
}
else {
/* ptT = NTAB + STAB[Nvt-3] + (ivB+1)*Nvt; */
idx = STAB[Nvt-3] + (ivB+1)*Nvt;
if (f_line_conection)
{ Tri_Stripe[iST] = Tri_Stripe[iST-1]; iST++;
/* WLH 25 Oct 97
Tri_Stripe[iST++] = Vert_f_Pol[NTAB[idx-1]+off];
*/
// System.out.println("5 Tri_Stripe[" + iST + "] from poly " + cpol);
Tri_Stripe[iST++] = polyToVert[cpol][NTAB[idx-1]];
f_line_conection = false;
}
for (ii=0; ii<((Nvt < 6) ? Nvt:6); ii++) {
/* WLH 25 Oct 97
Tri_Stripe[iST++] = Vert_f_Pol[NTAB[--idx]+off];
*/
// System.out.println("6 Tri_Stripe[" + iST + "] from poly " + cpol);
Tri_Stripe[iST++] = polyToVert[cpol][NTAB[--idx]];
}
}
vB = Tri_Stripe[iST-1];
vA = Tri_Stripe[iST-2];
cpol = npol;
while (true)
{
/* get_vertices_of_pol(cpol,Vt,Nvt) { */
/* WLH 25 Oct 97
Nvt = Vert_f_Pol [(j=cpol*7)+6];
off = j;
*/
Nvt = polyToVert[cpol].length;
/* } */
/* update_polygon(cpol) */
// System.out.println("2 vet_pol[" + cpol + "] = false");
vet_pol[cpol] = false;
/* WLH 25 Oct 97
for (ivA=0; ivA=0 && vB>=0) {
i=vA*9;
k=i+Pol_f_Vert [i+8];
j=vB*9;
m=j+Pol_f_Vert [j+8];
while (i>0 && j>0 && i=0 && vB>=0) {
int ilim = vertToPoly[vA].length;
k = 0;
int jlim = vertToPoly[vB].length;
m = 0;
// while (0=0 && vB>=0) {
i=vA*9;
k=i+Pol_f_Vert [i+8];
j=vB*9;
m=j+Pol_f_Vert [j+8];
while (i>0 && j>0 && i=0 && vB>=0) {
int ilim = vertToPoly[vA].length;
k = 0;
int jlim = vertToPoly[vB].length;
m = 0;
// while (0= 0 && tri[j] == null) {
// compare signatures of neighbors
int v1 = Tri[i][w];
int v2 = Tri[i][(w + 1) % 3];
int i1 = -1;
int i2 = -1;
int j1 = -1;
int j2 = -1;
for (int k=0; k<3; k++) {
if (tri[i][k] == v1) i1 = k;
if (tri[i][k] == v2) i2 = k;
if (Tri[j][k] == v1) j1 = k;
if (Tri[j][k] == v2) j2 = k;
}
tri[j] = new int[3];
tri[j][0] = Tri[j][0];
if ( ( (((i1 + 1) % 3) == i2) && (((j1 + 1) % 3) == j2) ) ||
( (((i2 + 1) % 3) == i1) && (((j2 + 1) % 3) == j1) ) ) {
tri[j][1] = Tri[j][2];
tri[j][2] = Tri[j][1];
}
else {
tri[j][1] = Tri[j][1];
tri[j][2] = Tri[j][2];
}
// add j to stack
sp++;
walk_stack[sp] = 0;
poly_stack[sp] = j;
}
else { // (j < 0 || tri[j] != null)
while (true) {
walk_stack[sp]++;
if (walk_stack[sp] < 3) {
break;
}
else {
sp--;
if (sp < 0) break;
}
} // end while (true)
} // end if (j < 0 || tri[j] != null)
if (sp < 0) break;
} // end while (true)
} // end if (tri[ii] == null)
} // end for (int ii=0; ii #"+index[i]);
}
// test valueToInterp function
System.out.println("\nvalueToInterp test:");
int[][] indices = new int[value[0].length][];
float[][] weights = new float[value[0].length][];
iSet3D.valueToInterp(value, indices, weights);
for (int i=0; i ["
+indices[i][0]+", "
+indices[i][1]+", "
+indices[i][2]+", "
+indices[i][3]+"]\tweight total: "
+(weights[i][0]+weights[i][1]
+weights[i][2]+weights[i][3]));
}
// test makeIsosurface function
System.out.println("\nmakeIsosurface test:");
float[] field = {100, 300, 320, 250, 80, 70, 135, 110, 105};
float[][] slice = new float[3][];
int[][][] polyvert = new int[1][][];
int[][][] vertpoly = new int[1][][];
iSet3D.makeIsosurface(288, field, null, slice, null, polyvert, vertpoly);
for (int i=0; i #0
221.0, 161.0, 300.0 --> #2
319.0, 158.0, 127.0 --> #5
215.0, 139.0, 287.0 --> #2
196.0, 285.0, 194.0 --> #6
valueToInterp test:
189.0, 166.0, 207.0 --> [0, 1, 2, 4] weight total: 1.0
221.0, 161.0, 300.0 --> [1, 2, 3, 8] weight total: 1.0
319.0, 158.0, 127.0 --> [4, 5, 6, 8] weight total: 0.9999999999999999
215.0, 139.0, 287.0 --> [1, 2, 3, 8] weight total: 1.0
196.0, 285.0, 194.0 --> [1, 5, 6, 7] weight total: 1.0
makeIsosurface test:
polygons:
#0: 237.843, 226.93, 291.817
(tri) 227.2, 236.6, 292.709
236.547, 236.937, 288.368
#1: 228.82, 222.3, 290.2
(quad) 182.418, 142.4, 313.273
225.127, 228.382, 291.291
172.733, 156.133, 316.267
#2: 225.127, 228.382, 291.291
(tri) 227.2, 236.6, 292.709
235.323, 222.262, 291.215
#3: 228.82, 222.3, 290.2
(quad) 182.418, 142.4, 313.273
235.323, 222.262, 291.215
198.33, 142.623, 315.637
#4: 234.88, 205.08, 313.24
(tri) 237.843, 226.93, 291.817
235.323, 222.262, 291.215
#5: 182.418, 142.4, 313.273
(tri) 210.886, 138.743, 348.743
198.33, 142.623, 315.637
#6: 234.88, 205.08, 313.24
(quad) 235.323, 222.262, 291.215
210.886, 138.743, 348.743
198.33, 142.623, 315.637
#7: 225.127, 228.382, 291.291
(quad) 227.2, 236.6, 292.709
172.733, 156.133, 316.267
180.059, 178.984, 318.497
#8: 234.88, 205.08, 313.24
(tri) 237.843, 226.93, 291.817
236.547, 236.937, 288.368
#9: 228.82, 222.3, 290.2
(tri) 225.127, 228.382, 291.291
235.323, 222.262, 291.215
#10: 237.843, 226.93, 291.817
(tri) 227.2, 236.6, 292.709
235.323, 222.262, 291.215
iris 46%
*/
}
Ads help maintain this website.