valentina/src/libs/vgeometry/vspline.cpp

877 lines
28 KiB
C++
Raw Normal View History

/************************************************************************
**
** @file vspline.cpp
** @author Roman Telezhynskyi <dismine(at)gmail.com>
** @date November 15, 2013
**
** @brief
** @copyright
** This source code is part of the Valentine project, a pattern making
** program, whose allow create and modeling patterns of clothing.
** Copyright (C) 2013-2015 Valentina project
** <https://bitbucket.org/dismine/valentina> All Rights Reserved.
**
** Valentina 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 3 of the License, or
** (at your option) any later version.
**
** Valentina 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 Valentina. If not, see <http://www.gnu.org/licenses/>.
**
*************************************************************************/
2013-08-05 10:37:56 +02:00
#include "vspline.h"
#include "vspline_p.h"
#include <QDebug>
#include <QPainterPath>
#include <QtCore/qmath.h>
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief VSpline default constructor
*/
VSpline::VSpline()
:VAbstractCurve(GOType::Spline), d(new VSplineData)
{}
2013-08-05 10:37:56 +02:00
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief VSpline constructor.
* @param spline spline from which the copy.
*/
VSpline::VSpline ( const VSpline & spline )
:VAbstractCurve(spline), d(spline.d)
{}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief VSpline constructor.
* @param p1 first point spline.
* @param p4 last point spline.
* @param angle1 angle from first point to first control point.
* @param angle2 angle from second point to second control point.
* @param kCurve coefficient of curvature spline.
* @param kAsm1 coefficient of length first control line.
* @param kAsm2 coefficient of length second control line.
*/
VSpline::VSpline (VPointF p1, VPointF p4, qreal angle1, qreal angle2, qreal kAsm1, qreal kAsm2, qreal kCurve,
quint32 idObject, Draw mode)
:VAbstractCurve(GOType::Spline, idObject, mode), d(new VSplineData(p1, p4, angle1, angle2, kAsm1, kAsm2, kCurve))
{
CreateName();
2013-08-05 10:37:56 +02:00
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief VSpline constructor.
* @param p1 first point spline.
* @param p2 first control point.
* @param p3 second control point.
* @param p4 second point spline.
*/
VSpline::VSpline (VPointF p1, QPointF p2, QPointF p3, VPointF p4, qreal kCurve, quint32 idObject, Draw mode)
:VAbstractCurve(GOType::Spline, idObject, mode), d(new VSplineData(p1, p2, p3, p4, kCurve))
{
CreateName();
}
2013-08-05 10:37:56 +02:00
//---------------------------------------------------------------------------------------------------------------------
VSpline::~VSpline()
{}
2013-08-05 10:37:56 +02:00
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief GetLength return length of spline.
* @return length.
*/
qreal VSpline::GetLength () const
{
return LengthBezier ( GetP1().toQPointF(), d->p2, d->p3, GetP4().toQPointF());
2013-08-05 10:37:56 +02:00
}
//---------------------------------------------------------------------------------------------------------------------
qreal VSpline::LengthT(qreal t) const
{
if (t < 0 || t > 1)
{
qDebug()<<"Wrong value t.";
return 0;
}
QLineF seg1_2 ( GetP1 ().toQPointF(), GetP2 () );
seg1_2.setLength(seg1_2.length () * t);
QPointF p12 = seg1_2.p2();
QLineF seg2_3 ( GetP2 (), GetP3 () );
seg2_3.setLength(seg2_3.length () * t);
QPointF p23 = seg2_3.p2();
QLineF seg12_23 ( p12, p23 );
seg12_23.setLength(seg12_23.length () * t);
QPointF p123 = seg12_23.p2();
QLineF seg3_4 ( GetP3 (), GetP4 ().toQPointF() );
seg3_4.setLength(seg3_4.length () * t);
QPointF p34 = seg3_4.p2();
QLineF seg23_34 ( p23, p34 );
seg23_34.setLength(seg23_34.length () * t);
QPointF p234 = seg23_34.p2();
QLineF seg123_234 ( p123, p234 );
seg123_234.setLength(seg123_234.length () * t);
QPointF p1234 = seg123_234.p2();
return LengthBezier ( GetP1().toQPointF(), p12, p123, p1234);
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief CutSpline cut spline. GetPointP1() of base spline will return first point for first spline, GetPointP4()
* of base spline will return forth point of second spline.
* @param length length first spline
* @param spl1p2 second point of first spline
* @param spl1p3 third point of first spline
* @param spl2p2 second point of second spline
* @param spl2p3 third point of second spline
* @return point of cutting. This point is forth point of first spline and first point of second spline.
*/
QPointF VSpline::CutSpline ( qreal length, QPointF &spl1p2, QPointF &spl1p3, QPointF &spl2p2, QPointF &spl2p3 ) const
{
//Always need return two splines, so we must correct wrong length.
if (length < GetLength()*0.02)
{
length = GetLength()*0.02;
}
else if ( length > GetLength()*0.98)
{
length = GetLength()*0.98;
}
// Very stupid way find correct value of t.
// Better first compare with t = 0.5. Find length of spline.
// If length larger, take t = 0.75 and so on.
// If length less, take t = 0.25 and so on.
qreal parT = 0;
qreal step = 0.001;
while (1)
{
parT = parT + step;
qreal splLength = LengthT(parT);
if (splLength >= length || parT > 1)
{
break;
}
}
QLineF seg1_2 ( GetP1 ().toQPointF(), GetP2 () );
seg1_2.setLength(seg1_2.length () * parT);
QPointF p12 = seg1_2.p2();
QLineF seg2_3 ( GetP2 (), GetP3 () );
seg2_3.setLength(seg2_3.length () * parT);
QPointF p23 = seg2_3.p2();
QLineF seg12_23 ( p12, p23 );
seg12_23.setLength(seg12_23.length () * parT);
QPointF p123 = seg12_23.p2();
QLineF seg3_4 ( GetP3 (), GetP4 ().toQPointF() );
seg3_4.setLength(seg3_4.length () * parT);
QPointF p34 = seg3_4.p2();
QLineF seg23_34 ( p23, p34 );
seg23_34.setLength(seg23_34.length () * parT);
QPointF p234 = seg23_34.p2();
QLineF seg123_234 ( p123, p234 );
seg123_234.setLength(seg123_234.length () * parT);
QPointF p1234 = seg123_234.p2();
spl1p2 = p12;
spl1p3 = p123;
spl2p2 = p234;
spl2p3 = p34;
return p1234;
}
2013-08-05 10:37:56 +02:00
//---------------------------------------------------------------------------------------------------------------------
QPointF VSpline::CutSpline(qreal length, VSpline &spl1, VSpline &spl2) const
{
QPointF spl1p2;
QPointF spl1p3;
QPointF spl2p2;
QPointF spl2p3;
QPointF cutPoint = CutSpline (length, spl1p2, spl1p3, spl2p2, spl2p3 );
spl1 = VSpline(GetP1(), spl1p2, spl1p3, cutPoint, GetKcurve());
spl2 = VSpline(cutPoint, spl2p2, spl2p3, GetP4(), GetKcurve());
return cutPoint;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief GetPoints return list with spline points.
* @return list of points.
*/
QVector<QPointF> VSpline::GetPoints () const
{
return GetPoints(GetP1().toQPointF(), d->p2, d->p3, GetP4().toQPointF());
2013-08-05 10:37:56 +02:00
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief GetPoints return list with spline points.
* @param p1 first spline point.
* @param p2 first control point.
* @param p3 second control point.
* @param p4 last spline point.
* @return list of points.
*/
QVector<QPointF> VSpline::GetPoints (const QPointF &p1, const QPointF &p2, const QPointF &p3, const QPointF &p4)
{
2013-08-05 10:37:56 +02:00
QVector<QPointF> pvector;
QVector<qreal> x;
QVector<qreal> y;
QVector<qreal>& wx = x;
QVector<qreal>& wy = y;
x.append ( p1.x () );
y.append ( p1.y () );
PointBezier_r ( p1.x (), p1.y (), p2.x (), p2.y (),
p3.x (), p3.y (), p4.x (), p4.y (), 0, wx, wy );
x.append ( p4.x () );
y.append ( p4.y () );
for ( qint32 i = 0; i < x.count(); ++i )
{
pvector.append( QPointF ( x.at(i), y.at(i)) );
2013-08-05 10:37:56 +02:00
}
return pvector;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief LengthBezier return spline length using 4 spline point.
* @param p1 first spline point
* @param p2 first control point.
* @param p3 second control point.
* @param p4 last spline point.
* @return length.
*/
qreal VSpline::LengthBezier ( const QPointF &p1, const QPointF &p2, const QPointF &p3, const QPointF &p4 )
{
2013-08-05 10:37:56 +02:00
QPainterPath splinePath;
QVector<QPointF> points = GetPoints (p1, p2, p3, p4);
splinePath.moveTo(points.at(0));
for (qint32 i = 1; i < points.count(); ++i)
{
splinePath.lineTo(points.at(i));
2013-08-05 10:37:56 +02:00
}
return splinePath.length();
2013-08-05 10:37:56 +02:00
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief PointBezier_r find spline point using four point of spline.
* @param x1 х coordinate first point.
* @param y1 у coordinate first point.
* @param x2 х coordinate first control point.
* @param y2 у coordinate first control point.
* @param x3 х coordinate second control point.
* @param y3 у coordinate second control point.
* @param x4 х coordinate last point.
* @param y4 у coordinate last point.
* @param level level of recursion. In the begin 0.
* @param px list х coordinat spline points.
* @param py list у coordinat spline points.
*/
2013-08-05 10:37:56 +02:00
void VSpline::PointBezier_r ( qreal x1, qreal y1, qreal x2, qreal y2,
qreal x3, qreal y3, qreal x4, qreal y4,
qint16 level, QVector<qreal> &px, QVector<qreal> &py)
{
if (px.size() >= 2)
{
for (int i=1; i < px.size(); ++i)
{
if (QPointF(px.at(i-1), py.at(i-1)) == QPointF(px.at(i), py.at(i)))
{
qDebug("All neighbors points in path must be unique.");
}
}
}
const double curve_collinearity_epsilon = 1e-30;
const double curve_angle_tolerance_epsilon = 0.01;
2013-08-05 10:37:56 +02:00
const double m_angle_tolerance = 0.0;
enum curve_recursion_limit_e { curve_recursion_limit = 32 };
const double m_cusp_limit = 0.0;
double m_approximation_scale = 1.0;
double m_distance_tolerance_square;
2013-08-05 10:37:56 +02:00
m_distance_tolerance_square = 0.5 / m_approximation_scale;
m_distance_tolerance_square *= m_distance_tolerance_square;
if (level > curve_recursion_limit)
2013-08-05 10:37:56 +02:00
{
return;
}
2013-08-05 10:37:56 +02:00
// Calculate all the mid-points of the line segments
//----------------------
2014-03-03 19:14:56 +01:00
const double x12 = (x1 + x2) / 2;
const double y12 = (y1 + y2) / 2;
const double x23 = (x2 + x3) / 2;
const double y23 = (y2 + y3) / 2;
const double x34 = (x3 + x4) / 2;
const double y34 = (y3 + y4) / 2;
const double x123 = (x12 + x23) / 2;
const double y123 = (y12 + y23) / 2;
const double x234 = (x23 + x34) / 2;
const double y234 = (y23 + y34) / 2;
const double x1234 = (x123 + x234) / 2;
const double y1234 = (y123 + y234) / 2;
2013-08-05 10:37:56 +02:00
// Try to approximate the full cubic curve by a single straight line
//------------------
2014-03-03 19:14:56 +01:00
const double dx = x4-x1;
const double dy = y4-y1;
double d2 = fabs((x2 - x4) * dy - (y2 - y4) * dx);
double d3 = fabs((x3 - x4) * dy - (y3 - y4) * dx);
switch ((static_cast<int>(d2 > curve_collinearity_epsilon) << 1) +
static_cast<int>(d3 > curve_collinearity_epsilon))
2013-08-05 10:37:56 +02:00
{
case 0:
2014-03-03 19:14:56 +01:00
{
// All collinear OR p1==p4
//----------------------
2014-03-03 19:14:56 +01:00
double k = dx*dx + dy*dy;
if (k < 0.000000001)
2013-08-05 10:37:56 +02:00
{
d2 = CalcSqDistance(x1, y1, x2, y2);
d3 = CalcSqDistance(x4, y4, x3, y3);
2013-08-05 10:37:56 +02:00
}
else
{
k = 1 / k;
2014-03-03 19:14:56 +01:00
{
const double da1 = x2 - x1;
const double da2 = y2 - y1;
d2 = k * (da1*dx + da2*dy);
}
{
const double da1 = x3 - x1;
const double da2 = y3 - y1;
d3 = k * (da1*dx + da2*dy);
}
if (d2 > 0 && d2 < 1 && d3 > 0 && d3 < 1)
{
// Simple collinear case, 1---2---3---4
// We can leave just two endpoints
return;
}
if (d2 <= 0)
{
d2 = CalcSqDistance(x2, y2, x1, y1);
}
else if (d2 >= 1)
{
d2 = CalcSqDistance(x2, y2, x4, y4);
}
else
{
d2 = CalcSqDistance(x2, y2, x1 + d2*dx, y1 + d2*dy);
}
if (d3 <= 0)
{
d3 = CalcSqDistance(x3, y3, x1, y1);
}
else if (d3 >= 1)
{
d3 = CalcSqDistance(x3, y3, x4, y4);
}
else
{
d3 = CalcSqDistance(x3, y3, x1 + d3*dx, y1 + d3*dy);
}
2013-08-05 10:37:56 +02:00
}
if (d2 > d3)
2013-08-05 10:37:56 +02:00
{
if (d2 < m_distance_tolerance_square)
{
px.append(x2);
py.append(y2);
return;
}
2013-08-05 10:37:56 +02:00
}
else
2013-08-05 10:37:56 +02:00
{
if (d3 < m_distance_tolerance_square)
2013-08-05 10:37:56 +02:00
{
px.append(x3);
py.append(y3);
return;
}
}
break;
2014-03-03 19:14:56 +01:00
}
case 1:
2014-03-03 19:14:56 +01:00
{
// p1,p2,p4 are collinear, p3 is significant
2013-08-05 10:37:56 +02:00
//----------------------
if (d3 * d3 <= m_distance_tolerance_square * (dx*dx + dy*dy))
2013-08-05 10:37:56 +02:00
{
if (m_angle_tolerance < curve_angle_tolerance_epsilon)
2013-08-05 10:37:56 +02:00
{
px.append(x23);
py.append(y23);
return;
}
// Angle Condition
//----------------------
2014-03-03 19:14:56 +01:00
double da1 = fabs(atan2(y4 - y3, x4 - x3) - atan2(y3 - y2, x3 - x2));
if (da1 >= M_PI)
{
da1 = 2*M_PI - da1;
}
if (da1 < m_angle_tolerance)
{
2013-08-05 10:37:56 +02:00
px.append(x2);
py.append(y2);
px.append(x3);
py.append(y3);
2013-08-05 10:37:56 +02:00
return;
}
if (m_cusp_limit > 0.0 || m_cusp_limit < 0.0)
{
if (da1 > m_cusp_limit)
{
px.append(x3);
py.append(y3);
return;
}
}
2013-08-05 10:37:56 +02:00
}
break;
2014-03-03 19:14:56 +01:00
}
case 2:
2014-03-03 19:14:56 +01:00
{
// p1,p3,p4 are collinear, p2 is significant
2013-08-05 10:37:56 +02:00
//----------------------
if (d2 * d2 <= m_distance_tolerance_square * (dx*dx + dy*dy))
2013-08-05 10:37:56 +02:00
{
if (m_angle_tolerance < curve_angle_tolerance_epsilon)
{
px.append(x23);
py.append(y23);
return;
}
// Angle Condition
2013-08-05 10:37:56 +02:00
//----------------------
2014-03-03 19:14:56 +01:00
double da1 = fabs(atan2(y3 - y2, x3 - x2) - atan2(y2 - y1, x2 - x1));
if (da1 >= M_PI)
2013-08-05 10:37:56 +02:00
{
da1 = 2*M_PI - da1;
}
if (da1 < m_angle_tolerance)
{
2013-08-05 10:37:56 +02:00
px.append(x2);
py.append(y2);
px.append(x3);
py.append(y3);
2013-08-05 10:37:56 +02:00
return;
}
if (m_cusp_limit > 0.0 || m_cusp_limit < 0.0)
2013-08-05 10:37:56 +02:00
{
if (da1 > m_cusp_limit)
{
px.append(x2);
py.append(y2);
return;
}
}
}
break;
2014-03-03 19:14:56 +01:00
}
case 3:
2014-03-03 19:14:56 +01:00
{
// Regular case
//-----------------
if ((d2 + d3)*(d2 + d3) <= m_distance_tolerance_square * (dx*dx + dy*dy))
{
// If the curvature doesn't exceed the distance_tolerance value
// we tend to finish subdivisions.
//----------------------
if (m_angle_tolerance < curve_angle_tolerance_epsilon)
{
px.append(x23);
py.append(y23);
return;
}
// Angle & Cusp Condition
//----------------------
2014-03-03 19:14:56 +01:00
const double k = atan2(y3 - y2, x3 - x2);
double da1 = fabs(k - atan2(y2 - y1, x2 - x1));
double da2 = fabs(atan2(y4 - y3, x4 - x3) - k);
if (da1 >= M_PI)
{
da1 = 2*M_PI - da1;
}
if (da2 >= M_PI)
{
da2 = 2*M_PI - da2;
}
if (da1 + da2 < m_angle_tolerance)
{
// Finally we can stop the recursion
//----------------------
px.append(x23);
py.append(y23);
2013-08-05 10:37:56 +02:00
return;
}
if (m_cusp_limit > 0.0 || m_cusp_limit < 0.0)
{
if (da1 > m_cusp_limit)
{
px.append(x2);
py.append(y2);
return;
}
if (da2 > m_cusp_limit)
{
px.append(x3);
py.append(y3);
return;
}
}
2013-08-05 10:37:56 +02:00
}
break;
2014-03-03 19:14:56 +01:00
}
default:
break;
2013-08-05 10:37:56 +02:00
}
2013-08-05 10:37:56 +02:00
// Continue subdivision
//----------------------
2013-08-20 12:26:02 +02:00
PointBezier_r(x1, y1, x12, y12, x123, y123, x1234, y1234, static_cast<qint16>(level + 1), px, py);
PointBezier_r(x1234, y1234, x234, y234, x34, y34, x4, y4, static_cast<qint16>(level + 1), px, py);
2013-08-05 10:37:56 +02:00
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief CalcSqDistance calculate squared distance.
* @param x1 х coordinate first point.
* @param y1 у coordinate first point.
* @param x2 х coordinate second point.
* @param y2 у coordinate second point.
* @return squared length.
*/
qreal VSpline::CalcSqDistance (qreal x1, qreal y1, qreal x2, qreal y2)
{
2013-08-05 10:37:56 +02:00
qreal dx = x2 - x1;
qreal dy = y2 - y1;
return dx * dx + dy * dy;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief CreateName create spline name.
*/
void VSpline::CreateName()
{
setName(SPL_ + QString("%1_%2").arg(this->GetP1().name(), this->GetP4().name()));
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief SplinePoints return list with spline points.
* @param p1 first spline point.
* @param p4 last spline point.
* @param angle1 angle from first point to first control point.
* @param angle2 angle from second point to second control point.
* @param kAsm1 coefficient of length first control line.
* @param kAsm2 coefficient of length second control line.
* @param kCurve coefficient of curvature spline.
* @return list with spline points.
*/
2014-05-11 20:15:32 +02:00
// cppcheck-suppress unusedFunction
QVector<QPointF> VSpline::SplinePoints(const QPointF &p1, const QPointF &p4, qreal angle1, qreal angle2, qreal kAsm1,
qreal kAsm2, qreal kCurve)
{
2013-08-28 10:55:11 +02:00
QLineF p1pX(p1.x(), p1.y(), p1.x() + 100, p1.y());
p1pX.setAngle( angle1 );
qreal L = 0, radius = 0, angle = 90;
radius = QLineF(QPointF(p1.x(), p4.y()), p4).length();
2013-08-28 10:55:11 +02:00
L = kCurve * radius * 4 / 3 * tan( angle * M_PI / 180.0 / 4 );
QLineF p1p2(p1.x(), p1.y(), p1.x() + L * kAsm1, p1.y());
p1p2.setAngle(angle1);
QLineF p4p3(p4.x(), p4.y(), p4.x() + L * kAsm2, p4.y());
p4p3.setAngle(angle2);
QPointF p2 = p1p2.p2();
QPointF p3 = p4p3.p2();
return GetPoints(p1, p2, p3, p4);
}
//---------------------------------------------------------------------------------------------------------------------
VSpline &VSpline::operator =(const VSpline &spline)
{
if ( &spline == this )
{
return *this;
}
2014-06-25 16:08:12 +02:00
VAbstractCurve::operator=(spline);
d = spline.d;
return *this;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief GetP1 return first spline point.
* @return first point.
*/
VPointF VSpline::GetP1() const
{
return d->p1;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief GetP2 return first control point.
* @return first control point.
*/
QPointF VSpline::GetP2() const
{
return d->p2;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief GetP3 return second control point.
* @return second control point.
*/
QPointF VSpline::GetP3() const
{
return d->p3;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief GetP4 return last spline point.
* @return остання точка сплайну.
*/
VPointF VSpline::GetP4() const
{
return d->p4;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief GetAngle1 return first angle control line.
* @return angle.
*/
qreal VSpline::GetStartAngle() const
{
return d->angle1;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief GetAngle2 return second angle control line.
* @return angle.
*/
qreal VSpline::GetEndAngle() const
{
return d->angle2;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief GetKasm1 return coefficient of length first control line.
* @return coefficient.
*/
qreal VSpline::GetKasm1() const
{
return d->kAsm1;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief GetKasm2 return coefficient of length second control line.
* @return coefficient.
*/
qreal VSpline::GetKasm2() const
{
return d->kAsm2;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief GetKcurve return coefficient of curvature spline.
* @return coefficient
*/
qreal VSpline::GetKcurve() const
{
return d->kCurve;
}
//---------------------------------------------------------------------------------------------------------------------
void VSpline::SetKcurve(qreal factor)
{
if (factor > 0)
{
d->kCurve = factor;
}
}
//---------------------------------------------------------------------------------------------------------------------
2015-03-02 18:54:33 +01:00
int VSpline::Sign(long double ld)
{
2015-03-02 18:11:43 +01:00
if (qAbs(ld)<0.00000000001)
{
return 0;
}
return (ld>0) ? 1 : -1;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief Cubic Cubic equation solution. Real coefficients case.
*
* This method use method Vieta-Cardano for eval cubic equations.
* Cubic equation write in form x3+a*x2+b*x+c=0.
*
* Output:
* 3 real roots -> then x is filled with them;
* 1 real + 2 complex -> x[0] is real, x[1] is real part of complex roots, x[2] - non-negative imaginary part.
*
* @param x solution array (size 3).
* @param a coefficient
* @param b coefficient
* @param c coefficient
* @return 3 - 3 real roots;
* 1 - 1 real root + 2 complex;
* 2 - 1 real root + complex roots imaginary part is zero (i.e. 2 real roots).
*/
2015-03-02 18:54:33 +01:00
qint32 VSpline::Cubic(QVector<qreal> &x, qreal a, qreal b, qreal c)
{
//To find cubic equation roots in the case of real coefficients, calculated at the beginning
const qreal q = (pow(a, 2) - 3*b)/9.;
const qreal r = (2*pow(a, 3) - 9*a*b + 27.*c)/54.;
if (pow(r, 2) < pow(q, 3))
{ // equation has three real roots, use formula Vieta
const qreal t = acos(r/sqrt(pow(q, 3)))/3.;
x.insert(0, -2.*sqrt(q)*cos(t)-a/3);
x.insert(1, -2.*sqrt(q)*cos(t + (2*M_2PI/3.)) - a/3.);
x.insert(2, -2.*sqrt(q)*cos(t - (2*M_2PI/3.)) - a/3.);
return(3);
}
else
{ // 1 real root + 2 complex
//Formula Cardano
const qreal aa = -Sign(r)*pow(fabs(r)+sqrt(pow(r, 2)-pow(q, 3)), 1./3.);
const qreal bb = Sign(aa) == 0 ? 0 : q/aa;
x.insert(0, aa+bb-a/3.); // Real root
x.insert(1, (-0.5)*(aa+bb)-a/3.); //Complex root
x.insert(2, (sqrt(3.)*0.5)*fabs(aa-bb)); // Complex root
if (qFuzzyCompare(x.at(2) + 1, 0. + 1))
{
return(2);
}
return(1);
}
}
//---------------------------------------------------------------------------------------------------------------------
QVector<qreal> VSpline::CalcT (qreal curveCoord1, qreal curveCoord2, qreal curveCoord3,
qreal curveCoord4, qreal pointCoord) const
{
const qreal a = -curveCoord1 + 3*curveCoord2 - 3*curveCoord3 + curveCoord4;
const qreal b = 3*curveCoord1 - 6*curveCoord2 + 3*curveCoord3;
const qreal c = -3*curveCoord1 + 3*curveCoord2;
const qreal d = -pointCoord + curveCoord1;
QVector<qreal> t = QVector<qreal>(3, -1);
Cubic(t, b/a, c/a, d/a);
QVector<qreal> retT;
for (int i=0; i < t.size(); ++i)
{
if ( t.at(i) >= 0 && t.at(i) <= 1 )
{
retT.append(t.at(i));
}
}
return retT;
}
//---------------------------------------------------------------------------------------------------------------------
/**
* @brief VSpline::ParamT calculate t coeffient that reprezent point on curve.
*
* Each point that belongs to Cubic Bézier curve can be shown by coefficient in interval [0; 1].
*
* @param pBt point on curve
* @return t coeffient that reprezent this point on curve. Return -1 if point doesn't belongs to curve.
*/
qreal VSpline::ParamT (const QPointF &pBt) const
{
QVector<qreal> ts;
// Calculate t coefficient for each axis
ts += CalcT (GetP1().toQPointF().x(), d->p2.x(), d->p3.x(), GetP4().toQPointF().x(), pBt.x());
ts += CalcT (GetP1().toQPointF().y(), d->p2.y(), d->p3.y(), GetP4().toQPointF().y(), pBt.y());
2015-03-02 18:11:43 +01:00
if (ts.isEmpty())
{
return -1; // We don't have candidates
}
qreal tx = -1;
qreal eps = 3; // Error calculation
// In morst case we will have 6 result in interval [0; 1].
// Here we try find closest to our point.
for (int i=0; i< ts.size(); ++i)
{
const qreal t = ts.at(i);
const QPointF p0 = GetP1().toQPointF();
const QPointF p1 = d->p2;
const QPointF p2 = d->p3;
const QPointF p3 = GetP4().toQPointF();
//The explicit form of the Cubic Bézier curve
const qreal pointX = pow(1-t, 3)*p0.x() + 3*pow(1-t, 2)*t*p1.x() + 3*(1-t)*pow(t, 2)*p2.x() + pow(t, 3)*p3.x();
const qreal pointY = pow(1-t, 3)*p0.y() + 3*pow(1-t, 2)*t*p1.y() + 3*(1-t)*pow(t, 2)*p2.y() + pow(t, 3)*p3.y();
const QLineF line(pBt, QPointF(pointX, pointY));
if (line.length() <= eps)
{
tx = t;
eps = line.length(); //Next point should be even closest
}
}
return tx;
}