/*
!=======================================================================
!
!  PROGRAM  PHASE-Viewer  (PHASE-Viewer 2014.01 ver.3.3.0)
!
!  Created on ----
!  AUTHOR(S): KOGA, Junichiro
!  File : MarchingCube.java
!  
!  Contact address :  Phase System Consortium
!                     E-mail: phase_system@nims.go.jp URL https://azuma.nims.go.jp
!
!
!   Since 2002, this program set had been intensively developed as a part of the following 
!  national projects supported by the Ministry of Education, Culture, Sports, Science and
!  Technology (MEXT) of Japan; "Frontier Simulation Software for Industrial Science
!  (FSIS)" from 2002 to 2005, "Revolutionary Simulation Software (RSS21)" from 2006 to
!  2008. "Research and Development of Innovative Simulation Software (RISS)" from 2008
!  to 2013. These projects is lead by the Center for Research on Innovative Simulation 
!  Software (CISS), the Institute of Industrial Science (IIS), the University of Tokyo.
!   Since 2013, this program set has been further developed centering on PHASE System
!  Consortium. 
!   The activity of development of this program set has been supervised by Takahisa Ohno.
!
!=======================================================================
 */

//======================================================================
//                                                                     #
// Software Name : BioStationViewer Version 4.00                       #
//                                                                     #
//      Class Name :  This file name                                   #
//                                                                     #
//                             Written by Akifumi Kato  2003/11/07     #
//                                                                     #
//      Contact address :  The University of Tokyo, FSIS project       #
//                                                                     #
//      "Quantum Molecular Inreraction Analysis"  Group.               #
//                                                                     #
//======================================================================

package ciss.phase_viewer.acviewer.scenegraphelements.abinitmp;

import java.util.Vector;

import javax.vecmath.Point3f;
import javax.vecmath.Vector3f;

import ciss.phase_viewer.atomcoord.VolumetricData;

/****************************************************************************
 * This Program original from AVS module
 * 
 * avsのモジュールのソースを元に等値面をマーティンキューブ法で 求めるJavaソース
 * 
 * 入力 x,y,zの座標 格子点上の値 求める等値面の値
 * 
 * 出力 等値面を構成する三角ポリゴンの座標
 * 
 * 使用方法 MarchingCube.create_isosurf(x,y,z,data,value,tri_pos,normal) float
 * x[],y[],z[],data[],value Vector tri_pos; Point3fのVectorでポリゴンの座標を返す。 Vector
 * normal Vector3fのVectorでポリゴンの頂点の法線を返す。
 * 
 ****************************************************************************/

/*------------------------------------------------------------------
 The algorithm to create isosurface from 3D uniform data volume
 is a marching cubes method from Lorensen et al
 The main routine to generate isosurface is "create_isosurf".
 Input : 3D field scalar byte uniform
 Output: Geom object for rendering routine
 --------------------------------------------------------------------*/

/*
 * define case # of pattern from Marching Cube method the file must be included
 * in a main AVS routine for fast cube intersections retrieval by looking up
 * into the well-defined pattern table
 * 
 * Date: 12/10/92
 */

public class MarchingCube {

    /**
     * 補間方法の指定
     */
    public static final int LOG = 0;
    public static final int LINEAR = 1;

    /*
     * define initial case table for Rotation operation
     * 
     * entry format per word [pattern no][8 corner node indices] 6 bits 8x3=24
     * bits a initial pattern table with size equal to 128 must be initialized
     * with initial values from the following reference pattern structures
     */
    private static final byte PLYEDGES1[][] = { { 1, 3, 6 } };
    private static final byte PLYEDGES2[][] = { { 1, 2, 6 }, { 2, 6, 7 } };
    private static final byte PLYEDGES3[][] = { { 1, 3, 6 }, { 0, 2, 5 } };
    private static final byte PLYEDGES4[][] = { { 1, 3, 6 }, { 5, 8, 10 } };
    private static final byte PLYEDGES5[][] = { { 2, 6, 9 }, { 2, 3, 6 },
            { 2, 9, 10 } };
    private static final byte PLYEDGES6[][] = { { 5, 8, 10 }, { 1, 2, 6 },
            { 2, 6, 7 } };
    private static final byte PLYEDGES7[][] = { { 0, 1, 4 }, { 5, 8, 10 },
            { 2, 3, 7 } };
    private static final byte PLYEDGES8[][] = { { 1, 9, 10 }, { 1, 2, 10 } };
    private static final byte PLYEDGES9[][] = { { 1, 4, 8 }, { 1, 3, 8 },
            { 3, 8, 10 }, { 3, 7, 10 } };
    private static final byte PLYEDGES10[][] = { { 0, 4, 6 }, { 0, 3, 6 },
            { 5, 8, 11 }, { 5, 7, 11 } };
    private static final byte PLYEDGES11[][] = { { 1, 3, 9 }, { 5, 8, 9 },
            { 3, 5, 9 }, { 3, 5, 7 } };
    private static final byte PLYEDGES12[][] = { { 0, 1, 4 }, { 2, 6, 9 },
            { 2, 9, 10 }, { 2, 3, 6 } };
    private static final byte PLYEDGES13[][] = { { 4, 8, 9 }, { 1, 3, 6 },
            { 0, 2, 5 }, { 7, 10, 11 } };
    private static final byte PLYEDGES14[][] = { { 3, 4, 6 }, { 3, 4, 10 },
            { 4, 8, 10 }, { 2, 3, 10 } };

    /* POLYGON NUMBER */
    private static final int PATNO = 256;
    private static final int PlyNo1 = PLYEDGES1.length;
    private static final int PlyNo2 = PLYEDGES2.length;
    private static final int PlyNo3 = PLYEDGES3.length;
    private static final int PlyNo4 = PLYEDGES4.length;
    private static final int PlyNo5 = PLYEDGES5.length;
    private static final int PlyNo6 = PLYEDGES6.length;
    private static final int PlyNo7 = PLYEDGES7.length;
    private static final int PlyNo8 = PLYEDGES8.length;
    private static final int PlyNo9 = PLYEDGES9.length;
    private static final int PlyNo10 = PLYEDGES10.length;
    private static final int PlyNo11 = PLYEDGES11.length;
    private static final int PlyNo12 = PLYEDGES12.length;
    private static final int PlyNo13 = PLYEDGES13.length;
    private static final int PlyNo14 = PLYEDGES14.length;

    /* CROSSEDGE table defines the bit pattern which shows edge crossing */
    private static final int CROSSEDGE[] = { 00000, 00112, 00306, 00157, 02552,
            03114, 02746, 02677, 03006, 02632, 04771, 01652, 03137, 07777,
            02534 };
    private static final int EDGETABLE[][] = { { -1, 0, 1, -1, 4, -1, -1, -1 },
            { 0, -1, -1, 2, -1, 5, -1, -1 }, { 1, -1, -1, 3, -1, -1, 6, -1 },
            { -1, 2, 3, -1, -1, -1, -1, 7 }, { 4, -1, -1, -1, -1, 8, 9, -1 },
            { -1, 5, -1, -1, 8, -1, -1, 10 }, { -1, -1, 6, -1, 9, -1, -1, 11 },
            { -1, -1, -1, 7, -1, 10, 11, -1 } };
    private static final int UEDGETABLE[][] = { { 0, 1 }, { 0, 2 }, { 1, 3 },
            { 2, 3 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }, { 4, 5 },
            { 4, 6 }, { 5, 7 }, { 6, 7 } };

