FabLab_Chemnitz/mightyscape-1.1-deprecated
FabLab_Chemnitz/mightyscape-1.1-deprecated
Archived
1
0
Fork 0
Code Issues Pull Requests Releases Activity

added ray tracing

Browse Source
This commit is contained in:
Mario Voigt 2021-10-24 19:19:52 +02:00
parent 1dc02d905e
commit 77fcd7a08a
25 changed files with 1447 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
extensions/fablabchemnitz/raytracing/desc_parser.py Normal file
View File

@ -0,0 +1,20 @@
import re
rgx_float = r"[-+]?(\d+([.,]\d*)?|[.,]\d+)([eE][-+]?\d+)?"
rgx_name = "[a-z,_]*"
optics_pattern = re.compile(
f"optics *: *(?P<material>{rgx_name})(: *(?P<num>{rgx_float}))?",
re.IGNORECASE | re.MULTILINE,
)
def get_optics_fields(string_: str):
fields = re.finditer(optics_pattern, string_)
return fields
def clear_description(desc: str) -> str:
"""Removes text corresponding to an optical property"""
new_desc = re.sub(optics_pattern, "", desc)
return new_desc

42
extensions/fablabchemnitz/raytracing/lens.inx Normal file
View File

@ -0,0 +1,42 @@
<?xml version="1.0" encoding="UTF-8"?>
<inkscape-extension xmlns="http://www.inkscape.org/namespace/inkscape/extension">
<name>Ray Tracing - Insert Lens Optics</name>
<id>fablabchemnitz.de.raytracing_insert_lens_optics</id>
<param name="focal_length" type="float" gui-text="Focal length:" min="-10000." max="10000." precision="3">100.</param>
<param name="focal_length_unit" type="optiongroup" appearance="combo" gui-text=" ">
<option value="mm">mm</option>
<option value="cm">cm</option>
<option value="m">m</option>
<option value="in">in</option>
</param>
<param name="diameter" type="float" gui-text="Diameter:" min="0" max="10000" precision="3">1</param>
<param name="diameter_unit" type="optiongroup" appearance="combo" gui-text=" ">
<option value="in">in</option>
<option value="mm">mm</option>
<option value="cm">cm</option>
<option value="m">m</option>
</param>
<param name="edge_thickness" type="float" gui-text="Edge thickness:" min="0" max="10000" precision="3">2</param>
<param name="edge_thickness_unit" type="optiongroup" appearance="combo" gui-text=" ">
<option value="mm">mm</option>
<option value="in">in</option>
<option value="cm">cm</option>
<option value="m">m</option>
</param>
<param name="optical_index" type="float" min="1." max="3." precision="4" gui-text="Optical index:">1.5168</param>
<param name="lens_type" type="optiongroup" appearance="combo" gui-text="Lens type:">
<option value="plano_con">Plano-concave/convex</option>
<option value="bi_con">Bi-concave/convex</option>
</param>
<effect>
<object-type>all</object-type>
<effects-menu>
<submenu name="FabLab Chemnitz">
<submenu name="Various"/>
</submenu>
</effects-menu>
</effect>
<script>
<command location="inx" interpreter="python">lens.py</command>
</script>
</inkscape-extension>

229
extensions/fablabchemnitz/raytracing/lens.py Normal file
View File

