第08章 - 空间关系谓词与 DE-9IM

// 反例：手写点是否在矩形内
boolean inside = p.x >= minX && p.x <= maxX && p.y >= minY && p.y <= maxY;

import org.locationtech.jts.geom.IntersectionMatrix;

// 1. 空矩阵（所有元素初始为 F，即 -1）
IntersectionMatrix im = new IntersectionMatrix();

// 2. 从 9 字符串构造
IntersectionMatrix im2 = new IntersectionMatrix("212101212");

// 3. 从另一个矩阵拷贝
IntersectionMatrix im3 = new IntersectionMatrix(im2);

// 读取某个格子（行 0=内部 1=边界 2=外部）
int dim = im2.get(0, 0);        // II 的维度，本例是 2

// 写某个格子
im.set(0, 0, 2);                // 设置 II = 2
im.setAtLeast(0, 0, 1);         // 取较大值（1 与 2 取 2）

// 转成 9 字符串
String s = im2.toString();      // "212101212"

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

public class RelateMatrixTable {

  private static final WKTReader READER = new WKTReader(new GeometryFactory());

  static void show(String label, String a, String b) throws Exception {
    Geometry ga = READER.read(a);
    Geometry gb = READER.read(b);
    IntersectionMatrix im = ga.relate(gb);
    System.out.printf("%-24s IM=%s  intersects=%s contains=%s covers=%s touches=%s crosses=%s overlaps=%s%n",
        label, im,
        ga.intersects(gb), ga.contains(gb), ga.covers(gb),
        ga.touches(gb), ga.crosses(gb), ga.overlaps(gb));
  }

  public static void main(String[] args) throws Exception {
    String A = "POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))";
    show("面-面 重叠", A, "POLYGON ((5 5, 15 5, 15 15, 5 15, 5 5))");
    show("面-面 共边", A, "POLYGON ((10 0, 20 0, 20 10, 10 10, 10 0))");
    show("面-面 共角", A, "POLYGON ((10 10, 20 10, 20 20, 10 20, 10 10))");
    show("面-面 内含", A, "POLYGON ((2 2, 4 2, 4 4, 2 4, 2 2))");
    show("面-面 相同", A, A);
    show("面-面 相离", A, "POLYGON ((20 20, 21 20, 21 21, 20 21, 20 20))");
    show("线-面 穿过", "LINESTRING (0 5, 20 5)", A);
    show("线-面 内部", "LINESTRING (2 2, 8 8)", A);
    show("线-面 沿边", "LINESTRING (0 0, 10 0)", A);
    show("线-线 交叉", "LINESTRING (0 0, 10 10)", "LINESTRING (0 10, 10 0)");
    show("线-线 端点接", "LINESTRING (0 0, 5 0)", "LINESTRING (5 0, 10 0)");
    show("线-线 重叠", "LINESTRING (0 0, 10 0)", "LINESTRING (5 0, 15 0)");
    show("点-面 边界上", "POINT (0 0)", A);
    show("点-面 内部", "POINT (5 5)", A);
  }
}

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

public class PatternUseCases {

  private static final WKTReader R = new WKTReader(new GeometryFactory());

  /** 打印矩阵，并测试一组候选模式。 */
  static void check(String label, String a, String b, String... patterns) throws Exception {
    Geometry ga = R.read(a);
    Geometry gb = R.read(b);
    IntersectionMatrix im = ga.relate(gb);
    StringBuilder sb = new StringBuilder("IM=" + im + "  ");
    for (String p : patterns) {
      sb.append(p).append("=").append(im.matches(p)).append("  ");
    }
    System.out.printf("%-22s %s%n", label, sb);
  }

