package nxm.sys.libg;

/** A class to create high precision contour lines from a 2D mesh.
    This code is Contour with doubles instead of floats.
    <pre>
    The general process is:
    1) Process the exterior grid squares first (left, right, top, bottom).
    2) Once an entry point is found, doLine processes that grid block.
    3) doLine figures out the exit point, which grid the line goes into next,
       where it exits that grid block, ... until the line goes out-of-bounds.
    4) doLine then returns the line segment along with header information
    3) As each grid square is checked, the chk array is updated, so the
       square does not need to be checked again.
    4) After exterior lines/squares are checked, interior ones are checked.
    </pre>

    Based on {@link Contour}
    @author Jeff Schoen, Jeanette McCall
 */
public class ContourHighPrecision {

  /** left */      public static final int P_L=1;
  /** top */       public static final int P_T=2;
  /** right */     public static final int P_R=3;
  /** bottom */    public static final int P_B=4;

  /** x start */   public double xstart;
  /** x delta */   public double xdelta;
  /** y start */   public double ystart;
  /** y delta */   public double ydelta;
  /** z start */   public double zstart;
  /** z delta */   public double zdelta;
  /** number x */  public int nx=0;
  /** number y */  public int ny=0;
  /** number z */  public int nz=0;
  /** number xy */ public int nxy=0;

  /** maximum number output */
  public int maxout=0;
  /** number output */
  public int nout;
  private double work[],z;
  /** whether this grid block has been checked */
  private int chk[];
  private int ix1,ix2,iy1,iy2;
  /** output buffer that holds contour line as x1,y1,x2,y2,... where each pair
   *  is a location where the value matches z (ex: horizon line just has z of 0)*/
  public double out[];
  /** For each line the first 6 values are npts, z, xmin, xmax, yim, ymax */
  private static final int NUM_HEADER_VALUES = 6;

  /** Constructor. */
  public ContourHighPrecision () {
  }

  /** sets the x-axis
      @param start x-axis start value
      @param delta x-axis delta value
      @param n     number of x values
   */
  public void setXAxis (double start, double delta, int n) {
    xstart=start; xdelta=delta; nx=n;
    setBounds(0,0,nx,ny);
  }
  /** sets the y-axis
      @param start y-axis start value
      @param delta y-axis delta value
      @param n     number of y values
  */
  public void setYAxis (double start, double delta, int n) {
    ystart=start; ydelta=delta; ny=n;
    setBounds(0,0,nx,ny);
  }
  /** sets the z-axis
      @param start z-axis start value
      @param delta z-axis delta value
      @param n     number of z values
   */
  public void setZAxis (double start, double delta, int n) {
    zstart=start; zdelta=delta; nz=n;
  }

  /** Sets the x-axis and y-axis bounds
      @param ix1 first x-coordinate
      @param iy1 first y-coordinate
      @param ix2 second x-coordinate
      @param iy2 second y-coordinate
   */
  public void setBounds (int ix1, int iy1, int ix2, int iy2) {
    this.ix1=Math.max(0,Math.min(nx,Math.min(ix1,ix2)));
    this.ix2=Math.max(0,Math.min(nx,Math.max(ix1,ix2)));
    this.iy1=Math.max(0,Math.min(ny,Math.min(iy1,iy2)));
    this.iy2=Math.max(0,Math.min(ny,Math.max(iy1,iy2)));
  }

  /** process this buffer of data and place values in contour.out
   * <pre>
   *  [npts 1st group,zmax,xstart,ystart,xmax,ymax,lon,lat,lon,lat,....,
   *  npts 2nd group,zmax,xstart,ystart,xmax,ymax,lon,lat,lon,lat,.....,
   *  ...
   *  ]
   *  </pre>
   *  With a separate set of points for each contour segment, that
   *  enters/exits the bounds of the plot.
   *
   *  @param buf the buffer of data to process
   *  @return number of valid values in contour.out
   */
  public int process (double[] buf) {
    int ix,iy,iz,ii;
    work = buf;
    if (nx*ny*4 > maxout) {
      maxout = nx*ny*4;
      out = new double[maxout];
    }
    if (nx*ny > nxy) {
      nxy = nx*ny;
      chk = new int[nxy];
    }
    nout=0;
    for (iz=0; iz<nz; iz++) {
      for (iy=iy1; iy<iy2; iy++) {
        ii = iy*nx + ix1;
        for (ix=ix1; ix<ix2; ix++) chk[ii++] = 0;
      }
      z = zstart + iz*zdelta;
      // check left and right edges
      if (ix2>ix1) for (iy=iy1; iy<iy2-1; iy++) {
        ii = iy*nx+ix1;
        if (chk[ii]==0 && (z>work[ii])^(z>work[ii+nx])) doLine(ix1,iy,P_L);
        ii += (ix2-ix1-1);
        if (chk[ii-1]==0 && (z>work[ii])^(z>work[ii+nx])) doLine(ix2-2,iy,P_R);
      }
      // check top and bottom edges
      if (iy2>iy1) for (ix=ix1; ix<ix2-1; ix++) {
        ii = ix+iy1*nx;
        if (chk[ii]==0 && (z>work[ii])^(z>work[ii+1])) doLine(ix,iy1,P_T);
        ii += nx*(iy2-iy1-1); //= ix + nx*iy2 -nx
        if (chk[ii-nx]==0 && (z>work[ii])^(z>work[ii+1])) doLine(ix,iy2-2,P_B);
      }
      // check top edge of all internal cells
      for (iy=iy1+1; iy<iy2-1; iy++) {
        ii = iy*nx+ix1+1;
        for (ix=ix1+1; ix<ix2-1; ix++,ii++) {
          if ((chk[ii]==0) &&
              ((z>work[ii])^(z>work[ii+1]))) doLine(ix,iy,P_T);
        }
      }
    }
    return nout;
  }