    /*
     * class PATSTRUC { int cubeindex; int patcode;
     * 
     * PATSTRUC(int c, int p) { cubeindex = c; patcode = p; } };
     * 
     * private static PATSTRUC PAT0[] = { new PATSTRUC(0x00,00076543210)};
     */

    /* PATTERN 0 */
    private static int PAT0[][] = { { 0x00, 00076543210 } };

    /* PATTERN 1 */
    private static int PAT1[][] = { { 0x04, 00176543210 },
            { 0x01, 00164752031 }, { 0x02, 00173625140 },
            { 0x08, 00157461302 }, { 0x10, 00126370415 },
            { 0x20, 00131207564 }, { 0x80, 00123016745 }, { 0x40, 00137152604 } };

    /* PATTERN 2 */
    private static int PAT2[][] = { { 0x0c, 00276543210 },
            { 0x44, 00237152604 }, { 0xc0, 00223016745 },
            { 0x88, 00262407351 }, { 0x05, 00264752031 },
            { 0x0a, 00257461302 }, { 0xa0, 00246025713 },
            { 0x50, 00275316420 }, { 0x03, 00232761054 },
            { 0x22, 00273625140 }, { 0x30, 00210325476 }, { 0x11, 00226370415 } };

    /* PATTERN 3 */
    private static int PAT3[][] = { { 0x06, 00376543210 },
            { 0x09, 00364752031 }, { 0x28, 00362407351 },
            { 0x82, 00320643175 }, { 0x41, 00337152604 },
            { 0x14, 00313570246 }, { 0x60, 00331207564 },
            { 0x90, 00323016745 }, { 0x12, 00326370415 },
            { 0x21, 00332761054 }, { 0x84, 00340516273 }, { 0x48, 00354107632 } };

    /* PATTERN 4 */
    private static int PAT4[][] = { { 0x24, 00476543210 },
            { 0x18, 00462407351 }, { 0x81, 00423016745 }, { 0x42, 00437152604 } };

    /* PATTERN 5 */
    private static int PAT5[][] = { { 0xc8, 00576543210 },
            { 0x8c, 00537152604 }, { 0x4c, 00523016745 },
            { 0xc4, 00562407351 }, { 0x0b, 00513570246 },
            { 0x0e, 00532761054 }, { 0x0d, 00520643175 },
            { 0x07, 00501452367 }, { 0x8a, 00531207564 },
            { 0x2a, 00515043726 }, { 0xa2, 00557461302 },
            { 0xa8, 00573625140 }, { 0x15, 00502134657 },
            { 0x51, 00540516273 }, { 0x54, 00564752031 },
            { 0x45, 00526370415 }, { 0x13, 00504261537 },
            { 0x31, 00545670123 }, { 0x32, 00551734062 },
            { 0x23, 00510325476 }, { 0xd0, 00567234501 },
            { 0x70, 00546025713 }, { 0xb0, 00554107632 }, { 0xe0, 00575316420 } };

    /* PATTERN 6 */
    private static int PAT6[][] = { { 0x2c, 00676543210 },
            { 0x46, 00637152604 }, { 0xc1, 00623016745 },
            { 0x98, 00662407351 }, { 0xa1, 00646025713 },
            { 0x26, 00604261537 }, { 0x4a, 00620643175 },
            { 0x85, 00664752031 }, { 0x43, 00645670123 },
            { 0x58, 00675316420 }, { 0x38, 00610325476 },
            { 0x91, 00651734062 }, { 0x1c, 00601452367 },
            { 0x64, 00640516273 }, { 0xc2, 00654107632 },
            { 0x89, 00615043726 }, { 0xa4, 00631207564 },
            { 0x25, 00613570246 }, { 0x83, 00632761054 },
            { 0x52, 00602134657 }, { 0x34, 00667234501 },
            { 0x19, 00626370415 }, { 0x62, 00673625140 }, { 0x1a, 00657461302 } };

    /* PATTERN 7 */
    private static int PAT7[][] = { { 0x29, 00776543210 },
            { 0x92, 00762407351 }, { 0x86, 00754107632 },
            { 0x68, 00740516273 }, { 0x94, 00701452367 },
            { 0x49, 00715043726 }, { 0x61, 00773625140 }, { 0x16, 00737152604 } };

    /* PATTERN 8 */
    private static int PAT8[][] = { { 0xcc, 01076543210 },
            { 0x0f, 01032761054 }, { 0x33, 01010325476 },
            { 0xf0, 01054107632 }, { 0xaa, 01031207564 }, { 0x55, 01026370415 } };

    /* PATTERN 9 */
    private static int PAT9[][] = { { 0xd4, 01176543210 },
            { 0xe8, 01137152604 }, { 0x8e, 01123016745 },
            { 0x4d, 01162407351 }, { 0x71, 01154107632 },
            { 0xb2, 01115043726 }, { 0x2b, 01101452367 }, { 0x17, 01140516273 } };

    /* PATTERN 10 */
    private static int PAT10[][] = { { 0xa5, 01276543210 },
            { 0x5a, 01237152604 }, { 0x66, 01254107632 },
            { 0xc3, 01215043726 }, { 0x99, 01201452367 }, { 0x3c, 01240516273 } };

    /* PATTERN 11 */
    private static int PAT11[][] = { { 0xe4, 01376543210 },
            { 0xca, 01337152604 }, { 0x8d, 01323016745 },
            { 0x5c, 01362407351 }, { 0x72, 01354107632 },
            { 0xa3, 01315043726 }, { 0x1b, 01301452367 },
            { 0x35, 01340516273 }, { 0x2e, 01331207564 },
            { 0x47, 01302134657 }, { 0xd1, 01346025713 }, { 0xb8, 01357461302 } };

    /* PATTERN 12 */
    private static int PAT12[][] = { { 0xc9, 01476543210 },
            { 0x9c, 01437152604 }, { 0x6c, 01423016745 },
            { 0xc6, 01462407351 }, { 0x56, 01464752031 },
            { 0x65, 01426370415 }, { 0x95, 01402134657 },
            { 0x59, 01440516273 }, { 0x1e, 01432761054 },
            { 0x2d, 01420643175 }, { 0x87, 01401452367 },
            { 0x4b, 01413570246 }, { 0xa9, 01473625140 },
            { 0x9a, 01431207564 }, { 0x6a, 01415043726 },
            { 0xa6, 01457461302 }, { 0xd2, 01467234501 },
            { 0x78, 01446025713 }, { 0xb4, 01454107632 },
            { 0xe1, 01475316420 }, { 0x39, 01445670123 },
            { 0x36, 01451734062 }, { 0x63, 01410325476 }, { 0x93, 01404261537 } };

    /* PATTERN 13 */
    private static int PAT13[][] = { { 0x96, 01576543210 },
            { 0x69, 01537152604 } };

    /* PATTERN 14 */
    private static int PAT14[][] = { { 0xd8, 01676543210 },
            { 0xac, 01637152604 }, { 0x4e, 01623016745 },
            { 0xc5, 01662407351 }, { 0x3a, 01615043726 },
            { 0x27, 01601452367 }, { 0x53, 01640516273 },
            { 0xb1, 01654107632 }, { 0x8b, 01613570246 },
            { 0xe2, 01675316420 }, { 0x1d, 01602134657 }, { 0x74, 01646025713 } };

