mirror of https://github.com/CGAL/cgal
4008 lines
71 KiB
C++
4008 lines
71 KiB
C++
|
|
#ifndef CGAL_POLYAP_FCT
|
|
#define CGAL_POLYAP_FCT
|
|
|
|
|
|
#include <assert.h>
|
|
#include <CGAL/basic.h>
|
|
|
|
CGAL_BEGIN_NAMESPACE
|
|
|
|
|
|
//////////////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Dynamic Programming
|
|
|
|
// Bounded-# (min-E) local error assessment
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_DP_bnp_lea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t n_pt_bound,
|
|
typename DistTraits::FT &approx_error,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
typename std::size_t n,m,i;
|
|
|
|
InputIterator c_n,c_i,c_m;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_error=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
if(n_pt_bound<2)
|
|
return result;
|
|
|
|
for(n=0,c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
if(n_pt_bound>n)
|
|
{
|
|
approx_error=0;
|
|
|
|
for(c_i=begin;c_i!=beyond;c_i++)
|
|
*result++=*c_i;
|
|
|
|
return result;
|
|
}
|
|
|
|
typename std::size_t N=n+1;
|
|
|
|
FT **ApproxError;
|
|
typename std::size_t **SplitPoints;
|
|
FT **Err;
|
|
|
|
ApproxError=new FT*[N];
|
|
assert(ApproxError);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
ApproxError[m]=new FT[n_pt_bound];
|
|
assert(ApproxError[m]);
|
|
}
|
|
|
|
SplitPoints=new typename std::size_t*[N];
|
|
assert(SplitPoints);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
SplitPoints[m]=new typename std::size_t[n_pt_bound];
|
|
assert(SplitPoints[m]);
|
|
}
|
|
|
|
Err=new FT*[N];
|
|
assert(Err);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
Err[m]=new FT[N];
|
|
assert(Err[m]);
|
|
}
|
|
|
|
Err[0][0]=FT(0);
|
|
Err[N-1][N-1]=FT(0);
|
|
for(n=1;n<N-1;n++)
|
|
Err[n][n]=Err[n][n-1]=Err[n][n+1]=Err[n-1][n]=Err[n+1][n]=FT(0);
|
|
|
|
for(n=0,c_n=begin;n<N-2;n++,c_n++)
|
|
for(m=n+2,c_m=c_n,c_m++,c_m++,c_m++;m<N;m++,c_m++)
|
|
Err[m][n]=Err[n][m]=error(c_n,c_m);
|
|
|
|
if(n_pt_bound==2)
|
|
{
|
|
approx_error=Err[0][N-1];
|
|
*result++=*begin;
|
|
c_n=beyond;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
for(n=1;n<N;n++)
|
|
ApproxError[n][1]=Err[n][0];
|
|
|
|
FT err_crt;
|
|
|
|
for(m=2;m<n_pt_bound;m++)
|
|
{
|
|
ApproxError[m][m]=FT(0);
|
|
SplitPoints[m][m]=m-1;
|
|
|
|
for(n=m+1;n<N;n++)
|
|
{
|
|
ApproxError[n][m]=ApproxError[n-1][m-1];
|
|
SplitPoints[n][m]=n-1;
|
|
for(i=n-2;i>=m;i--)
|
|
{
|
|
if(ApproxError[i][m-1]>Err[i][n])
|
|
err_crt=ApproxError[i][m-1];
|
|
else
|
|
err_crt=Err[i][n];
|
|
|
|
if(err_crt<ApproxError[n][m])
|
|
{
|
|
ApproxError[n][m]=err_crt;
|
|
SplitPoints[n][m]=i;
|
|
}
|
|
}
|
|
err_crt=Err[m-1][n];
|
|
if(err_crt<ApproxError[n][m])
|
|
{
|
|
ApproxError[n][m]=err_crt;
|
|
SplitPoints[n][m]=m-1;
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
approx_error=ApproxError[N-1][n_pt_bound-1];
|
|
|
|
c_i=beyond;
|
|
c_i--;
|
|
|
|
typename std::list<InputIterator> local_container;
|
|
typename std::list<InputIterator>::iterator cri;
|
|
|
|
local_container.push_front(c_i);
|
|
|
|
n=N-1;
|
|
|
|
typename std::size_t nn=N-1;
|
|
|
|
for(m=n_pt_bound-1;m>=2;m--)
|
|
{
|
|
n=SplitPoints[n][m];
|
|
|
|
for(typename std::size_t ii=nn;ii>n;ii--)
|
|
c_i--;
|
|
nn=n;
|
|
|
|
local_container.push_front(c_i);
|
|
}
|
|
|
|
local_container.push_front(begin);
|
|
|
|
for(cri=local_container.begin();cri!=local_container.end();cri++)
|
|
*result++=**cri;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] ApproxError[m];
|
|
delete[] ApproxError;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] SplitPoints[m];
|
|
delete[] SplitPoints;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Err[m];
|
|
delete[] Err;
|
|
|
|
return result;
|
|
};
|
|
|
|
|
|
// Bounded-E (min-#) local error assessment
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_DP_be_lea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t &approx_n_pt,
|
|
typename DistTraits::FT error_bound,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
typename std::size_t n,m,i;
|
|
FT err_crt;
|
|
InputIterator c_n,c_i,c_m;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_n_pt=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=1;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=2;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
for(n=0,c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
typename std::size_t N=n+1;
|
|
|
|
FT **ApproxError;
|
|
typename std::size_t **SplitPoints;
|
|
FT **Err;
|
|
|
|
ApproxError=new FT*[N];
|
|
assert(ApproxError);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
ApproxError[m]=new FT[N];
|
|
assert(ApproxError[m]);
|
|
}
|
|
|
|
Err=new FT*[N];
|
|
assert(Err);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
Err[m]=new FT[N];
|
|
assert(Err[m]);
|
|
}
|
|
|
|
Err[0][0]=Err[0][1]=Err[1][0]=FT(0);
|
|
|
|
if(N>2)
|
|
for(n=2,c_n=begin,c_n++,c_n++,c_n++;n<N;n++,c_n++)
|
|
ApproxError[n][1]=Err[n][0]=Err[0][n]=error(begin,c_n);
|
|
|
|
if(ApproxError[N-1][1]<=error_bound)
|
|
{
|
|
approx_n_pt=2;
|
|
|
|
c_i=beyond;
|
|
c_i--;
|
|
|
|
*result++=*begin;
|
|
*result=*c_i;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] ApproxError[m];
|
|
delete[] ApproxError;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Err[m];
|
|
delete[] Err;
|
|
|
|
return result;
|
|
}
|
|
|
|
SplitPoints=new typename std::size_t*[N];
|
|
assert(SplitPoints);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
SplitPoints[m]=new typename std::size_t[N];
|
|
assert(SplitPoints[m]);
|
|
}
|
|
|
|
Err[N-1][N-1]=FT(0);
|
|
for(n=1;n<N-1;n++)
|
|
Err[n][n]=Err[n][n-1]=Err[n][n+1]=Err[n-1][n]=Err[n+1][n]=FT(0);
|
|
|
|
for(n=1,c_n=begin,c_n++;n<N-2;n++,c_n++)
|
|
for(m=n+2,c_m=c_n,c_m++,c_m++,c_m++;m<N;m++,c_m++)
|
|
Err[m][n]=Err[n][m]=error(c_n,c_m);
|
|
|
|
for(m=2;m<N;m++)
|
|
{
|
|
ApproxError[m][m]=FT(0);
|
|
SplitPoints[m][m]=m-1;
|
|
|
|
for(n=m+1;n<N;n++)
|
|
{
|
|
ApproxError[n][m]=ApproxError[n-1][m-1];
|
|
SplitPoints[n][m]=n-1;
|
|
for(i=n-2;i>=m;i--)
|
|
{
|
|
if(ApproxError[i][m-1]>Err[i][n])
|
|
err_crt=ApproxError[i][m-1];
|
|
else
|
|
err_crt=Err[i][n];
|
|
|
|
if(err_crt<ApproxError[n][m])
|
|
{
|
|
ApproxError[n][m]=err_crt;
|
|
SplitPoints[n][m]=i;
|
|
}
|
|
}
|
|
err_crt=Err[m-1][n];
|
|
if(err_crt<ApproxError[n][m])
|
|
{
|
|
ApproxError[n][m]=err_crt;
|
|
SplitPoints[n][m]=m-1;
|
|
}
|
|
}
|
|
|
|
if(ApproxError[N-1][m]<=error_bound)
|
|
{
|
|
approx_n_pt=m+1;
|
|
break;
|
|
}
|
|
}
|
|
|
|
c_i=beyond;
|
|
c_i--;
|
|
|
|
typename std::list<InputIterator> local_container;
|
|
typename std::list<InputIterator>::iterator cri;
|
|
|
|
local_container.push_front(c_i);
|
|
|
|
n=N-1;
|
|
|
|
typename std::size_t nn=N-1;
|
|
|
|
for(m=approx_n_pt-1;m>=2;m--)
|
|
{
|
|
n=SplitPoints[n][m];
|
|
|
|
for(typename std::size_t ii=nn;ii>n;ii--)
|
|
c_i--;
|
|
|
|
nn=n;
|
|
local_container.push_front(c_i);
|
|
}
|
|
|
|
local_container.push_front(begin);
|
|
|
|
for(cri=local_container.begin();cri!=local_container.end();cri++)
|
|
*result++=**cri;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] ApproxError[m];
|
|
delete[] ApproxError;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] SplitPoints[m];
|
|
delete[] SplitPoints;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Err[m];
|
|
delete[] Err;
|
|
|
|
return result;
|
|
};
|
|
|
|
|
|
// Bounded-# (min-E) global error assessment
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_DP_bnp_gea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t n_pt_bound,
|
|
typename DistTraits::FT &approx_error,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
typename std::size_t n,m,i;
|
|
|
|
InputIterator c_n,c_i,c_m;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_error=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
if(n_pt_bound<2)
|
|
return result;
|
|
|
|
n=0;
|
|
for(c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
if(n_pt_bound>n)
|
|
{
|
|
approx_error=0;
|
|
|
|
for(c_i=begin;c_i!=beyond;c_i++)
|
|
*result++=*c_i;
|
|
|
|
return result;
|
|
}
|
|
|
|
typename std::size_t N=n+1;
|
|
|
|
FT **ApproxError;
|
|
typename std::size_t **SplitPoints;
|
|
FT **Err;
|
|
|
|
ApproxError=new FT*[N];
|
|
assert(ApproxError);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
ApproxError[m]=new FT[n_pt_bound];
|
|
assert(ApproxError[m]);
|
|
}
|
|
|
|
SplitPoints=new typename std::size_t*[N];
|
|
assert(SplitPoints);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
SplitPoints[m]=new typename std::size_t[n_pt_bound];
|
|
assert(SplitPoints[m]);
|
|
}
|
|
|
|
Err=new FT*[N];
|
|
assert(Err);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
Err[m]=new FT[N];
|
|
assert(Err[m]);
|
|
}
|
|
|
|
Err[0][0]=FT(0);
|
|
Err[N-1][N-1]=FT(0);
|
|
for(n=1;n<N-1;n++)
|
|
Err[n][n]=Err[n][n-1]=Err[n][n+1]=Err[n-1][n]=Err[n+1][n]=FT(0);
|
|
|
|
for(n=0,c_n=begin;n<N-2;n++,c_n++)
|
|
for(m=n+2,c_m=c_n,c_m++,c_m++,c_m++;m<N;m++,c_m++)
|
|
Err[m][n]=Err[n][m]=error(c_n,c_m);;
|
|
|
|
if(n_pt_bound==2)
|
|
{
|
|
approx_error=Err[0][N-1];
|
|
*result++=*begin;
|
|
c_n=beyond;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
for(n=1;n<N;n++)
|
|
ApproxError[n][1]=Err[n][0];
|
|
|
|
FT err_crt;
|
|
|
|
for(m=2;m<n_pt_bound;m++)
|
|
{
|
|
ApproxError[m][m]=FT(0);
|
|
SplitPoints[m][m]=m-1;
|
|
|
|
for(n=m+1;n<N;n++)
|
|
{
|
|
ApproxError[n][m]=ApproxError[n-1][m-1];
|
|
SplitPoints[n][m]=n-1;
|
|
for(i=n-2;i>=m;i--)
|
|
{
|
|
err_crt=ApproxError[i][m-1]+Err[i][n];
|
|
if(err_crt<ApproxError[n][m])
|
|
{
|
|
ApproxError[n][m]=err_crt;
|
|
SplitPoints[n][m]=i;
|
|
}
|
|
}
|
|
err_crt=Err[m-1][n];
|
|
if(err_crt<ApproxError[n][m])
|
|
{
|
|
ApproxError[n][m]=err_crt;
|
|
SplitPoints[n][m]=m-1;
|
|
}
|
|
}
|
|
}
|
|
|
|
approx_error=ApproxError[N-1][n_pt_bound-1];
|
|
|
|
c_i=beyond;
|
|
c_i--;
|
|
|
|
|
|
typename std::list<InputIterator> local_container;
|
|
typename std::list<InputIterator>::iterator cri;
|
|
|
|
local_container.push_front(c_i);
|
|
|
|
n=N-1;
|
|
|
|
typename std::size_t nn=N-1;
|
|
|
|
for(m=n_pt_bound-1;m>=2;m--)
|
|
{
|
|
n=SplitPoints[n][m];
|
|
|
|
for(typename std::size_t ii=nn;ii>n;ii--)
|
|
c_i--;
|
|
nn=n;
|
|
|
|
local_container.push_front(c_i);
|
|
}
|
|
|
|
local_container.push_front(begin);
|
|
|
|
for(cri=local_container.begin();cri!=local_container.end();cri++)
|
|
*result++=**cri;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] ApproxError[m];
|
|
delete[] ApproxError;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] SplitPoints[m];
|
|
delete[] SplitPoints;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Err[m];
|
|
delete[] Err;
|
|
|
|
return result;
|
|
};
|
|
|
|
|
|
|
|
// Bounded-E (min-#) global error assessment
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_DP_be_gea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t &approx_n_pt,
|
|
typename DistTraits::FT error_bound,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
typename std::size_t n,m,i;
|
|
FT err_crt;
|
|
InputIterator c_n,c_i,c_m;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_n_pt=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=1;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=2;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
for(n=0,c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
typename std::size_t N=n+1;
|
|
|
|
FT **ApproxError;
|
|
typename std::size_t **SplitPoints;
|
|
FT **Err;
|
|
|
|
ApproxError=new FT*[N];
|
|
assert(ApproxError);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
ApproxError[m]=new FT[N];
|
|
assert(ApproxError[m]);
|
|
}
|
|
|
|
Err=new FT*[N];
|
|
assert(Err);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
Err[m]=new FT[N];
|
|
assert(Err[m]);
|
|
}
|
|
|
|
Err[0][0]=Err[0][1]=Err[1][0]=FT(0);
|
|
|
|
if(N>2)
|
|
for(n=2,c_n=begin,c_n++,c_n++,c_n++;n<N;n++,c_n++)
|
|
ApproxError[n][1]=Err[n][0]=Err[0][n]=error(begin,c_n);
|
|
|
|
if(ApproxError[N-1][1]<=error_bound)
|
|
{
|
|
approx_n_pt=2;
|
|
|
|
c_i=beyond;
|
|
c_i--;
|
|
|
|
*result++=*begin;
|
|
*result=*c_i;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] ApproxError[m];
|
|
delete[] ApproxError;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Err[m];
|
|
delete[] Err;
|
|
|
|
return result;
|
|
}
|
|
|
|
SplitPoints=new typename std::size_t*[N];
|
|
assert(SplitPoints);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
SplitPoints[m]=new typename std::size_t[N];
|
|
assert(SplitPoints[m]);
|
|
}
|
|
|
|
Err[N-1][N-1]=FT(0);
|
|
for(n=1;n<N-1;n++)
|
|
Err[n][n]=Err[n][n-1]=Err[n][n+1]=Err[n-1][n]=Err[n+1][n]=FT(0);
|
|
|
|
for(n=1,c_n=begin,c_n++;n<N-2;n++,c_n++)
|
|
for(m=n+2,c_m=c_n,c_m++,c_m++,c_m++;m<N;m++,c_m++)
|
|
Err[m][n]=Err[n][m]=error(c_n,c_m);
|
|
|
|
for(m=2;m<N;m++)
|
|
{
|
|
ApproxError[m][m]=FT(0);
|
|
SplitPoints[m][m]=m-1;
|
|
|
|
for(n=m+1;n<N;n++)
|
|
{
|
|
ApproxError[n][m]=ApproxError[n-1][m-1];
|
|
SplitPoints[n][m]=n-1;
|
|
for(i=n-2;i>=m;i--)
|
|
{
|
|
err_crt=ApproxError[i][m-1]+Err[i][n];
|
|
if(err_crt<ApproxError[n][m])
|
|
{
|
|
ApproxError[n][m]=err_crt;
|
|
SplitPoints[n][m]=i;
|
|
}
|
|
}
|
|
err_crt=Err[m-1][n];
|
|
if(err_crt<ApproxError[n][m])
|
|
{
|
|
ApproxError[n][m]=err_crt;
|
|
SplitPoints[n][m]=m-1;
|
|
}
|
|
}
|
|
|
|
if(ApproxError[N-1][m]<=error_bound)
|
|
{
|
|
approx_n_pt=m+1;
|
|
break;
|
|
}
|
|
}
|
|
|
|
c_i=beyond;
|
|
c_i--;
|
|
|
|
typename std::list<InputIterator> local_container;
|
|
typename std::list<InputIterator>::iterator cri;
|
|
|
|
local_container.push_front(c_i);
|
|
|
|
n=N-1;
|
|
|
|
typename std::size_t nn=N-1;
|
|
|
|
for(m=approx_n_pt-1;m>=2;m--)
|
|
{
|
|
n=SplitPoints[n][m];
|
|
|
|
for(typename std::size_t ii=nn;ii>n;ii--)
|
|
c_i--;
|
|
|
|
nn=n;
|
|
local_container.push_front(c_i);
|
|
}
|
|
|
|
local_container.push_front(begin);
|
|
|
|
for(cri=local_container.begin();cri!=local_container.end();cri++)
|
|
*result++=**cri;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] ApproxError[m];
|
|
delete[] ApproxError;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] SplitPoints[m];
|
|
delete[] SplitPoints;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Err[m];
|
|
delete[] Err;
|
|
|
|
return result;
|
|
};
|
|
|
|
|
|
//////////////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Recursive Split
|
|
|
|
// Bounded-E local error assessment recursive implementation
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_RSr_be_lea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t &approx_n_pt,
|
|
typename DistTraits::FT error_bound,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
static bool fl;
|
|
FT crt_err;
|
|
InputIterator split_pt,temp_it;
|
|
|
|
if(fl==0)
|
|
{
|
|
if(begin==beyond)
|
|
{
|
|
approx_n_pt=0;
|
|
return result;
|
|
}
|
|
|
|
temp_it=begin;
|
|
temp_it++;
|
|
if(temp_it==beyond)
|
|
{
|
|
approx_n_pt=1;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
temp_it++;
|
|
if(temp_it==beyond)
|
|
{
|
|
approx_n_pt=2;
|
|
*result++=*begin;
|
|
temp_it--;
|
|
*result++=*temp_it;
|
|
return result;
|
|
}
|
|
|
|
fl=1;
|
|
approx_n_pt=1;
|
|
*result++=*begin; // save begin as the first point of the result
|
|
}
|
|
|
|
crt_err=error(begin,beyond,split_pt);
|
|
|
|
if(crt_err>error_bound)
|
|
{
|
|
temp_it=split_pt;
|
|
temp_it++;
|
|
polygonal_approximation_RSr_be_lea(begin,temp_it,approx_n_pt,error_bound,result,error);
|
|
polygonal_approximation_RSr_be_lea(split_pt,beyond,approx_n_pt,error_bound,result,error);
|
|
}
|
|
else
|
|
{
|
|
temp_it=beyond;
|
|
temp_it--;
|
|
*result++=*temp_it;
|
|
approx_n_pt++;
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
template<class InputIterator>
|
|
struct Sp_List
|
|
{ InputIterator begin;
|
|
InputIterator beyond;
|
|
InputIterator split_pt;
|
|
Sp_List *next;
|
|
};
|
|
|
|
template<class InputIterator,class FTp>
|
|
struct WaitList
|
|
{ FTp err;
|
|
Sp_List<InputIterator> *adr_el;
|
|
WaitList *next;
|
|
};
|
|
|
|
|
|
// Bounded-E local error assessment
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_RS_be_lea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t &approx_n_pt,
|
|
typename DistTraits::FT error_bound,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
Sp_List<InputIterator> *l_sp,*sp_cr,*sp_new;
|
|
InputIterator split_pt,c_n;
|
|
FT crt_err;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_n_pt=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=1;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=2;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
l_sp=new Sp_List<InputIterator>;
|
|
assert(l_sp);
|
|
|
|
l_sp->next=NULL;
|
|
l_sp->begin=begin;
|
|
l_sp->beyond=beyond;
|
|
|
|
crt_err=error(begin,beyond,split_pt);
|
|
|
|
l_sp->split_pt=split_pt;
|
|
|
|
sp_cr=l_sp;
|
|
|
|
approx_n_pt=2;
|
|
|
|
while(crt_err>error_bound)
|
|
{
|
|
sp_new=new Sp_List<InputIterator>;
|
|
assert(sp_new);
|
|
|
|
approx_n_pt++;
|
|
|
|
sp_new->next=sp_cr->next;
|
|
sp_cr->next=sp_new;
|
|
|
|
sp_new->beyond=sp_cr->beyond;
|
|
sp_new->begin=sp_cr->split_pt;
|
|
sp_cr->beyond=sp_cr->split_pt;
|
|
sp_cr->beyond++;
|
|
|
|
crt_err=error(sp_cr->begin,sp_cr->beyond,split_pt);
|
|
|
|
while(sp_cr->next && crt_err<=error_bound)
|
|
{
|
|
sp_cr=sp_cr->next;
|
|
crt_err=error(sp_cr->begin,sp_cr->beyond,split_pt);
|
|
}
|
|
|
|
sp_cr->split_pt=split_pt;
|
|
}
|
|
|
|
*result++=*begin;
|
|
|
|
while(l_sp!=NULL)
|
|
{
|
|
sp_cr=l_sp;
|
|
l_sp->beyond--;
|
|
*result++=*l_sp->beyond;
|
|
l_sp=l_sp->next;
|
|
delete sp_cr;
|
|
}
|
|
|
|
return result;
|
|
};
|
|
|
|
|
|
// Bounded-# local error assessment
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_RS_bnp_lea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t n_pt_bound,
|
|
typename DistTraits::FT &approx_error,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
InputIterator c_n;
|
|
typename std::size_t n=0;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_error=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
if(n_pt_bound<2)
|
|
return result;
|
|
|
|
if(n_pt_bound==2)
|
|
{
|
|
approx_error=error(begin,beyond,c_n);
|
|
*result++=*begin;
|
|
c_n=beyond;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
for(c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
if(n_pt_bound>n)
|
|
{
|
|
approx_error=0;
|
|
|
|
for(c_n=begin;c_n!=beyond;c_n++)
|
|
*result++=*c_n;
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
Sp_List<InputIterator> *l_sp,*sp_cr,*sp_new;
|
|
WaitList<InputIterator,FT> *l_wt,*wt_cr,*wt_new;
|
|
InputIterator split_pt;
|
|
|
|
l_sp=new Sp_List<InputIterator>;
|
|
assert(l_sp);
|
|
l_wt=new WaitList<InputIterator,FT>;
|
|
assert(l_wt);
|
|
|
|
l_sp->next=NULL;
|
|
l_sp->begin=begin;
|
|
l_sp->beyond=beyond;
|
|
|
|
l_wt->next=NULL;
|
|
l_wt->adr_el=l_sp;
|
|
l_wt->err=error(begin,beyond,split_pt);
|
|
|
|
l_sp->split_pt=split_pt;
|
|
|
|
typename std::size_t m=2;
|
|
|
|
while(m<n_pt_bound)
|
|
{
|
|
sp_cr=l_wt->adr_el;
|
|
|
|
sp_new=new Sp_List<InputIterator>;
|
|
assert(sp_new);
|
|
|
|
sp_new->next=sp_cr->next;
|
|
sp_cr->next=sp_new;
|
|
|
|
sp_new->beyond=sp_cr->beyond;
|
|
sp_new->begin=sp_cr->split_pt;
|
|
sp_cr->beyond=sp_cr->split_pt;
|
|
sp_cr->beyond++;
|
|
|
|
wt_new=new WaitList<InputIterator,FT>;
|
|
assert(wt_new);
|
|
|
|
wt_new->err=error(sp_cr->begin,sp_cr->beyond,split_pt);
|
|
sp_cr->split_pt=split_pt;
|
|
wt_new->adr_el=sp_cr;
|
|
|
|
wt_cr=l_wt;
|
|
while(wt_cr->next!=NULL && wt_cr->next->err>wt_new->err)
|
|
wt_cr=wt_cr->next;
|
|
wt_new->next=wt_cr->next;
|
|
wt_cr->next=wt_new;
|
|
|
|
wt_new=new WaitList<InputIterator,FT>;
|
|
assert(wt_new);
|
|
|
|
wt_new->err=error(sp_new->begin,sp_new->beyond,split_pt);
|
|
wt_new->adr_el=sp_new;
|
|
sp_new->split_pt=split_pt;
|
|
|
|
wt_cr=l_wt;
|
|
while(wt_cr->next!=NULL && wt_cr->next->err>wt_new->err)
|
|
wt_cr=wt_cr->next;
|
|
wt_new->next=wt_cr->next;
|
|
wt_cr->next=wt_new;
|
|
|
|
wt_cr=l_wt;
|
|
l_wt=l_wt->next;
|
|
delete wt_cr;
|
|
|
|
m++;
|
|
}
|
|
|
|
approx_error=l_wt->err;
|
|
|
|
*result++=*begin;
|
|
|
|
while(l_sp!=NULL)
|
|
{
|
|
sp_cr=l_sp;
|
|
l_sp->beyond--;
|
|
*result++=*l_sp->beyond;
|
|
l_sp=l_sp->next;
|
|
delete sp_cr;
|
|
}
|
|
|
|
while(l_wt!=NULL)
|
|
{
|
|
wt_cr=l_wt;
|
|
l_wt=l_wt->next;
|
|
delete wt_cr;
|
|
}
|
|
|
|
return result;
|
|
};
|
|
|
|
|
|
// Bounded-E global error assessment
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_RS_be_gea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t &approx_n_pt,
|
|
typename DistTraits::FT error_bound,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
Sp_List<InputIterator> *l_sp,*sp_cr,*sp_new;
|
|
WaitList<InputIterator,FT> *l_wt,*wt_cr,*wt_new;
|
|
InputIterator split_pt,c_n;
|
|
FT crt_err;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_n_pt=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=1;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=2;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
l_sp=new Sp_List<InputIterator>;
|
|
assert(l_sp);
|
|
l_wt=new WaitList<InputIterator,FT>;
|
|
assert(l_wt);
|
|
|
|
l_sp->next=NULL;
|
|
l_sp->begin=begin;
|
|
l_sp->beyond=beyond;
|
|
|
|
l_wt->next=NULL;
|
|
l_wt->adr_el=l_sp;
|
|
l_wt->err=error(begin,beyond,split_pt);
|
|
|
|
l_sp->split_pt=split_pt;
|
|
|
|
typename std::size_t m=2;
|
|
|
|
crt_err=l_wt->err;
|
|
|
|
while(crt_err>error_bound)
|
|
{
|
|
sp_cr=l_wt->adr_el;
|
|
|
|
crt_err-=l_wt->err;
|
|
|
|
sp_new=new Sp_List<InputIterator>;
|
|
assert(sp_new);
|
|
|
|
sp_new->next=sp_cr->next;
|
|
sp_cr->next=sp_new;
|
|
|
|
sp_new->beyond=sp_cr->beyond;
|
|
sp_new->begin=sp_cr->split_pt;
|
|
sp_cr->beyond=sp_cr->split_pt;
|
|
sp_cr->beyond++;
|
|
|
|
wt_new=new WaitList<InputIterator,FT>;
|
|
assert(wt_new);
|
|
|
|
wt_new->err=error(sp_cr->begin,sp_cr->beyond,split_pt);
|
|
crt_err+=wt_new->err;
|
|
sp_cr->split_pt=split_pt;
|
|
wt_new->adr_el=sp_cr;
|
|
|
|
wt_cr=l_wt;
|
|
while(wt_cr->next!=NULL && wt_cr->next->err>wt_new->err)
|
|
wt_cr=wt_cr->next;
|
|
wt_new->next=wt_cr->next;
|
|
wt_cr->next=wt_new;
|
|
|
|
wt_new=new WaitList<InputIterator,FT>;
|
|
assert(wt_new);
|
|
|
|
wt_new->err=error(sp_new->begin,sp_new->beyond,split_pt);
|
|
crt_err+=wt_new->err;
|
|
wt_new->adr_el=sp_new;
|
|
sp_new->split_pt=split_pt;
|
|
|
|
wt_cr=l_wt;
|
|
while(wt_cr->next!=NULL && wt_cr->next->err>wt_new->err)
|
|
wt_cr=wt_cr->next;
|
|
wt_new->next=wt_cr->next;
|
|
wt_cr->next=wt_new;
|
|
|
|
wt_cr=l_wt;
|
|
l_wt=l_wt->next;
|
|
delete wt_cr;
|
|
|
|
m++;
|
|
}
|
|
|
|
approx_n_pt=m;
|
|
|
|
*result++=*begin;
|
|
|
|
while(l_sp!=NULL)
|
|
{
|
|
sp_cr=l_sp;
|
|
l_sp->beyond--;
|
|
*result++=*l_sp->beyond;
|
|
l_sp=l_sp->next;
|
|
delete sp_cr;
|
|
}
|
|
|
|
while(l_wt!=NULL)
|
|
{
|
|
wt_cr=l_wt;
|
|
l_wt=l_wt->next;
|
|
delete wt_cr;
|
|
}
|
|
|
|
return result;
|
|
};
|
|
|
|
|
|
// Bounded-# global error assessment
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_RS_bnp_gea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t n_pt_bound,
|
|
typename DistTraits::FT &approx_error,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
InputIterator c_n;
|
|
typename std::size_t n=0;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_error=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
if(n_pt_bound<2)
|
|
return result;
|
|
|
|
if(n_pt_bound==2)
|
|
{
|
|
approx_error=error(begin,beyond,c_n);
|
|
*result++=*begin;
|
|
c_n=beyond;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
for(c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
if(n_pt_bound>n)
|
|
{
|
|
approx_error=0;
|
|
|
|
for(c_n=begin;c_n!=beyond;c_n++)
|
|
*result++=*c_n;
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
Sp_List<InputIterator> *l_sp,*sp_cr,*sp_new;
|
|
WaitList<InputIterator,FT> *l_wt,*wt_cr,*wt_new;
|
|
InputIterator split_pt;
|
|
FT crt_err;
|
|
|
|
l_sp=new Sp_List<InputIterator>;
|
|
assert(l_sp);
|
|
l_wt=new WaitList<InputIterator,FT>;
|
|
assert(l_wt);
|
|
|
|
l_sp->next=NULL;
|
|
l_sp->begin=begin;
|
|
l_sp->beyond=beyond;
|
|
|
|
l_wt->next=NULL;
|
|
l_wt->adr_el=l_sp;
|
|
l_wt->err=error(begin,beyond,split_pt);
|
|
|
|
l_sp->split_pt=split_pt;
|
|
|
|
typename std::size_t m=2;
|
|
|
|
crt_err=l_wt->err;
|
|
|
|
while(m<n_pt_bound)
|
|
{
|
|
sp_cr=l_wt->adr_el;
|
|
|
|
crt_err-=l_wt->err;
|
|
|
|
sp_new=new Sp_List<InputIterator>;
|
|
assert(sp_new);
|
|
|
|
sp_new->next=sp_cr->next;
|
|
sp_cr->next=sp_new;
|
|
|
|
sp_new->beyond=sp_cr->beyond;
|
|
sp_new->begin=sp_cr->split_pt;
|
|
sp_cr->beyond=sp_cr->split_pt;
|
|
sp_cr->beyond++;
|
|
|
|
wt_new=new WaitList<InputIterator,FT>;
|
|
assert(wt_new);
|
|
|
|
wt_new->err=error(sp_cr->begin,sp_cr->beyond,split_pt);
|
|
crt_err+=wt_new->err;
|
|
sp_cr->split_pt=split_pt;
|
|
wt_new->adr_el=sp_cr;
|
|
|
|
wt_cr=l_wt;
|
|
while(wt_cr->next!=NULL && wt_cr->next->err>wt_new->err)
|
|
wt_cr=wt_cr->next;
|
|
wt_new->next=wt_cr->next;
|
|
wt_cr->next=wt_new;
|
|
|
|
wt_new=new WaitList<InputIterator,FT>;
|
|
assert(wt_new);
|
|
|
|
wt_new->err=error(sp_new->begin,sp_new->beyond,split_pt);
|
|
crt_err+=wt_new->err;
|
|
wt_new->adr_el=sp_new;
|
|
sp_new->split_pt=split_pt;
|
|
|
|
wt_cr=l_wt;
|
|
while(wt_cr->next!=NULL && wt_cr->next->err>wt_new->err)
|
|
wt_cr=wt_cr->next;
|
|
wt_new->next=wt_cr->next;
|
|
wt_cr->next=wt_new;
|
|
|
|
wt_cr=l_wt;
|
|
l_wt=l_wt->next;
|
|
delete wt_cr;
|
|
|
|
m++;
|
|
}
|
|
|
|
approx_error=crt_err;
|
|
|
|
*result++=*begin;
|
|
|
|
while(l_sp!=NULL)
|
|
{
|
|
sp_cr=l_sp;
|
|
l_sp->beyond--;
|
|
*result++=*l_sp->beyond;
|
|
l_sp=l_sp->next;
|
|
delete sp_cr;
|
|
}
|
|
|
|
while(l_wt!=NULL)
|
|
{
|
|
wt_cr=l_wt;
|
|
l_wt=l_wt->next;
|
|
delete wt_cr;
|
|
}
|
|
|
|
return result;
|
|
};
|
|
|
|
|
|
///////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Graph Search
|
|
|
|
// Bounded-E local error assessment
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_GS_be_lea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t &approx_n_pt,
|
|
typename DistTraits::FT error_bound,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
typename std::size_t path;
|
|
typename std::size_t n,m,m_pt,i,n_vis,m_crt;
|
|
typename DistTraits::FT err_crt,max,vm,err_max;
|
|
InputIterator c_n,c_i,c_m;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_n_pt=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=1;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=2;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
for(n=0,c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
typename std::size_t N=n+1;
|
|
|
|
typedef struct cel
|
|
{ typename std::size_t pt_crt;
|
|
struct cel *next;
|
|
} CEL;
|
|
|
|
CEL *queue_b,*queue_e,*queue_crt;
|
|
queue_e=queue_b=queue_crt=NULL;
|
|
|
|
CEL *breadth_crt;
|
|
breadth_crt=NULL;
|
|
|
|
m_crt=N+1;
|
|
err_crt=error_bound+1;
|
|
|
|
FT **Graph;
|
|
CEL **visited;
|
|
|
|
visited=new CEL*[N];
|
|
assert(visited);
|
|
for(m=0;m<N;m++)
|
|
visited[m]=NULL;
|
|
|
|
Graph=new FT*[N];
|
|
assert(Graph);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
Graph[m]=new FT[N];
|
|
assert(Graph[m]);
|
|
}
|
|
|
|
Graph[0][0]=Graph[0][1]=Graph[1][0]=FT(0);
|
|
Graph[N-1][N-1]=FT(0);
|
|
|
|
for(n=1;n<N-1;n++)
|
|
Graph[n][n]=Graph[n][n-1]=Graph[n][n+1]=Graph[n-1][n]=Graph[n+1][n]=FT(0);
|
|
|
|
for(n=0,c_n=begin;n<N-2;n++,c_n++)
|
|
for(m=n+2,c_m=c_n,c_m++,c_m++,c_m++;m<N;m++,c_m++)
|
|
Graph[n][m]=Graph[m][n]=error(c_n,c_m);
|
|
|
|
if(Graph[0][N-1]<=error_bound)
|
|
{
|
|
approx_n_pt=2;
|
|
|
|
*result++=*begin;
|
|
c_n=beyond;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
if(Graph[0][N-1]<=error_bound)
|
|
Graph[0][N-1]=Graph[N-1][0]=error_bound+1;
|
|
|
|
queue_b=queue_e=new CEL;
|
|
assert(queue_b);
|
|
queue_b->next=NULL;
|
|
queue_b->pt_crt=0;
|
|
|
|
breadth_crt=new CEL;
|
|
assert(breadth_crt);
|
|
breadth_crt->pt_crt=0;
|
|
breadth_crt->next=NULL;
|
|
|
|
visited[0]=breadth_crt;
|
|
n_vis=1;
|
|
|
|
while(n_vis<N)
|
|
{
|
|
for(i=queue_b->pt_crt+1;i<N;i++)
|
|
{
|
|
if(Graph[queue_b->pt_crt][i]<=error_bound)
|
|
{
|
|
if(i==N-1)
|
|
{
|
|
max=Graph[queue_b->pt_crt][N-1];
|
|
breadth_crt=visited[queue_b->pt_crt];
|
|
m_pt=2;
|
|
while(breadth_crt->next)
|
|
{
|
|
m_pt++;
|
|
vm=Graph[breadth_crt->pt_crt][breadth_crt->next->pt_crt];
|
|
if(vm>max)
|
|
max=vm;
|
|
breadth_crt=breadth_crt->next;
|
|
}
|
|
if(m_pt<m_crt || (m_pt==m_crt && err_max>max))
|
|
{
|
|
m_crt=m_pt;
|
|
err_max=max;
|
|
path=queue_b->pt_crt;
|
|
}
|
|
}
|
|
|
|
if(visited[i]==0)
|
|
{
|
|
breadth_crt=new CEL;
|
|
assert(breadth_crt);
|
|
|
|
breadth_crt->pt_crt=i;
|
|
breadth_crt->next=visited[queue_b->pt_crt];
|
|
visited[i]=breadth_crt;
|
|
n_vis++;
|
|
|
|
queue_e->next=new CEL;
|
|
assert(queue_e->next);
|
|
|
|
queue_e=queue_e->next;
|
|
queue_e->next=NULL;
|
|
queue_e->pt_crt=i;
|
|
}
|
|
}
|
|
}
|
|
|
|
queue_crt=queue_b;
|
|
queue_b=queue_b->next;
|
|
delete queue_crt;
|
|
}
|
|
|
|
while(queue_b)
|
|
{
|
|
if(queue_b->pt_crt!=N-1)
|
|
if(Graph[queue_b->pt_crt][N-1]<=error_bound)
|
|
{
|
|
max=Graph[queue_b->pt_crt][N-1];
|
|
breadth_crt=visited[queue_b->pt_crt];
|
|
m_pt=2;
|
|
while(breadth_crt->next)
|
|
{
|
|
m_pt++;
|
|
vm=Graph[breadth_crt->pt_crt][breadth_crt->next->pt_crt];
|
|
if(vm>max)
|
|
max=vm;
|
|
breadth_crt=breadth_crt->next;
|
|
}
|
|
if(m_pt<m_crt || (m_pt==m_crt && err_max>max))
|
|
{
|
|
m_crt=m_pt;
|
|
err_max=max;
|
|
path=queue_b->pt_crt;
|
|
}
|
|
}
|
|
queue_crt=queue_b;
|
|
queue_b=queue_b->next;
|
|
delete queue_crt;
|
|
}
|
|
|
|
approx_n_pt=m_crt;
|
|
|
|
delete visited[N-1];
|
|
visited[N-1]=NULL;
|
|
|
|
breadth_crt=visited[path]->next;
|
|
|
|
delete visited[path];
|
|
visited[path]=NULL;
|
|
|
|
while(breadth_crt)
|
|
{
|
|
path=breadth_crt->pt_crt;
|
|
breadth_crt=breadth_crt->next;
|
|
|
|
delete visited[path];
|
|
visited[path]=NULL;
|
|
}
|
|
|
|
for(c_n=begin,i=0;c_n!=beyond;c_n++,i++)
|
|
if(visited[i]==NULL)
|
|
*result++=*c_n;
|
|
|
|
for(i=0;i<N;i++)
|
|
if(visited[i])
|
|
delete visited[i];
|
|
|
|
delete[] visited;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Graph[m];
|
|
delete[] Graph;
|
|
|
|
return result;
|
|
};
|
|
|
|
|
|
|
|
template<class FT>
|
|
void QuickSort(FT* p,typename std::size_t n) //quick sort
|
|
{
|
|
FT tmp,pivot;
|
|
|
|
if(n<2)
|
|
return;
|
|
|
|
if(n==2)
|
|
{
|
|
if(p[0]>p[1])
|
|
{
|
|
tmp=p[0];
|
|
p[0]=p[1];
|
|
p[1]=tmp;
|
|
}
|
|
return;
|
|
}
|
|
|
|
pivot=p[n/2];
|
|
std::size_t i=0,j=n-1;
|
|
|
|
while(i<j)
|
|
{
|
|
while(p[i]<pivot)
|
|
i++;
|
|
while(p[j]>pivot)
|
|
j--;
|
|
if(i<j)
|
|
{
|
|
tmp=p[i];
|
|
p[i]=p[j];
|
|
p[j]=tmp;
|
|
if(p[i]==p[j])
|
|
{
|
|
i++;
|
|
j--;
|
|
}
|
|
}
|
|
}
|
|
|
|
if(i==j)
|
|
{ if(p[i]>pivot)
|
|
i=i-1;
|
|
}
|
|
else if(j<i)
|
|
i=j;
|
|
QuickSort<FT>(p,i+1);
|
|
j=n-i-1;
|
|
QuickSort<FT>(&p[i+1],j);
|
|
}
|
|
|
|
|
|
// Bounded-# local error assessment
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_GS_bnp_lea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t n_pt_bound,
|
|
typename DistTraits::FT &approx_error,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
typename std::size_t path;
|
|
typename std::size_t n,m,m_pt,i,n_vis,m_crt;
|
|
typename DistTraits::FT err_crt,max,vm,err_max;
|
|
InputIterator c_n,c_i,c_m;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_error=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
if(n_pt_bound<2)
|
|
return result;
|
|
|
|
n=0;
|
|
for(c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
if(n_pt_bound>n)
|
|
{
|
|
approx_error=0;
|
|
|
|
for(c_n=begin;c_n!=beyond;c_n++)
|
|
*result++=*c_n;
|
|
|
|
return result;
|
|
}
|
|
|
|
typename std::size_t N=n+1;
|
|
|
|
typedef struct cel
|
|
{ typename std::size_t pt_crt;
|
|
struct cel *next;
|
|
} CEL;
|
|
|
|
CEL *queue_b,*queue_e,*queue_crt;
|
|
queue_e=queue_b=queue_crt=NULL;
|
|
|
|
CEL *breadth_crt;
|
|
breadth_crt=NULL;
|
|
|
|
FT **Graph;
|
|
CEL **visited;
|
|
|
|
visited=new CEL*[N];
|
|
assert(visited);
|
|
for(m=0;m<N;m++)
|
|
visited[m]=NULL;
|
|
|
|
Graph=new FT*[N];
|
|
assert(Graph);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
Graph[m]=new FT[N];
|
|
assert(Graph[m]);
|
|
}
|
|
|
|
Graph[0][0]=Graph[0][1]=Graph[1][0]=FT(0);
|
|
Graph[N-1][N-1]=FT(0);
|
|
|
|
for(n=1;n<N-1;n++)
|
|
Graph[n][n]=Graph[n][n-1]=Graph[n][n+1]=Graph[n-1][n]=Graph[n+1][n]=FT(0);
|
|
|
|
typename std::size_t N_err=(N-1)*(N-2)/2;
|
|
FT *List_err;
|
|
List_err=new FT[N_err];
|
|
assert(List_err);
|
|
|
|
typename std::size_t idx=0;
|
|
|
|
for(n=0,c_n=begin;n<N-2;n++,c_n++)
|
|
for(m=n+2,c_m=c_n,c_m++,c_m++,c_m++;m<N;m++,c_m++)
|
|
Graph[n][m]=Graph[m][n]=List_err[idx++]=error(c_n,c_m);
|
|
|
|
if(n_pt_bound==2)
|
|
{
|
|
approx_error=Graph[0][N-1];
|
|
*result++=*begin;
|
|
c_n=beyond;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
FT cr_err_bound;
|
|
bool fl_circ=false;
|
|
|
|
typename std::size_t N_crt,N_mf;
|
|
double delta_N,nc;
|
|
|
|
if(Graph[0][N-1]<=0)
|
|
fl_circ=true;
|
|
|
|
QuickSort<FT>(List_err,N_err);
|
|
|
|
delta_N=N_err/2.;
|
|
nc=0;
|
|
m=n_pt_bound+1;
|
|
|
|
while(delta_N>0.4)
|
|
{
|
|
if(m>n_pt_bound)
|
|
nc+=delta_N;
|
|
else
|
|
nc-=delta_N;
|
|
|
|
N_crt=(int)(nc+0.5);
|
|
if(N_crt<0)
|
|
N_crt=0;
|
|
if(N_crt>=N_err)
|
|
N_crt=N_err;
|
|
|
|
cr_err_bound=List_err[N_crt];
|
|
m_crt=N+1;
|
|
err_crt=cr_err_bound+1;
|
|
|
|
if(fl_circ)
|
|
Graph[0][N-1]=Graph[N-1][0]=cr_err_bound+1;
|
|
|
|
for(m=0;m<N;m++)
|
|
visited[m]=NULL;
|
|
|
|
queue_b=queue_e=new CEL;
|
|
assert(queue_b);
|
|
|
|
queue_b->next=NULL;
|
|
queue_b->pt_crt=0;
|
|
|
|
breadth_crt=new CEL;
|
|
assert(breadth_crt);
|
|
|
|
breadth_crt->pt_crt=0;
|
|
breadth_crt->next=NULL;
|
|
|
|
visited[0]=breadth_crt;
|
|
n_vis=1;
|
|
|
|
while(n_vis<N)
|
|
{
|
|
for(i=(queue_b->pt_crt+1);i<N;i++)
|
|
{
|
|
if(Graph[queue_b->pt_crt][i]<=cr_err_bound)
|
|
{
|
|
if(i==N-1)
|
|
{
|
|
max=Graph[queue_b->pt_crt][N-1];
|
|
breadth_crt=visited[queue_b->pt_crt];
|
|
m_pt=2;
|
|
while(breadth_crt->next)
|
|
{
|
|
m_pt++;
|
|
vm=Graph[breadth_crt->pt_crt][breadth_crt->next->pt_crt];
|
|
if(vm>max)
|
|
max=vm;
|
|
breadth_crt=breadth_crt->next;
|
|
}
|
|
if(m_pt<m_crt || (m_pt==m_crt && err_max>max))
|
|
{
|
|
m_crt=m_pt;
|
|
err_max=max;
|
|
path=queue_b->pt_crt;
|
|
}
|
|
}
|
|
|
|
if(visited[i]==0)
|
|
{
|
|
breadth_crt=new CEL;
|
|
assert(breadth_crt);
|
|
breadth_crt->pt_crt=i;
|
|
breadth_crt->next=visited[queue_b->pt_crt];
|
|
visited[i]=breadth_crt;
|
|
n_vis++;
|
|
queue_e->next=new CEL;
|
|
assert(queue_e->next);
|
|
queue_e=queue_e->next;
|
|
queue_e->next=NULL;
|
|
queue_e->pt_crt=i;
|
|
}
|
|
}
|
|
}
|
|
|
|
queue_crt=queue_b;
|
|
queue_b=queue_b->next;
|
|
delete queue_crt;
|
|
}
|
|
|
|
while(queue_b)
|
|
{
|
|
if(queue_b->pt_crt!=N-1)
|
|
if(Graph[queue_b->pt_crt][N-1]<=cr_err_bound)
|
|
{
|
|
max=Graph[queue_b->pt_crt][N-1];
|
|
breadth_crt=visited[queue_b->pt_crt];
|
|
m_pt=2;
|
|
while(breadth_crt->next)
|
|
{
|
|
m_pt++;
|
|
vm=Graph[breadth_crt->pt_crt][breadth_crt->next->pt_crt];
|
|
if(vm>max)
|
|
max=vm;
|
|
breadth_crt=breadth_crt->next;
|
|
}
|
|
if(m_pt<m_crt || (m_pt==m_crt && err_max>max))
|
|
{
|
|
m_crt=m_pt;
|
|
err_max=max;
|
|
path=queue_b->pt_crt;
|
|
}
|
|
}
|
|
queue_crt=queue_b;
|
|
queue_b=queue_b->next;
|
|
delete queue_crt;
|
|
}
|
|
|
|
m=m_crt;
|
|
|
|
delete visited[N-1];
|
|
visited[N-1]=NULL;
|
|
|
|
breadth_crt=visited[path]->next;
|
|
|
|
delete visited[path];
|
|
visited[path]=NULL;
|
|
|
|
while(breadth_crt)
|
|
{
|
|
path=breadth_crt->pt_crt;
|
|
breadth_crt=breadth_crt->next;
|
|
|
|
delete visited[path];
|
|
visited[path]=NULL;
|
|
}
|
|
|
|
for(i=0;i<N;i++)
|
|
if(visited[i])
|
|
delete visited[i];
|
|
|
|
approx_error=cr_err_bound;
|
|
|
|
delta_N/=2.;
|
|
|
|
if(m<=n_pt_bound)
|
|
N_mf=N_crt;
|
|
}
|
|
|
|
if(m!=n_pt_bound)
|
|
N_crt=N_mf;
|
|
|
|
cr_err_bound=List_err[N_crt];
|
|
m_crt=N+1;
|
|
err_crt=cr_err_bound+1;
|
|
|
|
if(fl_circ)
|
|
Graph[0][N-1]=Graph[N-1][0]=cr_err_bound+1;
|
|
|
|
for(m=0;m<N;m++)
|
|
visited[m]=NULL;
|
|
|
|
queue_b=queue_e=new CEL;
|
|
assert(queue_b);
|
|
|
|
queue_b->next=NULL;
|
|
queue_b->pt_crt=0;
|
|
|
|
breadth_crt=new CEL;
|
|
assert(breadth_crt);
|
|
|
|
breadth_crt->pt_crt=0;
|
|
breadth_crt->next=NULL;
|
|
|
|
visited[0]=breadth_crt;
|
|
n_vis=1;
|
|
|
|
while(n_vis<N)
|
|
{
|
|
for(i=queue_b->pt_crt+1;i<N;i++)
|
|
{
|
|
if(Graph[queue_b->pt_crt][i]<=cr_err_bound)
|
|
{
|
|
if(i==N-1)
|
|
{
|
|
max=Graph[queue_b->pt_crt][N-1];
|
|
breadth_crt=visited[queue_b->pt_crt];
|
|
m_pt=2;
|
|
while(breadth_crt->next)
|
|
{
|
|
m_pt++;
|
|
vm=Graph[breadth_crt->pt_crt][breadth_crt->next->pt_crt];
|
|
if(vm>max)
|
|
max=vm;
|
|
breadth_crt=breadth_crt->next;
|
|
}
|
|
if(m_pt<m_crt || (m_pt==m_crt && err_max>max))
|
|
{
|
|
m_crt=m_pt;
|
|
err_max=max;
|
|
path=queue_b->pt_crt;
|
|
}
|
|
}
|
|
|
|
if(visited[i]==0)
|
|
{
|
|
breadth_crt=new CEL;
|
|
assert(breadth_crt);
|
|
|
|
breadth_crt->pt_crt=i;
|
|
breadth_crt->next=visited[queue_b->pt_crt];
|
|
visited[i]=breadth_crt;
|
|
n_vis++;
|
|
|
|
queue_e->next=new CEL;
|
|
assert(queue_e->next);
|
|
|
|
queue_e=queue_e->next;
|
|
queue_e->next=NULL;
|
|
queue_e->pt_crt=i;
|
|
}
|
|
}
|
|
}
|
|
|
|
queue_crt=queue_b;
|
|
queue_b=queue_b->next;
|
|
delete queue_crt;
|
|
}
|
|
|
|
while(queue_b)
|
|
{
|
|
if(queue_b->pt_crt!=N-1)
|
|
if(Graph[queue_b->pt_crt][N-1]<=cr_err_bound)
|
|
{
|
|
max=Graph[queue_b->pt_crt][N-1];
|
|
breadth_crt=visited[queue_b->pt_crt];
|
|
m_pt=2;
|
|
while(breadth_crt->next)
|
|
{
|
|
m_pt++;
|
|
vm=Graph[breadth_crt->pt_crt][breadth_crt->next->pt_crt];
|
|
|
|
if(vm>max)
|
|
max=vm;
|
|
breadth_crt=breadth_crt->next;
|
|
}
|
|
if(m_pt<m_crt || (m_pt==m_crt && err_max>max))
|
|
{
|
|
m_crt=m_pt;
|
|
err_max=max;
|
|
path=queue_b->pt_crt;
|
|
}
|
|
}
|
|
queue_crt=queue_b;
|
|
queue_b=queue_b->next;
|
|
delete queue_crt;
|
|
}
|
|
|
|
m=m_crt;
|
|
n_pt_bound=m;
|
|
|
|
delete visited[N-1];
|
|
visited[N-1]=NULL;
|
|
|
|
breadth_crt=visited[path]->next;
|
|
|
|
delete visited[path];
|
|
visited[path]=NULL;
|
|
|
|
while(breadth_crt)
|
|
{
|
|
path=breadth_crt->pt_crt;
|
|
breadth_crt=breadth_crt->next;
|
|
|
|
delete visited[path];
|
|
visited[path]=NULL;
|
|
}
|
|
|
|
for(c_n=begin,i=0;c_n!=beyond;c_n++,i++)
|
|
if(visited[i]==NULL)
|
|
*result++=*c_n;
|
|
|
|
for(i=0;i<N;i++)
|
|
if(visited[i])
|
|
delete visited[i];
|
|
|
|
approx_error=cr_err_bound;
|
|
|
|
delete[] visited;
|
|
|
|
delete[] List_err;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Graph[m];
|
|
delete[] Graph;
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
// Bounded-E global error assessment
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_GS_be_gea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t &approx_n_pt,
|
|
typename DistTraits::FT error_bound,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
typename std::size_t n,m,i,j;
|
|
InputIterator c_n,c_i,c_m;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_n_pt=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=1;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=2;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
for(n=0,c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
typename std::size_t N=n+1;
|
|
|
|
FT **Graph;
|
|
|
|
Graph=new FT*[N];
|
|
assert(Graph);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
Graph[m]=new FT[N];
|
|
assert(Graph[m]);
|
|
}
|
|
|
|
Graph[0][0]=Graph[0][1]=Graph[1][0]=FT(0);
|
|
Graph[N-1][N-1]=FT(0);
|
|
|
|
for(n=1;n<N-1;n++)
|
|
Graph[n][n]=Graph[n][n-1]=Graph[n][n+1]=Graph[n-1][n]=Graph[n+1][n]=FT(0);
|
|
|
|
for(n=0,c_n=begin;n<N-2;n++,c_n++)
|
|
for(m=n+2,c_m=c_n,c_m++,c_m++,c_m++;m<N;m++,c_m++)
|
|
Graph[n][m]=Graph[m][n]=error(c_n,c_m);
|
|
|
|
if(Graph[0][N-1]<=error_bound)
|
|
{
|
|
approx_n_pt=2;
|
|
|
|
*result++=*begin;
|
|
c_n=beyond;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
if(Graph[0][N-1]<=error_bound)
|
|
Graph[0][N-1]=Graph[N-1][0]=error_bound+1;
|
|
|
|
bool fl_find=false;
|
|
FT err_crt;
|
|
FT* ApproxError;
|
|
typename std::size_t **SplitPoints;
|
|
|
|
SplitPoints=new typename std::size_t*[N+1];
|
|
assert(SplitPoints);
|
|
|
|
ApproxError=new FT[N+1];
|
|
assert(ApproxError);
|
|
|
|
m=2;
|
|
|
|
for(i=0;i<=N-1;i++)
|
|
ApproxError[i]=Graph[i][0];
|
|
|
|
while(!fl_find)
|
|
{
|
|
m++;
|
|
SplitPoints[m]=new typename std::size_t[N];
|
|
assert(SplitPoints[m]);
|
|
|
|
ApproxError[N-1]=ApproxError[N-2];
|
|
SplitPoints[m][N-1]=N-2;
|
|
if(ApproxError[N-1]<=error_bound)
|
|
fl_find=true;
|
|
else
|
|
for(j=N-3;j>=m-2;j--)
|
|
{
|
|
if(ApproxError[j]<=error_bound && Graph[j][N-1]<=error_bound)
|
|
{
|
|
err_crt=ApproxError[j]+Graph[j][N-1];
|
|
|
|
if(err_crt<ApproxError[N-1])
|
|
{
|
|
ApproxError[N-1]=err_crt;
|
|
SplitPoints[m][N-1]=j;
|
|
if(ApproxError[N-1]<=error_bound)
|
|
{
|
|
fl_find=true;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if(!fl_find)
|
|
for(i=N-2;i>=m-1;i--)
|
|
{
|
|
ApproxError[i]=ApproxError[i-1];
|
|
SplitPoints[m][i]=i-1;
|
|
for(j=i-2;j>=m-2;j--)
|
|
{
|
|
if(ApproxError[j]<=error_bound && Graph[j][i]<=error_bound)
|
|
{
|
|
err_crt=ApproxError[j]+Graph[j][i];
|
|
if(err_crt<ApproxError[i])
|
|
{
|
|
ApproxError[i]=err_crt;
|
|
SplitPoints[m][i]=j;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
approx_n_pt=m;
|
|
|
|
c_i=beyond;
|
|
c_i--;
|
|
|
|
typename std::list<InputIterator> local_container;
|
|
typename std::list<InputIterator>::iterator cri;
|
|
|
|
local_container.push_front(c_i);
|
|
|
|
n=N-1;
|
|
|
|
typename std::size_t nn=N-1;
|
|
|
|
for(m=approx_n_pt;m>2;m--)
|
|
{
|
|
n=SplitPoints[m][n];
|
|
|
|
for(typename std::size_t ii=nn;ii>n;ii--)
|
|
c_i--;
|
|
|
|
nn=n;
|
|
local_container.push_front(c_i);
|
|
}
|
|
|
|
local_container.push_front(begin);
|
|
|
|
for(cri=local_container.begin();cri!=local_container.end();cri++)
|
|
*result++=**cri;
|
|
|
|
delete[] ApproxError;
|
|
|
|
for(m=3;m<=approx_n_pt;m++)
|
|
delete[] SplitPoints[m];
|
|
delete[] SplitPoints;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Graph[m];
|
|
delete[] Graph;
|
|
|
|
return result;
|
|
};
|
|
|
|
|
|
// Bounded-# global error assessment
|
|
|
|
template<class InputIterator,class OutputIterator, class DistTraits>
|
|
OutputIterator polygonal_approximation_GS_bnp_gea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t n_pt_bound,
|
|
typename DistTraits::FT &approx_error,
|
|
OutputIterator result,
|
|
DistTraits error )
|
|
{
|
|
typedef typename DistTraits::FT FT;
|
|
|
|
typename std::size_t n,m,i,j;
|
|
InputIterator c_n,c_i,c_m;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_error=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
if(n_pt_bound<2)
|
|
return result;
|
|
|
|
n=0;
|
|
for(c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
if(n_pt_bound>n)
|
|
{
|
|
approx_error=0;
|
|
|
|
for(c_n=begin;c_n!=beyond;c_n++)
|
|
*result++=*c_n;
|
|
|
|
return result;
|
|
}
|
|
|
|
typename std::size_t N=n+1;
|
|
|
|
FT **Graph;
|
|
|
|
Graph=new FT*[N];
|
|
assert(Graph);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
Graph[m]=new FT[N];
|
|
assert(Graph[m]);
|
|
}
|
|
|
|
Graph[0][0]=Graph[0][1]=Graph[1][0]=FT(0);
|
|
Graph[N-1][N-1]=FT(0);
|
|
|
|
for(n=1;n<N-1;n++)
|
|
Graph[n][n]=Graph[n][n-1]=Graph[n][n+1]=Graph[n-1][n]=Graph[n+1][n]=FT(0);
|
|
|
|
for(n=0,c_n=begin;n<N-2;n++,c_n++)
|
|
for(m=n+2,c_m=c_n,c_m++,c_m++,c_m++;m<N;m++,c_m++)
|
|
Graph[n][m]=Graph[m][n]=error(c_n,c_m);
|
|
|
|
if(n_pt_bound==2)
|
|
{
|
|
approx_error=Graph[0][N-1];
|
|
*result++=*begin;
|
|
c_n=beyond;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
typename std::size_t mm;
|
|
|
|
FT err_crt;
|
|
FT* ApproxError;
|
|
typename std::size_t **SplitPoints;
|
|
|
|
SplitPoints=new typename std::size_t*[N+1];
|
|
assert(SplitPoints);
|
|
|
|
ApproxError=new FT[N+1];
|
|
assert(ApproxError);
|
|
|
|
for(i=0;i<=N-1;i++)
|
|
ApproxError[i]=Graph[i][0];
|
|
|
|
FT epsi_crt;
|
|
|
|
mm=2;
|
|
|
|
while(mm<n_pt_bound)
|
|
{
|
|
mm++;
|
|
SplitPoints[mm]=new typename std::size_t[N];
|
|
assert(SplitPoints[mm]);
|
|
|
|
epsi_crt=ApproxError[N-1]=ApproxError[N-2];
|
|
SplitPoints[mm][N-1]=N-2;
|
|
|
|
for(j=N-3;j>=mm-2;j--)
|
|
{
|
|
if(ApproxError[j]<=epsi_crt && Graph[j][N-1]<=epsi_crt)
|
|
{
|
|
err_crt=ApproxError[j]+Graph[j][N-1];
|
|
|
|
if(err_crt<ApproxError[N-1])
|
|
{
|
|
epsi_crt=ApproxError[N-1]=err_crt;
|
|
SplitPoints[mm][N-1]=j;
|
|
}
|
|
}
|
|
}
|
|
|
|
if(mm<n_pt_bound)
|
|
for(i=N-2;i>=mm-1;i--)
|
|
{
|
|
epsi_crt=ApproxError[i]=ApproxError[i-1];
|
|
SplitPoints[mm][i]=i-1;
|
|
for(j=i-2;j>=mm-2;j--)
|
|
{
|
|
if(ApproxError[j]<=epsi_crt && Graph[j][i]<=epsi_crt)
|
|
{
|
|
err_crt=ApproxError[j]+Graph[j][i];
|
|
if(err_crt<ApproxError[i])
|
|
{
|
|
epsi_crt=ApproxError[i]=err_crt;
|
|
SplitPoints[mm][i]=j;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
approx_error=ApproxError[N-1];
|
|
|
|
c_i=beyond;
|
|
c_i--;
|
|
|
|
typename std::list<InputIterator> local_container;
|
|
typename std::list<InputIterator>::iterator cri;
|
|
|
|
local_container.push_front(c_i);
|
|
|
|
n=N-1;
|
|
|
|
typename std::size_t nn=N-1;
|
|
|
|
for(m=n_pt_bound;m>2;m--)
|
|
{
|
|
n=SplitPoints[m][n];
|
|
|
|
for(typename std::size_t ii=nn;ii>n;ii--)
|
|
c_i--;
|
|
|
|
nn=n;
|
|
local_container.push_front(c_i);
|
|
}
|
|
|
|
local_container.push_front(begin);
|
|
|
|
for(cri=local_container.begin();cri!=local_container.end();cri++)
|
|
*result++=**cri;
|
|
|
|
delete[] ApproxError;
|
|
|
|
for(m=3;m<=mm;m++)
|
|
delete[] SplitPoints[m];
|
|
delete[] SplitPoints;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Graph[m];
|
|
delete[] Graph;
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
///////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Optimal algorithms for Squared Euclidean distance
|
|
|
|
template<class InputIterator, class FT>
|
|
bool slope_sign(InputIterator p,InputIterator q,InputIterator begin,InputIterator end)
|
|
{
|
|
FT sHomogX,sHomogY,deltaXpq,deltaYpq,slopeVal;
|
|
|
|
sHomogX=begin->y()-end->y();
|
|
sHomogY=end->x()-begin->x();
|
|
|
|
if(sHomogX==FT(0) && sHomogY==FT(0))
|
|
sHomogX=FT(1);
|
|
|
|
deltaXpq=q->x()-p->x();
|
|
deltaYpq=q->y()-p->y();
|
|
slopeVal=sHomogX*deltaXpq+sHomogY*deltaYpq;
|
|
|
|
return (slopeVal>=FT(0) );
|
|
}
|
|
|
|
|
|
// Bounded-E local error assessment; Based on the online convex hull (Toussaint)
|
|
|
|
template<class InputIterator,class OutputIterator, class FT>
|
|
OutputIterator polygonal_approximation_CHGS_be_lea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t &approx_n_pt,
|
|
FT error_bound,
|
|
OutputIterator result )
|
|
{
|
|
typename std::size_t path;
|
|
|
|
typename std::size_t n,m,m_pt,i,n_vis,m_crt;
|
|
FT err_crt,max,vm,err_max;
|
|
InputIterator c_n,c_i,c_m;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_n_pt=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=1;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=2;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
for(n=0,c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
typename std::size_t N=n+1;
|
|
|
|
|
|
// build the double stack
|
|
|
|
typename std::size_t top; //Top end of elt.
|
|
typename std::size_t bot; //Bottom end of elt.
|
|
InputIterator* elt; //Double ended queue storing a convex hull. dim=TWICE_HULL_MAX
|
|
|
|
typename std::size_t HULL_MAX;
|
|
typename std::size_t TWICE_HULL_MAX;
|
|
typename std::size_t THRICE_HULL_MAX;
|
|
|
|
bool topflag,botflag;
|
|
|
|
HULL_MAX=N;
|
|
TWICE_HULL_MAX=2*HULL_MAX;
|
|
THRICE_HULL_MAX=3*HULL_MAX;
|
|
|
|
elt=new InputIterator[TWICE_HULL_MAX];
|
|
assert(elt);
|
|
|
|
for (typename std::size_t j=0;j<TWICE_HULL_MAX;j++)
|
|
elt[j]=NULL;
|
|
|
|
typedef struct cel
|
|
{ typename std::size_t pt_crt;
|
|
struct cel *next;
|
|
} CEL;
|
|
|
|
CEL *queue_b,*queue_e,*queue_crt;
|
|
queue_e=queue_b=queue_crt=NULL;
|
|
|
|
CEL *breadth_crt;
|
|
breadth_crt=NULL;
|
|
|
|
m_crt=N+1;
|
|
err_crt=error_bound+1;
|
|
|
|
FT **Graph;
|
|
CEL **visited;
|
|
|
|
visited=new CEL*[N];
|
|
assert(visited);
|
|
|
|
for(m=0;m<N;m++)
|
|
visited[m]=NULL;
|
|
|
|
Graph=new FT*[N];
|
|
assert(Graph);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
Graph[m]=new FT[N];
|
|
assert(Graph[m]);
|
|
}
|
|
|
|
Graph[0][0]=Graph[0][1]=Graph[1][0]=FT(0);
|
|
|
|
Graph[N-1][N-1]=FT(0);
|
|
for(n=1;n<N-1;n++)
|
|
Graph[n][n]=Graph[n][n-1]=Graph[n][n+1]=Graph[n-1][n]=Graph[n+1][n]=FT(0);
|
|
|
|
for(n=0,c_n=begin;n<N-2;n++,c_n++)
|
|
{
|
|
top=0;
|
|
bot=0;
|
|
c_i=c_n;
|
|
c_i++;
|
|
|
|
elt[HULL_MAX]=c_n;
|
|
top=HULL_MAX+1;
|
|
elt[top]=c_i;
|
|
bot=HULL_MAX-1;
|
|
elt[bot]=c_i;
|
|
|
|
for(m=n+2,c_m=c_n,c_m++,c_m++;c_m!=beyond;m++,c_m++)
|
|
{
|
|
topflag=((elt[top]->x()-c_m->x())*(elt[top-1]->y()-c_m->y())) >= ((elt[top-1]->x()-c_m->x())*(elt[top]->y()-c_m->y())); // left_of
|
|
botflag=((elt[bot+1]->x()-c_m->x())*(elt[bot]->y()-c_m->y())) >= ((elt[bot]->x()-c_m->x())*(elt[bot+1]->y()-c_m->y()));
|
|
|
|
if(topflag || botflag) //If the new point is outside the hull.
|
|
{
|
|
while(topflag && (top>HULL_MAX)) //Pop points until convexity is ensured.
|
|
{
|
|
top--;
|
|
topflag=((elt[top]->x()-c_m->x())*(elt[top-1]->y()-c_m->y())) >= ((elt[top-1]->x()-c_m->x())*(elt[top]->y()-c_m->y()));
|
|
}
|
|
|
|
while(botflag && (bot<HULL_MAX))
|
|
{
|
|
bot++;
|
|
botflag=((elt[bot+1]->x()-c_m->x())*(elt[bot]->y()-c_m->y())) >= ((elt[bot]->x()-c_m->x())*(elt[bot+1]->y()-c_m->y()));
|
|
}
|
|
|
|
top++; //Then push the new point on the top and bottom
|
|
bot--; //of the queue.
|
|
elt[top]=c_m;
|
|
elt[bot]=c_m;
|
|
}
|
|
|
|
// find the furthest point
|
|
|
|
typename std::size_t mid,lo,m1,brk,m2,hi;
|
|
bool sbase, sbrk;
|
|
FT d1,d2,dx,dy,d,tr_y;
|
|
|
|
if((top-bot)>6) //If there are > 6 points on the hull
|
|
{ //(otherwise we will just look at them all.
|
|
lo=bot;
|
|
hi=top-1;
|
|
sbase = slope_sign<InputIterator,FT>(elt[hi],elt[lo],c_n,c_m); //The sign of the base edge.
|
|
|
|
do
|
|
{
|
|
brk = (lo + hi) / 2; //Binary search for an edge with opposite sign.
|
|
sbrk = slope_sign<InputIterator,FT>(elt[brk], elt[brk+1], c_n,c_m);
|
|
|
|
if(sbase == sbrk )
|
|
if(sbase == (slope_sign<InputIterator,FT>(elt[lo], elt[brk+1], c_n,c_m)) )
|
|
lo = brk + 1;
|
|
else
|
|
hi = brk;
|
|
}
|
|
while(sbase==sbrk);
|
|
|
|
m1=brk; //Now, the sign changes between the base edge and
|
|
while(lo<m1) // brk+1. Binary search for the extreme point.
|
|
{
|
|
mid=(lo+ m1)/2;
|
|
|
|
if(sbase==(slope_sign<InputIterator,FT>(elt[mid],elt[mid+1],c_n,c_m)) )
|
|
lo=mid+1;
|
|
else
|
|
m1=mid;
|
|
}
|
|
|
|
m2=brk; //The sign also chagnes between brk and the base.
|
|
while(m2<hi) //Binary search again.
|
|
{
|
|
mid=(m2+hi)/2;
|
|
|
|
if(sbase==slope_sign<InputIterator,FT>(elt[mid],elt[mid+1],c_n,c_m) )
|
|
hi = mid;
|
|
else
|
|
m2=mid+1;
|
|
}
|
|
|
|
dx=c_m->x()-c_n->x();
|
|
dy=c_m->y()-c_n->y();
|
|
d=dx*dx+dy*dy;
|
|
if(d==0)
|
|
{
|
|
d1=(elt[lo]->x()-c_n->x())*(elt[lo]->x()-c_n->x())+(elt[lo]->y()-c_n->y())*(elt[lo]->y()-c_n->y());
|
|
d2=(elt[hi]->x()-c_n->x())*(elt[hi]->x()-c_n->x())+(elt[hi]->y()-c_n->y())*(elt[hi]->y()-c_n->y());
|
|
}
|
|
else
|
|
{
|
|
//d1=squared_distance(begin,end,elt[lo]);
|
|
tr_y=(elt[lo]->y()-c_n->y())*dx-(elt[lo]->x()-c_n->x())*dy;
|
|
d1=tr_y*tr_y/d;
|
|
|
|
//d2=squared_distance(begin,end,elt[hi]);
|
|
tr_y=(elt[hi]->y()-c_n->y())*dx-(elt[hi]->x()-c_n->x())*dy;
|
|
d2=tr_y*tr_y/d;
|
|
}
|
|
|
|
if(d1>d2)
|
|
d2=d1;
|
|
}
|
|
else //Few points in hull--search by brute force.
|
|
{
|
|
d2=FT(0);
|
|
for(mid=bot;mid<top;mid++)
|
|
{
|
|
//d1=squared_distance(begin,end,elt[mid]);
|
|
dx=c_m->x()-c_n->x();
|
|
dy=c_m->y()-c_n->y();
|
|
d=dx*dx+dy*dy;
|
|
if(d==0)
|
|
{
|
|
d1=(elt[mid]->x()-c_n->x())*(elt[mid]->x()-c_n->x())+(elt[mid]->y()-c_n->y())*(elt[mid]->y()-c_n->y());
|
|
}
|
|
else
|
|
{
|
|
tr_y=(elt[mid]->y()-c_n->y())*dx-(elt[mid]->x()-c_n->x())*dy;
|
|
d1=tr_y*tr_y/d;
|
|
}
|
|
|
|
if(d1>d2)
|
|
d2=d1;
|
|
}
|
|
}
|
|
Graph[n][m]=Graph[m][n]=d2;
|
|
}
|
|
}
|
|
|
|
delete[] elt;
|
|
elt=NULL;
|
|
|
|
if(Graph[0][N-1]<=error_bound)
|
|
{
|
|
approx_n_pt=2;
|
|
|
|
*result++=*begin;
|
|
c_n=beyond;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
if(Graph[0][N-1]<=error_bound)
|
|
Graph[0][N-1]=Graph[N-1][0]=error_bound+1;
|
|
|
|
queue_b=queue_e=new CEL;
|
|
assert(queue_b);
|
|
|
|
queue_b->next=NULL;
|
|
queue_b->pt_crt=0;
|
|
|
|
breadth_crt=new CEL;
|
|
assert(breadth_crt);
|
|
breadth_crt->pt_crt=0;
|
|
breadth_crt->next=NULL;
|
|
|
|
visited[0]=breadth_crt;
|
|
n_vis=1;
|
|
|
|
while(n_vis<N)
|
|
{
|
|
for(i=queue_b->pt_crt+1;i<N;i++)
|
|
{
|
|
if(Graph[queue_b->pt_crt][i]<=error_bound)
|
|
{
|
|
if(i==N-1)
|
|
{
|
|
max=Graph[queue_b->pt_crt][N-1];
|
|
breadth_crt=visited[queue_b->pt_crt];
|
|
m_pt=2;
|
|
while(breadth_crt->next)
|
|
{
|
|
m_pt++;
|
|
vm=Graph[breadth_crt->pt_crt][breadth_crt->next->pt_crt];
|
|
if(vm>max)
|
|
max=vm;
|
|
breadth_crt=breadth_crt->next;
|
|
}
|
|
if(m_pt<m_crt || (m_pt==m_crt && err_max>max))
|
|
{
|
|
m_crt=m_pt;
|
|
err_max=max;
|
|
path=queue_b->pt_crt;
|
|
}
|
|
}
|
|
|
|
if(visited[i]==0)
|
|
{
|
|
breadth_crt=new CEL;
|
|
assert(breadth_crt);
|
|
|
|
breadth_crt->pt_crt=i;
|
|
breadth_crt->next=visited[queue_b->pt_crt];
|
|
visited[i]=breadth_crt;
|
|
n_vis++;
|
|
|
|
queue_e->next=new CEL;
|
|
assert(queue_e->next);
|
|
|
|
queue_e=queue_e->next;
|
|
queue_e->next=NULL;
|
|
queue_e->pt_crt=i;
|
|
}
|
|
}
|
|
}
|
|
|
|
queue_crt=queue_b;
|
|
queue_b=queue_b->next;
|
|
delete queue_crt;
|
|
}
|
|
|
|
while(queue_b)
|
|
{
|
|
if(queue_b->pt_crt!=N-1)
|
|
if(Graph[queue_b->pt_crt][N-1]<=error_bound)
|
|
{
|
|
max=Graph[queue_b->pt_crt][N-1];
|
|
breadth_crt=visited[queue_b->pt_crt];
|
|
m_pt=2;
|
|
while(breadth_crt->next)
|
|
{
|
|
m_pt++;
|
|
vm=Graph[breadth_crt->pt_crt][breadth_crt->next->pt_crt];
|
|
if(vm>max)
|
|
max=vm;
|
|
breadth_crt=breadth_crt->next;
|
|
}
|
|
if(m_pt<m_crt || (m_pt==m_crt && err_max>max))
|
|
{
|
|
m_crt=m_pt;
|
|
err_max=max;
|
|
path=queue_b->pt_crt;
|
|
}
|
|
}
|
|
queue_crt=queue_b;
|
|
queue_b=queue_b->next;
|
|
delete queue_crt;
|
|
}
|
|
|
|
approx_n_pt=m_crt;
|
|
|
|
delete visited[N-1];
|
|
visited[N-1]=NULL;
|
|
|
|
breadth_crt=visited[path]->next;
|
|
|
|
delete visited[path];
|
|
visited[path]=NULL;
|
|
|
|
while(breadth_crt)
|
|
{
|
|
path=breadth_crt->pt_crt;
|
|
breadth_crt=breadth_crt->next;
|
|
|
|
delete visited[path];
|
|
visited[path]=NULL;
|
|
}
|
|
|
|
for(c_n=begin,i=0;c_n!=beyond;c_n++,i++)
|
|
if(visited[i]==NULL)
|
|
*result++=*c_n;
|
|
|
|
for(i=0;i<N;i++)
|
|
if(visited[i])
|
|
delete visited[i];
|
|
|
|
delete[] visited;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Graph[m];
|
|
delete[] Graph;
|
|
|
|
return result;
|
|
};
|
|
|
|
|
|
// Bounded-# local error assessment; Based on the online convex hull (Toussaint)
|
|
|
|
template<class InputIterator,class OutputIterator, class FT>
|
|
OutputIterator polygonal_approximation_CHGS_bnp_lea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t n_pt_bound,
|
|
FT &approx_error,
|
|
OutputIterator result )
|
|
{
|
|
typename std::size_t path;
|
|
typename std::size_t n,m,i;
|
|
typename std::size_t m_pt,n_vis,m_crt;
|
|
FT err_crt,max,vm,err_max;
|
|
InputIterator c_n,c_i,c_m;
|
|
FT cr_err_bound;
|
|
bool fl_circ=false;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_error=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
if(n_pt_bound<2)
|
|
return result;
|
|
|
|
n=0;
|
|
for(c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
if(n_pt_bound>n)
|
|
{
|
|
approx_error=0;
|
|
|
|
for(c_n=begin;c_n!=beyond;c_n++)
|
|
*result++=*c_n;
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
typename std::size_t N=n+1;
|
|
typename std::size_t N_err=(N-1)*(N-2)/2;
|
|
typename std::size_t N_crt,N_mf;
|
|
double delta_N,nc;
|
|
|
|
|
|
// build the double stack
|
|
|
|
typename std::size_t top; //Top end of elt.
|
|
typename std::size_t bot; //Bottom end of elt.
|
|
InputIterator* elt; //Double ended queue storing a convex hull. dim=TWICE_HULL_MAX
|
|
|
|
typename std::size_t HULL_MAX;
|
|
typename std::size_t TWICE_HULL_MAX;
|
|
typename std::size_t THRICE_HULL_MAX;
|
|
|
|
bool topflag,botflag;
|
|
|
|
HULL_MAX=N;
|
|
TWICE_HULL_MAX=2*HULL_MAX;
|
|
THRICE_HULL_MAX=3*HULL_MAX;
|
|
|
|
elt=new InputIterator[TWICE_HULL_MAX];
|
|
assert(elt);
|
|
|
|
for (typename std::size_t j=0;j<TWICE_HULL_MAX;j++)
|
|
elt[j]=NULL;
|
|
|
|
typedef struct cel
|
|
{ typename std::size_t pt_crt;
|
|
struct cel *next;
|
|
} CEL;
|
|
|
|
CEL *queue_b,*queue_e,*queue_crt;
|
|
queue_e=queue_b=queue_crt=NULL;
|
|
|
|
CEL *breadth_crt;
|
|
breadth_crt=NULL;
|
|
|
|
CEL **visited;
|
|
|
|
visited=new CEL*[N];
|
|
assert(visited);
|
|
|
|
for(m=0;m<N;m++)
|
|
visited[m]=NULL;
|
|
|
|
FT *List_err;
|
|
FT **Graph;
|
|
|
|
List_err=new FT[N_err];
|
|
assert(List_err);
|
|
|
|
Graph=new FT*[N];
|
|
assert(Graph);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
Graph[m]=new FT[N];
|
|
assert(Graph[m]);
|
|
}
|
|
|
|
Graph[0][0]=Graph[0][1]=Graph[1][0]=FT(0);
|
|
Graph[N-1][N-1]=FT(0);
|
|
for(n=1;n<N-1;n++)
|
|
Graph[n][n]=Graph[n][n-1]=Graph[n][n+1]=Graph[n-1][n]=Graph[n+1][n]=FT(0);
|
|
|
|
typename std::size_t index1=0;
|
|
|
|
for(n=0,c_n=begin;n<N-2;n++,c_n++)
|
|
{
|
|
top=0;
|
|
bot=0;
|
|
c_i=c_n;
|
|
c_i++;
|
|
|
|
elt[HULL_MAX]=c_n;
|
|
top=HULL_MAX+1;
|
|
elt[top]=c_i;
|
|
bot=HULL_MAX-1;
|
|
elt[bot]=c_i;
|
|
|
|
for(m=n+2,c_m=c_n,c_m++,c_m++;c_m!=beyond;m++,c_m++)
|
|
{
|
|
topflag=((elt[top]->x()-c_m->x())*(elt[top-1]->y()-c_m->y())) >= ((elt[top-1]->x()-c_m->x())*(elt[top]->y()-c_m->y())); // left_of
|
|
botflag=((elt[bot+1]->x()-c_m->x())*(elt[bot]->y()-c_m->y())) >= ((elt[bot]->x()-c_m->x())*(elt[bot+1]->y()-c_m->y()));
|
|
|
|
if(topflag || botflag) //If the new point is outside the hull.
|
|
{
|
|
while(topflag && (top>HULL_MAX)) //Pop points until convexity is ensured.
|
|
{
|
|
top--;
|
|
topflag=((elt[top]->x()-c_m->x())*(elt[top-1]->y()-c_m->y())) >= ((elt[top-1]->x()-c_m->x())*(elt[top]->y()-c_m->y()));
|
|
}
|
|
|
|
while(botflag && (bot<HULL_MAX))
|
|
{
|
|
bot++;
|
|
botflag=((elt[bot+1]->x()-c_m->x())*(elt[bot]->y()-c_m->y())) >= ((elt[bot]->x()-c_m->x())*(elt[bot+1]->y()-c_m->y()));
|
|
}
|
|
|
|
top++; //Then push the new point on the top and bottom
|
|
bot--; //of the queue.
|
|
elt[top]=c_m;
|
|
elt[bot]=c_m;
|
|
}
|
|
|
|
// find the furthest point
|
|
|
|
typename std::size_t mid,lo,m1,brk,m2,hi;
|
|
bool sbase, sbrk;
|
|
FT d1,d2,dx,dy,d,tr_y;
|
|
|
|
if((top-bot)>6) //If there are > 6 points on the hull
|
|
{ //(otherwise we will just look at them all.
|
|
lo=bot;
|
|
hi=top-1;
|
|
sbase = slope_sign<InputIterator,FT>(elt[hi],elt[lo],c_n,c_m); //The sign of the base edge.
|
|
|
|
do
|
|
{
|
|
brk = (lo + hi) / 2; //Binary search for an edge with opposite sign.
|
|
sbrk = slope_sign<InputIterator,FT>(elt[brk], elt[brk+1], c_n,c_m);
|
|
|
|
if(sbase == sbrk )
|
|
if(sbase == (slope_sign<InputIterator,FT>(elt[lo], elt[brk+1], c_n,c_m)) )
|
|
lo = brk + 1;
|
|
else
|
|
hi = brk;
|
|
}
|
|
while(sbase==sbrk);
|
|
|
|
m1=brk; //Now, the sign changes between the base edge and
|
|
while(lo<m1) // brk+1. Binary search for the extreme point.
|
|
{
|
|
mid=(lo+ m1)/2;
|
|
|
|
if(sbase==(slope_sign<InputIterator,FT>(elt[mid],elt[mid+1],c_n,c_m)) )
|
|
lo=mid+1;
|
|
else
|
|
m1=mid;
|
|
}
|
|
|
|
m2=brk; //The sign also chagnes between brk and the base.
|
|
while(m2<hi) //Binary search again.
|
|
{
|
|
mid=(m2+hi)/2;
|
|
|
|
if(sbase==slope_sign<InputIterator,FT>(elt[mid],elt[mid+1],c_n,c_m) )
|
|
hi = mid;
|
|
else
|
|
m2=mid+1;
|
|
}
|
|
|
|
//d1=squared_distance(begin,end,elt[lo]);
|
|
dx=c_m->x()-c_n->x();
|
|
dy=c_m->y()-c_n->y();
|
|
d=dx*dx+dy*dy;
|
|
if(d==0)
|
|
{
|
|
d1=(elt[lo]->x()-c_n->x())*(elt[lo]->x()-c_n->x())+(elt[lo]->y()-c_n->y())*(elt[lo]->y()-c_n->y());
|
|
d2=(elt[hi]->x()-c_n->x())*(elt[hi]->x()-c_n->x())+(elt[hi]->y()-c_n->y())*(elt[hi]->y()-c_n->y());
|
|
}
|
|
else
|
|
{
|
|
tr_y=(elt[lo]->y()-c_n->y())*dx-(elt[lo]->x()-c_n->x())*dy;
|
|
d1=tr_y*tr_y/d;
|
|
|
|
//d2=squared_distance(begin,end,elt[hi]);
|
|
tr_y=(elt[hi]->y()-c_n->y())*dx-(elt[hi]->x()-c_n->x())*dy;
|
|
d2=tr_y*tr_y/d;
|
|
}
|
|
|
|
if(d1>d2)
|
|
d2=d1;
|
|
}
|
|
else //Few points in hull--search by brute force.
|
|
{
|
|
d2=FT(0);
|
|
for(mid=bot;mid<top;mid++)
|
|
{
|
|
//d1=squared_distance(begin,end,elt[mid]);
|
|
dx=c_m->x()-c_n->x();
|
|
dy=c_m->y()-c_n->y();
|
|
d=dx*dx+dy*dy;
|
|
if(d==0)
|
|
d1=(elt[mid]->x()-c_n->x())*(elt[mid]->x()-c_n->x())+(elt[mid]->y()-c_n->y())*(elt[mid]->y()-c_n->y());
|
|
else
|
|
tr_y=(elt[mid]->y()-c_n->y())*dx-(elt[mid]->x()-c_n->x())*dy;
|
|
d1=tr_y*tr_y/d;
|
|
|
|
if(d1>d2)
|
|
d2=d1;
|
|
}
|
|
}
|
|
|
|
List_err[index1]=d2;
|
|
index1++;
|
|
|
|
Graph[n][m]=Graph[m][n]=d2;
|
|
}
|
|
}
|
|
|
|
delete[] elt;
|
|
elt=NULL;
|
|
|
|
if(n_pt_bound==2)
|
|
{
|
|
approx_error=Graph[0][N-1];
|
|
*result++=*begin;
|
|
c_n=beyond;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
if(Graph[0][N-1]<=0)
|
|
fl_circ=true;
|
|
|
|
QuickSort<FT>(List_err,N_err);
|
|
|
|
delta_N=N_err/2.;
|
|
nc=0;
|
|
m=n_pt_bound+1;
|
|
|
|
while(delta_N>0.4)
|
|
{
|
|
if(m>n_pt_bound)
|
|
nc+=delta_N;
|
|
else
|
|
nc-=delta_N;
|
|
|
|
N_crt=(int)(nc+0.5);
|
|
if(N_crt<0)
|
|
N_crt=0;
|
|
if(N_crt>=N_err)
|
|
N_crt=N_err;
|
|
|
|
cr_err_bound=List_err[N_crt];
|
|
m_crt=N+1;
|
|
err_crt=cr_err_bound+1;
|
|
|
|
if(fl_circ)
|
|
Graph[0][N-1]=Graph[N-1][0]=cr_err_bound+1;
|
|
|
|
for(m=0;m<N;m++)
|
|
visited[m]=NULL;
|
|
|
|
queue_b=queue_e=new CEL;
|
|
assert(queue_b);
|
|
queue_b->next=NULL;
|
|
queue_b->pt_crt=0;
|
|
|
|
breadth_crt=new CEL;
|
|
assert(breadth_crt);
|
|
|
|
breadth_crt->pt_crt=0;
|
|
breadth_crt->next=NULL;
|
|
|
|
visited[0]=breadth_crt;
|
|
n_vis=1;
|
|
|
|
while(n_vis<N)
|
|
{
|
|
for(i=(queue_b->pt_crt+1);i<N;i++)
|
|
{
|
|
if(Graph[queue_b->pt_crt][i]<=cr_err_bound)
|
|
{
|
|
if(i==N-1)
|
|
{
|
|
max=Graph[queue_b->pt_crt][N-1];
|
|
breadth_crt=visited[queue_b->pt_crt];
|
|
m_pt=2;
|
|
while(breadth_crt->next)
|
|
{
|
|
m_pt++;
|
|
vm=Graph[breadth_crt->pt_crt][breadth_crt->next->pt_crt];
|
|
if(vm>max)
|
|
max=vm;
|
|
breadth_crt=breadth_crt->next;
|
|
}
|
|
if(m_pt<m_crt || (m_pt==m_crt && err_max>max))
|
|
{
|
|
m_crt=m_pt;
|
|
err_max=max;
|
|
path=queue_b->pt_crt;
|
|
}
|
|
}
|
|
|
|
if(visited[i]==0)
|
|
{
|
|
breadth_crt=new CEL;
|
|
assert(breadth_crt);
|
|
|
|
breadth_crt->pt_crt=i;
|
|
breadth_crt->next=visited[queue_b->pt_crt];
|
|
visited[i]=breadth_crt;
|
|
n_vis++;
|
|
|
|
queue_e->next=new CEL;
|
|
assert(queue_e->next);
|
|
|
|
queue_e=queue_e->next;
|
|
queue_e->next=NULL;
|
|
queue_e->pt_crt=i;
|
|
}
|
|
}
|
|
}
|
|
|
|
queue_crt=queue_b;
|
|
queue_b=queue_b->next;
|
|
delete queue_crt;
|
|
}
|
|
|
|
while(queue_b)
|
|
{
|
|
if(queue_b->pt_crt!=N-1)
|
|
if(Graph[queue_b->pt_crt][N-1]<=cr_err_bound)
|
|
{
|
|
max=Graph[queue_b->pt_crt][N-1];
|
|
breadth_crt=visited[queue_b->pt_crt];
|
|
m_pt=2;
|
|
while(breadth_crt->next)
|
|
{
|
|
m_pt++;
|
|
vm=Graph[breadth_crt->pt_crt][breadth_crt->next->pt_crt];
|
|
if(vm>max)
|
|
max=vm;
|
|
breadth_crt=breadth_crt->next;
|
|
}
|
|
if(m_pt<m_crt || (m_pt==m_crt && err_max>max))
|
|
{
|
|
m_crt=m_pt;
|
|
err_max=max;
|
|
path=queue_b->pt_crt;
|
|
}
|
|
}
|
|
queue_crt=queue_b;
|
|
queue_b=queue_b->next;
|
|
delete queue_crt;
|
|
}
|
|
|
|
m=m_crt;
|
|
|
|
delete visited[N-1];
|
|
visited[N-1]=NULL;
|
|
|
|
breadth_crt=visited[path]->next;
|
|
|
|
delete visited[path];
|
|
visited[path]=NULL;
|
|
|
|
while(breadth_crt)
|
|
{
|
|
path=breadth_crt->pt_crt;
|
|
breadth_crt=breadth_crt->next;
|
|
|
|
delete visited[path];
|
|
visited[path]=NULL;
|
|
}
|
|
|
|
for(i=0;i<N;i++)
|
|
if(visited[i])
|
|
delete visited[i];
|
|
|
|
approx_error=cr_err_bound;
|
|
|
|
delta_N/=2.;
|
|
|
|
if(m<=n_pt_bound)
|
|
N_mf=N_crt;
|
|
}
|
|
|
|
if(m!=n_pt_bound)
|
|
N_crt=N_mf;
|
|
|
|
cr_err_bound=List_err[N_crt];
|
|
m_crt=N+1;
|
|
err_crt=cr_err_bound+1;
|
|
|
|
if(fl_circ)
|
|
Graph[0][N-1]=Graph[N-1][0]=cr_err_bound+1;
|
|
|
|
for(m=0;m<N;m++)
|
|
visited[m]=NULL;
|
|
|
|
queue_b=queue_e=new CEL;
|
|
assert(queue_b);
|
|
|
|
queue_b->next=NULL;
|
|
queue_b->pt_crt=0;
|
|
|
|
breadth_crt=new CEL;
|
|
assert(breadth_crt);
|
|
breadth_crt->pt_crt=0;
|
|
breadth_crt->next=NULL;
|
|
|
|
visited[0]=breadth_crt;
|
|
n_vis=1;
|
|
|
|
while(n_vis<N)
|
|
{
|
|
for(i=queue_b->pt_crt+1;i<N;i++)
|
|
{
|
|
if(Graph[queue_b->pt_crt][i]<=cr_err_bound)
|
|
{
|
|
if(i==N-1)
|
|
{
|
|
max=Graph[queue_b->pt_crt][N-1];
|
|
breadth_crt=visited[queue_b->pt_crt];
|
|
m_pt=2;
|
|
while(breadth_crt->next)
|
|
{
|
|
m_pt++;
|
|
vm=Graph[breadth_crt->pt_crt][breadth_crt->next->pt_crt];
|
|
if(vm>max)
|
|
max=vm;
|
|
breadth_crt=breadth_crt->next;
|
|
}
|
|
if(m_pt<m_crt || (m_pt==m_crt && err_max>max))
|
|
{
|
|
m_crt=m_pt;
|
|
err_max=max;
|
|
path=queue_b->pt_crt;
|
|
}
|
|
}
|
|
|
|
if(visited[i]==0)
|
|
{
|
|
breadth_crt=new CEL;
|
|
assert(breadth_crt);
|
|
|
|
breadth_crt->pt_crt=i;
|
|
breadth_crt->next=visited[queue_b->pt_crt];
|
|
visited[i]=breadth_crt;
|
|
n_vis++;
|
|
|
|
queue_e->next=new CEL;
|
|
assert(queue_e->next);
|
|
|
|
queue_e=queue_e->next;
|
|
queue_e->next=NULL;
|
|
queue_e->pt_crt=i;
|
|
}
|
|
}
|
|
}
|
|
|
|
queue_crt=queue_b;
|
|
queue_b=queue_b->next;
|
|
delete queue_crt;
|
|
}
|
|
|
|
while(queue_b)
|
|
{
|
|
if(queue_b->pt_crt!=N-1)
|
|
if(Graph[queue_b->pt_crt][N-1]<=cr_err_bound)
|
|
{
|
|
max=Graph[queue_b->pt_crt][N-1];
|
|
breadth_crt=visited[queue_b->pt_crt];
|
|
m_pt=2;
|
|
while(breadth_crt->next)
|
|
{
|
|
m_pt++;
|
|
vm=Graph[breadth_crt->pt_crt][breadth_crt->next->pt_crt];
|
|
|
|
if(vm>max)
|
|
max=vm;
|
|
breadth_crt=breadth_crt->next;
|
|
}
|
|
if(m_pt<m_crt || (m_pt==m_crt && err_max>max))
|
|
{
|
|
m_crt=m_pt;
|
|
err_max=max;
|
|
path=queue_b->pt_crt;
|
|
}
|
|
}
|
|
queue_crt=queue_b;
|
|
queue_b=queue_b->next;
|
|
delete queue_crt;
|
|
}
|
|
|
|
m=m_crt;
|
|
n_pt_bound=m;
|
|
|
|
delete visited[N-1];
|
|
visited[N-1]=NULL;
|
|
|
|
breadth_crt=visited[path]->next;
|
|
|
|
delete visited[path];
|
|
visited[path]=NULL;
|
|
|
|
while(breadth_crt)
|
|
{
|
|
path=breadth_crt->pt_crt;
|
|
breadth_crt=breadth_crt->next;
|
|
|
|
delete visited[path];
|
|
visited[path]=NULL;
|
|
}
|
|
|
|
for(c_n=begin,i=0;c_n!=beyond;c_n++,i++)
|
|
if(visited[i]==NULL)
|
|
*result++=*c_n;
|
|
|
|
approx_error=cr_err_bound;
|
|
|
|
for(i=0;i<N;i++)
|
|
if(visited[i])
|
|
delete visited[i];
|
|
|
|
delete[] visited;
|
|
|
|
delete[] List_err;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Graph[m];
|
|
delete[] Graph;
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
// Bounded-E global error assessment; Based on the incremental technique (Perez and Vidal)
|
|
|
|
template<class InputIterator,class OutputIterator, class FT>
|
|
OutputIterator polygonal_approximation_ITGS_be_gea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t &approx_n_pt,
|
|
FT error_bound,
|
|
OutputIterator result )
|
|
{
|
|
typename std::size_t n,m,i,j;
|
|
InputIterator c_n,c_i,c_m;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_n_pt=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=1;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_n_pt=2;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
for(n=0,c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
typename std::size_t N=n+1;
|
|
|
|
FT x_sum,y_sum,x2_sum,y2_sum,xy_sum,xx;
|
|
FT nr_p;
|
|
FT ssvdErr,ssedErr,xk,yk,yValY,dir_coef,dif;
|
|
|
|
FT **Graph;
|
|
|
|
Graph=new FT*[N];
|
|
assert(Graph);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
Graph[m]=new FT[N];
|
|
assert(Graph[m]);
|
|
}
|
|
|
|
Graph[0][0]=Graph[0][1]=Graph[1][0]=FT(0);
|
|
|
|
Graph[N-1][N-1]=FT(0);
|
|
for(n=1;n<N-1;n++)
|
|
Graph[n][n]=Graph[n][n-1]=Graph[n][n+1]=Graph[n-1][n]=Graph[n+1][n]=FT(0);
|
|
|
|
for(n=0,c_n=begin;n<N-2;n++,c_n++)
|
|
{
|
|
x_sum=0;
|
|
y_sum=0;
|
|
x2_sum=0;
|
|
y2_sum=0;
|
|
xy_sum=0;
|
|
nr_p=FT(0);
|
|
|
|
for(m=n+2,c_m=c_n,c_m++;m<N;m++)
|
|
{
|
|
nr_p+=FT(1);
|
|
xk=c_m->x();
|
|
yk=c_m->y();
|
|
|
|
c_m++;
|
|
|
|
x_sum += xk;
|
|
y_sum += yk;
|
|
x2_sum += (xk * xk);
|
|
y2_sum += (yk * yk);
|
|
xy_sum += (xk * yk);
|
|
|
|
//Compute the direction coefficient and the y-value at the y-axis of approx. curve.
|
|
|
|
dif=(c_m->x()-c_n->x());
|
|
|
|
if(dif!=FT(0))
|
|
{
|
|
dir_coef=(c_m->y()-c_n->y())/dif;
|
|
yValY=c_n->y()-dir_coef*c_n->x();
|
|
|
|
ssvdErr=yValY*yValY*nr_p+FT(2)*yValY*dir_coef*x_sum-FT(2)*yValY*y_sum+dir_coef*dir_coef*x2_sum+y2_sum-FT(2)*dir_coef*xy_sum;
|
|
ssedErr=(FT(1)/(dir_coef*dir_coef+FT(1)))*ssvdErr;
|
|
}
|
|
else
|
|
{
|
|
if((c_m->y()-c_n->y())==FT(0)) // Closed curve case (*begin==*end)
|
|
{
|
|
xx=c_n->x();
|
|
yValY=c_n->y();
|
|
|
|
ssedErr=xx*xx*nr_p+yValY*yValY*nr_p-FT(2)*yValY*y_sum+y2_sum+x2_sum-FT(2)*xx*x_sum;
|
|
}
|
|
else
|
|
{
|
|
xx=c_n->x();
|
|
ssedErr=x2_sum-FT(2)*xx*x_sum+nr_p*xx*xx;
|
|
}
|
|
}
|
|
|
|
if(ssedErr<FT(0))
|
|
ssedErr=-ssedErr;
|
|
|
|
Graph[n][m]=Graph[m][n]=ssedErr;
|
|
}
|
|
}
|
|
|
|
if(Graph[0][N-1]<=error_bound)
|
|
{
|
|
approx_n_pt=2;
|
|
|
|
*result++=*begin;
|
|
c_n=beyond;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
bool fl_find=false;
|
|
FT err_crt;
|
|
FT* ApproxError;
|
|
typename std::size_t **SplitPoints;
|
|
|
|
SplitPoints=new typename std::size_t*[N+1];
|
|
assert(SplitPoints);
|
|
|
|
ApproxError=new FT[N+1];
|
|
assert(ApproxError);
|
|
|
|
m=2;
|
|
|
|
for(i=0;i<=N-1;i++)
|
|
ApproxError[i]=Graph[i][0];
|
|
|
|
while(!fl_find)
|
|
{
|
|
m++;
|
|
SplitPoints[m]=new typename std::size_t[N];
|
|
assert(SplitPoints[m]);
|
|
|
|
ApproxError[N-1]=ApproxError[N-2];
|
|
SplitPoints[m][N-1]=N-2;
|
|
if(ApproxError[N-1]<=error_bound)
|
|
fl_find=true;
|
|
else
|
|
for(j=N-3;j>=m-2;j--)
|
|
{
|
|
if(ApproxError[j]<=error_bound && Graph[j][N-1]<=error_bound)
|
|
{
|
|
err_crt=ApproxError[j]+Graph[j][N-1];
|
|
|
|
if(err_crt<ApproxError[N-1])
|
|
{
|
|
ApproxError[N-1]=err_crt;
|
|
SplitPoints[m][N-1]=j;
|
|
if(ApproxError[N-1]<=error_bound)
|
|
{
|
|
fl_find=true;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if(!fl_find)
|
|
for(i=N-2;i>=m-1;i--)
|
|
{
|
|
ApproxError[i]=ApproxError[i-1];
|
|
SplitPoints[m][i]=i-1;
|
|
for(j=i-2;j>=m-2;j--)
|
|
{
|
|
if(ApproxError[j]<=error_bound && Graph[j][i]<=error_bound)
|
|
{
|
|
err_crt=ApproxError[j]+Graph[j][i];
|
|
if(err_crt<ApproxError[i])
|
|
{
|
|
ApproxError[i]=err_crt;
|
|
SplitPoints[m][i]=j;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
approx_n_pt=m;
|
|
|
|
c_i=beyond;
|
|
c_i--;
|
|
|
|
typename std::list<InputIterator> local_container;
|
|
typename std::list<InputIterator>::iterator cri;
|
|
|
|
local_container.push_front(c_i);
|
|
|
|
n=N-1;
|
|
|
|
typename std::size_t nn=N-1;
|
|
|
|
for(m=approx_n_pt;m>2;m--)
|
|
{
|
|
n=SplitPoints[m][n];
|
|
|
|
for(typename std::size_t ii=nn;ii>n;ii--)
|
|
c_i--;
|
|
|
|
nn=n;
|
|
local_container.push_front(c_i);
|
|
}
|
|
|
|
local_container.push_front(begin);
|
|
|
|
for(cri=local_container.begin();cri!=local_container.end();cri++)
|
|
*result++=**cri;
|
|
|
|
delete[] ApproxError;
|
|
|
|
for(m=3;m<=approx_n_pt;m++)
|
|
delete[] SplitPoints[m];
|
|
delete[] SplitPoints;
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Graph[m];
|
|
delete[] Graph;
|
|
|
|
|
|
return result;
|
|
};
|
|
|
|
|
|
// Bounded-# global error assessment; Based on the incremental technique (Perez and Vidal)
|
|
|
|
template<class InputIterator,class OutputIterator, class FT>
|
|
OutputIterator polygonal_approximation_ITGS_bnp_gea( InputIterator begin,
|
|
InputIterator beyond,
|
|
typename std::size_t n_pt_bound,
|
|
FT &approx_error,
|
|
OutputIterator result )
|
|
{
|
|
typename std::size_t n,m,i,j;
|
|
InputIterator c_n,c_i,c_m;
|
|
|
|
if(begin==beyond)
|
|
{
|
|
approx_error=0;
|
|
return result;
|
|
}
|
|
|
|
c_n=begin;
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
return result;
|
|
}
|
|
|
|
c_n++;
|
|
if(c_n==beyond)
|
|
{
|
|
approx_error=0;
|
|
*result++=*begin;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
if(n_pt_bound<2)
|
|
return result;
|
|
|
|
n=0;
|
|
for(c_n=begin,c_n++;c_n!=beyond;c_n++)
|
|
n++;
|
|
|
|
if(n_pt_bound>n)
|
|
{
|
|
approx_error=0;
|
|
|
|
for(c_n=begin;c_n!=beyond;c_n++)
|
|
*result++=*c_n;
|
|
|
|
return result;
|
|
}
|
|
|
|
typename std::size_t N=n+1;
|
|
|
|
FT x_sum,y_sum,x2_sum,y2_sum,xy_sum,xx;
|
|
FT nr_p;
|
|
FT ssvdErr,ssedErr,xk,yk,yValY,dir_coef,dif;
|
|
|
|
FT **Graph;
|
|
|
|
Graph=new FT*[N];
|
|
assert(Graph);
|
|
for(m=0;m<N;m++)
|
|
{
|
|
Graph[m]=new FT[N];
|
|
assert(Graph[m]);
|
|
}
|
|
|
|
Graph[0][0]=Graph[0][1]=Graph[1][0]=FT(0);
|
|
|
|
Graph[N-1][N-1]=FT(0);
|
|
for(n=1;n<N-1;n++)
|
|
Graph[n][n]=Graph[n][n-1]=Graph[n][n+1]=Graph[n-1][n]=Graph[n+1][n]=FT(0);
|
|
|
|
for(n=0,c_n=begin;n<N-2;n++,c_n++)
|
|
{
|
|
x_sum=0;
|
|
y_sum=0;
|
|
x2_sum=0;
|
|
y2_sum=0;
|
|
xy_sum=0;
|
|
nr_p=FT(0);
|
|
|
|
for(m=n+2,c_m=c_n,c_m++;m<N;m++)
|
|
{
|
|
nr_p+=FT(1);
|
|
xk=c_m->x();
|
|
yk=c_m->y();
|
|
|
|
c_m++;
|
|
|
|
x_sum += xk;
|
|
y_sum += yk;
|
|
x2_sum += (xk * xk);
|
|
y2_sum += (yk * yk);
|
|
xy_sum += (xk * yk);
|
|
|
|
//Compute the direction coefficient and the y-value at the y-axis of approx. curve.
|
|
|
|
dif=(c_m->x()-c_n->x());
|
|
|
|
if(dif!=FT(0))
|
|
{
|
|
dir_coef=(c_m->y()-c_n->y())/dif;
|
|
yValY=c_n->y()-dir_coef*c_n->x();
|
|
|
|
ssvdErr=yValY*yValY*nr_p+FT(2)*yValY*dir_coef*x_sum-FT(2)*yValY*y_sum+dir_coef*dir_coef*x2_sum+y2_sum-FT(2)*dir_coef*xy_sum;
|
|
ssedErr=(FT(1)/(dir_coef*dir_coef+FT(1)))*ssvdErr;
|
|
}
|
|
else
|
|
{
|
|
if((c_m->y()-c_n->y())==FT(0)) // Closed curve case (*begin==*end)
|
|
{
|
|
xx=c_n->x();
|
|
yValY=c_n->y();
|
|
|
|
ssedErr=xx*xx*nr_p+yValY*yValY*nr_p-FT(2)*yValY*y_sum+y2_sum+x2_sum-FT(2)*xx*x_sum;
|
|
}
|
|
else
|
|
{
|
|
xx=c_n->x();
|
|
ssedErr=x2_sum-FT(2)*xx*x_sum+nr_p*xx*xx;
|
|
}
|
|
}
|
|
|
|
if(ssedErr<FT(0))
|
|
ssedErr=-ssedErr;
|
|
|
|
Graph[n][m]=Graph[m][n]=ssedErr;
|
|
}
|
|
}
|
|
|
|
if(n_pt_bound==2)
|
|
{
|
|
approx_error=Graph[0][N-1];
|
|
*result++=*begin;
|
|
c_n=beyond;
|
|
c_n--;
|
|
*result++=*c_n;
|
|
return result;
|
|
}
|
|
|
|
|
|
typename std::size_t mm;
|
|
|
|
FT err_crt;
|
|
FT* ApproxError;
|
|
typename std::size_t **SplitPoints;
|
|
|
|
SplitPoints=new typename std::size_t*[N+1];
|
|
assert(SplitPoints);
|
|
|
|
ApproxError=new FT[N+1];
|
|
assert(ApproxError);
|
|
|
|
for(i=0;i<=N-1;i++)
|
|
ApproxError[i]=Graph[i][0];
|
|
|
|
FT epsi_crt;
|
|
|
|
mm=2;
|
|
|
|
while(mm<n_pt_bound)
|
|
{
|
|
mm++;
|
|
SplitPoints[mm]=new typename std::size_t[N];
|
|
assert(SplitPoints[mm]);
|
|
|
|
epsi_crt=ApproxError[N-1]=ApproxError[N-2];
|
|
SplitPoints[mm][N-1]=N-2;
|
|
|
|
for(j=N-3;j>=mm-2;j--)
|
|
{
|
|
if(ApproxError[j]<=epsi_crt && Graph[j][N-1]<=epsi_crt)
|
|
{
|
|
err_crt=ApproxError[j]+Graph[j][N-1];
|
|
|
|
if(err_crt<ApproxError[N-1])
|
|
{
|
|
epsi_crt=ApproxError[N-1]=err_crt;
|
|
SplitPoints[mm][N-1]=j;
|
|
}
|
|
}
|
|
}
|
|
|
|
if(mm<n_pt_bound)
|
|
for(i=N-2;i>=mm-1;i--)
|
|
{
|
|
epsi_crt=ApproxError[i]=ApproxError[i-1];
|
|
SplitPoints[mm][i]=i-1;
|
|
for(j=i-2;j>=mm-2;j--)
|
|
{
|
|
if(ApproxError[j]<=epsi_crt && Graph[j][i]<=epsi_crt)
|
|
{
|
|
err_crt=ApproxError[j]+Graph[j][i];
|
|
if(err_crt<ApproxError[i])
|
|
{
|
|
epsi_crt=ApproxError[i]=err_crt;
|
|
SplitPoints[mm][i]=j;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
approx_error=ApproxError[N-1];
|
|
|
|
c_i=beyond;
|
|
c_i--;
|
|
|
|
typename std::list<InputIterator> local_container;
|
|
typename std::list<InputIterator>::iterator cri;
|
|
|
|
local_container.push_front(c_i);
|
|
|
|
n=N-1;
|
|
|
|
typename std::size_t nn=N-1;
|
|
|
|
for(m=n_pt_bound;m>2;m--)
|
|
{
|
|
n=SplitPoints[m][n];
|
|
|
|
for(typename std::size_t ii=nn;ii>n;ii--)
|
|
c_i--;
|
|
|
|
nn=n;
|
|
local_container.push_front(c_i);
|
|
}
|
|
|
|
local_container.push_front(begin);
|
|
|
|
for(cri=local_container.begin();cri!=local_container.end();cri++)
|
|
*result++=**cri;
|
|
|
|
delete[] ApproxError;
|
|
|
|
for(m=3;m<=mm;m++)
|
|
delete[] SplitPoints[m];
|
|
delete[] SplitPoints;
|
|
|
|
|
|
for(m=0;m<N;m++)
|
|
delete[] Graph[m];
|
|
delete[] Graph;
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
CGAL_END_NAMESPACE
|
|
|
|
#endif // CGAL_POLYAP_FCT
|