• Skip to content
  • Skip to link menu
KDE API Reference
  • KDE API Reference
  • kdeedu API Reference
  • KDE Home
  • Contact Us
 

kig

  • sources
  • kde-4.14
  • kdeedu
  • kig
  • objects
curve_imp.cc
Go to the documentation of this file.
1 // Copyright (C) 2002 Dominique Devriese <devriese@kde.org>
2 
3 // This program is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU General Public License
5 // as published by the Free Software Foundation; either version 2
6 // of the License, or (at your option) any later version.
7 
8 // This program is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU General Public License for more details.
12 
13 // You should have received a copy of the GNU General Public License
14 // along with this program; if not, write to the Free Software
15 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16 // 02110-1301, USA.
17 
18 #include "curve_imp.h"
19 #include "../misc/common.h"
20 #include "../misc/coordinate.h"
21 #include "../misc/kignumerics.h"
22 #include "../misc/equation.h"
23 #include "../kig/kig_document.h"
24 
25 #include <cmath>
26 
27 const ObjectImpType* CurveImp::stype()
28 {
29  static const ObjectImpType t(
30  Parent::stype(), "curve",
31  I18N_NOOP( "curve" ),
32  I18N_NOOP( "Select this curve" ),
33  I18N_NOOP( "Select curve %1" ),
34  I18N_NOOP( "Remove a Curve" ),
35  I18N_NOOP( "Add a Curve" ),
36  I18N_NOOP( "Move a Curve" ),
37  I18N_NOOP( "Attach to this curve" ),
38  I18N_NOOP( "Show a Curve" ),
39  I18N_NOOP( "Hide a Curve" )
40  );
41  return &t;
42 }
43 
44 Coordinate CurveImp::attachPoint() const
45 {
46  return Coordinate::invalidCoord();
47 }
48 
49 /*
50  * Generic algorithm for getParam()
51  */
52 
58 double CurveImp::getParamofmin( double a, double b,
59  const Coordinate& p,
60  const KigDocument& doc ) const
61 {
62  double epsilons = 1.e-08;
63  double epsilonl = 2.e-02;
64 
65  //assert( a < b && a >= 0. && b <= 1.0);
66  assert( a < b && a >= 0.);
67 
68  double r2 = ( sqrt( 5. ) - 1 ) / 2.; // golden ratio
69  double r1 = 1. - r2;
70 
71  double t2 = a + r2 * ( b - a );
72  double t1 = a + r1 * ( b - a );
73  Coordinate p1 = getPoint( t1, doc);
74  double f1 = (p1 - p).length();
75  Coordinate p2 = getPoint( t2, doc);
76  double f2 = (p2 - p).length();
77 
78  double fmin, tmin;
79  if (f1 < f2)
80  {
81  b = t2;
82  fmin = f1;
83  tmin = t1;
84  }
85  else
86  {
87  a = t1;
88  fmin = f2;
89  tmin = t2;
90  }
91 
92  while ( ( b - a ) > epsilons &&
93  ( (p1 - p2).length() > 0.4 * fmin
94  || (b - a) > epsilonl) &&
95  fmin > 1.e-8 )
96  {
97  if ( f1 < f2 )
98  {
99  t2 = t1;
100  t1 = a + r1*(b - a);
101  f2 = f1;
102  //p1 = getPoint( fmod( t1, 1. ), doc);
103  p1 = getPoint( t1, doc);
104  f1 = (p1 - p).length();
105  }
106  else
107  {
108  t1 = t2;
109  t2 = a + r2*(b - a);
110  f1 = f2;
111  // p2 = getPoint( fmod( t2, 1. ), doc);
112  p2 = getPoint( t2, doc);
113  f2 = (p2 - p).length();
114  }
115  if ( f1 < f2 )
116  {
117  b = t2;
118  fmin = f1;
119  tmin = t1;
120  }
121  else
122  {
123  a = t1;
124  fmin = f2;
125  tmin = t2;
126  }
127  }
128 
129  return(tmin);
130 }
131 
137 double CurveImp::getDist(double param, const Coordinate& p, const KigDocument& doc) const
138 {
139 // param = fmod( param, 1 );
140 // if( param < 0 ) param += 1.;
141  if( param < 0 ) param = 0.;
142  if( param > 1 ) param = 1.;
143  Coordinate p1 = getPoint( param, doc );
144  // i don't think the p1.valid() switch is really necessary, but I
145  // prefer to not take any chances :)
146  return p1.valid() ? ( p1 - p ).length() : +double_inf;
147 }
148 
149 double CurveImp::getParam( const Coordinate& p, const KigDocument& doc ) const
150 {
151  // this function ( and related functions like getInterval etc. ) is
152  // written by Franco Pasquarelli <pasqui@dmf.bs.unicatt.it>.
153  // I ( domi ) have adapted and documented it a bit.
154 
155  // mp: the following two lines are especially useful in conjunction to
156  // differential geometry constructions like tangent, center of curvature,...
157  // such constructions need to recover the param associated to a (constrained)
158  // PointImp, but do not have direct access to it since it is a parent of the
159  // calcer accociated to the ConstrainedPointType, whereas we only have the
160  // ObjectImps of the Curve and of the Point; in such case the only possibility
161  // consists in a call to getParam, which is unnecessarily heavy since the PointImp
162  // was itself computed previously using getPoint. So the param used in getPoint
163  // is cached in LocusImp, BezierImp, ... and then checked for validity here.
164 
165  if ( doc.mcachedparam >= 0. && doc.mcachedparam <= 1. &&
166  getPoint ( doc.mcachedparam, doc ) == p )
167  return doc.mcachedparam;
168 
169  // consider the function that returns the distance for a point at
170  // parameter x to the locus for a given parameter x. What we do
171  // here is look for the global minimum of this function. We do that
172  // by dividing the range ( [0,1] ) into N parts, and start looking
173  // for a local minimum from there on. If we find one, we keep it if
174  // it is the lowest of all the ones we've already found..
175 
176  const int N = 64;
177  const double incr = 1. / (double) N;
178 
179  // xm is the best parameter we've found so far, fxm is the distance
180  // to the locus from that point. We start with some
181  // pseudo-values.
182  // (mp) note that if the distance is actually increasing in the
183  // whole interval [0,1] this value will be returned in the end.
184  double xm = 0.;
185  double fxm = getDist( xm, p, doc );
186  double x1,x2;
187 
188  double mm[N+1];
189  mm[0] = fxm;
190 
191  for(int j = 1; j < N + 1; j++)
192  {
193  // [x1,x2] is the range we're currently considering..
194  x1 = j * incr;
195 
196  // check the range x1,x2 for the first local maximum..
197  double mm1 = getDist( x1, p, doc);
198  if( mm1 < fxm)
199  {
200  xm=x1;
201  fxm=mm1;
202  }
203  mm[j]=mm1;
204  }
205  if ( xm == 0.)
206  {
207  x1=0.;
208  x2=incr;
209  }
210  else if ( xm >= 1.)
211  {
212  x1=1.-incr;
213  x2=1.;
214  }
215  else
216  {
217  x1=xm-incr;
218  x2=xm+incr;
219  }
220  double xm1 = getParamofmin( x1, x2, p, doc);
221  double fxm1 = getDist( xm1, p, doc );
222  if( fxm1 < fxm )
223  {
224  // we found a new minimum, save it..
225  xm=xm1;
226  fxm=fxm1;
227  }
228  for(int j = 1; j < N - 1; j++)
229  {
230  if (mm[j] < mm[j-1] && mm[j] < mm[j+1])
231  {
232  if ( fxm > 2*mm[j]-mm[j-1] || fxm > 2*mm[j]-mm[j+1])
233  {
234  xm1 = getParamofmin( (j-1)*incr, (j+1)*incr, p, doc);
235  fxm1 = getDist( xm1, p, doc );
236  if( fxm1 < fxm )
237  {
238  // we found a new minimum, save it..
239  xm=xm1;
240  fxm=fxm1;
241  }
242  }
243  }
244  }
245  return xm;
246 }
247 
248 //This function is used to obtain a pseudo-random number using bitwise operators
249 //it probably should be moved elsewhere, or made completely local...
250 //
251 double CurveImp::revert (int n) const
252 {
253  int nl = n;
254  double b = 1.0;
255  double t = 0.0;
256 
257  assert (n > 0);
258  while (nl)
259  {
260  b /= 2;
261  if (nl & 1)
262  t += b;
263  nl >>= 1;
264  }
265  t += b/2 - b*(rand() / (RAND_MAX+1.0));
266  assert(t <1 && t >0);
267  return (t);
268 }
269 
270 QString CurveImp::cartesianEquationString (const KigDocument& doc ) const
271 {
272  EquationString ret = EquationString( "" );
273  const double threshold = 1e-10;
274  const int degmax=6;
275  const int N = (degmax + 1)*(degmax + 2)/2;
276  const int M = N - 1;
277  double rows[M][N];
278  double *matrix[M];
279  double sol[N];
280  double dist;
281  int deg, degx, degy, deglocus;
282  int i=0, j=0;
283  bool test=true;
284  double exp[degmax+2];
285  double x, y;
286  double xi, yi;
287  double t;
288  double f, fx, fy;
289  int m, n, k, idx;
290  int exchange[N];
291  bool needsign = false;
292  Coordinate point[M];
293  Coordinate p;
294  for(m = 0; m < M; m++) matrix[m] = rows[m];
295  for(deglocus = 1; deglocus <= degmax; ++deglocus)
296  {
297  n=((deglocus+1)*(deglocus+2))/2;
298  m=n-1;
299  k=1;
300  for(i=0; i<m; ++i)
301  {
302  do{
303  point[i] = getPoint(revert(k++), doc);
304  } while(point[i].valid() == false); // now we have a valid point
305  xi = point[i].x;
306  yi = point[i].y;
307  j=0;
308  matrix[i][j]=1;
309  for(deg=1; deg<=deglocus; ++deg)
310  {
311  for(degy=0; degy<=deg; ++degy)
312  {
313  j++;
314  degx=deg-degy;
315  if(degx==0)
316  {
317  matrix[i][j]=matrix[i][j-deg-1]*yi;
318  }
319  else
320  {
321  matrix[i][j]=matrix[i][j-deg]*xi;
322  }
323  }
324  }
325  }
326  assert(j==m);
327  GaussianElimination( matrix, m, n, exchange );
328  BackwardSubstitution( matrix, m, n, exchange, sol);
329  for(i=0; i<m; ++i)
330  {
331  test = true;
332  do{
333  t=revert(k++);
334  p = getPoint(t,doc);
335  } while(p.valid() == false);
336  x = p.x;
337  y = p.y;
338  exp[0]=1;
339  exp[1]=1;
340  fx=sol[1];
341  fy=sol[2];
342  f= sol[0];
343  //now we test if p is a point of the curve
344  idx=1;
345  for (deg=1; deg<=deglocus; ++deg){
346  for(degy=0; degy<=deg; ++degy)
347  {
348  degx = deg - degy;
349  if(degy!=deg)
350  {
351  exp[degy] *= x;
352  }
353  else
354  {
355  exp[degy] *= y;
356  exp[degy+1]=exp[degy];
357  }
358  if(deg!= deglocus)
359  {
360 // i_col(dx, dy) =
361 // d = dx+dy;
362 // icol = (d*(d+1)/2) + dy
363 //
364 // qui ci sarebbe i_col (degx+1, degy)
365 // quindi degnuovo = degvecchio + 1, ovvero
366 // il nuovo valore di d*(d+1)/2 e'
367 // (dv+2)*(dv+1)/2 = dv*(dv+1)/2 + (dv+1)
368 // quindi sarebbe giusto "deg+1+idx", se idx non fosse stato incrementato!
369  fx += (degx+1)*exp[degy]*sol[deg+idx+1];
370 // qui ci sarebbe i_col (degx, degy+1)
371  fy += (degy+1)*exp[degy]*sol[deg+idx+2];
372  }
373  f += exp[degy]*sol[idx++];
374  }
375  }
376  dist=fabs(f)/(fabs(fx) + fabs(fy));
377  if (dist > threshold){
378  test=false;
379  break;
380  }
381  }
382  if(test==true) // now we can costruct the cartesian equation of the locus
383  {
384  assert ( deglocus >= 0 && deglocus <= degmax );
385  for (deg = deglocus; deg > 0; --deg)
386  {
387  j = deg*(deg+1)/2;
388  for (degy = 0; degy <= deg; ++degy)
389  {
390  degx = deg - degy;
391  ret.addTerm( sol[j++], ret.xnym(degx,degy), needsign );
392  }
393  }
394  ret.addTerm( sol[0], "", needsign );
395  ret.append( " = 0" );
396  return ret;
397  }
398  }
399  ret = i18n("Possibly trascendental curve");
400  return ret;
401 }
402 
ObjectImpType
Instances of this class represent a certain ObjectImp type.
Definition: object_imp.h:95
QString::append
QString & append(QChar ch)
CurveImp::attachPoint
Coordinate attachPoint() const
Returns a reference point where to attach labels; when this returns an invalidCoord then the attachme...
Definition: curve_imp.cc:44
EquationString
Simple class that represents an equation.
Definition: equation.h:49
CurveImp::cartesianEquationString
QString cartesianEquationString(const KigDocument &w) const
Definition: curve_imp.cc:270
ObjectImp::stype
static const ObjectImpType * stype()
The ObjectImpType representing the base ObjectImp class.
Definition: object_imp.cc:284
Coordinate
The Coordinate class is the basic class representing a 2D location by its x and y components...
Definition: coordinate.h:33
CurveImp::getPoint
virtual const Coordinate getPoint(double param, const KigDocument &) const =0
GaussianElimination
bool GaussianElimination(double *matrix[], int numrows, int numcols, int exchange[])
Gaussian Elimination.
Definition: kignumerics.cpp:271
EquationString::addTerm
void addTerm(double coeff, const QString &unknowns, bool &needsign)
Definition: equation.cc:41
QString
CurveImp::getParamofmin
double getParamofmin(double a, double b, const Coordinate &p, const KigDocument &doc) const
This function calculates the parameter of the point that realizes the minimum in [a,b] of the distance between the points of the locus and the point of coordinate p, using the golden ration method.
Definition: curve_imp.cc:58
EquationString::xnym
const QString xnym(int n, int m) const
Definition: equation.cc:123
Coordinate::invalidCoord
static Coordinate invalidCoord()
Create an invalid Coordinate.
Definition: coordinate.cpp:171
double_inf
const double double_inf
Definition: common.cpp:509
CurveImp::stype
static const ObjectImpType * stype()
Returns the ObjectImpType representing the CurveImp type.
Definition: curve_imp.cc:27
KigDocument
KigDocument is the class holding the real data in a Kig document.
Definition: kig_document.h:36
KigDocument::mcachedparam
double mcachedparam
Definition: kig_document.h:71
BackwardSubstitution
void BackwardSubstitution(double *matrix[], int numrows, int numcols, int exchange[], double solution[])
Definition: kignumerics.cpp:340
Coordinate::x
double x
X Component.
Definition: coordinate.h:126
Coordinate::y
double y
Y Component.
Definition: coordinate.h:129
ObjectImp::valid
bool valid() const
Returns true if this is a valid ObjectImp.
Definition: object_imp.cc:41
CurveImp::getDist
double getDist(double param, const Coordinate &p, const KigDocument &doc) const
This function returns the distance between the point with parameter param and point p...
Definition: curve_imp.cc:137
curve_imp.h
CurveImp::getParam
virtual double getParam(const Coordinate &point, const KigDocument &) const
Definition: curve_imp.cc:149
Coordinate::valid
bool valid() const
Return whether this is a valid Coordinate.
Definition: coordinate.cpp:176
This file is part of the KDE documentation.
Documentation copyright © 1996-2020 The KDE developers.
Generated on Mon Jun 22 2020 13:12:05 by doxygen 1.8.7 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.

kig

Skip menu "kig"
  • Main Page
  • Namespace List
  • Namespace Members
  • Alphabetical List
  • Class List
  • Class Hierarchy
  • Class Members
  • File List
  • File Members

kdeedu API Reference

Skip menu "kdeedu API Reference"
  • Analitza
  •     lib
  • kalgebra
  • kalzium
  •   libscience
  • kanagram
  • kig
  •   lib
  • klettres
  • marble
  • parley
  • rocs
  •   App
  •   RocsCore
  •   VisualEditor
  •   stepcore

Search



Report problems with this website to our bug tracking system.
Contact the specific authors with questions and comments about the page contents.

KDE® and the K Desktop Environment® logo are registered trademarks of KDE e.V. | Legal