summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libdimension/Makefile.am1
-rw-r--r--libdimension/cylinder.c34
-rw-r--r--libdimension/dimension.h1
-rw-r--r--libdimension/dimension/geometry.h7
-rw-r--r--libdimension/dimension/polynomial.h51
-rw-r--r--libdimension/polynomial.c313
-rw-r--r--libdimension/sphere.c36
7 files changed, 411 insertions, 32 deletions
diff --git a/libdimension/Makefile.am b/libdimension/Makefile.am
index 36517ee..d444682 100644
--- a/libdimension/Makefile.am
+++ b/libdimension/Makefile.am
@@ -74,6 +74,7 @@ libdimension_la_SOURCES = $(nobase_include_HEADERS) \
platform.c \
platform.h \
point_light.c \
+ polynomial.c \
progress.c \
progress-impl.h \
prtree.c \
diff --git a/libdimension/cylinder.c b/libdimension/cylinder.c
index e9db50c..ba016cf 100644
--- a/libdimension/cylinder.c
+++ b/libdimension/cylinder.c
@@ -69,19 +69,27 @@ dmnsn_cylinder_intersection_fn(const dmnsn_object *cylinder, dmnsn_line line,
/* Solve (x0 + nx*t)^2 + (z0 + nz*t)^2
== (((r2 - r1)*(y0 + ny*t) + r1 + r2)/2)^2 */
- double a, b, c, t;
- a = l.n.x*l.n.x + l.n.z*l.n.z - l.n.y*l.n.y*(r2 - r1)*(r2 - r1)/4.0;
- b = 2.0*(l.n.x*l.x0.x + l.n.z*l.x0.z)
- - l.n.y*(r2 - r1)*(l.x0.y*(r2 - r1) + r2 + r1)/2.0;
- c = l.x0.x*l.x0.x + l.x0.z*l.x0.z
- - (l.x0.y*(r2 - r1) + r2 + r1)*(l.x0.y*(r2 - r1) + r2 + r1)/4;
-
- if (b*b - 4.0*a*c >= 0.0) {
- t = (-b - sqrt(b*b - 4.0*a*c))/(2.0*a);
- dmnsn_vector p = dmnsn_line_point(l, t);
-
- if (t < 0.0 || p.y <= -1.0 || p.y >= 1.0) {
- t = (-b + sqrt(b*b - 4.0*a*c))/(2.0*a);
+ double poly[3], x[2];
+ poly[2] = l.n.x*l.n.x + l.n.z*l.n.z - l.n.y*l.n.y*(r2 - r1)*(r2 - r1)/4.0;
+ poly[1] = 2.0*(l.n.x*l.x0.x + l.n.z*l.x0.z)
+ - l.n.y*(r2 - r1)*(l.x0.y*(r2 - r1) + r2 + r1)/2.0;
+ poly[0] = l.x0.x*l.x0.x + l.x0.z*l.x0.z
+ - (l.x0.y*(r2 - r1) + r2 + r1)*(l.x0.y*(r2 - r1) + r2 + r1)/4;
+
+ size_t n = dmnsn_solve_polynomial(poly, 2, x);
+
+ if (n > 0) {
+ double t = x[0];
+ dmnsn_vector p;
+ if (n == 2) {
+ t = dmnsn_min(t, x[1]);
+ p = dmnsn_line_point(l, t);
+
+ if (p.y <= -1.0 || p.y >= 1.0) {
+ t = dmnsn_max(x[0], x[1]);
+ p = dmnsn_line_point(l, t);
+ }
+ } else {
p = dmnsn_line_point(l, t);
}
diff --git a/libdimension/dimension.h b/libdimension/dimension.h
index c58eb0e..5cdc5b2 100644
--- a/libdimension/dimension.h
+++ b/libdimension/dimension.h
@@ -67,6 +67,7 @@ typedef void dmnsn_free_fn(void *ptr);
#include <dimension/progress.h>
#include <dimension/timer.h>
#include <dimension/geometry.h>
+#include <dimension/polynomial.h>
#include <dimension/color.h>
#include <dimension/canvas.h>
#include <dimension/gl.h>
diff --git a/libdimension/dimension/geometry.h b/libdimension/dimension/geometry.h
index 87e7ab7..4c4f8b8 100644
--- a/libdimension/dimension/geometry.h
+++ b/libdimension/dimension/geometry.h
@@ -96,6 +96,13 @@ dmnsn_degrees(double radians)
return radians*45.0/atan(1.0);
}
+DMNSN_INLINE int
+dmnsn_signbit(double n)
+{
+ /* Guarantee a 1 or 0 return, to allow testing two signs for equality */
+ return signbit(n) ? 1 : 0;
+}
+
/* Shorthand for vector/matrix construction */
DMNSN_INLINE dmnsn_vector
diff --git a/libdimension/dimension/polynomial.h b/libdimension/dimension/polynomial.h
new file mode 100644
index 0000000..52a44dd
--- /dev/null
+++ b/libdimension/dimension/polynomial.h
@@ -0,0 +1,51 @@
+/*************************************************************************
+ * Copyright (C) 2009-2010 Tavian Barnes <tavianator@gmail.com> *
+ * *
+ * This file is part of The Dimension Library. *
+ * *
+ * The Dimension Library is free software; you can redistribute it and/ *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version. *
+ * *
+ * The Dimension Library is distributed in the hope that it will be *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
+ * Lesser General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU Lesser General Public *
+ * License along with this program. If not, see *
+ * <http://www.gnu.org/licenses/>. *
+ *************************************************************************/
+
+/*
+ * Utility functions for working with and numerically solving polynomials.
+ * Polynomials are represented as simple arrays where the ith element is the
+ * coefficient on x^i. In general, we are only interested in positive roots.
+ */
+
+#ifndef DIMENSION_POLYNOMIAL_H
+#define DIMENSION_POLYNOMIAL_H
+
+#include <stddef.h>
+#include <stdio.h>
+
+DMNSN_INLINE double
+dmnsn_evaluate_polynomial(double poly[], size_t degree, double x)
+{
+ double ret = poly[degree];
+ ssize_t i;
+ for (i = degree - 1; i >= 0; --i) {
+ ret = ret*x + poly[i];
+ }
+ return ret;
+}
+
+/* Stores the non-negative roots of poly[] in x[], and returns the number of
+ such roots that were stored */
+size_t dmnsn_solve_polynomial(double poly[], size_t degree, double x[]);
+
+/* Helper function to print a polynomial */
+void dmnsn_print_polynomial(FILE *file, double poly[], size_t degree);
+
+#endif /* DIMENSION_POLYNOMIAL_H */
diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c
new file mode 100644
index 0000000..41313bb
--- /dev/null
+++ b/libdimension/polynomial.c
@@ -0,0 +1,313 @@
+/*************************************************************************
+ * Copyright (C) 2009-2010 Tavian Barnes <tavianator@gmail.com> *
+ * *
+ * This file is part of The Dimension Library. *
+ * *
+ * The Dimension Library is free software; you can redistribute it and/ *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version. *
+ * *
+ * The Dimension Library is distributed in the hope that it will be *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
+ * Lesser General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU Lesser General Public *
+ * License along with this program. If not, see *
+ * <http://www.gnu.org/licenses/>. *
+ *************************************************************************/
+
+#include "dimension.h"
+#include <math.h>
+
+/* Basic solving methods */
+
+static inline size_t
+dmnsn_solve_linear(double poly[2], double x[1])
+{
+ x[0] = -poly[0]/poly[1];
+ return (x[0] >= 0.0) ? 1 : 0;
+}
+
+static inline size_t
+dmnsn_solve_quadratic(double poly[3], double x[2])
+{
+ double disc = poly[1]*poly[1] - 4.0*poly[0]*poly[2];
+ if (disc >= 0.0) {
+ double s = sqrt(disc);
+ x[0] = (-poly[1] + s)/(2.0*poly[2]);
+ x[1] = (-poly[1] - s)/(2.0*poly[2]);
+ return (x[0] >= 0.0) ? ((x[1] >= 0.0) ? 2 : 1) : 0;
+ } else {
+ return 0;
+ }
+}
+
+static inline size_t
+dmnsn_zero_roots(double poly[], size_t *degree, double x[])
+{
+ size_t i = 0;
+ while (i <= *degree && poly[i] == 0.0) {
+ x[i] = 0.0;
+ ++i;
+ }
+
+ if (i > 0) {
+ for (size_t j = 0; j + i <= *degree; ++j) {
+ poly[j] = poly[j + i];
+ }
+ }
+
+ *degree -= i;
+ return i;
+}
+
+/* Get the real degree of a polynomial, ignoring leading zeros */
+static inline size_t
+dmnsn_real_degree(double poly[], size_t degree)
+{
+ for (ssize_t i = degree; i >= 0; ++i) {
+ if (poly[i] != 0.0) {
+ return i;
+ }
+ }
+
+ return 0;
+}
+
+/* Improve a root with Newton's method */
+static inline double
+dmnsn_improve_root(double poly[], size_t degree, double x)
+{
+ double error;
+ do {
+ /* Calculate the value of the polynomial and its derrivative at once */
+ double p = poly[degree], dp = 0.0;
+ for (ssize_t i = degree - 1; i >= 0; --i) {
+ dp = dp*x + p;
+ p = p*x + poly[i];
+ }
+
+ double dx = p/dp;
+ error = fabs(dx/x);
+ x -= dx;
+ } while (error > dmnsn_epsilon);
+
+ return x;
+}
+
+/* Use the method of bisection to find a root in a range that contains exactly
+ one root, counting multiplicity */
+static inline double
+dmnsn_bisect_root(double poly[], size_t degree, double min, double max)
+{
+ double evmin = dmnsn_evaluate_polynomial(poly, degree, min);
+ double evmax = dmnsn_evaluate_polynomial(poly, degree, max);
+
+ while (max - min > dmnsn_epsilon) {
+ double mid = (min + max)/2.0;
+ double evmid = dmnsn_evaluate_polynomial(poly, degree, mid);
+ if (dmnsn_signbit(evmid) == dmnsn_signbit(evmax)) {
+ max = mid;
+ evmax = evmid;
+ } else {
+ min = mid;
+ evmin = evmid;
+ }
+ }
+
+ return (min + max)/2.0;
+}
+
+/* Use synthetic division to eliminate the root `r' from poly[] */
+static inline void
+dmnsn_eliminate_root(double poly[], size_t *degree, double r)
+{
+ double rem = poly[*degree];
+ for (ssize_t i = *degree - 1; i >= 0; --i) {
+ double temp = poly[i];
+ poly[i] = rem;
+ rem = temp + r*rem;
+ }
+ --*degree;
+}
+
+/* Returns the number of sign changes between coefficients of `poly' */
+static inline size_t
+dmnsn_descartes_rule(double poly[], size_t degree)
+{
+ int lastsign = 0;
+ size_t i;
+ for (i = 0; i <= degree; ++i) {
+ if (poly[i] != 0.0) {
+ lastsign = dmnsn_signbit(poly[i]);
+ break;
+ }
+ }
+
+ size_t changes = 0;
+ for (++i; i <= degree; ++i) {
+ int sign = dmnsn_signbit(poly[i]);
+ if (poly[i] != 0.0 && sign != lastsign) {
+ lastsign = sign;
+ ++changes;
+ }
+ }
+
+ return changes;
+}
+
+/* Get the (n k) binomial coefficient */
+static double
+dmnsn_binom(size_t n, size_t k)
+{
+ if (k > n - k) {
+ k = n - k;
+ }
+
+ double ret = 1.0;
+ for (size_t i = 0; i < k; ++i) {
+ ret *= n - i;
+ ret /= i + 1;
+ }
+
+ return ret;
+}
+
+/* Find a range that contains a single root, with Uspensky's algorithm */
+static bool
+dmnsn_uspensky_bound(double poly[], size_t degree, double *min, double *max)
+{
+ size_t n = dmnsn_descartes_rule(poly, degree);
+ if (n == 0) {
+ return false;
+ } else if (n == 1) {
+ return true;
+ } else {
+ double a[degree];
+ /* a is the expanded form of poly(x + 1) */
+ for (size_t i = 0; i <= degree; ++i) {
+ a[i] = poly[i];
+ for (size_t j = i + 1; j <= degree; ++j) {
+ a[i] += dmnsn_binom(j, i)*poly[j];
+ }
+ }
+
+ if (a[0] == 0.0) {
+ *max = *min;
+ return true;
+ } else if (dmnsn_uspensky_bound(a, degree, min, max)) {
+ *min += 1.0;
+ *max += 1.0;
+ return true;
+ }
+
+ double b[degree];
+ /* b is the expanded form of (x + 1)^degree * poly(1/(x + 1)) */
+ for (size_t i = 0; i <= degree; ++i) {
+ b[i] = poly[degree - i];
+ for (size_t j = i + 1; j <= degree; ++j) {
+ b[i] += dmnsn_binom(j, i)*poly[degree - j];
+ }
+ }
+
+ if (dmnsn_uspensky_bound(b, degree, min, max)) {
+ double temp = *min;
+ *min = 1.0/(*max + 1.0);
+ *max = 1.0/(temp + 1.0);
+ return true;
+ }
+
+ return false;
+ }
+}
+
+/* Modified Uspensky's algorithm */
+size_t
+dmnsn_solve_polynomial(double poly[], size_t degree, double x[])
+{
+ /* Copy the polynomial so we can be destructive */
+ double p[degree + 1];
+ for (ssize_t i = degree; i >= 0; --i) {
+ p[i] = poly[i];
+ }
+
+ size_t i = 0; /* Index into x[] */
+
+ /* Eliminate simple zero roots */
+ i += dmnsn_zero_roots(p, &degree, x + i);
+ /* Account for leading zero coefficients */
+ degree = dmnsn_real_degree(p, degree);
+
+ while (degree > 2) {
+ /* Get a bound on the range of positive roots */
+ double min = +0.0, max = INFINITY;
+ if (dmnsn_uspensky_bound(p, degree, &min, &max)) {
+ if (isinf(max)) {
+ /* Replace an infinite upper bound with a finite one due to Cauchy */
+ max = 0.0;
+ for (size_t j = 0; j < degree; ++j) {
+ max = dmnsn_max(max, fabs(p[j]));
+ }
+ max /= fabs(p[degree]);
+ max += 1.0;
+ }
+
+ /* Bisect within the found range */
+ double r = dmnsn_bisect_root(p, degree, min, max);
+ r = dmnsn_improve_root(p, degree, r);
+
+ /* Use synthetic division to eliminate the root `r' */
+ dmnsn_eliminate_root(p, &degree, r);
+
+ /* Store the found root */
+ x[i] = r;
+ ++i;
+ } else {
+ break;
+ }
+
+ i += dmnsn_zero_roots(p, &degree, x + i);
+ degree = dmnsn_real_degree(p, degree);
+ }
+
+ switch (degree) {
+ case 1:
+ i += dmnsn_solve_linear(p, x + i);
+ break;
+
+ case 2:
+ i += dmnsn_solve_quadratic(p, x + i);
+ break;
+ }
+
+ return i;
+}
+
+void
+dmnsn_print_polynomial(FILE *file, double poly[], size_t degree)
+{
+ fprintf(file, "%g*x^%zu", poly[degree], degree);
+ for (size_t i = degree - 1; i > 1; --i) {
+ if (poly[i] > 0.0) {
+ fprintf(file, " + %g*x^%zu", poly[i], i);
+ } else if (poly[i] < 0.0) {
+ fprintf(file, " - %g*x^%zu", -poly[i], i);
+ }
+ }
+
+ if (poly[1] > 0.0) {
+ fprintf(file, " + %g*x", poly[1]);
+ } else if (poly[1] < 0.0) {
+ fprintf(file, " - %g*x", -poly[1]);
+ }
+
+ if (poly[0] > 0.0) {
+ fprintf(file, " + %g", poly[0]);
+ } else if (poly[0] < 0.0) {
+ fprintf(file, " - %g", -poly[0]);
+ }
+
+ fprintf(file, "\n");
+}
diff --git a/libdimension/sphere.c b/libdimension/sphere.c
index 98be4b5..490074b 100644
--- a/libdimension/sphere.c
+++ b/libdimension/sphere.c
@@ -53,29 +53,27 @@ dmnsn_sphere_intersection_fn(const dmnsn_object *sphere, dmnsn_line line,
dmnsn_line l = dmnsn_transform_line(sphere->trans_inv, line);
/* Solve (x0 + nx*t)^2 + (y0 + ny*t)^2 + (z0 + nz*t)^2 == 1 */
- double a, b, c, t;
- a = dmnsn_vector_dot(l.n, l.n);
- b = 2.0*dmnsn_vector_dot(l.n, l.x0);
- c = dmnsn_vector_dot(l.x0, l.x0) - 1.0;
+ double poly[3], x[2];
+ poly[2] = dmnsn_vector_dot(l.n, l.n);
+ poly[1] = 2.0*dmnsn_vector_dot(l.n, l.x0);
+ poly[0] = dmnsn_vector_dot(l.x0, l.x0) - 1.0;
- if (b*b - 4.0*a*c >= 0.0) {
- t = (-b - sqrt(b*b - 4.0*a*c))/(2.0*a);
- if (t < 0.0) {
- t = (-b + sqrt(b*b - 4.0*a*c))/(2.0*a);
- }
+ size_t n = dmnsn_solve_polynomial(poly, 2, x);
+ if (n == 0) {
+ return false;
+ } else {
+ double t = x[0];
+ if (n == 2)
+ t = dmnsn_min(t, x[1]);
- if (t >= 0.0) {
- intersection->ray = line;
- intersection->t = t;
- intersection->normal = dmnsn_transform_normal(sphere->trans,
+ intersection->ray = line;
+ intersection->t = t;
+ intersection->normal = dmnsn_transform_normal(sphere->trans,
dmnsn_line_point(l, t));
- intersection->texture = sphere->texture;
- intersection->interior = sphere->interior;
- return true;
- }
+ intersection->texture = sphere->texture;
+ intersection->interior = sphere->interior;
+ return true;
}
-
- return false;
}
/* Return whether a point is inside a sphere (x**2 + y**2 + z**2 < 1.0) */
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-23 13:53:16 +0000