Examples with DoubleWrapper net.jafama.DoubleWrapper used on opensource projects

Example 1 with DoubleWrapper

use of net.jafama.DoubleWrapper in project elki by elki-project.

the class DistanceFunctionVisualization method drawCosine.

/**
 * Visualizes Cosine and ArcCosine distance functions
 *
 * @param svgp SVG Plot
 * @param proj Visualization projection
 * @param mid mean vector
 * @param angle Opening angle in radians
 * @return path element
 */
public static Element drawCosine(SVGPlot svgp, Projection2D proj, NumberVector mid, double angle) {
    // Project origin
    double[] pointOfOrigin = proj.fastProjectDataToRenderSpace(new double[proj.getInputDimensionality()]);
    // direction of the selected Point
    double[] selPoint = proj.fastProjectDataToRenderSpace(mid);
    double[] range1, range2;
    {
        // Rotation plane:
        double[] p1 = proj.fastProjectRenderToDataSpace(selPoint[0] + 10, selPoint[1]);
        double[] p2 = proj.fastProjectRenderToDataSpace(selPoint[0], selPoint[1] + 10);
        double[] pm = mid.toArray();
        // Compute relative vectors
        minusEquals(p1, pm);
        minusEquals(p2, pm);
        // Scale p1 and p2 to unit length:
        timesEquals(p1, 1. / euclideanLength(p1));
        timesEquals(p2, 1. / euclideanLength(p2));
        {
            double test = scalarProduct(p1, p2);
            if (Math.abs(test) > 1E-10) {
                LoggingUtil.warning("Projection does not seem to be orthogonal?");
            }
        }
        // Project onto p1, p2:
        double l1 = scalarProduct(pm, p1), l2 = scalarProduct(pm, p2);
        // Rotate projection by + and - angle
        // Using sin(-x) = -sin(x) and cos(-x)=cos(x)
        // To return cosine
        final DoubleWrapper tmp = new DoubleWrapper();
        final double sangle = FastMath.sinAndCos(angle, tmp), cangle = tmp.value;
        double r11 = +cangle * l1 - sangle * l2, r12 = +sangle * l1 + cangle * l2;
        double r21 = +cangle * l1 + sangle * l2, r22 = -sangle * l1 + cangle * l2;
        // Build rotated vectors - remove projected component, add rotated
        // component:
        double[] r1 = copy(pm), r2 = copy(pm);
        plusTimesEquals(r1, p1, -l1 + r11);
        plusTimesEquals(r1, p2, -l2 + r12);
        plusTimesEquals(r2, p1, -l1 + r21);
        plusTimesEquals(r2, p2, -l2 + r22);
        // Project to render space:
        range1 = proj.fastProjectDataToRenderSpace(r1);
        range2 = proj.fastProjectDataToRenderSpace(r2);
    }
    // Continue lines to viewport.
    {
        CanvasSize viewport = proj.estimateViewport();
        minusEquals(range1, pointOfOrigin);
        minusEquals(range2, pointOfOrigin);
        timesEquals(range1, viewport.continueToMargin(pointOfOrigin, range1));
        timesEquals(range2, viewport.continueToMargin(pointOfOrigin, range2));
        plusEquals(range1, pointOfOrigin);
        plusEquals(range2, pointOfOrigin);
        // Go backwards into the other direction - the origin might not be in the
        // viewport!
        double[] start1 = minus(pointOfOrigin, range1);
        double[] start2 = minus(pointOfOrigin, range2);
        timesEquals(start1, viewport.continueToMargin(range1, start1));
        timesEquals(start2, viewport.continueToMargin(range2, start2));
        plusEquals(start1, range1);
        plusEquals(start2, range2);
        // TODO: add filled variant?
        SVGPath path = new SVGPath();
        path.moveTo(start1);
        path.lineTo(range1);
        path.moveTo(start2);
        path.lineTo(range2);
        return path.makeElement(svgp);
    }
}

Example 2 with DoubleWrapper

use of net.jafama.DoubleWrapper in project elki by elki-project.

the class SphereUtil method bearingRad.

/**
 * Compute the bearing from start to end.
 *
 * @param latS Start latitude, in radians
 * @param lngS Start longitude, in radians
 * @param latE End latitude, in radians
 * @param lngE End longitude, in radians
 * @return Bearing in degree
 */
public static double bearingRad(double latS, double lngS, double latE, double lngE) {
    // To return cosine
    final DoubleWrapper tmp = new DoubleWrapper();
    final double slatS = sinAndCos(latS, tmp), clatS = tmp.value;
    final double slatE = sinAndCos(latE, tmp), clatE = tmp.value;
    return atan2(-sin(lngS - lngE) * clatE, clatS * slatE - slatS * clatE * cos(lngS - lngE));
}

Example 3 with DoubleWrapper

use of net.jafama.DoubleWrapper in project elki by elki-project.

the class SphereUtil method latlngMinDistRadFull.

