Java tutorial
/* Copyright 2007-2009 361DEGRES This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program 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 General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ package ch.epfl.scapetoad; import com.vividsolutions.jts.geom.Coordinate; import com.vividsolutions.jts.geom.Envelope; import com.vividsolutions.jts.geom.GeometryFactory; import com.vividsolutions.jts.geom.LinearRing; import com.vividsolutions.jts.geom.Polygon; import com.vividsolutions.jump.feature.AttributeType; import com.vividsolutions.jump.feature.BasicFeature; import com.vividsolutions.jump.feature.FeatureDataset; import com.vividsolutions.jump.feature.FeatureSchema; import org.apache.commons.math.special.Erf; import org.apache.commons.math.MathException; /** * This class implements Gastner's diffusion algorithm. * @author Christian.Kaiser@91nord.com * @version v1.0.0, 2007-11-30 */ public class CartogramGastner { /** * The cartogram grid which contains the density and which we * will deform. */ CartogramGrid mGrid; /** * The cartogram bounding box. It is slightly larger than the * grid bounding box for computational reasons. */ Envelope mExtent; /** * Length of map in x direction. This is our grid size (number of cells). * The number should always be a power of 2. */ private int lx; /** * Length of map in y direction. This is our grid size (number of cells). * The number should always be a power of 2. */ private int ly; /** * Array for the density at time t = 0. */ private double[][] rho_0; /** * Array for the density at time t > 0. */ private double[][] rho; /** * Array for the velocity field in x direction at position (j,k). */ private double[][] gridvx; /** * Array for the velocity field in y direction at position (j,k). */ private double[][] gridvy; /** * Array for the x position at t > 0. x[j][k] is the x-coordinate for the * element that was at position (j,k) at time t = 0. */ private double[][] x; /** * Array for the y position at t > 0. y[j][k] is the y-coordinate for the * element that was at position (j,k) at time t = 0. */ private double[][] y; // Arrays for the velocity field at position (x[j][k], y[j][k]). private double[][] vx; private double[][] vy; // Definition of some other class wide variables. private double minpop = 0.0; private int nblurs = 0; private double[][] xappr; private double[][] yappr; /** * Some constants needed for the cartogram computation. */ private double CONVERGENCE = 1e-100; private double INFTY = 1e100; private double HINITIAL = 1e-4; private int IMAX = 50; private double MINH = 1e-5; private int MAXINTSTEPS = 3000; private double SIGMA = 0.1; private double SIGMAFAC = 1.2; private double TIMELIMIT = 1e8; private double TOLF = 1e-3; private double TOLINT = 1e-3; private double TOLX = 1e-3; /** * Attributes for the computation progress bar. */ public int mProgressStart; public int mProgressEnd; public String mProgressText; public CartogramWizard mCartogramWizard; /** * The constructor takes the cartogram grid which contains the * density values. */ public CartogramGastner(CartogramGrid grid) { mGrid = grid; } // CartogramGastner.<init> /** * Starts the cartogram computation using the given grid size. * @param gridSize the size of the grid used for computation. The * grid size must be a power of 2. */ public void compute(int gridSize) throws InterruptedException { // Store the grid size in the lx and ly attributes. lx = gridSize; ly = gridSize; this.initializeArrays(); this.computeInitialDensity(); FFT.coscosft(rho_0, 1, 1); boolean hasConverged = false; while (hasConverged == false) { if (Thread.interrupted()) { // Raise an InterruptedException. throw new InterruptedException("Computation has been interrupted by the user."); } hasConverged = this.integrateNonlinearVolterraEquation(); } this.projectCartogramGrid(); } // CartogramGastner.compute /** * Initializes the arrays using the grid size. */ private void initializeArrays() { rho_0 = new double[lx + 1][ly + 1]; rho = new double[lx + 1][ly + 1]; gridvx = new double[lx + 1][ly + 1]; gridvy = new double[lx + 1][ly + 1]; x = new double[lx + 1][ly + 1]; y = new double[lx + 1][ly + 1]; vx = new double[lx + 1][ly + 1]; vy = new double[lx + 1][ly + 1]; xappr = new double[lx + 1][ly + 1]; yappr = new double[lx + 1][ly + 1]; } // CartogramGastner.initializeArrays /** * Computes the initial density of the Gastner grid using the provided * cartogram grid. */ private void computeInitialDensity() { // Compute the cartogram extent. mExtent = this.cartogramExtent(mGrid.envelope(), lx, ly); // Compute the cell size in x and y direction. double cellSizeX = mExtent.getWidth() / lx; double cellSizeY = mExtent.getHeight() / ly; // Store the extent's minimum and maximum coordinates. double extentMinX = mExtent.getMinX() - (cellSizeX / 2); double extentMinY = mExtent.getMinY() - (cellSizeY / 2); double extentMaxX = mExtent.getMaxX() + (cellSizeX / 2); double extentMaxY = mExtent.getMaxY() + (cellSizeY / 2); // Let the cartogram grid fill in the density. mGrid.fillRegularDensityGrid(rho_0, extentMinX, extentMaxX, extentMinY, extentMaxY); // If there are 0 density values, introduce a small bias. double minimumDensity = rho_0[0][0]; double maximumDensity = rho_0[0][0]; for (int j = 0; j <= ly; j++) { for (int i = 0; i <= lx; i++) { if (rho_0[i][j] < minimumDensity) minimumDensity = rho_0[i][j]; if (rho_0[i][j] > maximumDensity) maximumDensity = rho_0[i][j]; } } if ((minimumDensity * 1000) < maximumDensity) { double bias = (maximumDensity / 1000) - minimumDensity; for (int j = 0; j <= ly; j++) for (int i = 0; i <= lx; i++) rho_0[i][j] += bias; } this.correctInitialDensityEdges(); } // CartogramGastner.computeInitialDensity /** * Fills the edges of the initial density grid correctly. */ private void correctInitialDensityEdges() { int i, j; rho_0[0][0] += rho_0[0][ly] + rho_0[lx][0] + rho_0[lx][ly]; for (i = 1; i < lx; i++) rho_0[i][0] += rho_0[i][ly]; for (j = 1; j < ly; j++) rho_0[0][j] += rho_0[lx][j]; for (i = 0; i < lx; i++) rho_0[i][ly] = rho_0[i][0]; for (j = 0; j <= ly; j++) rho_0[lx][j] = rho_0[0][j]; } // CartogramGastner.correctInitialDensityEdges /** * Computes the cartogram extent bounding box using the layer extent * and the grid size. * @param env the bounding box of the layers (or the cartogram grid). * @param gridX the grid size in x * @param gridY the grid size in y * @return the cartogram extent as Envelope */ private Envelope cartogramExtent(Envelope env, int gridX, int gridY) { double margin = 1.5; double minx, maxx, miny, maxy; if ((env.getWidth() / gridX) > (env.getHeight() / gridY)) { maxx = 0.5 * ((1 + margin) * env.getMaxX() + (1 - margin) * env.getMinX()); minx = 0.5 * ((1 - margin) * env.getMaxX() + (1 + margin) * env.getMinX()); maxy = 0.5 * (env.getMaxY() + env.getMinY() + (maxx - minx) * gridY / gridX); miny = 0.5 * (env.getMaxY() + env.getMinY() - (maxx - minx) * gridY / gridX); } else { maxy = 0.5 * ((1 + margin) * env.getMaxY() + (1 - margin) * env.getMinY()); miny = 0.5 * ((1 - margin) * env.getMaxY() + (1 + margin) * env.getMinY()); maxx = 0.5 * (env.getMaxX() + env.getMinX() + (maxy - miny) * gridX / gridY); minx = 0.5 * (env.getMaxX() + env.getMinX() - (maxy - miny) * gridX / gridY); } // Creating the Envelope Envelope bbox = new Envelope(minx, maxx, miny, maxy); return bbox; } // CartogramGastner.cartogramExtent /** * Integrates the non-linear Volterra equation. * @return true if the displacement field has converged, false otherwise. */ private boolean integrateNonlinearVolterraEquation() throws InterruptedException { boolean stepsize_ok; double h, maxchange = this.INFTY, t, vxplus, vyplus, xguess, yguess; int i, j, k; do { this.initcond(); this.nblurs++; //if (this.minpop < 0.0) // double sigmaVal = SIGMA * Math.pow(this.SIGMAFAC, this.nblurs); } while (this.minpop < 0.0); h = HINITIAL; t = 0; for (j = 0; j <= this.lx; j++) { for (k = 0; k <= this.ly; k++) { this.x[j][k] = j; this.y[j][k] = k; } } this.calculateVelocityField(0.0); for (j = 0; j <= this.lx; j++) { for (k = 0; k <= this.ly; k++) { vx[j][k] = gridvx[j][k]; vy[j][k] = gridvy[j][k]; } } i = 1; do { // Stop if the user has interrupted the process. if (Thread.interrupted()) { // Raise an InterruptedException. throw new InterruptedException("Computation has been interrupted by the user."); } stepsize_ok = true; this.calculateVelocityField(t + h); for (j = 0; j <= this.lx; j++) { for (k = 0; k <= this.ly; k++) { double xinterpol = this.x[j][k] + (h * this.vx[j][k]); double yinterpol = this.y[j][k] + (h * this.vy[j][k]); if (xinterpol < 0.0 || yinterpol < 0.0) { if (AppContext.DEBUG) System.out.println("[ERROR] Cartogram out of bounds !"); } vxplus = this.interpolateBilinear(this.gridvx, xinterpol, yinterpol); vyplus = this.interpolateBilinear(this.gridvy, xinterpol, yinterpol); xguess = this.x[j][k] + (0.5 * h * (this.vx[j][k] + vxplus)); yguess = this.y[j][k] + (0.5 * h * (this.vy[j][k] + vyplus)); double[] ptappr = new double[2]; ptappr[0] = this.xappr[j][k]; ptappr[1] = this.yappr[j][k]; boolean solving_ok = this.newt2(h, ptappr, xguess, yguess, j, k); this.xappr[j][k] = ptappr[0]; this.yappr[j][k] = ptappr[1]; if (solving_ok == false) return false; if (((xguess - this.xappr[j][k]) * (xguess - this.xappr[j][k])) + ((yguess - this.yappr[j][k]) * (yguess - this.yappr[j][k])) > this.TOLINT) { if (h < this.MINH) { //double sigmaVal = this.SIGMA * Math.pow( // this.SIGMAFAC, this.nblurs); this.nblurs++; return false; } h = h / 10; stepsize_ok = false; break; } } // for (k = 0; k <= this.ly; k++) } // for (j = 0; j <= this.lx; j++) if (!stepsize_ok) { continue; } else { t += h; maxchange = 0.0; for (j = 0; j <= this.lx; j++) { for (k = 0; k <= this.ly; k++) { if (((this.x[j][k] - this.xappr[j][k]) * (this.x[j][k] - this.xappr[j][k])) + ((this.y[j][k] - this.yappr[j][k]) * (this.y[j][k] - this.yappr[j][k])) > maxchange) { maxchange = ((this.x[j][k] - this.xappr[j][k]) * (this.x[j][k] - this.xappr[j][k])) + ((this.y[j][k] - this.yappr[j][k]) * (this.y[j][k] - this.yappr[j][k])); } this.x[j][k] = this.xappr[j][k]; this.y[j][k] = this.yappr[j][k]; this.vx[j][k] = this.interpolateBilinear(this.gridvx, this.xappr[j][k], this.yappr[j][k]); this.vy[j][k] = this.interpolateBilinear(this.gridvy, this.xappr[j][k], this.yappr[j][k]); } // for (k=0; k<=ly; k++) } // for (j = 0; j <= this.lx; j++) } h = 1.2 * h; int progress = mProgressEnd; if (i < 200) progress = mProgressStart + (i * ((mProgressEnd - mProgressStart) / 200)); mCartogramWizard.updateRunningStatus(progress, mProgressText, "Doing time step " + i); i++; } while (i < this.MAXINTSTEPS && t < this.TIMELIMIT && maxchange > this.CONVERGENCE); return true; } // CartogramGastner.integrateNonlinearVolterraEquation private void initcond() { double maxpop; int i, j; FFT.coscosft(rho_0, -1, -1); for (i = 0; i < this.lx; i++) { for (j = 0; j < this.ly; j++) { if (this.rho_0[i][j] < -1e10) { this.rho_0[i][j] = 0.0; } } } this.gaussianBlur(); minpop = this.rho_0[0][0]; maxpop = this.rho_0[0][0]; for (i = 0; i < this.lx; i++) { for (j = 0; j < this.ly; j++) { if (this.rho_0[i][j] < minpop) { minpop = rho_0[i][j]; } } } for (i = 0; i < this.lx; i++) { for (j = 0; j < this.ly; j++) { if (this.rho_0[i][j] > maxpop) { maxpop = this.rho_0[i][j]; } } } FFT.coscosft(this.rho_0, 1, 1); } // CartogramGastner.initcond /** * Performs au Gaussian blur on the density grid. */ public void gaussianBlur() { double[][][] blur = new double[1][this.lx][this.ly]; double[][][] conv = new double[1][this.lx][this.ly]; double[][][] pop = new double[1][this.lx][this.ly]; double[][] speqblur = new double[1][2 * this.lx]; double[][] speqconv = new double[1][2 * this.lx]; double[][] speqpop = new double[1][2 * this.lx]; int i, j, p, q; for (i = 1; i <= this.lx; i++) { for (j = 1; j <= this.ly; j++) { if (i > (this.lx / 2)) p = i - 1 - this.lx; else p = i - 1; if (j > (this.ly / 2)) q = j - 1 - this.ly; else q = j - 1; pop[0][i - 1][j - 1] = this.rho_0[i - 1][j - 1]; double erfDenominator = Math.sqrt(2.0) * (this.SIGMA * Math.pow(this.SIGMAFAC, this.nblurs)); double erfParam1 = this.erf((p + 0.5) / erfDenominator); double erfParam2 = this.erf((p - 0.5) / erfDenominator); double erfParam3 = this.erf((q + 0.5) / erfDenominator); double erfParam4 = this.erf((q - 0.5) / erfDenominator); conv[0][i - 1][j - 1] = 0.5 * (erfParam1 - erfParam2) * (erfParam3 - erfParam4) / (this.lx * this.ly); } } FFT.rlft3(pop, speqpop, 1, this.lx, this.ly, 1); FFT.rlft3(conv, speqconv, 1, this.lx, this.ly, 1); for (i = 1; i <= this.lx; i++) { for (j = 1; j <= (this.ly / 2); j++) { blur[0][i - 1][(2 * j) - 2] = pop[0][i - 1][(2 * j) - 2] * conv[0][i - 1][(2 * j) - 2] - pop[0][i - 1][(2 * j) - 1] * conv[0][i - 1][(2 * j) - 1]; blur[0][i - 1][(2 * j) - 1] = pop[0][i - 1][(2 * j) - 1] * conv[0][i - 1][(2 * j) - 2] + pop[0][i - 1][(2 * j) - 2] * conv[0][i - 1][(2 * j) - 1]; } } for (i = 1; i <= this.lx; i++) { speqblur[0][(2 * i) - 2] = speqpop[0][(2 * i) - 2] * speqconv[0][(2 * i) - 2] - speqpop[0][(2 * i) - 1] * speqconv[0][(2 * i) - 1]; speqblur[0][(2 * i) - 1] = speqpop[0][(2 * i) - 1] * speqconv[0][(2 * i) - 2] + speqpop[0][(2 * i) - 2] * speqconv[0][(2 * i) - 1]; } FFT.rlft3(blur, speqblur, 1, this.lx, this.ly, -1); for (i = 1; i <= this.lx; i++) { for (j = 1; j <= this.ly; j++) { rho_0[i - 1][j - 1] = blur[0][i - 1][j - 1]; } } } // CartogramGastner.gaussianBlur /** * Computes the velocity field at time t. * @param t the time. */ private void calculateVelocityField(double t) { int j, k; for (j = 0; j <= this.lx; j++) { for (k = 0; k <= this.ly; k++) { this.rho[j][k] = Math.exp(-1 * ((Math.PI * j / this.lx) * (Math.PI * j / this.lx) + (Math.PI * k / this.ly) * (Math.PI * k / this.ly)) * t) * this.rho_0[j][k]; } } for (j = 0; j <= this.lx; j++) { for (k = 0; k <= this.ly; k++) { this.gridvx[j][k] = -1 * (Math.PI * j / this.lx) * this.rho[j][k]; this.gridvy[j][k] = -1 * (Math.PI * k / this.ly) * this.rho[j][k]; } } FFT.coscosft(this.rho, -1, -1); FFT.sincosft(this.gridvx, -1, -1); FFT.cossinft(this.gridvy, -1, -1); for (j = 0; j <= this.lx; j++) { for (k = 0; k <= this.ly; k++) { this.gridvx[j][k] = (-1 * this.gridvx[j][k]) / this.rho[j][k]; this.gridvy[j][k] = (-1 * this.gridvy[j][k]) / this.rho[j][k]; } } } // CartogramGastner.calculateVelocityField /** * Bilinear interpolation in 2D. */ private double interpolateBilinear(double[][] arr, double x, double y) { int gaussx, gaussy; double deltax, deltay; if (x < 0.0 || y < 0.0) { return 0.0; } int xlen = arr.length; int ylen = arr[0].length; if (x >= xlen || y >= ylen) { return 0.0; } Double xobj = new Double(x); gaussx = xobj.intValue(); Double yobj = new Double(y); gaussy = yobj.intValue(); deltax = x - gaussx; deltay = y - gaussy; if (gaussx == this.lx && gaussy == this.ly) return arr[gaussx][gaussy]; if (gaussx == this.lx) return ((1 - deltay) * arr[gaussx][gaussy]) + (deltay * arr[gaussx][gaussy + 1]); if (gaussy == this.ly) return ((1 - deltax) * arr[gaussx][gaussy]) + (deltax * arr[gaussx + 1][gaussy]); return ((1 - deltax) * (1 - deltay) * arr[gaussx][gaussy]) + ((1 - deltax) * deltay * arr[gaussx][gaussy + 1]) + (deltax * (1 - deltay) * arr[gaussx + 1][gaussy]) + (deltax * deltay * arr[gaussx + 1][gaussy + 1]); } // CartogramGastner.interpolateBilinear public boolean newt2(double h, double[] ptappr, double xguess, double yguess, int j, int k) { double deltax, deltay, dfxdx, dfxdy, dfydx, dfydy, fx, fy; int gaussx, gaussxplus, gaussy, gaussyplus, i; double temp; Double tempobj = null; ptappr[0] = xguess; ptappr[1] = yguess; for (i = 1; i <= IMAX; i++) { temp = this.interpolateBilinear(this.gridvx, ptappr[0], ptappr[1]); fx = ptappr[0] - (0.5 * h * temp) - this.x[j][k] - (0.5 * h * this.vx[j][k]); temp = this.interpolateBilinear(this.gridvy, ptappr[0], ptappr[1]); fy = ptappr[1] - (0.5 * h * temp) - this.y[j][k] - (0.5 * h * this.vy[j][k]); tempobj = new Double(ptappr[0]); gaussx = tempobj.intValue(); tempobj = new Double(ptappr[1]); gaussy = tempobj.intValue(); if (gaussx == this.lx) gaussxplus = 0; else gaussxplus = gaussx + 1; if (gaussy == this.ly) gaussyplus = 0; else gaussyplus = gaussy + 1; deltax = this.x[j][k] - gaussx; deltay = this.y[j][k] - gaussy; dfxdx = 1 - 0.5 * h * ((1 - deltay) * (this.gridvx[gaussxplus][gaussy] - this.gridvx[gaussx][gaussy]) + deltay * (this.gridvx[gaussxplus][gaussyplus] - this.gridvx[gaussx][gaussyplus])); dfxdy = -0.5 * h * ((1 - deltax) * (this.gridvx[gaussx][gaussyplus] - this.gridvx[gaussx][gaussy]) + deltax * (this.gridvx[gaussxplus][gaussyplus] - this.gridvx[gaussxplus][gaussy])); dfydx = -0.5 * h * ((1 - deltay) * (this.gridvy[gaussxplus][gaussy] - this.gridvy[gaussx][gaussy]) + deltay * (this.gridvy[gaussxplus][gaussyplus] - this.gridvy[gaussx][gaussyplus])); dfydy = 1 - 0.5 * h * ((1 - deltax) * (this.gridvy[gaussx][gaussyplus] - this.gridvy[gaussx][gaussy]) + deltax * (this.gridvy[gaussxplus][gaussyplus] - this.gridvy[gaussxplus][gaussy])); if ((fx * fx + fy * fy) < this.TOLF) return true; deltax = (fy * dfxdy - fx * dfydy) / (dfxdx * dfydy - dfxdy * dfydx); deltay = (fx * dfydx - fy * dfxdx) / (dfxdx * dfydy - dfxdy * dfydx); if ((deltax * deltax + deltay * deltay) < this.TOLX) return true; ptappr[0] += deltax; ptappr[1] += deltay; } return false; } // CartogramGastner.newt2 /** * Our wrapper function for the Jakarta Commons Math Erf.erf function. * For values <= -4 or >= 4, we return -1 or 1 directly, without * computation. Erf.erf raises too often an exception for failing * convergence. */ public static double erf(double value) { if (value <= -4.0) return -1.0; if (value >= 4.0) return 1.0; double erf = 0.0; try { erf = Erf.erf(value); } catch (MathException mexc) { if (value < 0) return -1.0; else return 1.0; } return erf; } // CartogramGastner.erf /** * Applies the cartogram deformation to the cartogram grid. */ private void projectCartogramGrid() { // Project each point in the cartogram grid. double[][] x = mGrid.getXCoordinates(); double[][] y = mGrid.getYCoordinates(); int gridSizeX = x.length; int gridSizeY = x[0].length; int i, j; for (i = 0; i < gridSizeX; i++) { for (j = 0; j < gridSizeY; j++) { double[] projectedPoint = this.projectPoint(x[i][j], y[i][j]); x[i][j] = projectedPoint[0]; y[i][j] = projectedPoint[1]; } } } // CartogramGastner.projectCartogramGrid /** * Projects one point using the deformed grid. * @param x the x coordinate of the point to project. * @param y the y coordinate of the point to project. * @return a double array with the coordinates of the projected point. */ private double[] projectPoint(double x, double y) { double deltax, deltay, den, t, u, temp; long gaussx, gaussy; // Get the grid size and the cellsize. double cellSizeX = mExtent.getWidth() / lx; double cellSizeY = mExtent.getHeight() / ly; // Make a copy of the point coordinate. double px = x; double py = y; px = (px - mExtent.getMinX()) * lx / mExtent.getWidth(); py = (py - mExtent.getMinY()) * ly / mExtent.getHeight(); temp = Math.floor(px); gaussx = Math.round(temp); temp = Math.floor(py); gaussy = Math.round(temp); if (gaussx < 0 || gaussx > lx || gaussy < 0 || gaussy > ly) { System.out.println("[ERROR] Coordinate limits exceeded."); return null; } deltax = px - gaussx; deltay = py - gaussy; double ax = (1 - deltax) * this.x[(int) gaussx][(int) gaussy] + deltax * this.x[(int) (gaussx + 1)][(int) gaussy]; double ay = (1 - deltax) * this.y[(int) gaussx][(int) gaussy] + deltax * this.y[(int) (gaussx + 1)][(int) gaussy]; double bx = (1 - deltax) * this.x[(int) gaussx][(int) (gaussy + 1)] + deltax * this.x[(int) (gaussx + 1)][(int) (gaussy + 1)]; double by = (1 - deltax) * this.y[(int) gaussx][(int) (gaussy + 1)] + deltax * this.y[(int) (gaussx + 1)][(int) (gaussy + 1)]; double cx = (1 - deltay) * this.x[(int) gaussx][(int) gaussy] + deltay * this.x[(int) gaussx][(int) (gaussy + 1)]; double cy = (1 - deltay) * this.y[(int) gaussx][(int) gaussy] + deltay * this.y[(int) gaussx][(int) (gaussy + 1)]; double dx = (1 - deltay) * this.x[(int) (gaussx + 1)][(int) gaussy] + deltay * this.x[(int) (gaussx + 1)][(int) (gaussy + 1)]; double dy = (1 - deltay) * this.y[(int) (gaussx + 1)][(int) gaussy] + deltay * this.y[(int) (gaussx + 1)][(int) (gaussy + 1)]; den = (bx - ax) * (cy - dy) + (ay - by) * (cx - dx); if (Math.abs(den) < 1e-12) { double ix = (ax + bx + cx + dx) / 4; double iy = (ay + by + cy + dy) / 4; double meanpoint[] = new double[2]; meanpoint[0] = (ix * (lx / mExtent.getWidth())) + mExtent.getMinX(); meanpoint[1] = (iy * (ly / mExtent.getHeight())) + mExtent.getMinY(); return meanpoint; } t = ((cx - ax) * (cy - dy) + (ay - cy) * (cx - dx)) / den; u = ((bx - ax) * (cy - ay) + (ay - by) * (cx - ax)) / den; px = (1 - (ax + t * (bx - ax)) / lx) * mExtent.getMinX() + ((ax + t * (bx - ax)) / lx) * mExtent.getMaxX(); py = (1 - (ay + t * (by - ay)) / ly) * mExtent.getMinY() + ((ay + t * (by - ay)) / ly) * mExtent.getMaxY(); double point[] = new double[2]; point[0] = px; point[1] = py; return point; } // CartogramGastner.projectPointWithGrid /** * Writes the grid into the specified shape file. * @param shapefile the path to the shape file. */ public void writeToShapefile(String shapefile) { // Create a new Feature Schema for our shape file. FeatureSchema fs = new FeatureSchema(); // We add the following attributes to the Feature Schema: // cellId : a serial number starting at 1 // geom : the geometry (polygon) // i : the index of the cell in x direction // j : the index of the cell in y direction fs.addAttribute("cellId", AttributeType.INTEGER); fs.addAttribute("geom", AttributeType.GEOMETRY); fs.addAttribute("i", AttributeType.INTEGER); fs.addAttribute("j", AttributeType.INTEGER); fs.addAttribute("rho_0", AttributeType.DOUBLE); fs.addAttribute("rho", AttributeType.DOUBLE); // Create a new Geometry Factory for creating our geometries. GeometryFactory gf = new GeometryFactory(); // Create a new Feature Dataset in order to store our new Features. FeatureDataset fd = new FeatureDataset(fs); // Create one Feature for each cell. int i, j; int cellId = 0; for (j = 0; j < ly; j++) { for (i = 0; i < lx; i++) { cellId++; // Extract the coordinates for the cell polygon. Coordinate[] coords = new Coordinate[5]; coords[0] = new Coordinate(x[i][j], y[i][j]); coords[1] = new Coordinate(x[i][j + 1], y[i][j + 1]); coords[2] = new Coordinate(x[i + 1][j + 1], y[i + 1][j + 1]); coords[3] = new Coordinate(x[i + 1][j], y[i + 1][j]); coords[4] = coords[0]; // Create the polygon. LinearRing ring = gf.createLinearRing(coords); Polygon poly = gf.createPolygon(ring, null); // Create a new Feature. BasicFeature feat = new BasicFeature(fs); // Setting the Feature's attributes. feat.setAttribute("cellId", new Integer(cellId)); feat.setAttribute("geom", poly); feat.setAttribute("i", new Integer(i)); feat.setAttribute("j", new Integer(j)); feat.setAttribute("rho_0", new Double(rho_0[i][j])); feat.setAttribute("rho", new Double(rho[i][j])); // Add the Feature to the Feature Dataset. fd.add(feat); } } // Write the Feature Dataset to the Shape file. IOManager.writeShapefile(fd, shapefile); } // CartogramGrid.writeToShapefile } // CartogramGastner class FFT { public static void coscosft(double[][] y, int isign1, int isign2) { int lx = y.length - 1; int ly = y[0].length - 1; double temp[] = new double[lx + 1]; int i, j; for (i = 0; i <= lx; i++) { FFT.cosft(y[i], ly, isign2); } for (j = 0; j <= ly; j++) { for (i = 0; i <= lx; i++) { temp[i] = y[i][j]; } FFT.cosft(temp, lx, isign1); for (i = 0; i <= lx; i++) { y[i][j] = temp[i]; } } } public static void cosft(double z[], int n, int isign) { double theta, wi = 0.0, wpi, wpr, wr = 1.0, wtemp; double[] a; double sum, y1, y2; int j, n2; a = new double[n + 2]; for (j = 1; j <= (n + 1); j++) { a[j] = z[j - 1]; } theta = Math.PI / n; wtemp = Math.sin(0.5 * theta); wpr = -2.0 * wtemp * wtemp; wpi = Math.sin(theta); sum = 0.5 * (a[1] - a[n + 1]); a[1] = 0.5 * (a[1] + a[n + 1]); n2 = n + 2; for (j = 2; j <= (n / 2); j++) { wtemp = wr; wr = wr * wpr - wi * wpi + wr; wi = wi * wpr + wtemp * wpi + wi; y1 = 0.5 * (a[j] + a[n2 - j]); y2 = a[j] - a[n2 - j]; a[j] = y1 - wi * y2; a[n2 - j] = y1 + wi * y2; sum += wr * y2; } FFT.realft(a, n, 1); a[n + 1] = a[2]; a[2] = sum; for (j = 4; j <= n; j += 2) { sum += a[j]; a[j] = sum; } if (isign == 1) { for (j = 1; j <= (n + 1); j++) z[j - 1] = a[j]; } else if (isign == -1) { for (j = 1; j <= (n + 1); j++) z[j - 1] = 2.0 * a[j] / n; } } // FFT.cosft public static void cossinft(double[][] y, int isign1, int isign2) { int lx = y.length - 1; int ly = y[0].length - 1; double[] temp = new double[lx + 1]; int i, j; for (i = 0; i <= lx; i++) FFT.sinft(y[i], ly, isign2); for (j = 0; j <= ly; j++) { for (i = 0; i <= lx; i++) temp[i] = y[i][j]; FFT.cosft(temp, lx, isign1); for (i = 0; i <= lx; i++) y[i][j] = temp[i]; } } public static void four1(double[] data, int nn, int isign) { double theta, wi, wpi, wpr, wr, wtemp; double tempi, tempr; int i, istep, j, m, mmax, n; n = nn * 2; j = 1; for (i = 1; i < n; i += 2) { if (j > i) { tempr = data[j]; data[j] = data[i]; data[i] = tempr; tempr = data[j + 1]; data[j + 1] = data[i + 1]; data[i + 1] = tempr; } m = n / 2; while (m >= 2 && j > m) { j -= m; m = m / 2; } j += m; } mmax = 2; while (n > mmax) { istep = mmax * 2; theta = isign * (2 * Math.PI / mmax); wtemp = Math.sin(0.5 * theta); wpr = -2.0 * wtemp * wtemp; wpi = Math.sin(theta); wr = 1.0; wi = 0.0; for (m = 1; m < mmax; m += 2) { for (i = m; i <= n; i += istep) { j = i + mmax; tempr = wr * data[j] - wi * data[j + 1]; tempi = wr * data[j + 1] + wi * data[j]; data[j] = data[i] - tempr; data[j + 1] = data[i + 1] - tempi; data[i] += tempr; data[i + 1] += tempi; } wtemp = wr; wr = wtemp * wpr - wi * wpi + wr; wi = wi * wpr + wtemp * wpi + wi; } mmax = istep; } } public static void fourn(double[] data, int[] nn, int ndim, int isign) { int idim; int i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2; int ibit, k1, k2, n, nprev, nrem, ntot; double tempi, tempr; double theta, wi, wpi, wpr, wr, wtemp; double swaptemp; ntot = 1; for (idim = 1; idim <= ndim; idim++) { ntot *= nn[idim]; } nprev = 1; for (idim = ndim; idim >= 1; idim--) { n = nn[idim]; nrem = ntot / (n * nprev); ip1 = nprev * 2; ip2 = ip1 * n; ip3 = ip2 * nrem; i2rev = 1; for (i2 = 1; i2 <= ip2; i2 += ip1) { if (i2 < i2rev) { for (i1 = i2; i1 <= (i2 + ip1 - 2); i1 += 2) { for (i3 = i1; i3 <= ip3; i3 += ip2) { i3rev = i2rev + i3 - i2; swaptemp = data[i3 - 1]; data[i3 - 1] = data[i3rev - 1]; data[i3rev - 1] = swaptemp; swaptemp = data[i3]; data[i3] = data[i3rev]; data[i3rev] = swaptemp; } } } ibit = ip2 / 2; while (ibit >= ip1 && i2rev > ibit) { i2rev -= ibit; ibit = ibit / 2; } i2rev += ibit; } ifp1 = ip1; while (ifp1 < ip2) { ifp2 = ifp1 * 2; theta = 2 * isign * Math.PI / (ifp2 / ip1); wtemp = Math.sin(0.5 * theta); wpr = -2.0 * wtemp * wtemp; wpi = Math.sin(theta); wr = 1.0; wi = 0.0; for (i3 = 1; i3 <= ifp1; i3 += ip1) { for (i1 = i3; i1 <= (i3 + ip1 - 2); i1 += 2) { for (i2 = i1; i2 <= ip3; i2 += ifp2) { k1 = i2; k2 = k1 + ifp1; tempr = (wr * data[k2 - 1]) - (wi * data[k2]); tempi = (wr * data[k2]) + (wi * data[k2 - 1]); data[k2 - 1] = data[k1 - 1] - tempr; data[k2] = data[k1] - tempi; data[k1 - 1] += tempr; data[k1] += tempi; } } wtemp = wr; wr = wtemp * wpr - wi * wpi + wr; wi = wi * wpr + wtemp * wpi + wi; } ifp1 = ifp2; } nprev *= n; } } public static void realft(double[] data, int n, int isign) { double theta, wi, wpi, wpr, wr, wtemp; double c1 = 0.5, c2, h1i, h1r, h2i, h2r; int i, i1, i2, i3, i4, np3; theta = Math.PI / (n / 2); if (isign == 1) { c2 = -0.5; FFT.four1(data, (n / 2), 1); } else { c2 = 0.5; theta = -1.0 * theta; } wtemp = Math.sin(0.5 * theta); wpr = -2.0 * wtemp * wtemp; wpi = Math.sin(theta); wr = 1.0 + wpr; wi = wpi; np3 = n + 3; for (i = 2; i <= (n / 4); i++) { i1 = i + i - 1; i2 = 1 + i + i - 1; i3 = np3 - i2; i4 = 1 + i3; h1r = c1 * (data[i1] + data[i3]); h1i = c1 * (data[i2] - data[i4]); h2r = (-1.0 * c2) * (data[i2] + data[i4]); h2i = c2 * (data[i1] - data[i3]); data[i1] = h1r + wr * h2r - wi * h2i; data[i2] = h1i + wr * h2i + wi * h2r; data[i3] = h1r - wr * h2r + wi * h2i; data[i4] = (-1.0 * h1i) + wr * h2i + wi * h2r; wtemp = wr; wr = wr * wpr - wi * wpi + wr; wi = wi * wpr + wtemp * wpi + wi; } if (isign == 1) { h1r = data[1]; data[1] = h1r + data[2]; data[2] = h1r - data[2]; } else { h1r = data[1]; data[1] = c1 * (h1r + data[2]); data[2] = c1 * (h1r - data[2]); FFT.four1(data, (n / 2), -1); } } public static void rlft3(double[][][] data, double[][] speq, int nn1, int nn2, int nn3, int isign) { double theta, wi, wpi, wpr, wr, wtemp; double c1, c2, h1r, h1i, h2r, h2i; int i1, i2, i3, j1, j2, j3, ii3; int[] nn = new int[4]; c1 = 0.5; c2 = -0.5 * isign; theta = 2 * isign * (Math.PI / nn3); wtemp = Math.sin(0.5 * theta); wpr = -2.0 * wtemp * wtemp; wpi = Math.sin(theta); nn[1] = nn1; nn[2] = nn2; nn[3] = nn3 / 2; double[] datatemp = new double[(nn1 * nn2 * nn3)]; if (isign == 1) { j1 = 0; for (i1 = 0; i1 < nn1; i1++) { for (i2 = 0; i2 < nn2; i2++) { for (i3 = 0; i3 < nn3; i3++) { datatemp[j1] = data[i1][i2][i3]; j1++; } } } FFT.fourn(datatemp, nn, 3, isign); j1 = 0; for (i1 = 0; i1 < nn1; i1++) { for (i2 = 0; i2 < nn2; i2++) { for (i3 = 0; i3 < nn3; i3++) { data[i1][i2][i3] = datatemp[j1]; j1++; } } } for (i1 = 1; i1 <= nn1; i1++) { for (i2 = 1, j2 = 0; i2 <= nn2; i2++) { speq[i1 - 1][j2] = data[i1 - 1][i2 - 1][0]; j2++; speq[i1 - 1][j2] = data[i1 - 1][i2 - 1][1]; j2++; } } } for (i1 = 1; i1 <= nn1; i1++) { if (i1 != 1) j1 = nn1 - i1 + 2; else j1 = 1; wr = 1.0; wi = 0.0; for (ii3 = 1, i3 = 1; i3 <= ((nn3 / 4) + 1); i3++, ii3 += 2) { for (i2 = 1; i2 <= nn2; i2++) { if (i3 == 1) { if (i2 != 1) j2 = ((nn2 - i2) * 2) + 3; else j2 = 1; h1r = c1 * (data[i1 - 1][i2 - 1][0] + speq[j1 - 1][j2 - 1]); h1i = c1 * (data[i1 - 1][i2 - 1][1] - speq[j1 - 1][j2]); h2i = c2 * (data[i1 - 1][i2 - 1][0] - speq[j1 - 1][j2 - 1]); h2r = -1 * c2 * (data[i1 - 1][i2 - 1][1] + speq[j1 - 1][j2]); data[i1 - 1][i2 - 1][0] = h1r + h2r; data[i1 - 1][i2 - 1][1] = h1i + h2i; speq[j1 - 1][j2 - 1] = h1r - h2r; speq[j1 - 1][j2] = h2i - h1i; } else { if (i2 != 1) j2 = nn2 - i2 + 2; else j2 = 1; j3 = nn3 + 3 - (i3 * 2); h1r = c1 * (data[i1 - 1][i2 - 1][ii3 - 1] + data[j1 - 1][j2 - 1][j3 - 1]); h1i = c1 * (data[i1 - 1][i2 - 1][ii3] - data[j1 - 1][j2 - 1][j3]); h2i = c2 * (data[i1 - 1][i2 - 1][ii3 - 1] - data[j1 - 1][j2 - 1][j3 - 1]); h2r = -1 * c2 * (data[i1 - 1][i2 - 1][ii3] + data[j1 - 1][j2 - 1][j3]); data[i1 - 1][i2 - 1][ii3 - 1] = h1r + wr * h2r - wi * h2i; data[i1 - 1][i2 - 1][ii3] = h1i + wr * h2i + wi * h2r; data[j1 - 1][j2 - 1][j3 - 1] = h1r - wr * h2r + wi * h2i; data[j1 - 1][j2 - 1][j3] = -1 * h1i + wr * h2i + wi * h2r; } } wtemp = wr; wr = wtemp * wpr - wi * wpi + wr; wi = wi * wpr + wtemp * wpi + wi; } } if (isign == -1) { j1 = 0; for (i1 = 0; i1 < nn1; i1++) { for (i2 = 0; i2 < nn2; i2++) { for (i3 = 0; i3 < nn3; i3++) { datatemp[j1] = data[i1][i2][i3]; j1++; } } } FFT.fourn(datatemp, nn, 3, isign); j1 = 0; for (i1 = 0; i1 < nn1; i1++) { for (i2 = 0; i2 < nn2; i2++) { for (i3 = 0; i3 < nn3; i3++) { data[i1][i2][i3] = datatemp[j1]; j1++; } } } } } public static void sincosft(double[][] y, int isign1, int isign2) { int lx = y.length - 1; int ly = y[0].length - 1; double[] temp = new double[lx + 1]; int i, j; for (i = 0; i <= lx; i++) { FFT.cosft(y[i], ly, isign2); } for (j = 0; j <= ly; j++) { for (i = 0; i <= lx; i++) { temp[i] = y[i][j]; } FFT.sinft(temp, lx, isign1); for (i = 0; i <= lx; i++) { y[i][j] = temp[i]; } } } public static void sinft(double[] z, int n, int isign) { double theta, wi = 0.0, wpi, wpr, wr = 1.0, wtemp; double[] a; double sum, y1, y2; int j; int n2 = n + 2; a = new double[n + 1]; for (j = 1; j <= n; j++) { a[j] = z[j - 1]; } theta = Math.PI / n; wtemp = Math.sin(0.5 * theta); wpr = -2.0 * wtemp * wtemp; wpi = Math.sin(theta); a[1] = 0.0; for (j = 2; j <= ((n / 2) + 1); j++) { wtemp = wr; wr = wtemp * wpr - wi * wpi + wr; wi = wi * wpr + wtemp * wpi + wi; y1 = wi * (a[j] + a[n2 - j]); y2 = 0.5 * (a[j] - a[n2 - j]); a[j] = y1 + y2; a[n2 - j] = y1 - y2; } FFT.realft(a, n, 1); a[1] *= 0.5; sum = a[2] = 0.0; for (j = 1; j <= (n - 1); j += 2) { sum += a[j]; a[j] = a[j + 1]; a[j + 1] = sum; } if (isign == 1) { for (j = 1; j <= n; j++) { z[j - 1] = a[j]; } } else if (isign == -1) { for (j = 1; j <= n; j++) { z[j - 1] = 2.0 * a[j] / n; } } z[n] = 0.0; } }