@ -0,0 +1,229 @@
"""
Module to add a lens object in the document
"""
from math import cos, pi, sin, sqrt, acos, tan
import inkex
class Lens(inkex.GenerateExtension):
"""
Produces a PathElement corresponding to the shape of the lens calculated
from user parameters.
"""
@property
def style(self):
return {
"stroke": "#000000",
"fill": "#b7c2dd",
"stroke-linejoin": "round",
"stroke-width": str(self.svg.unittouu("1px")),
}
def add_arguments(self, pars):
pars.add_argument("--focal_length", type=float, default=100.0)
pars.add_argument("--focal_length_unit", default="mm")
pars.add_argument("--diameter", type=float, default=1.0)
pars.add_argument("--diameter_unit", default="in")
pars.add_argument("--edge_thickness", type=float, default=2.0)
pars.add_argument("--edge_thickness_unit", default="mm")
pars.add_argument("--optical_index", type=float, default=1.5168)
pars.add_argument("--lens_type", default="plano_con")
def to_document_units(self, value: float, unit: str) -> float:
c1x, c1y, c2x, c2y = self.svg.get_viewbox()
document_width = self.svg.unittouu(self.document.getroot().get("width"))
scale_factor = (c2x - c1x) / document_width
return self.svg.unittouu(str(value) + unit) * scale_factor
def generate(self):
opts = self.options
d = self.to_document_units(opts.diameter, opts.diameter_unit)
f = self.to_document_units(opts.focal_length, opts.focal_length_unit)
e = self.to_document_units(opts.edge_thickness, opts.edge_thickness_unit)
optical_index = opts.optical_index
lens_path = []
if opts.lens_type == "plano_con":
# Radius of curvature from Lensmaker's equation
roc = (optical_index - 1) * abs(f)
if 2 * roc < d:
inkex.utils.errormsg(
"Focal length is too short or diameter is too large."
)
return None
elif (roc ** 2 - (d / 2) ** 2) ** 0.5 - roc < -e and f < 0:
inkex.utils.errormsg("Edge thickness is too small.")
return None
else:
sweep = 1 if f < 0 else 0
lens_path = arc_to_path(
[-d / 2, 0], [roc, roc, 0.0, 0, sweep, +d / 2, 0]
)
lens_path += [
[[+d / 2, 0], [+d / 2, 0], [+d / 2, -e]],
[[+d / 2, -e], [+d / 2, -e], [-d / 2, -e]],
[[+d / 2, -e], [-d / 2, -e], [-d / 2, +e]],
]
# no need to close the path correctly as it's done after
elif opts.lens_type == "bi_con":
roc = (
(optical_index - 1) * abs(f) * (1 + (1 - e / f / optical_index) ** 0.5)
)
if 2 * roc < d:
inkex.utils.errormsg(
"Focal length is too short or diameter is too large."
)
return None
elif (roc ** 2 - (d / 2) ** 2) ** 0.5 - roc < -e / 2 and f < 0:
inkex.utils.errormsg("Edge thickness is too small.")
return None
else:
sweep = 1 if f < 0 else 0
lens_path = arc_to_path(
[-d / 2, 0], [roc, roc, 0.0, 0, sweep, +d / 2, 0]
)
lens_path += [
[[+d / 2, 0], [+d / 2, 0], [+d / 2, -e]],
[[+d / 2, -e], [+d / 2, -e], [+d / 2, -e]],
]
lens_path += arc_to_path(
[+d / 2, -e], [roc, roc, 0.0, 0, sweep, -d / 2, -e]
)
lens_path += [
[[-d / 2, -e], [-d / 2, -e], [-d / 2, 0]],
[[-d / 2, -e], [-d / 2, 0], [-d / 2, 0]],
]
lens = inkex.PathElement()
lens.style = self.style
closed_path = inkex.Path(inkex.CubicSuperPath([lens_path]))
closed_path.close()
lens.path = closed_path.transform(inkex.Transform("rotate(90)"))
lens.desc = (
f"L{opts.focal_length}{opts.focal_length_unit}\n"
f"optics:glass:{optical_index:.4f}"
)
yield lens
def arc_to_path(point, params):
"""Approximates an arc with cubic bezier segments.
Arguments:
point: Starting point (absolute coords)
params: Arcs parameters as per
https://www.w3.org/TR/SVG/paths.html#PathDataEllipticalArcCommands
Returns a list of triplets of points : [control_point_before, node, control_point_after]
(first and last returned triplets are [p1, p1, *] and [*, p2, p2])
"""
A = point[:]
rx, ry, theta, long_flag, sweep_flag, x2, y2 = params[:]
theta = theta * pi / 180.0
B = [x2, y2]
# Degenerate ellipse
if rx == 0 or ry == 0 or A == B:
return [[A[:], A[:], A[:]], [B[:], B[:], B[:]]]
# turn coordinates so that the ellipse morph into a *unit circle* (not 0-centered)
mat = mat_prod(
(rot_mat(theta), [[1.0 / rx, 0.0], [0.0, 1.0 / ry]], rot_mat(-theta))
)
apply_mat(mat, A)
apply_mat(mat, B)
k = [-(B[1] - A[1]), B[0] - A[0]]
d = k[0] * k[0] + k[1] * k[1]
k[0] /= sqrt(d)
k[1] /= sqrt(d)
d = sqrt(max(0, 1 - d / 4.0))
# k is the unit normal to AB vector, pointing to center O
# d is distance from center to AB segment (distance from O to the midpoint of AB)
# for the last line, remember this is a unit circle, and kd vector is orthogonal to AB (Pythagorean thm)
if (
long_flag == sweep_flag
): # top-right ellipse in SVG example https://www.w3.org/TR/SVG/images/paths/arcs02.svg
d *= -1
O = [(B[0] + A[0]) / 2.0 + d * k[0], (B[1] + A[1]) / 2.0 + d * k[1]]
OA = [A[0] - O[0], A[1] - O[1]]
OB = [B[0] - O[0], B[1] - O[1]]
start = acos(OA[0] / norm(OA))
if OA[1] < 0:
start *= -1
end = acos(OB[0] / norm(OB))
if OB[1] < 0:
end *= -1
# start and end are the angles from center of the circle to A and to B respectively
if sweep_flag and start > end:
end += 2 * pi
if (not sweep_flag) and start < end:
end -= 2 * pi
nb_sectors = int(abs(start - end) * 2 / pi) + 1
d_theta = (end - start) / nb_sectors
v = 4 * tan(d_theta / 4.0) / 3.0
# I would use v = tan(d_theta/2)*4*(sqrt(2)-1)/3 ?
p = []
for i in range(0, nb_sectors + 1, 1):
angle = start + i * d_theta
v1 = [
O[0] + cos(angle) - (-v) * sin(angle),
O[1] + sin(angle) + (-v) * cos(angle),
]
pt = [O[0] + cos(angle), O[1] + sin(angle)]
v2 = [O[0] + cos(angle) - v * sin(angle), O[1] + sin(angle) + v * cos(angle)]
p.append([v1, pt, v2])
p[0][0] = p[0][1][:]
p[-1][2] = p[-1][1][:]
# go back to the original coordinate system
mat = mat_prod((rot_mat(theta), [[rx, 0], [0, ry]], rot_mat(-theta)))
for pts in p:
apply_mat(mat, pts[0])
apply_mat(mat, pts[1])
apply_mat(mat, pts[2])
return p
def mat_prod(m_list):
"""Get the product of the mat"""
prod = m_list[0]
for mat in m_list[1:]:
a00 = prod[0][0] * mat[0][0] + prod[0][1] * mat[1][0]
a01 = prod[0][0] * mat[0][1] + prod[0][1] * mat[1][1]
a10 = prod[1][0] * mat[0][0] + prod[1][1] * mat[1][0]
a11 = prod[1][0] * mat[0][1] + prod[1][1] * mat[1][1]
prod = [[a00, a01], [a10, a11]]
return prod
def rot_mat(theta):
"""Rotate the mat"""
return [[cos(theta), -sin(theta)], [sin(theta), cos(theta)]]
def apply_mat(mat, point):
"""Apply the given mat"""
x = mat[0][0] * point[0] + mat[0][1] * point[1]
y = mat[1][0] * point[0] + mat[1][1] * point[1]
point[0] = x
point[1] = y
def norm(point):
"""Normalise"""
return sqrt(point[0] * point[0] + point[1] * point[1])
if __name__ == "__main__":
Lens().run()

20
extensions/fablabchemnitz/raytracing/meta.json Normal file
View File

