- performance optimization (without caching it's faster)

- more robust affine calculation (with linear solver)
- binary data loading instead of ASCII

build.gradle

@@ -1,10 +1,11 @@
 apply plugin: 'java-library'
 apply plugin: 'eclipse'
 
+sourceCompatibility = 1.5
+targetCompatibility = 1.5
+
 eclipse {
   jdt {
-    sourceCompatibility = 1.5
-    targetCompatibility = 1.5
     javaRuntimeName = "J2SE-1.5"
   }
 }
@@ -15,8 +16,8 @@ repositories {
 
 dependencies {
     api 'org.locationtech.jts:jts-core:1.15.0'
-	api 'org.slf4j:slf4j-api:1.7.25'
 
     testImplementation 'junit:junit:4.12'
 }
 
+

src/main/java/com/sinergise/geometry/crs/sitrans96/SiTrans96.java

@@ -20,6 +20,12 @@
  * 
  */
 public class SiTrans96 {
+	static final TriangleProvider TIN_D48GK;
+	static final TriangleProvider TIN_D96TM;
+	static {
+		TIN_D48GK = TriangularTransformationTinFromFileLoader.loadTin();
+		TIN_D96TM = TIN_D48GK.inverse();
+	}
 
 	public static class D48gkToD96tmTriangular implements SimpleTransform {
 		public Coordinate transformPoint(Coordinate pos) {
@@ -33,14 +39,6 @@ public Coordinate transformPoint(Coordinate pos) {
 		}
 	}
 
-	static final TriangleProvider TIN_D48GK;
-	static final TriangleProvider TIN_D96TM;
-	static {
-		TriangleFromTinProvider srcProvider = TriangularTransformationTinFromFileLoader.loadTin();
-		TIN_D48GK = new CachingTriangleProvider(srcProvider);
-		TIN_D96TM = new CachingTriangleProvider(srcProvider.inverse());
-	}
-
 	public static final SimpleTransform D48GK_TO_D96TM_TRIANGULAR = new D48gkToD96tmTriangular();
 	public static final SimpleTransform D96TM_TO_D48GK_TRIANGULAR = new D96tmToD48gkTriangular();
 
@@ -59,8 +57,4 @@ public static Coordinate d48gk_to_d96tm(Coordinate pos) {
 	public static Coordinate d96tm_to_d48gk(Coordinate pos) {
 		return D96TM_TO_D48GK_TRIANGULAR.transformPoint(pos);
 	}
-
-	public static void main(String[] args) {
-		System.out.println(d48gk_to_d96tm(350394, 142349));
-	}
 }

src/main/java/com/sinergise/geometry/crs/sitrans96/TriangleProvider.java

@@ -1,7 +1,73 @@
 package com.sinergise.geometry.crs.sitrans96;
 
+import java.util.Arrays;
+
 import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.Envelope;
+import org.locationtech.jts.triangulate.IncrementalDelaunayTriangulator;
+import org.locationtech.jts.triangulate.quadedge.QuadEdge;
+import org.locationtech.jts.triangulate.quadedge.QuadEdgeSubdivision;
+import org.locationtech.jts.triangulate.quadedge.TriangleVisitor;
+
+public class TriangleProvider {
+	QuadEdgeSubdivision delaunay;
+	TriangularTransformationPoint[] points;
+
+	public TriangleProvider(TriangularTransformationPoint[] pointData) {
+		this.points = pointData;
+		this.delaunay = buildTriangles(pointData);
+	}
+
+	public TriangleProvider inverse() {
+		TriangularTransformationPoint[] invPoints = new TriangularTransformationPoint[points.length];
+		for (int i = 0; i < invPoints.length; i++) {
+			invPoints[i] = points[i].reverse();
+		}
+		return new TriangleProvider(invPoints);
+	}
+
+	public TriangularTransformationTriangle getTriangleAt(Coordinate pos) {
+		QuadEdge e = delaunay.locate(pos);
+		if (e != null) {
+			TriangularTransformationTriangle data = (TriangularTransformationTriangle) e.getData();
+			if (data != null) {
+				return data;
+			}
+		}
+		throw new IllegalArgumentException("Coordinate out of bounds " + pos + " ENV: " + delaunay.getEnvelope());
+	}
+
+	private static final QuadEdgeSubdivision buildTriangles(TriangularTransformationPoint[] points) {
+		QuadEdgeSubdivision subdiv = new QuadEdgeSubdivision(aggregateEnvelope(points), 0);
+
+		IncrementalDelaunayTriangulator triangulator = new IncrementalDelaunayTriangulator(subdiv);
+		triangulator.insertSites(Arrays.asList(points));
+
+		buildTriangleObjects(subdiv);
+
+		return subdiv;
+	}
+
+	private static Envelope aggregateEnvelope(TriangularTransformationPoint[] points) {
+		Envelope e = new Envelope();
+		for (TriangularTransformationPoint tp : points) {
+			e.expandToInclude(tp.src());
+		}
+		return e;
+	}
 
-public interface TriangleProvider {
-	TriangularTransformationTriangle getTriangleAt(Coordinate pos);
+	private static void buildTriangleObjects(QuadEdgeSubdivision subdiv) {
+		TriangleVisitor triVisitor = new TriangleVisitor() {
+			public void visit(QuadEdge[] triEdges) {
+				TriangularTransformationTriangle tri = new TriangularTransformationTriangle(//
+						(TriangularTransformationPoint) triEdges[0].orig(),
+						(TriangularTransformationPoint) triEdges[1].orig(),
+						(TriangularTransformationPoint) triEdges[2].orig());
+				triEdges[0].setData(tri);
+				triEdges[1].setData(tri);
+				triEdges[2].setData(tri);
+			}
+		};
+		subdiv.visitTriangles(triVisitor, false);
+	}
 }

src/main/java/com/sinergise/geometry/crs/sitrans96/TriangularTransformationTinFromFileLoader.java

@@ -1,53 +1,42 @@
 package com.sinergise.geometry.crs.sitrans96;
 
-import static java.lang.Double.parseDouble;
-
 import java.io.IOException;
 import java.io.InputStream;
-import java.io.InputStreamReader;
-import java.io.LineNumberReader;
+import java.nio.ByteBuffer;
+import java.nio.DoubleBuffer;
+import java.nio.channels.Channels;
 import java.util.ArrayList;
 
 import org.locationtech.jts.geom.Coordinate;
 
 import com.sinergise.geometry.crs.util.Util;
 
 public class TriangularTransformationTinFromFileLoader {
-	private static final String FNAME_TRANS_POINTS = "/GK2TM_VVT4.csv";
+	private static final String FNAME_TRANS_POINTS_BIN = "/GK2TM_VVT4.bin";
 
-	public static TriangleFromTinProvider loadTin() {
-		return new TriangleFromTinProvider(loadPointDataFromFile());
+	public static TriangleProvider loadTin() {
+		return new TriangleProvider(loadPointDataFromFile());
 	}
 
 	private static TriangularTransformationPoint[] loadPointDataFromFile() {
 		ArrayList<TriangularTransformationPoint> ret = new ArrayList<TriangularTransformationPoint>(1000);
-
-		LineNumberReader rdr = null;
 		InputStream is = null;
-		try {
-			is = TriangularTransformationTinFromFileLoader.class.getResourceAsStream(FNAME_TRANS_POINTS);
-			rdr = new LineNumberReader(new InputStreamReader(is));
 
-			String curLine;
-			while ((curLine = rdr.readLine()) != null) {
-				ret.add(readPointFromLine(curLine));
+		try {
+			is = TriangularTransformationTinFromFileLoader.class.getResourceAsStream(FNAME_TRANS_POINTS_BIN);
+			ByteBuffer buf = ByteBuffer.allocate(is.available());
+			Channels.newChannel(is).read(buf);
+			buf.rewind();
+			DoubleBuffer dbuf = buf.asDoubleBuffer();
+			while (dbuf.hasRemaining()) {
+				ret.add(new TriangularTransformationPoint(new Coordinate(dbuf.get(), dbuf.get()),
+						new Coordinate(dbuf.get(), dbuf.get())));
 			}
-
 			return ret.toArray(new TriangularTransformationPoint[ret.size()]);
-		} catch(IOException e) {
+		} catch (IOException e) {
 			throw new RuntimeException(e);
 		} finally {
-			Util.closeSilent(rdr, is);
+			Util.closeSilent(is);
 		}
 	}
-
-	private static TriangularTransformationPoint readPointFromLine(String curLine) {
-		String[] parts = curLine.trim().split("\\s+");
-
-		double gkE = parseDouble(parts[3]);
-		double gkN = parseDouble(parts[4]);
-		double d96e = parseDouble(parts[1]);
-		double d96n = parseDouble(parts[2]);
-		return new TriangularTransformationPoint(new Coordinate(gkE, gkN), new Coordinate(d96e, d96n));
-	}
 }

src/main/java/com/sinergise/geometry/crs/sitrans96/TriangularTransformationTriangle.java

@@ -3,7 +3,6 @@
 import org.locationtech.jts.geom.Coordinate;
 
 import com.sinergise.geometry.crs.util.AffineTransform2D;
-import com.sinergise.geometry.crs.util.GeomUtil;
 import com.sinergise.geometry.crs.util.SimpleTransform;
 
 public class TriangularTransformationTriangle implements SimpleTransform {
@@ -14,7 +13,7 @@ public class TriangularTransformationTriangle implements SimpleTransform {
 	final AffineTransform2D affine;
 
 	public TriangularTransformationTriangle(TriangularTransformationPoint a, TriangularTransformationPoint b, TriangularTransformationPoint c) {
-		this(a, b, c, GeomUtil.triangleAffine(//
+		this(a, b, c, AffineTransform2D.triangleAffine(//
 				a.src(), a.tgt(), //
 				b.src(), b.tgt(), //
 				c.src(), c.tgt())
@@ -28,10 +27,6 @@ public TriangularTransformationTriangle(TriangularTransformationPoint a, Triangu
 		this.affine = affine;
 	}
 
-	public boolean contains(final Coordinate pos) {
-		return GeomUtil.triangleContains(a.src(), b.src(), c.src(), pos);
-	}
-
 	public Coordinate transformPoint(Coordinate input) {
 		return affine.transformPoint(input);
 	}

src/main/java/com/sinergise/geometry/crs/util/AffineTransform2D.java

@@ -1,6 +1,7 @@
 package com.sinergise.geometry.crs.util;
 
 import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.math.Matrix;
 
 public final class AffineTransform2D implements SimpleTransform {
 	final double m0;
@@ -30,4 +31,49 @@ public double x(final double x, final double y) {
 	public double y(final double x, final double y) {
 		return m1 * x + m3 * y + m5;
 	}
+
+	public static AffineTransform2D triangleAffine(Coordinate src1, Coordinate tgt1, Coordinate src2, Coordinate tgt2, Coordinate src3, Coordinate tgt3) {
+		double sx1 = src1.x;
+		double sy1 = src1.y;
+		double sx2 = src2.x;
+		double sy2 = src2.y;
+		double sx3 = src3.x;
+		double sy3 = src3.y;
+
+		double tx1 = tgt1.x;
+		double ty1 = tgt1.y;
+		double tx2 = tgt2.x;
+		double ty2 = tgt2.y;
+		double tx3 = tgt3.x;
+		double ty3 = tgt3.y;
+
+
+		double sumSx2 = sx1 * sx1 + sx2 * sx2 + sx3 * sx3;
+		double sumSy2 = sy1 * sy1 + sy2 * sy2 + sy3 * sy3;
+		double sumSxy = sx1 * sy1 + sx2 * sy2 + sx3 * sy3;
+		double sumSx = sx1 + sx2 + sx3;
+		double sumSy = sy1 + sy2 + sy3;
+		double[][] Ax = new double[][] {
+			{sumSx2, sumSxy, sumSx},
+			{sumSxy, sumSy2, sumSy},
+			{sumSx, sumSy, 3}
+		};
+		double[][] Ay = new double[3][3];
+		for (int i = 0; i < 3; i++) {
+			System.arraycopy(Ax[i], 0, Ay[i], 0, 3);
+		}
+		double[] Bxx = new double[] {
+			tx1 * sx1 + tx2 * sx2 + tx3 * sx3,
+			tx1 * sy1 + tx2 * sy2 + tx3 * sy3,
+			tx1 + tx2 + tx3
+		};
+		double[] Bxy = new double[] {
+			ty1 * sx1 + ty2 * sx2 + ty3 * sx3,
+			ty1 * sy1 + ty2 * sy2 + ty3 * sy3,
+			ty1 + ty2 + ty3
+		};
+		double[] mx = Matrix.solve(Ax, Bxx);
+		double[] my = Matrix.solve(Ay, Bxy);
+		return new AffineTransform2D(mx[0], my[0], mx[1], my[1], mx[2], my[2]);
+	}
 }