QGpCoreTools/Curve.h

Go to the documentation of this file.

/***************************************************************************
**
**  This file is part of QGpCoreTools.
**
**  This library is free software; you can redistribute it and/or
**  modify it under the terms of the GNU Lesser General Public
**  License as published by the Free Software Foundation; either
**  version 2.1 of the License, or (at your option) any later version.
**
**  This file 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 Lesser General Public
**  License for more details.
**
**  You should have received a copy of the GNU Lesser General Public
**  License along with this library; if not, write to the Free Software
**  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
**
**  See http://www.geopsy.org for more information.
**
**  Created : 2006-09-07
**  Authors :
**    Marc Wathelet
**    Marc Wathelet (LGIT, Grenoble, France)
**
***************************************************************************/

#ifndef CURVE_H
#define CURVE_H

#include <math.h>

#include "Global.h"
#include "Trace.h"
#include "Point2D.h"
#include "SmoothingParameters.h"
#include "CurvePointOptions.h"
#include "Translations.h"
#include "CoreApplication.h"

#include "QGpCoreToolsDLLExport.h"

namespace QGpCoreTools {

enum SamplingOption {LinearScale=0, LogScale=1, InversedScale=2, Interpole=4, Function=8};
Q_DECLARE_FLAGS(SamplingOptions, SamplingOption)

template <class pointType>
class QGPCORETOOLS_EXPORT Curve : private QVector<pointType>
{
  TRANSLATIONS("Curve");
public:
  Curve() {}
  Curve(int n) : QVector<pointType>(n) {}
  Curve(const Curve<pointType>& o) : QVector<pointType>(o) {}
  Curve(const QVector<pointType>& o) : QVector<pointType>(o) {}

  void setFunction();
  void unique();

  // Vector functions
  bool operator==(const Curve<pointType>& o) const {return QVector<pointType>::operator==(o);}
  int count() const {return QVector<pointType>::count();}
  bool isEmpty() const {return QVector<pointType>::isEmpty();}
  void resize(int n) {QVector<pointType>::resize(n);}
  int indexAfter(double valX) const;
  int indexOf(const pointType& p) const {return QVector<pointType>::indexOf(p);}
  pointType at(double x) const;
  const pointType& at(int index) const {return QVector<pointType>::at(index);}
  pointType& operator[](int index) {return QVector<pointType>::operator[](index);}
  const pointType& operator[](int index) const {return QVector<pointType>::at(index);}
  const pointType& first() const {return QVector<pointType>::first();}
  const pointType& last() const {return QVector<pointType>::last();}
  pointType& first() {return QVector<pointType>::first();}
  pointType& last() {return QVector<pointType>::last();}
  void remove(int index) {QVector<pointType>::remove(index);}
  void insert(int index, const pointType& p) {QVector<pointType>::insert(index, p);}
  void insert(const pointType& p);
  void clear() {QVector<pointType>::clear();}
  void prepend(const pointType& p) {QVector<pointType>::prepend(p);}
  void append(const pointType& p) {QVector<pointType>::append(p);}
  void append(const Curve<pointType>& c) {QVector<pointType>::operator+=(c);}
  typedef pointType * iterator;
  typedef const pointType * const_iterator;
  const_iterator begin() const {return QVector<pointType>::begin();}
  const_iterator end() const {return QVector<pointType>::end();}
  iterator begin() {return QVector<pointType>::begin();}
  iterator end() {return QVector<pointType>::end();}

  bool sameSampling(const Curve<pointType>& o) const;

  void xLog10();
  void xExp10();
  void xMultiply(double factor, SamplingOptions options=Function);
  void xInverse(SamplingOptions options=Function);

  void ySetValue(double value, const CurvePointOptions * pointOptions=0);
  void ySetMinimumValue(double value, const CurvePointOptions * pointOptions=0);
  void ySetMaximumValue(double value, const CurvePointOptions * pointOptions=0);
  void ySum(double value, const CurvePointOptions * pointOptions=0);
  void ySum(const Curve<pointType>& o, const CurvePointOptions * pointOptions=0);
  void yMultiply(double factor, const CurvePointOptions * pointOptions=0);
  void yMultiply(const Curve<pointType>& o, const CurvePointOptions * pointOptions=0);
  void ySquare(const CurvePointOptions * pointOptions=0);
  void ySqrt(const CurvePointOptions * pointOptions=0);
  void yInverse(const CurvePointOptions * pointOptions=0);
  void yLog10(const CurvePointOptions * pointOptions=0);
  void yExp10(const CurvePointOptions * pointOptions=0);

