summaryrefslogtreecommitdiffstats
path: root/kig/misc/cubic-common.cc
diff options
context:
space:
mode:
authorMichele Calgaro <michele.calgaro@yahoo.it>2020-12-11 21:51:45 +0900
committerMichele Calgaro <michele.calgaro@yahoo.it>2020-12-11 21:51:45 +0900
commit65a9f54e1051ee8ab936975e78dcb7117b265ef5 (patch)
tree4eb0dc3c3c34f8e744e70a6478cce1b5d18d9eca /kig/misc/cubic-common.cc
parent2f2f921072921d62baf5cffc44fb53ce516afebe (diff)
downloadtdeedu-65a9f54e1051ee8ab936975e78dcb7117b265ef5.tar.gz
tdeedu-65a9f54e1051ee8ab936975e78dcb7117b265ef5.zip
Renaming of files in preparation for code style tools.
Signed-off-by: Michele Calgaro <michele.calgaro@yahoo.it>
Diffstat (limited to 'kig/misc/cubic-common.cc')
-rw-r--r--kig/misc/cubic-common.cc527
1 files changed, 0 insertions, 527 deletions
diff --git a/kig/misc/cubic-common.cc b/kig/misc/cubic-common.cc
deleted file mode 100644
index 029f1194..00000000
--- a/kig/misc/cubic-common.cc
+++ /dev/null
@@ -1,527 +0,0 @@
-// Copyright (C) 2003 Dominique Devriese <devriese@kde.org>
-
-// This program is free software; you can redistribute it and/or
-// modify it under the terms of the GNU General Public License
-// as published by the Free Software Foundation; either version 2
-// of the License, or (at your option) any later version.
-
-// This program 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 General Public License for more details.
-
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
-// 02110-1301, USA.
-
-#include <config.h>
-
-#include "cubic-common.h"
-#include "kignumerics.h"
-#include "kigtransform.h"
-
-#ifdef HAVE_IEEEFP_H
-#include <ieeefp.h>
-#endif
-
-/*
- * coefficients of the cartesian equation for cubics
- */
-
-CubicCartesianData::CubicCartesianData()
-{
- std::fill( coeffs, coeffs + 10, 0 );
-}
-
-CubicCartesianData::CubicCartesianData(
- const double incoeffs[10] )
-{
- std::copy( incoeffs, incoeffs + 10, coeffs );
-}
-
-const CubicCartesianData calcCubicThroughPoints (
- const std::vector<Coordinate>& points )
-{
- // points is a vector of at most 9 points through which the cubic is
- // constrained.
- // this routine should compute the coefficients in the cartesian equation
- // they are defined up to a multiplicative factor.
- // since we don't know (in advance) which one of them is nonzero, we
- // simply keep all 10 parameters, obtaining a 9x10 linear system which
- // we solve using gaussian elimination with complete pivoting
- // If there are too few, then we choose some cool way to fill in the
- // empty parts in the matrix according to the LinearConstraints
- // given..
-
- // 9 rows, 10 columns..
- double row0[10];
- double row1[10];
- double row2[10];
- double row3[10];
- double row4[10];
- double row5[10];
- double row6[10];
- double row7[10];
- double row8[10];
- double *matrix[9] = {row0, row1, row2, row3, row4, row5, row6, row7, row8};
- double solution[10];
- int scambio[10];
-
- int numpoints = points.size();
- int numconstraints = 9;
-
- // fill in the matrix elements
- for ( int i = 0; i < numpoints; ++i )
- {
- double xi = points[i].x;
- double yi = points[i].y;
- matrix[i][0] = 1.0;
- matrix[i][1] = xi;
- matrix[i][2] = yi;
- matrix[i][3] = xi*xi;
- matrix[i][4] = xi*yi;
- matrix[i][5] = yi*yi;
- matrix[i][6] = xi*xi*xi;
- matrix[i][7] = xi*xi*yi;
- matrix[i][8] = xi*yi*yi;
- matrix[i][9] = yi*yi*yi;
- }
-
- for ( int i = 0; i < numconstraints; i++ )
- {
- if (numpoints >= 9) break; // don't add constraints if we have enough
- for (int j = 0; j < 10; ++j) matrix[numpoints][j] = 0.0;
- bool addedconstraint = true;
- switch (i)
- {
- case 0:
- matrix[numpoints][7] = 1.0;
- matrix[numpoints][8] = -1.0;
- break;
- case 1:
- matrix[numpoints][7] = 1.0;
- break;
- case 2:
- matrix[numpoints][9] = 1.0;
- break;
- case 3:
- matrix[numpoints][4] = 1.0;
- break;
- case 4:
- matrix[numpoints][5] = 1.0;
- break;
- case 5:
- matrix[numpoints][3] = 1.0;
- break;
- case 6:
- matrix[numpoints][1] = 1.0;
- break;
-
- default:
- addedconstraint = false;
- break;
- }
-
- if (addedconstraint) ++numpoints;
- }
-
- if ( ! GaussianElimination( matrix, numpoints, 10, scambio ) )
- return CubicCartesianData::invalidData();
- // fine della fase di eliminazione
- BackwardSubstitution( matrix, numpoints, 10, scambio, solution );
-
- // now solution should contain the correct coefficients..
- return CubicCartesianData( solution );
-}
-
-const CubicCartesianData calcCubicCuspThroughPoints (
- const std::vector<Coordinate>& points )
-{
- // points is a vector of at most 4 points through which the cubic is
- // constrained. Moreover the cubic is required to have a cusp at the
- // origin.
-
- // 9 rows, 10 columns..
- double row0[10];
- double row1[10];
- double row2[10];
- double row3[10];
- double row4[10];
- double row5[10];
- double row6[10];
- double row7[10];
- double row8[10];
- double *matrix[9] = {row0, row1, row2, row3, row4, row5, row6, row7, row8};
- double solution[10];
- int scambio[10];
-
- int numpoints = points.size();
- int numconstraints = 9;
-
- // fill in the matrix elements
- for ( int i = 0; i < numpoints; ++i )
- {
- double xi = points[i].x;
- double yi = points[i].y;
- matrix[i][0] = 1.0;
- matrix[i][1] = xi;
- matrix[i][2] = yi;
- matrix[i][3] = xi*xi;
- matrix[i][4] = xi*yi;
- matrix[i][5] = yi*yi;
- matrix[i][6] = xi*xi*xi;
- matrix[i][7] = xi*xi*yi;
- matrix[i][8] = xi*yi*yi;
- matrix[i][9] = yi*yi*yi;
- }
-
- for ( int i = 0; i < numconstraints; i++ )
- {
- if (numpoints >= 9) break; // don't add constraints if we have enough
- for (int j = 0; j < 10; ++j) matrix[numpoints][j] = 0.0;
- bool addedconstraint = true;
- switch (i)
- {
- case 0:
- matrix[numpoints][0] = 1.0; // through the origin
- break;
- case 1:
- matrix[numpoints][1] = 1.0;
- break;
- case 2:
- matrix[numpoints][2] = 1.0; // no first degree term
- break;
- case 3:
- matrix[numpoints][3] = 1.0; // a011 (x^2 coeff) = 0
- break;
- case 4:
- matrix[numpoints][4] = 1.0; // a012 (xy coeff) = 0
- break;
- case 5:
- matrix[numpoints][7] = 1.0;
- matrix[numpoints][8] = -1.0;
- break;
- case 6:
- matrix[numpoints][7] = 1.0;
- break;
- case 7:
- matrix[numpoints][9] = 1.0;
- break;
- case 8:
- matrix[numpoints][6] = 1.0;
- break;
-
- default:
- addedconstraint = false;
- break;
- }
-
- if (addedconstraint) ++numpoints;
- }
-
- if ( ! GaussianElimination( matrix, numpoints, 10, scambio ) )
- return CubicCartesianData::invalidData();
- // fine della fase di eliminazione
- BackwardSubstitution( matrix, numpoints, 10, scambio, solution );
-
- // now solution should contain the correct coefficients..
- return CubicCartesianData( solution );
-}
-
-const CubicCartesianData calcCubicNodeThroughPoints (
- const std::vector<Coordinate>& points )
-{
- // points is a vector of at most 6 points through which the cubic is
- // constrained. Moreover the cubic is required to have a node at the
- // origin.
-
- // 9 rows, 10 columns..
- double row0[10];
- double row1[10];
- double row2[10];
- double row3[10];
- double row4[10];
- double row5[10];
- double row6[10];
- double row7[10];
- double row8[10];
- double *matrix[9] = {row0, row1, row2, row3, row4, row5, row6, row7, row8};
- double solution[10];
- int scambio[10];
-
- int numpoints = points.size();
- int numconstraints = 9;
-
- // fill in the matrix elements
- for ( int i = 0; i < numpoints; ++i )
- {
- double xi = points[i].x;
- double yi = points[i].y;
- matrix[i][0] = 1.0;
- matrix[i][1] = xi;
- matrix[i][2] = yi;
- matrix[i][3] = xi*xi;
- matrix[i][4] = xi*yi;
- matrix[i][5] = yi*yi;
- matrix[i][6] = xi*xi*xi;
- matrix[i][7] = xi*xi*yi;
- matrix[i][8] = xi*yi*yi;
- matrix[i][9] = yi*yi*yi;
- }
-
- for ( int i = 0; i < numconstraints; i++ )
- {
- if (numpoints >= 9) break; // don't add constraints if we have enough
- for (int j = 0; j < 10; ++j) matrix[numpoints][j] = 0.0;
- bool addedconstraint = true;
- switch (i)
- {
- case 0:
- matrix[numpoints][0] = 1.0;
- break;
- case 1:
- matrix[numpoints][1] = 1.0;
- break;
- case 2:
- matrix[numpoints][2] = 1.0;
- break;
- case 3:
- matrix[numpoints][7] = 1.0;
- matrix[numpoints][8] = -1.0;
- break;
- case 4:
- matrix[numpoints][7] = 1.0;
- break;
- case 5:
- matrix[numpoints][9] = 1.0;
- break;
- case 6:
- matrix[numpoints][4] = 1.0;
- break;
- case 7:
- matrix[numpoints][5] = 1.0;
- break;
- case 8:
- matrix[numpoints][3] = 1.0;
- break;
-
- default:
- addedconstraint = false;
- break;
- }
-
- if (addedconstraint) ++numpoints;
- }
-
- if ( ! GaussianElimination( matrix, numpoints, 10, scambio ) )
- return CubicCartesianData::invalidData();
- // fine della fase di eliminazione
- BackwardSubstitution( matrix, numpoints, 10, scambio, solution );
-
- // now solution should contain the correct coefficients..
- return CubicCartesianData( solution );
-}
-
-/*
- * computation of the y value corresponding to some x value
- */
-
-double calcCubicYvalue ( double x, double ymin, double ymax, int root,
- CubicCartesianData data, bool& valid,
- int &numroots )
-{
- valid = true;
-
- // compute the third degree polinomial:
- double a000 = data.coeffs[0];
- double a001 = data.coeffs[1];
- double a002 = data.coeffs[2];
- double a011 = data.coeffs[3];
- double a012 = data.coeffs[4];
- double a022 = data.coeffs[5];
- double a111 = data.coeffs[6];
- double a112 = data.coeffs[7];
- double a122 = data.coeffs[8];
- double a222 = data.coeffs[9];
-
- // first the y^3 coefficient, it coming only from a222:
- double a = a222;
- // next the y^2 coefficient (from a122 and a022):
- double b = a122*x + a022;
- // next the y coefficient (from a112, a012 and a002):
- double c = a112*x*x + a012*x + a002;
- // finally the constant coefficient (from a111, a011, a001 and a000):
- double d = a111*x*x*x + a011*x*x + a001*x + a000;
-
- return calcCubicRoot ( ymin, ymax, a, b, c, d, root, valid, numroots );
-}
-
-const Coordinate calcCubicLineIntersect( const CubicCartesianData& cu,
- const LineData& l,
- int root, bool& valid )
-{
- assert( root == 1 || root == 2 || root == 3 );
-
- double a, b, c, d;
- calcCubicLineRestriction ( cu, l.a, l.b-l.a, a, b, c, d );
- int numroots;
- double param =
- calcCubicRoot ( -1e10, 1e10, a, b, c, d, root, valid, numroots );
- return l.a + param*(l.b - l.a);
-}
-
-/*
- * calculate the cubic polynomial resulting from the restriction
- * of a cubic to a line (defined by two "Coordinates": a point and a
- * direction)
- */
-
-void calcCubicLineRestriction ( CubicCartesianData data,
- Coordinate p, Coordinate v,
- double& a, double& b, double& c, double& d )
-{
- a = b = c = d = 0;
-
- double a000 = data.coeffs[0];
- double a001 = data.coeffs[1];
- double a002 = data.coeffs[2];
- double a011 = data.coeffs[3];
- double a012 = data.coeffs[4];
- double a022 = data.coeffs[5];
- double a111 = data.coeffs[6];
- double a112 = data.coeffs[7];
- double a122 = data.coeffs[8];
- double a222 = data.coeffs[9];
-
- // zero degree term
- d += a000;
-
- // first degree terms
- d += a001*p.x + a002*p.y;
- c += a001*v.x + a002*v.y;
-
- // second degree terms
- d += a011*p.x*p.x + a012*p.x*p.y + a022*p.y*p.y;
- c += 2*a011*p.x*v.x + a012*(p.x*v.y + v.x*p.y) + 2*a022*p.y*v.y;
- b += a011*v.x*v.x + a012*v.x*v.y + a022*v.y*v.y;
-
- // third degree terms: a111 x^3 + a222 y^3
- d += a111*p.x*p.x*p.x + a222*p.y*p.y*p.y;
- c += 3*(a111*p.x*p.x*v.x + a222*p.y*p.y*v.y);
- b += 3*(a111*p.x*v.x*v.x + a222*p.y*v.y*v.y);
- a += a111*v.x*v.x*v.x + a222*v.y*v.y*v.y;
-
- // third degree terms: a112 x^2 y + a122 x y^2
- d += a112*p.x*p.x*p.y + a122*p.x*p.y*p.y;
- c += a112*(p.x*p.x*v.y + 2*p.x*v.x*p.y) + a122*(v.x*p.y*p.y + 2*p.x*p.y*v.y);
- b += a112*(v.x*v.x*p.y + 2*v.x*p.x*v.y) + a122*(p.x*v.y*v.y + 2*v.x*v.y*p.y);
- a += a112*v.x*v.x*v.y + a122*v.x*v.y*v.y;
-}
-
-
-const CubicCartesianData calcCubicTransformation (
- const CubicCartesianData& data,
- const Transformation& t, bool& valid )
-{
- double a[3][3][3];
- double b[3][3][3];
- CubicCartesianData dataout;
-
- int icount = 0;
- for (int i=0; i < 3; i++)
- {
- for (int j=i; j < 3; j++)
- {
- for (int k=j; k < 3; k++)
- {
- a[i][j][k] = data.coeffs[icount++];
- if ( i < k )
- {
- if ( i == j ) // case aiik
- {
- a[i][i][k] /= 3.;
- a[i][k][i] = a[k][i][i] = a[i][i][k];
- }
- else if ( j == k ) // case aijj
- {
- a[i][j][j] /= 3.;
- a[j][i][j] = a[j][j][i] = a[i][j][j];
- }
- else // case aijk (i<j<k)
- {
- a[i][j][k] /= 6.;
- a[i][k][j] = a[j][i][k] = a[j][k][i] =
- a[k][i][j] = a[k][j][i] = a[i][j][k];
- }
- }
- }
- }
- }
-
- Transformation ti = t.inverse( valid );
- if ( ! valid ) return dataout;
-
- for (int i = 0; i < 3; i++)
- {
- for (int j = 0; j < 3; j++)
- {
- for (int k = 0; k < 3; k++)
- {
- b[i][j][k] = 0.;
- for (int ii = 0; ii < 3; ii++)
- {
- for (int jj = 0; jj < 3; jj++)
- {
- for (int kk = 0; kk < 3; kk++)
- {
- b[i][j][k] += a[ii][jj][kk]*ti.data( ii, i )*ti.data( jj, j )*ti.data( kk, k );
- }
- }
- }
- }
- }
- }
-
-// assert (fabs(b[0][1][2] - b[1][2][0]) < 1e-8); // test a couple of cases
-// assert (fabs(b[0][1][1] - b[1][1][0]) < 1e-8);
-
- // apparently, the above assertions are wrong ( due to rounding
- // errors, Maurizio and I hope :) ), so since the symmetry is not
- // present, we just take the sum of the parts of the matrix elements
- // that should be symmetric, instead of taking one of them, and
- // multiplying it..
-
- dataout.coeffs[0] = b[0][0][0];
- dataout.coeffs[1] = b[0][0][1] + b[0][1][0] + b[1][0][0];
- dataout.coeffs[2] = b[0][0][2] + b[0][2][0] + b[2][0][0];
- dataout.coeffs[3] = b[0][1][1] + b[1][0][1] + b[1][1][0];
- dataout.coeffs[4] = b[0][1][2] + b[0][2][1] + b[1][2][0] + b[1][0][2] + b[2][1][0] + b[2][0][1];
- dataout.coeffs[5] = b[0][2][2] + b[2][0][2] + b[2][2][0];
- dataout.coeffs[6] = b[1][1][1];
- dataout.coeffs[7] = b[1][1][2] + b[1][2][1] + b[2][1][1];
- dataout.coeffs[8] = b[1][2][2] + b[2][1][2] + b[2][2][1];
- dataout.coeffs[9] = b[2][2][2];
-
- return dataout;
-}
-
-bool operator==( const CubicCartesianData& lhs, const CubicCartesianData& rhs )
-{
- for ( int i = 0; i < 10; ++i )
- if ( lhs.coeffs[i] != rhs.coeffs[i] )
- return false;
- return true;
-}
-
-CubicCartesianData CubicCartesianData::invalidData()
-{
- CubicCartesianData ret;
- ret.coeffs[0] = double_inf;
- return ret;
-}
-
-bool CubicCartesianData::valid() const
-{
- return finite( coeffs[0] );
-}