@ -0,0 +1,20 @@
[
{
"name": "Ray Tracing - <various>",
"id": "fablabchemnitz.de.raytracing.<various>",
"path": "raytracing",
"original_name": "<various>",
"original_id": "damienBloch/inkscape-raytracing/<various>",
"license": "GNU GPL v3",
"license_url": "https://github.com/damienBloch/inkscape-raytracing/blob/master/LICENSE",
"comment": "",
"source_url": "https://gitea.fablabchemnitz.de/FabLab_Chemnitz/mightyscape-1.X/src/branch/master/extensions/fablabchemnitz/raytracing",
"fork_url": "https://github.com/damienBloch/inkscape-raytracing",
"documentation_url": "https://stadtfabrikanten.org/display/IFM/Ray+Tracing",
"inkscape_gallery_url": null,
"main_authors": [
"github.com/damienBloch",
"github.com/vmario89"
]
}
]

4
extensions/fablabchemnitz/raytracing/raytracing/__init__.py Normal file
View File

@ -0,0 +1,4 @@
from .optical_object import *
from .ray import *
from .vector import *
from .world import *

2
extensions/fablabchemnitz/raytracing/raytracing/geometry/__init__.py Normal file
View File

@ -0,0 +1,2 @@
from .cubic_bezier import CubicBezier
from .geometric_object import GeometricObject, CompoundGeometricObject, AABBox

210
extensions/fablabchemnitz/raytracing/raytracing/geometry/cubic_bezier.py Normal file
View File

@ -0,0 +1,210 @@
"""
Module for handling objects composed of cubic bezier curves
"""
from __future__ import annotations
import math
from dataclasses import dataclass
from functools import cached_property
import numpy
from .geometric_object import AABBox, GeometricObject, GeometryError
from ..ray import Ray
from ..shade import ShadeRec
from ..vector import Vector, UnitVector
@dataclass(frozen=True)
class CubicBezier:
r"""
Cubic bezier segment defined as
.. math::
\vec{X}(s) = (1-s)^3 \vec{p_0} + 3 s (1-s)^2 \vec{p_1}
+ 3 s^2 (1-s) \vec{p_2} + s^3 \vec{p_3}
for :math:`0 \le s \le 1`
"""
p0: Vector
p1: Vector
p2: Vector
p3: Vector
def eval(self, s) -> Vector:
return (
(1 - s) ** 3 * self.p0
+ 3 * s * (1 - s) ** 2 * self.p1
+ 3 * s ** 2 * (1 - s) * self.p2
+ s ** 3 * self.p3
)
@cached_property
def aabbox(self) -> AABBox:
# The box is slightly larger than the minimal box.
# It prevents the box to have a zero dimension if the object is a line
# aligned with vertical or horizontal.
lower_left = Vector(
min(self.p0.x, self.p1.x, self.p2.x, self.p3.x) - 1e-6,
min(self.p0.y, self.p1.y, self.p2.y, self.p3.y) - 1e-6,
)
upper_right = Vector(
max(self.p0.x, self.p1.x, self.p2.x, self.p3.x) + 1e-6,
max(self.p0.y, self.p1.y, self.p2.y, self.p3.y) + 1e-6,
)
return AABBox(lower_left, upper_right)
def tangent(self, s: float) -> UnitVector:
"""Returns the tangent at the curve at curvilinear coordinate s"""
diff_1 = (
-3 * (self.p0 - 3 * self.p1 + 3 * self.p2 - self.p3) * s ** 2
+ 6 * (self.p0 - 2 * self.p1 + self.p2) * s
- 3 * (self.p0 - self.p1)
)
# If the first derivative is not zero, it is parallel to the tangent
if diff_1.norm() > 1e-8:
return diff_1.normalize()
# but is the first derivative is zero, we need to get the second order
else:
diff_2 = -6 * (self.p0 - 3 * self.p1 + 3 * self.p2 - self.p3) * s + 6 * (
self.p0 - 2 * self.p1 + self.p2
)
if diff_2.norm() > 1e-8:
return diff_2.normalize()
else: # and even to the 3rd derivative if necessary
diff_3 = -6 * (self.p0 - 3 * self.p1 + 3 * self.p2 - self.p3)
return diff_3.normalize()
def normal(self, s: float) -> UnitVector:
"""Returns a vector normal at the curve at curvilinear coordinate s"""
return self.tangent(s).orthogonal()
def intersection_beam(self, ray: Ray) -> list[tuple[float, float]]:
r"""
Returns all couples :math:`(s, t)` such that there exist
:math:`\vec{X}` satisfying
.. math::
\vec{X} = (1-s)^3 \vec{p_0} + 3 s (1-s)^2 \vec{p_1}
+ 3 s^2 (1-s) \vec{p_2} + s^3 \vec{p_3}
and
.. math::
\vec{X} = \vec{o} + t \vec{d}
with :math:`0 \lq s \lq 1` and :math:`t >= 0`
"""
a = ray.direction.orthogonal()
a0 = a * (self.p0 - ray.origin)
a1 = -3 * a * (self.p0 - self.p1)
a2 = 3 * a * (self.p0 - 2 * self.p1 + self.p2)
a3 = a * (-self.p0 + 3 * self.p1 - 3 * self.p2 + self.p3)
roots = cubic_real_roots(a0, a1, a2, a3)
intersection_points = [self.eval(s) for s in roots]
travel = [(X - ray.origin) * ray.direction for X in intersection_points]
def valid_domain(s, t):
return 0 <= s <= 1 and t > Ray.min_travel
return [(s, t) for (s, t) in zip(roots, travel) if valid_domain(s, t)]
def num_hits(self, ray: Ray) -> int:
if self.aabbox.hit(ray):
return len(self.intersection_beam(ray))
else:
return 0
def hit(self, ray: Ray) -> ShadeRec:
"""
Returns a shade with the information for the first intersection
of a beam with the bezier segment
"""
shade = ShadeRec() # default no hit
if self.aabbox.hit(ray):
intersect_params = self.intersection_beam(ray)
travel_dist = [t for (__, t) in intersect_params]
if len(travel_dist) > 0: # otherwise error with np.argmin
shade.normal = True
first_hit = numpy.argmin(travel_dist)
shade.travel_dist = travel_dist[first_hit]
shade.local_hit_point = ray.origin + shade.travel_dist * ray.direction
shade.normal = self.normal(intersect_params[first_hit][0])
shade.set_normal_same_side(ray.origin)
return shade
def is_inside(self, ray: Ray) -> bool:
raise GeometryError(f"Can't define an inside for {self}.")
def cubic_real_roots(d: float, c: float, b: float, a: float) -> list[float]:
"""
Returns the real roots X of a cubic polynomial defined as
.. math::
a X^3 + b X^2 + c X + d = 0
"""
# For more information see:
# https://en.wikipedia.org/wiki/Cubic_equation#General_cubic_formula
if not is_almost_zero(a): # true cubic equation
p = (3 * a * c - b ** 2) / 3 / a ** 2
q = (2 * b ** 3 - 9 * a * b * c + 27 * a ** 2 * d) / 27 / a ** 3
if is_almost_zero(p):
t = [numpy.cbrt(-q)]
else:
discr = -(4 * p ** 3 + 27 * q ** 2)
if is_almost_zero(discr):
if is_almost_zero(q):
t = [0]
else:
t = [3 * q / p, -3 * q / 2 / p]
elif discr < 0:
t = [
numpy.cbrt(-q / 2 + numpy.sqrt(-discr / 108))
+ numpy.cbrt(-q / 2 - numpy.sqrt(-discr / 108))
]
else:
t = [
2
* numpy.sqrt(-p / 3)
* numpy.cos(
1 / 3 * numpy.arccos(3 * q / 2 / p * numpy.sqrt(-3 / p))
- 2 * numpy.pi * k / 3
)
for k in range(3)
]
return [x - b / 3 / a for x in t]
else:
return quadratic_roots(b, c, d)
def quadratic_roots(a: float, b: float, c: float) -> list[float]:
if not is_almost_zero(a):
discr = b ** 2 - 4 * a * c
if discr > 0:
return [
(-b + numpy.sqrt(discr)) / 2 / a,
(-b - numpy.sqrt(discr)) / 2 / a,
]
elif is_almost_zero(discr):
return [-b / 2 / a]
else:
return []
else:
return linear_root(b, c)
def linear_root(a: float, b: float) -> list[float]:
if is_almost_zero(a): # No solutions for 0*X+b=0
return [] # Ignore infinite solutions for a=b=0
else:
return [-b / a]
def is_almost_zero(x: float) -> bool:
return math.isclose(x, 0, abs_tol=1e-8)