    /*-------------------------------ba9 876 543 210 */
    private static final int XTRACED = 01122; /* 001 001 010 010 */
    private static final int YTRACED = 00461; /* 000 100 110 001 */
    private static final int ZTRACED = 00017; /* 000 000 001 111 */

    private int BOUND(float x, float t) {
        return (x >= t ? 1 : 0);
    }

    private int EXTRACT(int x, int n) {
        // debug
        // System.out.println("EXTRACT x "+x+" n "+n+" r "+((x>>(n*3))&07));
        return (x >> (n * 3)) & 07;
    }

    private int EXTRACT1(int x, int n) {
        // debug
        // System.out.println("EXTRACT1 x "+x+" n "+n+" r "+((x>>n)&01));
        return (x >> n) & 01;
    }

    /*------------------------------------------------------------------*/
    /* one pattern table with 128 entries is defined */
    /* IAC code change - static */
    private int pattab[] = new int[PATNO];

    /* pre-computed dim0,dim1,dim2,dim01,dim012 */
    /* IAC code change - static */
    private int dim0, dim1, dim2, dim01, dim012, dim_01;

    /*------------------------------------------------------------------*/

    /*
     * Main routine to generate isosurface geometric object by Marching method
     * Note: the input field is a uniform one
     */

    /* STATIC GLOBAL VARIABLES */
    /* IAC code change - static */private float xp[], yp[], zp[];
    /* IAC code change - static */private Point3f b0[], f0[], m0[], b[], f[],
            m[], tbf[];
    /* IAC code change - static */private Vector3f b0n[], f0n[], m0n[], bn[],
            fn[], mn[], tbfn[];
    /* IAC code change - static */private int pno[] = new int[15];
    /* IAC code change - static */private int xtraced, ytraced, ztraced,
            traced;

    /**
     * 値の補間をLogで行うかLinearで行うかの指定 0:log, 1:linear
     */
    private int logOrLinear;

    /*------------------------------------------------------------------*/
    /**
     * ノードの作成
     * 
     * 
     * @return なし
     * @param x
     *            x座標
     * @param y
     *            y座標
     * @param z
     *            z座標
     * @param nx
     *            x座標の数
     * @param ny
     *            y座標の数
     * @param nz
     *            z座標の数
     * @param data
     *            格子上の値
     * @param value
     *            等値面の値
     * @param tri_pos
     *            等値面のポリゴンの座標
     * @param ll
     *            0:log, 1:linear
     */
    public void create_isosurf(float x[], float y[], float z[], int nx, int ny,
            int nz, float data[], float value, Vector tri_pos, Vector normal,
            int ll)
    // , float data2[])
    {

        /*------------------------------------------------------------------*/
        /* look up pattern code prototype */
        int typeno[] = new int[15];
        int pcode, ptype;
        /* miscellaneous local variable declarations */
        int zoffset0, zoffset1;
        int yoffset0, yoffset1;
        /*------------------------------------------------------------------*/

        logOrLinear = ll;
        init_pattab();

        /* get dimension for x, y and z and pre-compute multiply */
        dim0 = nx;
        dim1 = ny;
        dim2 = nz;
        dim01 = dim0 * dim1;
        dim012 = dim01 * dim2;

        try {
            b0 = new Point3f[2 * dim01 - dim0 - dim1];
            for (int i = 0; i < b0.length; i++)
                b0[i] = new Point3f();
            f0 = new Point3f[2 * dim01 - dim0 - dim1];
            for (int i = 0; i < f0.length; i++)
                f0[i] = new Point3f();
            m0 = new Point3f[dim01];
            for (int i = 0; i < m0.length; i++)
                m0[i] = new Point3f();
            b0n = new Vector3f[2 * dim01 - dim0 - dim1];
            for (int i = 0; i < b0n.length; i++)
                b0n[i] = new Vector3f();
            f0n = new Vector3f[2 * dim01 - dim0 - dim1];
            for (int i = 0; i < f0n.length; i++)
                f0n[i] = new Vector3f();
            m0n = new Vector3f[dim01];
            for (int i = 0; i < m0n.length; i++)
                m0n[i] = new Vector3f();
        } catch (Exception e) {
            System.err.println("MarchingCube/create_isosurf: " + e);
            return;
        }

        /* xp -->start of all x coord under dimension */
        xp = x;
        /* yp -->start of all y coord under dimension */
        yp = y;
        /* zp -->start of all z coord under dimension */
        zp = z;

        for (int i = 0; i < 15; i++)
            pno[i] = 0;

        /* for all cubes match pattern code table to generate triangles */
        zoffset1 = 0;
        b = b0;
        m = m0;
        f = f0;
        bn = b0n;
        mn = m0n;
        fn = f0n;
        for (int k = 0, ztraced = 0; k < dim2 - 1; k++, ztraced = ZTRACED) {
            // @@zoffset1= dim01*5;
            // @@for (k=5,ztraced=0;k<6;k++,ztraced=ZTRACED) {
            // if (k==5) break; /*@@@@@ for debug */
            zoffset0 = zoffset1;
            zoffset1 = zoffset0 + dim01;
            /* deal with only one z plane cubes */
            read1plane(data, value, zoffset0, k, tri_pos, normal);
            // debaug用 ,data2);

            /* redirect front pointer to back pointer */
            tbf = f;
            f = b;
            b = tbf;
            tbfn = fn;
            fn = bn;
            bn = tbfn;
        }
        // for (int i=0;i<15;i++) System.out.println(i+" "+ pno[i]);
    }

