commit 36d7f037c3dbefe7a20282057d1b8a580d5892e8 Author: svxf Date: Mon Apr 1 16:02:59 2019 +0400 init diff --git a/Pipfile b/Pipfile new file mode 100644 index 0000000..f583c09 --- /dev/null +++ b/Pipfile @@ -0,0 +1,14 @@ +[[source]] +url = "https://pypi.org/simple" +verify_ssl = true +name = "pypi" + +[packages] +numpy = "*" +pytest = "*" +matplotlib = "*" + +[dev-packages] + +[requires] +python_version = "3.7" diff --git a/Pipfile.lock b/Pipfile.lock new file mode 100644 index 0000000..0c70eb0 --- /dev/null +++ b/Pipfile.lock @@ -0,0 +1,174 @@ +{ + "_meta": { + "hash": { + "sha256": "38985b134d3f8d2a44b9c4d413de1240fd302ad89956c4c8265719c3445ee91b" + }, + "pipfile-spec": 6, + "requires": { + "python_version": "3.7" + }, + "sources": [ + { + "name": "pypi", + "url": "https://pypi.org/simple", + "verify_ssl": true + } + ] + }, + "default": { + "atomicwrites": { + "hashes": [ + "sha256:03472c30eb2c5d1ba9227e4c2ca66ab8287fbfbbda3888aa93dc2e28fc6811b4", + "sha256:75a9445bac02d8d058d5e1fe689654ba5a6556a1dfd8ce6ec55a0ed79866cfa6" + ], + "version": "==1.3.0" + }, + "attrs": { + "hashes": [ + "sha256:69c0dbf2ed392de1cb5ec704444b08a5ef81680a61cb899dc08127123af36a79", + "sha256:f0b870f674851ecbfbbbd364d6b5cbdff9dcedbc7f3f5e18a6891057f21fe399" + ], + "version": "==19.1.0" + }, + "cycler": { + "hashes": [ + "sha256:1d8a5ae1ff6c5cf9b93e8811e581232ad8920aeec647c37316ceac982b08cb2d", + "sha256:cd7b2d1018258d7247a71425e9f26463dfb444d411c39569972f4ce586b0c9d8" + ], + "version": "==0.10.0" + }, + "kiwisolver": { + "hashes": [ + "sha256:0ee4ed8b3ae8f5f712b0aa9ebd2858b5b232f1b9a96b0943dceb34df2a223bc3", + "sha256:0f7f532f3c94e99545a29f4c3f05637f4d2713e7fd91b4dd8abfc18340b86cd5", + "sha256:1a078f5dd7e99317098f0e0d490257fd0349d79363e8c923d5bb76428f318421", + "sha256:1aa0b55a0eb1bd3fa82e704f44fb8f16e26702af1a073cc5030eea399e617b56", + "sha256:2874060b91e131ceeff00574b7c2140749c9355817a4ed498e82a4ffa308ecbc", + "sha256:379d97783ba8d2934d52221c833407f20ca287b36d949b4bba6c75274bcf6363", + "sha256:3b791ddf2aefc56382aadc26ea5b352e86a2921e4e85c31c1f770f527eb06ce4", + "sha256:4329008a167fac233e398e8a600d1b91539dc33c5a3eadee84c0d4b04d4494fa", + "sha256:45813e0873bbb679334a161b28cb9606d9665e70561fd6caa8863e279b5e464b", + "sha256:53a5b27e6b5717bdc0125338a822605084054c80f382051fb945d2c0e6899a20", + "sha256:574f24b9805cb1c72d02b9f7749aa0cc0b81aa82571be5201aa1453190390ae5", + "sha256:66f82819ff47fa67a11540da96966fb9245504b7f496034f534b81cacf333861", + "sha256:79e5fe3ccd5144ae80777e12973027bd2f4f5e3ae8eb286cabe787bed9780138", + "sha256:83410258eb886f3456714eea4d4304db3a1fc8624623fc3f38a487ab36c0f653", + "sha256:8b6a7b596ce1d2a6d93c3562f1178ebd3b7bb445b3b0dd33b09f9255e312a965", + "sha256:9576cb63897fbfa69df60f994082c3f4b8e6adb49cccb60efb2a80a208e6f996", + "sha256:95a25d9f3449046ecbe9065be8f8380c03c56081bc5d41fe0fb964aaa30b2195", + "sha256:a424f048bebc4476620e77f3e4d1f282920cef9bc376ba16d0b8fe97eec87cde", + "sha256:aaec1cfd94f4f3e9a25e144d5b0ed1eb8a9596ec36d7318a504d813412563a85", + "sha256:acb673eecbae089ea3be3dcf75bfe45fc8d4dcdc951e27d8691887963cf421c7", + "sha256:b15bc8d2c2848a4a7c04f76c9b3dc3561e95d4dabc6b4f24bfabe5fd81a0b14f", + "sha256:b1c240d565e977d80c0083404c01e4d59c5772c977fae2c483f100567f50847b", + "sha256:c595693de998461bcd49b8d20568c8870b3209b8ea323b2a7b0ea86d85864694", + "sha256:ce3be5d520b4d2c3e5eeb4cd2ef62b9b9ab8ac6b6fedbaa0e39cdb6f50644278", + "sha256:e0f910f84b35c36a3513b96d816e6442ae138862257ae18a0019d2fc67b041dc", + "sha256:ea36e19ac0a483eea239320aef0bd40702404ff8c7e42179a2d9d36c5afcb55c", + "sha256:efabbcd4f406b532206b8801058c8bab9e79645b9880329253ae3322b7b02cd5", + "sha256:f923406e6b32c86309261b8195e24e18b6a8801df0cfc7814ac44017bfcb3939" + ], + "version": "==1.0.1" + }, + "matplotlib": { + "hashes": [ + "sha256:1ae6549976b6ceb6ee426272a28c0fc9715b3e3669694d560c8f661c5b39e2c5", + "sha256:4d4250bf508dd07cca3b43888097f873cadb66eec6ac63dbbfb798798ec07af2", + "sha256:53af2e01d7f1700ed2b64a9091bc865360c9c4032f625451c4589a826854c787", + "sha256:63e498067d32d627111cd1162cae1621f1221f9d4c6a9745dd7233f29de581b6", + "sha256:7169a34971e398dd58e87e173f97366fd88a3fa80852704530433eb224a8ca57", + "sha256:91c54d6bb9eeaaff965656c5ea6cbdcbf780bad8462ac99b30b451548194746f", + "sha256:aeef177647bb3fccfe09065481989d7dfc5ac59e9367d6a00a3481062cf651e4", + "sha256:cf8ae10559a78aee0409ede1e9d4fda03895433eeafe609dd9ed67e45f552db0", + "sha256:d51d0889d1c4d51c51a9822265c0494ea3e70a52bdd88358e0863daca46fa23a", + "sha256:de5ccd3500247f85fe4f9fad90f80a8bd397e4f110a4c33fabf95f07403e8372", + "sha256:e1d33589e32f482d0a7d1957bf473d43341115d40d33f578dad44432e47df7b7", + "sha256:e8d1939262aa6b36d0c51f50a50a43a04b9618d20db31e6c0192b1463067aeef", + "sha256:e918d51b1fda82a65fdf52d2f3914b2246481cc2a9cd10e223e6be6078916ff3" + ], + "index": "pypi", + "version": "==3.0.3" + }, + "more-itertools": { + "hashes": [ + "sha256:0125e8f60e9e031347105eb1682cef932f5e97d7b9a1a28d9bf00c22a5daef40", + "sha256:590044e3942351a1bdb1de960b739ff4ce277960f2425ad4509446dbace8d9d1" + ], + "markers": "python_version > '2.7'", + "version": "==6.0.0" + }, + "numpy": { + "hashes": [ + "sha256:1980f8d84548d74921685f68096911585fee393975f53797614b34d4f409b6da", + "sha256:22752cd809272671b273bb86df0f505f505a12368a3a5fc0aa811c7ece4dfd5c", + "sha256:23cc40313036cffd5d1873ef3ce2e949bdee0646c5d6f375bf7ee4f368db2511", + "sha256:2b0b118ff547fecabc247a2668f48f48b3b1f7d63676ebc5be7352a5fd9e85a5", + "sha256:3a0bd1edf64f6a911427b608a894111f9fcdb25284f724016f34a84c9a3a6ea9", + "sha256:3f25f6c7b0d000017e5ac55977a3999b0b1a74491eacb3c1aa716f0e01f6dcd1", + "sha256:4061c79ac2230594a7419151028e808239450e676c39e58302ad296232e3c2e8", + "sha256:560ceaa24f971ab37dede7ba030fc5d8fa173305d94365f814d9523ffd5d5916", + "sha256:62be044cd58da2a947b7e7b2252a10b42920df9520fc3d39f5c4c70d5460b8ba", + "sha256:6c692e3879dde0b67a9dc78f9bfb6f61c666b4562fd8619632d7043fb5b691b0", + "sha256:6f65e37b5a331df950ef6ff03bd4136b3c0bbcf44d4b8e99135d68a537711b5a", + "sha256:7a78cc4ddb253a55971115f8320a7ce28fd23a065fc33166d601f51760eecfa9", + "sha256:80a41edf64a3626e729a62df7dd278474fc1726836552b67a8c6396fd7e86760", + "sha256:893f4d75255f25a7b8516feb5766c6b63c54780323b9bd4bc51cdd7efc943c73", + "sha256:972ea92f9c1b54cc1c1a3d8508e326c0114aaf0f34996772a30f3f52b73b942f", + "sha256:9f1d4865436f794accdabadc57a8395bd3faa755449b4f65b88b7df65ae05f89", + "sha256:9f4cd7832b35e736b739be03b55875706c8c3e5fe334a06210f1a61e5c2c8ca5", + "sha256:adab43bf657488300d3aeeb8030d7f024fcc86e3a9b8848741ea2ea903e56610", + "sha256:bd2834d496ba9b1bdda3a6cf3de4dc0d4a0e7be306335940402ec95132ad063d", + "sha256:d20c0360940f30003a23c0adae2fe50a0a04f3e48dc05c298493b51fd6280197", + "sha256:d3b3ed87061d2314ff3659bb73896e622252da52558f2380f12c421fbdee3d89", + "sha256:dc235bf29a406dfda5790d01b998a1c01d7d37f449128c0b1b7d1c89a84fae8b", + "sha256:fb3c83554f39f48f3fa3123b9c24aecf681b1c289f9334f8215c1d3c8e2f6e5b" + ], + "index": "pypi", + "version": "==1.16.2" + }, + "pluggy": { + "hashes": [ + "sha256:19ecf9ce9db2fce065a7a0586e07cfb4ac8614fe96edf628a264b1c70116cf8f", + "sha256:84d306a647cc805219916e62aab89caa97a33a1dd8c342e87a37f91073cd4746" + ], + "version": "==0.9.0" + }, + "py": { + "hashes": [ + "sha256:64f65755aee5b381cea27766a3a147c3f15b9b6b9ac88676de66ba2ae36793fa", + "sha256:dc639b046a6e2cff5bbe40194ad65936d6ba360b52b3c3fe1d08a82dd50b5e53" + ], + "version": "==1.8.0" + }, + "pyparsing": { + "hashes": [ + "sha256:66c9268862641abcac4a96ba74506e594c884e3f57690a696d21ad8210ed667a", + "sha256:f6c5ef0d7480ad048c054c37632c67fca55299990fff127850181659eea33fc3" + ], + "version": "==2.3.1" + }, + "pytest": { + "hashes": [ + "sha256:592eaa2c33fae68c7d75aacf042efc9f77b27c08a6224a4f59beab8d9a420523", + "sha256:ad3ad5c450284819ecde191a654c09b0ec72257a2c711b9633d677c71c9850c4" + ], + "index": "pypi", + "version": "==4.3.1" + }, + "python-dateutil": { + "hashes": [ + "sha256:7e6584c74aeed623791615e26efd690f29817a27c73085b78e4bad02493df2fb", + "sha256:c89805f6f4d64db21ed966fda138f8a5ed7a4fdbc1a8ee329ce1b74e3c74da9e" + ], + "version": "==2.8.0" + }, + "six": { + "hashes": [ + "sha256:3350809f0555b11f552448330d0b52d5f24c91a322ea4a15ef22629740f3761c", + "sha256:d16a0141ec1a18405cd4ce8b4613101da75da0e9a7aec5bdd4fa804d0e0eba73" + ], + "version": "==1.12.0" + } + }, + "develop": {} +} diff --git a/movie/frame000.png b/movie/frame000.png new file mode 100644 index 0000000..0515f06 Binary files /dev/null and b/movie/frame000.png differ diff --git a/movie/frame001.png b/movie/frame001.png new file mode 100644 index 0000000..a60bc6d Binary files /dev/null and b/movie/frame001.png differ diff --git a/movie/frame002.png b/movie/frame002.png new file mode 100644 index 0000000..5995768 Binary files /dev/null and b/movie/frame002.png differ diff --git a/movie/frame003.png b/movie/frame003.png new file mode 100644 index 0000000..28bd820 Binary files /dev/null and b/movie/frame003.png differ diff --git a/movie/frame004.png b/movie/frame004.png new file mode 100644 index 0000000..d3412ad Binary files /dev/null and b/movie/frame004.png differ diff --git a/movie/frame005.png b/movie/frame005.png new file mode 100644 index 0000000..6e6b426 Binary files /dev/null and b/movie/frame005.png differ diff --git a/movie/frame006.png b/movie/frame006.png new file mode 100644 index 0000000..bcf82e6 Binary files /dev/null and b/movie/frame006.png differ diff --git a/movie/frame007.png b/movie/frame007.png new file mode 100644 index 0000000..a746674 Binary files /dev/null and b/movie/frame007.png differ diff --git a/movie/frame008.png b/movie/frame008.png new file mode 100644 index 0000000..4a58f78 Binary files /dev/null and b/movie/frame008.png differ diff --git a/movie/frame009.png b/movie/frame009.png new file mode 100644 index 0000000..77a5186 Binary files /dev/null and b/movie/frame009.png differ diff --git a/movie/frame010.png b/movie/frame010.png new file mode 100644 index 0000000..4af5d5f Binary files /dev/null and b/movie/frame010.png differ diff --git a/movie/frame011.png b/movie/frame011.png new file mode 100644 index 0000000..73b8913 Binary files /dev/null and b/movie/frame011.png differ diff --git a/movie/frame012.png b/movie/frame012.png new file mode 100644 index 0000000..5514206 Binary files /dev/null and b/movie/frame012.png differ diff --git a/movie/frame013.png b/movie/frame013.png new file mode 100644 index 0000000..5cee0f3 Binary files /dev/null and b/movie/frame013.png differ diff --git a/movie/frame014.png b/movie/frame014.png new file mode 100644 index 0000000..2c295ae Binary files /dev/null and b/movie/frame014.png differ diff --git a/movie/frame015.png b/movie/frame015.png new file mode 100644 index 0000000..0515f06 Binary files /dev/null and b/movie/frame015.png differ diff --git a/src/__init__.py b/src/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/geometry/__init__.py b/src/geometry/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/geometry/affine.py b/src/geometry/affine.py new file mode 100644 index 0000000..42b8188 --- /dev/null +++ b/src/geometry/affine.py @@ -0,0 +1,48 @@ +from dataclasses import dataclass + +import numpy as np + +from geometry.vector import Vector + + +@dataclass +class AffineTransform: + A: np.ndarray + A_extended: np.ndarray + + inv_A: np.ndarray + inv_A_extended: np.ndarray + + def __init__(self, A: np.ndarray, b: Vector): + n = len(A) + + assert A.shape == (n, n) + assert b.a.shape == (n, 1) + + self.A = np.array(A, dtype=float) + self.A.setflags(write=False) + + self.A_extended = np.zeros((n + 1, n + 1), dtype=float) + self.A_extended[:-1, :-1] = self.A + self.A_extended[:-1, -1:] = b.a + self.A_extended[-1, -1] = 1 + self.A_extended.setflags(write=False) + + self.inv_A = np.linalg.inv(self.A) + self.inv_A_extended = np.linalg.inv(self.A_extended) + + def apply_shift(self, v: Vector) -> Vector: + v_extended = np.vstack((v.a, np.eye(1))) + result_extended = self.A_extended @ v_extended + return Vector(result_extended[:-1, :]) + + def apply_no_shift(self, v: Vector) -> Vector: + return Vector(self.A @ v.a) + + def invert_shift(self, v: Vector) -> Vector: + v_extended = np.vstack((v.a, np.eye(1))) + result_extended = self.inv_A_extended @ v_extended + return Vector(result_extended[:-1, :]) + + def invert_no_shift(self, v: Vector) -> Vector: + return Vector(self.inv_A @ v.a) diff --git a/src/geometry/ellipse.py b/src/geometry/ellipse.py new file mode 100644 index 0000000..f939b8e --- /dev/null +++ b/src/geometry/ellipse.py @@ -0,0 +1,167 @@ +from dataclasses import dataclass +from typing import Optional + +import numpy as np + +from geometry.ray import Ray +from geometry.surface import Surface +from geometry.line_drawers import draw_line_3d +from geometry.vector import Vector +from geometry.rotations import rotate2d, rotate3d, rotate_rodrigues +from util import double_equality, minpos, solve_quadratic + + +@dataclass +class Ellipse(Surface): + refraction_index: float + f1: Vector + f2: Vector + distance: float + + def __post_init__(self): + if abs(self.f1 - self.f2) > self.distance: + raise ValueError() + + def distance_to(self, v): + return abs(self.f1 - v) + abs(self.f2 - v) + + def __and__(self, v: Vector) -> bool: + return double_equality(self.distance_to(v), self.distance) + + def __contains__(self, v: Vector) -> bool: + return self & v or self.distance_to(v) < self.distance + + def normal_at(self, v: Vector) -> Vector: + d1 = (v - self.f1).normalize() + d2 = (v - self.f2).normalize() + return (d1 + d2).normalize() + + def intersect(self, ray: Ray) -> Optional[float]: + D = self.distance + b1 = 2 * (ray.e @ (ray.r0 - self.f1)) + b2 = 2 * (ray.e @ (ray.r0 - self.f2)) + c1 = (ray.r0 - self.f1) @ (ray.r0 - self.f1) + c2 = (ray.r0 - self.f2) @ (ray.r0 - self.f2) + + a = 4 * D ** 2 - (b1 - b2) ** 2 + b = 2 * D ** 2 * (b1 + b2) - 2 * b1 * c1 + 2 * b1 * c2 + 2 * b2 * c1 - 2 * b2 * c2 + c = 2 * D ** 2 * (c1 + c2) - (c1 - c2) ** 2 - D ** 4 + + return minpos(solve_quadratic(a, b, c)) + + def draw2d(self, plt, **kwargs): + r0 = (self.f1 + self.f2) / 2 + e = Vector.fromargs(1, 0) + + xs = [] + ys = [] + + for alpha in np.linspace(0, 2 * np.pi, 200): + r = Ray(r0=r0, e=rotate2d(e, alpha)) + p = r.prolong(self.intersect(r)) + xs.append(p.a[0, 0]) + ys.append(p.a[1, 0]) + + # n = self.normal_at(p) + # p2 = p + n + # plt.plot( + # [p.a[0, 0], p2.a[0, 0]], + # [p.a[1, 0], p2.a[1, 0]], + # ) + + plt.plot(xs, ys, **kwargs) + + def draw3d_old(self, axis, **kwargs): + r0 = (self.f1 + self.f2) / 2 + e = Vector.fromargs(0, 0, 1) + + N1 = 20 + N2 = 100 + + for alpha in np.linspace(0, 2 * np.pi, N1): + xs = [] + ys = [] + zs = [] + + for beta in np.linspace(0, 2 * np.pi, N2): + d = Vector.fromargs( + np.sin(alpha) * np.cos(beta), + np.sin(alpha) * np.sin(beta), + np.cos(alpha), + ) + + r = Ray(r0=r0, e=d) + p = r.prolong(self.intersect(r)) + + xs.append(p.a[0, 0]) + ys.append(p.a[1, 0]) + zs.append(p.a[2, 0]) + + axis.plot(xs, ys, zs, color='blue', **kwargs) + + xs = [] + ys = [] + zs = [] + + for beta in np.linspace(0, 2 * np.pi, N2): + d = Vector.fromargs( + np.sin(beta) * np.cos(alpha), + np.sin(beta) * np.sin(alpha), + np.cos(beta), + ) + + r = Ray(r0=r0, e=d) + p = r.prolong(self.intersect(r)) + + xs.append(p.a[0, 0]) + ys.append(p.a[1, 0]) + zs.append(p.a[2, 0]) + + axis.plot(xs, ys, zs, color='blue', **kwargs) + + def draw3d(self, axis, **kwargs): + r0 = (self.f1 + self.f2) / 2 + a1 = (self.f2 - r0).normalize() + + a2_candidates = [ + a1 | Vector.fromargs(1, 0, 0), + a1 | Vector.fromargs(0, 1, 0), + a1 | Vector.fromargs(0, 0, 1), + ] + a2 = a2_candidates[np.argmax(np.abs(a2_candidates))] + a2 = a2.normalize() + a3 = (a1 | a2).normalize() + + draw_line_3d(axis, r0, r0 + a1, lw=5, zorder=0, color='red') + draw_line_3d(axis, r0, r0 + a2, lw=5, zorder=0, color='green') + draw_line_3d(axis, r0, r0 + a3, lw=5, zorder=0, color='blue') + + N1 = 8 + N2 = 100 + + for lat in np.linspace(0, np.pi, N1): + vs = [] + e = rotate_rodrigues(a1, a2, lat) + + for long in np.linspace(0, 2 * np.pi, N2): + r = Ray(r0=r0, e=rotate_rodrigues(e, a1, long)) + p = r.prolong(self.intersect(r)) + + vs.append(p.a) + + V = np.hstack(vs) + axis.plot(V[0, :], V[1, :], V[2, :], color='navy', **kwargs) + + for long in np.linspace(0, np.pi, N1): + vs = [] + e = rotate_rodrigues(a2, a1, long) + a3_e = rotate_rodrigues(a3, a1, long) + + for lat in np.linspace(0, 2 * np.pi, N2): + r = Ray(r0=r0, e=rotate_rodrigues(e, a3_e, lat)) + p = r.prolong(self.intersect(r)) + + vs.append(p.a) + + V = np.hstack(vs) + axis.plot(V[0, :], V[1, :], V[2, :], color='indigo', **kwargs) diff --git a/src/geometry/line_drawers.py b/src/geometry/line_drawers.py new file mode 100644 index 0000000..03e93e3 --- /dev/null +++ b/src/geometry/line_drawers.py @@ -0,0 +1,55 @@ +from dataclasses import dataclass +from typing import List + +import numpy as np +from matplotlib import pyplot as plt + +from geometry.affine import AffineTransform +from geometry.vector import Vector + + +def draw_line_2d(axis, v1: Vector, v2: Vector, **kwargs): + axis.plot( + [v1.a[0, 0], v2.a[0, 0]], + [v1.a[1, 0], v2.a[1, 0]], + **kwargs) + + +def draw_line_3d(axis, v1: Vector, v2: Vector, **kwargs): + axis.plot( + [v1.a[0, 0], v2.a[0, 0]], + [v1.a[1, 0], v2.a[1, 0]], + [v1.a[2, 0], v2.a[2, 0]], + **kwargs) + + +class Painter: + axis: plt.Axes + + def __init__(self, axis): + self.axis = axis + + def plot3d(self, vs: List[Vector], **kwargs): + V = np.hstack([v.a for v in vs]) + + self.axis.plot(V[0, :], V[1, :], V[2, :], **kwargs) + + def line3d(self, a: Vector, b: Vector, **kwargs): + self.plot3d([a, b], **kwargs) + + +class TransformedPainter: + painter: Painter + transform: AffineTransform + + def __init__(self, painter, transform): + self.painter = painter + self.transform = transform + + def plot3d(self, vs: List[Vector], **kwargs): + vs_ = list(map(self.transform.invert_shift, vs)) + + self.painter.plot3d(vs_, **kwargs) + + def line3d(self, a: Vector, b: Vector, **kwargs): + self.plot3d([a, b], **kwargs) diff --git a/src/geometry/plane.py b/src/geometry/plane.py new file mode 100644 index 0000000..b5f78cb --- /dev/null +++ b/src/geometry/plane.py @@ -0,0 +1,35 @@ +from dataclasses import dataclass +from typing import Optional + +from ray import Ray +from surface import Surface +from util import minpos +from vector import Vector + + +@dataclass +class Plane(Surface): + refraction_index: float + r0: Vector + normal: Vector + + def __post_init__(self): + self.normal = self.normal.normalize() + + def __and__(self, v: Vector) -> bool: + return self.normal % (self.r0 - v) + + def __contains__(self, v: Vector) -> bool: + proj = (self.r0 - v) @ self.normal + return self & v or proj > 0 + + def normal_at(self, r: Vector) -> Vector: + return self.normal + + def intersect(self, ray: Ray) -> Optional[float]: + if self.normal % ray.e: + return None + + t = (self.normal @ (self.r0 - ray.r0)) / (self.normal @ ray.e) + + return minpos([t]) diff --git a/src/geometry/ray.py b/src/geometry/ray.py new file mode 100644 index 0000000..ec1a704 --- /dev/null +++ b/src/geometry/ray.py @@ -0,0 +1,15 @@ +from dataclasses import dataclass + +from geometry.vector import Vector + + +@dataclass +class Ray: + r0: Vector + e: Vector + + def __post_init__(self): + self.e = self.e.normalize() + + def prolong(self, t: float) -> Vector: + return self.r0 + self.e * t diff --git a/src/geometry/reflections.py b/src/geometry/reflections.py new file mode 100644 index 0000000..08e3c46 --- /dev/null +++ b/src/geometry/reflections.py @@ -0,0 +1,20 @@ +import numpy as np + +from geometry.vector import Vector + + +def reflect(direction: Vector, normal: Vector) -> Vector: + reflected = direction - 2 * (direction @ normal) * normal + return reflected.normalize() + + +def refract(direction: Vector, normal: Vector, n1: float, n2: float) -> Vector: + if n1 > n2: + return refract(direction, -1 * normal, n2, n1) + + proj = normal @ (n1 * direction) + k = (n2 ** 2 - n1 ** 2) / proj ** 2 + 1 + refracted = (1 / n2) * (n1 * direction + normal * proj + * (np.sqrt(k) - 1)) + + return refracted.normalize() diff --git a/src/geometry/rotations.py b/src/geometry/rotations.py new file mode 100644 index 0000000..8a4fe8c --- /dev/null +++ b/src/geometry/rotations.py @@ -0,0 +1,45 @@ +import numpy as np + +from geometry.vector import Vector + + +def rotate2d(v: Vector, alpha: float) -> Vector: + R = np.array([ + [np.cos(alpha), -np.sin(alpha)], + [np.sin(alpha), np.cos(alpha)], + ]) + + return Vector(R @ v.a) + + +def rotate3d(v: Vector, alpha: float, axis: int) -> Vector: + R = np.eye(3) + + if axis == 0: + R[1, 1] = np.cos(alpha) + R[1, 2] = -np.sin(alpha) + R[2, 1] = np.sin(alpha) + R[2, 2] = np.cos(alpha) + + elif axis == 1: + R[0, 0] = np.cos(alpha) + R[0, 2] = np.sin(alpha) + R[2, 0] = -np.sin(alpha) + R[2, 2] = np.cos(alpha) + + elif axis == 2: + R[0, 0] = np.cos(alpha) + R[0, 1] = -np.sin(alpha) + R[1, 0] = np.sin(alpha) + R[1, 1] = np.cos(alpha) + + else: + raise ValueError + + return Vector(R @ v.a) + + +def rotate_rodrigues(v: Vector, k: Vector, alpha: float) -> Vector: + k = k.normalize() + + return v * np.cos(alpha) + (k | v) * np.sin(alpha) + k * (k @ v) * (1 - np.cos(alpha)) diff --git a/src/geometry/sphere.py b/src/geometry/sphere.py new file mode 100644 index 0000000..87e36c3 --- /dev/null +++ b/src/geometry/sphere.py @@ -0,0 +1,66 @@ +from dataclasses import dataclass +from typing import Optional + +import numpy as np + +from geometry.line_drawers import Painter +from geometry.ray import Ray +from geometry.rotations import rotate_rodrigues +from geometry.surface import Surface +from util import double_equality, minpos, solve_quadratic +from geometry.vector import Vector + + +@dataclass +class Sphere(Surface): + refraction_index: float + r0: Vector + radius: float + + def __and__(self, v: Vector) -> bool: + return double_equality(abs(self.r0 - v), self.radius) + + def __contains__(self, v: Vector) -> bool: + return self & v or abs(self.r0 - v) < self.radius + + def normal_at(self, v: Vector) -> Vector: + return (v - self.r0).normalize() + + def intersect(self, ray: Ray) -> Optional[float]: + b = -2 * ((self.r0 - ray.r0) @ ray.e) + c = (self.r0 - ray.r0) @ (self.r0 - ray.r0) - self.radius ** 2 + + return minpos(solve_quadratic(1, b, c)) + + def draw3d(self, painter: Painter): + a1 = Vector.fromargs(0, 0, 1) + a2 = Vector.fromargs(1, 0, 0) + a3 = Vector.fromargs(0, 1, 0) + + painter.line3d(self.r0, self.r0 + a1, lw=5, zorder=0, color='red') + painter.line3d(self.r0, self.r0 + a2, lw=5, zorder=0, color='green') + painter.line3d(self.r0, self.r0 + a3, lw=5, zorder=0, color='blue') + + N1 = 8 + N2 = 100 + + for lat in np.linspace(0, np.pi, N1): + vs = [] + e = rotate_rodrigues(a1, a2, lat) + + for long in np.linspace(0, 2 * np.pi, N2): + r = Ray(r0=self.r0, e=rotate_rodrigues(e, a1, long)) + vs.append(r.prolong(self.intersect(r))) + + painter.plot3d(vs, color='navy') + + for long in np.linspace(0, np.pi, N1): + vs = [] + e = rotate_rodrigues(a2, a1, long) + a3_e = rotate_rodrigues(a3, a1, long) + + for lat in np.linspace(0, 2 * np.pi, N2): + r = Ray(r0=self.r0, e=rotate_rodrigues(e, a3_e, lat)) + vs.append(r.prolong(self.intersect(r))) + + painter.plot3d(vs, color='indigo') diff --git a/src/geometry/surface.py b/src/geometry/surface.py new file mode 100644 index 0000000..6cb00b3 --- /dev/null +++ b/src/geometry/surface.py @@ -0,0 +1,31 @@ +from dataclasses import dataclass +from typing import Optional + +from geometry.ray import Ray +from geometry.vector import Vector + + +@dataclass +class Surface: + refraction_index: float + + def __and__(self, v: Vector) -> bool: + """ + Does this point lie on the surface itself? + """ + raise NotImplementedError() + + def __contains__(self, v: Vector) -> bool: + """ + Is this point _inside_ the region of space, separated by this surface? + """ + raise NotImplementedError() + + def normal_at(self, v: Vector) -> Vector: + raise NotImplementedError() + + def intersect(self, ray: Ray) -> Optional[float]: + raise NotImplementedError() + + def draw2d(self, plt, **kwargs): + raise NotImplementedError() diff --git a/src/geometry/transformed_surface.py b/src/geometry/transformed_surface.py new file mode 100644 index 0000000..7ab13eb --- /dev/null +++ b/src/geometry/transformed_surface.py @@ -0,0 +1,53 @@ +from dataclasses import dataclass +from typing import Optional + +import numpy as np + +from geometry.affine import AffineTransform +from geometry.ray import Ray +from geometry.surface import Surface +from geometry.line_drawers import draw_line_3d, Painter, TransformedPainter +from geometry.vector import Vector +from geometry.rotations import rotate2d, rotate3d, rotate_rodrigues +from util import double_equality, minpos, solve_quadratic + + +@dataclass +class TransformedSurface(Surface): + refraction_index: float + subsurface: Surface + transform: AffineTransform + + def __contains__(self, v: Vector) -> bool: + return self.transform.apply_shift(v) in self.subsurface + + def __and__(self, v: Vector) -> bool: + return self.subsurface & self.transform.apply_shift(v) + + def normal_at(self, v: Vector) -> Vector: + v_ = self.transform.apply_shift(v) + result_ = self.subsurface.normal_at(v_) + result = self.transform.invert_no_shift(result_) + return result.normalize() + + def intersect(self, ray: Ray) -> Optional[float]: + e_ = self.transform.apply_no_shift(ray.e) + length_ratio = abs(e_) + + ray_ = Ray( + r0=self.transform.apply_shift(ray.r0), + e=e_.normalize(), + ) + + t_ = self.subsurface.intersect(ray_) + if t_ is None: + return None + + t = t_ / length_ratio + + return t + + def draw3d(self, painter: Painter): + transformed_painter = TransformedPainter(painter, self.transform) + + self.subsurface.draw3d(transformed_painter) diff --git a/src/geometry/vector.py b/src/geometry/vector.py new file mode 100644 index 0000000..eb03574 --- /dev/null +++ b/src/geometry/vector.py @@ -0,0 +1,104 @@ +from __future__ import annotations + +import numpy as np + +from util import iszero + + +class Vector: + def __init__(self, a): + self.a: np.ndarray = np.array(a, dtype=float) + self.a.setflags(write=False) + + assert len(self.a.shape) == 2 + + @staticmethod + def fromiterable(iterable): + a = np.zeros((len(iterable), 1), dtype=float) + + for i, j in enumerate(iterable): + a[i, 0] = j + + return Vector(a) + + @staticmethod + def fromargs(*args): + return Vector.fromiterable(args) + + @staticmethod + def random(): + return Vector.fromiterable(np.random.uniform(low=-10, high=10, size=3)) + + def __add__(self, other: Vector) -> Vector: + return Vector(self.a + other.a) + + def __sub__(self, other: Vector) -> Vector: + return Vector(self.a - other.a) + + def __matmul__(self, other: Vector) -> float: + """ + Dot product + """ + # noinspection PyTypeChecker + return np.sum(self.a * other.a) + + def __or__(self, other: Vector) -> Vector: + """ + Cross product + """ + + if self.a.shape != other.a.shape or self.a.shape != (3, 1): + raise ValueError() + + return Vector.fromargs( + self.a[1, 0] * other.a[2, 0] - self.a[2, 0] * other.a[1, 0], + self.a[2, 0] * other.a[0, 0] - self.a[0, 0] * other.a[2, 0], + self.a[0, 0] * other.a[1, 0] - self.a[1, 0] * other.a[0, 0], + ) + + def __abs__(self) -> float: + """ + Norm of self + """ + return np.sqrt(self @ self) + + def __eq__(self, other: Vector): + if self.a.shape != other.a.shape: + return False + + return (self - other).iszero() + + def __str__(self): + return str(self.a) + + def __repr__(self): + return f'Vector({repr(self.a)})' + + def __mul__(self, other): + return Vector(self.a * other) + + def __rmul__(self, other): + return Vector(other * self.a) + + def __truediv__(self, other): + return Vector(self.a / other) + + def __mod__(self, other): + """ + Is self perpendicular to other? + """ + return iszero(self @ other) + + def __neg__(self): + return -1 * self + + def iszero(self): + return iszero(self @ self) + + def normalize(self): + if self.iszero(): + raise ValueError() + return self / abs(self) + + def __hash__(self): + return hash(self.a) diff --git a/src/rendering/_3d/__init__.py b/src/rendering/_3d/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/rendering/_3d/render.py b/src/rendering/_3d/render.py new file mode 100644 index 0000000..3b29deb --- /dev/null +++ b/src/rendering/_3d/render.py @@ -0,0 +1,80 @@ +import numpy as np +from matplotlib import pyplot as plt +# noinspection PyUnresolvedReferences +from mpl_toolkits.mplot3d import Axes3D + +from geometry.ellipse import Ellipse +from geometry.line_drawers import draw_line_3d +from geometry.ray import Ray +from geometry.rotations import rotate_rodrigues +from geometry.vector import Vector +from rendering import raytracing +from rendering.raytracing import RaytracingSettings + +figure = plt.figure(figsize=(8, 8), dpi=200) +axis = figure.add_subplot(111, projection='3d') +axis.view_init(elev=10, azim=-80) + +surfaces = [] + +ellipse = Ellipse( + refraction_index=1.5, + f1=Vector.fromargs(-2, 2, 1), + f2=Vector.fromargs(1, 0, -2), + distance=5.8, +) +surfaces.append(ellipse) +ellipse.draw3d(axis) + +axis.set_xlim(-4, 4) +axis.set_ylim(-4, 4) +axis.set_zlim(-4, 4) + +axis.set_xlabel('x') +axis.set_ylabel('y') +axis.set_zlabel('z') + +axis.grid(True) + +settings = RaytracingSettings( + STARTING_LW=2, + STARTING_HUE=0, + + REFLECT_HUE_SHIFT=0, + REFLECT_LW_DECAY=0.7, + + REFRACT_HUE_SHIFT=1 / 8, + REFRACT_LW_DECAY=0.8, + + MAX_DEPTH=3 +) + +for alpha in np.linspace(-np.pi / 32, np.pi / 32, 4): + for beta in np.linspace(-np.pi / 32, np.pi / 32, 4): + e = rotate_rodrigues( + Vector.fromargs(-1, 0, 0), + Vector.fromargs(0, 1, 0), + alpha, + ) + e = rotate_rodrigues( + e, + Vector.fromargs(0, 0, 1), + beta, + ) + + raytracing.raytrace( + Ray( + r0=Vector.fromargs(4, 0, 0), + e=e + ), + settings, + axis, + surfaces, + draw_line_3d + ) + +figure.show() + +for i, azim in enumerate(np.linspace(0, 360, 16)): + axis.view_init(elev=10, azim=azim) + figure.savefig(f'movie/frame{i:03d}.png') diff --git a/src/rendering/_3d/render_transformed.py b/src/rendering/_3d/render_transformed.py new file mode 100644 index 0000000..27cdaaf --- /dev/null +++ b/src/rendering/_3d/render_transformed.py @@ -0,0 +1,111 @@ +import numpy as np +from matplotlib import pyplot as plt +# noinspection PyUnresolvedReferences +from mpl_toolkits.mplot3d import Axes3D + +from geometry.affine import AffineTransform +from geometry.line_drawers import draw_line_3d, Painter +from geometry.ray import Ray +from geometry.rotations import rotate_rodrigues +from geometry.sphere import Sphere +from geometry.transformed_surface import TransformedSurface +from geometry.vector import Vector +from rendering import raytracing +from rendering.raytracing import RaytracingSettings + +figure = plt.figure(figsize=(8, 8), dpi=200) +axis = figure.add_subplot(111, projection='3d') +axis.view_init(elev=10, azim=-80) + +painter = Painter(axis) + +surfaces = [] + +# sphere = Sphere( +# refraction_index=1.5, +# r0=Vector.fromargs(0, 0, 1), +# radius=2, +# ) +# surfaces.append(sphere) +# sphere.draw3d(painter) + +np.random.seed(5) +transform_matrix = np.random.normal(loc=1, scale=2, size=(3, 3)) +transform_shift = np.random.normal(loc=0, scale=2, size=(3, 1)) + +transform = AffineTransform( + # A=np.array([ + # [2, 0, 0], + # [0, 0.5, 0], + # [0, 0, 1], + # ]), + A=transform_matrix, + # b=Vector.fromargs(0, 0, 0), + b=Vector(transform_shift), +) +sphere = Sphere( + refraction_index=1.5, + r0=Vector.fromargs(0, 0, 1), + radius=1, +) +transformed_sphere = TransformedSurface( + refraction_index=sphere.refraction_index, + subsurface=sphere, + transform=transform +) +# sphere.draw3d(painter) +transformed_sphere.draw3d(painter) +surfaces.append(transformed_sphere) + +axis.set_xlim(-4, 4) +axis.set_ylim(-4, 4) +axis.set_zlim(-4, 4) + +axis.set_xlabel('x') +axis.set_ylabel('y') +axis.set_zlabel('z') + +axis.grid(True) + +settings = RaytracingSettings( + STARTING_LW=2, + STARTING_HUE=0, + + REFLECT_HUE_SHIFT=0, + REFLECT_LW_DECAY=0.7, + + REFRACT_HUE_SHIFT=1 / 8, + REFRACT_LW_DECAY=0.8, + + MAX_DEPTH=3 +) + +for alpha in np.linspace(-np.pi / 32, np.pi / 32, 4): + for beta in np.linspace(-np.pi / 32, np.pi / 32, 4): + e = rotate_rodrigues( + Vector.fromargs(-1, 0, 0), + Vector.fromargs(0, 1, 0), + alpha, + ) + e = rotate_rodrigues( + e, + Vector.fromargs(0, 0, 1), + beta, + ) + + raytracing.raytrace( + Ray( + r0=Vector.fromargs(4, 0, 0), + e=e + ), + settings, + axis, + surfaces, + draw_line_3d + ) + +figure.show() + +for i, azim in enumerate(np.linspace(0, 360, 16)): + axis.view_init(elev=10, azim=azim) + figure.savefig(f'movie/frame{i:03d}.png') diff --git a/src/rendering/__init__.py b/src/rendering/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/rendering/raytracing.py b/src/rendering/raytracing.py new file mode 100644 index 0000000..be6e678 --- /dev/null +++ b/src/rendering/raytracing.py @@ -0,0 +1,110 @@ +import dataclasses +from typing import List + +from matplotlib.colors import hsv_to_rgb + +from geometry.ray import Ray +from geometry.reflections import reflect, refract +from geometry.surface import Surface +from util import iszero + + +@dataclasses.dataclass +class RaytracingSettings: + STARTING_LW: float + REFLECT_LW_DECAY: float + REFRACT_LW_DECAY: float + + STARTING_HUE: float + REFLECT_HUE_SHIFT: float + REFRACT_HUE_SHIFT: float + + MAX_DEPTH: int + + lw: float = None + hue: float = None + depth: float = None + n: float = 1 + + def __post_init__(self): + if self.lw is None: + self.lw = self.STARTING_LW + + if self.hue is None: + self.hue = self.STARTING_HUE + + if self.depth is None: + self.depth = self.MAX_DEPTH + + def copy(self): + return dataclasses.replace(self) + + def reflect(self): + return dataclasses.replace( + self, + lw=self.lw * self.REFLECT_LW_DECAY, + hue=(self.hue + self.REFLECT_HUE_SHIFT) % 1, + depth=self.depth - 1, + ) + + def refract(self, new_n): + return dataclasses.replace( + self, + lw=self.lw * self.REFRACT_LW_DECAY, + hue=(self.hue + self.REFRACT_HUE_SHIFT) % 1, + depth=self.depth - 1, + n=new_n, + ) + + +def raytrace(ray: Ray, + settings: RaytracingSettings, + axis, + surfaces: List[Surface], + line_drawer): + if settings.depth == 0: + return + + min_t = None + surface = None + + for cur_surface in surfaces: + t = cur_surface.intersect(ray) + if t is None or iszero(t): + continue + + if min_t is None or t < min_t: + min_t = t + surface = cur_surface + + if min_t is None: + line_drawer(axis, ray.r0, ray.prolong(1), + color=hsv_to_rgb([settings.hue, settings.lw / settings.STARTING_LW, 1]), + zorder=settings.depth + 5, + linestyle='--', + lw=settings.lw) + return + + intersection = ray.prolong(min_t) + line_drawer(axis, ray.r0, intersection, + color=hsv_to_rgb([settings.hue, settings.lw / settings.STARTING_LW, 1]), + zorder=settings.depth + 5, + linestyle='-', + lw=settings.lw) + + if ray.r0 not in surface: + # moves into + normal = surface.normal_at(intersection) + new_n = settings.n * surface.refraction_index + else: + # moves outside + normal = -surface.normal_at(intersection) + new_n = settings.n / surface.refraction_index + + # reflection + reflected = reflect(ray.e, normal) + raytrace(Ray(intersection, reflected), settings.reflect(), axis, surfaces, line_drawer) + + # refraction + refracted = refract(ray.e, normal, settings.n, new_n) + raytrace(Ray(intersection, refracted), settings.refract(new_n), axis, surfaces, line_drawer) diff --git a/src/rendering/render.py b/src/rendering/render.py new file mode 100644 index 0000000..bfba130 --- /dev/null +++ b/src/rendering/render.py @@ -0,0 +1,111 @@ +from typing import Set + +import numpy as np +from matplotlib import pyplot as plt + +from plane import Plane +from ray import Ray +from reflections import reflect, refract +from sphere import Sphere +from vector import Vector +from util import iszero + + +def plot_v(v1: Vector, v2: Vector, *args, **kwargs): + plt.plot( + [v1.a[0, 0], v2.a[0, 0]], + [v1.a[1, 0], v2.a[1, 0]], + *args, **kwargs) + + +intersectables = [] +f = plt.figure(figsize=(6, 6)) +ax = plt.gca() + +# declare a plane +plane = Plane( + r0=Vector.fromargs(1, 0), + n=Vector.fromargs(-2, 1), +) +intersectables.append((plane, 1.5)) +plt.plot([-1, 3], [-4, 4], label='плоскость') + +# declare a sphere +sphere = Sphere( + r0=Vector.fromargs(0, 2), + radius=0.8, +) +intersectables.append((sphere, 2)) +circle = plt.Circle( + (sphere.r0.a[0, 0], sphere.r0.a[1, 0]), + sphere.radius, color='blue', fill=False) +ax.add_artist(circle) + + +def raytrace( + ray: Ray, + current_n: float, + already_intersected: Set, + depth: int): + if depth == 0: + return + + min_t = None + surface = None + surface_n = None + + for intersectable, intersectable_n in intersectables: + t = intersectable.intersect(ray) + if t is None or iszero(t): + continue + + if min_t is None or t < min_t: + min_t = t + surface = intersectable + surface_n = intersectable_n + + if min_t is None: + plot_v(ray.r0, ray.prolong(5), color='red') + return + + intersection = ray.prolong(min_t) + plot_v(ray.r0, intersection, color='orange') + + normal = surface.normal_at(intersection) + new_n = surface_n + new_already_intersected = set(already_intersected) + + if id(surface) not in already_intersected: + # moves_into + new_already_intersected.add(id(surface)) + else: + # moves outside + new_already_intersected.remove(id(surface)) + normal = -1 * normal + new_n = 1 + + # reflection + reflected = reflect(ray.e, normal) + raytrace(Ray(intersection, reflected), current_n, already_intersected, depth - 1) + + # refraction + refracted = refract(ray.e, normal, current_n, new_n) + raytrace(Ray(intersection, refracted), new_n, new_already_intersected, depth - 1) + + +raytrace(Ray( + r0=Vector.fromargs(-4, 2), + e=Vector.fromargs(2, -0.5) +), 1, set(), 10) + +# for i in range(30): +# trace_ray(Ray( +# r0=Vector.fromargs(-4, -2), +# e=Vector.fromargs(np.cos(i / 20 / np.pi), np.sin(i / 20 / np.pi)), +# )) + +plt.xlim(-5, 5) +plt.ylim(-5, 5) +# plt.legend() +plt.grid(True) +plt.show() diff --git a/src/rendering/render2.py b/src/rendering/render2.py new file mode 100644 index 0000000..d2ef9d9 --- /dev/null +++ b/src/rendering/render2.py @@ -0,0 +1,125 @@ +from typing import Set + +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.colors import hsv_to_rgb + +from plane import Plane +from ray import Ray +from reflections import reflect, refract +from sphere import Sphere +from vector import Vector +from util import iszero + + +def draw_line(v1: Vector, v2: Vector, **kwargs): + plt.plot( + [v1.a[0, 0], v2.a[0, 0]], + [v1.a[1, 0], v2.a[1, 0]], + **kwargs) + + +def draw_point(v: Vector, **kwargs): + plt.gca().add_artist( + plt.Circle( + (v.a[0, 0], v.a[1, 0]), 0.1, + fill=False, **kwargs)) + + +f = plt.figure(figsize=(8, 8)) + +surfaces = [] + +np.random.seed(5) + +for i in range(5): + sphere = Sphere( + refraction_index=3, + r0=Vector.fromiterable(np.random.uniform(-2, 2, size=2)), + radius=np.random.uniform(0.5, 2), + ) + surfaces.append(sphere) + + plt.gca().add_artist(plt.Circle( + (sphere.r0.a[0, 0], sphere.r0.a[1, 0]), + sphere.radius, color='blue', fill=False, zorder=200)) + plt.gca().add_artist(plt.Circle( + (sphere.r0.a[0, 0], sphere.r0.a[1, 0]), + sphere.radius, color='blue', fill=True, zorder=0, alpha=0.05)) + + +def raytrace( + ray: Ray, + hue: float, + lw: float, + current_n: float, + depth: int): + if depth == 0: + return + + min_t = None + surface = None + + for cur_surface in surfaces: + t = cur_surface.intersect(ray) + if t is None or iszero(t): + continue + + if min_t is None or t < min_t: + min_t = t + surface = cur_surface + + if min_t is None: + draw_line(ray.r0, ray.prolong(5), color=hsv_to_rgb([hue, 1, 1]), zorder=depth + 5, lw=lw) + return + + intersection = ray.prolong(min_t) + draw_line(ray.r0, intersection, color=hsv_to_rgb([hue, 1, 1]), zorder=depth + 5, lw=lw) + + # normal = surface.normal_at(intersection) + # new_n = current_n * + + if ray.r0 not in surface: + # moves into + normal = surface.normal_at(intersection) + new_n = current_n * surface.refraction_index + else: + # moves outside + normal = -1 * surface.normal_at(intersection) + new_n = current_n / surface.refraction_index + + # reflection + reflected = reflect(ray.e, normal) + raytrace( + Ray(intersection, reflected), + hue=hue, + lw=lw * 1, + current_n=current_n, + depth=depth - 1) + + # refraction + refracted = refract(ray.e, normal, current_n, new_n) + raytrace( + Ray(intersection, refracted), + hue=(hue + 0.15) % 1, + lw=lw * 0.7, + current_n=new_n, + depth=depth - 1) + + +raytrace(Ray( + r0=Vector.fromargs(-3.5, 0), + e=Vector.fromargs(2, -0.5) +), hue=0, lw=5, current_n=1, depth=8) + +# for i in range(30): +# trace_ray(Ray( +# r0=Vector.fromargs(-4, -2), +# e=Vector.fromargs(np.cos(i / 20 / np.pi), np.sin(i / 20 / np.pi)), +# )) + +plt.xlim(-4, 4) +plt.ylim(-4, 4) +# plt.legend() +plt.grid(True, zorder=1) +plt.show() diff --git a/src/rendering/render3.py b/src/rendering/render3.py new file mode 100644 index 0000000..d7ee018 --- /dev/null +++ b/src/rendering/render3.py @@ -0,0 +1,131 @@ +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.colors import hsv_to_rgb + +from ellipse import Ellipse +from ray import Ray +from reflections import reflect, refract +from util import iszero +from vector import Vector + +MAX_DEPTH = 10 + +STARTING_LW = 5 + +REFRACT_HUE_SHIFT = 1 / 8 + +LW_DECAY = 0.8 + + +def draw_line(v1: Vector, v2: Vector, **kwargs): + plt.plot( + [v1.a[0, 0], v2.a[0, 0]], + [v1.a[1, 0], v2.a[1, 0]], + **kwargs) + + +def draw_point(v: Vector, **kwargs): + plt.gca().add_artist( + plt.Circle( + (v.a[0, 0], v.a[1, 0]), 0.1, + fill=False, **kwargs)) + + +f = plt.figure(figsize=(8, 8)) + +surfaces = [] + +np.random.seed(2) + +# ellipse = Ellipse( +# refraction_index=1.3, +# f1=Vector.fromargs(1, 0), +# f2=Vector.fromargs(-1, 0), +# distance=3, +# ) +# ellipse.draw2d(plt) +# surfaces.append(ellipse) + +for i in range(5): + f1 = Vector.fromiterable(np.random.uniform(-3, 3, size=2)) + f2 = Vector.fromiterable(np.random.uniform(-3, 3, size=2)) + distance = abs(f1 - f2) + np.random.uniform(0, 2) + ellipse = Ellipse( + refraction_index=2, + f1=f1, + f2=f2, + distance=distance, + ) + ellipse.draw2d(plt, color='blue') + surfaces.append(ellipse) + + +def raytrace( + ray: Ray, + hue: float, + lw: float, + current_n: float, + depth: int): + if depth == 0: + return + + min_t = None + surface = None + + for cur_surface in surfaces: + t = cur_surface.intersect(ray) + if t is None or iszero(t): + continue + + if min_t is None or t < min_t: + min_t = t + surface = cur_surface + + if min_t is None: + draw_line(ray.r0, ray.prolong(5), color=hsv_to_rgb([hue, 1, 1]), zorder=depth + 5, lw=lw) + return + + intersection = ray.prolong(min_t) + draw_line(ray.r0, intersection, color=hsv_to_rgb([hue, 1, 1]), zorder=depth + 5, lw=lw) + + # normal = surface.normal_at(intersection) + # new_n = current_n * + + if ray.r0 not in surface: + # moves into + normal = surface.normal_at(intersection) + new_n = current_n * surface.refraction_index + else: + # moves outside + normal = -1 * surface.normal_at(intersection) + new_n = current_n / surface.refraction_index + + # reflection + reflected = reflect(ray.e, normal) + raytrace( + Ray(intersection, reflected), + hue=hue, + lw=lw * LW_DECAY, + current_n=current_n, + depth=depth - 1) + + # refraction + refracted = refract(ray.e, normal, current_n, new_n) + raytrace( + Ray(intersection, refracted), + hue=(hue + REFRACT_HUE_SHIFT) % 1, + lw=lw * LW_DECAY, + current_n=new_n, + depth=depth - 1) + + +raytrace(Ray( + r0=Vector.fromargs(-3, 3), + e=Vector.fromargs(1, -1) +), hue=0, lw=STARTING_LW, current_n=1, depth=MAX_DEPTH) + +plt.xlim(-4, 4) +plt.ylim(-4, 4) +# plt.legend() +plt.grid(True, zorder=1) +plt.show() diff --git a/src/rendering/render4.py b/src/rendering/render4.py new file mode 100644 index 0000000..14192d5 --- /dev/null +++ b/src/rendering/render4.py @@ -0,0 +1,136 @@ +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.colors import hsv_to_rgb + +from geometry.ellipse import Ellipse +from geometry.ray import Ray +from geometry.reflections import reflect, refract +from geometry.vector import Vector +from geometry.rotations import rotate2d +from util import iszero + +MAX_DEPTH = 2 + +STARTING_LW = 3 + +REFRACT_HUE_SHIFT = 1 / 8 + +LW_DECAY = 0.7 + + +def draw_line(v1: Vector, v2: Vector, **kwargs): + plt.plot( + [v1.a[0, 0], v2.a[0, 0]], + [v1.a[1, 0], v2.a[1, 0]], + **kwargs) + + +def draw_point(v: Vector, **kwargs): + plt.gca().add_artist( + plt.Circle( + (v.a[0, 0], v.a[1, 0]), 0.1, + fill=False, **kwargs)) + + +f = plt.figure(figsize=(8, 8), dpi=200) + +surfaces = [] + +np.random.seed(2) + + +# sphere = Sphere( +# refraction_index=1.3, +# r0=Vector.fromargs(-1, 2), +# radius=1, +# ) +# # sphere.draw2d(plt) +# surfaces.append(sphere) + + +def raytrace( + ray: Ray, + hue: float, + lw: float, + current_n: float, + depth: int): + if depth == 0: + return + + min_t = None + surface = None + + for cur_surface in surfaces: + t = cur_surface.intersect(ray) + if t is None or iszero(t): + continue + + if min_t is None or t < min_t: + min_t = t + surface = cur_surface + + if min_t is None: + draw_line(ray.r0, ray.prolong(10), + color=hsv_to_rgb([hue, lw / STARTING_LW, 1]), + zorder=depth + 5, lw=lw) + return + + intersection = ray.prolong(min_t) + draw_line(ray.r0, intersection, + color=hsv_to_rgb([hue, lw / STARTING_LW, 1]), + zorder=depth + 5, lw=lw) + + # normal = surface.normal_at(intersection) + # new_n = current_n * + + if ray.r0 not in surface: + # moves into + normal = surface.normal_at(intersection) + new_n = current_n * surface.refraction_index + else: + # moves outside + normal = -1 * surface.normal_at(intersection) + new_n = current_n / surface.refraction_index + + # reflection + reflected = reflect(ray.e, normal) + raytrace( + Ray(intersection, reflected), + hue=hue, + lw=lw * LW_DECAY, + current_n=current_n, + depth=depth - 1) + + # refraction + refracted = refract(ray.e, normal, current_n, new_n) + raytrace( + Ray(intersection, refracted), + hue=(hue + REFRACT_HUE_SHIFT) % 1, + lw=lw * LW_DECAY, + current_n=new_n, + depth=depth - 1) + + +ellipse = Ellipse( + refraction_index=1, + f1=Vector.fromargs(0, -2), + f2=Vector.fromargs(2, 2), + distance=5, +) +ellipse.draw2d(plt, color='blue') +surfaces.append(ellipse) + +r0 = Vector.fromargs(-3, 0) +e = Vector.fromargs(1, 0) + +for alpha in np.linspace(-np.pi / 4, np.pi / 4, 31): + raytrace(Ray( + r0=r0, + e=rotate2d(e, alpha) + ), hue=0, lw=STARTING_LW, current_n=1, depth=MAX_DEPTH) + +plt.xlim(-4, 4) +plt.ylim(-4, 4) +# plt.legend() +plt.grid(True, zorder=1) +plt.show() diff --git a/src/rendering/render5.py b/src/rendering/render5.py new file mode 100644 index 0000000..0cd9c8c --- /dev/null +++ b/src/rendering/render5.py @@ -0,0 +1,142 @@ +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.colors import hsv_to_rgb + +from ellipse import Ellipse +from ray import Ray +from reflections import reflect, refract +from util import iszero +from vector import Vector, rotate2d + +MAX_DEPTH = 6 + +STARTING_LW = 5 + +REFRACT_HUE_SHIFT = 1 / 8 + +LW_DECAY = 0.7 + + +def draw_line(v1: Vector, v2: Vector, **kwargs): + plt.plot( + [v1.a[0, 0], v2.a[0, 0]], + [v1.a[1, 0], v2.a[1, 0]], + **kwargs) + + +def draw_point(v: Vector, **kwargs): + plt.gca().add_artist( + plt.Circle( + (v.a[0, 0], v.a[1, 0]), 0.1, + fill=False, **kwargs)) + + +f = plt.figure(figsize=(8, 8), dpi=200) + +surfaces = [] + +np.random.seed(2) + + +# ellipse = Ellipse( +# refraction_index=1.3, +# f1=Vector.fromargs(1, 0), +# f2=Vector.fromargs(-1, 0), +# distance=3, +# ) +# ellipse.draw2d(plt) +# surfaces.append(ellipse) + + +def raytrace( + ray: Ray, + hue: float, + lw: float, + current_n: float, + depth: int): + if depth == 0: + return + + min_t = None + surface = None + + for cur_surface in surfaces: + t = cur_surface.intersect(ray) + if t is None or iszero(t): + continue + + if min_t is None or t < min_t: + min_t = t + surface = cur_surface + + if min_t is None: + draw_line(ray.r0, ray.prolong(10), + color=hsv_to_rgb([hue, lw / STARTING_LW, 1]), + zorder=depth + 5, lw=lw) + return + + intersection = ray.prolong(min_t) + draw_line(ray.r0, intersection, + color=hsv_to_rgb([hue, lw / STARTING_LW, 1]), + zorder=depth + 5, lw=lw) + + if ray.r0 not in surface: + # moves into + normal = surface.normal_at(intersection) + new_n = current_n * surface.refraction_index + else: + # moves outside + normal = -1 * surface.normal_at(intersection) + new_n = current_n / surface.refraction_index + + # reflection + reflected = reflect(ray.e, normal) + raytrace( + Ray(intersection, reflected), + hue=hue, + lw=lw * LW_DECAY, + current_n=current_n, + depth=depth - 1) + + # refraction + refracted = refract(ray.e, normal, current_n, new_n) + raytrace( + Ray(intersection, refracted), + hue=(hue + REFRACT_HUE_SHIFT) % 1, + lw=lw * LW_DECAY, + current_n=new_n, + depth=depth - 1) + + +ellipse1 = Ellipse( + refraction_index=3, + f1=Vector.fromargs(2, -2), + f2=Vector.fromargs(3, 2), + distance=5, +) +ellipse1.draw2d(plt, color='blue') +surfaces.append(ellipse1) + +ellipse2 = Ellipse( + refraction_index=3, + f1=Vector.fromargs(-1, -2), + f2=Vector.fromargs(-1, 2), + distance=4.2, +) +ellipse2.draw2d(plt, color='blue') +surfaces.append(ellipse2) + +r0 = Vector.fromargs(-3, 0) +e = Vector.fromargs(1, 0) + +for alpha in np.linspace(-np.pi / 4, np.pi / 4, 9): + raytrace(Ray( + r0=r0, + e=rotate2d(e, alpha) + ), hue=0, lw=STARTING_LW, current_n=1, depth=MAX_DEPTH) + +plt.xlim(-4, 4) +plt.ylim(-4, 4) +# plt.legend() +plt.grid(True, zorder=1) +plt.show() diff --git a/src/test.py b/src/test.py new file mode 100644 index 0000000..7fd96cf --- /dev/null +++ b/src/test.py @@ -0,0 +1,18 @@ +import numpy as np + +from plane import Plane +from vector import Vector + +if __name__ == '__main__': + a = Vector.fromargs(1, 2, 3) + b = Vector.fromiterable(np.eye(3)[:, 0]) + print(a @ b) + print(abs(a)) + print(a / 2) + + p = Plane( + Vector.fromargs(0, 0, 0), + Vector.fromargs(0, 1, 0), + ) + + print(p & Vector.fromargs(1, 1, 0)) diff --git a/src/util.py b/src/util.py new file mode 100644 index 0000000..1afb9ef --- /dev/null +++ b/src/util.py @@ -0,0 +1,46 @@ +from __future__ import annotations + +from typing import List, Optional + +import numpy as np + +EPS = 1e-6 + + +def double_equality(a, b) -> bool: + return abs(a - b) < EPS + + +def iszero(x) -> bool: + return double_equality(np.sum(x), 0) + + +def minpos_criterion(x: float) -> bool: + return x is not None and not iszero(x) and x > 0 + + +def minpos(args) -> Optional[float]: + args = list(filter(minpos_criterion, args)) + + if not args: + return None + + return min(args) + + +def solve_quadratic(a, b, c) -> List[float]: + if iszero(a): + if iszero(b): + return [] + + return [-c / b] + + D = b ** 2 - 4 * a * c + + if iszero(D): + return [-b / 2 / a] + + if D < 0: + return [] + + return [(-b - np.sqrt(D)) / 2 / a, (-b + np.sqrt(D)) / 2 / a] diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_intersections.py b/tests/test_intersections.py new file mode 100644 index 0000000..d5a7fa1 --- /dev/null +++ b/tests/test_intersections.py @@ -0,0 +1,30 @@ +import numpy as np + +from intersections import intersect_plane_ray +from plane import Plane +from ray import Ray +from vector import Vector + + +def test_intersections(): + r = Ray( + Vector.fromargs(0, 0, 1), + Vector.fromargs(-1, -1, -1), + ) + p = Plane( + Vector.fromargs(0, 0, 0), + Vector.fromargs(0, 0, 1), + ) + + assert intersect_plane_ray(p, r) == Vector.fromargs(-1, -1, 0) + + np.random.seed(0) + for i in range(10): + p = Plane(Vector.random(), Vector.random()) + r = Ray(Vector.random(), Vector.random()) + + intersection = intersect_plane_ray(p, r) + + print(p, r, intersection) + + assert p & intersection diff --git a/tests/test_plane.py b/tests/test_plane.py new file mode 100644 index 0000000..411cc64 --- /dev/null +++ b/tests/test_plane.py @@ -0,0 +1,37 @@ +import numpy as np + +from plane import Plane +from ray import Ray +from vector import Vector + + +def test_plane(): + p = Plane( + Vector.fromargs(0, 0, 0), + Vector.fromargs(0, 0, 1), + ) + assert (p & Vector.fromargs(0, 0, 0)) + assert (p & Vector.fromargs(1, 1, 0)) + assert not (p & Vector.fromargs(0, 0, -1)) + + p = Plane( + Vector.fromargs(1, 1), + Vector.fromargs(1, 1), + ) + + assert not (p & Vector.fromargs(0, 0)) + assert (p & Vector.fromargs(2, 0)) + + +def test_reflection(): + p = Plane( + Vector.fromargs(0, 0), + Vector.fromargs(0, 1), + ) + + r = Ray( + Vector.fromargs(-1, 1), + Vector.fromargs(1, -1), + ) + + assert p.intersect(r) == np.sqrt(2) diff --git a/tests/test_sphere.py b/tests/test_sphere.py new file mode 100644 index 0000000..e123920 --- /dev/null +++ b/tests/test_sphere.py @@ -0,0 +1,22 @@ +import numpy as np + +from ray import Ray +from sphere import Sphere +from vector import Vector + + +def test_intersect(): + sphere = Sphere( + r0=Vector.fromargs(0, 0), + radius=1, + ) + + assert sphere.intersect(Ray( + r0=Vector.fromargs(-2, 0), + e=Vector.fromargs(1, 0), + )) == 1 + + assert sphere.intersect(Ray( + r0=Vector.fromargs(-1, -1), + e=Vector.fromargs(1, 1), + )) == np.sqrt(2) - 1 diff --git a/tests/test_vector.py b/tests/test_vector.py new file mode 100644 index 0000000..2bb8946 --- /dev/null +++ b/tests/test_vector.py @@ -0,0 +1,7 @@ +from vector import Vector + + +def test_vector(): + v = Vector.fromargs(1, 2) + print(v) + print(-v)