  void line(const pointType& p1, const pointType& p2);
  void resample(int n, double min, double max, SamplingOptions options,
                double valX, double valY, const CurvePointOptions * pointOptions=0);
  void resample(int n, double min, double max, SamplingOptions options);
  void resample(const Curve<pointType>& xModel, SamplingOptions options);
  void resample(const QVector<double>& xModel, SamplingOptions options);
  void cut(double min, double max, SamplingOptions options);
  void smooth(double min, double max, const SmoothingParameters& param);
  void average(const Curve<pointType>& o);

  bool leastSquare(double& a, double& b, const CurvePointOptions * pointOptions=0) const;
  double azimuth(const CurvePointOptions * pointOptions=0) const;
  QVector<double> project(const CurvePointOptions * pointOptions=0) const;

  void swapXY(const CurvePointOptions * pointOptions=0);
  Curve<pointType> derivative(const CurvePointOptions * pointOptions=0) const;

  void setValid(bool v);
  void removeInvalid();
  int minimumX(int startIndex=0) const;
  int maximumX(int startIndex=0) const;
  int minimumY(int startIndex=0, const CurvePointOptions * pointOptions=0) const;
  int maximumY(int startIndex=0, const CurvePointOptions * pointOptions=0) const;
  int closestMax(int index);

  QVector<double> xVector() const;
  QString toString(int precision=6, char format='g') const;

