/**
* \mainpage
* \section Usage
* @code
* PowFit2D [-v] [-p maxdegree] inputfilename > outputfile
* @endcode
* The @c inputfilename is a list of function values on
* a finite discrete set of 2D Cartesian coordinates as described in PowFit2D::readf().
* The standard output contains the fitting coefficients as described in
* PowFit2D::display(). If one wants to delete the comment lines that are dispersed
* in the standard output, one can use a pipe like
* @code
* PowFit2D [-v] [-p maxdegree] inputfilename | fgrep -v '#' | awk '{print $1 $2 $3}' > outputfile
* @endcode
* One can extract the original input and the fit by using the -v option
* and looking at the lines
* between FITTED START and FITTED END:
* @code
* PowFit2D -v inputfilename | sed -n '/FITTED START/,/FITTED END/p' | sed '/FITTED/d' | tr "#" " " > outputfile
* @endcode
* and this output is immediately useful for gnuplot(1)'s splot command.
* @code
* gnuplot
* gnuplot> splot "outputfile" title "original", "outputfile" using 1:2:4 title "fitted"
* gnuplot> quit
* @endcode
*
* \subsection Command line options
* See PowFit2D::main() .
*
* \section Description
* See PowFit2D::fit() for the fitting formul.
*
* \section Compilation
* A prerequisite is that the LAPACK library of netlib.org
* has been installed.
* @code
* g++ -o PowFit2D -I PowFit2D.cpp -L -llapack
* @endcode
* The makefile entry is
* @code
* LIBFLAGS=-L
* INCFLAGS=-I
* PowFit2D: PowFit2D.cpp
* $(CXX) -o $@ $(INCFLAGS) $< $(LIBFLAGS) -llapack
* @endcode
* IF the GSL library is available, this ought be indicated by
* setting preprocessor variable HAVE_GSL, then the LAPACK library is not needed
* @code
* g++ -o PowFit2D -I -DHAVE_GSL PowFit2D.cpp -L... -lgslcblas -lgsl
* @endcode
*
* @author Richard J. Mathar
* @todo compute error bars on the fitting coefficients
* @since 2007-07-12
* Richard J. Mathar homepage.
* @see
* \cite{GuoOptik117}
*/
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#ifdef HAVE_GSL
#include
#include
#include
#include
#else /* HAVE_GSL */
extern "C" {
/** The f2c header file f2c.h.
* A prerequisite to include clapack.h below.
* This is obtained from www.netlib.org/f2c.
*/
#include
}
extern "C" {
/** The clapack header file \c clapack.h.
* This is obtained from www.netlib.org/clapack.
*/
#include
}
#endif /* HAVE_GSL */
using namespace std ;
/** The class contains a list of 2D points \f$(x_k,y_k)\f$, \f$k=1,..,N\f$, and scalar
* functions values \f$f_k\f$ of these, and performs a least squares fit on a polynomial
* basis in the sense of \f$f_k\approx \sum_{i,j} \alpha_{i,j} x_k^i y_k^j\f$.
* @since 2007-07-12
*/
class PowFit2D {
public:
/** Ctor.
* This initializes the PowFit2D::x, PowFit2D::y and PowFit2D::f vectors with
* the triples of floating point numbers that are to be fitted.
* @param[in] fname the file name of the input file with the tabulated (x,y,f) triples.
* @since 2007-07-12
*/
PowFit2D(char *fname, const bool verbos=false) : verb(verbos)
{
readf(fname,verbos) ;
}
/** Do the fitting.
* Minimize
* \f$ \sum_k [f(x_k,y_k)-f_k]^2\rightarrow \min \f$.
* with the ansatz
* \f$ f(x,y)= \sum_{p=0}^n \sum_{i=0}^p \alpha_{i,p-i} x^i y^{p-i} \f$.
* which leads to the linear system of equations for the unknown \f$ \alpha_{i,p-i}\f$
* \f$ \sum_{p,i}\alpha_{i,p-i} \sum_k x_k^iy_k^{p-i} x_k^l y_k^m = \sum_k f_k x_k^l y_k^m \f$.
* @param[in] maxP the maximum power in the fitting procedure.
* The default of -1 means that it is chosen automatically to be the maximum
* allowed by the number of 2D points that have been read in.
* @since 2007-07-12
*/
void fit(int maxP = -1)
{
/* the degrees of freedom as determined by the number of points.
*/
const int deg = x.size() ;
if( deg <=0)
{
cerr << "Number of points " << deg << " insufficient\n";
exit(EXIT_FAILURE) ;
}
/* determine maxP automatically. deg=0 does not allow fitting.
* deg=1,2 means maxP=0, deg=3..5 means maxP=1 etc.
*/
if ( maxP < 0 )
{
maxP = 0 ;
while ( triang(maxP+2) <= deg)
maxP++ ;
}
else
{
while ( triang(maxP+1) > deg)
maxP-- ;
}
maxHybrP = maxP ;
/* the number of fitting degrees of freedom
* We keep this an overdetermined system, n<=deg. The case of n=deg
* means the fit is actually an interpolation (ie, exact).
*/
int n=triang(maxP+1) ;
if (verb)
{
cout << "# maximum total power " << maxP << " for " << deg << " input points; "
<< n << " fitting degrees\n" ;
}
#ifdef HAVE_GSL
/* the n by n matrix of the linear system of equations */
gsl_matrix *A=gsl_matrix_alloc(n,n) ;
gsl_matrix_set_all(A,0.) ;
/* the RHS of the linear system of equations */
gsl_vector *rhs=gsl_vector_alloc(n) ;
gsl_vector_set_all(rhs,0.) ;
/* fill the matrix and the right hand side
*/
for(int row=0; rowZernike.txt.
* - For lines that are not comment lines, the first non-negative integer is the power \f$i\f$
* in the x coordinate, the next
* non-negative integer is the power \f$j\f$ in the y coordinate, and the following number
* is the expansion coefficient \f$\alpha_{i,j}\f$ in \f$f\approx \sum_{i,j}\alpha_{i,j}x^iy^j\f$.
* The rest of the line is an expression summarizing this in
* more human-readable fashion.
* @warning this only makes sense if PowFit2D.fit() had been called before.
* @parameter[in] verb verbosity level
* @since 2007-07-12
*/
void display(const bool verb=false)
{
/* Display the main result, which is the expansion coefficients
* and their associated powers in the two cartesian coordinates.
*/
for(int i=0; i < alpha.size() ; i++)
{
int powx, powy ;
strIndxInv(i,powx,powy) ;
cout << powx << " " << powy << " " << alpha[i] << " "
<< alpha[i] << "*x^" << powx << "*y^" << powy << endl ;
}
/* If verbose output is demanded, also show the fit of the quality
* by providing the original triples together with their fitted
* approximation by the bi-variate polynomial and the error between
* the fit and the input value.
*/
if ( verb)
{
cout << "# POWFIT2D FITTED START\n" ;
for(int i=0; i < x.size() ; i++)
{
const double fi = fitted(x[i],y[i]) ;
cout << "# " << x[i] << " " << y[i] << " " << f[i] << " " << fi << " " << fi-f[i] << endl ;
}
cout << "# POWFIT2D FITTED END\n" ;
}
/* Finally convert the expansion coefficients in the Cartesian
* coordinates to expansion coefficients in the spherical system.
*/
cout << "# POWFIT2D SPHERICAL START x=r*cos phi, y=r*sin phi\n" ;
/* loop over the powers r^p
*/
for(int p=0; p <= maxHybrP ; p++)
{
/* loop over the sin(m*phi) with m<0 and the cos(m*phi) with m>=0.
*/
for(int m= p; m >= -p ; m--)
{
double co = 0. ;
/* gather the contribution of all mixed powers individually.
* The powx+powy=p that contribute with x^powx*y^powy are all in a consecutiave
* range of straightened indices.
*/
for(int s= triang(p) ; s < triang(p+1) ; s++)
{
int powx, powy ;
strIndxInv(s,powx,powy) ;
const int degfact = ( m==0) ? 1 : 2 ;
co += degfact*alpha[s]*cosFlatn(powx,powy,m) ;
}
cout << "# " << p << " " << m << " " << co << " " << co << "*r^" << p ;
if ( m > 0)
cout << "*cos(" << m << "*phi)" << endl ;
else if ( m < 0 )
cout << "*sin(" << -m << "*phi)" << endl ;
else
cout << endl ;
}
}
cout << "# POWFIT2D SPHERICAL END\n" ;
}
protected:
private:
/** list of x coordinates
*/
vector x ;
/** list of y coordinates, one per x coordinate
*/
vector y ;
/** list of function values, one per (x,y) pair.
*/
vector f ;
/** verbosity rised if true.
*/
bool verb ;
/** list of expansion coefficients, 1D version (straightened)
*/
vector alpha ;
/** the maximum combined power of the x and y coordinates for the fit
*/
int maxHybrP ;
/** classify the input line as comment or non-comment.
*/
static bool isComment(const string inpli, const regex_t & preg)
{
/* debugging
* cout << __FILE__ << " " << __LINE__ << " " << inpli << endl ;
*/
if ( regexec(&preg,inpli.c_str(),0,0,0) == 0)
return true ;
else
return false ;
}
/** The binomial C(a,b).
* \latexonly
* \cite[3.1.2]{AS}
* \endlatexonly
* @param[in] a the non-negative numerator of the binomial
* @param[in] b the non-negative denomiator, less or equal to a
* @return the binomial coefficient \f$C(a,b)=a!/b!/(a-b)!\f$
* @since 2007-07-12
* @see OEIS A007318
*/
static double binomial(int a, int b)
{
#ifdef HAVE_GSL
return gsl_sf_choose(unsigned(a),unsigned(b)) ;
#else
if ( a<0 || b<0 || b>a)
return 0. ;
if (a-b < b)
b = a-b ;
if ( b == 0 )
return 1. ;
double resul = a ;
while ( b-- > 1)
resul *= double(a-b+1)/double(b) ;
return a ;
#endif /* HAVE_GSL */
}
#if 0
/** sign function.
* @param[in] x the value of the argument
* @return 1 if x>=0, -1 if x<0.
* @since 2007-07-12
*/
inline static int sign(const double x)
{
return (x> 0.) ? 1 (x==0. ? 0 : -1) ;
}
#endif
/** The contribution of \f$ \cos^l(\varphi) \sin^k(\varphi)\f$ to \f$\cos(m \varphi)\f$
* or \f$\sin(m\varphi)\f$. With the Euler formula
* \latexonly
* \cite[4.3]{AS}
* \endlatexonly
* and the binomial expansion
* \latexonly
* \cite[3.1.3]{AS}
* \endlatexonly
* we have
* \f$\cos^l\varphi \sin^k\varphi=(e^{i\varphi}+e^{-i\varphi})^l/(2^l) (e^{i\varphi}-e^{-i\varphi})^k/(2i)^k\f$
* \f$\sim\frac{1}{2^{k+l}}(-)^{\langle k/2\rangle}\sum_{s=0}^l{l \choose s}\sum_{t=0}^k{k \choose t}(-)^t e^{i(2s-l+k-2t)\varphi}\f$
* to be divided by \f$i\f$ if \f$k\f$ is odd. The algorithm is simply to match
* the integer in the exponential with \f$m\f$.
* @param[in] l the power of the cosine
* @param[in] k the power of the sine
* @param[in] m the cycle number of the cosine (m>=0) or the sine (m<0)
* for which the coefficient is created.
* @return the parameter beta in the expansion cos^l phi sin^k phi = .... beta*cos(m*phi)
* or = .. beta*sin(|m|*phi)
* @since 2007-07-12
*/
static double cosFlatn(const int l, const int k, int m)
{
double resul=0. ;
/* for nonzero result, an even k must meet a cosine expansion coeffieicnt
* and an odd k must meet a sine expansion coefficient.
*/
if ( m >=0 && k % 2 == 0 || m < 0 && k % 2)
{
m = abs(m) ;
/* k/2 if k is even , floor(k/2) if k is odd
*/
const int khalf = k /2 ;
if ( (m+l+k) % 2 ==0)
{
for(int s=0 ; s <= l ; s++)
{
const int t= s-(m+l-k)/2 ;
if ( t>=0 && t <= k)
if ( t % 2 )
resul -= binomial(l,s)*binomial(k,t) ;
else
resul += binomial(l,s)*binomial(k,t) ;
}
if ( khalf % 2)
resul *= -1. ;
resul /= pow(2.,l+k) ;
}
}
return resul ;
}
/** Calculate a value of the fitting function.
* @param[in] x Cartesian x coordinate of the 2D point
* @param[in] y Cartesian y coordinate of the 2D point
* @return the mixed polynomial (the fit) at the 2D point.
* @since 2007-07-12
*/
double fitted(const double x, const double y)
{
double f=0. ;
for(int i=0; i < alpha.size() ; i++)
{
int powx, powy ;
strIndxInv(i,powx,powy) ;
f += alpha[i]*pow(x,powx)*pow(y,powy) ;
}
return f ;
}
/** Read the (x,y,f) triples from an ASCII file.
* The syntax of the lines in the ASCII file is as follows:
* - Lines that start with optional white space followed by a hash (#) or
* containing only white space are ignored (a.k.a. comment lines)
* - All other lines start with three floating point numbers separated
* by white space. Letters and numbers in the rest of these lines are
* ignored. The three floating point numbers are the x Cartesian coordinate,
* then the y Cartesian coordinate, then the function value f.
* The count of lines of this form implies the number of input triples.
* @param[in] fname the UNIX file name
* @param[in] verb if true turns on verbose commenting on the inputs
* @since 2007-07-12
*/
void readf(char *fname, bool verb=false)
{
ifstream inf(fname) ;
/* complain if the file is not readable.
*/
if ( !inf)
{
cerr << "Cannot open " << fname << endl ;
exit(EXIT_FAILURE) ;
}
/* a regular expression characterizing comment lines
*/
regex_t preg;
if ( int rege = regcomp(&preg,"(^[[:space:]]*#)|(^[[:space:]]*$)",REG_NOSUB|REG_EXTENDED) )
{
char errbuf[257] ;
regerror(rege,&preg,errbuf,257) ;
cerr << __FILE__ << " " << __LINE__ << " " << errbuf << endl ;
}
/* read the file line by line until EOF.
*/
while( !inf.eof() )
{
/* one line of the input file
*/
string inpli ;
getline(inf,inpli) ;
if ( !isComment(inpli,preg) )
{
double vals[3] ;
istringstream s(inpli) ;
s >> vals[0] >> vals[1] >> vals[2] ;
if ( verb)
cout << "# " << x.size() << " " << vals[0] << " " << vals[1] << " " << vals[2] << endl ;
x.push_back(vals[0]) ;
y.push_back(vals[1]) ;
f.push_back(vals[2]) ;
}
}
regfree(&preg) ;
}
/** Compute triangular number.
* OEIS A000217
* \latexonly
* \cite{EIS}
* \endlatexonly
* @param[in] i
* @return the i'th triangular number, i(i+1)/2.
* @since 2007-07-12
*/
int triang( const int i)
{
return i*(i+1)/2 ;
}
/** Derive a straightened single index from a pair index. The result is
* constructed by reading the infinite 2x2 array of integer indices
* in the first quadrant along anti-diagonals.
* @param[in] i first index of the pair, >=0
* @param[in] j second index of the pair, >=0
* @return a unique number larger or equal to zero. (0,0) is mapped on 0,
* (0,1) on 1, (1,0) on 2, (0,2) on 3, (1,1) on 4, (2,0) on 5 etc. So a group
* of common i+j produces a sequential list of integers, and the smallest i
* are maped on the lower single indices.
* @since 2007-07-12
*/
int strIndx( const int i, const int j)
{
return triang(i+j)+i ;
}
/** Decompose a straightened single index into a pair of indices.
* This is the inverse function of strIndx().
* @param[in] n the straightened single index.
* @param[out] i first index of the pair, >=0
* @param[out] j second index of the pair, >=0
* @since 2007-07-12
*/
void strIndxInv( const int n, int &i, int & j)
{
/* the triangular number that starts a sequence of common sum i+j
*/
int p=0 ;
while ( triang(p) <= n)
p++ ;
p-- ;
i = n-triang(p) ;
j=p-i ;
}
} ;
/** Print a usage information.
* @param[in] argv0 the command name
* @since 2007-07-12
*/
void usage(char *argv0)
{
cout << "usage : " << argv0 << " [-v] [-p maxdegree] inptfile\n" ;
}
/** The main executable.
* @param[in] -v the option \c -v can be used to produce more
* verbose output
* @param[in] -p the option \c -p followed by a non-negative integer indicates
* that the maximum (hybrid) degree of the powers in the fit should be reduced
* to this integer. If not used, the program will calculate a default, which
* is to take the maximum fitting degree which still matches a least-squares
* algorithm. The default means that the degree is 0 for 1 Cartesian input point,
* the degree is 1 for 2 or 3 input points, the degree is 2 for 4 to 6 input
* points etc. This can be seen as always including a full set of azimuthal
* parameters \f$m\f$ for each \f$\sim r^p\cos(m\varphi)\f$ and \f$\sim r^p\sin(m\varphi)\f$,
* \f$0\le m\le p\f$, such that no directional preference is found in the
* algorithm.
* @param[in] fname the single command line argument is the ASCII file name
* with input lines as described in PowFit2D::readf().
* @return 0 if the fitting was apparently successful, -1 if wrong options
* were used or the input file couldn't be opened.
* @since 2007-07-12
*/
int main(int argc, char *argv[])
{
/* command line option.
*/
int c ;
/* Maximum radial power in the fit. Initialized with negative value
* means we automate the calculation.
*/
int maxP = -1 ;
/* verbosity
*/
bool verb=false ;
/* read command line options
*/
while( (c=getopt(argc,argv,"vp:")) != -1 )
{
switch(c)
{
case 'v':
verb=true ;
break ;
case 'p':
maxP = atoi(optarg) ;
break ;
case '?':
usage(argv[0]) ;
return -1 ;
break ;
}
}
if ( optind >= argc)
{
usage(argv[0]) ;
cerr << "Missing file name argument\n" ;
return -1 ;
}
if (verb)
{
time_t now=time(NULL) ;
cout << "# " << getenv("USER") << " " << ctime(&now) << "# " ;
for(c=0 ; c < argc; c++)
cout << argv[c] << " " ;
cout << endl ;
}
/* read the 2D points to be fitted from the ASCII File
*/
PowFit2D fitProbl(argv[optind],verb) ;
/* do the fitting.
*/
fitProbl.fit(maxP) ;
/* show the results
*/
fitProbl.display(verb) ;
return 0 ;
}