151
extensions/fablabchemnitz/raytracing/raytracing/geometry/geometric_object.py Normal file
View File

@ -0,0 +1,151 @@
from __future__ import annotations
import functools
from dataclasses import dataclass
from typing import Protocol, Iterable, TypeVar, Generic
import numpy
from ..ray import Ray
from ..shade import ShadeRec
from ..vector import Vector
class GeometricObject(Protocol):
"""Protocol for a geometric object (line, rectangle, circle, ...)"""
def hit(self, ray: Ray) -> ShadeRec:
"""Tests if a collision between a beam and the object occurred
Returns a shade that contains the information about the collision in
case it happened.
"""
raise NotImplementedError
def num_hits(self, ray: Ray) -> int:
"""Returns the number of times a beam intersect the object boundary"""
raise NotImplementedError
@property
def aabbox(self) -> AABBox:
"""Computes an axis aligned bounding box for the object"""
raise NotImplementedError
def is_inside(self, ray: Ray) -> bool:
"""Indicates if a ray is inside or outside the object
It is not possible to define an inside for every object, for example if it is
not closed. In this case it should raise a GeometryError.
"""
raise NotImplementedError
class GeometryError(RuntimeError):
pass
T = TypeVar("T", bound=GeometricObject)
@dataclass(frozen=True)
class CompoundGeometricObject(GeometricObject, Generic[T]):
sub_objects: tuple[T, ...]
def __init__(self, sub_objects: Iterable[T]):
object.__setattr__(self, "sub_objects", tuple(sub_objects))
def __iter__(self) -> Iterable[T]:
return iter(self.sub_objects)
def __getitem__(self, item) -> T:
return self.sub_objects[item]
@functools.cached_property
def aabbox(self):
sub_boxes = (sub.aabbox for sub in self.sub_objects)
return AABBox.englobing(sub_boxes)
def hit(self, ray: Ray) -> ShadeRec:
"""
Returns a shade with the information for the first intersection
of a beam with one of the object composing the composite object
"""
result = ShadeRec()
if self.aabbox.hit(ray):
result = find_first_hit(ray, self.sub_objects)
result.hit_geometry = self
return result
def is_inside(self, ray: Ray) -> bool:
# A ray is inside an object if it intersect its boundary an odd
# number of times
return (self.num_hits(ray) % 2) == 1
def num_hits(self, ray: Ray) -> int:
if self.aabbox.hit(ray):
return sum([obj.num_hits(ray) for obj in self.sub_objects])
else:
return 0
def find_first_hit(ray: Ray, objects: Iterable[GeometricObject]) -> ShadeRec:
result = ShadeRec()
for obj in objects:
shade = obj.hit(ray)
if Ray.min_travel < shade.travel_dist < result.travel_dist:
result = shade
return result
@dataclass(frozen=True)
class AABBox:
"""
Implements an axis-aligned bounding box
This is used to accelerate the intersection between a beam and an object.
If the beam doesn't hit the bounding box, it is not necessary to do
expensive intersection calculations with the object.
"""
lower_left: Vector
upper_right: Vector
@classmethod
def englobing(cls, aabboxes: Iterable[AABBox]) -> AABBox:
return functools.reduce(cls.englobing_two, aabboxes)
@classmethod
def englobing_two(cls, b1: AABBox, b2: AABBox) -> AABBox:
union_lower_left = Vector(
min(b1.lower_left.x, b2.lower_left.x),
min(b1.lower_left.y, b2.lower_left.y),
)
union_upper_right = Vector(
max(b1.upper_right.x, b2.upper_right.x),
max(b1.upper_right.y, b2.upper_right.y),
)
return AABBox(union_lower_left, union_upper_right)
def hit(self, ray: Ray) -> bool:
"""Tests if a beam intersects the bounding box"""
# This algorithm uses the properties of IEEE floating-point
# arithmetic to correctly handle cases where the ray travels
# parallel to a coordinate axis.
# See Williams et al. "An efficient and robust ray-box intersection
# algorithm" for more details.
p0 = numpy.array([self.lower_left.x, self.lower_left.y])
p1 = numpy.array([self.upper_right.x, self.upper_right.y])
direction = numpy.array([ray.direction.x, ray.direction.y])
origin = numpy.array([ray.origin.x, ray.origin.y])
# The implementation safely handles the case where an element
# of ray.direction is zero. Warning for floating point error
# can be ignored for this step.
with numpy.errstate(invalid="ignore", divide="ignore"):
a = 1 / direction
t_min = (numpy.where(a >= 0, p0, p1) - origin) * a
t_max = (numpy.where(a >= 0, p1, p0) - origin) * a
t0 = numpy.max(t_min)
t1 = numpy.min(t_max)
return (t0 < t1) and (t1 > Ray.min_travel)

