/* Copyright (c) 2010, Chris Want
   Algorithm adapted from the applet 'base26' by Karsten Schmidt
   Please see BSD style licence at the end of this file */

class IsoSurface {
  boolean zUp=false;
  color col = color(128, 128, 128);
  float opMax = 128.0, opMin = 0.0;

  ScalarField scalars;
  int activeScalars;
  
  // space dimensions
  int gridx, gridy;
  float isoStep;
  float isoVal;
  int w2;
  // vertice cache
  float[] tempX,tempY,tempZ;
  float[] cacheX,cacheY,cacheZ;
  // current list of nodes/voxels in the space
  
  // internally used lookup tables by the marching square algorithm
  // see comments in renderSlice() method for more detail
  float[][] _sliceData;
  boolean[][] _edgeFlags;
  
  int[][] _edges = {{0, 1}, {1, 2}, {2, 3}, {3, 0}};
  int[][] _offsets = {{0, 0}, {1, 0}, {1, 1}, {0, 1}};
  int[][] _lines = {
    {-1, -1, -1, -1, -1},
    {0, 3, -1, -1, -1},
    {0, 1, -1, -1, -1},
    {3, 1, -1, -1, -1},
    {1, 2, -1, -1, -1},
    {1, 2, 0, 3, -1},
    {0, 2, -1, -1, -1},
    {3, 2, -1, -1, -1},
    {3, 2, -1, -1, -1},
    {0, 2, -1, -1, -1},
    {3, 2, 0, 2, -1},
    {1, 2, -1, -1, -1},
    {3, 1, -1, -1, -1},
    {0, 1, -1, -1, -1},
    {0, 3, -1, -1, -1},
    {-1, -1, -1, -1, -1}
  };
  // vertex counter
  int vCount;
  
  IsoSurface() {
    scalars = new ScalarField();

    initialize();
  }

  IsoSurface(ScalarField myScalars) {
    setScalars(myScalars);
  }

  IsoSurface(ScalarField myScalars, int arr) {
    setScalars(myScalars, arr);
  }

  IsoSurface(String filename) {
    setScalars(new ScalarField(filename));
  }

  void setScalars(ScalarField myScalars) {
    setScalars(myScalars, myScalars.getActiveScalars());
  }
  
  void setScalars(ScalarField myScalars, int arr) {
    scalars = myScalars;
    activeScalars = arr;

    initialize();
  }

  void initialize() {
    int oldActive = scalars.getActiveScalars();
    scalars.setActiveScalars(activeScalars);

    gridx = scalars.xRes;
    gridy = max(scalars.yRes,scalars.zRes);
    _sliceData = new float[gridx][gridy];
    _edgeFlags = new boolean[gridx][gridy];

    cacheX=null;
    cacheY=null;
    cacheZ=null;
    
    w2 = (scalars.xRes >> 1);

    isoVal  = (scalars.getScalMin() + scalars.getScalMax()) / 2.0;
    isoStep = (scalars.getScalMax() - scalars.getScalMin()) / 100.0;

    scalars.setActiveScalars(oldActive);
  }
  
  // compute new surface hull
  // s=step size
  // isoV = threshold value to be used for finding surface
  void increaseIso(int times) {
    isoVal += isoStep*times;
    if (isoVal > scalars.getScalMax(activeScalars)) {
      isoVal = scalars.getScalMax(activeScalars);
    }
  }

  void decreaseIso(int times) {
    isoVal -= isoStep*times;
    if (isoVal < scalars.getScalMin(activeScalars)) {
      isoVal = scalars.getScalMin(activeScalars);
    }
  }

  void update() {
    update(isoVal);
  }

  void updateUnit(float unit) {
    update(scalars.getScalMax(activeScalars) * unit + 
           scalars.getScalMin(activeScalars) * (1.0 - unit));
  }
  
  void update(float isoV) {
    update(isoV, false);
  }
  
  void updateInvalid() {
    update(0.5, true);
  }

  void update(float isoV, boolean invalidMode) {

    isoVal = isoV;
    vCount=0;
    cacheX=null;
    cacheY=null;
    cacheZ=null;

    if (invalidMode && !scalars.hasInvalid) return;
    
    int oldActive = scalars.getActiveScalars();
    scalars.setActiveScalars(activeScalars);

    tempX = new float[8*(gridx-1)*(gridy-1)];
    tempY = new float[8*(gridx-1)*(gridy-1)];
    tempZ = new float[8*(gridx-1)*(gridy-1)];
    for(int i=0; i < scalars.zRes; ++i) {
      computeSlice(true,i,isoV, invalidMode);
    }
    if (scalars.zRes > 1) {
      for(int i=0; i < scalars.yRes; ++i) {
        computeSlice(false,i,isoV, invalidMode);
      }
    }
    tempX=null;
    tempY=null;
    tempZ=null;

    scalars.setActiveScalars(oldActive);
    
  }
  
  void zUpOn() {
    zUp = true;
  }

  void zUpOff() {
    zUp = false;
  }

  void setColor(int r, int g, int b) {
    col = color(r, g, b);
  }