  void operator+=(const Curve<pointType>& o);
  void operator*=(const Curve<pointType>& o);
private:
  inline static bool lessThan(const pointType& p1, const pointType& p2);
  double sumX() const;
  double sumY(const CurvePointOptions * pointOptions=0) const;
  double sumXY(const CurvePointOptions * pointOptions=0) const;
  double sumX2() const;
  double sumY2(const CurvePointOptions * pointOptions=0) const;
};

template <class pointType>
void Curve<pointType>::setFunction()
{
  if(!isEmpty()) {
    qSort(begin(), end(), lessThan);
    // In case of lot of double entries, this is function was quite long with a remove()
    Curve<pointType> f;
    f.reserve(f.count());
    f.append(at(0));
    for(int i=1; i<count(); i++) {
      if(at(i-1).x()!=at(i).x()) {
        f.append(at(i));
      }
    }
    f.squeeze();
    // This is quite fast with implicit sharing
    *this=f;
  }
}

template <class pointType>
void Curve<pointType>::unique()
{
  ::QGpCoreTools::unique(*this);
}

template <class pointType>
inline bool Curve<pointType>::lessThan(const pointType& p1, const pointType& p2)
{
  return p1.x()<p2.x();
}

template <class pointType>
int Curve<pointType>::indexAfter (double valX) const
{
  TRACE;
  ASSERT(count()>1);
  if(valX <= at(0).x()) return 0;
  int lasti=count()-1;
  if(valX > at(lasti).x()) return lasti+1;
  int i=1;
  while(i < count()) i=i << 1;
  int step2=i >> 1;
  while(step2>0) {
    if(i>lasti) i-=step2;
    else if(valX <= at(i).x()) {
      if(valX > at(i-1).x()) break;
      i-=step2;
    } else i+=step2;
    step2=step2 >> 1;
  }
  return i;
}

template <class pointType>
pointType Curve<pointType>::at(double x) const
{
  if(count()<2) {
    if(count()==1) {
      return at(0);
    } else {
      return pointType();
    }
  } else {
    int i=indexAfter(x);
    if(i==0) i++;                // Force interpolation if outside range
    else if(i==count()) i--;
    pointType p;
    p.interpole(x, at(i-1), at(i));
    return p;
  }
}

template <class pointType>
void Curve<pointType>::insert(const pointType& p)
{
  TRACE;
  if(count()<2) {
    if(isEmpty()) {
      append(p);
    } else if(p.x() < at(0).x()) {
      prepend(p);
    } else if(p.x() > at(0).x()) {
      append(p);
    } else {
      (*this)[0]=p;
    }
  } else {
    int i=indexAfter(p.x());
    if(i<count() && at(i).x()==p.x()) {
      (*this)[i]=p;
    } else {
      if(at(p.x()).y()!=p.y()) {
        insert(i, p);
      }
    }
  }
}

template <class pointType>
void Curve<pointType>::resample(int n, double min, double max, SamplingOptions options, double valX, double valY,
                                const CurvePointOptions * pointOptions)
{
  TRACE;
  // Transforms X scale according to options (log and inv)
  if(options & InversedScale) {
    xInverse(options);
    valX=1.0/valX;
  }
  if(options & LogScale) {
    xLog10();
    valX=::log10(valX);
    if(options & Function) {
      min=::log10(min);
      max=::log10(max);
    }
  }
  if(options & Function) {
    if(min < first().x())
      min=first().x();
    if(max > last().x())
      max=last().x();
  }
  int currentN=count();

  // Calculate the total "distance" curve
  Curve<Point2D> distances(currentN);
  distances[ 0 ].setX(at(0).x());
  distances[ 0 ].setY(0.0);
  for(int i=1; i < currentN; i++ ) {
    const pointType& p1=at(i - 1);
    const pointType& p2=at(i);
    distances[ i ].setX(p2.x());
    double deltaX=(p2.x() - p1.x())/valX;
    double deltaY=(p2.y(pointOptions) - p1.y(pointOptions))/valY;
    distances[ i ].setY(distances[ i - 1 ].y() + ::sqrt(deltaX * deltaX + deltaY * deltaY));
  }
  distances.setFunction();

  // Calculate the distance of min and max, and constant step
  double distMin, distMax;
  if(options & Function) {
    Point2D p;
    int i;
    i=distances.indexAfter(min);
    p.interpole(min, distances.at(i-1), distances.at(i));
    distMin=p.y();
    i=distances.indexAfter(max);
    p.interpole(max, distances.at(i-1), distances.at(i));
    distMax=p.y();
  } else {
    distMin=0;
    distMax=distances[ currentN - 1 ].y();
  }
  double cstStep=(distMax - distMin)/(currentN - 1);

  // Map x as a function of distance, function monotoneous, then no need to setFunction()
  distances.swapXY();

  // Intepolate x for constant step distance
  QVector<double> x(n);
  double dist=distMin;
  Point2D p;
  for(int i=0; i < n; i++, dist += cstStep) {
    int iDist=distances.indexAfter(dist  );
    p.interpole(dist, distances.at(iDist-1), distances.at(iDist));
    x[ i ]=p.y();
  }
  resample(x, options);

  // Re-Transforms axis according options (log and inv)
  if(options & LogScale)
    xExp10();
  if(options & InversedScale)
    xInverse(options);
}

template <class pointType>
void Curve<pointType>::resample(int n, double min, double max, SamplingOptions options)
{
  TRACE;
  ASSERT(options & Function);
  // Transforms X scale according to options (log and inv)
  if(options & InversedScale) {
    xInverse(options);
  }
  if(options & LogScale) {
    xLog10();
    min=::log10(min);
    max=::log10(max);
  }
  // Calculate the constant step
  double cstStep=(max - min)/(n - 1);
  QVector<double> x(n);
  for(int i=0;i < n;i++ )
    x[ i ]=min + (double)i*cstStep;
  resample(x, options);

  // Re-Transforms axis according options (log and inv)
  if(options & LogScale)
    xExp10();
  if(options & InversedScale)
    xInverse(options);
}

template <class pointType>
void Curve<pointType>::resample(const Curve<pointType>& xModel, SamplingOptions options)
{
  TRACE;
  ASSERT(options & Function);
  QVector<double> x(xModel.count());
  for(int i=0;i < xModel.count();i++ )
    x[ i ]=xModel.at(i).x();
  resample(x, options);
}

template <class pointType>
void Curve<pointType>::resample(const QVector<double>& xModel, SamplingOptions options)
{
  TRACE;
  ASSERT(options & Function);
  setFunction(); // Make sure values are sorted
  int nCurve=count();
  int nModel=xModel.count();
  if(nCurve<2 || nModel<2) return;
  Curve<pointType> interpolated(nModel);
  // Insert values inside xModel
  int iCurve=1;
  double min=at(0).x()*(1+(at(0).x()<0 ? 1e-15 : -1e-15));
  double max=at(nCurve-1).x()*(1+(at(nCurve-1).x()<0 ? -1e-15 : 1e-15));
  for(int i =0; i<nModel;i++) {
    double x=xModel.at(i);
    while(iCurve<nCurve-1 && at(iCurve).x()<x) iCurve++;
    interpolated[i].interpole(x, at(iCurve-1), at(iCurve));
    if(x<min || x>max || !at(iCurve-1).isValid() || !at(iCurve).isValid()) interpolated[i].setValid(false);
  }
  // Insert all values outside xModel
  if(nModel>1) {
    min=xModel.at(0)*(1+(xModel.at(0)<0 ? 1e-15 : -1e-15));
    max=xModel.at(nModel-1)*(1+(xModel.at(nModel-1)<0 ? -1e-15 : 1e-15));
    for(int i=0;i<nCurve;i++) {
      const pointType& p=at(i);
      if(p.x()<min || max<p.x()) {
        interpolated.append(p);
      }
    }
  }
  interpolated.setFunction();
  *this=interpolated;
}

template <class pointType>
void Curve<pointType>::cut(double min, double max, SamplingOptions options)
{
  TRACE;
  ASSERT(options & Function);
  // Assert that count()>1 is made by indexAfter()
  // Transforms X scale according to options (log and inv)
  if(options & InversedScale) {
    xInverse(options);
  }
  if(min < first().x())
    min=first().x();
  if(max > last().x())
    max=last().x();
  // Find the first point
  int ix=indexAfter(min);
  if(ix>0) {
    if(fabs(min - at(ix-1).x()) < 1e-15) {
      ix--;
    } else if((options & Interpole) && fabs(min - at(ix).x()) > 1e-15) {
      ix--;
      (*this)[ ix ].interpole(min, at(ix), at(ix+1));
    }
  }
  // Remove everything before
  for(int i=0; i<ix; i++) remove(0);
  // Find the last point
  ix=indexAfter(max);
  if(ix < count()) {
    if(fabs(max - at(ix).x()) < 1e-15) {
      ix++;
    } else if((options & Interpole) && fabs(max - at(ix-1).x()) > 1e-15) {
      (*this)[ ix ].interpole(max, at(ix-1), at(ix));
      ix++;
    }
  }
  // Remove everything after
  for(int i=count()-1; i>=ix; i--) remove(i);
  // Re-Transforms axis according options (log and inv)
  if(options & InversedScale)
    xInverse(options);
}

template <class pointType>
void Curve<pointType>::smooth(double min, double max, const SmoothingParameters& param)
{
  TRACE;
  
}

template <class pointType>
void Curve<pointType>::average(const Curve<pointType>& o)
{
  TRACE;
  if(o.count()<2) return;
  if(count()==0) {
    *this=o;
    return;
  }
  // Set the same sampling for both curves
  QVector<double> x=xVector();
  x+=o.xVector();
  qSort(x);
  ::QGpCoreTools::unique(x);
  resample(x, Function);
  Curve<pointType> oResamp=o;
  oResamp.resample(x, Function);
  // Take the average of the two curves
  for(int i=x.count()-1; i>=0; i--) {
    operator[](i).average(oResamp.at(i));
  }
  setFunction();
}

template <class pointType>
void Curve<pointType>::setValid(bool v)
{
  TRACE;
  iterator it;
  for(it=begin(); it!=end(); ++it) {
    it->setValid(v);
  }
}

template <class pointType>
void Curve<pointType>::removeInvalid()
{
  TRACE;
  for(int i=count() - 1; i>=0; i-- ) {
    if( !at(i).isValid()) {
      remove(i);
    }
  }
}

template <class pointType>
void Curve<pointType>::swapXY(const CurvePointOptions * pointOptions)
{
  TRACE;
  iterator it;
  for(it=begin(); it!=end(); ++it) {
    double tmp=it->x();
    it->setX(it->y(pointOptions));
    it->setY(tmp, pointOptions);
  }
}

template <class pointType>
Curve<pointType> Curve<pointType>::derivative(const CurvePointOptions * pointOptions) const
{
  TRACE;
  int n=count()-1;
  ASSERT(n>1);
  Curve<pointType> d(n-1);
  for(int i=1; i<n; i++ ) {
    double v=(at(i).y(pointOptions)-at(i-1).y(pointOptions))
               /(at(i).x()-at(i-1).x());
    v+=(at(i+1).y(pointOptions)-at(i).y(pointOptions))
      /(at(i+1).x()-at(i).x());
    v*=2.0;
    pointType& p=d[i-1];
    p.setX(at(i).x());
    p.setY(v, pointOptions);
  }
  return d;
}

template <class pointType>
void Curve<pointType>::xMultiply(double factor, SamplingOptions options)
{
  TRACE;
  iterator it;
  for(it=begin(); it!=end(); ++it) {
    pointType& p=*it;
    p.setX(p.x()*factor);
  }
  if(factor<0.0 && options & Function) setFunction();
}

template <class pointType>
void Curve<pointType>::xInverse(SamplingOptions options)
{
  TRACE;
  iterator it;
  for(it=begin(); it!=end(); ++it) {
    pointType& p=*it;
    p.setX(1.0/p.x());
  }
  if(options & Function) setFunction();
}

template <class pointType>
void Curve<pointType>::xExp10()
{
  TRACE;
  iterator it;
  for(it=begin(); it!=end(); ++it) {
    pointType& p=*it;
    p.setX(pow(10, p.x()));
  }
}

template <class pointType>
void Curve<pointType>::xLog10()
{
  TRACE;
  iterator it;
  for(it=begin(); it!=end(); ++it) {
    pointType& p=*it;
    if(p.x()>0.0) p.setX(::log10(p.x()));
  }
}

template <class pointType>
void Curve<pointType>::ySetValue(double v, const CurvePointOptions * pointOptions)
{
  int n=count();
  for(int i=0; i<n; i++) {
    (*this)[i].setY(v, pointOptions);
  }
}

template <class pointType>
void Curve<pointType>::ySetMinimumValue(double v, const CurvePointOptions * pointOptions)
{
  int n=count();
  for(int i=0; i<n; i++) {
    pointType& p=(*this)[i];
    if(p.y(pointOptions)<v) {
      p.setY(v, pointOptions);
    }
  }
}

template <class pointType>
void Curve<pointType>::ySetMaximumValue(double v, const CurvePointOptions * pointOptions)
{
  int n=count();
  for(int i=0; i<n; i++) {
    pointType& p=(*this)[i];
    if(p.y(pointOptions)>v) {
      p.setY(v, pointOptions);
    }
  }
}

template <class pointType>
void Curve<pointType>::ySum(double value, const CurvePointOptions * pointOptions)
{
  TRACE;
  iterator it;
  for(it=begin(); it!=end(); ++it) {
    pointType& p=*it;
    p.setY(p.y(pointOptions)+value, pointOptions);
  }
}

template <class pointType>
void Curve<pointType>::ySum(const Curve<pointType>& o, const CurvePointOptions * pointOptions)
{
  if(isEmpty()) {
    *this=o;
    return;
  }
  int n=count();
  if(n!=o.count()) {
    App::stream() << tr("Curve::ySum(): uncompatible sampling.") << endl;
    return;
  }
  for(int i=0; i<n; i++) {
    pointType& p=(*this)[i];
    p.setY(p.y(pointOptions)+o.at(i).y(pointOptions), pointOptions);
  }
}

template <class pointType>
void Curve<pointType>::yMultiply(double factor, const CurvePointOptions * pointOptions)
{
  TRACE;
  iterator it;
  for(it=begin(); it!=end(); ++it) {
    pointType& p=*it;
    p.setY(p.y(pointOptions)*factor, pointOptions);
  }
}

template <class pointType>
void Curve<pointType>::yMultiply(const Curve<pointType>& o, const CurvePointOptions * pointOptions)
{
  if(isEmpty()) {
    *this=o;
    return;
  }
  int n=count();
  if(n!=o.count()) {
    App::stream() << tr("Curve::yMultiply(): uncompatible sampling.") << endl;
    return;
  }
  for(int i=0; i<n; i++) {
    pointType& p=(*this)[i];
    p.setY(p.y(pointOptions)*o.at(i).y(pointOptions), pointOptions);
  }
}

template <class pointType>
void Curve<pointType>::ySquare(const CurvePointOptions * pointOptions)
{
  int n=count();
  for(int i=0; i<n; i++) {
    pointType& p=(*this)[i];
    double v=p.y(pointOptions);
    p.setY(v*v, pointOptions);
  }
}

template <class pointType>
void Curve<pointType>::ySqrt(const CurvePointOptions * pointOptions)
{
  int n=count();
  for(int i=0; i<n; i++) {
    pointType& p=(*this)[i];
    double v=p.y(pointOptions);
    if(v>=0.0) p.setY(::sqrt(v), pointOptions);
  }
}

template <class pointType>
void Curve<pointType>::yInverse(const CurvePointOptions * pointOptions)
{
  TRACE;
  iterator it;
  for(it=begin(); it!=end(); ++it) {
    pointType& p=*it;
    double v=p.y(pointOptions);
    if(v!=0.0) p.setY(1.0/v, pointOptions);
  }
}

template <class pointType>
void Curve<pointType>::yLog10(const CurvePointOptions * pointOptions)
{
  TRACE;
  iterator it;
  for(it=begin(); it!=end(); ++it) {
    pointType& p=*it;
    double v=p.y(pointOptions);
    if(v>0.0) p.setY(::log10(v), pointOptions);
  }
}

template <class pointType>
void Curve<pointType>::yExp10(const CurvePointOptions * pointOptions)
{
  TRACE;
  iterator it;
  for(it=begin(); it!=end(); ++it) {
    pointType& p=*it;
    p.setY(pow(10, p.y(pointOptions)), pointOptions);
  }
}

template <class pointType>
void Curve<pointType>::line(const pointType& p1, const pointType& p2)
{
  TRACE;
  resize(2);
  (*this)[0]=p1;
  (*this)[1]=p2;
}

template <class pointType>
QVector<double> Curve<pointType>::xVector() const
{
  int n=count();
  QVector<double> x(n);
  for(int i=0; i<n; i++ ) {
    x[i]=at(i).x();
  }
  return x;
}

template <class pointType>
int Curve<pointType>::minimumX(int startIndex) const
{
  int n=count();
  ASSERT(n>0 && startIndex<n);
  if(startIndex<0) startIndex=0;
  const pointType * p=&at(startIndex);
  int iMin=startIndex;
  for(int i=startIndex+1; i<n; i++ ) {
    if(at(i).x() < p->x()) {
      p=&at(i);
      iMin=i;
    }
  }
  return iMin;
}

template <class pointType>
int Curve<pointType>::maximumX(int startIndex) const
{
  int n=count();
  ASSERT(n>0 && startIndex<n);
  if(startIndex<0) startIndex=0;
  const pointType * p=&at(startIndex);
  int iMax=startIndex;
  for(int i=startIndex+1; i<n; i++ ) {
    if(at(i).x() > p->x()) {
      p=&at(i);
      iMax=i;
    }
  }
  return iMax;
}

template <class pointType>
int Curve<pointType>::minimumY(int startIndex, const CurvePointOptions * pointOptions) const
{
  int n=count();
  ASSERT(n>0 && startIndex<n);
  if(startIndex<0) startIndex=0;
  const pointType * p=&at(startIndex);
  int iMin=startIndex;
  for(int i=startIndex+1; i<n; i++ ) {
    if(at(i).y(pointOptions) < p->y(pointOptions)) {
      p=&at(i);
      iMin=i;
    }
  }
  return iMin;
}

template <class pointType>
int Curve<pointType>::maximumY(int startIndex, const CurvePointOptions * pointOptions) const
{
  int n=count();
  ASSERT(n>0 && startIndex<n);
  if(startIndex<0) startIndex=0;
  const pointType * p=&at(startIndex);
  int iMax=startIndex;
  for(int i=startIndex+1; i<n; i++ ) {
    if(at(i).y(pointOptions) > p->y(pointOptions)) {
      p=&at(i);
      iMax=i;
    }
  }
  return iMax;
}

template <class pointType>
int Curve<pointType>::closestMax(int index)
{
  TRACE;
  int n=count();
  int di;

  if(index==0) {
    di=1;
    index=1;
  }
  else if(index>=n) {
    index=n;
    di=-1;
  } else if(at(index).y()>at(index-1).y()) di=1;
  else if(at(index).y()<at(index-1).y()) di=-1;
  else return index;

  if(di>0) {
    while(index<n) {
      if(at(index).y()<at(index-1).y()) {
        index--;
        break;
      }
      index+=di;
    }
    if(index==n) index--;
  } else {
    while(index>0) {
      if(at(index).y()>at(index-1).y()) break;
      index+=di;
    }
  }
  return index;
}

template <class pointType>
QString Curve<pointType>::toString(int precision, char format) const
{
  TRACE;
  QString s;
  for(const_iterator it=begin(); it!=end(); ++it) {
    s+=it->toString(precision, format);
    s+="\n";
  }
  return s;
}

template <class pointType>
double Curve<pointType>::sumX() const
{
  TRACE;
  double sum=0.0;
  const_iterator it;
  for(it=begin();it!=end();++it) sum+=it->x();
  return sum;
}

template <class pointType>
double Curve<pointType>::sumY(const CurvePointOptions * pointOptions) const
{
  TRACE;
  double sum=0.0;
  const_iterator it;
  for(it=begin();it!=end();++it) sum+=it->y(pointOptions);
  return sum;
}

template <class pointType>
double Curve<pointType>::sumXY(const CurvePointOptions * pointOptions) const
{
  TRACE;
  double sum=0.0;
  const_iterator it;
  for(it=begin();it!=end();++it) sum+=it->x()*it->y(pointOptions);
  return sum;
}

template <class pointType>
double Curve<pointType>::sumX2() const
{
  TRACE;
  double sum=0.0;
  const_iterator it;
  for(it=begin();it!=end();++it) sum+=it->x()*it->x();
  return sum;
}

template <class pointType>
double Curve<pointType>::sumY2(const CurvePointOptions * pointOptions) const
{
  TRACE;
  double sum=0.0;
  const_iterator it;
  for(it=begin();it!=end();++it) sum+=it->y()*it->y(pointOptions);
  return sum;
}

template <class pointType>
bool Curve<pointType>::leastSquare(double& a, double& b, const CurvePointOptions * pointOptions) const
{
  TRACE;
  int n=count();
  if(n>1) {
    double invn=1.0/(double)n;
    double sx=sumX();
    double sy=sumY(pointOptions);
    double sxy=sumXY(pointOptions);
    double sx2=sumX2();
    printf("sum x=%lg sum x*sum x=%lg nsum2 x=%lg\n", sx, sx*sx, n*sx2);
    double denom=sx2-invn*sx*sx;
    if(denom!=0.0) {
      a=(sxy-invn*sy*sx)/denom;
      b=(invn*sy*sx2-invn*sxy*sx)/denom;
      return true;
    } else {
      a=0.0;
      b=sx*invn;
    }
  }
  return false;
}

template <class pointType>
double Curve<pointType>::azimuth(const CurvePointOptions * pointOptions) const
{
  TRACE;
  int n=count();
  if(n>1) {
    double a, b;
    if(leastSquare(a, b, pointOptions)) {
      return atan(a);
    } else {
      return 0.5*M_PI;
    }
  }
  return 0.0;
}

template <class pointType>
QVector<double> Curve<pointType>::project(const CurvePointOptions * pointOptions) const
{
  TRACE;
  int n=count();
  if(!n) return QVector<double>();
  QVector<double> projections(n);
  if(n>1) {
    double a, b;
    double min=1e99;
    if(leastSquare(a, b, pointOptions)) {
      for(int i=0; i<n; i++) {
        const pointType& p=at(i);
        double denom=a*a+1.0;
        double cp=-p.x()-a*p.y(pointOptions);
        double xp=(-cp-a*b)/denom;
        double yp=(-a*cp+b)/denom-b;
        double d=::sqrt(xp*xp+yp*yp);
        projections[i]=(xp>=0)? d : -d;
        if(projections[i]<min) min=projections[i];
      }
    } else {
      for(int i=0; i<n; i++) {
        const pointType& p=at(i);
        projections[i]=p.y(pointOptions);
        if(projections[i]<min) min=projections[i];
      }
    }
    // set min to abcsisse=0
    for(int i=0; i<n;i++) projections[i]-=min;
  } else {
    projections[0]=0;
  }
  return projections;
}

template <class pointType>
bool Curve<pointType>::sameSampling(const Curve<pointType>& o) const
{
  TRACE;
  int n=count();
  if(n!=o.count()) {
    return false;
  }
  for(int i=0; i<n; i++) {
    if(at(i).x()!=o.at(i).x()) {
      return false;
    }
  }
  return true;
}

} // namespace QGpCoreTools

Q_DECLARE_OPERATORS_FOR_FLAGS(QGpCoreTools::SamplingOptions)

#endif // CURVE_H