5
extensions/fablabchemnitz/raytracing/raytracing/material/__init__.py Normal file
View File

@ -0,0 +1,5 @@
from .optic_material import *
from .beamdump import *
from .mirror import *
from .beamsplitter import *
from .glass import *

15
extensions/fablabchemnitz/raytracing/raytracing/material/beamdump.py Normal file
View File

@ -0,0 +1,15 @@
from __future__ import annotations
from .optic_material import OpticMaterial
from ..ray import Ray
from ..shade import ShadeRec
class BeamDump(OpticMaterial):
"""Material absorbing all beams that hit it"""
def __repr__(self):
return "BeamDump()"
def generated_beams(self, ray: Ray, shade: ShadeRec) -> list[Ray]:
return list()

27
extensions/fablabchemnitz/raytracing/raytracing/material/beamsplitter.py Normal file
View File

@ -0,0 +1,27 @@
from typing import List
import numpy as np
from ..ray import Ray
from ..shade import ShadeRec
from .optic_material import OpticMaterial
class BeamSplitter(OpticMaterial):
"""
Material producing two beams after collision. One is reflected and
the other is transmitted.
"""
def __init__(self):
super().__init__()
def __repr__(self):
return "Mirror()"
def generated_beams(self, ray: Ray, shade: ShadeRec) -> List[Ray]:
o, d = shade.local_hit_point, ray.direction
n = shade.normal
reflected_ray = Ray(o, d - 2 * np.dot(d, n) * n)
transmitted_ray = Ray(o, d)
return [reflected_ray, transmitted_ray]

39
extensions/fablabchemnitz/raytracing/raytracing/material/glass.py Normal file
View File

@ -0,0 +1,39 @@
from typing import List
import numpy as np
from ..ray import Ray
from ..shade import ShadeRec
from .optic_material import OpticMaterial
class Glass(OpticMaterial):
"""Material that transmits and bends beams hitting it"""
def __init__(self, optical_index):
self._optical_index = optical_index
@property
def optical_index(self):
return self._optical_index
def __repr__(self):
return f"Glass({self._optical_index})"
def generated_beams(self, ray: Ray, shade: ShadeRec) -> List[Ray]:
o, d = shade.local_hit_point, ray.direction
n = shade.normal
if shade.hit_geometry.is_inside(ray):
n_1, n_2 = self.optical_index, 1
else:
n_1, n_2 = 1, self.optical_index
r = n_1 / n_2
c1 = -np.dot(d, n)
u = 1 - r ** 2 * (1 - c1 ** 2)
if u < 0: # total internal reflection
reflected_ray = Ray(o, d - 2 * np.dot(d, n) * n)
return [reflected_ray]
else: # refraction
c2 = np.sqrt(u)
transmitted_ray = Ray(o, r * d + (r * c1 - c2) * n)
return [transmitted_ray]

20
extensions/fablabchemnitz/raytracing/raytracing/material/mirror.py Normal file
View File

@ -0,0 +1,20 @@
from typing import List
import numpy
from ..ray import Ray
from ..shade import ShadeRec
from .optic_material import OpticMaterial
class Mirror(OpticMaterial):
"""Material reflecting beams that hit it"""
def __repr__(self):
return "Mirror()"
def generated_beams(self, ray: Ray, shade: ShadeRec) -> List[Ray]:
o, d = shade.local_hit_point, ray.direction
n = shade.normal
reflected_ray = Ray(o, d - 2 * numpy.dot(d, n) * n)
return [reflected_ray]

19
extensions/fablabchemnitz/raytracing/raytracing/material/optic_material.py Normal file
View File

@ -0,0 +1,19 @@
from abc import abstractmethod
from typing import Protocol, List
from ..ray import Ray
from ..shade import ShadeRec
class OpticMaterial(Protocol):
"""Protocol for an optical material"""
@abstractmethod
def generated_beams(self, ray: Ray, shade: ShadeRec) -> List[Ray]:
"""Compute the beams generated after intersection of a beam with this
material
Returns list of new beam seeds to start from after the intersection
of a beam and an object.
"""
raise NotImplementedError

10
extensions/fablabchemnitz/raytracing/raytracing/optical_object.py Normal file
View File

@ -0,0 +1,10 @@
from dataclasses import dataclass
from .geometry import GeometricObject
from .material import OpticMaterial
@dataclass(frozen=True)
class OpticalObject:
geometry: GeometricObject
material: OpticMaterial

18
extensions/fablabchemnitz/raytracing/raytracing/ray.py Normal file
View File

@ -0,0 +1,18 @@
from dataclasses import dataclass
from typing import ClassVar
from .vector import UnitVector, Vector
@dataclass(frozen=True)
class Ray:
"""This class implements a 2D line with an origin point and a direction."""
origin: Vector
direction: UnitVector
travel: float = 0
# If a beam hits an object before having traveled a minimum distance
# from its origin, the collision is ignored. This prevents infinite
# collision in case the origin of a beam is on the surface of an object
min_travel: ClassVar[float] = 1e-7

35
extensions/fablabchemnitz/raytracing/raytracing/shade.py Normal file
View File

