From 940b4d1848e8c70ab7642901a68594e8016caffc Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Sat, 27 Apr 2024 18:51:28 +0200 Subject: Adding upstream version 1:7.0.4. Signed-off-by: Daniel Baumann --- basegfx/source/workbench/Makefile | 34 + basegfx/source/workbench/bezierclip.cxx | 2004 +++++++++++++++++++++++++++++++ basegfx/source/workbench/bezierclip.hxx | 87 ++ basegfx/source/workbench/convexhull.cxx | 200 +++ basegfx/source/workbench/gauss.hxx | 170 +++ 5 files changed, 2495 insertions(+) create mode 100644 basegfx/source/workbench/Makefile create mode 100644 basegfx/source/workbench/bezierclip.cxx create mode 100644 basegfx/source/workbench/bezierclip.hxx create mode 100644 basegfx/source/workbench/convexhull.cxx create mode 100644 basegfx/source/workbench/gauss.hxx (limited to 'basegfx/source/workbench') diff --git a/basegfx/source/workbench/Makefile b/basegfx/source/workbench/Makefile new file mode 100644 index 000000000..21dfc1400 --- /dev/null +++ b/basegfx/source/workbench/Makefile @@ -0,0 +1,34 @@ +# +# This file is part of the LibreOffice project. +# +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. +# +# This file incorporates work covered by the following license notice: +# +# Licensed to the Apache Software Foundation (ASF) under one or more +# contributor license agreements. See the NOTICE file distributed +# with this work for additional information regarding copyright +# ownership. The ASF licenses this file to you under the Apache +# License, Version 2.0 (the "License"); you may not use this file +# except in compliance with the License. You may obtain a copy of +# the License at http://www.apache.org/licenses/LICENSE-2.0 . +# + +# Testbuild + +#test : bezierclip.cxx convexhull.cxx +# g++ -Wall -g \ +# -I. -I. -I../inc -I./inc -I./unx/inc -I./unxlngi4/inc -I. -I/develop4/update/SRX644/unxlngi4/inc.m4/stl -I/develop4/update/SRX644/unxlngi4/inc.m4/external -I/develop4/update/SRX644/unxlngi4/inc.m4 -I/develop4/update/SRX644/src.m4/solenv/unxlngi4/inc -I/net/grande/develop6/update/dev/gcc_3.0.1_linux_libc2.11_turbolinux/include -I/develop4/update/SRX644/src.m4/solenv/inc -I/develop4/update/SRX644/unxlngi4/inc.m4/stl -I/net/grande.germany/develop6/update/dev/gcc_3.0.1_linux_libc2.11_turbolinux/redhat60/usr/include -I/net/grande.germany/develop6/update/dev/gcc_3.0.1_linux_libc2.11_turbolinux/redhat60/usr/include/X11 -I/develop4/update/SRX644/src.m4/res -I/net/grande/develop6/update/dev/Linux_JDK_1.4.0/include -I/net/grande/develop6/update/dev/Linux_JDK_1.4.0/include/linux -I. -I./res -I. \ +# -include preinclude.h -D_USE_NAMESPACE -DGLIBC=2 -D_USE_NAMESPACE=1 -D_DEBUG_RUNTIME \ +# bezierclip.cxx convexhull.cxx -o bezierclip + +prog : bezierclip.cxx convexhull.cxx + g++ -Wall -g bezierclip.cxx convexhull.cxx -o bezierclip + +test : testconvexhull.cxx + g++ -Wall -g testconvexhull.cxx -o testhull + +.cxx.o: + g++ -c $(LOCALDEFINES) $(CCFLAGS) $< diff --git a/basegfx/source/workbench/bezierclip.cxx b/basegfx/source/workbench/bezierclip.cxx new file mode 100644 index 000000000..1f16ed37c --- /dev/null +++ b/basegfx/source/workbench/bezierclip.cxx @@ -0,0 +1,2004 @@ +/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ +/* + * This file is part of the LibreOffice project. + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * This file incorporates work covered by the following license notice: + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed + * with this work for additional information regarding copyright + * ownership. The ASF licenses this file to you under the Apache + * License, Version 2.0 (the "License"); you may not use this file + * except in compliance with the License. You may obtain a copy of + * the License at http://www.apache.org/licenses/LICENSE-2.0 . + */ + +#include +#include +#include +#include + +#include + +#include +#include + +// what to test +#define WITH_ASSERTIONS +//#define WITH_CONVEXHULL_TEST +//#define WITH_MULTISUBDIVIDE_TEST +//#define WITH_FATLINE_TEST +//#define WITH_CALCFOCUS_TEST +//#define WITH_SAFEPARAMBASE_TEST +//#define WITH_SAFEPARAMS_TEST +//#define WITH_SAFEPARAM_DETAILED_TEST +//#define WITH_SAFEFOCUSPARAM_CALCFOCUS +//#define WITH_SAFEFOCUSPARAM_TEST +//#define WITH_SAFEFOCUSPARAM_DETAILED_TEST +#define WITH_BEZIERCLIP_TEST + +/* Implementation of the so-called 'Fat-Line Bezier Clipping Algorithm' by Sederberg et al. + * + * Actual reference is: T. W. Sederberg and T Nishita: Curve + * intersection using Bezier clipping. In Computer Aided Design, 22 + * (9), 1990, pp. 538--549 + */ + +/* Misc helper + * =========== + */ +int fallFac( int n, int k ) +{ +#ifdef WITH_ASSERTIONS + assert(n>=k); // "For factorials, n must be greater or equal k" + assert(n>=0); // "For factorials, n must be positive" + assert(k>=0); // "For factorials, k must be positive" +#endif + + int res( 1 ); + + while( k-- && n ) res *= n--; + + return res; +} + +int fac( int n ) +{ + return fallFac(n, n); +} + +/* Bezier fat line clipping part + * ============================= + */ + +void Impl_calcFatLine( FatLine& line, const Bezier& c ) +{ + // Prepare normalized implicit line + // ================================ + + // calculate vector orthogonal to p1-p4: + line.a = -(c.p0.y - c.p3.y); + line.b = (c.p0.x - c.p3.x); + + // normalize + const double len( sqrt( line.a*line.a + line.b*line.b ) ); + if( !tolZero(len) ) + { + line.a /= len; + line.b /= len; + } + + line.c = -(line.a*c.p0.x + line.b*c.p0.y); + + // Determine bounding fat line from it + // =================================== + + // calc control point distances + const double dP2( calcLineDistance(line.a, line.b, line.c, c.p1.x, c.p1.y ) ); + const double dP3( calcLineDistance(line.a, line.b, line.c, c.p2.x, c.p2.y ) ); + + // calc approximate bounding lines to curve (tight bounds are + // possible here, but more expensive to calculate and thus not + // worth the overhead) + if( dP2 * dP3 > 0.0 ) + { + line.dMin = 3.0/4.0 * std::min(0.0, std::min(dP2, dP3)); + line.dMax = 3.0/4.0 * std::max(0.0, std::max(dP2, dP3)); + } + else + { + line.dMin = 4.0/9.0 * std::min(0.0, std::min(dP2, dP3)); + line.dMax = 4.0/9.0 * std::max(0.0, std::max(dP2, dP3)); + } +} + +void Impl_calcBounds( Point2D& leftTop, + Point2D& rightBottom, + const Bezier& c1 ) +{ + leftTop.x = std::min( c1.p0.x, std::min( c1.p1.x, std::min( c1.p2.x, c1.p3.x ) ) ); + leftTop.y = std::min( c1.p0.y, std::min( c1.p1.y, std::min( c1.p2.y, c1.p3.y ) ) ); + rightBottom.x = std::max( c1.p0.x, std::max( c1.p1.x, std::max( c1.p2.x, c1.p3.x ) ) ); + rightBottom.y = std::max( c1.p0.y, std::max( c1.p1.y, std::max( c1.p2.y, c1.p3.y ) ) ); +} + +bool Impl_doBBoxIntersect( const Bezier& c1, + const Bezier& c2 ) +{ + // calc rectangular boxes from c1 and c2 + Point2D lt1; + Point2D rb1; + Point2D lt2; + Point2D rb2; + + Impl_calcBounds( lt1, rb1, c1 ); + Impl_calcBounds( lt2, rb2, c2 ); + + if( std::min(rb1.x, rb2.x) < std::max(lt1.x, lt2.x) || + std::min(rb1.y, rb2.y) < std::max(lt1.y, lt2.y) ) + { + return false; + } + else + { + return true; + } +} + +/* calculates two t's for the given bernstein control polygon: the first is + * the intersection of the min value line with the convex hull from + * the left, the second is the intersection of the max value line with + * the convex hull from the right. + */ +bool Impl_calcSafeParams( double& t1, + double& t2, + const Polygon2D& rPoly, + double lowerYBound, + double upperYBound ) +{ + // need the convex hull of the control polygon, as this is + // guaranteed to completely bound the curve + Polygon2D convHull( convexHull(rPoly) ); + + // init min and max buffers + t1 = 0.0 ; + double currLowerT( 1.0 ); + + t2 = 1.0; + double currHigherT( 0.0 ); + + if( convHull.size() <= 1 ) + return false; // only one point? Then we're done with clipping + + /* now, clip against lower and higher bounds */ + Point2D p0; + Point2D p1; + + bool bIntersection( false ); + + for( Polygon2D::size_type i=0; i part2 is input curve, part1 is empty (input.p0, that is) + part1.p0.x = part1.p1.x = part1.p2.x = part1.p3.x = input.p0.x; + part1.p0.y = part1.p1.y = part1.p2.y = part1.p3.y = input.p0.y; + part2 = input; + } + else if( tolEqual(t, 1.0) ) + { + // t is one -> part1 is input curve, part2 is empty (input.p3, that is) + part1 = input; + part2.p0.x = part2.p1.x = part2.p2.x = part2.p3.x = input.p3.x; + part2.p0.y = part2.p1.y = part2.p2.y = part2.p3.y = input.p3.y; + } + else + { + part1.p0.x = input.p0.x; part1.p0.y = input.p0.y; + part1.p1.x = (1.0 - t)*part1.p0.x + t*input.p1.x; part1.p1.y = (1.0 - t)*part1.p0.y + t*input.p1.y; + const double Hx ( (1.0 - t)*input.p1.x + t*input.p2.x ), Hy ( (1.0 - t)*input.p1.y + t*input.p2.y ); + part1.p2.x = (1.0 - t)*part1.p1.x + t*Hx; part1.p2.y = (1.0 - t)*part1.p1.y + t*Hy; + part2.p3.x = input.p3.x; part2.p3.y = input.p3.y; + part2.p2.x = (1.0 - t)*input.p2.x + t*input.p3.x; part2.p2.y = (1.0 - t)*input.p2.y + t*input.p3.y; + part2.p1.x = (1.0 - t)*Hx + t*part2.p2.x; part2.p1.y = (1.0 - t)*Hy + t*part2.p2.y; + part2.p0.x = (1.0 - t)*part1.p2.x + t*part2.p1.x; part2.p0.y = (1.0 - t)*part1.p2.y + t*part2.p1.y; + part1.p3.x = part2.p0.x; part1.p3.y = part2.p0.y; + } +} + +void printCurvesWithSafeRange( const Bezier& c1, const Bezier& c2, double t1_c1, double t2_c1, + const Bezier& c2_part, const FatLine& bounds_c2 ) +{ + static int offset = 0; + + cout << "# safe param range testing" << endl + << "plot [t=0.0:1.0] "; + + // clip safe ranges off c1 + Bezier c1_part1; + Bezier c1_part2; + Bezier c1_part3; + + // subdivide at t1_c1 + Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1_c1 ); + // subdivide at t2_c1 + Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, t2_c1 ); + + // output remaining segment (c1_part1) + + cout << "bez(" + << c1.p0.x+offset << "," + << c1.p1.x+offset << "," + << c1.p2.x+offset << "," + << c1.p3.x+offset << ",t),bez(" + << c1.p0.y << "," + << c1.p1.y << "," + << c1.p2.y << "," + << c1.p3.y << ",t), bez(" + << c2.p0.x+offset << "," + << c2.p1.x+offset << "," + << c2.p2.x+offset << "," + << c2.p3.x+offset << ",t),bez(" + << c2.p0.y << "," + << c2.p1.y << "," + << c2.p2.y << "," + << c2.p3.y << ",t), " +#if 1 + << "bez(" + << c1_part1.p0.x+offset << "," + << c1_part1.p1.x+offset << "," + << c1_part1.p2.x+offset << "," + << c1_part1.p3.x+offset << ",t),bez(" + << c1_part1.p0.y << "," + << c1_part1.p1.y << "," + << c1_part1.p2.y << "," + << c1_part1.p3.y << ",t), " +#endif +#if 1 + << "bez(" + << c2_part.p0.x+offset << "," + << c2_part.p1.x+offset << "," + << c2_part.p2.x+offset << "," + << c2_part.p3.x+offset << ",t),bez(" + << c2_part.p0.y << "," + << c2_part.p1.y << "," + << c2_part.p2.y << "," + << c2_part.p3.y << ",t), " +#endif + << "linex(" + << bounds_c2.a << "," + << bounds_c2.b << "," + << bounds_c2.c << ",t)+" << offset << ", liney(" + << bounds_c2.a << "," + << bounds_c2.b << "," + << bounds_c2.c << ",t) title \"fat line (center)\", linex(" + << bounds_c2.a << "," + << bounds_c2.b << "," + << bounds_c2.c-bounds_c2.dMin << ",t)+" << offset << ", liney(" + << bounds_c2.a << "," + << bounds_c2.b << "," + << bounds_c2.c-bounds_c2.dMin << ",t) title \"fat line (min) \", linex(" + << bounds_c2.a << "," + << bounds_c2.b << "," + << bounds_c2.c-bounds_c2.dMax << ",t)+" << offset << ", liney(" + << bounds_c2.a << "," + << bounds_c2.b << "," + << bounds_c2.c-bounds_c2.dMax << ",t) title \"fat line (max) \"" << endl; + + offset += 1; +} + +void printResultWithFinalCurves( const Bezier& c1, const Bezier& c1_part, + const Bezier& c2, const Bezier& c2_part, + double t1_c1, double t2_c1 ) +{ + static int offset = 0; + + cout << "# final result" << endl + << "plot [t=0.0:1.0] "; + + cout << "bez(" + << c1.p0.x+offset << "," + << c1.p1.x+offset << "," + << c1.p2.x+offset << "," + << c1.p3.x+offset << ",t),bez(" + << c1.p0.y << "," + << c1.p1.y << "," + << c1.p2.y << "," + << c1.p3.y << ",t), bez(" + << c1_part.p0.x+offset << "," + << c1_part.p1.x+offset << "," + << c1_part.p2.x+offset << "," + << c1_part.p3.x+offset << ",t),bez(" + << c1_part.p0.y << "," + << c1_part.p1.y << "," + << c1_part.p2.y << "," + << c1_part.p3.y << ",t), " + << " pointmarkx(bez(" + << c1.p0.x+offset << "," + << c1.p1.x+offset << "," + << c1.p2.x+offset << "," + << c1.p3.x+offset << "," + << t1_c1 << "),t), " + << " pointmarky(bez(" + << c1.p0.y << "," + << c1.p1.y << "," + << c1.p2.y << "," + << c1.p3.y << "," + << t1_c1 << "),t), " + << " pointmarkx(bez(" + << c1.p0.x+offset << "," + << c1.p1.x+offset << "," + << c1.p2.x+offset << "," + << c1.p3.x+offset << "," + << t2_c1 << "),t), " + << " pointmarky(bez(" + << c1.p0.y << "," + << c1.p1.y << "," + << c1.p2.y << "," + << c1.p3.y << "," + << t2_c1 << "),t), " + + << "bez(" + << c2.p0.x+offset << "," + << c2.p1.x+offset << "," + << c2.p2.x+offset << "," + << c2.p3.x+offset << ",t),bez(" + << c2.p0.y << "," + << c2.p1.y << "," + << c2.p2.y << "," + << c2.p3.y << ",t), " + << "bez(" + << c2_part.p0.x+offset << "," + << c2_part.p1.x+offset << "," + << c2_part.p2.x+offset << "," + << c2_part.p3.x+offset << ",t),bez(" + << c2_part.p0.y << "," + << c2_part.p1.y << "," + << c2_part.p2.y << "," + << c2_part.p3.y << ",t)" << endl; + + offset += 1; +} + +/** determine parameter ranges [0,t1) and (t2,1] on c1, where c1 is guaranteed to lie outside c2. + Returns false, if the two curves don't even intersect. + + @param t1 + Range [0,t1) on c1 is guaranteed to lie outside c2 + + @param t2 + Range (t2,1] on c1 is guaranteed to lie outside c2 + + @param c1_orig + Original curve c1 + + @param c1_part + Subdivided current part of c1 + + @param c2_orig + Original curve c2 + + @param c2_part + Subdivided current part of c2 + */ +bool Impl_calcClipRange( double& t1, + double& t2, + const Bezier& c1_orig, + const Bezier& c1_part, + const Bezier& c2_orig, + const Bezier& c2_part ) +{ + // TODO: Maybe also check fat line orthogonal to P0P3, having P0 + // and P3 as the extremal points + + if( Impl_doBBoxIntersect(c1_part, c2_part) ) + { + // Calculate fat lines around c1 + FatLine bounds_c2; + + // must use the subdivided version of c2, since the fat line + // algorithm works implicitly with the convex hull bounding + // box. + Impl_calcFatLine(bounds_c2, c2_part); + + // determine clip positions on c2. Can use original c1 (which + // is necessary anyway, to get the t's on the original curve), + // since the distance calculations work directly in the + // Bernstein polynomial parameter domain. + if( Impl_calcSafeParams_clip( t1, t2, bounds_c2, + calcLineDistance( bounds_c2.a, + bounds_c2.b, + bounds_c2.c, + c1_orig.p0.x, + c1_orig.p0.y ), + calcLineDistance( bounds_c2.a, + bounds_c2.b, + bounds_c2.c, + c1_orig.p1.x, + c1_orig.p1.y ), + calcLineDistance( bounds_c2.a, + bounds_c2.b, + bounds_c2.c, + c1_orig.p2.x, + c1_orig.p2.y ), + calcLineDistance( bounds_c2.a, + bounds_c2.b, + bounds_c2.c, + c1_orig.p3.x, + c1_orig.p3.y ) ) ) + { + //printCurvesWithSafeRange(c1_orig, c2_orig, t1, t2, c2_part, bounds_c2); + + // they do intersect + return true; + } + } + + // they don't intersect: nothing to do + return false; +} + +/* Tangent intersection part + * ========================= + */ + +void Impl_calcFocus( Bezier& res, const Bezier& c ) +{ + // arbitrary small value, for now + // TODO: find meaningful value + const double minPivotValue( 1.0e-20 ); + + Point2D::value_type fMatrix[6]; + Point2D::value_type fRes[2]; + + // calc new curve from hodograph, c and linear blend + + // Coefficients for derivative of c are (C_i=n(C_{i+1} - C_i)): + + // 3(P1 - P0), 3(P2 - P1), 3(P3 - P2) (bezier curve of degree 2) + + // The hodograph is then (bezier curve of 2nd degree is P0(1-t)^2 + 2P1(1-t)t + P2t^2): + + // 3(P1 - P0)(1-t)^2 + 6(P2 - P1)(1-t)t + 3(P3 - P2)t^2 + + // rotate by 90 degrees: x=-y, y=x and you get the normal vector function N(t): + + // x(t) = -(3(P1.y - P0.y)(1-t)^2 + 6(P2.y - P1.y)(1-t)t + 3(P3.y - P2.y)t^2) + // y(t) = 3(P1.x - P0.x)(1-t)^2 + 6(P2.x - P1.x)(1-t)t + 3(P3.x - P2.x)t^2 + + // Now, the focus curve is defined to be F(t)=P(t) + c(t)N(t), + // where P(t) is the original curve, and c(t)=c0(1-t) + c1 t + + // This results in the following expression for F(t): + + // x(t) = P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 - + // (c0(1-t) + c1 t)(3(P1.y - P0.y)(1-t)^2 + 6(P2.y - P1.y)(1-t)t + 3(P3.y - P2.y)t^2) + + // y(t) = P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1.t)t^2 + P3.y t^3 + + // (c0(1-t) + c1 t)(3(P1.x - P0.x)(1-t)^2 + 6(P2.x - P1.x)(1-t)t + 3(P3.x - P2.x)t^2) + + // As a heuristic, we set F(0)=F(1) (thus, the curve is closed and _tends_ to be small): + + // For F(0), the following results: + + // x(0) = P0.x - c0 3(P1.y - P0.y) + // y(0) = P0.y + c0 3(P1.x - P0.x) + + // For F(1), the following results: + + // x(1) = P3.x - c1 3(P3.y - P2.y) + // y(1) = P3.y + c1 3(P3.x - P2.x) + + // Reorder, collect and substitute into F(0)=F(1): + + // P0.x - c0 3(P1.y - P0.y) = P3.x - c1 3(P3.y - P2.y) + // P0.y + c0 3(P1.x - P0.x) = P3.y + c1 3(P3.x - P2.x) + + // which yields + + // (P0.y - P1.y)c0 + (P3.y - P2.y)c1 = (P3.x - P0.x)/3 + // (P1.x - P0.x)c0 + (P2.x - P3.x)c1 = (P3.y - P0.y)/3 + + // so, this is what we calculate here (determine c0 and c1): + fMatrix[0] = c.p1.x - c.p0.x; + fMatrix[1] = c.p2.x - c.p3.x; + fMatrix[2] = (c.p3.y - c.p0.y)/3.0; + fMatrix[3] = c.p0.y - c.p1.y; + fMatrix[4] = c.p3.y - c.p2.y; + fMatrix[5] = (c.p3.x - c.p0.x)/3.0; + + // TODO: determine meaningful value for + if( !solve(fMatrix, 2, 3, fRes, minPivotValue) ) + { + // TODO: generate meaningful values here + // singular or nearly singular system -- use arbitrary + // values for res + fRes[0] = 0.0; + fRes[1] = 1.0; + + cerr << "Matrix singular!" << endl; + } + + // now, the reordered and per-coefficient collected focus curve is + // the following third degree bezier curve F(t): + + // x(t) = P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 - + // (c0(1-t) + c1 t)(3(P1.y - P0.y)(1-t)^2 + 6(P2.y - P1.y)(1-t)t + 3(P3.y - P2.y)t^2) + // = P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 - + // (3c0P1.y(1-t)^3 - 3c0P0.y(1-t)^3 + 6c0P2.y(1-t)^2t - 6c0P1.y(1-t)^2t + + // 3c0P3.y(1-t)t^2 - 3c0P2.y(1-t)t^2 + + // 3c1P1.y(1-t)^2t - 3c1P0.y(1-t)^2t + 6c1P2.y(1-t)t^2 - 6c1P1.y(1-t)t^2 + + // 3c1P3.yt^3 - 3c1P2.yt^3) + // = (P0.x - 3 c0 P1.y + 3 c0 P0.y)(1-t)^3 + + // 3(P1.x - c1 P1.y + c1 P0.y - 2 c0 P2.y + 2 c0 P1.y)(1-t)^2t + + // 3(P2.x - 2 c1 P2.y + 2 c1 P1.y - c0 P3.y + c0 P2.y)(1-t)t^2 + + // (P3.x - 3 c1 P3.y + 3 c1 P2.y)t^3 + // = (P0.x - 3 c0(P1.y - P0.y))(1-t)^3 + + // 3(P1.x - c1(P1.y - P0.y) - 2c0(P2.y - P1.y))(1-t)^2t + + // 3(P2.x - 2 c1(P2.y - P1.y) - c0(P3.y - P2.y))(1-t)t^2 + + // (P3.x - 3 c1(P3.y - P2.y))t^3 + + // y(t) = P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1-t)t^2 + P3.y t^3 + + // (c0(1-t) + c1 t)(3(P1.x - P0.x)(1-t)^2 + 6(P2.x - P1.x)(1-t)t + 3(P3.x - P2.x)t^2) + // = P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1-t)t^2 + P3.y t^3 + + // 3c0(P1.x - P0.x)(1-t)^3 + 6c0(P2.x - P1.x)(1-t)^2t + 3c0(P3.x - P2.x)(1-t)t^2 + + // 3c1(P1.x - P0.x)(1-t)^2t + 6c1(P2.x - P1.x)(1-t)t^2 + 3c1(P3.x - P2.x)t^3 + // = (P0.y + 3 c0 (P1.x - P0.x))(1-t)^3 + + // 3(P1.y + 2 c0 (P2.x - P1.x) + c1 (P1.x - P0.x))(1-t)^2t + + // 3(P2.y + c0 (P3.x - P2.x) + 2 c1 (P2.x - P1.x))(1-t)t^2 + + // (P3.y + 3 c1 (P3.x - P2.x))t^3 + + // Therefore, the coefficients F0 to F3 of the focus curve are: + + // F0.x = (P0.x - 3 c0(P1.y - P0.y)) F0.y = (P0.y + 3 c0 (P1.x - P0.x)) + // F1.x = (P1.x - c1(P1.y - P0.y) - 2c0(P2.y - P1.y)) F1.y = (P1.y + 2 c0 (P2.x - P1.x) + c1 (P1.x - P0.x)) + // F2.x = (P2.x - 2 c1(P2.y - P1.y) - c0(P3.y - P2.y)) F2.y = (P2.y + c0 (P3.x - P2.x) + 2 c1 (P2.x - P1.x)) + // F3.x = (P3.x - 3 c1(P3.y - P2.y)) F3.y = (P3.y + 3 c1 (P3.x - P2.x)) + + res.p0.x = c.p0.x - 3*fRes[0]*(c.p1.y - c.p0.y); + res.p1.x = c.p1.x - fRes[1]*(c.p1.y - c.p0.y) - 2*fRes[0]*(c.p2.y - c.p1.y); + res.p2.x = c.p2.x - 2*fRes[1]*(c.p2.y - c.p1.y) - fRes[0]*(c.p3.y - c.p2.y); + res.p3.x = c.p3.x - 3*fRes[1]*(c.p3.y - c.p2.y); + + res.p0.y = c.p0.y + 3*fRes[0]*(c.p1.x - c.p0.x); + res.p1.y = c.p1.y + 2*fRes[0]*(c.p2.x - c.p1.x) + fRes[1]*(c.p1.x - c.p0.x); + res.p2.y = c.p2.y + fRes[0]*(c.p3.x - c.p2.x) + 2*fRes[1]*(c.p2.x - c.p1.x); + res.p3.y = c.p3.y + 3*fRes[1]*(c.p3.x - c.p2.x); +} + +bool Impl_calcSafeParams_focus( double& t1, + double& t2, + const Bezier& curve, + const Bezier& focus ) +{ + // now, we want to determine which normals of the original curve + // P(t) intersect with the focus curve F(t). The condition for + // this statement is P'(t)(P(t) - F) = 0, i.e. hodograph P'(t) and + // line through P(t) and F are perpendicular. + // If you expand this equation, you end up with something like + + // (\sum_{i=0}^n (P_i - F)B_i^n(t))^T (\sum_{j=0}^{n-1} n(P_{j+1} - P_j)B_j^{n-1}(t)) + + // Multiplying that out (as the scalar product is linear, we can + // extract some terms) yields: + + // (P_i - F)^T n(P_{j+1} - P_j) B_i^n(t)B_j^{n-1}(t) + ... + + // If we combine the B_i^n(t)B_j^{n-1}(t) product, we arrive at a + // Bernstein polynomial of degree 2n-1, as + + // \binom{n}{i}(1-t)^{n-i}t^i) \binom{n-1}{j}(1-t)^{n-1-j}t^j) = + // \binom{n}{i}\binom{n-1}{j}(1-t)^{2n-1-i-j}t^{i+j} + + // Thus, with the defining equation for a 2n-1 degree Bernstein + // polynomial + + // \sum_{i=0}^{2n-1} d_i B_i^{2n-1}(t) + + // the d_i are calculated as follows: + + // d_i = \sum_{j+k=i, j\in\{0,...,n\}, k\in\{0,...,n-1\}} \frac{\binom{n}{j}\binom{n-1}{k}}{\binom{2n-1}{i}} n (P_{k+1} - P_k)^T(P_j - F) + + // Okay, but F is now not a single point, but itself a curve + // F(u). Thus, for every value of u, we get a different 2n-1 + // bezier curve from the above equation. Therefore, we have a + // tensor product bezier patch, with the following defining + // equation: + + // d(t,u) = \sum_{i=0}^{2n-1} \sum_{j=0}^m B_i^{2n-1}(t) B_j^{m}(u) d_{ij}, where + // d_{ij} = \sum_{k+l=i, l\in\{0,...,n\}, k\in\{0,...,n-1\}} \frac{\binom{n}{l}\binom{n-1}{k}}{\binom{2n-1}{i}} n (P_{k+1} - P_k)^T(P_l - F_j) + + // as above, only that now F is one of the focus' control points. + + // Note the difference in the binomial coefficients to the + // reference paper, these formulas most probably contained a typo. + + // To determine, where D(t,u) is _not_ zero (these are the parts + // of the curve that don't share normals with the focus and can + // thus be safely clipped away), we project D(u,t) onto the + // (d(t,u), t) plane, determine the convex hull there and proceed + // as for the curve intersection part (projection is orthogonal to + // u axis, thus simply throw away u coordinate). + + // \fallfac are so-called falling factorials (see Concrete + // Mathematics, p. 47 for a definition). + + // now, for tensor product bezier curves, the convex hull property + // holds, too. Thus, we simply project the control points (t_{ij}, + // u_{ij}, d_{ij}) onto the (t,d) plane and calculate the + // intersections of the convex hull with the t axis, as for the + // bezier clipping case. + + // calc polygon of control points (t_{ij}, d_{ij}): + + const int n( 3 ); // cubic bezier curves, as a matter of fact + const int i_card( 2*n ); + const int j_card( n + 1 ); + const int k_max( n-1 ); + Polygon2D controlPolygon( i_card*j_card ); // vector of (t_{ij}, d_{ij}) in row-major order + + int i, j, k, l; // variable notation from formulas above and Sederberg article + Point2D::value_type d; + for( i=0; i=0 && k<=n-1); // k \in {0,...,n-1} + assert(l>=0 && l<=n); // l \in {0,...,n} + + // TODO: find, document and assert proper limits for n and int's max_val. + // This becomes important should anybody wants to use + // this code for higher-than-cubic beziers + d += static_cast(fallFac(n,l)*fallFac(n-1,k)*fac(i)) / + static_cast(fac(l)*fac(k) * fallFac(2*n-1,i)) * n * + ( (curve[k+1].x - curve[k].x)*(curve[l].x - focus[j].x) + // dot product here + (curve[k+1].y - curve[k].y)*(curve[l].y - focus[j].y) ); + } + + // Note that the t_{ij} values are evenly spaced on the + // [0,1] interval, thus t_{ij}=i/(2n-1) + controlPolygon[ i*j_card + j ] = Point2D( i/(2.0*n-1.0), d ); + } + } + +#ifndef WITH_SAFEFOCUSPARAM_DETAILED_TEST + + // calc safe parameter range, to determine [0,t1] and [t2,1] where + // no zero crossing is guaranteed. + return Impl_calcSafeParams( t1, t2, controlPolygon, 0.0, 0.0 ); + +#else + bool bRet( Impl_calcSafeParams( t1, t2, controlPolygon, 0.0, 0.0 ) ); + + Polygon2D convHull( convexHull( controlPolygon ) ); + + cout << "# convex hull testing (focus)" << endl + << "plot [t=0:1] "; + cout << "'-' using ($1):($2) title \"control polygon\" with lp, " + << "'-' using ($1):($2) title \"convex hull\" with lp" << endl; + + unsigned int count; + for( count=0; count void Impl_applySafeRanges_rec( std::back_insert_iterator< std::vector< std::pair > >& result, + double delta, + const Functor& safeRangeFunctor, + int recursionLevel, + const Bezier& c1_orig, + const Bezier& c1_part, + double last_t1_c1, + double last_t2_c1, + const Bezier& c2_orig, + const Bezier& c2_part, + double last_t1_c2, + double last_t2_c2 ) +{ + // check end condition + // =================== + + // TODO: tidy up recursion handling. maybe put everything in a + // struct and swap that here at method entry + + // TODO: Implement limit on recursion depth. Should that limit be + // reached, chances are that we're on a higher-order tangency. For + // this case, AW proposed to take the middle of the current + // interval, and to correct both curve's tangents at that new + // endpoint to be equal. That virtually generates a first-order + // tangency, and justifies to return a single intersection + // point. Otherwise, inside/outside test might fail here. + + for( int i=0; i last_t2_c2 - last_t1_c2 ) + { + // subdivide c1 + // ============ + + intervalMiddle = last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0; + + // subdivide at the middle of the interval (as + // we're not subdividing on the original + // curve, this simply amounts to subdivision + // at 0.5) + Impl_deCasteljauAt( part1, part2, c1_part, 0.5 ); + + // and descend recursively with swapped curves + Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1, + c2_orig, c2_part, last_t1_c2, last_t2_c2, + c1_orig, part1, last_t1_c1, intervalMiddle ); + + Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1, + c2_orig, c2_part, last_t1_c2, last_t2_c2, + c1_orig, part2, intervalMiddle, last_t2_c1 ); + } + else + { + // subdivide c2 + // ============ + + intervalMiddle = last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0; + + // subdivide at the middle of the interval (as + // we're not subdividing on the original + // curve, this simply amounts to subdivision + // at 0.5) + Impl_deCasteljauAt( part1, part2, c2_part, 0.5 ); + + // and descend recursively with swapped curves + Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1, + c2_orig, part1, last_t1_c2, intervalMiddle, + c1_orig, c1_part, last_t1_c1, last_t2_c1 ); + + Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1, + c2_orig, part2, intervalMiddle, last_t2_c2, + c1_orig, c1_part, last_t1_c1, last_t2_c1 ); + } + } + else + { + // apply calculated clip + // ===================== + + // clip safe ranges off c1_orig + Bezier c1_part1; + Bezier c1_part2; + Bezier c1_part3; + + // subdivide at t1_c1 + Impl_deCasteljauAt( c1_part1, c1_part2, c1_orig, t1_c1 ); + + // subdivide at t2_c1. As we're working on + // c1_part2 now, we have to adapt t2_c1 since + // we're no longer in the original parameter + // interval. This is based on the following + // assumption: t2_new = (t2-t1)/(1-t1), which + // relates the t2 value into the new parameter + // range [0,1] of c1_part2. + Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2_c1-t1_c1)/(1.0-t1_c1) ); + + // descend with swapped curves and c1_part1 as the + // remaining (middle) segment + Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1, + c2_orig, c2_part, last_t1_c2, last_t2_c2, + c1_orig, c1_part1, t1_c1, t2_c1 ); + } + } + } + } +} + +struct ClipBezierFunctor +{ + bool operator()( double& t1_c1, + double& t2_c1, + const Bezier& c1_orig, + const Bezier& c1_part, + const Bezier& c2_orig, + const Bezier& c2_part ) const + { + return Impl_calcClipRange( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part ); + } +}; + +struct BezierTangencyFunctor +{ + bool operator()( double& t1_c1, + double& t2_c1, + const Bezier& c1_orig, + const Bezier& c1_part, + const Bezier& c2_orig, + const Bezier& c2_part ) const + { + // calc focus curve of c2 + Bezier focus; + Impl_calcFocus(focus, c2_part); // need to use subdivided c2 + // here, as the whole curve is + // used for focus calculation + + // determine safe range on c1_orig + bool bRet( Impl_calcSafeParams_focus( t1_c1, t2_c1, + c1_orig, // use orig curve here, need t's on original curve + focus ) ); + + cerr << "range: " << t2_c1 - t1_c1 << ", ret: " << bRet << endl; + + return bRet; + } +}; + +/** Perform a bezier clip (curve against curve) + + @param result + Output iterator where the final t values are added to. This + iterator will remain empty, if there are no intersections. + + @param delta + Maximal allowed distance to true intersection (measured in the + original curve's coordinate system) + */ +void clipBezier( std::back_insert_iterator< std::vector< std::pair > >& result, + double delta, + const Bezier& c1, + const Bezier& c2 ) +{ +#if 0 + // first of all, determine list of collinear normals. Collinear + // normals typically separate two intersections, thus, subdivide + // at all collinear normal's t values beforehand. This will cater + // for tangent intersections, where two or more intersections are + // infinitesimally close together. + + // TODO: evaluate effects of higher-than-second-order + // tangencies. Sederberg et al. state that collinear normal + // algorithm then degrades quickly. + + std::vector< std::pair > results; + std::back_insert_iterator< std::vector< std::pair > > ii(results); + + Impl_calcCollinearNormals( ii, delta, 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 ); + + // As Sederberg's collinear normal theorem is only sufficient, not + // necessary for two intersections left and right, we've to test + // all segments generated by the collinear normal algorithm + // against each other. In other words, if the two curves are both + // divided in a left and a right part, the collinear normal + // theorem does _not_ state that the left part of curve 1 does not + // e.g. intersect with the right part of curve 2. + + // divide c1 and c2 at collinear normal intersection points + std::vector< Bezier > c1_segments( results.size()+1 ); + std::vector< Bezier > c2_segments( results.size()+1 ); + Bezier c1_remainder( c1 ); + Bezier c2_remainder( c2 ); + unsigned int i; + for( i=0; i 1 ) + { + Bezier c1( someCurves[0] ); + Bezier c2( someCurves[1] ); + + c1.p0.x += safeParams2_xOffset; + c1.p1.x += safeParams2_xOffset; + c1.p2.x += safeParams2_xOffset; + c1.p3.x += safeParams2_xOffset; + c2.p0.x += safeParams2_xOffset; + c2.p1.x += safeParams2_xOffset; + c2.p2.x += safeParams2_xOffset; + c2.p3.x += safeParams2_xOffset; + + double t1, t2; + + // output happens here + Impl_calcClipRange(t1, t2, c1, c1, c2, c2); + } +#endif + +#ifdef WITH_SAFEFOCUSPARAM_TEST + // test safe parameter range from focus algorithm + const double safeParamsFocus_xOffset( curr_Offset ); + curr_Offset += 20; + cout << "# safe param range from focus testing" << endl + << "plot [t=0.0:1.0] "; + for( i=0; i 1 ) + { + Bezier c1( someCurves[0] ); + Bezier c2( someCurves[1] ); + + c1.p0.x += safeParams3_xOffset; + c1.p1.x += safeParams3_xOffset; + c1.p2.x += safeParams3_xOffset; + c1.p3.x += safeParams3_xOffset; + c2.p0.x += safeParams3_xOffset; + c2.p1.x += safeParams3_xOffset; + c2.p2.x += safeParams3_xOffset; + c2.p3.x += safeParams3_xOffset; + + double t1, t2; + + Bezier focus; +#ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS + Impl_calcFocus( focus, c2 ); +#else + focus = c2; +#endif + + // determine safe range on c1, output happens here + Impl_calcSafeParams_focus( t1, t2, + c1, focus ); + } +#endif + +#ifdef WITH_BEZIERCLIP_TEST + std::vector< std::pair > result; + std::back_insert_iterator< std::vector< std::pair > > ii(result); + + // test full bezier clipping + const double bezierClip_xOffset( curr_Offset ); + cout << endl << endl << "# bezier clip testing" << endl + << "plot [t=0:1] "; + for( i=0; i + +struct Point2D +{ + typedef double value_type; + Point2D( double _x, double _y ) : x(_x), y(_y) {} + Point2D() : x(), y() {} + double x; + double y; +}; + +struct Bezier +{ + Point2D p0; + Point2D p1; + Point2D p2; + Point2D p3; + + Point2D& operator[]( int i ) { return reinterpret_cast(this)[i]; } + const Point2D& operator[]( int i ) const { return reinterpret_cast(this)[i]; } +}; + +struct FatLine +{ + // line L through p1 and p4 in normalized implicit form + double a; + double b; + double c; + + // the upper and lower distance from this line + double dMin; + double dMax; +}; + +template DataType calcLineDistance( const DataType& a, + const DataType& b, + const DataType& c, + const DataType& x, + const DataType& y ) +{ + return a*x + b*y + c; +} + +typedef std::vector< Point2D > Polygon2D; + +/* little abs template */ +template NumType absval( NumType x ) +{ + return x<0 ? -x : x; +} + +Polygon2D convexHull( const Polygon2D& rPoly ); + +// TODO: find proper epsilon here (try std::numeric_limits::epsilon()?)! +#define DBL_EPSILON 1.0e-100 + +/* little approximate comparisons */ +template bool tolZero( NumType n ) { return fabs(n) < DBL_EPSILON; } +template bool tolEqual( NumType n1, NumType n2 ) { return tolZero(n1-n2); } +template bool tolLessEqual( NumType n1, NumType n2 ) { return tolEqual(n1,n2) || n1 bool tolGreaterEqual( NumType n1, NumType n2 ) { return tolEqual(n1,n2) || n1>n2; } + +#endif // INCLUDED_BASEGFX_SOURCE_WORKBENCH_BEZIERCLIP_HXX + +/* vim:set shiftwidth=4 softtabstop=4 expandtab: */ diff --git a/basegfx/source/workbench/convexhull.cxx b/basegfx/source/workbench/convexhull.cxx new file mode 100644 index 000000000..ddcb04046 --- /dev/null +++ b/basegfx/source/workbench/convexhull.cxx @@ -0,0 +1,200 @@ +/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ +/* + * This file is part of the LibreOffice project. + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * This file incorporates work covered by the following license notice: + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed + * with this work for additional information regarding copyright + * ownership. The ASF licenses this file to you under the Apache + * License, Version 2.0 (the "License"); you may not use this file + * except in compliance with the License. You may obtain a copy of + * the License at http://www.apache.org/licenses/LICENSE-2.0 . + */ + +#include +#include + +#include + +/* Implements the theta function from Sedgewick: Algorithms in XXX, chapter 24 */ +template double theta( const PointType& p1, const PointType& p2 ) +{ + typename PointType::value_type dx, dy, ax, ay; + double t; + + dx = p2.x - p1.x; ax = absval( dx ); + dy = p2.y - p1.y; ay = absval( dy ); + t = (ax+ay == 0) ? 0 : (double) dy/(ax+ay); + if( dx < 0 ) + t = 2-t; + else if( dy < 0 ) + t = 4+t; + + return t*90.0; +} + +/* Model of LessThanComparable for theta sort. + * Uses the theta function from Sedgewick: Algorithms in XXX, chapter 24 + */ +template class ThetaCompare +{ +public: + explicit ThetaCompare( const PointType& rRefPoint ) : maRefPoint( rRefPoint ) {} + + bool operator() ( const PointType& p1, const PointType& p2 ) + { + // return true, if p1 precedes p2 in the angle relative to maRefPoint + return theta(maRefPoint, p1) < theta(maRefPoint, p2); + } + + double operator() ( const PointType& p ) const + { + return theta(maRefPoint, p); + } + +private: + PointType maRefPoint; +}; + +/* Implementation of the predicate 'counter-clockwise' for three points, from Sedgewick: Algorithms in XXX, chapter 24 */ +template typename PointType::value_type ccw( const PointType& p0, const PointType& p1, const PointType& p2, CmpFunctor& thetaCmp ) +{ + typename PointType::value_type dx1, dx2, dy1, dy2; + typename PointType::value_type theta0( thetaCmp(p0) ); + typename PointType::value_type theta1( thetaCmp(p1) ); + typename PointType::value_type theta2( thetaCmp(p2) ); + +#if 0 + if( theta0 == theta1 || + theta0 == theta2 || + theta1 == theta2 ) + { + // cannot reliably compare, as at least two points are + // theta-equal. See bug description below + return 0; + } +#endif + + dx1 = p1.x - p0.x; dy1 = p1.y - p0.y; + dx2 = p2.x - p0.x; dy2 = p2.y - p0.y; + + if( dx1*dy2 > dy1*dx2 ) + return +1; + + if( dx1*dy2 < dy1*dx2 ) + return -1; + + if( (dx1*dx2 < 0) || (dy1*dy2 < 0) ) + return -1; + + if( (dx1*dx1 + dy1*dy1) < (dx2*dx2 + dy2*dy2) ) + return +1; + + return 0; +} + +/* + Bug + === + + Sometimes, the resulting polygon is not the convex hull (see below + for an edge configuration to reproduce that problem) + + Problem analysis: + ================= + + The root cause of this bug is the fact that the second part of + the algorithm (the 'wrapping' of the point set) relies on the + previous theta sorting. Namely, it is required that the + generated point ordering, when interpreted as a polygon, is not + self-intersecting. This property, although, cannot be + guaranteed due to limited floating point accuracy. For example, + for two points very close together, and at the same time very + far away from the theta reference point p1, can appear on the + same theta value (because floating point accuracy is limited), + although on different rays to p1 when inspected locally. + + Example: + + / + P3 / + |\ / + | / + |/ \ + P2 \ + \ + ...____________\ + P1 + + Here, P2 and P3 are theta-equal relative to P1, but the local + ccw measure always says 'left turn'. Thus, the convex hull is + wrong at this place. + + Solution: + ========= + + If two points are theta-equal and checked via ccw, ccw must + also classify them as 'equal'. Thus, the second stage of the + convex hull algorithm sorts the first one out, effectively + reducing a cluster of theta-equal points to only one. This + single point can then be treated correctly. +*/ + +/* Implementation of Graham's convex hull algorithm, see Sedgewick: Algorithms in XXX, chapter 25 */ +Polygon2D convexHull( const Polygon2D& rPoly ) +{ + const Polygon2D::size_type N( rPoly.size() ); + Polygon2D result( N + 1 ); + std::copy(rPoly.begin(), rPoly.end(), result.begin()+1 ); + Polygon2D::size_type min, i; + + // determine safe point on hull (smallest y value) + for( min=1, i=2; i<=N; ++i ) + { + if( result[i].y < result[min].y ) + min = i; + } + + // determine safe point on hull (largest x value) + for( i=1; i<=N; ++i ) + { + if( result[i].y == result[min].y && + result[i].x > result[min].x ) + { + min = i; + } + } + + // TODO: add inner elimination optimization from Sedgewick: Algorithms in XXX, chapter 25 + // TODO: use radix sort instead of std::sort(), calc theta only once and store + + // setup first point and sort + std::swap( result[1], result[min] ); + ThetaCompare cmpFunc(result[1]); + std::sort( result.begin()+1, result.end(), cmpFunc ); + + // setup sentinel + result[0] = result[N]; + + // generate convex hull + Polygon2D::size_type M; + for( M=3, i=4; i<=N; ++i ) + { + while( ccw(result[M], result[M-1], result[i], cmpFunc) >= 0 ) + --M; + + ++M; + std::swap( result[M], result[i] ); + } + + // copy range [1,M] to output + return Polygon2D( result.begin()+1, result.begin()+M+1 ); +} + +/* vim:set shiftwidth=4 softtabstop=4 expandtab: */ diff --git a/basegfx/source/workbench/gauss.hxx b/basegfx/source/workbench/gauss.hxx new file mode 100644 index 000000000..fc352fe7e --- /dev/null +++ b/basegfx/source/workbench/gauss.hxx @@ -0,0 +1,170 @@ +/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ +/* + * This file is part of the LibreOffice project. + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * This file incorporates work covered by the following license notice: + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed + * with this work for additional information regarding copyright + * ownership. The ASF licenses this file to you under the Apache + * License, Version 2.0 (the "License"); you may not use this file + * except in compliance with the License. You may obtain a copy of + * the License at http://www.apache.org/licenses/LICENSE-2.0 . + */ + +/** This method eliminates elements below main diagonal in the given + matrix by gaussian elimination. + + @param matrix + The matrix to operate on. Last column is the result vector (right + hand side of the linear equation). After successful termination, + the matrix is upper triangular. The matrix is expected to be in + row major order. + + @param rows + Number of rows in matrix + + @param cols + Number of columns in matrix + + @param minPivot + If the pivot element gets lesser than minPivot, this method fails, + otherwise, elimination succeeds and true is returned. + + @return true, if elimination succeeded. + */ + +#ifndef INCLUDED_BASEGFX_SOURCE_WORKBENCH_GAUSS_HXX +#define INCLUDED_BASEGFX_SOURCE_WORKBENCH_GAUSS_HXX + +template +bool eliminate( Matrix& matrix, + int rows, + int cols, + const BaseType& minPivot ) +{ + BaseType temp; + + /* i, j, k *must* be signed, when looping like: j>=0 ! */ + /* eliminate below main diagonal */ + for(int i=0; i fabs(matrix[ max*cols + i ]) ) + max = j; + + /* check pivot value */ + if( fabs(matrix[ max*cols + i ]) < minPivot ) + return false; /* pivot too small! */ + + /* interchange rows 'max' and 'i' */ + for(int k=0; k=i; --k) + matrix[ j*cols + k ] -= matrix[ i*cols + k ] * + matrix[ j*cols + i ] / matrix[ i*cols + i ]; + } + + /* everything went well */ + return true; +} + +/** Retrieve solution vector of linear system by substituting backwards. + + This operation _relies_ on the previous successful + application of eliminate()! + + @param matrix + Matrix in upper diagonal form, as e.g. generated by eliminate() + + @param rows + Number of rows in matrix + + @param cols + Number of columns in matrix + + @param result + Result vector. Given matrix must have space for one column (rows entries). + + @return true, if back substitution was possible (i.e. no division + by zero occurred). + */ +template +bool substitute( const Matrix& matrix, + int rows, + int cols, + Vector& result ) +{ + BaseType temp; + + /* j, k *must* be signed, when looping like: j>=0 ! */ + /* substitute backwards */ + for(int j=rows-1; j>=0; --j) + { + temp = 0.0; + for(int k=j+1; k +bool solve( Matrix& matrix, + int rows, + int cols, + Vector& result, + BaseType minPivot ) +{ + if( eliminate(matrix, rows, cols, minPivot) ) + return substitute(matrix, rows, cols, result); + + return false; +} + +#endif // INCLUDED_BASEGFX_SOURCE_WORKBENCH_GAUSS_HXX + +/* vim:set shiftwidth=4 softtabstop=4 expandtab: */ -- cgit v1.2.3