/* * Conrec.java * * Created on 5 August 2001, 15:03 * * Copyright (c) 1996-1997 Nicholas Yue * * This software is copyrighted by Nicholas Yue. This code is base on the work of * Paul D. Bourke CONREC.F routine * * The authors hereby grant permission to use, copy, and distribute this * software and its documentation for any purpose, provided that existing * copyright notices are retained in all copies and that this notice is included * verbatim in any distributions. Additionally, the authors grant permission to * modify this software and its documentation for any purpose, provided that * such modifications are not distributed without the explicit consent of the * authors and that existing copyright notices are retained in all copies. Some * of the algorithms implemented by this software are patented, observe all * applicable patent law. * * IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO ANY PARTY FOR * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT * OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY DERIVATIVES THEREOF, * EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * * THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING, * BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A * PARTICULAR PURPOSE, AND NON-INFRINGEMENT. THIS SOFTWARE IS PROVIDED ON AN * "AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE NO OBLIGATION TO PROVIDE * MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. */ /** * Conrec a straightforward method of contouring some surface represented a regular * triangular mesh. * * Ported from the C++ code by Nicholas Yue (see above copyright notice). * @see http://paulbourke.net/papers/conrec for full description * of code and original C++ source. * * @author Bradley White * @version 1.0 */ public class Conrec { private double [] h = new double [5]; private int [] sh = new int [5]; private double [] xh = new double [5]; private double [] yh = new double [5]; // Object that knows how to draw the contour private Render render = null; /** Creates new Conrec */ public Conrec(Render render) throws Exception { if (render == null){ throw new Exception ("Render null"); } this.render = render; } /** * contour is a contouring subroutine for rectangularily spaced data * * It emits calls to a line drawing subroutine supplied by the user * which draws a contour map corresponding to real*4data on a randomly * spaced rectangular grid. The coordinates emitted are in the same * units given in the x() and y() arrays. * * Any number of contour levels may be specified but they must be * in order of increasing value. * * * @param d - matrix of data to contour * @param ilb,iub,jlb,jub - index bounds of data matrix * * The following two, one dimensional arrays (x and y) contain the horizontal and * vertical coordinates of each sample points. * @param x - data matrix column coordinates * @param y - data matrix row coordinates * @param nc - number of contour levels * @param z - contour levels in increasing order. * */ public void contour(double [][] d, int ilb, int iub, int jlb, int jub, double [] x, double [] y, int nc, double [] z) { int m1; int m2; int m3; int case_value; double dmin; double dmax; double x1 = 0.0; double x2 = 0.0; double y1 = 0.0; double y2 = 0.0; int i,j,k,m; // The indexing of im and jm should be noted as it has to start from zero // unlike the fortran counter part int [] im = {0,1,1,0}; int [] jm = {0,0,1,1}; // Note that castab is arranged differently from the FORTRAN code because // Fortran and C/C++ arrays are transposed of each other, in this case // it is more tricky as castab is in 3 dimension int [][][] castab= { { {0,0,8},{0,2,5},{7,6,9} }, { {0,3,4},{1,3,1},{4,3,0} }, { {9,6,7},{5,2,0},{8,0,0} } }; for (j=(jub-1);j>=jlb;j--) { for (i=ilb;i<=iub-1;i++) { double temp1,temp2; temp1 = Math.min(d[i][j],d[i][j+1]); temp2 = Math.min(d[i+1][j],d[i+1][j+1]); dmin = Math.min(temp1,temp2); temp1 = Math.max(d[i][j],d[i][j+1]); temp2 = Math.max(d[i+1][j],d[i+1][j+1]); dmax = Math.max(temp1,temp2); if (dmax>=z[0]&&dmin<=z[nc-1]) { for (k=0;k=dmin&&z[k]<=dmax) { for (m=4;m>=0;m--) { if (m>0) { // The indexing of im and jm should be noted as it has to // start from zero h[m] = d[i+im[m-1]][j+jm[m-1]]-z[k]; xh[m] = x[i+im[m-1]]; yh[m] = y[j+jm[m-1]]; } else { h[0] = 0.25*(h[1]+h[2]+h[3]+h[4]); xh[0]=0.5*(x[i]+x[i+1]); yh[0]=0.5*(y[j]+y[j+1]); } if (h[m]>0.0) { sh[m] = 1; } else if (h[m]<0.0) { sh[m] = -1; } else sh[m] = 0; } // // Note: at this stage the relative heights of the corners and the // centre are in the h array, and the corresponding coordinates are // in the xh and yh arrays. The centre of the box is indexed by 0 // and the 4 corners by 1 to 4 as shown below. // Each triangle is then indexed by the parameter m, and the 3 // vertices of each triangle are indexed by parameters m1,m2,and // m3. // It is assumed that the centre of the box is always vertex 2 // though this isimportant only when all 3 vertices lie exactly on // the same contour level, in which case only the side of the box // is drawn. // // // vertex 4 +-------------------+ vertex 3 // | \ / | // | \ m-3 / | // | \ / | // | \ / | // | m=2 X m=2 | the centre is vertex 0 // | / \ | // | / \ | // | / m=1 \ | // | / \ | // vertex 1 +-------------------+ vertex 2 // // // // Scan each triangle in the box // for (m=1;m<=4;m++) { m1 = m; m2 = 0; if (m!=4) { m3 = m+1; } else { m3 = 1; } case_value = castab[sh[m1]+1][sh[m2]+1][sh[m3]+1]; if (case_value!=0) { switch (case_value) { case 1: // Line between vertices 1 and 2 x1=xh[m1]; y1=yh[m1]; x2=xh[m2]; y2=yh[m2]; break; case 2: // Line between vertices 2 and 3 x1=xh[m2]; y1=yh[m2]; x2=xh[m3]; y2=yh[m3]; break; case 3: // Line between vertices 3 and 1 x1=xh[m3]; y1=yh[m3]; x2=xh[m1]; y2=yh[m1]; break; case 4: // Line between vertex 1 and side 2-3 x1=xh[m1]; y1=yh[m1]; x2=xsect(m2,m3); y2=ysect(m2,m3); break; case 5: // Line between vertex 2 and side 3-1 x1=xh[m2]; y1=yh[m2]; x2=xsect(m3,m1); y2=ysect(m3,m1); break; case 6: // Line between vertex 3 and side 1-2 x1=xh[m3]; y1=yh[m3]; x2=xsect(m1,m2); y2=ysect(m1,m2); break; case 7: // Line between sides 1-2 and 2-3 x1=xsect(m1,m2); y1=ysect(m1,m2); x2=xsect(m2,m3); y2=ysect(m2,m3); break; case 8: // Line between sides 2-3 and 3-1 x1=xsect(m2,m3); y1=ysect(m2,m3); x2=xsect(m3,m1); y2=ysect(m3,m1); break; case 9: // Line between sides 3-1 and 1-2 x1=xsect(m3,m1); y1=ysect(m3,m1); x2=xsect(m1,m2); y2=ysect(m1,m2); break; default: break; } // Put your processing code here and comment out the printf //printf("%f %f %f %f %f\n",x1,y1,x2,y2,z[k]); render.drawContour(x1,y1,x2,y2,z[k]); } } } } } } } } private double xsect(int p1, int p2){ return (h[p2]*xh[p1]-h[p1]*xh[p2])/(h[p2]-h[p1]); } private double ysect(int p1, int p2){ return (h[p2]*yh[p1]-h[p1]*yh[p2])/(h[p2]-h[p1]); } }