@ -0,0 +1,35 @@
import numpy as np
from typing import Optional
class ShadeRec(object):
"""
This object contains the information needed to process the collision
between a ray and an object.
"""
def __init__(self):
self.hit_an_object: bool = False
self.local_hit_point: Optional[np.ndarray] = None
self.normal: Optional[np.ndarray] = None
self.travel_dist: float = np.inf
from .geometry import GeometricObject
self.hit_geometry: Optional[GeometricObject] = None
def __repr__(self):
return (
f"ShadeRec({self.hit_an_object}, {self.local_hit_point}, "
f"{self.normal}, {self.travel_dist})"
)
def set_normal_same_side(self, point: np.ndarray):
if self.normal is None:
raise RuntimeError("Can't find normal orientation if not already defined.")
elif self.local_hit_point is None:
raise RuntimeError(
"Can't find normal orientation if hit point not defined."
)
elif np.dot(self.normal, self.local_hit_point - point) > 0:
self.normal = -self.normal

75
extensions/fablabchemnitz/raytracing/raytracing/vector.py Normal file
View File

@ -0,0 +1,75 @@
from __future__ import annotations
from dataclasses import dataclass, field
from functools import singledispatchmethod
from math import sqrt
from numbers import Real
@dataclass(frozen=True)
class Vector:
x: float = field()
y: float = field()
def orthogonal(self) -> Vector:
"""Return a vector obtained by a pi/2 rotation"""
return UnitVector(-self.y, self.x)
@singledispatchmethod
def __mul__(self, other):
raise NotImplementedError
@__mul__.register
def _(self, other: Real):
return Vector(self.x * other, self.y * other)
@singledispatchmethod
def __rmul__(self, other):
raise NotImplementedError(type(other))
@__rmul__.register
def _(self, other: Real):
return Vector(self.x * other, self.y * other)
@singledispatchmethod
def __add__(self, other) -> Vector:
raise NotImplementedError
@singledispatchmethod
def __sub__(self, other) -> Vector:
raise NotImplementedError
def __neg__(self) -> Vector:
return Vector(-self.x, -self.y)
def norm(self):
return sqrt(self * self)
def normalize(self) -> UnitVector:
return UnitVector(self.x, self.y)
@dataclass(frozen=True)
class UnitVector(Vector):
def __init__(self, x, y):
norm = sqrt(x ** 2 + y ** 2)
super().__init__(x / norm, y / norm)
@Vector.__add__.register
def _(self, other: Vector):
return Vector(self.x + other.x, self.y + other.y)
@Vector.__sub__.register
def _(self, other: Vector):
return Vector(self.x - other.x, self.y - other.y)
@Vector.__mul__.register
def _(self, other: Vector) -> float:
return self.x * other.x + self.y * other.y
@Vector.__rmul__.register
def _(self, other: Vector) -> float:
return self.x * other.x + self.y * other.y

93
extensions/fablabchemnitz/raytracing/raytracing/world.py Normal file
View File

@ -0,0 +1,93 @@
"""
Module to describe and interact with a scene composed of various optical
objects
"""
from __future__ import annotations
import warnings
from dataclasses import dataclass, field
from typing import Optional, List, NamedTuple, Iterable, Tuple
from .geometry import GeometricObject
from .material import OpticMaterial, BeamDump
from .ray import Ray
from .shade import ShadeRec
class OpticalObject(NamedTuple):
geometry: GeometricObject
material: OpticMaterial
@dataclass
class World:
"""Stores a scene and computes the interaction with a ray"""
objects: Optional[list[OpticalObject]] = field(default_factory=list)
# default recursion depth can be changed, but should not exceed
# system recursion limit.
max_recursion_depth: Optional[int] = 500
def add(self, obj: OpticalObject):
self.objects.append(obj)
def __iter__(self) -> Iterable[OpticalObject]:
return iter(self.objects)
@property
def num_objects(self) -> int:
return len(self.objects)
def first_hit(self, ray: Ray) -> Tuple[ShadeRec, OpticMaterial]:
"""
Returns the information about the first collision of the beam
with an object.
:return: A shade for the collision geometric information and the
material of the object hit.
"""
result = ShadeRec()
material = BeamDump()
for obj in self:
shade = obj.geometry.hit(ray)
if Ray.min_travel < shade.travel_dist < result.travel_dist:
result = shade
material = obj.material
return result, material
def propagate_beams(self, seed):
return self._propagate_beams([[seed]], 0)
def _propagate_beams(self, beams: List[List[Ray]], depth) -> List[List[Ray]]:
"""Computes the propagation of beams in the system
:return: List of all the beam paths generated by these seeds.
It is stored as
[path0[Ray0, Ray1, ...], path1[...], ...].
Each path is a list of successive rays having each traveled a
given distance.
:raise: warning if recursion depth hits a limit.
"""
if depth >= self.max_recursion_depth:
err_msg = (
f"Maximal recursion depth exceeded ({self.max_recursion_depth})."
"It is likely that not all beams have been rendered."
)
warnings.warn(err_msg)
return beams
else:
new_beams = list()
for index, beam in enumerate(beams):
ray = beam[-1]
if ray.travel <= 0:
shade, material = self.first_hit(ray)
new_seeds = material.generated_beams(ray, shade)
beams[index][-1] = Ray(ray.origin, ray.direction, shade.travel_dist)
if len(new_seeds) == 0:
new_beams.append(beams[index])
for seed in new_seeds:
generated_beams = self._propagate_beams([[seed]], depth + 1)
for new_beam in generated_beams:
new_beams.append(beams[index] + new_beam)
return new_beams

17
extensions/fablabchemnitz/raytracing/render.inx Normal file
View File

@ -0,0 +1,17 @@
<?xml version="1.0" encoding="UTF-8"?>
<inkscape-extension xmlns="http://www.inkscape.org/namespace/inkscape/extension">
<name>Ray Tracing - Render</name>
<id>fablabchemnitz.de.raytracing.render</id>
<effect>
<object-type>all</object-type>
<effects-menu>
<submenu name="FabLab Chemnitz">
<submenu name="Various"/>
</submenu>
</effects-menu>
</effect>
<script>
<command location="inx" interpreter="python">render.py</command>
</script>
</inkscape-extension>

