//-----------------------------------------------------------------------------
// tstrokedeformations.cpp
//
// Revision History:
// 16/07/2002 m_lengthOfDeformation is almost TConsts::epsilon
//-----------------------------------------------------------------------------
#include "tstrokedeformations.h"
#include "tcurveutil.h"
#include "tmathutil.h"
#include "tstroke.h"
using namespace std;
//=============================================================================
namespace
{
const double nine_inv = 1.0 / 9.0;
/*!
_____| r
/ |\ inner
__/ \___|
| r
outer
{ 1 if r <= r
{ inner
{
bowl potential(r)= { cos( (r - r ) / (r - r ) * pi_2 ) if r < r <= r
{ i o i i o
{
{ 0 if r > r
o
*/
struct bowlPotential {
double m_radiusInner;
double m_radiusOuter;
bowlPotential(double radiusInner,
double radiusOuter)
: m_radiusInner(radiusInner), m_radiusOuter(radiusOuter)
{
assert(m_radiusInner > 0);
assert(m_radiusOuter >= m_radiusInner);
}
virtual double value(double radiusToTest)
{
assert(radiusToTest >= 0);
if (radiusToTest <= m_radiusInner)
return 1.0;
if (radiusToTest > m_radiusOuter)
return 0.0;
return 0.5 * (1.0 + cos((radiusToTest - m_radiusInner) / (m_radiusOuter - m_radiusInner) * TConsts::pi));
}
virtual double gradient(double radiusToTest)
{
assert(radiusToTest >= 0);
if (radiusToTest <= m_radiusInner || radiusToTest > m_radiusOuter)
return 0.0;
double den = TConsts::pi / (m_radiusOuter - m_radiusInner);
return -0.5 * den * sin(den * (radiusToTest - m_radiusInner));
}
virtual ~bowlPotential() {}
};
double wyvillPotential(double r, double R)
{
if (0.0 == R)
return 0.0;
if (0 > r || r >= R)
return 0.0;
double ratio2 = sq(r / R);
double ratio4 = sq(ratio2);
double ratio6 = ratio2 * ratio4;
return 1.0 + (17.0 * ratio4 - (4.0 * ratio6 + 22.0 * ratio2)) * nine_inv;
}
/*
*/
double derivateOfWyvillPotential(double r, double R)
{
if (0.0 == R)
return 0.0;
if (0 > r || r > R)
return 0.0;
double inv_of_R = 1.0 / R;
double ratio = r / R;
double ratio2 = sq(ratio);
double ratio3 = ratio * ratio2;
double ratio5 = ratio3 * ratio2;
return inv_of_R * nine_inv * (68.0 * ratio3 - (24.0 * ratio5 + 66.0 * ratio));
}
const double c_maxLenghtOfGaussian = 3.0;
/*
*/
double gaussianPotential(double x)
{
return exp(-sq(x));
}
/*
*/
double derivateOfGaussianPotential(double x)
{
return -2 * x * exp(-sq(x));
}
/*
Check if vector distance is in segment.
bool pointProjectionIsInSegment( const TPointD& p, const TSegment& seg )
{
TPointD
a ( p - seg.getP0()),
b ( seg.getP1() - seg.getP0());
double b2 = b*b;
if( ! isAlmostZero(b2) )
{
b2 = a * b / b2;
if ( 0 <= b2 && b2 <= 1.0 ) return true;
}
return false;
}
*/
double normOfGradientOfWyvillPotential(double r, double R)
{
TPointD grad;
grad.x = 1.0;
grad.y = derivateOfWyvillPotential(r, R);
return norm(grad);
}
} // end of unnamed namespace
//=============================================================================
struct TStrokePointDeformation::Imp {
TPointD m_circleCenter;
double m_circleRadius;
TPointD *m_vect;
bowlPotential *m_potential;
Imp(const TPointD ¢er,
double radius)
: m_circleCenter(center), m_circleRadius(radius), m_vect(0)
{
m_potential = new bowlPotential(0.3 * m_circleRadius, m_circleRadius);
}
Imp(const TPointD &vect,
const TPointD ¢er,
double radius)
: m_circleCenter(center), m_circleRadius(radius), m_vect(new TPointD(vect))
{
m_potential = new bowlPotential(0.3 * m_circleRadius, m_circleRadius);
}
~Imp()
{
delete m_vect;
delete m_potential;
}
};
TStrokePointDeformation::TStrokePointDeformation(const TPointD ¢er,
double radius)
: m_imp(new Imp(center, radius))
{
}
//-----------------------------------------------------------------------------
TStrokePointDeformation::TStrokePointDeformation(const TPointD &vect,
const TPointD ¢er,
double radius)
: m_imp(new Imp(vect, center, radius))
{
}
//-----------------------------------------------------------------------------
TStrokePointDeformation::~TStrokePointDeformation()
{
}
//-----------------------------------------------------------------------------
TThickPoint TStrokePointDeformation::getDisplacementForControlPoint(const TStroke &stroke, UINT n) const
{
// riferimento ad un punto ciccione della stroke
TPointD pntOfStroke(convert(stroke.getControlPoint(n)));
double d = tdistance(pntOfStroke, m_imp->m_circleCenter);
if (m_imp->m_vect)
return m_imp->m_potential->value(d) * TThickPoint(*m_imp->m_vect, 0);
else {
double outVal = m_imp->m_potential->value(d);
return TThickPoint(outVal, outVal, 0);
}
}
TThickPoint TStrokePointDeformation::getDisplacementForControlPointLen(const TStroke &stroke, double cpLen) const
{
assert(0);
return TThickPoint();
}
//-----------------------------------------------------------------------------
TThickPoint TStrokePointDeformation::getDisplacement(const TStroke &stroke, double w) const
{
// riferimento ad un punto ciccione della stroke
TThickPoint
thickPnt = m_imp->m_vect ? stroke.getControlPointAtParameter(w) : stroke.getThickPoint(w);
assert(thickPnt != TConsts::natp);
TPointD pntOfStroke(convert(thickPnt));
double d = tdistance(pntOfStroke, m_imp->m_circleCenter);
if (m_imp->m_vect)
return m_imp->m_potential->value(d) * TThickPoint(*m_imp->m_vect, 0);
else {
double outVal = m_imp->m_potential->value(d);
return TThickPoint(outVal, outVal, 0);
}
}
//-----------------------------------------------------------------------------
double TStrokePointDeformation::getDelta(const TStroke &stroke, double w) const
{
// reference to a thickpoint
TThickPoint
thickPnt = m_imp->m_vect ? stroke.getControlPointAtParameter(w) : stroke.getThickPoint(w);
assert(thickPnt != TConsts::natp);
TPointD
pntOfStroke = convert(thickPnt);
double d = tdistance(pntOfStroke, m_imp->m_circleCenter);
return m_imp->m_potential->gradient(d);
}
//-----------------------------------------------------------------------------
double TStrokePointDeformation::getMaxDiff() const
{
return 0.005;
}
//=============================================================================
TStrokeParamDeformation::TStrokeParamDeformation(const TStroke *ref,
const TPointD &vect,
double s,
double l)
: m_pRef(ref), m_startParameter(s), m_lengthOfDeformation(l), m_vect(new TPointD(vect))
{
assert(m_lengthOfDeformation >= 0);
if (isAlmostZero(m_lengthOfDeformation))
m_lengthOfDeformation = TConsts::epsilon;
}
//-----------------------------------------------------------------------------
TStrokeParamDeformation::TStrokeParamDeformation(const TStroke *ref,
double s,
double l)
: m_pRef(ref), m_startParameter(s), m_lengthOfDeformation(l), m_vect(0)
{
assert(m_lengthOfDeformation >= 0);
if (isAlmostZero(m_lengthOfDeformation))
m_lengthOfDeformation = TConsts::epsilon;
}
//-----------------------------------------------------------------------------
TThickPoint TStrokeParamDeformation::getDisplacementForControlPoint(const TStroke &stroke, UINT n) const
{
// potenziale exp^(-x^2) limitato tra [-c_maxLenghtOfGaussian,c_maxLenghtOfGaussian]
double diff = stroke.getLengthAtControlPoint(n);
diff = diff - m_startParameter;
if (fabs(diff) <= m_lengthOfDeformation) {
// modulo il vettore spostamento in funzione del potenziale
// il termine (1.0/m_lengthOfDeformation) * 3 scala
// il punto in diff su un sistema di coordinate
// normalizzato, associato con la curva exp^(-x^2)
diff *= (1.0 / m_lengthOfDeformation) * c_maxLenghtOfGaussian;
if (m_vect)
return gaussianPotential(diff) * TThickPoint(*m_vect, 0);
else {
double outVal = gaussianPotential(diff);
return TThickPoint(outVal, outVal, 0);
}
}
return TThickPoint();
}
TThickPoint TStrokeParamDeformation::getDisplacementForControlPointLen(const TStroke &stroke, double cpLenDiff) const
{
// potenziale exp^(-x^2) limitato tra [-c_maxLenghtOfGaussian,c_maxLenghtOfGaussian]
//double diff = stroke.getLengthAtControlPoint(n);
//double diff = cpLen;
double diff = cpLenDiff;
if (fabs(diff) <= m_lengthOfDeformation) {
// modulo il vettore spostamento in funzione del potenziale
// il termine (1.0/m_lengthOfDeformation) * 3 scala
// il punto in diff su un sistema di coordinate
// normalizzato, associato con la curva exp^(-x^2)
diff *= (1.0 / m_lengthOfDeformation) * c_maxLenghtOfGaussian;
if (m_vect)
return gaussianPotential(diff) * TThickPoint(*m_vect, 0);
else {
double outVal = gaussianPotential(diff);
return TThickPoint(outVal, outVal, 0);
}
}
return TThickPoint();
}
//-----------------------------------------------------------------------------
TThickPoint TStrokeParamDeformation::getDisplacement(const TStroke &stroke, double w) const
{
// potenziale exp^(-x^2) limitato tra [-c_maxLenghtOfGaussian,c_maxLenghtOfGaussian]
double diff = stroke.getLength(w);
diff = diff - m_startParameter;
if (fabs(diff) <= m_lengthOfDeformation) {
// modulo il vettore spostamento in funzione del potenziale
// il termine (1.0/m_lengthOfDeformation) * 3 scala
// il punto in diff su un sistema di coordinate
// normalizzato, associato con la curva exp^(-x^2)
diff *= (1.0 / m_lengthOfDeformation) * c_maxLenghtOfGaussian;
if (m_vect)
return gaussianPotential(diff) * TThickPoint(*m_vect, 0);
else {
double outVal = gaussianPotential(diff);
return TThickPoint(outVal, outVal, 0);
}
}
return TThickPoint();
}
//-----------------------------------------------------------------------------
double TStrokeParamDeformation::getDelta(const TStroke &stroke, double w) const
{
// potenziale exp^(-x^2) limitato tra [-c_maxLenghtOfGaussian,c_maxLenghtOfGaussian]
double diff = stroke.getLength(w);
diff = diff - m_startParameter;
if (fabs(diff) <= m_lengthOfDeformation) {
// modulo il vettore spostamento in funzione del potenziale
// il termine (1.0/m_lengthOfDeformation) * 3 scala
// il punto in diff su un sistema di coordinate
// normalizzato, associato con la curva exp^(-x^2)
diff *= (1.0 / m_lengthOfDeformation) * c_maxLenghtOfGaussian;
return derivateOfGaussianPotential(diff);
}
return 0;
}
//-----------------------------------------------------------------------------
double TStrokeParamDeformation::getMaxDiff() const
{
return 0.09;
}
//-----------------------------------------------------------------------------
TStrokeParamDeformation::~TStrokeParamDeformation()
{
delete m_vect;
}
//=============================================================================
//=============================================================================
/*
*/
TStrokeBenderDeformation::TStrokeBenderDeformation(const TStroke *ref,
double s,
double l)
: m_pRef(ref), m_startLength(s), m_lengthOfDeformation(l), m_vect(0), m_versus(INNER)
{
assert(m_lengthOfDeformation >= 0);
if (isAlmostZero(m_lengthOfDeformation))
m_lengthOfDeformation = TConsts::epsilon;
}
//-----------------------------------------------------------------------------
TStrokeBenderDeformation::TStrokeBenderDeformation(const TStroke *ref,
const TPointD &vect,
double angle,
double s,
int innerOrOuter,
double l)
: m_pRef(ref), m_startLength(s), m_lengthOfDeformation(l), m_vect(new TPointD(vect)), m_versus(innerOrOuter), m_angle(angle)
{
assert(m_lengthOfDeformation >= 0);
if (isAlmostZero(m_lengthOfDeformation))
m_lengthOfDeformation = TConsts::epsilon;
}
//-----------------------------------------------------------------------------
TStrokeBenderDeformation::~TStrokeBenderDeformation()
{
delete m_vect;
}
//-----------------------------------------------------------------------------
double TStrokeBenderDeformation::getMaxDiff() const
{
// return 0.09; OK per gaussiana
return 0.4;
}
//-----------------------------------------------------------------------------
TThickPoint TStrokeBenderDeformation::getDisplacementForControlPoint(const TStroke &s, UINT n) const
{
// potenziale exp^(-x^2) limitato tra [-c_maxLenghtOfGaussian,c_maxLenghtOfGaussian]
double strokeLengthAtParameter = s.getLengthAtControlPoint(n);
double diff = strokeLengthAtParameter - m_startLength;
if (m_vect) {
double outVal = 0;
if (fabs(diff) <= m_lengthOfDeformation && m_versus == INNER) {
diff *= (1.0 / m_lengthOfDeformation) * c_maxLenghtOfGaussian;
outVal = gaussianPotential(diff);
} else if (m_versus == OUTER) {
double valForGaussian =
-c_maxLenghtOfGaussian + 2 * c_maxLenghtOfGaussian / m_lengthOfDeformation * strokeLengthAtParameter;
outVal = 1.0 - gaussianPotential(valForGaussian);
}
TPointD cp = convert(s.getControlPoint(n));
TPointD p = cp;
TRotation rot(*m_vect, outVal * rad2degree(m_angle));
p = rot * p;
return TThickPoint(p - cp, 0.0);
}
return TThickPoint();
}
TThickPoint TStrokeBenderDeformation::getDisplacementForControlPointLen(const TStroke &stroke, double cpLen) const
{
assert(0);
return TThickPoint();
}
//-----------------------------------------------------------------------------
TThickPoint TStrokeBenderDeformation::getDisplacement(const TStroke &s, double w) const
{
// potenziale exp^(-x^2) limitato tra [-c_maxLenghtOfGaussian,c_maxLenghtOfGaussian]
double strokeLengthAtParameter = s.getLength(w);
double diff = strokeLengthAtParameter - m_startLength;
if (m_vect) {
double outVal = 0.0;
if (fabs(diff) <= m_lengthOfDeformation)
if (m_versus == INNER) {
diff *= (1.0 / m_lengthOfDeformation) * c_maxLenghtOfGaussian;
outVal = gaussianPotential(diff);
} else if (m_versus == OUTER) {
double valForGaussian =
-c_maxLenghtOfGaussian + 2 * c_maxLenghtOfGaussian / m_lengthOfDeformation * strokeLengthAtParameter;
outVal = 1.0 - gaussianPotential(valForGaussian);
}
TPointD cp = convert(s.getControlPointAtParameter(w));
TPointD p = cp;
TRotation rot(*m_vect, outVal * rad2degree(m_angle));
p = rot * p;
return TThickPoint(p - cp, 0.0);
}
return TThickPoint();
}
//-----------------------------------------------------------------------------
double TStrokeBenderDeformation::getDelta(const TStroke &s, double w) const
{
double totalLenght = s.getLength();
if (totalLenght != 0) {
double val = s.getLength(w) / totalLenght * TConsts::pi * 10.0;
return sin(val);
}
return 0;
}
//=============================================================================
TStrokeTwirlDeformation::TStrokeTwirlDeformation(const TPointD ¢er, double radius)
: m_center(center), m_innerRadius2(sq(radius)), m_vectorOfMovement(TPointD())
{
m_outerRadius = 1.25 * radius;
}
//-----------------------------------------------------------------------------
TStrokeTwirlDeformation::TStrokeTwirlDeformation(const TPointD ¢er, double radius, const TPointD &v)
: m_center(center), m_innerRadius2(sq(radius)), m_vectorOfMovement(v)
{
m_outerRadius = 1.25 * radius;
}
//-----------------------------------------------------------------------------
TStrokeTwirlDeformation::~TStrokeTwirlDeformation()
{
}
//-----------------------------------------------------------------------------
TThickPoint TStrokeTwirlDeformation::getDisplacement(const TStroke &stroke, double s) const
{
double outVal = 0;
double distance2 = tdistance2(convert(stroke.getControlPointAtParameter(s)), m_center);
if (distance2 <= m_innerRadius2)
outVal = wyvillPotential(distance2, m_innerRadius2);
return TThickPoint(m_vectorOfMovement * outVal, 0);
}
//-----------------------------------------------------------------------------
double TStrokeTwirlDeformation::getDelta(const TStroke &stroke, double s) const
{
/*
vector<DoublePair> vres;
if(intersect( stroke, m_center, m_outerRadius, vres))
{
double totalLenght = stroke.getLength();
if(totalLenght != 0)
{
double val = stroke.getLength(s)/totalLenght * TConsts::pi *11.0;
return sin(val);
}
}
*/
return 0;
}
//-----------------------------------------------------------------------------
double TStrokeTwirlDeformation::getMaxDiff() const
{
return 0.4;
}
//-----------------------------------------------------------------------------
// End Of File
//-----------------------------------------------------------------------------
//=============================================================================
TStrokeThicknessDeformation::TStrokeThicknessDeformation(const TStroke *ref,
const TPointD &vect,
double s,
double l,
double versus)
: m_lengthOfDeformation(l), m_startParameter(s), m_versus(versus), m_vect(new TPointD(vect)), m_pRef(ref)
{
assert(m_lengthOfDeformation >= 0);
if (isAlmostZero(m_lengthOfDeformation))
m_lengthOfDeformation = TConsts::epsilon;
}
//-----------------------------------------------------------------------------
TStrokeThicknessDeformation::TStrokeThicknessDeformation(const TStroke *ref,
double s,
double l)
: m_lengthOfDeformation(l), m_startParameter(s), m_vect(0), m_pRef(ref)
{
assert(m_lengthOfDeformation >= 0);
if (isAlmostZero(m_lengthOfDeformation))
m_lengthOfDeformation = TConsts::epsilon;
}
//-----------------------------------------------------------------------------
TThickPoint TStrokeThicknessDeformation::getDisplacementForControlPoint(const TStroke &stroke, UINT n) const
{
// potenziale exp^(-x^2) limitato tra [-c_maxLenghtOfGaussian,c_maxLenghtOfGaussian]
double diff = stroke.getLengthAtControlPoint(n);
diff = diff - m_startParameter;
if (fabs(diff) <= m_lengthOfDeformation) {
// modulo il vettore spostamento in funzione del potenziale
// il termine (1.0/m_lengthOfDeformation) * 3 scala
// il punto in diff su un sistema di coordinate
// normalizzato, associato con la curva exp^(-x^2)
diff *= (1.0 / m_lengthOfDeformation) * c_maxLenghtOfGaussian;
TThickPoint delta;
if (m_vect) {
//tsign(m_vect->y) * 0.1
delta = TThickPoint(0, 0, m_versus * norm(*m_vect) * gaussianPotential(diff));
} else {
double outVal = gaussianPotential(diff);
delta = TThickPoint(0, 0, outVal);
}
//TThickPoint cp = stroke.getControlPoint(n);
//if(cp.thick + delta.thick<0.001) delta.thick = 0.001-cp.thick;
return delta;
}
return TThickPoint();
}
TThickPoint TStrokeThicknessDeformation::getDisplacementForControlPointLen(const TStroke &stroke, double diff) const
{
if (fabs(diff) <= m_lengthOfDeformation) {
// modulo il vettore spostamento in funzione del potenziale
// il termine (1.0/m_lengthOfDeformation) * 3 scala
// il punto in diff su un sistema di coordinate
// normalizzato, associato con la curva exp^(-x^2)
diff *= (1.0 / m_lengthOfDeformation) * c_maxLenghtOfGaussian;
TThickPoint delta;
if (m_vect) {
//tsign(m_vect->y) * 0.1
delta = TThickPoint(0, 0, m_versus * norm(*m_vect) * gaussianPotential(diff));
} else {
double outVal = gaussianPotential(diff);
delta = TThickPoint(0, 0, outVal);
}
//TThickPoint cp = stroke.getControlPoint(n);
//if(cp.thick + delta.thick<0.001) delta.thick = 0.001-cp.thick;
return delta;
}
return TThickPoint();
}
TThickPoint TStrokeThicknessDeformation::getDisplacement(const TStroke &stroke, double w) const
{
// potenziale exp^(-x^2) limitato tra [-c_maxLenghtOfGaussian,c_maxLenghtOfGaussian]
double diff = stroke.getLength(w);
diff = diff - m_startParameter;
if (fabs(diff) <= m_lengthOfDeformation) {
// modulo il vettore spostamento in funzione del potenziale
// il termine (1.0/m_lengthOfDeformation) * 3 scala
// il punto in diff su un sistema di coordinate
// normalizzato, associato con la curva exp^(-x^2)
diff *= (1.0 / m_lengthOfDeformation) * c_maxLenghtOfGaussian;
if (m_vect)
return gaussianPotential(diff) * TThickPoint(*m_vect, 0);
else {
double outVal = gaussianPotential(diff);
return TThickPoint(0, 0, outVal);
}
}
return TThickPoint();
}
//-----------------------------------------------------------------------------
double TStrokeThicknessDeformation::getDelta(const TStroke &stroke, double w) const
{
// potenziale exp^(-x^2) limitato tra [-c_maxLenghtOfGaussian,c_maxLenghtOfGaussian]
double diff = stroke.getLength(w);
diff = diff - m_startParameter;
if (fabs(diff) <= m_lengthOfDeformation) {
// modulo il vettore spostamento in funzione del potenziale
// il termine (1.0/m_lengthOfDeformation) * 3 scala
// il punto in diff su un sistema di coordinate
// normalizzato, associato con la curva exp^(-x^2)
diff *= (1.0 / m_lengthOfDeformation) * c_maxLenghtOfGaussian;
return derivateOfGaussianPotential(diff);
}
return 0;
}
//-----------------------------------------------------------------------------
double TStrokeThicknessDeformation::getMaxDiff() const
{
return 0.09;
}
//-----------------------------------------------------------------------------
TStrokeThicknessDeformation::~TStrokeThicknessDeformation()
{
delete m_vect;
}
//=============================================================================
TPointDeformation::TPointDeformation(const TStroke *stroke,
const TPointD ¢er,
double radius)
: m_strokeRef(stroke), m_center(center), m_radius(radius)
{
assert(m_strokeRef);
}
//-----------------------------------------------------------------------------
TPointDeformation::TPointDeformation()
{
}
//-----------------------------------------------------------------------------
TPointDeformation::~TPointDeformation()
{
}
//-----------------------------------------------------------------------------
TThickPoint TPointDeformation::getDisplacement(double s) const
{
// riferimento ad un punto ciccione della stroke
TThickPoint
thickPnt = m_strokeRef->getPointAtLength(s);
assert(thickPnt != TConsts::natp);
TPointD pntOfStroke(convert(thickPnt));
double d = tdistance(pntOfStroke, m_center);
double outVal = wyvillPotential(d, m_radius);
return TThickPoint(outVal, outVal, 0);
}
//-----------------------------------------------------------------------------
// return ratio: (number of Control Point)/length
double TPointDeformation::getCPDensity(double s) const
{
// reference to a thickpoint
TThickPoint
thickPnt = m_strokeRef->getThickPointAtLength(s);
assert(thickPnt != TConsts::natp);
TPointD
pntOfStroke = convert(thickPnt);
double d = tdistance(pntOfStroke, m_center);
return normOfGradientOfWyvillPotential(d, m_radius);
}
//-----------------------------------------------------------------------------
double TPointDeformation::getCPCountInRange(double s0, double s1) const
{
if (s1 < s0)
swap(s1, s0);
double step = (s1 - s0) * 0.1;
double cpCount = 0.0;
for (double s = s0; s < s1; s += step)
cpCount += getCPDensity(s);
cpCount += getCPDensity(s1);
return cpCount;
}
//-----------------------------------------------------------------------------
double TPointDeformation::getMinCPDensity() const
{
return 0.3;
}
//-----------------------------------------------------------------------------
// End Of File
//-----------------------------------------------------------------------------