  public static void main(String[] args) throws Exception {
    String A = "POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))";
    String edge = "POLYGON ((10 0, 20 0, 20 10, 10 10, 10 0))";      // 共边
    String corner = "POLYGON ((10 10, 20 10, 20 20, 10 20, 10 10))"; // 共角

    // 1. 面与面相邻：共边 or 共角
    check("AA 共边", A, edge, "F***T****", "FF*F1****", "FT*******");
    check("AA 共角", A, corner, "F***T****", "FF*F0****", "FT*******");

    // 2. 包含 vs 覆盖（点在边界上）
    check("点落在角上", A, "POINT (10 10)", "T*****FF*", "F**T*****");
    System.out.println("covers  = " + R.read(A).covers(R.read("POINT (10 10)")));
    System.out.println("contains= " + R.read(A).contains(R.read("POINT (10 10)")));

    // 3. 线沿面边界（典型“贴边”数据）
    check("线沿边界", A, "LINESTRING (0 0, 10 0)", "F**T*****", "T*****FF*");

    // 4. 严格包含（不碰边界）
    check("严格包含", A, "POLYGON ((2 2, 4 2, 4 4, 2 4, 2 2))", "T**FF*FF*", "T*****FF*");
  }
}

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

WKTReader r = new WKTReader(new GeometryFactory());

// 自相交的蝴蝶结：无效，但“简单性”判定要看具体类型
Geometry bow = r.read("POLYGON ((0 0, 10 10, 10 0, 0 10, 0 0))");
System.out.println("bow isValid  = " + bow.isValid());    // false
System.out.println("bow isSimple = " + bow.isSimple());   // 面的 isSimple 定义与 isValid 不同

// 自相交的线：无效？不，线没有无效概念，但“不简单”
Geometry zigzag = r.read("LINESTRING (0 0, 10 10, 10 0, 0 10)");
System.out.println("zigzag isValid  = " + zigzag.isValid());   // true
System.out.println("zigzag isSimple = " + zigzag.isSimple());  // false

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

WKTReader r = new WKTReader(new GeometryFactory());

Geometry bow = r.read("POLYGON ((0 0, 10 10, 10 0, 0 10, 0 0))");  // 自相交
Geometry sq  = r.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");

System.out.println("intersects = " + bow.intersects(sq));
System.out.println("contains   = " + bow.contains(sq));
System.out.println("relate     = " + bow.relate(sq));

// modules/core/src/main/java/org/locationtech/jts/geom/GeometryRelate.java
public static String RELATE_PROPERTY_NAME = "jts.relate";
public static String RELATE_PROPERTY_VALUE_NG = "ng";
public static String RELATE_PROPERTY_VALUE_OLD = "old";
public static boolean RELATE_NG_DEFAULT = false;   // 默认还是旧引擎

static {
  setRelateImpl(System.getProperty(RELATE_PROPERTY_NAME));
}

static boolean intersects(Geometry a, Geometry b) {
  if (isRelateNG) {
    return RelateNG.relate(a, b, RelatePredicate.intersects());
  }
  // 旧引擎：GeometryCollection 手工展开成笛卡尔积
  if (a.isGeometryCollection() || b.isGeometryCollection()) {
    for (int i = 0 ; i < a.getNumGeometries() ; i++) {
      for (int j = 0 ; j < b.getNumGeometries() ; j++) {
        if (a.getGeometryN(i).intersects(b.getGeometryN(j))) {
          return true;
        }
      }
    }
    return false;
  }
  return RelateOp.relate(a, b).isIntersects();
}

static IntersectionMatrix relate(Geometry a, Geometry b) {
  if (isRelateNG) {
    return RelateNG.relate(a, b);
  }
  Geometry.checkNotGeometryCollection(a);
  Geometry.checkNotGeometryCollection(b);
  return RelateOp.relate(a, b);
}

Geometry poly = r.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");
Geometry gc = r.read("GEOMETRYCOLLECTION (POINT (5 5), POINT (20 20))");

System.out.println("intersects(gc) = " + poly.intersects(gc));   // true

