diff --git a/libs/geom/include/psemek/geom/intersection.hpp b/libs/geom/include/psemek/geom/intersection.hpp
index 3a697e8f..49d9dc94 100644
--- a/libs/geom/include/psemek/geom/intersection.hpp
+++ b/libs/geom/include/psemek/geom/intersection.hpp
@@ -6,6 +6,7 @@
 #include <psemek/geom/ray.hpp>
 #include <psemek/geom/box.hpp>
 #include <psemek/geom/sphere.hpp>
+#include <psemek/geom/cylinder.hpp>
 #include <psemek/geom/orientation.hpp>
 #include <psemek/geom/contains.hpp>
 #include <psemek/geom/gauss.hpp>
@@ -222,10 +223,42 @@ namespace psemek::geom
 		return {ir->first, ir->second};
 	}
 
+	template <typename T, std::size_t N>
+	interval<T> intersection(ray<T, N> const & r, cylinder<T, N> const & c)
+	{
+		auto Z2 = length_sqr(c.axis);
+		auto CZ = dot(c.center - r.origin, c.axis);
+		auto DZ = dot(r.direction, c.axis);
+		if (DZ == 0)
+			return {};
+
+		auto prmin = (-Z2 + CZ) / DZ;
+		auto prmax = ( Z2 + CZ) / DZ;
+		if (DZ < 0)
+			std::swap(prmin, prmax);
+
+		auto n = normalized(c.axis);
+		auto q = r.origin - c.center;
+		auto ir = solve_quadratic(length_sqr(r.direction) - sqr(dot(r.direction, n)),
+			2 * dot(q, r.direction) - 2 * dot(q, n) * dot(r.direction, n),
+			length_sqr(q) - sqr(dot(q, n)) - sqr(c.radius));
+
+		if (!ir)
+			return {};
+
+		return {std::max(ir->first, prmin), std::min(ir->second, prmax)};
+	}
+
 	template <typename T, std::size_t N>
 	bool intersect(ray<T, N> const & r, sphere<T, N> const & s)
 	{
 		return intersection(r, s).empty();
 	}
 
+	template <typename T, std::size_t N>
+	bool intersect(ray<T, N> const & r, cylinder<T, N> const & c)
+	{
+		return intersection(r, c).empty();
+	}
+
 }