289
extensions/fablabchemnitz/raytracing/render.py Normal file
View File

@ -0,0 +1,289 @@
"""
Extension for rendering beams in 2D optics with Inkscape
"""
from __future__ import annotations
from dataclasses import dataclass
from functools import singledispatchmethod
from typing import Iterable, Optional, Final
import inkex
from inkex.paths import Line, Move
import raytracing.material
from desc_parser import get_optics_fields
from raytracing import Vector, World, OpticalObject, Ray
from raytracing.geometry import CubicBezier, CompoundGeometricObject, GeometricObject
from utils import pairwise
@dataclass
class BeamSeed:
ray: Optional[Ray] = None
parent: Optional[inkex.ShapeElement] = None
def get_unlinked_copy(clone: inkex.Use) -> Optional[inkex.ShapeElement]:
"""Creates a copy of the original with all transformations applied"""
copy = clone.href.copy()
copy.transform = clone.composed_transform() * copy.transform
copy.style = clone.style + copy.style
copy.getparent = clone.getparent
return copy
def get_or_create_beam_layer(parent_layer: inkex.Layer) -> inkex.Layer:
for child in parent_layer:
if isinstance(child, inkex.Layer):
if child.get("inkscape:label") == "generated_beams":
return child
new_layer = parent_layer.add(inkex.Layer())
new_layer.label = "generated_beams"
return new_layer
def plot_beam(beam: list[Ray], node: inkex.ShapeElement, layer: inkex.Layer):
path = inkex.Path()
if beam:
path += [Move(beam[0].origin.x, beam[0].origin.y)]
for ray in beam:
p1 = ray.origin + ray.travel * ray.direction
path += [Line(p1.x, p1.y)]
element = layer.add(inkex.PathElement())
# Need to convert to path to get the correct style for inkex.Use
element.style = node.to_path_element().style
element.path = path
class Raytracing(inkex.EffectExtension):
"""Extension to renders the beams present in the document"""
# Ray tracing is only implemented for the following inkex primitives
filter_primitives: Final = (
inkex.PathElement,
inkex.Line,
inkex.Polyline,
inkex.Polygon,
inkex.Rectangle,
inkex.Ellipse,
inkex.Circle,
)
def __init__(self):
super().__init__()
self.world = World()
self.beam_seeds: list[BeamSeed] = list()
def effect(self) -> None:
"""
Loads the objects and outputs a svg with the beams after propagation
"""
# Can't set the border earlier because self.svg is not yet defined
self.document_border = self.get_document_borders_as_beamdump()
self.world.add(self.document_border)
filter_ = self.filter_primitives + (inkex.Group, inkex.Use)
for obj in self.svg.selection.filter(filter_):
self.add(obj)
if self.beam_seeds:
for seed in self.beam_seeds:
if self.is_inside_document(seed.ray):
generated = self.world.propagate_beams(seed.ray)
for beam in generated:
try:
new_layer = get_or_create_beam_layer(
get_containing_layer(seed.parent)
)
plot_beam(beam, seed.parent, new_layer)
except LayerError as e:
inkex.utils.errormsg(f"{e} It will be ignored.")
@singledispatchmethod
def add(self, obj):
pass
@add.register
def _(self, group: inkex.Group):
for child in group:
self.add(child)
@add.register
def _(self, clone: inkex.Use):
copy = get_unlinked_copy(clone)
self.add(copy)
for type in filter_primitives:
@add.register(type)
def _(self, obj):
"""
Extracts properties and adds the object to the ray tracing data
structure
"""
material = get_material(obj)
if material:
if isinstance(material, BeamSeed):
for ray in get_beams(obj):
self.beam_seeds.append(BeamSeed(ray, parent=obj))
else:
geometry = get_geometry(obj)
opt_obj = OpticalObject(geometry, material)
self.world.add(opt_obj)
def get_document_borders_as_beamdump(self) -> OpticalObject:
"""
Adds a beam blocking contour on the borders of the document to
prevent the beams from going to infinity
"""
c1x, c1y, c2x, c2y = self.svg.get_viewbox()
contour_geometry = CompoundGeometricObject(
(
CubicBezier(
Vector(c1x, c1y),
Vector(c1x, c1y),
Vector(c2x, c1y),
Vector(c2x, c1y),
),
CubicBezier(
Vector(c2x, c1y),
Vector(c2x, c1y),
Vector(c2x, c2y),
Vector(c2x, c2y),
),
CubicBezier(
Vector(c2x, c2y),
Vector(c2x, c2y),
Vector(c1x, c2y),
Vector(c1x, c2y),
),
CubicBezier(
Vector(c1x, c2y),
Vector(c1x, c2y),
Vector(c1x, c1y),
Vector(c1x, c1y),
),
)
)
return OpticalObject(contour_geometry, raytracing.material.BeamDump())
def is_inside_document(self, ray: Ray) -> bool:
return self.document_border.geometry.is_inside(ray)
def get_material(
obj: inkex.ShapeElement,
) -> Optional[raytracing.material.OpticMaterial | BeamSeed]:
"""Extracts the optical material of an object from its description"""
desc = obj.desc
if desc is None:
desc = ""
materials = get_materials_from_description(desc)
if len(materials) == 0:
return None
if len(materials) > 1:
raise_err_num_materials(obj)
elif len(materials) == 1:
return materials[0]
def get_materials_from_description(
desc: str,
) -> list[raytracing.material.OpticMaterial | BeamSeed]:
"""Run through the description to extract the material properties"""
materials = list()
class_alias = dict(
beam_dump=raytracing.material.BeamDump,
mirror=raytracing.material.Mirror,
beam_splitter=raytracing.material.BeamSplitter,
glass=raytracing.material.Glass,
beam=BeamSeed,
)
for match in get_optics_fields(desc):
material_type = match.group("material")
prop_str = match.group("num")
if material_type in class_alias:
if material_type == "glass" and prop_str is not None:
optical_index = float(prop_str)
materials.append(class_alias[material_type](optical_index))
else:
materials.append(class_alias[material_type]())
return materials
def raise_err_num_materials(obj):
inkex.utils.errormsg(
f"The element {obj.get_id()} has more than one optical material and will be"
f" ignored:\n{obj.desc}\n"
)
def get_geometry(obj: inkex.ShapeElement) -> GeometricObject:
"""
Converts the geometry of inkscape elements to a form suitable for the
ray tracing module
"""
# Treats all objects as cubic Bezier curves. This treatment is exact
# for most primitives except circles and ellipses that are only
# approximated by Bezier curves.
# TODO: implement exact representation for ellipses
path = get_absolute_path(obj)
composite_bezier = convert_to_composite_bezier(path)
return composite_bezier
def get_absolute_path(obj: inkex.ShapeElement) -> inkex.CubicSuperPath:
path = obj.to_path_element().path.to_absolute()
transformed_path = path.transform(obj.composed_transform())
return transformed_path.to_superpath()
def get_beams(element: inkex.ShapeElement) -> Iterable[Ray]:
"""
Returns a beam with origin at the endpoint of the path and tangent to
the path
"""
bezier_path = convert_to_composite_bezier(get_absolute_path(element))
for sub_path in bezier_path:
last_segment = sub_path[-1]
endpoint = last_segment.eval(1)
tangent = -last_segment.tangent(1)
yield Ray(endpoint, tangent)
def convert_to_composite_bezier(
superpath: inkex.CubicSuperPath,
) -> CompoundGeometricObject:
"""
Converts a superpath with a representation
[Subpath0[handle0_0, point0, handle0_1], ...], ...]
to a representation of consecutive bezier segments of the form
CompositeCubicBezier([CubicBezierPath[CubicBezier[point0, handle0_1,
handle1_0, point1], ...], ...]).
"""
composite_bezier = list()
for subpath in superpath:
bezier_path = list()
for (__, p0, p1), (p2, p3, __) in pairwise(subpath):
bezier = CubicBezier(Vector(*p0), Vector(*p1), Vector(*p2), Vector(*p3))
bezier_path.append(bezier)
composite_bezier.append(CompoundGeometricObject(bezier_path))
return CompoundGeometricObject(composite_bezier)
def get_containing_layer(obj: inkex.BaseElement) -> inkex.Layer:
try:
return obj.ancestors().filter(inkex.Layer)[0]
except IndexError:
raise LayerError(f"Object '{obj.get_id()}' is not inside a layer.")
class LayerError(RuntimeError):
pass
if __name__ == "__main__":
Raytracing().run()

