增加了求点到贝塞尔曲线最小距离的工具类

dev-VirtualTexture
yang.yongquan 2022-11-21 16:01:27 +08:00
parent 598e0ce7a9
commit 2e2d89b80c
5 changed files with 568 additions and 7 deletions

View File

@ -112,6 +112,7 @@
<ClCompile Include="src\Renderer\Painting\BvhTree.cpp" />
<ClCompile Include="src\Renderer\Camera.cpp" />
<ClCompile Include="src\Renderer\Painting\CubicBezier.cpp" />
<ClCompile Include="src\Renderer\Painting\CubicBezierSignedDistance.cpp" />
<ClCompile Include="src\Renderer\Painting\CubicMonotonization.cpp" />
<ClCompile Include="src\Renderer\Light.cpp" />
<ClCompile Include="src\Renderer\Painting\Element.cpp" />
@ -138,8 +139,6 @@
<ItemGroup>
<None Include="darkstyle.qss" />
<None Include="lightstyle.qss" />
<None Include="msvc_make.bat" />
<None Include="qt.conf" />
<None Include="Shaders\depth_init.comp" />
<None Include="Shaders\depth_mipmap.comp" />
<None Include="Shaders\final.frag" />
@ -165,6 +164,7 @@
<ClInclude Include="src\Editor\LayerStyle.h" />
<ClInclude Include="src\Editor\LayerWrapper.h" />
<QtMoc Include="src\Editor\PreviewWindow.h" />
<ClInclude Include="src\Renderer\Painting\CubicBezierSignedDistance.h" />
<ClInclude Include="src\Renderer\Painting\Element.h" />
<ClInclude Include="src\Renderer\Painting\LineTree.h" />
<ClInclude Include="src\Renderer\Painting\Painting.h" />

View File

@ -144,6 +144,9 @@
<ClCompile Include="src\Renderer\Painting\LineTree.cpp">
<Filter>Source Files\Renderer\Painting</Filter>
</ClCompile>
<ClCompile Include="src\Renderer\Painting\CubicBezierSignedDistance.cpp">
<Filter>Source Files\Renderer\Painting</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<QtMoc Include="src\Renderer\RendererGLWidget.h">
@ -229,10 +232,6 @@
<None Include="lightstyle.qss">
<Filter>Resource Files</Filter>
</None>
<None Include="msvc_make.bat">
<Filter>Source Files</Filter>
</None>
<None Include="qt.conf" />
</ItemGroup>
<ItemGroup>
<QtUic Include="EditorWidget.ui">
@ -309,6 +308,9 @@
<ClInclude Include="src\Renderer\Painting\LineTree.h">
<Filter>Header Files\Renderer\Painting</Filter>
</ClInclude>
<ClInclude Include="src\Renderer\Painting\CubicBezierSignedDistance.h">
<Filter>Header Files\Renderer\Painting</Filter>
</ClInclude>
</ItemGroup>
<ItemGroup>
<QtRcc Include="MainWindow.qrc">

View File

@ -0,0 +1,522 @@
#include "CubicBezierSignedDistance.h"
#include <glm/gtc/matrix_transform.hpp>
#include <algorithm>
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){
// /* no absolute term: y(y^3 + py + q) = 0 */
// dvec3 cubic_coeffs;
// cubic_coeffs[0] = q;
// cubic_coeffs[1] = p;
// cubic_coeffs[2] = 0.;
// num = solve_cubic(cubic_coeffs, s.xyz);
// s[num] = 0.;
// num++;
//}
{
/* solve the resolvent cubic ... */
dvec3 cubic_coeffs, Sxyz(s.x, s.y, s.z);
cubic_coeffs[0] = 1.0 / 2. * r * p - 1.0 / 8. * q * q;
cubic_coeffs[1] = -r;
cubic_coeffs[2] = -1.0 / 2. * p;
solve_cubic(cubic_coeffs, Sxyz);
s.x = Sxyz.x, s.y = Sxyz.y, s.z = Sxyz.z;
/* ... and take the one real solution ... */
double z = s[0];
/* ... to build two quadric equations */
double u = z * z - r;
double v = 2. * z - p;
if (u > -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;
}

View File

@ -0,0 +1,36 @@
#pragma once
#include <glm/vec2.hpp>
#include <glm/vec3.hpp>
#include <glm/vec4.hpp>
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);
};
}

View File

@ -6,6 +6,7 @@
#include <algorithm>
#include <memory>
#include <iostream>
namespace Renderer
{
@ -46,7 +47,7 @@ namespace Renderer
virtual void monotonization(std::vector <LinePtr>& 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;