    /* read 1st cubic plane with positions and normals */
    /**
     * xy平面基準のノードの作成
     * 
     * 
     * @return なし
     * @param data
     *            格子上の値
     * @param lev
     *            等値面の値
     * @param zoff
     * @param z
     * @param tri_pos
     *            等値面のポリゴンの座標
     */
    private void read1plane(float data[], float lev, int zoff, int z,
            Vector tri_pos, Vector normal)
    // , float data2[])
    {
        /*
         * ------------------------------------local variables
         * define---------------------
         */
        float vdata1 = 0, vdata2 = 0;
        int pcode;
        int ptype;
        int cedge;
        int uedge;
        int c1, c2;
        Point3f addressv;
        Point3f edge_array[];
        Vector3f addressn;
        Vector3f normal_array[];
        byte plyedges[][] = null;
        byte PlyNo;
        int edge_cnvt[] = new int[12];
        Point3f v0 = new Point3f();
        Point3f v1 = new Point3f();
        Vector3f n0 = new Vector3f();
        Vector3f n1 = new Vector3f();
        int yoffset1, yoffset0, zoffset1, zoffset0;

        /*
         * ------------------------------------z & y offset
         * initialization-----------------
         */
        zoffset0 = zoff;
        zoffset1 = zoffset0 + dim01;
        yoffset1 = 0;

        edge_array = new Point3f[12];
        normal_array = new Vector3f[12];
        for (int i = 0; i < edge_array.length; i++) {
            edge_array[i] = new Point3f();
            normal_array[i] = new Vector3f();
        }

        /*-----------------------------------begin to find isosurface----------------------*/
        /* loop for y direction with yExt=dim1-1 */
        for (int i = 0, ytraced = ztraced; i < dim1 - 1; i++, ytraced = ztraced, ytraced |= YTRACED) {
            // if (i==5) break; /*@@@@@ for debug */
            // @@yoffset1=dim0*22;
            // @@for (int
            // i=22,ytraced=ztraced;i<23;i++,ytraced=ztraced,ytraced|=YTRACED) {
            yoffset0 = yoffset1;
            yoffset1 = yoffset0 + dim0;
            /*
             * ----------------------loop for x direction with
             * xExt=dim0-1-------------
             */
            for (int j = 0, xtraced = ytraced; j < dim0 - 1; j++, xtraced = ytraced, xtraced |= XTRACED) {

                /*
                 * debug用 find_offset(data,0,zoffset0,yoffset0,j,i,z,v0,n0);
                 * System.out.println(j+" "+i+" "+z+" "+
                 * data[j+yoffset0+zoffset0]+" "+
                 * data2[j+yoffset0+zoffset0]+" "+ n0.x+" "+n0.y+" "+n0.z);
                 */

                // if (j==5) break; /*@@@@@ for debug */
                // @@for (int
                // j=7,xtraced=ytraced;j<8;j++,xtraced=ytraced,xtraced|=XTRACED)
                // {
                /* get pattern code from the current cubic */

                // int code=0,tcode,sum1=0;
                // int z1=z0+dim01;
                // int y1=y0+dim0;
                // /* 7 */
                // tcode=BOUND(d[z1+y1+x+1],lev);
                // code|=tcode; code<<=1; sum1+=tcode;
                // /* 6 */
                // tcode=BOUND(d[z1+y1+x],lev);
                // code|=tcode; code<<=1; sum1+=tcode;
                // /* 5 */
                // tcode=BOUND(d[z1+y0+x+1],lev);
                // code|=tcode; code<<=1; sum1+=tcode;
                // /* 4 */
                // tcode=BOUND(d[z1+y0+x],lev);
                // code|=tcode; code<<=1; sum1+=tcode;
                // /* 3 */
                // tcode=BOUND(d[z0+y1+x+1],lev);
                // code|=tcode; code<<=1; sum1+=tcode;
                // /* 2 */
                // tcode=BOUND(d[z0+y1+x],lev);
                // code|=tcode; code<<=1; sum1+=tcode;
                // /* 1 */
                // tcode=BOUND(d[z0+y0+x+1],lev);
                // code|=tcode; code<<=1; sum1+=tcode;
                // /* 0 */
                // tcode=BOUND(d[z0+y0+x],lev);
                // code|=tcode; sum1+=tcode;

                pcode = getpcode(data, zoffset0, yoffset0, j, lev);
                ptype = pcode;
                ptype = (ptype >> 24) & 077;
                pno[ptype]++;
                if (ptype > 14) {
                    System.out.println("error ");
                    continue;
                }

                /*
                 * debug if ((j==17||j==18) && i==10 && z==12) { int
                 * xn=j,yn=yoffset0,zn=zoffset0;
                 * System.out.println("x,y,z       "+data[j+yoffset0+zoffset0]);
                 * System
                 * .out.println("x+1,y,z     "+data[j+1+yoffset0+zoffset0]);
                 * System
                 * .out.println("x,y+1,z     "+data[j+yoffset0+dim0+zoffset0]);
                 * System
                 * .out.println("x+1,y+1,z   "+data[j+1+yoffset0+dim0+zoffset0
                 * ]);
                 * System.out.println("x,y,z+1     "+data[j+yoffset0+zoffset0
                 * +dim01]);
                 * System.out.println("x+1,y,z+1   "+data[j+1+yoffset0+
                 * zoffset0+dim01]);
                 * System.out.println("x,y+1,z+1   "+data[j+yoffset0
                 * +dim0+zoffset0+dim01]);
                 * System.out.println("x+1,y+1,z+1 "+data
                 * [j+1+yoffset0+dim0+zoffset0+dim01]);
                 * System.out.println("ptype  "+ptype); }
                 */
                /* get cross edge code for ptype */
                cedge = CROSSEDGE[ptype];
                if (cedge == 0)
                    continue; // 接点無しなので次のセルへ

                // debug
                // System.out.println("cedge "+cedge+" ptype "+ptype+
                // " pcode "+pcode);

                /*
                 * -------------------calculate interseted position
                 * --------------------
                 */
                for (int k = 0; k < 12; cedge >>= 1, k++) {
                    /*
                     * only consider the edges that are intersected by
                     * isosurface
                     */
                    if ((cedge & 1) == 1) {

                        /* convert real corner index for these two end points */
                        c1 = EXTRACT(pcode, UEDGETABLE[k][0]);
                        c2 = EXTRACT(pcode, UEDGETABLE[k][1]);

                        /* decide real edge that is intersected */
                        uedge = EDGETABLE[c1][c2];

                        /* ready for retrieval later */
                        edge_cnvt[k] = uedge;
                        // debug
                        // System.out.println("k "+k+" uedge "+uedge+" c1 "+c1+" c2 "+c2);

                        if (EXTRACT1(xtraced, uedge) != 0) {
                            /*
                             * this is a shared edge and only get data from
                             * cache
                             */
                            // debug
                            // System.out.println("get uedge "+uedge);
                            switch (uedge) {
                            case 0:
                                copyPoint3f(edge_array[uedge], f[yoffset0
                                        + yoffset0 - i + j]);
                                copyVector3f(normal_array[uedge], fn[yoffset0
                                        + yoffset0 - i + j]);
                                break;
                            case 1:
                                copyPoint3f(edge_array[uedge], f[dim0 - 1
                                        + yoffset0 + yoffset0 - i + j]);
                                copyVector3f(normal_array[uedge], fn[dim0 - 1
                                        + yoffset0 + yoffset0 - i + j]);
                                break;
                            case 2:
                                copyPoint3f(edge_array[uedge], f[dim0 - 1
                                        + yoffset0 + yoffset0 - i + j + 1]);
                                copyVector3f(normal_array[uedge], fn[dim0 - 1
                                        + yoffset0 + yoffset0 - i + j + 1]);
                                break;
                            case 3:
                                copyPoint3f(edge_array[uedge], f[yoffset0
                                        + dim0 + yoffset0 + dim0 - i - 1 + j]);
                                copyVector3f(normal_array[uedge], fn[yoffset0
                                        + dim0 + yoffset0 + dim0 - i - 1 + j]);
                                break;
                            case 4:
                                copyPoint3f(edge_array[uedge], m[yoffset0 + j]);
                                copyVector3f(normal_array[uedge], m0n[yoffset0
                                        + j]);
                                break;
                            case 5:
                                copyPoint3f(edge_array[uedge], m[yoffset0 + j
                                        + 1]);
                                copyVector3f(normal_array[uedge], m0n[yoffset0
                                        + j + 1]);
                                break;
                            case 6:
                                copyPoint3f(edge_array[uedge], m[yoffset1 + j]);
                                copyVector3f(normal_array[uedge], mn[yoffset1
                                        + j]);
                                break;
                            case 8:
                                copyPoint3f(edge_array[uedge], b[yoffset0
                                        + yoffset0 - i + j]);
                                copyVector3f(normal_array[uedge], bn[yoffset0
                                        + yoffset0 - i + j]);
                                break;
                            case 9:
                                copyPoint3f(edge_array[uedge], b[dim0 - 1
                                        + yoffset0 + yoffset0 - i + j]);
                                copyVector3f(normal_array[uedge], bn[dim0 - 1
                                        + yoffset0 + yoffset0 - i + j]);
                                break;
                            case 7:
                            case 10:
                            case 11:
                            default:
                                break;
                            }

                        }
                        /* v,n must be generated and copy them onto bmf */
                        else {
                            vdata1 = find_offset(data, c1, zoffset0, yoffset0,
                                    j, i, z, v0, n0);
                            vdata2 = find_offset(data, c2, zoffset0, yoffset0,
                                    j, i, z, v1, n1);
                            /* here v0,n0,vdata1 & v1,n1,vdata2 are available */
                            /* interpolate these two end points into the edge */
                            interpolate(v0, v1, vdata1, vdata2, lev,
                                    edge_array[uedge]);
                            interpolate(n0, n1, vdata1, vdata2, lev,
                                    normal_array[uedge]);

                            // debug
                            // System.out.println("u "+uedge+
                            // " edge "+edge_array[uedge]);

                            /* record it into cache */
                            switch (uedge) {
                            case 0:
                                copyPoint3f(f[yoffset0 + yoffset0 - i + j],
                                        edge_array[uedge]);
                                copyVector3f(fn[yoffset0 + yoffset0 - i + j],
                                        normal_array[uedge]);
                                break;
                            case 1:
                                copyPoint3f(f[dim0 - 1 + yoffset0 + yoffset0
                                        - i + j], edge_array[uedge]);
                                copyVector3f(fn[dim0 - 1 + yoffset0 + yoffset0
                                        - i + j], normal_array[uedge]);
                                break;
                            case 2:
                                copyPoint3f(f[dim0 - 1 + yoffset0 + yoffset0
                                        - i + j + 1], edge_array[uedge]);
                                copyVector3f(fn[dim0 - 1 + yoffset0 + yoffset0
                                        - i + j + 1], normal_array[uedge]);
                                break;
                            case 3:
                                copyPoint3f(f[yoffset0 + dim0 + yoffset0 + dim0
                                        - i - 1 + j], edge_array[uedge]);
                                copyVector3f(fn[yoffset0 + dim0 + yoffset0
                                        + dim0 - i - 1 + j],
                                        normal_array[uedge]);
                                break;
                            case 4:
                                copyPoint3f(m[yoffset0 + j], edge_array[uedge]);
                                copyVector3f(mn[yoffset0 + j],
                                        normal_array[uedge]);
                                break;
                            case 5:
                                copyPoint3f(m[yoffset0 + j + 1],
                                        edge_array[uedge]);
                                copyVector3f(mn[yoffset0 + j + 1],
                                        normal_array[uedge]);
                                break;
                            case 6:
                                copyPoint3f(m[yoffset1 + j], edge_array[uedge]);
                                copyVector3f(mn[yoffset1 + j],
                                        normal_array[uedge]);
                                break;
                            case 8:
                                copyPoint3f(b[yoffset0 + yoffset0 - i + j],
                                        edge_array[uedge]);
                                copyVector3f(bn[yoffset0 + yoffset0 - i + j],
                                        normal_array[uedge]);
                                break;
                            case 9:
                                copyPoint3f(b[dim0 - 1 + yoffset0 + yoffset0
                                        - i + j], edge_array[uedge]);
                                copyVector3f(bn[dim0 - 1 + yoffset0 + yoffset0
                                        - i + j], normal_array[uedge]);
                                break;
                            case 7:
                                copyPoint3f(m[yoffset1 + j + 1],
                                        edge_array[uedge]);
                                copyVector3f(mn[yoffset1 + j + 1],
                                        normal_array[uedge]);
                                break;
                            case 10:
                                copyPoint3f(b[dim0 - 1 + yoffset0 + yoffset0
                                        - i + j + 1], edge_array[uedge]);
                                copyVector3f(bn[dim0 - 1 + yoffset0 + yoffset0
                                        - i + j + 1], normal_array[uedge]);
                                break;
                            case 11:
                                copyPoint3f(b[yoffset0 + dim0 + yoffset0 + dim0
                                        - i - 1 + j], edge_array[uedge]);
                                copyVector3f(bn[yoffset0 + dim0 + yoffset0
                                        + dim0 - i - 1 + j],
                                        normal_array[uedge]);
                                break;
                            default:
                                break;
                            }

                            /*
                             * assign address pointer to be ready for fill
                             * coords and norms (8) 4------------5 (4) /| (5) /|
                             * 0------------1| || (0) || edge index assignment
                             * (0..11) ||(9) ||(10) (1) || (2) || corner index
                             * assignment (0..7) |6___(11)____|7 2/(6)_
                             * ______3/(7) (3)
                             */
                        }

                    }
                    /* unsetting edge doesn't care !! */
                }
                /* complete all these 12 edge checkings */
                /*
                 * PlyNo : decides how many triangles to generate plyedges :
                 * defines all triangle edge index both variables must be
                 * pointered to appropriate data structure for different type
                 */
                /*-----------------------------------------------------------*/
                switch (ptype) {
                case 1:
                    PlyNo = (byte) PlyNo1;
                    plyedges = PLYEDGES1;
                    break;
                case 2:
                    PlyNo = (byte) PlyNo2;
                    plyedges = PLYEDGES2;
                    break;
                case 3:
                    PlyNo = (byte) PlyNo3;
                    plyedges = PLYEDGES3;
                    break;
                case 4:
                    PlyNo = (byte) PlyNo4;
                    plyedges = PLYEDGES4;
                    break;
                case 5:
                    PlyNo = (byte) PlyNo5;
                    plyedges = PLYEDGES5;
                    break;
                case 6:
                    PlyNo = (byte) PlyNo6;
                    plyedges = PLYEDGES6;
                    break;
                case 7:
                    PlyNo = (byte) PlyNo7;
                    plyedges = PLYEDGES7;
                    break;
                case 8:
                    PlyNo = (byte) PlyNo8;
                    plyedges = PLYEDGES8;
                    break;
                case 9:
                    PlyNo = (byte) PlyNo9;
                    plyedges = PLYEDGES9;
                    break;
                case 10:
                    PlyNo = (byte) PlyNo10;
                    plyedges = PLYEDGES10;
                    break;
                case 11:
                    PlyNo = (byte) PlyNo11;
                    plyedges = PLYEDGES11;
                    break;
                case 12:
                    PlyNo = (byte) PlyNo12;
                    plyedges = PLYEDGES12;
                    break;
                case 13:
                    PlyNo = (byte) PlyNo13;
                    plyedges = PLYEDGES13;
                    break;
                case 14:
                    PlyNo = (byte) PlyNo14;
                    plyedges = PLYEDGES14;
                    break;
                default:
                    PlyNo = 0;
                    break;
                }
                /*------------------------------------------------------------*/
                /*
                 * ---------------here actually generate geometry object for
                 * triangles --------------
                 */

                float ignore = VolumetricData.meaninglessValue;
                int z1 = zoffset0 + dim01;
                int y1 = yoffset0 + dim0;
                int x = j;
                int z0 = zoffset0;
                int y0 = yoffset0;
                if (data[z1 + y1 + x + 1] == ignore
                        || data[z1 + y1 + x] == ignore
                        || data[z1 + y0 + x + 1] == ignore
                        || data[z1 + y0 + x] == ignore
                        || data[z0 + y1 + x + 1] == ignore
                        || data[z0 + y1 + x] == ignore
                        || data[z0 + y0 + x + 1] == ignore
                        || data[z0 + y0 + x] == ignore)
                    continue;

                for (int k = 0; k < PlyNo; k++) {
                    /* convert to correct edge lists */
                    addPolygon(plyedges[k][0], plyedges[k][1], plyedges[k][2],
                            lev, edge_cnvt, edge_array, normal_array, tri_pos,
                            normal);
                }

                // marching cubeの穴が空くところをふさぐ
                // ptype == 3のとき頂点0,1,2,3で接するセルのパターンが3または10の時と
                // ptype == 12のとき頂点0,1,2,3で接するセルのパターンが6の時とき
                // 0,1,2,3をポリゴンに加える
                boolean hanten1 = hanten;
                if (hanten && ptype == 3 || ptype == 6 || ptype == 7) {
                    // ptype == 10|| // @@@@
                    // ptype== 12 || ptype == 13) { // @@@
                    int xn = j, yn = yoffset0, zn = zoffset0;
                    int voffset = xn + yn + zn;
                    int e0;
                    int e1;
                    int e2;

                    if (ptype == 3) {
                        e0 = edge_cnvt[0];
                        e1 = edge_cnvt[1];
                        e2 = edge_cnvt[2];
                    } else {
                        e0 = edge_cnvt[2];
                        e1 = edge_cnvt[5];
                        e2 = edge_cnvt[7];
                    }

                    // System.out.println("orignal xn "+xn+" yn "+yn+" zn "+zn);

                    if (sameValue(edge_array[e0].x, edge_array[e1].x)
                            && sameValue(edge_array[e1].x, edge_array[e2].x)
                            && sameValue(edge_array[e0].x, edge_array[e2].x)) {
                        if (sameValue(edge_array[e1].x, xp[voffset]))
                            xn--;
                        else
                            xn++;
                    } else if (sameValue(edge_array[e0].y, edge_array[e1].y)
                            && sameValue(edge_array[e1].y, edge_array[e2].y)
                            && sameValue(edge_array[e0].y, edge_array[e2].y)) {
                        if (sameValue(edge_array[e1].y, yp[voffset]))
                            yn = yoffset0 - dim0;
                        else
                            yn = yoffset0 + dim0;
                    } else {
                        if (sameValue(edge_array[e1].z, zp[voffset]))
                            zn = zoffset0 - dim01;
                        else
                            zn = zoffset0 + dim01;
                    }
                    // System.out.println("\n"+j+" "+i+" "+z+
                    // " x "+xp[voffset]+" y "+yp[voffset]+" z "+zp[voffset]+" ptype "+ptype);
                    // System.out.println("xn "+xn+" yn "+yn+" zn "+zn);
                    // System.out.println("e0 "+edge_array[e0]);
                    // System.out.println("e1 "+edge_array[e1]);
                    // System.out.println("e2 "+edge_array[e2]);

                    int next_ptype = getpcode(data, zn, yn, xn, lev);
                    next_ptype = (next_ptype >> 24) & 077;

                    // System.out.println("next_ptype "+next_ptype);
                    if (!hanten
                            && (next_ptype == 3 || next_ptype == 6
                                    || next_ptype == 7 || next_ptype == 10
                                    || next_ptype == 12 || next_ptype == 13)) {

                        if (ptype == 3) {
                            addPolygon(0, 1, 3, lev, edge_cnvt, edge_array,
                                    normal_array, tri_pos, normal);
                            addPolygon(0, 3, 2, lev, edge_cnvt, edge_array,
                                    normal_array, tri_pos, normal);
                        } else {
                            addPolygon(2, 5, 7, lev, edge_cnvt, edge_array,
                                    normal_array, tri_pos, normal);
                            addPolygon(7, 5, 10, lev, edge_cnvt, edge_array,
                                    normal_array, tri_pos, normal);
                        }
                        // System.out.println("patch "+j+" "+i+" "+z);
                    }
                }
            }
            /* end of j */
        }
        /* end of i */
    }

