namespace H {
    public class Edge {
        public Vector3[] normals;
        public Vector3[] positions;
    }
    public class GridCell {
        public Edge[] edges;
        public float[] densities;
    }
    public class HermiteData {
        public GridCell[,,] cells;
    }

    private CMS.H.HermiteData CreateHermiteData(int resolution, Sample sample) {
        CMS.H.HermiteData h = new CMS.H.HermiteData();
        h.cells = new H.GridCell[resolution,resolution,resolution];

        for(int x = 0; x < resolution; x++) {
            for(int y = 0; y < resolution; y++) {
                for(int z = 0; z < resolution; z++) {
                    CMS.H.GridCell cell = new CMS.H.GridCell();
                    h.cells[x,y,z] = cell;
                    Vector3 min = new Vector3(x, y, z);
                    int corners = 0;

                    cell.densities = new float[8];
                    for(int i = 0; i < 8; i++) {
                        cell.densities[i] = sample(min + CMS.S.L.OFFSETS[i]);
                        if(cell.densities[i] > 0) corners |= (int)Mathf.Pow(2, i);
                    }


                    if(corners == 0 || corners == 255) {
                        continue;
                    }
                    cell.edges = new CMS.H.Edge[12];

                    for(int i = 0; i < 12; i++) {
                        cell.edges[i] = new CMS.H.Edge();
                        int c1 = S.L.edgevmap[i][0];
                        int c2 = S.L.edgevmap[i][1];

                        int m1 = cell.densities[c1] > 0 ? 1 : 0; //(corners >> c1) & 1;
                        int m2 = cell.densities[c2] > 0 ? 1 : 0; //(corners >> c2) & 1;

                        if((m1 == 0 && m2 == 0) || (m1 == 1 && m2 == 1)) {
                            continue;
                        }


                        Vector3 p1 = min + S.L.OFFSETS[c1];
                        Vector3 p2 = min + S.L.OFFSETS[c2];
                        Vector3 p = ApproximateZeroCrossingPosition(p1, p2, sample);

                        cell.edges[i].positions = new Vector3[1];
                        cell.edges[i].normals = new Vector3[1];

                        cell.edges[i].positions[0] = p;
                        cell.edges[i].normals[0] = CalculateSurfaceNormal(p, sample);
                    }
                }
            }
        }

        return h;
    }
    public static Vector3 CalculateSurfaceNormal(Vector3 p, SE.IsosurfaceAlgorithm.Sample sample) {
         const float H = 0.001f;
         float dx = sample(p + new Vector3(H, 0, 0)) - sample(p - new Vector3(H, 0, 0));
         float dy = sample(p + new Vector3(0, H, 0)) - sample(p - new Vector3(0, H, 0));
         float dz = sample(p + new Vector3(0, 0, H)) - sample(p - new Vector3(0, 0, H));

        return new Vector3(dx, dy, dz).normalized;
    }
    public static Vector3 ApproximateZeroCrossingPosition(Vector3 p0, Vector3 p1, SE.IsosurfaceAlgorithm.Sample sample)
    {
        // approximate the zero crossing by finding the min value along the edge
        float minValue = 100000;
        float t = 0;
        float currentT = 0;
        const int steps = 8;
        const float increment = 1f / (float)steps;
        while (currentT <= 1f)
        {
            Vector3 p = p0 + ((p1 - p0) * currentT);
            float density = Mathf.Abs(sample(p));
            if (density < minValue)
            {
                minValue = density;
                t = currentT;
            }

            currentT += increment;
        }

        return p0 + ((p1 - p0) * t);
    }