diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..f0bce6b --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,3 @@ +[build-system] +requires = ["setuptools>=61", "wheel", "Cython>=3.0"] +build-backend = "setuptools.build_meta" diff --git a/setup.py b/setup.py index 8c9cebe..4b49a5d 100644 --- a/setup.py +++ b/setup.py @@ -1,6 +1,8 @@ from setuptools import setup from setuptools.extension import Extension +from Cython.Build import cythonize + setup( name = "tps", version = "0.2", @@ -15,13 +17,21 @@ "Operating System :: OS Independent", "Programming Language :: C", "Programming Language :: C++", - "Programming Language :: Python :: 2.5", - "Programming Language :: Python :: 2.6", - "Programming Language :: Python :: 2.7", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.12", "Topic :: Scientific/Engineering", ], test_suite = 'tps.test.test_suite', - ext_modules=[ - Extension("tps._tps", ["tps/_tps.cpp", "tps/thinplatespline.cpp"], libraries=["stdc++"]), - ], + python_requires=">=3.8", + ext_modules=cythonize( + [ + Extension( + "tps._tps", + ["tps/_tps.pyx", "tps/thinplatespline.cpp"], + libraries=["stdc++"], + language="c++", + ), + ], + language_level=3, + ), ) diff --git a/tps/__init__.py b/tps/__init__.py index cf14886..b9fd013 100644 --- a/tps/__init__.py +++ b/tps/__init__.py @@ -18,6 +18,11 @@ # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER # DEALINGS IN THE SOFTWARE. -from _tps import TPS, TPSError, from_control_points +from importlib.util import find_spec -__all__ = ['TPS', 'TPSError', 'from_control_points'] \ No newline at end of file +if find_spec("tps._tps") is not None: + from ._tps import TPS, TPSError, from_control_points +else: + from ._tps_py import TPS, TPSError, from_control_points + +__all__ = ['TPS', 'TPSError', 'from_control_points'] diff --git a/tps/_tps_py.py b/tps/_tps_py.py new file mode 100644 index 0000000..579ce7f --- /dev/null +++ b/tps/_tps_py.py @@ -0,0 +1,169 @@ +# Copyright (c) 2011, Omniscale GmbH & Co. KG +# +# Permission is hereby granted, free of charge, to any person obtaining a +# copy of this software and associated documentation files (the "Software"), +# to deal in the Software without restriction, including without limitation +# the rights to use, copy, modify, merge, publish, distribute, sublicense, +# and/or sell copies of the Software, and to permit persons to whom the +# Software is furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +# DEALINGS IN THE SOFTWARE. + +from __future__ import annotations + +import math +from typing import Iterable, List, Sequence, Tuple + +Point = Tuple[float, float, float, float] + + +class TPSError(Exception): + pass + + +def _kernel(r: float) -> float: + if r == 0.0: + return 0.0 + r2 = r * r + return r2 * math.log(r2) + + +def _solve_linear_system(matrix: List[List[float]], vector: List[float]) -> List[float]: + size = len(matrix) + augmented = [row[:] + [vector[index]] for index, row in enumerate(matrix)] + + for col in range(size): + pivot_row = max(range(col, size), key=lambda r: abs(augmented[r][col])) + pivot_value = augmented[pivot_row][col] + if abs(pivot_value) < 1e-12: + raise TPSError("could not solve thin plate spline") + if pivot_row != col: + augmented[col], augmented[pivot_row] = augmented[pivot_row], augmented[col] + + pivot_value = augmented[col][col] + for j in range(col, size + 1): + augmented[col][j] /= pivot_value + + for row in range(size): + if row == col: + continue + factor = augmented[row][col] + if factor == 0.0: + continue + for j in range(col, size + 1): + augmented[row][j] -= factor * augmented[col][j] + + return [augmented[row][size] for row in range(size)] + + +class TPS: + """ + Thin Plate Spline computation class. + """ + + def __init__(self, points: Iterable[Point] | None = None) -> None: + self._points: List[Point] = [] + self._solved = False + self._weights_x: List[float] = [] + self._weights_y: List[float] = [] + self._affine_x: List[float] = [] + self._affine_y: List[float] = [] + if points: + for point in points: + self.add(*point) + + def add(self, src_x: float, src_y: float, dst_x: float, dst_y: float) -> None: + self._points.append((src_x, src_y, dst_x, dst_y)) + self._solved = False + + def solve(self) -> None: + if not self._points: + raise TPSError("could not solve thin plate spline") + + count = len(self._points) + size = count + 3 + matrix = [[0.0 for _ in range(size)] for _ in range(size)] + vector_x = [0.0 for _ in range(size)] + vector_y = [0.0 for _ in range(size)] + + for i, (src_x_i, src_y_i, dst_x_i, dst_y_i) in enumerate(self._points): + vector_x[i] = dst_x_i + vector_y[i] = dst_y_i + matrix[i][count] = 1.0 + matrix[i][count + 1] = src_x_i + matrix[i][count + 2] = src_y_i + matrix[count][i] = 1.0 + matrix[count + 1][i] = src_x_i + matrix[count + 2][i] = src_y_i + for j, (src_x_j, src_y_j, _, _) in enumerate(self._points): + dx = src_x_i - src_x_j + dy = src_y_i - src_y_j + matrix[i][j] = _kernel(math.hypot(dx, dy)) + matrix[i][i] += 1e-12 + + params_x = _solve_linear_system(matrix, vector_x) + params_y = _solve_linear_system(matrix, vector_y) + + self._weights_x = params_x[:count] + self._weights_y = params_y[:count] + self._affine_x = params_x[count:] + self._affine_y = params_y[count:] + self._solved = True + + def transform(self, src_x: float, src_y: float) -> Tuple[float, float]: + if not self._points: + raise TPSError("could not solve thin plate spline") + + if len(self._points) == 1: + _, _, dst_x, dst_y = self._points[0] + return float(dst_x), float(dst_y) + + if len(self._points) == 2: + (src_x_1, src_y_1, dst_x_1, dst_y_1), (src_x_2, src_y_2, dst_x_2, dst_y_2) = ( + self._points + ) + dx = src_x_2 - src_x_1 + dy = src_y_2 - src_y_1 + denom = dx * dx + dy * dy + if denom == 0.0: + return float(dst_x_1), float(dst_y_1) + t = ((src_x - src_x_1) * dx + (src_y - src_y_1) * dy) / denom + value_x = dst_x_1 + t * (dst_x_2 - dst_x_1) + value_y = dst_y_1 + t * (dst_y_2 - dst_y_1) + return round(value_x, 12), round(value_y, 12) + + if not self._solved: + self.solve() + + value_x = self._affine_x[0] + self._affine_x[1] * src_x + self._affine_x[2] * src_y + value_y = self._affine_y[0] + self._affine_y[1] * src_x + self._affine_y[2] * src_y + + for (point_x, point_y, _, _), weight_x, weight_y in zip( + self._points, + self._weights_x, + self._weights_y, + ): + value = _kernel(math.hypot(src_x - point_x, src_y - point_y)) + value_x += weight_x * value + value_y += weight_y * value + + return round(value_x, 12), round(value_y, 12) + + +def from_control_points(points: Sequence[Point], backwards: bool = False) -> TPS: + tps = TPS() + for point in points: + if backwards: + tps.add(point[2], point[3], point[0], point[1]) + else: + tps.add(*point) + return tps diff --git a/tps/test/__init__.py b/tps/test/__init__.py index 5a76bb9..96530bc 100644 --- a/tps/test/__init__.py +++ b/tps/test/__init__.py @@ -1,9 +1,8 @@ from unittest import TestSuite -import test_tps +from . import test_tps def test_suite(): suite = TestSuite() suite.addTest(test_tps.test_suite()) return suite - diff --git a/tps/test/test_tps.py b/tps/test/test_tps.py index fe89d84..64572f7 100644 --- a/tps/test/test_tps.py +++ b/tps/test/test_tps.py @@ -25,17 +25,17 @@ class TestTPS(unittest.TestCase): def test_init_from_list(self): t = TPS([(0, 0, 50, 50), (10, 10, 100, 100)]) - self.assertEquals(t.transform(4, 5), (72.5, 72.5)) + self.assertEqual(t.transform(4, 5), (72.5, 72.5)) t.add(0, 10, 70, 100) - self.assertEquals(t.transform(4, 5), (72.0, 75.0)) + self.assertEqual(t.transform(4, 5), (72.0, 75.0)) def test_simple(self): t = TPS() t.add(0, 0, 50, 50) t.add(10, 10, 100, 100) - self.assertEquals(t.transform(4, 5), (72.5, 72.5)) + self.assertEqual(t.transform(4, 5), (72.5, 72.5)) t.add(0, 10, 70, 100) - self.assertEquals(t.transform(4, 5), (72.0, 75.0)) + self.assertEqual(t.transform(4, 5), (72.0, 75.0)) def test_no_points(self): try: @@ -52,7 +52,7 @@ def test_from_control_points_list(self): (10, 10, 100, 100), (0, 10, 70, 100)]) - self.assertEquals(t.transform(4, 5), (72.0, 75.0)) + self.assertEqual(t.transform(4, 5), (72.0, 75.0)) def test_from_control_points_list_backwards(self): t = from_control_points([ @@ -62,11 +62,11 @@ def test_from_control_points_list_backwards(self): backwards=True) results = t.transform(72, 75) - self.assertAlmostEquals(results[0], 4.0) - self.assertAlmostEquals(results[1], 5.0) + self.assertAlmostEqual(results[0], 4.0) + self.assertAlmostEqual(results[1], 5.0) def test_suite(): return unittest.TestLoader().loadTestsFromTestCase(TestTPS) if __name__ == '__main__': - unittest.main() \ No newline at end of file + unittest.main()