[mkgmap-dev] [PATCH] contours v3

5 Jul 2009

This version allows to store the contours in separate IMG files using 
the --dem-separate-img option.
The name (8-digit number) of the IMG file will be the original mapname + 
10000000.

I also replaced orthodromicDistance by quickDistance (I only use the 
distance calculation to find out if I should refine a contour line 
segment by splitting it).

Best wishes
Christian

Index: src/uk/me/parabola/mkgmap/Option.java
===================================================================

--- src/uk/me/parabola/mkgmap/Option.java	(Revision 1078)
+++ src/uk/me/parabola/mkgmap/Option.java	(Arbeitskopie)
@@ -34,7 +34,7 @@
 		}
 	}
 
-	protected Option(String option, String value) {
+	public Option(String option, String value) {
 		this.option = option;
 		this.value = value;
 	}
Index: src/uk/me/parabola/mkgmap/CommandArgsReader.java
===================================================================
--- src/uk/me/parabola/mkgmap/CommandArgsReader.java	(Revision 1078)
+++ src/uk/me/parabola/mkgmap/CommandArgsReader.java	(Arbeitskopie)
@@ -129,7 +129,7 @@
 	 * @param option The option name.
 	 * @param value Its value.
 	 */
-	private void addOption(String option, String value) {
+	public void addOption(String option, String value) {
 		CommandOption opt = new CommandOption(option, value);
 		addOption(opt);
 	}
@@ -181,13 +181,21 @@
 		return arglist.iterator();
 	}
 
+        public ArgList getArgList() {
+	        return arglist;
+	}
+
+        public EnhancedProperties getArgs() {
+	    return args;
+	}
+
 	/**
 	 * Read a config file that contains more options.  When the number of
 	 * options becomes large it is more convenient to place them in a file.
 	 *
 	 * @param filename The filename to obtain options from.
 	 */
-	private void readConfigFile(String filename) {
+	public void readConfigFile(String filename) {
 		Options opts = new Options(new OptionProcessor() {
 			public void processOption(Option opt) {
 				log.debug("incoming opt", opt.getOption(), opt.getValue());
@@ -206,14 +214,14 @@
 	 * the argument to be processed in order.  Options can be intersperced with
 	 * filenames.  The options take effect where they appear.
 	 */
-	interface ArgType {
+	public interface ArgType {
 		public abstract void processArg();
 	}
 
 	/**
 	 * A filename.
 	 */
-	class Filename implements ArgType {
+	public class Filename implements ArgType {
 		private final String name;
 		private boolean useFilenameAsMapname = true;
 
@@ -270,7 +278,7 @@
 	/**
 	 * An option argument.  A key value pair.
 	 */
-	class CommandOption implements ArgType {
+	public class CommandOption implements ArgType {
 		private final Option option;
 
 		private CommandOption(Option option) {
@@ -297,7 +305,7 @@
 	/**
 	 * The arguments are held in this list.
 	 */
-	class ArgList implements Iterable<CommandArgsReader.ArgType> {
+	public class ArgList implements Iterable<CommandArgsReader.ArgType> {
 		private final List<ArgType> alist;
 
 		private int filenameCount;
Index: src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java
===================================================================
--- src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java	(Revision 1078)
+++ src/uk/me/parabola/mkgmap/reader/MapperBasedMapDataSource.java	(Arbeitskopie)
@@ -27,12 +27,14 @@
 import uk.me.parabola.imgfmt.app.Area;
 import uk.me.parabola.imgfmt.app.Coord;
 import uk.me.parabola.imgfmt.app.trergn.Overview;
+import uk.me.parabola.mkgmap.general.LoadableMapDataSource;
 import uk.me.parabola.mkgmap.general.MapDataSource;
 import uk.me.parabola.mkgmap.general.MapDetails;
 import uk.me.parabola.mkgmap.general.MapLine;
 import uk.me.parabola.mkgmap.general.MapPoint;
 import uk.me.parabola.mkgmap.general.MapShape;
 import uk.me.parabola.mkgmap.general.RoadNetwork;
+import uk.me.parabola.mkgmap.reader.dem.DEM;
 import uk.me.parabola.util.Configurable;
 import uk.me.parabola.util.EnhancedProperties;
 
@@ -118,6 +120,11 @@
 		return configProps;
 	}
 
+        public MapDetails getMapper()
+        {
+	        return mapper;
+	}
+
 	/**
 	 * We add the background polygons if the map is not transparent.
 	 */
@@ -146,5 +153,8 @@
 			// the overview section.
 			mapper.addShape(background);
 		}
+		if (getConfig().getProperty("contours", false)) {		    
+		    DEM.createContours((LoadableMapDataSource) this, getConfig());
+		}
 	}
 }
Index: src/uk/me/parabola/mkgmap/reader/osm/Way.java
===================================================================
--- src/uk/me/parabola/mkgmap/reader/osm/Way.java	(Revision 1078)
+++ src/uk/me/parabola/mkgmap/reader/osm/Way.java	(Arbeitskopie)
@@ -30,12 +30,18 @@
  */
 public class Way extends Element {
 
-	private final List<Coord> points = new ArrayList<Coord>();
+        private final List<Coord> points;
 
 	public Way(long id) {
+                points = new ArrayList<Coord>();
 		setId(id);
 	}
 
+        public Way(long id, List<Coord> points) {
+	        this.points = points;
+		setId(id);
+	}
+
 	/**
 	 * Get the points that make up the way.  We attempt to re-order the segments
 	 * and return a list of points that traces the route of the way.

/*
 * Copyright (C) 2009 Christian Gawron
 * 
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License version 2 as
 *  published by the Free Software Foundation.
 * 
 *  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.
 * 
 * 
 * Author: Christian Gawron
 * Create date: 03-Jul-2009
 */
package uk.me.parabola.mkgmap.reader.dem;

import java.io.*;
import java.nio.channels.FileChannel;
import java.nio.MappedByteBuffer;

import java.util.Vector;
import java.util.List;
import java.util.ArrayList;
import java.util.Locale;

import com.sun.media.jai.codec.*;
import javax.media.jai.*;
import java.awt.image.renderable.ParameterBlock;
import java.awt.image.Raster;
import java.awt.image.DataBuffer;

import uk.me.parabola.imgfmt.Utils;
import uk.me.parabola.imgfmt.ExitException;
import uk.me.parabola.imgfmt.FormatException;
import uk.me.parabola.imgfmt.FileSystemParam;
import uk.me.parabola.imgfmt.app.Area;
import uk.me.parabola.imgfmt.app.Coord;
import uk.me.parabola.imgfmt.app.map.Map;

import uk.me.parabola.mkgmap.build.MapBuilder;
import uk.me.parabola.mkgmap.general.MapDetails;
import uk.me.parabola.mkgmap.general.MapPoint;
import uk.me.parabola.mkgmap.general.MapLine;
import uk.me.parabola.mkgmap.general.LevelInfo;
import uk.me.parabola.mkgmap.general.LoadableMapDataSource;
import uk.me.parabola.mkgmap.reader.MapperBasedMapDataSource;
import uk.me.parabola.mkgmap.reader.osm.OsmConverter;
import uk.me.parabola.mkgmap.reader.osm.Style;
import uk.me.parabola.mkgmap.reader.osm.Way;
import uk.me.parabola.mkgmap.osmstyle.StyleImpl;
import uk.me.parabola.mkgmap.osmstyle.StyledConverter;
import uk.me.parabola.mkgmap.osmstyle.eval.SyntaxException;
import uk.me.parabola.util.EnhancedProperties;

import static java.nio.channels.FileChannel.MapMode.READ_ONLY;

/**
 * Create contour lines using an algorithm similar to that described in 
 * <a href="http://mapcontext.com/autocarto/proceedings/auto-carto-5/pdf/an-adaptive-grid-contouring-algorithm.pdf">An Adaptive Grid Contouring Algorithm</a> by Downing and Zoraster. 
 */
public abstract class DEM 
{
    final static double epsilon = 1e-9;
    final static double delta = 1.5;
    final static int maxPoints=200000;
    final static double minDist = 15;
    final static double maxDist = 21;
    final static double step = 0.01;
    final static double semiMajorAxis = 6378137.0;
    final static double inverseFlattening = 298.257223563;

    static int M=1200;
    static int N=M;
    static double res = 1.0/N;
    static int id = -1;

    short values[] = null;
    int lat;
    int lon;
    
    abstract double ele(int x, int y);
    abstract void read(int minLon, int minLat, int maxLon, int maxLat);
    abstract void serializeCopyRight(Writer out) throws IOException;

    public static void createContours(LoadableMapDataSource mapData, EnhancedProperties config)
    {
	Area bounds = mapData.getBounds();

	double minLat = Utils.toDegrees(bounds.getMinLat());
	double minLon = Utils.toDegrees(bounds.getMinLong());
	double maxLat = Utils.toDegrees(bounds.getMaxLat());
	double maxLon = Utils.toDegrees(bounds.getMaxLong());
	
	System.out.printf("bounds: %f %f %f %f\n", minLat, minLon, maxLat, maxLon);
	DEM data;
	String demType = config.getProperty("dem-type", "CGIAR");

	String dataPath;
	if (demType.equals("ASTER")) {
	    dataPath = config.getProperty("dem-path", "ASTER");
	    data = new ASTERGeoTiff(dataPath, minLat, minLon, maxLat, maxLon);
	}
	else if (demType.equals("CGIAR")) {
	    dataPath = config.getProperty("dem-path", "CGIAR");
	    data = new CGIARGeoTiff(dataPath, minLat, minLon, maxLat, maxLon);
	}
	else {
	    dataPath = config.getProperty("dem-path", "SRTM");
	    data = new HGTDEM(dataPath, minLat, minLon, maxLat, maxLon);
	}
	
	Isolines lines = data.new Isolines(data, minLat, minLon, maxLat, maxLon);
	int increment = config.getProperty("dem-increment", 10);

	double minHeight = lines.getMinHeight();
	double maxHeight = lines.getMaxHeight();
	int maxLevels = config.getProperty("dem-maxlevels", 100);
	while ((maxHeight-minHeight)/increment > maxLevels)
	    increment *= 2;

	String loc = config.getProperty("style-file");
	if (loc == null)
	    loc = config.getProperty("map-features");
	String name = config.getProperty("style");
	
	if (loc == null && name == null)
	    name = "default";

	LoadableMapDataSource dest = mapData;
	if (config.getProperty("dem-separate-img", false)) {
	    dest = new DEMMapDataSource(mapData, config);
	}

	
	OsmConverter converter;
	try {
	    Style style = new StyleImpl(loc, name);
	    style.applyOptionOverride(config);
	    
	    converter = new StyledConverter(style, ((MapperBasedMapDataSource) dest).getMapper(), config);
	} catch (SyntaxException e) {
	    System.err.println("Error in style: " + e.getMessage());
	    throw new ExitException("Could not open style " + name);
	} catch (FileNotFoundException e) {
	    String name1 = (name != null)? name: loc;
	    throw new ExitException("Could not open style " + name1);
	}
	
	for (int level=0; level<maxHeight; level+=increment) {
	    if (level < minHeight) continue;

	    // create isolines
	    lines.addLevel(level);
	
    	    for (Isolines.Isoline line : lines.isolines) {
		Way way = new Way(id--, line.points);
		way.addTag("contour", "elevation");
		way.addTag("ele", String.format("%d", (int) line.level));
		converter.convertWay(way);
	    }
	    lines.isolines.clear();
	}
	
	if (config.getProperty("dem-separate-img", false)) {
	    MapBuilder builder = new MapBuilder();
	    builder.config(config);

	    FileSystemParam params = new FileSystemParam();
	    params.setMapDescription("contour lines");
	    long mapName = Integer.valueOf(config.getProperty("mapname", "63240000"));
	    try {
		Map map = Map.createMap(String.format("%08d", mapName+10000000), params);
		builder.makeMap(map, dest);
		map.close();
	    }
	    catch (Exception ex) {
		throw new RuntimeException(ex);
	    }
	}
    }

    public static class HGTDEM extends DEM
    {
	MappedByteBuffer buffer = null;

	public HGTDEM(String dataPath, double minLat, double minLon, double maxLat, double maxLon)
	{
	    this.lat = (int) minLat;
	    this.lon = (int) minLon;
	    if (maxLat > lat+1 || maxLon > lon+1)
		throw new RuntimeException("Area too large (must not span more than one SRTM file)");

	    String northSouth = lat < 0 ? "S" : "N";
	    String eastWest = lon > 0 ? "E" : "W";
	    String fileName = String.format("%s/%s%02d%s%03d.hgt", dataPath, 
					    northSouth, lat < 0 ? -lat : lat, 
					    eastWest, lon < 0 ? -lon : lon);
	    try {
		FileInputStream is = new FileInputStream(fileName);
		buffer = is.getChannel().map(READ_ONLY, 0, 2*(M+1)*(M+1));
	    }
	    catch (Exception e) {
		throw new RuntimeException(e);
	    }
	}

	public  void read(int minLon, int minLat, int maxLon, int maxLat)
	{
	}

	public double ele(int x, int y)
	{
	    return buffer.getShort(2*((M-y)*(M+1)+x))+delta;
	}

	public void serializeCopyRight(Writer out) throws IOException
	{
	    out.write("  <copyright>\n");
	    out.write("  Contour lines generated from DEM data by NASA\n");
	    out.write("  </copyright>\n");
	}
    }

    public static class CGIARGeoTiff extends DEM
    {
	Raster raster;
	String fileName;
	int minLat, minLon, maxLat, maxLon;
	PlanarImage image;

	public CGIARGeoTiff(String dataPath, double minLat, double minLon, double maxLat, double maxLon)
	{
	    this.lat = ((int) (minLat/5))*5;
	    this.lon = ((int) (minLon/5))*5;
	    if (maxLat > lat+5 || maxLon > lon+5)
		throw new RuntimeException("Area too large (must not span more than one CGIAR GeoTIFF)");

	    int tileX, tileY;
	    tileX = (180 + lon)/5 + 1;
	    tileY = (60 - lat)/5;
	    this.fileName = String.format("%s/srtm_%02d_%02d.tif", dataPath, tileX, tileY);
	    init();
	}

	public void serializeCopyRight(Writer out) throws IOException
	{
	    out.write("  <copyright>\n");
	    out.write("  Contour lines generated from improved SRTM data by CIAT-CSI (see http://srtm.csi.cgiar.org)\n");
	    out.write("  </copyright>\n");
	}
	
	public  void read(int minLon, int minLat, int maxLon, int maxLat)
	{
	    this.minLon = minLon;
	    this.minLat = minLat;
	    this.maxLon = maxLon;
	    this.maxLat = maxLat;
	    raster = image.getData(new java.awt.Rectangle(minLon, 6000-maxLat-1, maxLon-minLon+1, maxLat-minLat+1));
	    System.out.printf("read: %d %d %d %d\n", minLon, 6000-maxLat-1, maxLon-minLon+1, maxLat-minLat+1);
	}

	void init()
	{
	    System.out.printf("CGIAR GeoTIFF: %s\n", fileName);
	    N = 6000;
	    M = 6000;
	    res = 5.0/M;
	    
	    try {
		SeekableStream s = new FileSeekableStream(fileName);
		ParameterBlock pb = new ParameterBlock();
		pb.add(s);
		
		TIFFDecodeParam param = new TIFFDecodeParam();
		pb.add(param);
		
		RenderedOp op = JAI.create("tiff", pb);
		image = op.createInstance();
		System.out.printf("Image: %d %d %d %d\n", image.getWidth(), image.getHeight(), 
				  image.getNumXTiles(), image.getNumYTiles());
	    }
	    catch (Exception e) {
		throw new RuntimeException(e);
	    }
	}

	public double ele(int x, int y)
	{
	    try {
		int elevation = raster.getPixel(x, 6000-y-1, (int[])null)[0];
		return elevation+delta;
	    }
	    catch (ArrayIndexOutOfBoundsException ex) {
		System.out.printf("ele: (%d, %d) (%d, %d, %d, %d)  %s\n", 
				  x, 6000-y-1, 
				  raster.getMinX(), raster.getMinY(), 
				  raster.getWidth(), raster.getHeight(), ex.toString());
		throw ex;
	    }
	}
    }

    public static class ASTERGeoTiff extends DEM
    {
	Raster raster;
	String fileName;
	int minLat, minLon, maxLat, maxLon;
	PlanarImage image;

	public ASTERGeoTiff(String dataPath, double minLat, double minLon, double maxLat, double maxLon)
	{
	    this.lat = (int) minLat;
	    this.lon = (int) minLon;
	    if (maxLat > lat+1 || maxLon > lon+1)
		throw new RuntimeException("Area too large (must not span more than one ASTER GeoTIFF)");

	    String northSouth = lat < 0 ? "S" : "N";
	    String eastWest = lon > 0 ? "E" : "W";
	    fileName = String.format("%s/ASTGTM_%s%02d%s%03d_dem.tif", dataPath, 
				     northSouth, lat < 0 ? -lat : lat, 
				     eastWest, lon < 0 ? -lon : lon);
	    init();
	}

	public void serializeCopyRight(Writer out) throws IOException
	{
	    out.write("  <copyright>\n");
	    out.write("  Contour lines generated from DGM data by ASTER (see https://wist.echo.nasa.gov/~wist/api/imswelcome)\n");
	    out.write("  </copyright>\n");
	}
	
	public  void read(int minLon, int minLat, int maxLon, int maxLat)
	{
	    this.minLon = minLon;
	    this.minLat = minLat;
	    this.maxLon = maxLon;
	    this.maxLat = maxLat;
	    raster = image.getData(new java.awt.Rectangle(minLon, 3601-maxLat-1, maxLon-minLon+1, maxLat-minLat+1));
	    System.out.printf("read: %d %d %d %d\n", minLon, 3601-maxLat-1, maxLon-minLon+1, maxLat-minLat+1);
	}

	void init()
	{
	    System.out.printf("ASTER GeoTIFF: %s\n", fileName);
	    N = 3600;
	    M = 3600;
	    res = 1.0/M;
	    
	    try {
		SeekableStream s = new FileSeekableStream(fileName);
		ParameterBlock pb = new ParameterBlock();
		pb.add(s);
		
		TIFFDecodeParam param = new TIFFDecodeParam();
		pb.add(param);
		
		RenderedOp op = JAI.create("tiff", pb);
		image = op.createInstance();
		System.out.printf("Image: %d %d %d %d\n", image.getWidth(), image.getHeight(), 
				  image.getNumXTiles(), image.getNumYTiles());
	    }
	    catch (Exception e) {
		throw new RuntimeException(e);
	    }
	}

	public double ele(int x, int y)
	{
	    try {
		int elevation = raster.getPixel(x, 3601-y-1, (int[])null)[0];
		return elevation+delta;
	    }
	    catch (ArrayIndexOutOfBoundsException ex) {
		System.out.printf("ele: (%d, %d) (%d, %d, %d, %d)  %s\n", 
				  x, 3601-y-1, 
				  raster.getMinX(), raster.getMinY(), 
				  raster.getWidth(), raster.getHeight(), ex.toString());
		throw ex;
	    }
	}
    }

    private static void debug(String format, Object ... args)
    {
	//System.out.printf(format + "\n", args);
    }


    int lastXi = -1;
    int lastYi = -1;

    final static int bcInv[][]= {
	{ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, 
	{ -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1, 0, 0, 0, 0},
	{ 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0}, 
	{ 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, 
	{ 0, 0, 0, 0, -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1}, 
	{ 0, 0, 0, 0, 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1}, 
	{ -3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{ 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0}, 
	{ 9, -9, 9, -9, 6, 3, -3, -6, 6, -6, -3, 3, 4, 2, 1, 2}, 
	{ -6, 6, -6, 6, -4, -2, 2, 4, -3, 3, 3, -3, -2, -1, -1, -2}, 
	{ 2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{ 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0}, 
	{ -6, 6, -6, 6, -3, -3, 3, 3, -4, 4, 2, -2, -2, -2, -1, -1}, 
	{ 4, -4, 4, -4, 2, 2, -2, -2, 2, -2, -2, 2, 1, 1, 1, 1}
    };

    static int lastId=1000000000;
    static double lastX = 0;
    static double lastY = 0;

    static int edge[] = new int[2];
    
    static int off0[][] = {{ 0,  0},
			   { 0,  0},
			   { 0,  1},
			   { 1,  1}};
    
    static int off1[][] = {{ 1,  0},
			   { 0,  1},
			   { 1,  1},
			   { 1,  0}};
    
    static int brd[] = { 1, 2, 4, 8 };
    static int inv[] = { 4, 8, 1, 2 };
    
    static int rev[] = { 2, 3, 0, 1 };
    
    static int mov[][] = {{ 0, -1},
			  {-1,  0},
			  { 0,  1},
			  { 1,  0}};


    double bc[][] = new double[4][4];
    double bc_y[] = new double[4];
    double bc_y1[] = new double[4];
    double bc_y2[] = new double[4];
    double bc_y12[] = new double[4];
    double bc_Coeff[] = new double[16];
    double bc_x[] = new double[16];

    private void recalculateCoefficients(int xi, int yi)
    {
	double v00, vp0, v0p, vpp;
	double vm0, v0m, vpm, vmp, vmm, vmP, vPm; 
	double vP0, v0P, vPp, vpP, vPP;

	v00 = ele(xi, yi);
	v0p = ele(xi, yi+1);
	vpp = ele(xi+1, yi+1);
	vp0 = ele(xi+1, yi);

	vm0 = ele(xi-1, yi);
	v0m = ele(xi, yi-1);
	vmp = ele(xi-1, yi+1);
	vpm = ele(xi+1, yi-1);
	vmm = ele(xi-1, yi-1);
	vmP = ele(xi+2, yi-1);
	vPm = ele(xi-1, yi+2);

	vP0 = ele(xi+2, yi);
	v0P = ele(xi, yi+2);
	vPp = ele(xi+2, yi+1);
	vpP = ele(xi+1, yi+2);
	vPP = ele(xi+2, yi+2);

	bc_y[0] = v00;
	bc_y[1] = vp0;
	bc_y[2] = vpp;
	bc_y[3] = v0p;

	bc_y1[0] = (vp0-vm0)/2;
	bc_y1[1] = (vP0-v00)/2;
	bc_y1[2] = (vPp-v0p)/2;
	bc_y1[3] = (vpp-vmp)/2;

	bc_y2[0] = (v0p-v0m)/2;
	bc_y2[1] = (vpp-vpm)/2;
	bc_y2[2] = (vpP-vp0)/2;
	bc_y2[3] = (v0P-v00)/2;

	bc_y12[0] = (vpp - vpm - vmp + vmm) / 4;
	bc_y12[0] = (vPp - vPm - v0p + v0m) / 4;
	bc_y12[2] = (vPP - vP0 - v0P + v00) / 4;
	bc_y12[0] = (vpP - vp0 - vmP + vm0) / 4;

	int j, i;
	double s;

	for (i=0; i<4; i++) 
	{
	    bc_x[i]=bc_y[i];
	    bc_x[i+4]=bc_y1[i];
	    bc_x[i+8]=bc_y2[i];
	    bc_x[i+12]=bc_y12[i];
	}

	for (i=0; i<16; i++) 
	{
	    s = 0;
	    for (int k=0; k<16; k++) s += bcInv[i][k]*bc_x[k];
	    bc_Coeff[i] = s;
	}

	int l = 0;
	for (i=0; i<4; i++)
	    for (j=0; j<4; j++) 
		bc[i][j] = bc_Coeff[l++];
    }

    public double gradient(double lat, double lon, double[] grad)
    {
	grad[0] = 0;
	grad[1] = 0;

	double x = (lon-this.lon)/res;
	double y = (lat-this.lat)/res;
	
	int xi = (int) x;
	int yi = (int) y;

	if (lastXi != xi || lastYi != yi)
	{
	    debug("new Cell for interpolation: %d %d", xi, yi);
	    recalculateCoefficients(xi, yi);
	    lastXi = xi;
	    lastYi = yi;
	}

	double t = x - xi;
	double u = y - yi;

	if  (xi < 0 || xi > N+1 || yi < 0 || yi > N+1)
	    throw new IndexOutOfBoundsException(String.format("(%f, %f)->(%d, %d)", lat, lon, xi, yi));
	
	double val = 0;
	for (int i=3; i>=0; i--) 
	{
	    val = t*val + ((bc[i][3]*u + bc[i][2])*u + bc[i][1])*u + bc[i][0];
	    grad[0] = u*grad[0] + (3*bc[3][i]*t + 2*bc[2][i])*t + bc[1][i];
	    grad[1] = t*grad[1] + (3*bc[i][3]*t + 2*bc[i][2])*t + bc[i][1];
	}

	return val;
    }

    public double elevation(double lat, double lon)
    {
	double x = (lon-this.lon)/res;
	double y = (lat-this.lat)/res;
	
	int xi = (int) x;
	int yi = (int) y;

	if (lastXi != xi || lastYi != yi)
	{
	    debug("new Cell for interpolation: %d %d", xi, yi);
	    recalculateCoefficients(xi, yi);
	    lastXi = xi;
	    lastYi = yi;
	}

	double t = x - xi;
	double u = y - yi;

	if  (xi < 0 || xi > N+1 || yi < 0 || yi > N+1)
	    throw new IndexOutOfBoundsException(String.format("(%f, %f)->(%d, %d)", lat, lon, xi, yi));
	
	double val = 0;
	for (int i=3; i>=0; i--) 
	{
	    val = t*val + ((bc[i][3]*u + bc[i][2])*u + bc[i][1])*u + bc[i][0];
	}

	return val;
    }

    public double elevation(int x, int y)
    {
	if (x < 0 || x > N || y < 0 || y > N)
	    throw new IndexOutOfBoundsException(String.format("elevation: %d %d", x, y));
	return ele(x, y);
    }


    class Isolines
    {
	DEM data;
	int minX;
	int maxX; 
	int minY; 
	int maxY;

	double min;
	double max;

	Vector<Isoline> isolines = new Vector<Isoline>();

	class Isoline
	{
	    int id;
	    Vector<Coord> points;
	    double level;
	    boolean isClosed;
	    
	    private Isoline(double level)
	    {
		this.level = level;
		isClosed = false;
		id = lastId++;
		points = new Vector<Coord>();
	    }

	    private class Edge implements Brent.Function
	    {
		double x0, y0, x1, y1;
		Edge(double x0, double y0, double x1, double y1)
		{
		    this.x0 = x0;
		    this.y0 = y0;
		    this.x1 = x1;
		    this.y1 = y1;
		} 
		
		public double eval(double d)
		{
		    double f = data.elevation(x0 + d * (x1-x0), y0 + d * (y1-y0)) - level; 
		    //System.out.printf("evalEdge: %f %f\n", d, f); 
		    return f;
		}
	    }

	    private class FN implements Brent.Function
	    {
		double x0, y0;
		double dx, dy;

		public void setParameter(double x0, double y0, double dx, double dy)
		{
		    this.x0 = x0;
		    this.y0 = y0;
		    this.dx = dx;
		    this.dy = dy;
		}

		public double eval(double t) 
		{ 
		    double f = data.elevation(y0+t*dy, x0+t*dx) - level; 
		    return f; 
		}
	    }
	    
	    private FN fn = new FN();

	    double grad[] = new double[2];
	    double px[] = new double[4];
	    double py[] = new double[4];
	    int edges[] = new int[4];

	    boolean addCell(Position p, int direction)
	    {
		debug("addCell: %f %d %d %d %d", level, p.ix, p.iy, p.edge, direction);
		
		int c = 0;
		for (int k=0; k<4; k++)
		{
		    if (k == p.edge)
			continue;

		    int x0 = p.ix + off0[k][0];
		    int y0 = p.iy + off0[k][1];
		    int x1 = p.ix + off1[k][0];
		    int y1 = p.iy + off1[k][1];
		    
		    double l0 = elevation(x0, y0) - level;
		    double l1 = elevation(x1, y1) - level;
		    
		    if (Math.abs(l1) < epsilon || l0*l1 < 0)
		    {
			edges[c] = k;

			Brent.Function f = new Edge(data.lat + y0 * DEM.res, data.lon + x0 * DEM.res, data.lat + y1 * DEM.res, data.lon + x1 * DEM.res);
			double f0 = elevation(x0, y0) - level;
			double f1 = elevation(x1, y1) - level;
			double delta;
			
			if (Math.abs(1) < epsilon)
			{
			    delta = 1;
			}
			else if (Math.abs(f0) < epsilon)
			    throw new RuntimeException("implementation error!");
			else
			    delta = Brent.zero(f, epsilon, 1-epsilon);
			
			px[c] = data.lon + (x0+delta*(x1-x0))*DEM.res;
			py[c] = data.lat + (y0+delta*(y1-y0))*DEM.res;
			c++;
		    }
		}


		if (c == 1)
		{
		    p.edge = edges[0];

		    double px0 = p.x;
		    double py0 = p.y;
		    p.x = px[0];
		    p.y = py[0];
		    double px1 = p.x;
		    double py1 = p.y;
		    
		    double xMin = data.lon + p.ix * DEM.res;
		    double xMax = xMin + DEM.res;
		    double yMin = data.lat + p.iy * DEM.res;
		    double yMax = yMin + DEM.res;

		    refineAdaptively(xMin, yMin, xMax, yMax, px0, py0, px1, py1, direction, maxDist);

		    addPoint(p.x, p.y, direction);
		    p.moveCell();
		    return true;
		}
		else
		{
		    debug("addCellByStepping: %d", c);
		    return addCellByStepping(p, direction, c, edges, px, py);
		}		
	    }

	    private void refineAdaptively(double xMin, double yMin, double xMax, double yMax, 
					  double x0, double y0, double x1, double y1, 
					  int direction, double maxDist)
	    {
		double dist = quickDistance(x0, y0, x1, y1);
		if (dist > maxDist)
		{
		    double dx = x1-x0;
		    double dy = y1-y0;
		    double xm, ym, f0, f1, t0, t1, t;
		    Brent.Function f;
		    
		    xm = x0 + 0.5*dx;
		    ym = y0 + 0.5*dy;
		    double n = Math.sqrt(dx*dx+dy*dy);
		    fn.setParameter(xm, ym, -dy/n, dx/n);
		    f = fn;
		    t0 = -0.05*res;
		    t1 = 0.05*res;
		    f0 = f.eval(t0);
		    f1 = f.eval(t1);
		
		    int count = 0;
		    while (f0 * f1 > 0 && count++ < 20) {
			if ((count & 1) > 0)
			    t0 -= 0.05*res;
			else
			    t1 += 0.05*res; 
			f0 = f.eval(t0);
			f1 = f.eval(t1);
			debug("refine: %f %f %f %f", t0, t1, f0, f1);
		    }
			
		    if (f0 * f1 < 0)
		    {
			t = Brent.zero(f, t0, t1);
			xm -= t*dy;
			ym += t*dx;
		    }
		    else
		    {
			debug("refine failed: %f %f %f %f", t0, t1, f0, f1);
			if (false) throw new RuntimeException(String.format("refine failed: %f %f %f %f %f %f %f %f", xMin, yMin, xMax, yMax, x0, y0, x1, y1));
			return;
		    }

		    if (xm > xMin && xm < xMax && ym > yMin && ym < yMax)
			refineAdaptively(xMin, yMin, xMax, yMax, x0, y0, xm, ym, direction, maxDist*1.1);
		    addPoint(xm, ym, direction);
		    if (xm > xMin && xm < xMax && ym > yMin && ym < yMax)
			refineAdaptively(xMin, yMin, xMax, yMax, xm, ym, x1, y1, direction, maxDist*1.1);
		}
	    }

	    boolean addCellByStepping(Position p, int direction, int numEdges, int[] edges, double[] px, double[] py)
	    {
		debug("addCellByStepping: %f %d %d %d %d", level, p.ix, p.iy, p.edge, direction);


		double xMin = data.lon + p.ix * DEM.res;
		double xMax = xMin + DEM.res;
		double yMin = data.lat + p.iy * DEM.res;
		double yMax = yMin + DEM.res;

		double dt, t0, t1;
		double f0, f1;
		boolean edgeHit = false;
		double h[] = new double[4];
		
		int dir;
		double n2 = Math.sqrt(1.0/(grad[0]*grad[0] + grad[1]*grad[1]));
		double dx;
		double dy;

		int count=0;
		int iMin = -1;
		double md = 5000;
		
		for (int i=0; i<numEdges; i++) {
		    gradient(p.y, p.x, grad);
		    double dist = quickDistance(p.x, p.y, px[i], py[i]);
		    debug("distance %d: %f", i, dist);
		    
		    if (dist < md && (visited[p.iy*(N+1)+p.ix] & brd[edges[i]]) == 0) {
			md = dist;
			iMin = i;
		    }
		}

		p.edge = edges[iMin];
		
		double px0 = p.x;
		double py0 = p.y;
		p.x = px[iMin];
		p.y = py[iMin];
		double px1 = p.x;
		double py1 = p.y;
		
		xMin = data.lon + p.ix * DEM.res;
		xMax = xMin + DEM.res;
		yMin = data.lat + p.iy * DEM.res;
		yMax = yMin + DEM.res;
		
		refineAdaptively(xMin, yMin, xMax, yMax, px0, py0, px1, py1, direction, maxDist);
		
		addPoint(p.x, p.y, direction);
		p.moveCell();
		return true;
	    }

	    private void addPoint(double x, double y, int direction)
	    {
		double dist = quickDistance(x, y, lastX, lastY);
		debug("addPoint: %f %f %f", x, y, dist);
		
		if (dist > minDist)
		{
		    if (direction > 0)
			points.add(0, new Coord(y, x));
		    else
			points.add(points.size(), new Coord(y, x));
		    lastX = x;
		    lastY = y;
		}
	    }

	    private void close()
	    {
		points.add(points.size(), points.get(0));
		isClosed = true;
	    }

	    private void addMove(int x0, int y0, int x1, int y1, int direction)
	    {
		double x;
		double y;
		
		if (x0 < minX || x0 >= maxX || y0 < minY || y0 >= maxY)
		    return;
		double l0 = data.elevation(x0, y0);
		double l1 = data.elevation(x1, y1);
		debug("addMove: %d %d %d %d %.2f %.2f %d", x0, y0, x1, y1, l0, l1, direction);
		
		if ((l0 < level && l1 < level) || (l0 > level && l1 > level))
		    throw new RuntimeException("implementation error");
		
		if (l0 == level)
		{
		    x = data.lon + x0 * DEM.res;
		    y = data.lat + y0 * DEM.res;
		}
		else
		{
		    double delta = (l0-level)/(l0-l1);
		    x = data.lon + (x0 + delta * (x1-x0)) * DEM.res;
		    y = data.lat + (y0 + delta * (y1-y0)) * DEM.res;
		}
		double dist = quickDistance(x, y, lastX, lastY);
		debug("levels: %d %d %f, %d %d %f: %f %f %f", x0, y0, l0, x1, y1, l1, x, y, dist);
		
		if (dist > 1)
		{
		    if (direction > 0)
			points.add(0, new Coord(y, x));
		    else
			points.add(points.size(), new Coord(y, x));
		    lastX = x;
		    lastY = y;
		}
	    }
	}
	 
	public Isolines(DEM data, int minX, int maxX, int minY, int maxY)
	{
	    this.data = data;
	    this.minX = minX;
	    this.maxX = maxX;
	    this.minY = minY;
	    this.maxY = maxY;
	    
	    init();
	}

	public Isolines(DEM data, double  minLat, double minLon, double maxLat, double maxLon)
	{
	    System.out.printf("init: %f %f %f %f\n", minLat, minLon, maxLat, maxLon);

	    this.data = data;
	    this.minX = (int) ((minLon-data.lon)/data.res);
	    this.minY = (int) ((minLat-data.lat)/data.res);
	    this.maxX = (int) ((maxLon-data.lon)/data.res);
	    this.maxY = (int) ((maxLat-data.lat)/data.res);
	    
	    init();
	}

	private void init()
	{
	    System.out.printf("init: %d %d %d %d\n", minX, minY, maxX, maxY);
	    data.read(minX-2, minY-2, maxX+2, maxY+2); 
	    // we need some overlap for bicubic interpolation
	    max = -1000;
	    min = 10000;
	    for (int i=minX; i<maxX; i++)
		for (int j=minY; j<maxY; j++)
		{
		    if (data.elevation(i, j) < min) min = data.elevation(i, j);
		    if (data.elevation(i, j) > max) max = data.elevation(i, j);
		}
	
	    debug("min: %f, max: %f\n", min, max);
	}
	
	double getMinHeight()
	{
	    return min;
	}

	double getMaxHeight()
	{
	    return max;
	}

	public void addLevels(int start, int increment)
	{
	    for (int level=start; level<max; level+=increment)
		addLevel(level);
	}

	private class Edge implements Brent.Function
	{
	    double x0, y0, x1, y1, level;
	    Edge(double x0, double y0, double x1, double y1, double level)
	    {
		this.x0 = x0;
		this.y0 = y0;
		this.x1 = x1;
		this.y1 = y1;
		this.level = level;
	    } 
	    
	    public double eval(double d)
	    {
		double f = data.elevation(x0 + d * (x1-x0), y0 + d * (y1-y0)) - level; 
		//System.out.printf("evalEdge: %f %f\n", d, f); 
		return f;
	    }
	}

	private class Position
	{
	    int ix, iy;
	    double x, y;
	    int edge;

	    Position(int ix, int iy, double x, double y, int edge)
	    {
		this.ix = ix;
		this.iy = iy;
		this.x = x;
		this.y = y;
		this.edge = edge;
	    } 

	    Position(Position p)
	    {
		this.ix = p.ix;
		this.iy = p.iy;
		this.x = p.x;
		this.y = p.y;
		this.edge = p.edge;
	    } 

	    void markEdge()
	    {
		debug("marking edge: %d %d %d %d", ix, iy, edge, brd[edge]);
		visited[iy*(N+1)+ix] |= brd[edge];
	    }

	    void moveCell()
	    {
		markEdge();
		ix += mov[edge][0];
		iy += mov[edge][1];
		edge = rev[edge];
		markEdge();
	    }
	}

	byte visited[] = new byte[(N+1)*(N+1)];

	public void addLevel(double level)
	{
	    if (level < min || level > max)
		return;
	    
	    System.out.printf("addLevel: %f\n", level);
	    java.util.Arrays.fill(visited, (byte) 0);

	    for (int y=minY; y<maxY; y++)
	    {
		for (int x=minX; x<maxX; x++)
		{
		    byte v = 0;
		    // Mark the borders of the cell, represented by the four points (i, j), (i+1, j), (i, j+1), (i+1, j+1), 
		    // which are intersected by the contour. The values are:
		    // 1: top
		    // 2: left
		    // 4: bottom
		    // 8: right
		    if (data.elevation(x, y) >= level)
		    {
			if (data.elevation(x+1, y) < level) { v |= 1; }
			if (data.elevation(x, y+1) < level) { v |= 2; }
		    }
		    else
		    {
			if (data.elevation(x+1, y) > level) { v |= 1; }
			if (data.elevation(x, y+1) > level) { v |= 2; }
		    }

		    int k=-1;
		   
		    if ((v&1) > 0 && (visited[y*(N+1)+x]&1) == 0)
		    {
			k=0;
		    }
		    else if ((v&2) > 0 && (visited[y*(N+1)+x]&2) == 0)
		    {
			k=1;
		    }

		    if (k>=0)
		    {
			int x0 = x + off0[k][0];
			int y0 = y + off0[k][1];
			int x1 = x + off1[k][0];
			int y1 = y + off1[k][1];

			try {
			    Brent.Function f = new Edge(data.lat + y0 * DEM.res, data.lon + x0 * DEM.res, 
							data.lat + y1 * DEM.res, data.lon + x1 * DEM.res, 
							level);
			    double f0 = elevation(x0, y0) - level;
			    double f1 = elevation(x1, y1) - level;
			    double delta;
			    if (Math.abs(f0) < epsilon)
			    {
				delta = 0;
			    }
			    else if (Math.abs(f1) < epsilon)
				continue;
			    else
				delta = Brent.zero(f, 0, 1-epsilon);
			    
			    Position p = new Position(x, y, data.lon + (x0+delta*(x1-x0))*DEM.res, data.lat + (y0+delta*(y1-y0))*DEM.res, k);
			    p.markEdge();
			    isolines.add(traceByStepping(level, p));
			}
			catch (RuntimeException ex)
			{
			    debug("error: %s", ex.toString());
			    ex.printStackTrace();
			    return;
			}
		    }
		}
	    }
	}

	private Isoline traceByStepping(double level, Position p)
	{
	    debug("traceByStepping: starting contour %f %d %d %f %f %d", level, p.ix, p.iy, p.x, p.y, p.edge);
	    int direction = 1;
	    int n = 0;
	    Position startP = new Position(p);
	    boolean foundEnd = false;

	    Isoline line = new Isoline(level);

	    while (true)
	    {
		debug("traceByStepping: %f %d %d %f %f %d", level, p.ix, p.iy, p.x, p.y, p.edge);
		visited[p.iy*(N+1)+p.ix] |= brd[p.edge];

		if (n>0 && p.ix == startP.ix && p.iy == startP.iy && quickDistance(p.x, p.y, startP.x, startP.y) < 5)
		{
		    debug("closed curve!");
		    line.close();
		    break;
		}
		else if (p.ix < minX || p.iy < minY || p.ix >= maxX || p.iy >= maxY)
		{
		    if (foundEnd) // did we already reach one end?
		    {
			debug("second border reached!");
			break;
		    }
		    else
		    {
			debug("border reached!");
			foundEnd = true;
			n = 0;
			direction *= -1;
			p = new Position(startP);
			p.moveCell();
			continue;
		    }
		}
		n++;
		if (!line.addCell(p, direction) || line.points.size() > maxPoints)
		{
		    debug("ending contour");
		    isolines.add(line);
		    return line; 
		}
	    }

	    return line;
	}

    }

    public static double quickDistance(double long1, double lat1, double long2, double lat2) 
    {
	double latDiff;
	if (lat1 < lat2)
	    latDiff = lat2 - lat1;
	else
	    latDiff = lat1 - lat2;	
	if (latDiff > 90)
	    latDiff -= 180;
	
	double longDiff;
	if (long1 < long2)
	    longDiff = long2 - long1;
	else
	    longDiff = long1 - long2;
	if (longDiff > 180)
	    longDiff -= 360;
	
	// scale longDiff by cosine of average latitude
	longDiff *= Math.cos(Math.PI / 180 * Math.abs((lat1 + lat2) / 2));
	
	double distDegSq = (latDiff * latDiff) + (longDiff * longDiff);

	return 40075000 * Math.sqrt(distDegSq) / 360;
    }

    /**
     * Returns the orthodromic distance between two geographic coordinates in WGS84 datum.
     * The orthodromic distance is the shortest distance between two points
     * on a sphere's surface. The orthodromic path is always on a great circle.
     *
     * @param  x1 Longitude of first  point (in degrees).
     * @param  y1 Latitude  of first  point (in degrees).
     * @param  x2 Longitude of second point (in degrees).
     * @param  y2 Latitude  of second point (in degrees).
     * @return The orthodromic distance (in meters).
     *
     */
    public static double orthodromicDistance(double x1, double y1, double x2, double y2) 
    {
        x1 = Math.toRadians(x1);
        y1 = Math.toRadians(y1);
        x2 = Math.toRadians(x2);
        y2 = Math.toRadians(y2);
        final int MAX_ITERATIONS = 100;
        final double EPS = 5.0E-14;
        final double F = 1 / inverseFlattening;
        final double R = 1 - F;
        double tu1 = R * Math.sin(y1) / Math.cos(y1);
        double tu2 = R * Math.sin(y2) / Math.cos(y2);
        double cu1 = 1 / Math.sqrt(tu1 * tu1 + 1);
        double cu2 = 1 / Math.sqrt(tu2 * tu2 + 1);
        double su1 = cu1 * tu1;
        double s = cu1 * cu2;
        double baz = s * tu2;
        double faz = baz * tu1;
        double x = x2 - x1;
        for (int i = 0; i < MAX_ITERATIONS; i++) 
	{
            final double sx = Math.sin(x);
            final double cx = Math.cos(x);
            tu1 = cu2 * sx;
            tu2 = baz - su1 * cu2 * cx;
            final double sy = Math.sqrt(tu1*tu1 + tu2*tu2);;
            final double cy = s * cx + faz;
            final double y = Math.atan2(sy, cy);
            final double SA = s * sx / sy;
            final double c2a = 1 - SA * SA;
            double cz = faz + faz;
            if (c2a > 0) {
                cz = -cz / c2a + cy;
            }
            double e = cz * cz * 2 - 1;
            double c = ((-3 * c2a + 4) * F + 4) * c2a * F / 16;
            double d = x;
            x = ((e * cy * c + cz) * sy * c + y) * SA;
            x = (1 - c) * x * F + x2 - x1;
            if (Math.abs(d - x) <= EPS) {
                x = Math.sqrt((1 / (R * R) - 1) * c2a + 1) + 1;
                x = (x - 2) / x;
                c = 1 - x;
                c = (x * x / 4 + 1) / c;
                d = (0.375 * x * x - 1) * x;
                x = e * cy;
                s = 1 - 2 * e;
                s = ((((sy * sy * 4 - 3) * s * cz * d / 6 - x) * d / 4 + cz) * sy * d + y) * c * R * semiMajorAxis;
                return s;
            }
        }
        final double LEPS = 1.0E-10;
        if (Math.abs(x1 - x2) <= LEPS && Math.abs(y1 - y2) <= LEPS) {
            return 0;
        }
        if (Math.abs(y1) <= LEPS && Math.abs(y2) <= LEPS) {
            return Math.abs(x1 - x2) * semiMajorAxis;
        }
        throw new ArithmeticException("no convergence");
    }

    private static class DEMMapDataSource extends MapperBasedMapDataSource implements LoadableMapDataSource
    {
	LoadableMapDataSource parent;
	List<String> copyright = new ArrayList<String>();

	DEMMapDataSource(LoadableMapDataSource parent, EnhancedProperties props)
	{
	    this.parent = parent;
	    config(props);
	}

	public boolean isFileSupported(String name)
	{
	    return false;
	}

	public void load(String name)
	    throws FileNotFoundException, FormatException
	{
	    throw new FormatException("load not supported");
	}

	public LevelInfo[] mapLevels()
	{
	    return parent.mapLevels();
	}

	public String[] copyrightMessages()
	{
	    return copyright.toArray(new String[1]);
	}
    }
}

/*
 * Copyright (C) 2009 Christian Gawron
 * 
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License version 2 as
 *  published by the Free Software Foundation.
 * 
 *  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.
 * 
 * 
 * Author: Christian Gawron
 * Create date: 03-Jul-2009
 */
package uk.me.parabola.mkgmap.reader.dem;

/**
 * Find zero of a function using Brent's method.
 */
public class Brent
{
    
    public interface Function
    {
	public double eval(double x);
    }


    static double epsilon = 3.0e-10;

    public static double zero(Function f, double x1, double x2) 
    {
	return zero(f, x1, x2, 1e-8, 100);
    }

    public static double zero(Function f, double x1, double x2, double tol, int maxit)
    {
	double a=x1, b=x2, c=x2, d=0, e=0, min1, min2;
	double fa=f.eval(a), fb=f.eval(b), fc;
	double p, q, r, s, tolerance, xm;

	if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
	    throw new ArithmeticException("Root must be bracketed");

	fc=fb;
	for (int iterations=0; iterations < maxit; iterations++) 
	{
	    if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
		c=a;
		fc=fa;
		e=d=b-a;
	    }
	    if (Math.abs(fc) < Math.abs(fb)) {
		a=b;
		b=c;
		c=a;
		fa=fb;
		fb=fc;
		fc=fa;
	    }
	    tolerance=2.0*epsilon*Math.abs(b)+0.5*tol;
	    xm=0.5*(c-b);
	    if (Math.abs(xm) <= tolerance || fb == 0.0) return b;
	    if (Math.abs(e) >= tolerance && Math.abs(fa) > Math.abs(fb)) {
		s=fb/fa;
		if (a == c) {
		    p=2.0*xm*s;
		    q=1.0-s;
		} else {
		    q=fa/fc;
		    r=fb/fc;
		    p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
		    q=(q-1.0)*(r-1.0)*(s-1.0);
		}
		if (p > 0.0) q = -q;
		p=Math.abs(p);
		min1=3.0*xm*q-Math.abs(tolerance*q);
		min2=Math.abs(e*q);
		if (2.0*p < (min1 < min2 ? min1 : min2)) {
		    e=d;
		    d=p/q;
		} else {
		    d=xm;
		    e=d;
		}
	    } else {
		d=xm;
		e=d;
	    }
	    a=b;
	    fa=fb;
	    if (Math.abs(d) > tolerance)
		b += d;
	    else
		b += xm >= 0 ? tolerance : -tolerance;
	    fb=f.eval(b);
	}
	throw new ArithmeticException("Maximum number of iterations exceeded");
    }
}

    