    /**
     * ポリゴンの追加
     * 
     * @param int i1 頂点番号
     * @param int i2 頂点番号
     * @param int i3 頂点番号
     * @param float v 値
     * @param int [] edge_cnvt 頂点変換テーブル
     * @param Point3f
     *            [] edge_array 頂点の座標テーブル
     * @param Vector3f
     *            [] normal_array 法線のテーブル
     * @param vector
     *            tri_pos 等値面を構成するポリゴンの頂点座標
     * @param vector
     *            normal 等値面の頂点の法線
     */
    private void addPolygon(int i1, int i2, int i3, float v, int[] edge_cnvt,
            Point3f[] edge_array, Vector3f[] normal_array, Vector tri_pos,
            Vector normal) {
        int e1, e2, e3;
        e1 = edge_cnvt[i1];
        e2 = edge_cnvt[i2];
        e3 = edge_cnvt[i3];

        //

        // メッシュの切れ目で変なことになるのを防ぐためこの処理を追加
        // こんなことはなさそうなのでコメントにする 2003/2/18
        /*
         * float d1 = edge_array[e1].distance(edge_array[e2]); float d2 =
         * edge_array[e2].distance(edge_array[e3]); float d3 =
         * edge_array[e3].distance(edge_array[e1]); float d = (xp[1] - xp[0])
         * *2; if (d1 > d ||d2 > d ||d3 > d) {
         * System.out.println("i1 i2 i3 "+i1+" "+i2+" "+i3+" ");
         * System.out.println("e1 e2 e3 "+e1+" "+e2+" "+e3+" ");
         * System.out.println("d1 d2 d3 "+d1+" "+d2+" "+d3+" "); for (int i=0;
         * i<edge_array.length; i++) System.out.println(i+" "+edge_array[i]);
         * System.out.println(""); //return; }
         */

        Vector3f n1 = new Vector3f(normal_array[e1]);
        Vector3f n2 = new Vector3f(normal_array[e2]);
        Vector3f n3 = new Vector3f(normal_array[e3]);

        if (v > 0f) {
            n1.negate();
            n2.negate();
            n3.negate();
        }

        // ここでポリゴンの向きをそろえる
        // これをしないとうまくシェーディングされない
        // 平面の法線と頂点の法線ベクトルの向きを比べてポリゴンの順番を調整する
        Vector3f vv1 = new Vector3f(edge_array[e1]);
        Vector3f vv2 = new Vector3f(edge_array[e2]);
        Vector3f vv3 = new Vector3f(edge_array[e3]);

        vv3.sub(vv2);
        vv2.sub(vv1);
        vv1.cross(vv3, vv2); // 平面の法線

        /*
         * とりあえず元に戻しておく Vector3f base; if (!isEdge(edge_array[e1])) base = n1;
         * else if (!isEdge(edge_array[e2])) base = n2; else base = n3;
         */

        // 面に平行なものは境界付近のせいか異常なものとして
        // 無視して、その他の法線で判断する。
        int count = 0;
        if (vv1.angle(n1) < (float) Math.PI * 0.4f)
            count++;
        if (vv1.angle(n2) < (float) Math.PI * 0.4f)
            count++;
        if (vv1.angle(n3) < (float) Math.PI * 0.4f)
            count++;

        if (vv1.angle(n1) > (float) Math.PI * 0.6f)
            count--;
        if (vv1.angle(n2) > (float) Math.PI * 0.6f)
            count--;
        if (vv1.angle(n3) > (float) Math.PI * 0.6f)
            count--;

        // if (vv1.angle(n1) > (float) Math.PI/2.0f) {
        if (count >= 1) {
            tri_pos.addElement(new Point3f(edge_array[e1]));
            tri_pos.addElement(new Point3f(edge_array[e3]));
            tri_pos.addElement(new Point3f(edge_array[e2]));
            normal.addElement(n1);
            normal.addElement(n3);
            normal.addElement(n2);
        } else {
            tri_pos.addElement(new Point3f(edge_array[e1]));
            tri_pos.addElement(new Point3f(edge_array[e2]));
            tri_pos.addElement(new Point3f(edge_array[e3]));
            normal.addElement(n1);
            normal.addElement(n2);
            normal.addElement(n3);
        }
    }

