Имеется замкнутый контур G, представленный системой n точек. Точки делят кривую на неравные по длине дуги части. Нужно найти новую систему точек, которая делит кривую G на равные по длине дуги части.

Воспользовался линейной интерполяцией y1[m]=(y[k+1]-y[k])/(x[k+1]-x[k])*(x[m]-x[k])+y[k].

Но, на картинке наблюдается "заусенец", при многократном переразбиении меняющегося контура G,
то есть что-то не работает в нижеприведенном коде, не могу найти что:

void interp(std::vector<point> & G){
  size_t n = G.size();
  double L = 0;
  std::vector<double> Delta_L(n+1);
  std::vector<point> G1(n);
  // Вспомагательный массив Delta_L
  Delta_L[0] = 0;
  for( size_t m=1; m<n; ++m){
       L += sqrt( pow(G[m].x-G[m-1].x,2.0) + pow(G[m].y-G[m-1].y,2.0) );
       Delta_L[m] = L;
  }
  L += sqrt( pow(G[0].x-G[n-1].x,2.0) + pow(G[0].y-G[n-1].y,2.0) );
  Delta_L[n] = L;
  // Равномерный шаг
  double dl = L / n;
  // Первые точки совпадают
  G1[0].x = G[0].x;
  G1[0].y = G[0].y;
  // Вычисляем оставшиеся точки 1,2...n-1
  for( size_t m=1; m<n; ++m ){
       for( size_t k=0; k<n; ++k ){
             //ищю между какими x[k]..x[k+1] должен быть x[m]
            if( Delta_L[k]<m*dl && m*dl<=Delta_L[k+1] ){
                double temp = ( m*dl-Delta_L[k] ) /
                          sqrt(1+pow((G[k+1].y-G[k].y)/(G[k+1].x-G[k].x),2));
                if( G[k].x == G[k+1].x )
                    G1[m].x = G[k].x;
                else if( G[k].x < G[k+1].x )
                    G1[m].x = temp + G[k].x;
                else
                    G1[m].x = - temp + G[k].x;
                G1[m].y = ( G[k+1].y-G[k].y ) / ( G[k+1].x-G[k].x ) *
                      ( G1[m].x-G[k].x ) + G[k].y;
            }
       }
  }
  G.swap(G1);
}