/**
 * Point to rectangle minimum distance.
 *
 * Previous version, only around for reference.
 *
 * Complexity:
 * <ul>
 * <li>Trivial cases (on longitude slice): no trigonometric functions.</li>
 * <li>Cross-track case: 10+2 trig</li>
 * <li>Corner case: 10+3 trig, 1 sqrt</li>
 * </ul>
 *
 * Reference:
 * <p>
 * Erich Schubert, Arthur Zimek and Hans-Peter Kriegel<br />
 * Geodetic Distance Queries on R-Trees for Indexing Geographic Data<br />
 * 13th Int. Symposium on Advances in Spatial and Temporal Databases
 * </p>
 *
 * @param plat Latitude of query point.
 * @param plng Longitude of query point.
 * @param rminlat Min latitude of rectangle.
 * @param rminlng Min longitude of rectangle.
 * @param rmaxlat Max latitude of rectangle.
 * @param rmaxlng Max longitude of rectangle.
 * @return Distance in radians
 */
@// 
Reference(// 
authors = "Erich Schubert, Arthur Zimek and Hans-Peter Kriegel", // 
title = "Geodetic Distance Queries on R-Trees for Indexing Geographic Data", // 
booktitle = "13th Int. Symposium on Advances in Spatial and Temporal Databases", url = "http://dx.doi.org/10.1007/978-3-642-40235-7_9")
public static double latlngMinDistRadFull(double plat, double plng, double rminlat, double rminlng, double rmaxlat, double rmaxlng) {
    // Degenerate rectangles:
    if ((rminlat >= rmaxlat) && (rminlng >= rmaxlng)) {
        return haversineFormulaRad(rminlat, rminlng, plat, plng);
    }
    // The simplest case is when the query point is in the same "slice":
    if (rminlng <= plng && plng <= rmaxlng) {
        return // Inside
        (rminlat <= plat && plat <= rmaxlat) ? // Inside
        0 : // S, N
        (plat < rminlat) ? rminlat - plat : plat - rmaxlat;
    }
    // Determine whether going east or west is shorter.
    double lngE = rminlng - plng;
    lngE = lngE >= 0 ? lngE : lngE + TWOPI;
    // we keep this negative!
    double lngW = plng - rmaxlng;
    lngW = lngW >= 0 ? lngW : lngW + TWOPI;
    // Compute sine and cosine values we will certainly need below:
    // To return cosine
    final DoubleWrapper tmp = new DoubleWrapper();
    final double slatQ = sinAndCos(plat, tmp), clatQ = tmp.value;
    final double slatN = sinAndCos(rmaxlat, tmp), clatN = tmp.value;
    final double slatS = sinAndCos(rminlat, tmp), clatS = tmp.value;
    // Head east, to min edge:
    if (lngE <= lngW) {
        final double slngD = sinAndCos(lngE, tmp), clngD = tmp.value;
        // Bearing to south
        // atan2(slngD * clatS, clatQ * slatS - slatQ * clatS * clngD);
        // Bearing from south
        final double bs = atan2(slngD * clatQ, clatS * slatQ - slatS * clatQ * clngD);
        // Bearing to north
        // atan2(slngD * clatN, clatQ * slatN - slatQ * clatN * clngD);
        // Bearing from north
        final double bn = atan2(slngD * clatQ, clatN * slatQ - slatN * clatQ * clngD);
        if (bs < HALFPI && bn > HALFPI) {
            // Radians from south pole = abs(ATD)
            final double radFromS = -HALFPI - plat;
            // Cross-track-distance to longitude line.
            return asin(sin(radFromS) * -slngD);
        }
        if (bs - HALFPI < HALFPI - bn) {
            // Haversine to north corner.
            final double slatN2 = sin((plat - rmaxlat) * .5);
            final double slon = sin(lngE * .5);
            final double aN = slatN2 * slatN2 + slon * slon * clatQ * clatN;
            return 2 * atan2(sqrt(aN), sqrt(1 - aN));
        } else {
            // Haversine to south corner.
            final double slatS2 = sin((plat - rminlat) * .5);
            final double slon = sin(lngE * .5);
            final double aS = slatS2 * slatS2 + slon * slon * clatQ * clatS;
            return 2 * atan2(sqrt(aS), sqrt(1 - aS));
        }
    } else {
        // Head west, to max edge
        final double slngD = -sinAndCos(lngW, tmp), clngD = tmp.value;
        // Bearing to south
        // atan2(slngD * clatS, clatQ * slatS - slatQ * clatS * clngD);
        // Bearing from south
        final double bs = atan2(slngD * clatQ, clatS * slatQ - slatS * clatQ * clngD);
        // Bearing to north
        // atan2(slngD * clatN, clatQ * slatN - slatQ * clatN * clngD);
        // Bearing from north
        final double bn = atan2(slngD * clatQ, clatN * slatQ - slatN * clatQ * clngD);
        if (bs > -HALFPI && bn < -HALFPI) {
            // Radians from south = abs(ATD) = distance from pole
            final double radFromS = -HALFPI - plat;
            // Cross-track-distance to longitude line.
            return asin(sin(radFromS) * slngD);
        }
        if (-HALFPI - bs < bn + HALFPI) {
            // Haversine to north corner.
            final double slatN2 = sin((plat - rmaxlat) * .5);
            final double slon = sin(lngW * .5);
            final double aN = slatN2 * slatN2 + slon * slon * clatQ * clatN;
            return 2 * atan2(sqrt(aN), sqrt(1 - aN));
        } else {
            // Haversine to south corner.
            final double slatS2 = sin((plat - rminlat) * .5);
            final double slon = sin(lngW * .5);
            final double aS = slatS2 * slatS2 + slon * slon * clatQ * clatS;
            return 2 * atan2(sqrt(aS), sqrt(1 - aS));
        }
    }
}

