diff --git a/ArchitectureColoredPainting/ArchitectureColoredPainting.vcxproj b/ArchitectureColoredPainting/ArchitectureColoredPainting.vcxproj
index 807742f..493a418 100644
--- a/ArchitectureColoredPainting/ArchitectureColoredPainting.vcxproj
+++ b/ArchitectureColoredPainting/ArchitectureColoredPainting.vcxproj
@@ -112,6 +112,7 @@
+
@@ -138,8 +139,6 @@
-
-
@@ -165,6 +164,7 @@
+
diff --git a/ArchitectureColoredPainting/ArchitectureColoredPainting.vcxproj.filters b/ArchitectureColoredPainting/ArchitectureColoredPainting.vcxproj.filters
index d687ce6..b28ae3f 100644
--- a/ArchitectureColoredPainting/ArchitectureColoredPainting.vcxproj.filters
+++ b/ArchitectureColoredPainting/ArchitectureColoredPainting.vcxproj.filters
@@ -144,6 +144,9 @@
Source Files\Renderer\Painting
+
+ Source Files\Renderer\Painting
+
@@ -229,10 +232,6 @@
Resource Files
-
- Source Files
-
-
@@ -309,6 +308,9 @@
Header Files\Renderer\Painting
+
+ Header Files\Renderer\Painting
+
diff --git a/ArchitectureColoredPainting/src/Renderer/Painting/CubicBezierSignedDistance.cpp b/ArchitectureColoredPainting/src/Renderer/Painting/CubicBezierSignedDistance.cpp
new file mode 100644
index 0000000..d7ed172
--- /dev/null
+++ b/ArchitectureColoredPainting/src/Renderer/Painting/CubicBezierSignedDistance.cpp
@@ -0,0 +1,522 @@
+#include "CubicBezierSignedDistance.h"
+#include
+#include
+
+using Renderer::CubicBezierSignedDistance;
+using std::min;
+using std::max;
+using glm::dvec2;
+using glm::dvec3;
+using glm::dvec4;
+using glm::dot;
+using glm::sign;
+using glm::clamp;
+
+const double CubicBezierSignedDistance::eps = .000005;
+const double CubicBezierSignedDistance::zoom = 1.;
+const double CubicBezierSignedDistance::dot_size = .005;
+const dvec3 CubicBezierSignedDistance::point_col = dvec3(1, 1, 0);
+int CubicBezierSignedDistance::halley_iterations = 8;
+
+double CubicBezierSignedDistance::cubic_bezier_dis(dvec2 uv, dvec2 p0, dvec2 p1, dvec2 p2, dvec2 p3) {
+
+ //switch points when near to end point to minimize numerical error
+ //only needed when control point(s) very far away
+#if 0
+ dvec2 mid_curve = parametric_cub_bezier(.5, p0, p1, p2, p3);
+ dvec2 mid_points = (p0 + p3) / 2.;
+
+ dvec2 tang = mid_curve - mid_points;
+ dvec2 nor = dvec2(tang.y, -tang.x);
+
+ if (sign(dot(nor, uv - mid_curve)) != sign(dot(nor, p0 - mid_curve))) {
+ dvec2 tmp = p0;
+ p0 = p3;
+ p3 = tmp;
+
+ tmp = p2;
+ p2 = p1;
+ p1 = tmp;
+ }
+#endif
+
+ dvec2 a3 = (-p0 + 3. * p1 - 3. * p2 + p3);
+ dvec2 a2 = (3. * p0 - 6. * p1 + 3. * p2);
+ dvec2 a1 = (-3. * p0 + 3. * p1);
+ dvec2 a0 = p0 - uv;
+
+ //compute polynomial describing distance to current pixel dependent on a parameter t
+ double bc6 = dot(a3, a3);
+ double bc5 = 2. * dot(a3, a2);
+ double bc4 = dot(a2, a2) + 2. * dot(a1, a3);
+ double bc3 = 2. * (dot(a1, a2) + dot(a0, a3));
+ double bc2 = dot(a1, a1) + 2. * dot(a0, a2);
+ double bc1 = 2. * dot(a0, a1);
+ double bc0 = dot(a0, a0);
+
+ bc5 /= bc6;
+ bc4 /= bc6;
+ bc3 /= bc6;
+ bc2 /= bc6;
+ bc1 /= bc6;
+ bc0 /= bc6;
+
+ //compute derivatives of this polynomial
+
+ double b0 = bc1 / 6.;
+ double b1 = 2. * bc2 / 6.;
+ double b2 = 3. * bc3 / 6.;
+ double b3 = 4. * bc4 / 6.;
+ double b4 = 5. * bc5 / 6.;
+
+ dvec4 c1 = dvec4(b1, 2. * b2, 3. * b3, 4. * b4) / 5.;
+ dvec3 c2 = dvec3(c1[1], 2. * c1[2], 3. * c1[3]) / 4.;
+ dvec2 c3 = dvec2(c2[1], 2. * c2[2]) / 3.;
+ double c4 = c3[1] / 2.;
+
+ dvec4 roots_drv = dvec4(1e38);
+
+ int num_roots_drv = solve_quartic(c1, roots_drv);
+ sort_roots4(roots_drv);
+
+ double ub = upper_bound_lagrange5(b0, b1, b2, b3, b4);
+ double lb = lower_bound_lagrange5(b0, b1, b2, b3, b4);
+
+ dvec3 a = dvec3(1e38);
+ dvec3 b = dvec3(1e38);
+
+ dvec3 roots = dvec3(1e38);
+
+ int num_roots = 0;
+
+ //compute root isolating intervals by roots of derivative and outer root bounds
+ //only roots going form - to + considered, because only those result in a minimum
+ if (num_roots_drv == 4) {
+ if (eval_poly5(b0, b1, b2, b3, b4, roots_drv[0]) > 0.) {
+ a[0] = lb;
+ b[0] = roots_drv[0];
+ num_roots = 1;
+ }
+
+ if (sign(eval_poly5(b0, b1, b2, b3, b4, roots_drv[1])) != sign(eval_poly5(b0, b1, b2, b3, b4, roots_drv[2]))) {
+ if (num_roots == 0) {
+ a[0] = roots_drv[1];
+ b[0] = roots_drv[2];
+ num_roots = 1;
+ }
+ else {
+ a[1] = roots_drv[1];
+ b[1] = roots_drv[2];
+ num_roots = 2;
+ }
+ }
+
+ if (eval_poly5(b0, b1, b2, b3, b4, roots_drv[3]) < 0.) {
+ if (num_roots == 0) {
+ a[0] = roots_drv[3];
+ b[0] = ub;
+ num_roots = 1;
+ }
+ else if (num_roots == 1) {
+ a[1] = roots_drv[3];
+ b[1] = ub;
+ num_roots = 2;
+ }
+ else {
+ a[2] = roots_drv[3];
+ b[2] = ub;
+ num_roots = 3;
+ }
+ }
+ }
+ else {
+ if (num_roots_drv == 2) {
+ if (eval_poly5(b0, b1, b2, b3, b4, roots_drv[0]) < 0.) {
+ num_roots = 1;
+ a[0] = roots_drv[1];
+ b[0] = ub;
+ }
+ else if (eval_poly5(b0, b1, b2, b3, b4, roots_drv[1]) > 0.) {
+ num_roots = 1;
+ a[0] = lb;
+ b[0] = roots_drv[0];
+ }
+ else {
+ num_roots = 2;
+
+ a[0] = lb;
+ b[0] = roots_drv[0];
+
+ a[1] = roots_drv[1];
+ b[1] = ub;
+ }
+
+ }
+ else {//num_roots_drv==0
+ dvec3 roots_snd_drv = dvec3(1e38);
+ int num_roots_snd_drv = solve_cubic(c2, roots_snd_drv);
+
+ dvec2 roots_trd_drv = dvec2(1e38);
+ int num_roots_trd_drv = solve_quadric(c3, roots_trd_drv);
+ num_roots = 1;
+
+ a[0] = lb;
+ b[0] = ub;
+ }
+
+ //further subdivide intervals to guarantee convergence of halley's method
+ //by using roots of further derivatives
+ dvec3 roots_snd_drv = dvec3(1e38);
+ int num_roots_snd_drv = solve_cubic(c2, roots_snd_drv);
+ sort_roots3(roots_snd_drv);
+
+ int num_roots_trd_drv = 0;
+ dvec2 roots_trd_drv = dvec2(1e38);
+
+ if (num_roots_snd_drv != 3) {
+ num_roots_trd_drv = solve_quadric(c3, roots_trd_drv);
+ }
+
+ for (int i = 0; i < 3; i++) {
+ if (i < num_roots) {
+ for (int j = 0; j < 3; j += 2) {
+ if (j < num_roots_snd_drv) {
+ if (a[i] < roots_snd_drv[j] && b[i] > roots_snd_drv[j]) {
+ if (eval_poly5(b0, b1, b2, b3, b4, roots_snd_drv[j]) > 0.) {
+ b[i] = roots_snd_drv[j];
+ }
+ else {
+ a[i] = roots_snd_drv[j];
+ }
+ }
+ }
+ }
+ for (int j = 0; j < 2; j++) {
+ if (j < num_roots_trd_drv) {
+ if (a[i] < roots_trd_drv[j] && b[i] > roots_trd_drv[j]) {
+ if (eval_poly5(b0, b1, b2, b3, b4, roots_trd_drv[j]) > 0.) {
+ b[i] = roots_trd_drv[j];
+ }
+ else {
+ a[i] = roots_trd_drv[j];
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ double d0 = 1e38;
+
+ //compute roots with halley's method
+
+ for (int i = 0; i < 3; i++) {
+ if (i < num_roots) {
+ roots[i] = .5 * (a[i] + b[i]);
+
+ for (int j = 0; j < halley_iterations; j++) {
+ roots[i] = halley_iteration5(b0, b1, b2, b3, b4, roots[i]);
+ }
+
+
+ //compute squared distance to nearest point on curve
+ roots[i] = clamp(roots[i], 0., 1.);
+ dvec2 to_curve = uv - parametric_cub_bezier(roots[i], p0, p1, p2, p3);
+ d0 = min(d0, dot(to_curve, to_curve));
+ }
+ }
+
+ return sqrt(d0);
+}
+
+dvec2 CubicBezierSignedDistance::parametric_cub_bezier(double t, dvec2 p0, dvec2 p1, dvec2 p2, dvec2 p3) {
+ dvec2 a0 = (-p0 + 3. * p1 - 3. * p2 + p3);
+ dvec2 a1 = (3. * p0 - 6. * p1 + 3. * p2);
+ dvec2 a2 = (-3. * p0 + 3. * p1);
+ dvec2 a3 = p0;
+
+ return (((a0 * t) + a1) * t + a2) * t + a3;
+}
+
+int CubicBezierSignedDistance::solve_quartic(dvec4 coeffs, dvec4& s) {
+
+ double a = coeffs[3];
+ double b = coeffs[2];
+ double c = coeffs[1];
+ double d = coeffs[0];
+
+ /* substitute x = y - A/4 to eliminate cubic term:
+ x^4 + px^2 + qx + r = 0 */
+
+ double sq_a = a * a;
+ double p = -3. / 8. * sq_a + b;
+ double q = 1. / 8. * sq_a * a - 1. / 2. * a * b + c;
+ double r = -3. / 256. * sq_a * sq_a + 1. / 16. * sq_a * b - 1. / 4. * a * c + d;
+
+ int num;
+
+ /* doesn't seem to happen for me */
+ //if(abs(r) -eps) {
+ u = sqrt(abs(u));
+ }
+ else {
+ return 0;
+ }
+
+ if (v > -eps) {
+ v = sqrt(abs(v));
+ }
+ else {
+ return 0;
+ }
+
+ dvec2 quad_coeffs, Sxy(s.x, s.y);
+
+ quad_coeffs[0] = z - u;
+ quad_coeffs[1] = q < 0. ? -v : v;
+
+ num = solve_quadric(quad_coeffs, Sxy);
+ s.x = Sxy.x, s.y = Sxy.y;
+
+ quad_coeffs[0] = z + u;
+ quad_coeffs[1] = q < 0. ? v : -v;
+
+ dvec2 tmp = dvec2(1e38);
+ int old_num = num;
+
+ num += solve_quadric(quad_coeffs, tmp);
+ if (old_num != num) {
+ if (old_num == 0) {
+ s[0] = tmp[0];
+ s[1] = tmp[1];
+ }
+ else {//old_num == 2
+ s[2] = tmp[0];
+ s[3] = tmp[1];
+ }
+ }
+ }
+
+ /* resubstitute */
+
+ double sub = 1. / 4. * a;
+
+ /* single halley iteration to fix cancellation */
+ for (int i = 0; i < 4; i += 2) {
+ if (i < num) {
+ s[i] -= sub;
+ s[i] = halley_iteration4(coeffs, s[i]);
+
+ s[i + 1] -= sub;
+ s[i + 1] = halley_iteration4(coeffs, s[i + 1]);
+ }
+ }
+
+ return num;
+}
+
+int CubicBezierSignedDistance::solve_cubic(dvec3 coeffs, dvec3& r) {
+
+ double a = coeffs[2];
+ double b = coeffs[1];
+ double c = coeffs[0];
+
+ double p = b - a * a / 3.0;
+ double q = a * (2.0 * a * a - 9.0 * b) / 27.0 + c;
+ double p3 = p * p * p;
+ double d = q * q + 4.0 * p3 / 27.0;
+ double offset = -a / 3.0;
+ if (d >= 0.0) { // Single solution
+ double z = sqrt(d);
+ double u = (-q + z) / 2.0;
+ double v = (-q - z) / 2.0;
+ u = sign(u) * pow(abs(u), 1.0 / 3.0);
+ v = sign(v) * pow(abs(v), 1.0 / 3.0);
+ r[0] = offset + u + v;
+
+ //Single newton iteration to account for cancellation
+ double f = ((r[0] + a) * r[0] + b) * r[0] + c;
+ double f1 = (3. * r[0] + 2. * a) * r[0] + b;
+
+ r[0] -= f / f1;
+
+ return 1;
+ }
+ double u = sqrt(-p / 3.0);
+ double v = acos(-sqrt(-27.0 / p3) * q / 2.0) / 3.0;
+ double m = cos(v), n = sin(v) * 1.732050808;
+
+ //Single newton iteration to account for cancellation
+ //(once for every root)
+ r[0] = offset + u * (m + m);
+ r[1] = offset - u * (n + m);
+ r[2] = offset + u * (n - m);
+
+ dvec3 f = ((r + a) * r + b) * r + c;
+ dvec3 f1 = (3. * r + 2. * a) * r + b;
+
+ r -= f / f1;
+
+ return 3;
+}
+
+int CubicBezierSignedDistance::solve_quadric(dvec2 coeffs, dvec2& roots) {
+
+ // normal form: x^2 + px + q = 0
+ double p = coeffs[1] / 2.;
+ double q = coeffs[0];
+
+ double D = p * p - q;
+
+ if (D < 0.) {
+ return 0;
+ }
+ else if (D > 0.) {
+ roots[0] = -sqrt(D) - p;
+ roots[1] = sqrt(D) - p;
+
+ return 2;
+ }
+}
+
+void CubicBezierSignedDistance::sort_roots3(dvec3& roots) {
+ dvec3 tmp;
+
+ tmp[0] = min(roots[0], min(roots[1], roots[2]));
+ tmp[1] = max(roots[0], min(roots[1], roots[2]));
+ tmp[2] = max(roots[0], max(roots[1], roots[2]));
+
+ roots = tmp;
+}
+
+void CubicBezierSignedDistance::sort_roots4(dvec4& roots) {
+ dvec4 tmp;
+
+ dvec2 min1_2 = min(dvec2(roots.x, roots.z), dvec2(roots.y, roots.w));
+ dvec2 max1_2 = max(dvec2(roots.x, roots.z), dvec2(roots.y, roots.w));
+
+ double maxmin = max(min1_2.x, min1_2.y);
+ double minmax = min(max1_2.x, max1_2.y);
+
+ tmp[0] = min(min1_2.x, min1_2.y);
+ tmp[1] = min(maxmin, minmax);
+ tmp[2] = max(minmax, maxmin);
+ tmp[3] = max(max1_2.x, max1_2.y);
+
+ roots = tmp;
+}
+
+double CubicBezierSignedDistance::lower_bound_lagrange5(double a0, double a1, double a2, double a3, double a4) {
+
+ dvec4 coeffs1 = dvec4(-a0, a1, -a2, a3);
+
+ dvec4 neg1 = max(-coeffs1, dvec4(0));
+ double neg2 = max(-a4, 0.);
+
+ const dvec4 indizes1 = dvec4(0, 1, 2, 3);
+ const double indizes2 = 4.;
+
+ dvec4 bounds1 = pow(neg1, 1. / (5. - indizes1));
+ double bounds2 = pow(neg2, 1. / (5. - indizes2));
+
+ dvec2 min1_2 = min(dvec2(bounds1.x, bounds1.z), dvec2(bounds1.y, bounds1.w));
+ dvec2 max1_2 = max(dvec2(bounds1.x, bounds1.z), dvec2(bounds1.y, bounds1.w));
+
+ double maxmin = max(min1_2.x, min1_2.y);
+ double minmax = min(max1_2.x, max1_2.y);
+
+ double max3 = max(max1_2.x, max1_2.y);
+
+ double max_max = max(max3, bounds2);
+ double max_max2 = max(min(max3, bounds2), max(minmax, maxmin));
+
+ return -max_max - max_max2;
+}
+
+double CubicBezierSignedDistance::upper_bound_lagrange5(double a0, double a1, double a2, double a3, double a4) {
+
+ dvec4 coeffs1 = dvec4(a0, a1, a2, a3);
+
+ dvec4 neg1 = max(-coeffs1, dvec4(0));
+ double neg2 = max(-a4, 0.);
+
+ const dvec4 indizes1 = dvec4(0, 1, 2, 3);
+ const double indizes2 = 4.;
+
+ dvec4 bounds1 = pow(neg1, 1. / (5. - indizes1));
+ double bounds2 = pow(neg2, 1. / (5. - indizes2));
+
+ dvec2 min1_2 = min(dvec2(bounds1.x, bounds1.z), dvec2(bounds1.y, bounds1.w));
+ dvec2 max1_2 = max(dvec2(bounds1.x, bounds1.z), dvec2(bounds1.y, bounds1.w));
+
+ double maxmin = max(min1_2.x, min1_2.y);
+ double minmax = min(max1_2.x, max1_2.y);
+
+ double max3 = max(max1_2.x, max1_2.y);
+
+ double max_max = max(max3, bounds2);
+ double max_max2 = max(min(max3, bounds2), max(minmax, maxmin));
+
+ return max_max + max_max2;
+}
+
+double CubicBezierSignedDistance::halley_iteration4(dvec4 coeffs, double x) {
+
+ double f = (((x + coeffs[3]) * x + coeffs[2]) * x + coeffs[1]) * x + coeffs[0];
+ double f1 = ((4. * x + 3. * coeffs[3]) * x + 2. * coeffs[2]) * x + coeffs[1];
+ double f2 = (12. * x + 6. * coeffs[3]) * x + 2. * coeffs[2];
+
+ return x - (2. * f * f1) / (2. * f1 * f1 - f * f2);
+}
+
+double CubicBezierSignedDistance::halley_iteration5(double a0, double a1, double a2, double a3, double a4, double x) {
+
+ double f = ((((x + a4) * x + a3) * x + a2) * x + a1) * x + a0;
+ double f1 = (((5. * x + 4. * a4) * x + 3. * a3) * x + 2. * a2) * x + a1;
+ double f2 = ((20. * x + 12. * a4) * x + 6. * a3) * x + 2. * a2;
+
+ return x - (2. * f * f1) / (2. * f1 * f1 - f * f2);
+}
+
+double CubicBezierSignedDistance::eval_poly5(double a0, double a1, double a2, double a3, double a4, double x) {
+
+ double f = ((((x + a4) * x + a3) * x + a2) * x + a1) * x + a0;
+
+ return f;
+}
\ No newline at end of file
diff --git a/ArchitectureColoredPainting/src/Renderer/Painting/CubicBezierSignedDistance.h b/ArchitectureColoredPainting/src/Renderer/Painting/CubicBezierSignedDistance.h
new file mode 100644
index 0000000..ef3c0da
--- /dev/null
+++ b/ArchitectureColoredPainting/src/Renderer/Painting/CubicBezierSignedDistance.h
@@ -0,0 +1,36 @@
+#pragma once
+#include
+#include
+#include
+
+using glm::dvec2;
+using glm::dvec3;
+using glm::dvec4;
+
+namespace Renderer {
+ class CubicBezierSignedDistance
+ {
+ private:
+ static const double eps;
+ static const double zoom;
+ static const double dot_size;
+ static const dvec3 point_col;
+ static int halley_iterations;
+
+ public:
+ static int solve_cubic(dvec3 coeffs, dvec3& r);
+ static int solve_quartic(dvec4 coeffs, dvec4& s);
+ static int solve_quadric(dvec2 coeffs, dvec2& roots);
+ static void sort_roots3(dvec3& roots);
+ static void sort_roots4(dvec4& roots);
+ static double lower_bound_lagrange5(double a0, double a1, double a2, double a3, double a4);
+ static double upper_bound_lagrange5(double a0, double a1, double a2, double a3, double a4);
+ static double halley_iteration4(dvec4 coeffs, double x);
+ static double halley_iteration5(double a0, double a1, double a2, double a3, double a4, double x);
+ static dvec2 parametric_cub_bezier(double t, dvec2 p0, dvec2 p1, dvec2 p2, dvec2 p3);
+ static double eval_poly5(double a0, double a1, double a2, double a3, double a4, double x);
+ static double cubic_bezier_dis(dvec2 uv, dvec2 p0, dvec2 p1, dvec2 p2, dvec2 p3);
+
+ };
+}
+
diff --git a/ArchitectureColoredPainting/src/Renderer/Painting/Line.h b/ArchitectureColoredPainting/src/Renderer/Painting/Line.h
index c0677cb..9fa4618 100644
--- a/ArchitectureColoredPainting/src/Renderer/Painting/Line.h
+++ b/ArchitectureColoredPainting/src/Renderer/Painting/Line.h
@@ -6,6 +6,7 @@
#include
#include
#include
+
namespace Renderer
{
@@ -46,7 +47,7 @@ namespace Renderer
virtual void monotonization(std::vector & res) = 0;
virtual int judgeBoundIntersection(double xy, double l, double r, bool isY) = 0;
virtual bool judgeIntersection(QVector4D bound);
- virtual double minDistanceToLine(Point uv);
+ // virtual double minDistanceToLine(Point uv);
void upwardOffsetForY(bool isStart, double offset);
bool isLineContained(QVector4D bound);
Point operator[](int index) const;