    private boolean sameValue(float v1, float v2) {
        return (Math.abs(v1 - v2) < 1e-5f);
    }

    /**
     * 指定された点が格子の端の点かどうかを返す
     * 
     */
    private boolean isEdge(Point3f p) {
        /*
         * if (p.x == xp[0] || p.x == xp[xp.length-1] || p.y == yp[0] || p.y ==
         * yp[yp.length-1] || p.z == zp[0] || p.z == zp[zp.length-1]) {
         */
        if (sameValue(p.x, xp[0]) || sameValue(p.x, xp[xp.length - 1])
                || sameValue(p.y, yp[0]) || sameValue(p.y, yp[yp.length - 1])
                || sameValue(p.z, zp[0]) || sameValue(p.z, zp[zp.length - 1])) {
            System.out.println(p);
            return true;
        } else {
            return false;
        }
    }

    /**
     * 面の法線を計算する。
     * 
     * @param Point3f
     *            [] vertices 頂点の座標
     * @param int n vertices 頂点の開始番号
     */
    private Vector3f calcFaceNormal(Point3f[] vertices, int n) {
        float ax = vertices[n + 2].x - vertices[n + 1].x;
        float ay = vertices[n + 2].y - vertices[n + 1].y;
        float az = vertices[n + 2].z - vertices[n + 1].z;
        float bx = vertices[n + 0].x - vertices[n + 1].x;
        float by = vertices[n + 0].y - vertices[n + 1].y;
        float bz = vertices[n + 0].z - vertices[n + 1].z;

        float nx = ay * bz - az * by;
        float ny = az * bx - ax * bz;
        float nz = ax * by - ay * bx;

        Vector3f normal = new Vector3f((float) nx, (float) ny, (float) nz);
        normal.normalize();

        // DEBUG
        // System.out.println("nx=" + normal.x + ", ny=" + normal.y + ", nz=" +
        // normal.z);

        return normal;
    }