Example 4 with DoubleWrapper

use of net.jafama.DoubleWrapper in project elki by elki-project.

the class SphereUtil method cosineOrHaversineRad.

/**
 * Use cosine or haversine dynamically.
 *
 * Complexity: 4-5 trigonometric functions, 1 sqrt.
 *
 * @param lat1 Latitude of first point in degree
 * @param lon1 Longitude of first point in degree
 * @param lat2 Latitude of second point in degree
 * @param lon2 Longitude of second point in degree
 * @return Distance on unit sphere
 */
public static double cosineOrHaversineRad(double lat1, double lon1, double lat2, double lon2) {
    if (lat1 == lat2 && lon1 == lon2) {
        return 0.;
    }
    // To return cosine
    final DoubleWrapper tmp = new DoubleWrapper();
    final double slat1 = sinAndCos(lat1, tmp), clat1 = tmp.value;
    final double slat2 = sinAndCos(lat2, tmp), clat2 = tmp.value;
    final double dlat = lat1 - lat2, dlon = lon1 - lon2;
    // Use cosine, cheaper:
    if (dlat > 0.01 || dlat < -0.01 || dlon > 0.01 || dlat < -0.01) {
        final double a = slat1 * slat2 + clat1 * clat2 * cos(dlon);
        return a < .9999_9999_9999_999 ? acos(a) : 0;
    }
    // Haversine formula, higher precision at < 1 meters
    final double slat = sin(dlat * .5), slon = sin(dlon * .5);
    return 2 * asin(sqrt(slat * slat + slon * slon * clat1 * clat2));
}

Example 5 with DoubleWrapper

use of net.jafama.DoubleWrapper in project elki by elki-project.

the class SphereUtil method crossTrackDistanceRad.

/**
 * Compute the cross-track distance.
 *
 * XTD = asin(sin(dist_SQ)*sin(crs_SQ-crs_SE))
 *
 * @param lat1 Latitude of starting point.
 * @param lon1 Longitude of starting point.
 * @param lat2 Latitude of destination point.
 * @param lon2 Longitude of destination point.
 * @param latQ Latitude of query point.
 * @param lonQ Longitude of query point.
 * @return Cross-track distance in km. May be negative - this gives the side.
 */
public static double crossTrackDistanceRad(double lat1, double lon1, double lat2, double lon2, double latQ, double lonQ) {
    final double dlon12 = lon2 - lon1;
    final double dlon1Q = lonQ - lon1;
    final double dlat1Q = latQ - lat1;
    // Compute trigonometric functions only once.
    // To return cosine
    final DoubleWrapper tmp = new DoubleWrapper();
    final double slat1 = sinAndCos(lat1, tmp), clat1 = tmp.value;
    final double slatQ = sinAndCos(latQ, tmp), clatQ = tmp.value;
    final double slat2 = sinAndCos(lat2, tmp), clat2 = tmp.value;
    // Compute the course
    // y = sin(dlon) * cos(lat2)
    final double sdlon12 = sinAndCos(dlon12, tmp), cdlon12 = tmp.value;
    final double sdlon1Q = sinAndCos(dlon1Q, tmp), cdlon1Q = tmp.value;
    final double yE = sdlon12 * clat2;
    final double yQ = sdlon1Q * clatQ;
    // x = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dlon)
    final double xE = clat1 * slat2 - slat1 * clat2 * cdlon12;
    final double xQ = clat1 * slatQ - slat1 * clatQ * cdlon1Q;
    // Calculate cross-track distance
    // Haversine formula, higher precision at < 1 meters but maybe issues at
    // antipodal points - we do not yet multiply with the radius!
    final double slat = sin(dlat1Q * .5);
    final double slon = sin(dlon1Q * .5);
    final double a = slat * slat + slon * slon * clat1 * clatQ;
    if (a > 0.9999_9999_9999_999 || a < -0.9999_9999_9999_999 || a == 0.) {
        return 0.;
    }
    final double crs12 = atan2(yE, xE);
    final double crs1Q = atan2(yQ, xQ);
    return asin(sqrt(a) * sqrt(1 - a) * 2 * sin(crs1Q - crs12));
// final double angDist1Q = a < 1 ? 2 * asin(sqrt(a)) : 0;
// return asin(sin(angDist1Q) * sin(crs1Q - crs12));
}