  /** State array that determines exit side.
   *  This has 1 row for each side: s1/P_L s2/P_T s3 P_R s4/P_B<br>
   *  Each row is symmetrical, so only the first 8 values are listed.
   *  <pre>
   *  1) If the corner values on one side of the side are greater than z and
   *     the and the other side of the side has both values less than z, then
   *     just exit on the opposite side.
   *     Example: We are on s1 and z1 and z2 are less than z, and z3 and z4
   *     are greater than z, then we exit on s3.
   *  2) If (z - corner_z) of 3 corners is of one sign and is the opposite sign
   *     for the other corner we exit on the side that is next to the corner
   *     with the different z-corner sign (i.e. 3 positive 1 negative, or
   *     3 negative 1 positive)
   *     Example: We are on s2 and z1 is greater than z and the other corners
   *     are less than z, so we exit on s1
   *
   *       z1      s2       z2
   *         +-------------+
   *         |             |
   *         |             |
   *         |             |
   *      s1 |             | s3
   *         |             |
   *         |             |
   *         |             |
   *         +-------------+
   *       z4      s4       z3
   *  </pre>
   */
  private static final int statearray[] = {
    0,4,0,3,0,2,0,2,
    0,0,0,0,3,1,4,1,
    0,0,4,1,2,4,0,0,
    0,1,3,0,0,3,2,0 };

  private void doLine (int ix, int iy, int inside) {
    double z1;
    double z2;
    double z3;
    double z4;
    double rr; // ratio down or across grid block you need to go to have value=z
    int side = inside;
    int state, newside;
    int ii = ix+nx*iy;
    boolean inbounds=true;
    int n = nout+NUM_HEADER_VALUES;

    while (nout<maxout && inbounds && chk[ii]!=side) {
      chk[ii]=side;
      z1=work[ii];
      z2=work[ii+1];
      z3=work[ii+1+nx];
      z4=work[ii+nx];
      // calculate entry point
      if (n==nout+NUM_HEADER_VALUES) switch (side) {
        case P_L:
          rr = (z-z1)/(z4-z1); // ratio down left side
          out[n++]=xstart+(ix)*xdelta;
          out[n++]=ystart+(iy+rr)*ydelta;
          break;
        case P_T:
          rr = (z-z1)/(z2-z1); // ratio across top
          out[n++]=xstart+(ix+rr)*xdelta;
          out[n++]=ystart+(iy)*ydelta;
          break;
        case P_R:
          rr = (z-z2)/(z3-z2); // ratio down right side
          out[n++]=xstart+(ix+1)*xdelta;
          out[n++]=ystart+(iy+rr)*ydelta;
          break;
        case P_B:
          rr = (z-z4)/(z3-z4); // ratio across bottom
          out[n++]=xstart+(ix+rr)*xdelta;
          out[n++]=ystart+(iy+1)*ydelta;
          break;
      }
      // determine relationship of Z to corners
      state = 1;
      if (z>z1) state+=8;
      if (z>z2) state+=4;
      if (z>z3) state+=2;
      if (z>z4) state+=1;
      // the state array is symmetrical
      if (state>8) state = 17-state;
      // decide where to exit
      newside = statearray[ state-1 + (side-1)*8 ];
      // calculate exit point
      switch (newside) {
        case P_L:
          rr = (z-z1)/(z4-z1); // ratio down left side
          out[n++]=xstart+(ix)*xdelta;
          out[n++]=ystart+(iy+rr)*ydelta;
          side = P_R; inbounds = (--ix >= ix1);
          break;
        case P_T:
          rr = (z-z1)/(z2-z1); // ratio across top
          out[n++]=xstart+(ix+rr)*xdelta;
          out[n++]=ystart+(iy)*ydelta;
          side = P_B; inbounds = (--iy >= iy1);
          break;
        case P_R:
          rr = (z-z2)/(z3-z2); // ratio down right side
          out[n++]=xstart+(ix+1)*xdelta;
          out[n++]=ystart+(iy+rr)*ydelta;
          side = P_L; inbounds = (++ix < ix2-1);
          break;
        case P_B:
          rr = (z-z4)/(z3-z4); // ratio across bottom
          out[n++]=xstart+(ix+rr)*xdelta;
          out[n++]=ystart+(iy+1)*ydelta;
          side = P_T; inbounds = (++iy < iy2-1);
          break;
      }
      ii = ix+nx*iy;
    }
    // populate trace segment (i.e. line) header
    int npts =  (n-nout-NUM_HEADER_VALUES)/2;
    if (npts<=1) return;
    out[nout+0] = npts;
    out[nout+1] = z;
    if (xdelta > 0) {
      out[nout+2] = xstart;
      out[nout+4] = xstart+(nx-1)*xdelta;
    } else {
      out[nout+4] = xstart;
      out[nout+2] = xstart+(nx-1)*xdelta;
    }
    if (ydelta > 0) {
      out[nout+3] = ystart;
      out[nout+5] = ystart+(ny-1)*ydelta;
    } else {
      out[nout+5] = ystart;
      out[nout+3] = ystart+(ny-1)*ydelta;
    }
    nout = n;
  }

}
