src/2geom/numeric/fitting-model.h

   1 /*
   2  * Fitting Models for Geom Types
   3  *
   4  * Authors:
   5  *      Marco Cecchetti <mrcekets at gmail.com>
   6  *
   7  * Copyright 2008  authors
   8  *
   9  * This library is free software; you can redistribute it and/or
  10  * modify it either under the terms of the GNU Lesser General Public
  11  * License version 2.1 as published by the Free Software Foundation
  12  * (the "LGPL") or, at your option, under the terms of the Mozilla
  13  * Public License Version 1.1 (the "MPL"). If you do not alter this
  14  * notice, a recipient may use your version of this file under either
  15  * the MPL or the LGPL.
  16  *
  17  * You should have received a copy of the LGPL along with this library
  18  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
  19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  20  * You should have received a copy of the MPL along with this library
  21  * in the file COPYING-MPL-1.1
  22  *
  23  * The contents of this file are subject to the Mozilla Public License
  24  * Version 1.1 (the "License"); you may not use this file except in
  25  * compliance with the License. You may obtain a copy of the License at
  26  * http://www.mozilla.org/MPL/
  27  *
  28  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  29  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  30  * the specific language governing rights and limitations.
  31  */
  32
  33
  34 #ifndef _NL_FITTING_MODEL_H_
  35 #define _NL_FITTING_MODEL_H_
  36
  37
  38 #include <2geom/d2.h>
  39 #include <2geom/sbasis.h>
  40 #include <2geom/bezier.h>
  41 #include <2geom/bezier-curve.h>
  42 #include <2geom/poly.h>
  43 #include <2geom/ellipse.h>
  44 #include <2geom/circle.h>
  45 #include <2geom/utils.h>
  46
  47
  48 namespace Geom { namespace NL {
  49
  50 /*
  51  * A model is an abstraction for an expression dependent from a parameter where
  52  * the coefficients of this expression are the unknowns of the fitting problem.
  53  * For a ceratain number of parameter values we know the related values
  54  * the expression evaluates to: from each parameter value we get a row of
  55  * the matrix of the fitting problem, from each expression value we get
  56  * the related constant term.
  57  * Example: given the model a*x^2 + b*x + c = 0; from x = 1 we get
  58  * the equation a + b + c = 0, in this example the constant term is always
  59  * the same for each parameter value.
  60  *
  61  * A model is required to implement 3 methods:
  62  *
  63  *  - size : returns the number of unknown coefficients that appear in
  64  *           the expression of the fitting problem;
  65  *  - feed : its input is a parameter value and the related expression value,
  66  *           it generates a matrix row and a new entry of the constant vector
  67  *           of the fitting problem;
  68  *  - instance : it has an input parameter represented by the raw vector
  69  *               solution of the fitting problem and an output parameter
  70  *               of type InstanceType that return a specific object that is
  71  *               generated using the fitting problem solution, in the example
  72  *               above the object could be a Poly type.
  73  */
  74
  75 /*
  76  *   completely unknown models must inherit from this template class;
  77  *   example: the model a*x^2 + b*x + c = 0 to be solved wrt a, b, c;
  78  *   example: the model A(t) = known_sample_value_at(t) to be solved wrt
  79  *       the coefficients of the curve A(t) expressed in S-Basis form;
  80  *   parameter type: the type of x and t variable in the examples above;
  81  *   value type:     the type of the known sample values (in the first example
  82  *                   is constant )
  83  *   instance type:  the type of the objects produced by using
  84  *                   the fitting raw data solution
  85  */
  86 template< typename ParameterType, typename ValueType, typename InstanceType >
  87 class LinearFittingModel
  88 {
  89   public:
  90     typedef ParameterType       parameter_type;
  91     typedef ValueType           value_type;
  92     typedef InstanceType        instance_type;
  93
  94     static const bool WITH_FIXED_TERMS = false;
  95
  96     /*
  97      * a LinearFittingModel must implement the following methods:
  98      *
  99      * void feed( VectorView & vector,
 100      *            parameter_type const& sample_parameter ) const;
 101      *
 102      * size_t size() const;
 103      *
 104      * void instance(instance_type &, raw_type const& raw_data) const;
 105      *
 106      */
 107 };
 108
 109
 110 /*
 111  *   partially known models must inherit from this template class
 112  *   example: the model a*x^2 + 2*x + c = 0 to be solved wrt a and c
 113  */
 114 template< typename ParameterType, typename ValueType, typename InstanceType >
 115 class LinearFittingModelWithFixedTerms
 116 {
 117   public:
 118     typedef ParameterType       parameter_type;
 119     typedef ValueType           value_type;
 120     typedef InstanceType        instance_type;
 121
 122     static const bool WITH_FIXED_TERMS = true;
 123
 124     /*
 125      * a LinearFittingModelWithFixedTerms must implement the following methods:
 126      *
 127      * void feed( VectorView & vector,
 128      *            value_type & fixed_term,
 129      *            parameter_type const& sample_parameter ) const;
 130      *
 131      * size_t size() const;
 132      *
 133      * void instance(instance_type &, raw_type const& raw_data) const;
 134      *
 135      */
 136
 137
 138 };
 139
 140
 141 // incomplete model, it can be inherited to make up different kinds of
 142 // instance type; the raw data is a vector of coefficients of a polynomial
 143 // rapresented in standard power basis
 144 template< typename InstanceType >
 145 class LFMPowerBasis
 146     : public LinearFittingModel<double, double, InstanceType>
 147 {
 148   public:
 149     LFMPowerBasis(size_t degree)
 150         : m_size(degree + 1)
 151     {
 152     }
 153
 154     void feed( VectorView & coeff, double sample_parameter ) const
 155     {
 156         coeff[0] = 1;
 157         double x_i = 1;
 158         for (size_t i = 1; i < coeff.size(); ++i)
 159         {
 160           x_i *= sample_parameter;
 161           coeff[i] = x_i;
 162         }
 163     }
 164
 165     size_t size() const
 166     {
 167         return m_size;
 168     }
 169
 170   private:
 171     size_t m_size;
 172 };
 173
 174
 175 // this model generates Geom::Poly objects
 176 class LFMPoly
 177     : public LFMPowerBasis<Poly>
 178 {
 179   public:
 180     LFMPoly(size_t degree)
 181         : LFMPowerBasis<Poly>(degree)
 182     {
 183     }
 184
 185     void instance(Poly & poly, ConstVectorView const& raw_data) const
 186     {
 187         poly.clear();
 188         poly.resize(size());
 189         for (size_t i = 0; i < raw_data.size(); ++i)
 190         {
 191             poly[i] =  raw_data[i];
 192         }
 193     }
 194 };
 195
 196
 197 // incomplete model, it can be inherited to make up different kinds of
 198 // instance type; the raw data is a vector of coefficients of a polynomial
 199 // rapresented in standard power basis with leading term coefficient equal to 1
 200 template< typename InstanceType >
 201 class LFMNormalizedPowerBasis
 202     : public LinearFittingModelWithFixedTerms<double, double, InstanceType>
 203 {
 204   public:
 205     LFMNormalizedPowerBasis(size_t _degree)
 206         : m_model( _degree - 1)
 207     {
 208         assert(_degree > 0);
 209     }
 210
 211
 212     void feed( VectorView & coeff,
 213                double & known_term,
 214                double sample_parameter ) const
 215     {
 216         m_model.feed(coeff, sample_parameter);
 217         known_term = coeff[m_model.size()-1] * sample_parameter;
 218     }
 219
 220     size_t size() const
 221     {
 222         return m_model.size();
 223     }
 224
 225   private:
 226     LFMPowerBasis<InstanceType> m_model;
 227 };
 228
 229
 230 // incomplete model, it can be inherited to make up different kinds of
 231 // instance type; the raw data is a vector of coefficients of the equation
 232 // of an ellipse curve
 233 template< typename InstanceType >
 234 class LFMEllipseEquation
 235     : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>
 236 {
 237   public:
 238     void feed( VectorView & coeff, double & fixed_term, Point const& p ) const
 239     {
 240         coeff[0] = p[X] * p[Y];
 241         coeff[1] = p[Y] * p[Y];
 242         coeff[2] = p[X];
 243         coeff[3] = p[Y];
 244         coeff[4] = 1;
 245         fixed_term = p[X] * p[X];
 246     }
 247
 248     size_t size() const
 249     {
 250         return 5;
 251     }
 252 };
 253
 254
 255 // this model generates Ellipse curves
 256 class LFMEllipse
 257     : public LFMEllipseEquation<Ellipse>
 258 {
 259   public:
 260     void instance(Ellipse & e, ConstVectorView const& coeff) const
 261     {
 262         e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);
 263     }
 264 };
 265
 266
 267 // incomplete model, it can be inherited to make up different kinds of
 268 // instance type; the raw data is a vector of coefficients of the equation
 269 // of a circle curve
 270 template< typename InstanceType >
 271 class LFMCircleEquation
 272     : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>
 273 {
 274   public:
 275     void feed( VectorView & coeff, double & fixed_term, Point const& p ) const
 276     {
 277         coeff[0] = p[X];
 278         coeff[1] = p[Y];
 279         coeff[2] = 1;
 280         fixed_term = p[X] * p[X] + p[Y] * p[Y];
 281     }
 282
 283     size_t size() const
 284     {
 285         return 3;
 286     }
 287 };
 288
 289
 290 // this model generates Ellipse curves
 291 class LFMCircle
 292     : public LFMCircleEquation<Circle>
 293 {
 294   public:
 295     void instance(Circle & c, ConstVectorView const& coeff) const
 296     {
 297         c.set(1, coeff[0], coeff[1], coeff[2]);
 298     }
 299 };
 300
 301
 302 // this model generates SBasis objects
 303 class LFMSBasis
 304     : public LinearFittingModel<double, double, SBasis>
 305 {
 306   public:
 307     LFMSBasis( size_t _order )
 308         : m_size( 2*(_order+1) ),
 309           m_order(_order)
 310     {
 311     }
 312
 313     void feed( VectorView & coeff, double t ) const
 314     {
 315         double u0 = 1-t;
 316         double u1 = t;
 317         double s = u0 * u1;
 318         coeff[0] = u0;
 319         coeff[1] = u1;
 320         for (size_t i = 2; i < size(); i+=2)
 321         {
 322             u0 *= s;
 323             u1 *= s;
 324             coeff[i] = u0;
 325             coeff[i+1] = u1;
 326         }
 327     }
 328
 329     size_t size() const
 330     {
 331         return m_size;
 332     }
 333
 334     void instance(SBasis & sb, ConstVectorView const& raw_data) const
 335     {
 336         sb.clear();
 337         sb.resize(m_order+1);
 338         for (unsigned int i = 0, k = 0; i < raw_data.size(); i+=2, ++k)
 339         {
 340             sb[k][0] = raw_data[i];
 341             sb[k][1] = raw_data[i+1];
 342         }
 343     }
 344
 345   private:
 346     size_t m_size;
 347     size_t m_order;
 348 };
 349
 350
 351 // this model generates D2<SBasis> objects
 352 class LFMD2SBasis
 353     : public LinearFittingModel< double, Point, D2<SBasis> >
 354 {
 355   public:
 356     LFMD2SBasis( size_t _order )
 357         : mosb(_order)
 358     {
 359     }
 360
 361     void feed( VectorView & coeff, double t ) const
 362     {
 363         mosb.feed(coeff, t);
 364     }
 365
 366     size_t size() const
 367     {
 368         return mosb.size();
 369     }
 370
 371     void instance(D2<SBasis> & d2sb, ConstMatrixView const& raw_data) const
 372     {
 373         mosb.instance(d2sb[X], raw_data.column_const_view(X));
 374         mosb.instance(d2sb[Y], raw_data.column_const_view(Y));
 375     }
 376
 377   private:
 378     LFMSBasis mosb;
 379 };
 380
 381
 382 // this model generates Bezier objects
 383 class LFMBezier
 384     : public LinearFittingModel<double, double, Bezier>
 385 {
 386   public:
 387     LFMBezier( size_t _order )
 388         : m_size(_order + 1),
 389           m_order(_order)
 390     {
 391         binomial_coefficients(m_bc, m_order);
 392     }
 393
 394     void feed( VectorView & coeff, double t ) const
 395     {
 396         double s = 1;
 397         for (size_t i = 0; i < size(); ++i)
 398         {
 399             coeff[i] = s * m_bc[i];
 400             s *= t;
 401         }
 402         double u = 1-t;
 403         s = 1;
 404         for (size_t i = size()-1; i > 0; --i)
 405         {
 406             coeff[i] *= s;
 407             s *= u;
 408         }
 409         coeff[0] *= s;
 410     }
 411
 412     size_t size() const
 413     {
 414         return m_size;
 415     }
 416
 417     void instance(Bezier & b, ConstVectorView const& raw_data) const
 418     {
 419         assert(b.size() == raw_data.size());
 420         for (unsigned int i = 0; i < raw_data.size(); ++i)
 421         {
 422             b[i] = raw_data[i];
 423         }
 424     }
 425
 426   private:
 427     size_t m_size;
 428     size_t m_order;
 429     std::vector<size_t> m_bc;
 430 };
 431
 432
 433 // this model generates Bezier curves
 434 template< unsigned int N >
 435 class LFMBezierCurve
 436     : public LinearFittingModel< double, Point, BezierCurve<N> >
 437 {
 438   public:
 439     LFMBezierCurve( size_t _order )
 440         : mob(_order)
 441     {
 442     }
 443
 444     void feed( VectorView & coeff, double t ) const
 445     {
 446         mob.feed(coeff, t);
 447     }
 448
 449     size_t size() const
 450     {
 451         return mob.size();
 452     }
 453
 454     void instance(BezierCurve<N> & bc, ConstMatrixView const& raw_data) const
 455     {
 456         Bezier bx(size()-1);
 457         Bezier by(size()-1);
 458         mob.instance(bx, raw_data.column_const_view(X));
 459         mob.instance(by, raw_data.column_const_view(Y));
 460         bc = BezierCurve<N>(bx, by);
 461     }
 462
 463   private:
 464     LFMBezier mob;
 465 };
 466
 467 }  // end namespace NL
 468 }  // end namespace Geom
 469
 470
 471 #endif // _NL_FITTING_MODEL_H_
 472
 473
 474 /*
 475   Local Variables:
 476   mode:c++
 477   c-file-style:"stroustrup"
 478   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 479   indent-tabs-mode:nil
 480   fill-column:99
 481   End:
 482 */
 483 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :