第18章 - 精度模型与数值鲁棒性

public class DoubleBasics {
    public static void main(String[] args) {
        System.out.println(0.1 + 0.2 == 0.3);   // false
        System.out.println(1e16 + 1.0 == 1e16); // true，1 被吞掉
        System.out.println(1e16 + 1.0 - 1e16);  // 0.0
        System.out.println((0.1 + 0.2) - 0.3);  // 5.55e-17
    }
}

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.PrecisionModel;

public class PrecisionModelDemo {
    public static void main(String[] args) {
        // 1) 浮点模型：不做任何舍入
        PrecisionModel floating = new PrecisionModel(PrecisionModel.FLOATING);
        System.out.println(floating.isFloating());            // true
        System.out.println(floating.gridSize());              // NaN
        System.out.println(floating.makePrecise(1.23456789)); // 1.23456789

        // 2) 固定模型：保留 3 位小数
        PrecisionModel pm1000 = new PrecisionModel(1000);
        System.out.println(pm1000.getScale());                // 1000.0
        System.out.println(pm1000.gridSize());                // 0.001
        System.out.println(pm1000.makePrecise(1.23456789));   // 1.235

        // 3) 负 scale：直接指定 1 米网格
        PrecisionModel grid1m = new PrecisionModel(-1);
        System.out.println(grid1m.getScale());                // 1.0
        System.out.println(grid1m.gridSize());                // 1.0
        System.out.println(grid1m.makePrecise(1.6));          // 2.0
        System.out.println(grid1m.makePrecise(1.4));          // 1.0

        // 4) 单精度模型
        PrecisionModel single = new PrecisionModel(PrecisionModel.FLOATING_SINGLE);
        System.out.println(single.makePrecise(1.23456789));   // 1.2345679

        // 5) 对 Coordinate 就地舍入
        Coordinate c = new Coordinate(1.23456789, 9.87654321);
        pm1000.makePrecise(c);
        System.out.println(c);                                // (1.235, 9.877, NaN)
    }
}

import org.locationtech.jts.geom.PrecisionModel;

public class RoundingDemo {
    public static void main(String[] args) {
        PrecisionModel pm = new PrecisionModel(-1); // 1 米网格
        System.out.println(pm.makePrecise(0.5));    // 1.0
        System.out.println(pm.makePrecise(1.5));    // 2.0
        System.out.println(pm.makePrecise(2.5));    // 3.0
        System.out.println(pm.makePrecise(-0.5));   // 0.0
        System.out.println(pm.makePrecise(-1.5));   // -1.0
        System.out.println(pm.makePrecise(-2.5));   // -2.0

        // NaN 与 Inf 原样返回
        System.out.println(pm.makePrecise(Double.NaN));               // NaN
        System.out.println(pm.makePrecise(Double.POSITIVE_INFINITY)); // Infinity
    }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.precision.GeometryPrecisionReducer;

public class PrecisionReduceDemo {
    public static void main(String[] args) throws Exception {
        WKTReader r = new WKTReader();
        Geometry g = r.read(
            "POLYGON ((0.123456 0.654321, 10.987654 0.111111,"
          + " 10.222222 10.333333, 0.444444 10.555555, 0.123456 0.654321))");

        PrecisionModel pm = new PrecisionModel(-0.001); // 0.001 米网格

        // 1) 拓扑约简：结果有效，坐标钉在网格上
        Geometry topo = GeometryPrecisionReducer.reduce(g, pm);
        System.out.println(topo.isValid());      // true
        System.out.println(topo.getNumPoints()); // 5（重复点被移除）
        System.out.println(topo);

        // 2) 逐点约简：只改坐标，不做拓扑修复
        Geometry pointwise = GeometryPrecisionReducer.reducePointwise(g, pm);
        System.out.println(pointwise.isValid()); // 本例输入凸多边形，仍为 true

        // 3) 改变精度模型元数据
        GeometryPrecisionReducer reducer = new GeometryPrecisionReducer(pm);
        reducer.setChangePrecisionModel(true);
        Geometry withPm = reducer.reduce(g);
        System.out.println(withPm.getPrecisionModel().getType()); // FIXED

        // 4) 输入无效时拓扑约简可能抛异常
        Geometry bad = r.read("POLYGON ((0 0, 10 10, 10 0, 0 10, 0 0))");
        try {
            GeometryPrecisionReducer.reduce(bad, new PrecisionModel(-0.1));
            System.out.println("未抛异常：无效性可能被精度约简消除");
        } catch (IllegalArgumentException ex) {
            System.out.println("约简失败: " + ex.getMessage());
        }
    }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.precision.GeometryPrecisionReducer;

public class CollapseDemo {
    public static void main(String[] args) throws Exception {
        WKTReader r = new WKTReader();
        Geometry thin = r.read("POLYGON ((0 0, 10 0, 10 0.0004, 0 0.0004, 0 0))");
        PrecisionModel pm = new PrecisionModel(-1); // 1 米网格

        Geometry removed = GeometryPrecisionReducer.reduce(thin, pm);
        System.out.println(removed.isEmpty());      // 通常为 true，塌缩部件被移除

        Geometry kept = GeometryPrecisionReducer.reduceKeepCollapsed(thin, pm);
        System.out.println(kept.isEmpty());         // 取决于输入与版本，可能为空
        System.out.println(kept.getGeometryType()); // 可能保留为低维/退化几何
        System.out.println(kept.isValid());         // 可能为 false（退化几何）
    }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.overlayng.PrecisionUtil;

public class PrecisionUtilDemo {
    public static void main(String[] args) throws Exception {
        System.out.println(PrecisionUtil.safeScale(135.0));     // 1.0E11
        System.out.println(PrecisionUtil.safeScale(3.85e7));    // 1.0E6
        System.out.println(PrecisionUtil.MAX_ROBUST_DP_DIGITS); // 14

        WKTReader r = new WKTReader();
        Geometry g = r.read("LINESTRING (38400000 2500000, 38410000 2501000)");
        double scale = PrecisionUtil.safeScale(g);
        PrecisionModel pm = new PrecisionModel(scale);
        System.out.println(pm.getScale());   // 1.0E6
        System.out.println(pm.gridSize());   // 1.0E-6

        // 慢但能保住全部已有精度（仅演示，生产慎用）
        double inherent = PrecisionUtil.inherentScale(g);
        System.out.println(inherent);        // 1.0（本例坐标只有整数）
    }
}

import java.util.ArrayList;
import java.util.List;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.overlayng.OverlayNG;
import org.locationtech.jts.operation.overlayng.OverlayNGRobust;
import org.locationtech.jts.operation.overlayng.PrecisionUtil;

public class OverlayPrecisionDemo {
    public static void main(String[] args) throws Exception {
        WKTReader r = new WKTReader();
        Geometry a = r.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");
        Geometry b = r.read("POLYGON ((5 -1, 15 -1, 15 5, 5 5, 5 -1))");

        // 1) 纯浮点叠加（最常用）
        Geometry g1 = OverlayNG.overlay(a, b, OverlayNG.INTERSECTION);
        System.out.println(g1.getArea());          // 25.0

        // 2) 固定精度叠加：结果吸附到 0.001 网格
        PrecisionModel pm = new PrecisionModel(-0.001);
        Geometry g2 = OverlayNG.overlay(a, b, OverlayNG.INTERSECTION, pm);
        System.out.println(g2.getArea());          // 25.0（本例恰好不变）

        // 3) 稳健叠加：生产推荐入口
        Geometry g3 = OverlayNGRobust.overlay(a, b, OverlayNG.INTERSECTION);
        System.out.println(g3.getArea());          // 25.0

        // 4) 稳健判定比例尺（自动）
        PrecisionModel auto = PrecisionUtil.robustPM(a, b);
        System.out.println(auto.getScale());

        // 5) 稳健一元合并
        List<Geometry> parts = new ArrayList<Geometry>();
        parts.add(a);
        parts.add(b);
        Geometry merged = OverlayNGRobust.union(parts);
        System.out.println(merged.getArea());      // 175.0
    }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.overlayng.OverlayNG;

public class StrictModeDemo {
    public static void main(String[] args) throws Exception {
        WKTReader r = new WKTReader();
        Geometry a = r.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");
        // 一条穿过面的线：面与面求交时，线不属于交集；
        // 但面与线求交的非严格结果会保留线
        Geometry line = r.read("LINESTRING (-5 5, 15 5)");

        Geometry normal = OverlayNG.overlay(a, line, OverlayNG.INTERSECTION);
        System.out.println(normal.getGeometryType()); // LineString（线在面内部的部分）

        OverlayNG ov = new OverlayNG(a, line, new PrecisionModel(-0.001),
                OverlayNG.INTERSECTION);
        ov.setStrictMode(true);
        Geometry strict = ov.getResult();
        System.out.println(strict.getGeometryType());
        // 输入维度不同时，strict 仍保留低维交集；同维时才会丢弃塌缩线
    }
}

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.geom.TopologyException;
import org.locationtech.jts.io.WKTWriter;
import org.locationtech.jts.operation.overlayng.OverlayNG;
import org.locationtech.jts.operation.overlayng.OverlayNGRobust;
import org.locationtech.jts.operation.overlayng.PrecisionUtil;

public class RobustOverlay {
    public static Geometry overlay(Geometry g0, Geometry g1, int opCode) {
        try {
            return OverlayNG.overlay(g0, g1, opCode);
        } catch (TopologyException ex) {
            Coordinate pt = ex.getCoordinate();
            System.out.println("浮点叠加失败，位置=" + pt + "，消息=" + ex.getMessage());
            // 第一级回退：固定精度 + snap-rounding
            try {
                PrecisionModel pm = new PrecisionModel(PrecisionUtil.safeScale(g0, g1));
                return OverlayNG.overlay(g0, g1, opCode, pm);
            } catch (RuntimeException ex2) {
                // 第二级回退：完整稳健链（吸附 + snap-rounding）
                try {
                    return OverlayNGRobust.overlay(g0, g1, opCode);
                } catch (TopologyException ex3) {
                    WKTWriter w = new WKTWriter();
                    System.out.println("g0 = " + w.write(g0));
                    System.out.println("g1 = " + w.write(g1));
                    throw ex3;
                }
            }
        }
    }
}

import org.locationtech.jts.algorithm.CGAlgorithmsDD;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.math.DD;

public class DDDemo {
    public static void main(String[] args) {
        // 1) 普通 double 会丢掉 1，DD 不会
        double d = 1e16;
        System.out.println(d + 1.0 - d);                       // 0.0
        DD exact = DD.valueOf(d).add(1.0).subtract(d);
        System.out.println(exact);                             // 1
        System.out.println(exact.doubleValue());               // 1.0

        // 2) 行列式符号
        DD det = DD.determinant(1.0, 2.0, 2.0, 4.0);
        System.out.println(det.signum());                      // 0（共线）

        // 3) 方向判断：使用 DD 版本
        Coordinate p1 = new Coordinate(0, 0);
        Coordinate p2 = new Coordinate(10, 10);
        Coordinate q = new Coordinate(5, 5.000000000000001);
        System.out.println(CGAlgorithmsDD.orientationIndex(p1, p2, q)); // 1

        // 4) 精确交点（两条线平行时返回 null）
        Coordinate ip = CGAlgorithmsDD.intersection(
            new Coordinate(0, 0), new Coordinate(10, 10),
            new Coordinate(0, 10), new Coordinate(10, 0));
        System.out.println(ip);                                // (5, 5, NaN)
    }
}

import org.locationtech.jts.algorithm.RobustLineIntersector;
import org.locationtech.jts.geom.Coordinate;

public class IntersectorDemo {
    public static void main(String[] args) {
        RobustLineIntersector li = new RobustLineIntersector();

        // 1) 正常相交
        li.computeIntersection(new Coordinate(0, 0), new Coordinate(10, 10),
                               new Coordinate(0, 10), new Coordinate(10, 0));
        System.out.println(li.hasIntersection());      // true
        System.out.println(li.getIntersectionNum());   // 1
        System.out.println(li.getIntersection(0));     // (5, 5, NaN)
        System.out.println(li.isProper());             // true

        // 2) 共线重叠：交点数为 2
        li.computeIntersection(new Coordinate(0, 0), new Coordinate(10, 0),
                               new Coordinate(5, 0), new Coordinate(15, 0));
        System.out.println(li.getIntersectionNum());   // 2
        System.out.println(li.isProper());             // false

        // 3) 端点接触
        li.computeIntersection(new Coordinate(0, 0), new Coordinate(10, 0),
                               new Coordinate(10, 0), new Coordinate(10, 10));
        System.out.println(li.getIntersectionNum());   // 1
        System.out.println(li.isInteriorIntersection()); // false
    }
}

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.precision.CommonBitsOp;
import org.locationtech.jts.precision.CommonBitsRemover;

public class CommonBitsDemo {
    public static void main(String[] args) throws Exception {
        WKTReader r = new WKTReader();
        Geometry a = r.read("POLYGON ((10000000.1 20000000.1,"
                          + " 10000010.1 20000000.1, 10000010.1 20000010.1,"
                          + " 10000000.1 20000010.1, 10000000.1 20000000.1))");
        Geometry b = r.read("POLYGON ((10000005.05 20000005.05,"
                          + " 10000015.05 20000005.05, 10000015.05 20000015.05,"
                          + " 10000005.05 20000015.05, 10000005.05 20000005.05))");

        // 1) 直接看公共位
        CommonBitsRemover remover = new CommonBitsRemover();
        remover.add(a);
        remover.add(b);
        Coordinate common = remover.getCommonCoordinate();
        System.out.println("公共基准 = " + common);

        // 2) 用 CommonBitsOp 做交集
        CommonBitsOp op = new CommonBitsOp(true); // true = 结果回移到原精度
        Geometry inter = op.intersection(a, b);
        System.out.println(inter.getArea());      // 约 25

        // 3) 回移后的结果不做舍入，可能无效，必须检查
        System.out.println(inter.isValid());      // true（本例简单，通常有效）
    }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.precision.EnhancedPrecisionOp;

public class EnhancedPrecisionDemo {
    public static void main(String[] args) throws Exception {
        WKTReader r = new WKTReader();
        Geometry a = r.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");
        Geometry b = r.read("POLYGON ((5 5, 15 5, 15 15, 5 15, 5 5))");

        // 与普通 intersection 同签名，内部在失败时用增强精度重试
        Geometry inter = EnhancedPrecisionOp.intersection(a, b);
        System.out.println(inter.getArea()); // 25.0
    }
}

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.operation.valid.IsValidOp;

public class InvalidOrdinateDemo {
    public static void main(String[] args) {
        Coordinate ok = new Coordinate(116.4, 39.9);
        Coordinate nan = new Coordinate(Double.NaN, 39.9);
        Coordinate inf = new Coordinate(116.4, Double.POSITIVE_INFINITY);
        Coordinate badZ = new Coordinate(116.4, 39.9, Double.NaN);

        System.out.println(IsValidOp.isValid(ok));   // true
        System.out.println(IsValidOp.isValid(nan));  // false
        System.out.println(IsValidOp.isValid(inf));  // false
        System.out.println(badZ.isValid());          // true：z 为 NaN 不影响判定
    }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;

public class WktInfDemo {
    public static void main(String[] args) throws Exception {
        WKTReader reader = new WKTReader();
        // JTS 1.21 起可解析；1.20 及以前会抛 ParseException
        Geometry g = reader.read("LINESTRING (0 0, Inf 10)");
        System.out.println(g.getNumPoints());  // 2
        System.out.println(g.isValid());       // false（错误码 10）
        System.out.println(g.getEnvelopeInternal());
        // Envelope 可能包含 Infinity，后续求交/索引行为不可依赖
    }
}

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.LineString;
import org.locationtech.jts.geom.util.GeometryFixer;

public class CleanInvalidDemo {
    public static void main(String[] args) {
        GeometryFactory gf = new GeometryFactory();
        LineString dirty = gf.createLineString(new Coordinate[] {
            new Coordinate(0, 0),
            new Coordinate(Double.NaN, 5),
            new Coordinate(10, 10)
        });
        System.out.println(dirty.getNumPoints());       // 3
        Geometry cleaned = GeometryFixer.fix(dirty);
        System.out.println(cleaned.getNumPoints());     // 2
        System.out.println(cleaned.isValid());          // true
    }
}

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.CoordinateSequence;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryCollection;

public final class OrdinateCheck {
    private OrdinateCheck() {}