try {
  System.out.println(poly.relate(gc));
} catch (Exception e) {
  System.out.println("relate(gc) 异常: " + e.getMessage());
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.IntersectionMatrix;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.relateng.RelateNG;
import org.locationtech.jts.operation.relateng.RelatePredicate;

public class RelateNGDemo {

  public static void main(String[] args) throws Exception {
    WKTReader r = new WKTReader(new GeometryFactory());
    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))");

    // 1. 用谓词
    System.out.println("intersects = "
        + RelateNG.relate(a, b, RelatePredicate.intersects()));
    System.out.println("contains   = "
        + RelateNG.relate(a, b, RelatePredicate.contains()));
    System.out.println("touches    = "
        + RelateNG.relate(a, b, RelatePredicate.touches()));

    // 2. 用 DE-9IM 模式
    System.out.println("pattern T*T***T** = "
        + RelateNG.relate(a, b, "T*T***T**"));

    // 3. 取矩阵
    IntersectionMatrix im = RelateNG.relate(a, b);
    System.out.println("IM = " + im + "  isOverlaps = " + im.isOverlaps(2, 2));
  }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.relateng.RelateNG;
import org.locationtech.jts.operation.relateng.RelatePredicate;

public class PreparedRelateDemo {

  public static void main(String[] args) throws Exception {
    WKTReader r = new WKTReader(new GeometryFactory());
    Geometry parcel = r.read(
        "POLYGON ((0 0, 100 0, 100 100, 0 100, 0 0))");   // 一个大多边形

    // 预编译：内部会建空间索引
    RelateNG prepared = RelateNG.prepare(parcel);

    // 对大量小几何重复判定
    int inside = 0, touching = 0, outside = 0;
    for (int i = 0; i < 1000; i++) {
      Geometry pt = r.read("POINT (" + (i % 120) + " " + (i % 120) + ")");
      if (prepared.evaluate(pt, RelatePredicate.contains())) {
        inside++;
      }
      else if (prepared.evaluate(pt, RelatePredicate.covers())) {
        touching++;
      }
      else {
        outside++;
      }
    }
    System.out.println("inside=" + inside + " touching=" + touching + " outside=" + outside);
  }
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.prep.PreparedGeometry;
import org.locationtech.jts.geom.prep.PreparedGeometryFactory;
import org.locationtech.jts.io.WKTReader;

public class PreparedBench {

  public static void main(String[] args) throws Exception {
    WKTReader r = new WKTReader(new GeometryFactory());

    Geometry region = r.read(
        "POLYGON ((0 0, 1000 0, 1000 1000, 0 1000, 0 0))");

    Geometry[] points = new Geometry[20000];
    for (int i = 0; i < points.length; i++) {
      points[i] = r.read("POINT (" + (i % 1100) + " " + (i % 1100) + ")");
    }

    // 普通方式
    long t1 = System.nanoTime();
    int c1 = 0;
    for (Geometry g : points) {
      if (region.contains(g)) c1++;
    }
    long t2 = System.nanoTime();

    // Prepared 方式
    PreparedGeometry pg = PreparedGeometryFactory.prepare(region);
    long t3 = System.nanoTime();
    int c2 = 0;
    for (Geometry g : points) {
      if (pg.contains(g)) c2++;
    }
    long t4 = System.nanoTime();

    System.out.printf("普通: %d 次命中, %.1f ms%n", c1, (t2 - t1) / 1e6);
    System.out.printf("Prepared: %d 次命中, %.1f ms%n", c2, (t4 - t3) / 1e6);
    System.out.println("结果一致 = " + (c1 == c2));
  }
}

// 用 SineStarFactory 造一个复杂多边形
Geometry star = SineStarFactory.create(new Coordinate(0, 0), 200, 400, 20, 0.5);
Geometry complex = star.buffer(2);        // 621 个顶点

PreparedGeometry pg = PreparedGeometryFactory.prepare(complex);
// 5 万个点、2 万条线分别用普通方式与 Prepared 方式判定

// 会加速
pg.intersects(b);
pg.contains(b);
pg.covers(b);

// 不一定加速（可能回退）
pg.touches(b);
pg.crosses(b);

import org.locationtech.jts.algorithm.BoundaryNodeRule;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.BoundaryOp;
import org.locationtech.jts.operation.relateng.RelateNG;
import org.locationtech.jts.operation.relateng.RelatePredicate;

public class BoundaryRuleDemo {