  void setColor(color c) {
    col = c;
  }
  
  void setMinOpacity(float val) {
    opMin = val;
  }

  void setMaxOpacity(float val) {
    opMax = val;
  }
  
  void render() {
    render(g);
  }

  // render current contents of vertex cache
  void render(PGraphics pg) {
    float a, b, opacity;
    if (zUp) {
      if (scalars.zRes > 1) {
        a = (opMax - opMin) / (scalars.zMax - scalars.zMin);
        b = opMin - scalars.zMin * a;
      }
      else {
        a = 0;
        b = opMax;
      }
    }
    else {
      a = (opMax - opMin) / (scalars.yMax - scalars.yMin);
      b = opMin - scalars.yMin * a;
    }
    pg.beginShape(LINES);
    for(int i=0; i<vCount; i++) {
      if (zUp) {
        opacity = a*cacheZ[i] + b;
      }
      else {
        opacity = a*cacheY[i] + b;
      }
      pg.stroke(col, opacity);
      pg.vertex(cacheX[i],cacheY[i],cacheZ[i]);
    }
    pg.endShape();
  }
  
  void computeSlice(boolean zSlice, int currSlice, float isoV) {
    computeSlice(zSlice, currSlice, isoV, false);
  }
    // marching square implementation based on code/descriptions by paul bourke
  void computeSlice(boolean zSlice, int currSlice, float isoV, boolean invalidMode) {
    // those values are described below
    int[] startP;
    int[] endP;
    int startOffsetX;
    int startOffsetY;
    int endOffsetX;
    int endOffsetY;
    int corners;
    int[] currEdges;
    int yy,n;
    float v1;
    float dV;
    float lerpF;
    float vx1,vy1,vx2,vy2,vz;
    int tempCount = 0;
    float xMin, xSlope, yMin, ySlope, zMin, zSlope;
    
    if (invalidMode) isoV = 0.5;
    
    for (int x = 0; x < gridx; x++) {
      for (int y = 0; y < gridy; y++) {
          // reset all edges to unused
          _edgeFlags[x][y] = false;
      }
    }

    // main slice rendering pass
    // re-process values in _sliceData array and determine if
    // the field strength at each point is less than the given
    // threshold value. we have found an intersection with the 
    // boundary surface, if the current point's value is below,
    // but one of it's neighbour is >= threshold value
    xMin = scalars.xMin;
    xSlope = (scalars.xMax - scalars.xMin) / (scalars.xRes-1);
    
    if (zSlice) {
      yMin = scalars.yMin;
      ySlope = (scalars.yMax - scalars.yMin) / (scalars.yRes-1);
      zMin = scalars.zMin;
      if (scalars.zRes > 1) {
        zSlope = (scalars.zMax - scalars.zMin) / (scalars.zRes-1);
      }
      else {
        zSlope = 0.0;
      }
      for (int i=0; i < scalars.xRes; ++i) {
        for (int j=0; j < scalars.yRes; ++j) {
          if (!invalidMode) {
            _sliceData[i][j] = scalars.getData(i, j, currSlice);
          }
          else {
            if (scalars.isInvalid(i, j, currSlice)) {
              _sliceData[i][j] = 1.0;
            }
            else {
              _sliceData[i][j] = 0.0;
            }
          }
        }
      }
    }
    else {
      yMin = scalars.zMin;
      ySlope = (scalars.zMax - scalars.zMin) / (scalars.zRes - 1);
      zMin = scalars.yMin;
      zSlope = (scalars.yMax - scalars.yMin) / (scalars.yRes - 1);

      for (int i=0; i < scalars.xRes; ++i) {
        for (int j=0; j < scalars.zRes; ++j) {
          if (!invalidMode) {
            _sliceData[i][j] = scalars.getData(i, currSlice, j);
          }
          else {
            if (scalars.isInvalid(i, currSlice, j)) {
              _sliceData[i][j] = 1.0;
            }
            else {
              _sliceData[i][j] = 0.0;
            }
          }
        } 
      }
    }
    vz  = zMin + currSlice * zSlope;

    //beginShape(LINES);
    int gX1=scalars.xRes-1;
    int gY1=(zSlice ? scalars.yRes-1 : scalars.zRes-1);
    for (int y = 0; y < gY1; y++) {
      yy=y+1;
      for (int x = 0,xx=1; x < gX1; x++) {
        // check if this edge hasn't already been processed
        if (!_edgeFlags[x][y]) {
          corners = 0;
          // compare all 4 corners of the current grid square/cell with
          // the iso threshold value and set respective codes
          // 1=bottom left
          // 2=bottom right
          // 4=top right
          // 8=top left
          if (_sliceData[x][y] < isoV) {
            corners |= 1;
          }
          if (_sliceData[xx][y] < isoV) {
            corners |= 2;
          }
          if (_sliceData[xx][yy] < isoV) {
            corners |= 4;
          }
          if (_sliceData[x][yy] < isoV) {
            corners |= 8;
          }
          
          // if the result of corners=0 the entire square is within the surface area
          // if it is 1+2+4+8=15 the square is totally outside
          // in both cases no further steps are needed and we can skip that next part
          
          if (corners > 0 && corners < 0x0f) {

            n = 0;
            // here the algorithm makes heavy use of the lookup tables (defined above)
            // to compute the exact position of the surface line(s) in the current cell and at same
            // time take care of all possible symmetries.
            // note that there're 2 special cases where there will be 2 lines crossing the cell
            // this only happens when diagonally opposite points are within the surface and the other 2 outside
            currEdges=_lines[corners];
            // a value of -1 in the currEdges array means no more edges to check
            while (currEdges[n] != -1) {
              // retrieve edges involved
              int[] edge1 = _edges[currEdges[n++]];
              int[] edge2 = _edges[currEdges[n++]];
              
              // get indexes of start/end point relative to current grid xy position
              startP = _offsets[edge1[0]];
              endP = _offsets[edge1[1]];
              // get absolut grid position of those points by adding current xy's
              startOffsetX = x + startP[0];
              startOffsetY = y + startP[1];
              endOffsetX = x + endP[0];
              endOffsetY = y + endP[1];
              if (!invalidMode && scalars.hasInvalid) {
                float bad = scalars.getInvalidValue();
                if ( (_sliceData[startOffsetX][startOffsetY] == bad) ||
                     (_sliceData[endOffsetX][endOffsetY] == bad) ) {
                       _edgeFlags[x][y] = true;
                       continue;
                }
              }
              // get field value for the start point
              v1 = _sliceData[startOffsetX][startOffsetY];
              // compute difference in field strength to end point
              dV = _sliceData[endOffsetX][endOffsetY] - v1;
              // get linear interpolation factor based difference in field strengths
              lerpF = (dV != 0) ? (isoV - v1) / dV : 0.5;
              // compute interpolated position of start vertex
              vx1 = startOffsetX + lerpF * (endOffsetX - startOffsetX);
              vy1 = startOffsetY + lerpF * (endOffsetY - startOffsetY);
              
              // now do the same thing for the end point of the line
              startP = _offsets[edge2[0]];
              endP = _offsets[edge2[1]];
              startOffsetX = x + startP[0];
              startOffsetY = y + startP[1];
              endOffsetX = x + endP[0];
              endOffsetY = y + endP[1];
              v1 = _sliceData[startOffsetX][startOffsetY];
              dV = _sliceData[endOffsetX][endOffsetY] - v1;
              lerpF = (dV != 0) ? (isoV - v1) / dV : 0.5;
              vx2 = startOffsetX + lerpF * (endOffsetX - startOffsetX);
              vy2 = startOffsetY + lerpF * (endOffsetY - startOffsetY);
              
              vx1 = xMin + vx1 * xSlope;
              vx2 = xMin + vx2 * xSlope;
              vy1 = yMin + vy1 * ySlope;
              vy2 = yMin + vy2 * ySlope;

              if (zSlice) {
                tempX[tempCount]=vx1;
                tempY[tempCount]=vy1;
                tempZ[tempCount++]=vz;
                tempX[tempCount]=vx2;
                tempY[tempCount]=vy2;
                tempZ[tempCount++]=vz;
              }
              else {
                tempX[tempCount]=vx1;
                tempY[tempCount]=vz;
                tempZ[tempCount++]=vy1;
                tempX[tempCount]=vx2;
                tempY[tempCount]=vz;
                tempZ[tempCount++]=vy2;
              }

              // set edge as already processed
              _edgeFlags[x][y] = true;
            }
          }
        }
        xx++;
      }
    }
    if (tempCount > 0) {
      float[] newX = new float[vCount+tempCount];
      float[] newY = new float[vCount+tempCount];
      float[] newZ = new float[vCount+tempCount];

      for (int i=0; i<vCount; ++i) {
        newX[i] = cacheX[i];
        newY[i] = cacheY[i];
        newZ[i] = cacheZ[i];
      }
      for (int i=0; i<tempCount; ++i) {
        newX[vCount+i] = tempX[i];
        newY[vCount+i] = tempY[i];
        newZ[vCount+i] = tempZ[i];
      }
      cacheX = newX;
      cacheY = newY;
      cacheZ = newZ;
      vCount += tempCount;
    }
  }

  void destroy() {
    scalars = null;
    _sliceData = null;
    _edgeFlags = null; 

    cacheX = null;
    cacheY = null;
    cacheZ = null;

    
  }
}

/* Copyright (c) 2010, Chris Want, Research Support Group,
   AICT, University of Alberta. All rights reserved.

  Redistribution and use in source and binary forms, with or without 
  modification, are permitted provided that the following conditions are met:

  1) Redistributions of source code must retain the above copyright 
     notice, this list of conditions and the following disclaimer.
  2) Redistributions in binary form must reproduce the above copyright 
     notice, this list of conditions and the following disclaimer in the 
     documentation and/or other materials provided with the distribution.

  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF  
  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF 
  THE POSSIBILITY OF SUCH DAMAGE.

  Contributors: Chris Want (University of Alberta),
  Karsten Schmidt
*/