/*
* ContourGenerator.java
*
* Copyright (c) 2002-2015 Alexei Drummond, Andrew Rambaut and Marc Suchard
*
* This file is part of BEAST.
* See the NOTICE file distributed with this work for additional
* information regarding copyright ownership and licensing.
*
* BEAST is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* BEAST is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with BEAST; if not, write to the
* Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
* Boston, MA 02110-1301 USA
*/
package dr.geo.contouring;
import dr.geo.contouring.ContourAttrib;
import java.util.*;
/**
* <p> An object used to generate a list of contour lines
* or paths from a set of gridded three dimensional data.
* </p>
*
* <p> Based on contour_plot.c from NeXTcontour1.4 by Thomas H. Pulliam,
* pulliam@rft29.nas.nasa.gov, MS 202A-1 NASA Ames Research Center,
* Moffett Field, CA 94035.
* I don't know how the original Fortran code looked like or where it came from,
* other than that NeXTcontour1.4 is based on Pieter Bunings' PLOT3D package
* for Computational Fluid Dynamics.
* </p>
*
* <p> Ported from C to Java by Joseph A. Huwaldt, November 16, 2000. </p>
*
* <p> Modified by: Joseph A. Huwaldt </p>
*
* @author Joseph A. Huwaldt Date: November 11, 2000
* @version November 23, 2000
*
* @author Marc Suchard
*
**/
public class ContourGenerator {
// Debug flag.
private static final boolean DEBUG = false;
// Error messages.
private static final String kCancelMsg = "Method ContourGenerator.getContours() canceled by user.";
private static final String kInconsistantArrMsg = "Inconsistant array sizes.";
private static final String kArrSizeMsg = "Data arrays must have more than one row or column.";
private static final String kNegLogDataMsg = "Function data must be > 0 for logarithmic intervals.";
// Path buffer size.
private static final int kBufSize = 1000;
// The minimum number of points allowed in a contour path.
private static final int kMinNumPoints = 3;
// A list of contour paths.
private List pathList = new ArrayList();
// A flag to indicate that the contours have been computed or not.
private boolean cCalculated = false;
// Data arrays used for generating the contours.
private double[][] xArray, yArray, funcArray;
// Data arrays used when generating contours for 1D X & Y arrays.
private double[] xArr1D, yArr1D;
// Array of contour attributes, one for each contour level.
private ContourAttrib[] cAttr;
// The fraction of the task that is completed.
private float fracComplete = 0;
/**
* Used to indicate that the user wishes to cancel the calculation
* of contours.
**/
private boolean isCanceled = false;
// Variables in the original FORTRAN program.
private double[] pathbufxt, pathbufyt;
private int[] pathbufia;
private int lnstrt; // lnstrt=1 indicates starting a new line.
private int ignext;
private int icont; // Current contour level index.
private double cont; // The current contour level.
private int iss, iee, jss, jee; // i & j start and end index values.
private int ima; // ima tells which boundary region we are on.
private int iae; // Index to last element in the IA list.
private int ibeg, jbeg;
private int gi, gj; // Indexes into data arrays.
private double fij; // Data value at i,j in data array.
private int idir; // Indicates current direction.
private int np=0; // Number of points in current contour line.
private double wx=0, wy=0; // Starting point of a contour line.
/**
* Construct a ContourGenerator object using the specified data arrays
* and the specified attribute array. This constructor allows you
* to use data on an uneven X, Y grid.
*
* @param xArr 2D array containing the grid x coordinate data.
* @param yArr 2D array containing the grid y coordinate data.
* @param fArr 2D array containing the grid function (z) data.
* @param cAttr Array containing attributes of the contour levels.
**/
public ContourGenerator(double[][] xArr, double[][] yArr, double[][] fArr, ContourAttrib[] cAttr) {
// Make sure input data is reasonable.
if (yArr.length != xArr.length || yArr.length != fArr.length)
throw new IllegalArgumentException(kInconsistantArrMsg);
if (yArr[0].length != xArr[0].length || yArr[0].length != fArr[0].length)
throw new IllegalArgumentException(kInconsistantArrMsg);
if (xArr.length <= 1 || xArr[0].length <= 1)
throw new IllegalArgumentException(kArrSizeMsg);
this.cAttr = cAttr;
xArray = xArr;
yArray = yArr;
funcArray = fArr;
}
/**
* Construct a ContourGenerator object using the specified data arrays
* and the specified attribute array. This constructor allows you
* to use data on an evenly spaced grid where "X" values are invarient
* with "Y" and "Y" values are invarient with "X". This often occures
* where the data is on an evenly spaced cartesian grid.
*
* @param xArr 1D array containing the grid x coordinate data.
* @param yArr 1D array containing the grid y coordinate data.
* @param fArr 2D array containing the grid function (z) data.
* @param cAttr Array containing attributes of the contour levels.
**/
public ContourGenerator(double[] xArr, double[] yArr, double[][] fArr, ContourAttrib[] cAttr) {
// Make sure input data is reasonable.
if (yArr.length != fArr.length || xArr.length != fArr[0].length)
throw new IllegalArgumentException(kInconsistantArrMsg);
if (xArr.length <= 1)
throw new IllegalArgumentException(kArrSizeMsg);
this.cAttr = cAttr;
xArr1D = xArr;
yArr1D = yArr;
funcArray = fArr;
}
/**
* Construct a ContourGenerator object using the specified data arrays.
* Contour attributes, including the interval, are generated
* automatically. This constructor allows you to use data on an
* uneven X, Y grid.
*
* @param xArr 2D array containing the grid x coordinate data.
* @param yArr 2D array containing the grid y coordinate data.
* @param fArr 2D array containing the grid function (z) data.
* @param nc The number of contour levels to generate.
* @param logInterval Uses a logarithmic contour interval if true, and
* uses a linear interval if false.
**/
public ContourGenerator(double[][] xArr, double[][] yArr, double[][] fArr,
int nc, boolean logInterval) {
// Make sure input data is reasonable.
if (yArr.length != xArr.length || yArr.length != fArr.length)
throw new IllegalArgumentException(kInconsistantArrMsg);
if (yArr[0].length != xArr[0].length || yArr[0].length != fArr[0].length)
throw new IllegalArgumentException(kInconsistantArrMsg);
if (xArr.length <= 1 || xArr[0].length <= 1)
throw new IllegalArgumentException(kArrSizeMsg);
xArray = xArr;
yArray = yArr;
funcArray = fArr;
if (logInterval)
findLogIntervals(nc);
else
findLinearIntervals(nc);
}
/**
* Construct a ContourGenerator object using the specified data arrays.
* Contour attributes, including the interval, are generated
* automatically. This constructor allows you
* to use data on an evenly spaced grid where "X" values are invarient
* with "Y" and "Y" values are invarient with "X". This often occures
* where the data is on an evenly spaced cartesian grid.
*
* @param xArr 1D array containing the grid x coordinate data.
* @param yArr 1D array containing the grid y coordinate data.
* @param fArr 2D array containing the grid function (z) data.
* @param nc The number of contour levels to generate.
* @param logInterval Uses a logarithmic contour interval if true, and
* uses a linear interval if false.
**/
public ContourGenerator(double[] xArr, double[] yArr, double[][] fArr,
int nc, boolean logInterval) {
// Make sure input data is reasonable.
if (yArr.length != fArr.length || xArr.length != fArr[0].length)
throw new IllegalArgumentException(kInconsistantArrMsg);
if (xArr.length <= 1)
throw new IllegalArgumentException(kArrSizeMsg);
xArr1D = xArr;
yArr1D = yArr;
funcArray = fArr;
if (logInterval)
findLogIntervals(nc);
else
findLinearIntervals(nc);
}
/**
* Generate the contour paths and return them as an array
* of ContourPath objects. If there is a lot of data, this method
* method may take a long time, so be patient. Progress can be
* checked from another thread by calling "getProgress()".
*
* @return An array of contour path objects.
* @throws InterruptedException if the user cancels this process
* (by calling "cancel()" from another thread).
**/
public ContourPath[] getContours() throws InterruptedException {
if (!cCalculated) {
isCanceled = false;
pathList.clear();
// Go off an compute the contour paths.
computeContours();
// Now turn loose all our data arrays to be garbage collected.
cAttr = null;
xArray = yArray = funcArray = null;
xArr1D = yArr1D = null;
// Set our "done" flags.
cCalculated = true;
fracComplete = 1;
}
// Turn our pathList into an array and return the array.
int size = pathList.size();
ContourPath[] arr = new ContourPath[size];
for (int i=0; i < size; ++i)
arr[i] = (ContourPath)pathList.get(i);
return arr;
}
/**
* Returns true if the contour generation process is done. False if it is not.
**/
public boolean done() {
return cCalculated;
}
/**
* Call this method to cancel the generation of contours.
**/
public void cancel() {
isCanceled = true;
}
/**
* Returns the progress of the currently executing contour generation
* process: 0.0 (just starting) to 1.0 (done).
**/
public float getProgress() {
return fracComplete;
}
/**
* Find contour intervals that are linearly spaced through the data.
**/
private void findLinearIntervals(int nc) {
// Find min and max Z values.
double zMin = Double.MAX_VALUE;
double zMax = -zMin;
int ni = funcArray.length;
for (int i=0; i < ni; ++i) {
int nj = funcArray[i].length;
for (int j=0; j < nj; ++j) {
double zVal = funcArray[i][j];
zMin = Math.min(zMin, zVal);
zMax = Math.max(zMax, zVal);
}
}
// Allocate memory for contour attribute array.
cAttr = new ContourAttrib[nc];
// Determine contour levels.
double delta = (zMax-zMin)/(nc+1);
for (int i=0; i < nc; i++) {
cAttr[i] = new ContourAttrib( zMin + (i+1)*delta );
if (DEBUG)
System.out.println("level[" + i + "] = " + (zMin + (i+1)*delta));
}
}
/**
* Find contour intervals that are logarithmically spaced through the data.
**/
private void findLogIntervals(int nc) {
// Find min and max Z values.
double zMin = Double.MAX_VALUE;
double zMax = -zMin;
int ni = funcArray.length;
for (int i=0; i < ni; ++i) {
int nj = funcArray[i].length;
for (int j=0; j < nj; ++j) {
double zVal = funcArray[i][j];
zMin = Math.min(zMin, zVal);
zMax = Math.max(zMax, zVal);
}
}
if (zMin < 0)
throw new IllegalArgumentException(kNegLogDataMsg);
// Allocate memory for contour attribute array.
cAttr = new ContourAttrib[nc];
// Determine contour levels.
double temp = Math.log(zMin);
double delta = (Math.log(zMax) - temp)/(nc+1);
for (int i=0; i < nc; i++)
cAttr[i] = new ContourAttrib( Math.exp(temp + (i+1)*delta) );
}
/**
* Computes contour lines for gridded data and stores information about
* those contours. The result of this routine is a list of contour lines
* or paths.
**/
private void computeContours() throws InterruptedException {
int ncont = cAttr.length; // Number of contour levels.
// Find the number of data points in "I" and "J" directions.
int nx=0, ny=0;
if (xArray != null) {
ny = xArray.length;
nx = xArray[0].length;
} else {
nx = xArr1D.length;
ny = yArr1D.length;
}
// Allocate temporary storage space for path buffers.
pathbufxt = new double[kBufSize];
pathbufyt = new double[kBufSize];
pathbufia = new int[kBufSize*3];
// lnstrt=1 (line start) means we're starting a new line.
lnstrt = 1;
ignext = 0;
// Loop through each contour level.
for (icont = 0; icont < ncont; ++icont) {
// Check to see if the user has canceled.
if (isCanceled)
throw new InterruptedException(kCancelMsg);
// Begin working on this contour level.
cont = cAttr[icont].getLevel();
iss = 1;
iee = nx;
jss = 1;
jee = ny;
boolean subDivFlg = false;
/*L110*/ do {
// Find where function increases through the contour level.
FlagContourPassings();
boolean L10flg = false;
/*L210*/ do {
if (!L10flg) {
/* Search along the boundaries for contour line starts.
* IMA tells which boundary of the region we're on.
*/
ima = 1;
ibeg = iss - 1;
jbeg = jss;
}
/*L6*/ imaLoop:
do {
if (!L10flg) {
boolean imb = false;
boolean doneFlg = false;
do {
switch(ima) {
case 1:
++ibeg;
if (ibeg == iee)
ima = 2;
break;
case 2:
++jbeg;
if (jbeg == jee)
ima = 3;
break;
case 3:
--ibeg;
if (ibeg == iss)
ima = 4;
break;
case 4:
--jbeg;
if (jbeg == jss)
ima = 5;
break;
case 5:
continue imaLoop;
}
if (funcArray[jbeg -1][ibeg -1] <= cont) {
imb = true;
doneFlg = false;
} else if (imb == true)
doneFlg = true;
} while (!doneFlg);
// Got a start point.
gi = ibeg; // x index of starting point.
gj = jbeg; // y index of starting point.
fij = funcArray[jbeg -1][ibeg -1]; // z value of starting point.
// Round the corner if necessary.
/* Look different directions to see which way the contour line
* went:
* 4
* 1-|-3
* 2
*/
switch (ima) {
case 1:
Routine_L21();
break;
case 2:
if (gj != jss) {
if (!Routine_L31())
Routine_L21();
} else
Routine_L21();
break;
case 3:
if (gi != iee) {
if (!Routine_L41())
Routine_L21();
} else {
if (!Routine_L31())
Routine_L21();
}
break;
case 4:
if (gj != jee) {
if (!Routine_L51())
Routine_L21();
} else {
if (!Routine_L41())
Routine_L21();
}
break;
case 5:
if (!Routine_L51())
Routine_L21();
break;
}
} // end if(!L10flg)
// This is the end of a contour line. After this, we'll start a
// new line.
L10flg = false;
/*L90*/ lnstrt = 1; // Contour line start flag.
ignext = 0;
accumContour(np, icont, pathbufxt, pathbufyt, cAttr[icont]);
// If we're not done looking along the boundaries,
// go look there some more.
} while (ima != 5);
// Otherwise, get the next start out of IA.
/*L91*/ if (iae != 0) {
int ntmp3 = iae;
for (int iia = 1; iia <= ntmp3; ++iia) {
if (pathbufia[iia -1] != 0) {
// This is how we start in the middle of the region, using IA.
gi = pathbufia[iia - 1]/1000;
gj = pathbufia[iia - 1] - gi*1000;
fij = funcArray[gj -1][gi -1];
pathbufia[iia - 1] = 0;
Routine_L21();
L10flg = true;
break;
}
}
}
} while ( L10flg );
/* And if there are no more of these, we're done with this region.
* If we've subdivided, update the region pointers and go back for more.
*/
subDivFlg = false;
if (iee == nx) {
if (jee != ny) {
jss = jee;
jee = ny;
subDivFlg = true;
}
} else {
iss = iee;
iee = nx;
subDivFlg = true;
}
} while (subDivFlg);
// Update progress information.
fracComplete = (float)(icont+1)/(float)(ncont);
// Loop back for the next contour level.
} // Next icont
// Turn loose temporary arrays used to generate contours.
pathbufxt = null;
pathbufyt = null;
pathbufia = null;
}
/**
* Flag points in IA where the the function increases through the contour
* level, not including the boundaries. This is so we have a list of at least
* one point on each contour line that doesn't intersect a boundary.
**/
private void FlagContourPassings() {
iae = 0;
int ntmp2 = jee - 1;
for (int j=jss + 1; j <= ntmp2; ++j) {
boolean imb = false;
int iaend = iae;
int ntmp3 = iee;
for (int i=iss; i <= ntmp3; ++i) {
if (funcArray[j -1][i -1] <= cont)
imb = true;
else if (imb == true) {
++iae;
pathbufia[iae - 1] = i*1000 + j;
imb = false;
/* Check if the IA array is full. If so, the subdividing
* algorithm goes like this: if we've marked at least one
* J row, drop back to the last completed J and call that
* the region. If we haven't even finished one J row, our
* region just extends to this I location.
*/
if (iae == kBufSize*3) {
if (j > jss + 1) {
iae = iaend;
jee = j;
} else {
// Compute minimum.
jee = Math.min(j+1, jee);
iee = i;
}
// Break out of i & j loops.
return;
}
}
} // Next i
} // Next j
}
/**
* This function represents the block of code in the original
* FORTRAN program that comes after line 21.
**/
private void Routine_L21() {
while (true) {
--gi;
if (gi < iss)
return; // Goto L90.
idir = 1;
if (funcArray[gj -1][gi -1] <= cont) {
// Wipe this point out of IA if it's in the list.
/*L52*/ if (iae != 0) {
int ij = gi*1000 + gj + 1000;
int ntmp3 = iae;
for (int iia = 1; iia <= ntmp3; ++iia) {
if (pathbufia[iia - 1] == ij) {
pathbufia[iia - 1] = 0;
break;
}
}
}
doInterpolation();
return; // Goto L90.
}
fij = funcArray[gj -1][gi -1];
if (Routine_L31()) return; // Goto L90
}
}
/**
* This function represents the block of code in the original
* FORTRAN program that comes after line 31.
**/
private boolean Routine_L31() {
--gj;
if (gj < jss)
return true;
idir = 2;
if (funcArray[gj -1][gi -1] <= cont) {
doInterpolation();
return true;
}
fij = funcArray[gj -1][gi -1];
return (Routine_L41());
}
/**
* This function represents the block of code in the original
* FORTRAN program that comes after line 41.
**/
private boolean Routine_L41() {
++gi;
if (gi > iee)
return true;
idir = 3;
if (funcArray[gj -1][gi -1] <= cont) {
doInterpolation();
return true;
}
fij = funcArray[gj -1][gi -1];
return (Routine_L51());
}
/**
* This function represents the block of code in the original
* FORTRAN program that comes after line 51.
**/
private boolean Routine_L51() {
++gj;
idir = 4;
if (gj > jee)
return true;
if (funcArray[gj -1][gi -1] <= cont) {
doInterpolation();
return true;
}
fij = funcArray[gj -1][gi -1];
return false;
}
/**
* Do interpolation for X, Y coordinates.
*
* This function represents the block of code in the original
* FORTRAN program that comes after line 60.
**/
private void doInterpolation() {
// Do interpolation for X,Y coordinates.
double func = funcArray[gj -1][gi -1];
double xyf = (cont - func)/(fij - func);
/* This tests for a contour point coinciding with a grid point. In this case
* the contour routine comes up with the same physical coordinate twice. If
* If we don't trap it, it can (in some cases significantly) increase the
* number of points in a contour line. Also, if this happens on the first
* point in a line, the second point could be misinterpreted as the end of a
* (circling) contour line.
*/
if (xyf == 0)
++ignext;
double wxx=0, wyy=0;
double xVal=0, yVal=0;
if (xArray != null) {
// We have 2D arrays for the X & Y grid points.
xVal = xArray[gj -1][gi -1];
yVal = yArray[gj -1][gi -1];
switch (idir) {
case 1: // East
wxx = xVal + xyf*(xArray[gj -1][gi + 1 -1] - xVal);
wyy = yVal + xyf*(yArray[gj -1][gi + 1 -1] - yVal);
break;
case 2: // North
wxx = xVal + xyf*(xArray[gj + 1 -1][gi -1] - xVal);
wyy = yVal + xyf*(yArray[gj + 1 -1][gi -1] - yVal);
break;
case 3: // West
wxx = xVal + xyf*(xArray[gj -1][gi - 1 -1] - xVal);
wyy = yVal + xyf*(yArray[gj -1][gi - 1 -1] - yVal);
break;
case 4: // South
wxx = xVal + xyf*(xArray[gj - 1 -1][gi -1] - xVal);
wyy = yVal + xyf*(yArray[gj - 1 -1][gi -1] - yVal);
break;
}
} else {
// We have 1D arrays for the X & Y grid points.
xVal = xArr1D[gi -1];
yVal = yArr1D[gj -1];
switch (idir) {
case 1: // East
wxx = xVal + xyf*(xArr1D[gi + 1 -1] - xVal);
wyy = yVal;
break;
case 2: // North
wxx = xVal;
wyy = yVal + xyf*(yArr1D[gj + 1 -1] - yVal);
break;
case 3: // West
wxx = xVal + xyf*(xArr1D[gi - 1 -1] - xVal);
wyy = yVal;
break;
case 4: // South
wxx = xVal;
wyy = yVal + xyf*(yArr1D[gj - 1 -1] - yVal);
break;
}
}
if (DEBUG) {
System.out.println("i, j = " + gi + "," + gj);
System.out.println("cont = " + (float)cont + ", fij = " + (float)fij +
", func = " + (float)func + ", xyf = " + (float)xyf);
System.out.println("xVal = " + (float)xVal + ", yVal = " + (float)yVal);
System.out.println("wxx = " + (float)wxx + ", wyy = " + (float)wyy);
}
// Figure out what to do with this point.
if (lnstrt == 1) {
// This is the 1st point in the contour line.
np = 1;
pathbufxt[np -1] = wxx;
pathbufyt[np -1] = wyy;
// Save starting point as wx, wy.
wx = wxx;
wy = wyy;
// Clear the first point flag, we've got one now.
lnstrt = 0;
} else {
boolean skipFlg = false;
// Second point and after comes here.
// Add a point to this line. Check for duplicate point first.
if (ignext == 2) {
if (wxx == pathbufxt[np -1] && wyy == pathbufyt[np -1]) {
ignext = 0;
skipFlg = true;
} else
ignext = 1;
}
if (!skipFlg) {
// Increment # of points in contour.
++np;
pathbufxt[np -1] = wxx;
pathbufyt[np -1] = wyy;
// See if the temporary array xt, yt are full.
if (np == kBufSize) {
accumContour(np, icont, pathbufxt, pathbufyt, cAttr[icont]);
// Last point becomes 1st point to continue.
pathbufxt[0] = pathbufxt[np -1];
pathbufyt[0] = pathbufyt[np -1];
np =1;
}
// Check to see if we're back to the intial point.
if (wxx == wx && wyy == wy)
return;
}
}
// Search for the next point on this line.
/*L67*/ switch(idir) {
case 1:
++gi;
if (!Routine_L51())
Routine_L21();
break;
case 2:
++gj;
Routine_L21();
break;
case 3:
--gi;
if (!Routine_L31())
Routine_L21();
break;
case 4:
--gj;
if (!Routine_L41())
Routine_L21();
break;
}
return;
}
/**
* Accumulate contour paths, as they are generated, into
* an overall list of contours.
*
* @param np The number of points in the contour path buffers.
* @param icont The index to the current contour level.
* @param x,y Buffers containing x & y coordinates of contour points.
* @param cAttr The attributes for this particular contour level.
**/
private void accumContour(int np, int icont, double[] x, double[] y, ContourAttrib cAttr) {
// To few points for a contour line.
if (np < kMinNumPoints) return;
// Copy over coordinate points from buffers to their own arrays.
double[] xArr = new double[np];
double[] yArr = new double[np];
System.arraycopy(x, 0, xArr, 0, np);
System.arraycopy(y, 0, yArr, 0, np);
// Create a new contour path and add it to the list.
ContourPath path = new ContourPath(cAttr, icont, xArr, yArr);
pathList.add(path);
}
}