    public static boolean allFinite(Geometry g) {
        if (g.isEmpty()) return true;
        if (g instanceof GeometryCollection) {
            for (int i = 0; i < g.getNumGeometries(); i++) {
                if (!allFinite(g.getGeometryN(i))) return false;
            }
            return true;
        }
        CoordinateSequence seq = g.getCoordinateSequence();
        for (int i = 0; i < seq.size(); i++) {
            double x = seq.getX(i);
            double y = seq.getY(i);
            if (Double.isNaN(x) || Double.isInfinite(x)) return false;
            if (Double.isNaN(y) || Double.isInfinite(y)) return false;
        }
        return true;
    }

    public static void main(String[] args) {
        System.out.println(allFinite(new Coordinate(1, 2)));
        // 注意：Coordinate 不是 Geometry，这里仅演示方法骨架
    }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.geom.util.AffineTransformation;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.overlayng.PrecisionUtil;
import org.locationtech.jts.precision.GeometryPrecisionReducer;

public class ToleranceDesignDemo {
    public static void main(String[] args) throws Exception {
        WKTReader r = new WKTReader();
        Geometry g = r.read("LINESTRING (0.0001 0.0002, 10.0003 10.0004)");

        // 经纬度数据：先算安全比例尺，避免“凭感觉”设网格
        double safe = PrecisionUtil.safeScale(g);
        System.out.println(safe);   // 1.0E11（度制下可取到 1e-11 度）

        // 业务需要的是 1e-7 度（约 1 厘米），用 precisionModel(-1e-7) 表达网格
        PrecisionModel pm = new PrecisionModel(-1.0e-7);
        Geometry snapped = GeometryPrecisionReducer.reducePointwise(g, pm);
        System.out.println(snapped);

        // 大坐标数据：先平移到局部原点，再舍入（减少高位占用）
        AffineTransformation shift =
            AffineTransformation.translationInstance(-38400000.0, -2500000.0);
        Geometry local = shift.transform(g);
        System.out.println(local);
    }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.PrecisionModel;
import org.locationtech.jts.geom.TopologyException;
import org.locationtech.jts.operation.overlayng.OverlayNG;
import org.locationtech.jts.operation.overlayng.OverlayNGRobust;
import org.locationtech.jts.operation.overlayng.PrecisionUtil;
import org.locationtech.jts.operation.valid.IsValidOp;

public final class RobustOps {
    private RobustOps() {}

    public enum Mode { FLOAT, FIXED, ROBUST }

    public static Geometry overlay(Geometry g0, Geometry g1, int opCode, Mode mode) {
        if (!IsValidOp.isValid(g0) || !IsValidOp.isValid(g1)) {
            throw new IllegalArgumentException("输入几何无效，请先修复");
        }
        switch (mode) {
            case FIXED:
                PrecisionModel pm = new PrecisionModel(PrecisionUtil.safeScale(g0, g1));
                return OverlayNG.overlay(g0, g1, opCode, pm);
            case ROBUST:
                return OverlayNGRobust.overlay(g0, g1, opCode);
            case FLOAT:
            default:
                return OverlayNG.overlay(g0, g1, opCode);
        }
    }

    public static boolean safeIsValid(Geometry g) {
        try {
            return g != null && !g.isEmpty() && IsValidOp.isValid(g);
        } catch (RuntimeException ex) {
            return false;
        }
    }
}