Changeset 25577 in osm for applications


Ignore:
Timestamp:
2011-03-12T11:04:59+01:00 (14 years ago)
Author:
glebius
Message:
  • Use Newton's method to solve reverse equation lat=f(y) in the Scanex's projection instead of dichotomy.
  • To speed up convergence we cache previous solution and use it as start for Newton's method.
  • If function gets out of bounds, we start with new random value between bounds.

Testing showed that usually it requires 4 - 5 iterations to achieve
0.000001 precision.

Submitted by: Andrey Boltenkov <anboltenkov yandex.ru>

File:
1 edited

Legend:

Unmodified
Added
Removed
  • applications/viewer/jmapviewer/src/org/openstreetmap/gui/jmapviewer/tilesources/ScanexTileSource.java

    r25369 r25577  
    11package org.openstreetmap.gui.jmapviewer.tilesources;
     2
     3import java.util.Random;
    24
    35import org.openstreetmap.gui.jmapviewer.OsmMercator;
     
    8183
    8284    /*
    83      * DIRTY HACK ALERT!
    84      *
    85      * Until I can't solve the equation, use dihotomy :(
     85     * To solve inverse formula latitude = f(y) we use
     86     * Newton's method. We cache previous calculated latitude,
     87     * because new one is usually close to the old one. In case
     88     * if solution gets out of bounds, we reset to a new random
     89     * value.
    8690     */
     91    private double cached_lat = 0;
     92
    8793    @Override
    8894    public double tileYToLat(int y, int zoom) {
    89         double lat = 0;
    90         double minl = OsmMercator.MIN_LAT;
    91         double maxl = OsmMercator.MAX_LAT;
    92         double c;
     95        Random r= new Random();
     96        double lat0, lat;
    9397
    94         for (int i=0; i < 60; i++) {
    95             c = latToTileY(lat, zoom);
    96             if (c < y) {
    97                 maxl = lat;
    98                 lat -= (lat - minl)/2;
    99             } else {
    100                 minl = lat;
    101                 lat += (maxl - lat)/2;
     98        lat = cached_lat;
     99        do {
     100            lat0 = lat;
     101            lat = lat - Math.toDegrees(NextTerm(Math.toRadians(lat), y, zoom));
     102            if (lat > OsmMercator.MAX_LAT || lat < OsmMercator.MIN_LAT) {
     103                lat = OsmMercator.MIN_LAT +
     104                    (double )r.nextInt((int )(OsmMercator.MAX_LAT -
     105                    OsmMercator.MIN_LAT));
    102106            }
    103         }
     107        } while ((Math.abs(lat0 - lat) > 0.000001));
    104108
    105         return lat;
     109        cached_lat = lat;
     110
     111        return (lat);
     112    }
     113
     114    /* Next term in Newton's polynomial */
     115    private double NextTerm(double lat, double y, int zoom) {
     116        double sinl=Math.sin(lat);
     117        double cosl=Math.cos(lat);
     118        double ec, f, df;
     119
     120        zoom = (int )Math.pow(2.0, zoom - 1);
     121        ec = Math.exp((1 - y/zoom)*Math.PI);
     122
     123        f = (Math.tan(Math.PI/4+lat/2) -
     124            ec * Math.pow(Math.tan(Math.PI/4 + Math.asin(E * sinl)/2), E));
     125        df = 1/(1 - sinl) - ec * E * cosl/((1 - E * sinl) *
     126            (Math.sqrt (1 - E * E * sinl * sinl)));
     127
     128        return (f/df);
    106129    }
    107130
Note: See TracChangeset for help on using the changeset viewer.