//-----------------------------------------------------------------------------

//  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()

{

	delete m_imp;

}

//-----------------------------------------------------------------------------

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;</doublepair>

  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

//-----------------------------------------------------------------------------