    /* get vert,norm and data for interesting cubic index */
    private float find_offset(float d[], int c, int z0, int y0, int x, int y,
            int z, Point3f v, Vector3f n) {
        /*------------------------------------------------------------*/
        int offset = 0;
        int z1 = z0 + dim01;
        int y1 = y0 + dim0;
        float Gx, Gy, Gz, Nx, Ny, Nz, G, H;
        int X = x;
        int Y = y;
        int Z = z;
        /*------------------------------------------------------------*/

        /* adjust X,Y and Z and the offset to get data from */
        switch (c) {
        case 0:
            offset = z0 + y0 + x;
            break;
        case 1:
            offset = z0 + y0 + x + 1;
            X++;
            break;
        case 2:
            offset = z0 + y1 + x;
            Y++;
            break;
        case 3:
            offset = z0 + y1 + x + 1;
            X++;
            Y++;
            break;
        case 4:
            offset = z1 + y0 + x;
            Z++;
            break;
        case 5:
            offset = z1 + y0 + x + 1;
            Z++;
            X++;
            break;
        case 6:
            offset = z1 + y1 + x;
            Z++;
            Y++;
            break;
        case 7:
            offset = z1 + y1 + x + 1;
            Z++;
            Y++;
            X++;
            break;
        default:
            break;
        }

        /* get real coordinates for vertex v */
        // System.out.println("offset: "+offset+" xp.length: "+xp.length);
        v.x = xp[offset];
        v.y = yp[offset];
        v.z = zp[offset];

        // 境界を0とすると法線の向きがうまく計算されない場合がある。
        // ナノの場合穴があいてしまう。
        // try{
        if (X == 0)
            n.x = (d[offset + 1] - d[offset]) * 2;
        else if (X == (dim0 - 1))
            n.x = (d[offset] - d[offset - 1]) * 2;
        else
            n.x = d[offset + 1] - d[offset - 1];

        if (Y == 0)
            n.y = (d[offset + dim0] - d[offset]) * 2;
        else if (Y == (dim1 - 1))
            n.y = (d[offset] - d[offset - dim0]) * 2;
        else
            n.y = d[offset + dim0] - d[offset - dim0];

        if (Z == 0)
            n.z = (d[offset + dim01] - d[offset]) * 2;
        else if (Z == (dim2 - 1))
            n.z = (d[offset] - d[offset - dim01]) * 2;
        else
            n.z = d[offset + dim01] - d[offset - dim01];
        // } catch(Exception exc){}
        /*
         * 2003/03/30 にコメント if (X==0 || X==(dim0-1)) n.x=0.0f; else
         * n.x=d[offset+1]-d[offset-1]; if (Y==0 || Y==(dim1-1)) n.y=0.0f; else
         * n.y=d[offset+dim0]-d[offset-dim0]; if (Z==0 || Z==(dim2-1)) n.z=0.0f;
         * else n.z=d[offset+dim01]-d[offset-dim01];
         */

        /*
         * if (X==0) n.x= -1; else if (X==(dim0-1)) n.x=1; else
         * n.x=d[offset+1]-d[offset-1];
         * 
         * if (Y==0) n.y= -1; else if (Y==(dim1-1)) n.y= 1; else
         * n.y=d[offset+dim0]-d[offset-dim0];
         * 
         * if (Z==0) n.z= -1; else if (Z==(dim2-1)) n.z= 1; else
         * n.z=d[offset+dim01]-d[offset-dim01];
         * 
         * if (n.x == 1 || n.x == -1) { n.y = 0; n.z = 0; return(d[offset]); }
         * if (n.y == 1 || n.y == -1) { n.x = 0; n.z = 0; return(d[offset]); }
         * if (n.z == 1 || n.z == -1) { n.x = 0; n.y = 0; return(d[offset]); }
         */

        Gx = n.x;
        Gy = n.y;
        Gz = n.z;

        /*
         * if ((Gx= *n)<0.0) Gx=0.0-Gx; if ((Gy= *(n+1))<0.0) Gy=0.0-Gy; if
         * ((Gz= *(n+2))<0.0) Gz=0.0-Gz;
         */

        G = Gx * Gx + Gy * Gy + Gz * Gz;
        if (G > 0.0) {
            H = (float) Math.sqrt((double) G);
            n.x = Gx / H;
            n.y = Gy / H;
            n.z = Gz / H;
        }

        /* get vertice data */
        return (d[offset]);

    }