25
extensions/fablabchemnitz/raytracing/set_material.inx Normal file
View File

@ -0,0 +1,25 @@
<?xml version="1.0" encoding="UTF-8"?>
<inkscape-extension xmlns="http://www.inkscape.org/namespace/inkscape/extension">
<name>Ray Tracing - Set Lens Material</name>
<id>fablabchemnitz.de.raytracing.set_lens_material</id>
<param name="optical_material" type="optiongroup" appearance="combo" gui-text="Select material:">
<option value="none">None</option>
<option value="beam">Beam</option>
<option value="mirror">Mirror</option>
<option value="beam_dump">Beam dump</option>
<option value="beam_splitter">Beam splitter</option>
<option value="glass">Glass</option>
</param>
<param name="optical_index" type="float" min="1.0000" max="3.0000" precision="4" gui-text="Optical index:" indent="2">1.5168</param>
<effect>
<object-type>all</object-type>
<effects-menu>
<submenu name="FabLab Chemnitz">
<submenu name="Various"/>
</submenu>
</effects-menu>
</effect>
<script>
<command location="inx" interpreter="python">set_material.py</command>
</script>
</inkscape-extension>

67
extensions/fablabchemnitz/raytracing/set_material.py Normal file
View File

@ -0,0 +1,67 @@
from __future__ import annotations
from functools import singledispatchmethod
from typing import Final
import inkex
from desc_parser import clear_description
class SetMaterial(inkex.EffectExtension):
"""Writes the chosen optical property in the element description"""
# only change the description for these objects
filter_primitives: Final = (
inkex.PathElement,
inkex.Line,
inkex.Polyline,
inkex.Polygon,
inkex.Rectangle,
inkex.Ellipse,
inkex.Circle,
)
def __init__(self):
super().__init__()
def add_arguments(self, pars):
pars.add_argument("--optical_material", default="none", help="Name of the optical material to convert the selection to.")
pars.add_argument("--optical_index", type=float, default=1.5168)
def effect(self) -> None:
filter_ = self.filter_primitives + (inkex.Group,)
for obj in self.svg.selection.filter(filter_):
self.update_description(obj)
@singledispatchmethod
def update_description(self, arg):
pass
@update_description.register
def _(self, group: inkex.Group):
for obj in group:
self.update_description(obj)
for type in filter_primitives:
@update_description.register(type)
def _(self, obj):
"""
In the description of the element, replaces the optical properties
with the new one.
"""
desc = obj.desc
if desc is None:
desc = ""
new_desc = clear_description(desc)
if desc != "" and desc[-1] != "\n":
desc += "\n"
material_name = self.options.optical_material
if material_name is not None:
new_desc += f"optics:{material_name}"
if material_name == "glass":
new_desc += f":{self.options.optical_index:.4f}"
obj.desc = new_desc
if __name__ == "__main__":
SetMaterial().run()

4
extensions/fablabchemnitz/raytracing/test.py Normal file
View File

@ -0,0 +1,4 @@
from raytracing.geometry import CubicBezier
p = CubicBezier((1,1),(1,1),(1,1),(1,1))
print(p)

11
extensions/fablabchemnitz/raytracing/utils.py Normal file
View File

@ -0,0 +1,11 @@
import itertools
from typing import TypeVar, Iterator, Tuple
T = TypeVar("T")
def pairwise(iterable: Iterator[T]) -> Iterator[Tuple[T, T]]:
"""s -> (s0,s1), (s1,s2), (s2, s3), ..."""
a, b = itertools.tee(iterable)
next(b, None)
return zip(a, b)