  public static void main(String[] args) throws Exception {
    WKTReader r = new WKTReader(new GeometryFactory());

    // 两段首尾相接的线，中间点 (1 0) 被两条线共享
    Geometry mls = r.read("MULTILINESTRING ((0 0, 1 0), (1 0, 2 0))");
    Geometry mid = r.read("POINT (1 0)");
    Geometry end = r.read("POINT (0 0)");

    System.out.println("默认边界       = " + mls.getBoundary());
    System.out.println("ENDPOINT 边界  = "
        + BoundaryOp.getBoundary(mls, BoundaryNodeRule.ENDPOINT_BOUNDARY_RULE));

    for (BoundaryNodeRule rule : new BoundaryNodeRule[] {
        BoundaryNodeRule.MOD2_BOUNDARY_RULE,
        BoundaryNodeRule.ENDPOINT_BOUNDARY_RULE,
        BoundaryNodeRule.MULTIVALENT_ENDPOINT_BOUNDARY_RULE,
        BoundaryNodeRule.MONOVALENT_ENDPOINT_BOUNDARY_RULE }) {

      System.out.printf("%-42s mid.touches=%s end.touches=%s%n", rule,
          RelateNG.relate(mls, mid, RelatePredicate.touches(), rule),
          RelateNG.relate(mls, end, RelatePredicate.touches(), rule));
    }
  }
}

Geometry ring = r.read("LINESTRING (0 0, 10 0, 10 10, 0 10, 0 0)");
System.out.println("闭合线边界 = " + ring.getBoundary());
System.out.println("闭合线 touches 起点 = " + ring.touches(r.read("POINT (0 0)")));

Geometry poly = r.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");
Geometry empty = r.read("POINT EMPTY");

System.out.println("intersects = " + poly.intersects(empty));
System.out.println("disjoint   = " + poly.disjoint(empty));
System.out.println("contains   = " + poly.contains(empty));
System.out.println("covers     = " + poly.covers(empty));
System.out.println("relate     = " + RelateNG.relate(poly, empty));

public boolean safeContains(Geometry a, Geometry b) {
  if (a == null || b == null || a.isEmpty() || b.isEmpty()) {
    return false;   // 按业务约定：空几何一律不满足包含
  }
  return a.contains(b);
}

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.IntersectionMatrix;
import org.locationtech.jts.io.WKTReader;
import org.locationtech.jts.operation.relateng.RelateNG;
import org.locationtech.jts.operation.relateng.RelatePredicate;

WKTReader r = new WKTReader(new GeometryFactory());
Geometry poly = r.read("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))");

// 一个在内、一个在外的点集
Geometry gc = r.read("GEOMETRYCOLLECTION (POINT (5 5), POINT (20 20))");
// 两个点都在内
Geometry gcIn = r.read("GEOMETRYCOLLECTION (POINT (5 5), POINT (6 6))");

try {
  poly.relate(gc);
} catch (IllegalArgumentException e) {
  System.out.println("旧引擎 relate(GC): " + e.getMessage());
}

System.out.println("旧引擎 intersects(GC) = " + poly.intersects(gc));

IntersectionMatrix im = RelateNG.relate(poly, gc);
System.out.println("RelateNG IM(GC)      = " + im);
System.out.println("RelateNG intersects  = "
    + RelateNG.relate(poly, gc, RelatePredicate.intersects()));
System.out.println("RelateNG contains    = "
    + RelateNG.relate(poly, gc, RelatePredicate.contains()));

System.out.println("全在内的 GC IM       = " + RelateNG.relate(poly, gcIn));
System.out.println("全在内的 GC contains = "
    + RelateNG.relate(poly, gcIn, RelatePredicate.contains()));