    /* interpolate data for coord */
    private void interpolate(Point3f v1, Point3f v2, float d1, float d2,
            float lev, Point3f vert) {
        float p1[] = new float[3];
        float p2[] = new float[3];
        float c[] = new float[3];

        // debug
        // System.out.println("d1 "+d1+" d2 "+d2);
        // System.out.println("v1 "+v1);
        // System.out.println("v2 "+v2);

        p1[0] = v1.x;
        p1[1] = v1.y;
        p1[2] = v1.z;
        p2[0] = v2.x;
        p2[1] = v2.y;
        p2[2] = v2.z;

        if (logOrLinear == LOG) {
            CrossPoint.get(d1, d2, lev, p1, p2, c);
        } else {
            CrossPointLinear.get(d1, d2, lev, p1, p2, c);
        }

        vert.x = c[0];
        vert.y = c[1];
        vert.z = c[2];
    }

    /* interpolate data for coord */
    private void interpolate(Vector3f v1, Vector3f v2, float d1, float d2,
            float lev, Vector3f vert) {
        float p1[] = new float[3];
        float p2[] = new float[3];
        float c[] = new float[3];

        // debug
        // System.out.println("d1 "+d1+" d2 "+d2);
        // System.out.println("v1 "+v1);
        // System.out.println("v2 "+v2);

        p1[0] = v1.x;
        p1[1] = v1.y;
        p1[2] = v1.z;
        p2[0] = v2.x;
        p2[1] = v2.y;
        p2[2] = v2.z;

        if (logOrLinear == LOG) {
            CrossPoint.get(d1, d2, lev, p1, p2, c);
        } else {
            CrossPointLinear.get(d1, d2, lev, p1, p2, c);
        }

        vert.x = c[0];
        vert.y = c[1];
        vert.z = c[2];
    }

    /*
     * get pattern code for some cube the eight pos for this cube is arranged as
     * the following :
     * 
     * 
     * z1+y0+x z1+y0+x+1 (4)X=============X(5) /| /| (0)X=============X(1)
     * z0+y0+x z0+y0+x+1 ======> bit 7 6 5 4 3 2 1 0 | | | | -------------------
     * | |(6) | |(7) 0 0 1 1 0 1 0 0 |/z1+y1+x=====|/z1+y1+x+1 LOOK UP PATTAB
     * (2)X=============X(3) || z0+y1+x z0+y1+x+1 (int)(6 bits)(8X3 bits) | |
     * pattern type | result after rotations
     */
    private int getpcode(float d[], int z0, int y0, int x, float lev) {
        try {
            int code = 0, tcode, sum1 = 0;
            int z1 = z0 + dim01;
            int y1 = y0 + dim0;
            /* 7 */
            tcode = BOUND(d[z1 + y1 + x + 1], lev);
            code |= tcode;
            code <<= 1;
            sum1 += tcode;
            /* 6 */
            tcode = BOUND(d[z1 + y1 + x], lev);
            code |= tcode;
            code <<= 1;
            sum1 += tcode;
            /* 5 */
            tcode = BOUND(d[z1 + y0 + x + 1], lev);
            code |= tcode;
            code <<= 1;
            sum1 += tcode;
            /* 4 */
            tcode = BOUND(d[z1 + y0 + x], lev);
            code |= tcode;
            code <<= 1;
            sum1 += tcode;
            /* 3 */
            tcode = BOUND(d[z0 + y1 + x + 1], lev);
            code |= tcode;
            code <<= 1;
            sum1 += tcode;
            /* 2 */
            tcode = BOUND(d[z0 + y1 + x], lev);
            code |= tcode;
            code <<= 1;
            sum1 += tcode;
            /* 1 */
            tcode = BOUND(d[z0 + y0 + x + 1], lev);
            code |= tcode;
            code <<= 1;
            sum1 += tcode;
            /* 0 */
            tcode = BOUND(d[z0 + y0 + x], lev);
            code |= tcode;
            sum1 += tcode;
            /* here sum1 accumulates # of 1 in code */

            hanten = false;
            /* if # of 1 greater than four then inverse it */
            if (sum1 > 4) {
                code = 255 - code;
                // debug
                // System.out.println("sum1 "+sum1+" "+code+" "+pattab[code]);
                hanten = true;
            }

            /* return code in the pattern table */
            return (pattab[code]);

        } catch (Exception e) {
            // 　s.denで0.001で描こうとするとエラーになる。なぜだろう？　2003/03/30
            // e.printStackTrace();
            return (0);
        }

    }

    // 反転したかどうかのフラグ
    private boolean hanten;

    private void init_pattab() {
        int i;

        /* set initial pattern value */
        for (i = 0; i <= PATNO - 1; i++) {
            pattab[i] = 01700000000; /* only for checking empty pattern use */
        }
        /* fill special pattern values from predefined ones */
        for (i = 0; i < PAT0.length; i++)
            pattab[PAT0[i][0]] = PAT0[i][1];

        for (i = 0; i < PAT1.length; i++)
            pattab[PAT1[i][0]] = PAT1[i][1];

        for (i = 0; i < PAT2.length; i++)
            pattab[PAT2[i][0]] = PAT2[i][1];

        for (i = 0; i < PAT3.length; i++)
            pattab[PAT3[i][0]] = PAT3[i][1];

        for (i = 0; i < PAT4.length; i++)
            pattab[PAT4[i][0]] = PAT4[i][1];

        for (i = 0; i < PAT5.length; i++)
            pattab[PAT5[i][0]] = PAT5[i][1];

        for (i = 0; i < PAT6.length; i++)
            pattab[PAT6[i][0]] = PAT6[i][1];

        for (i = 0; i < PAT7.length; i++)
            pattab[PAT7[i][0]] = PAT7[i][1];

        for (i = 0; i < PAT8.length; i++)
            pattab[PAT8[i][0]] = PAT8[i][1];

        for (i = 0; i < PAT9.length; i++)
            pattab[PAT9[i][0]] = PAT9[i][1];

        for (i = 0; i < PAT10.length; i++)
            pattab[PAT10[i][0]] = PAT10[i][1];

        for (i = 0; i < PAT11.length; i++)
            pattab[PAT11[i][0]] = PAT11[i][1];

        for (i = 0; i < PAT12.length; i++)
            pattab[PAT12[i][0]] = PAT12[i][1];

        for (i = 0; i < PAT13.length; i++)
            pattab[PAT13[i][0]] = PAT13[i][1];

        for (i = 0; i < PAT14.length; i++)
            pattab[PAT14[i][0]] = PAT14[i][1];
    }

    /**
     * Point3f の値のコピー
     */
    private void copyPoint3f(Point3f to, Point3f from) {
        to.x = from.x;
        to.y = from.y;
        to.z = from.z;
    }

    /**
     * Vector3f の値のコピー
     */
    private void copyVector3f(Vector3f to, Vector3f from) {
        to.x = from.x;
        to.y = from.y;
        to.z = from.z;
    }
}
