Back to index

texmacs  1.0.7.15
poly_line.cpp
Go to the documentation of this file.
00001 
00002 /******************************************************************************
00003 * MODULE     : poly_line.cpp
00004 * DESCRIPTION: Poly lines
00005 * COPYRIGHT  : (C) 2012  Joris van der Hoeven
00006 *******************************************************************************
00007 * This software falls under the GNU general public license version 3 or later.
00008 * It comes WITHOUT ANY WARRANTY WHATSOEVER. For details, see the file LICENSE
00009 * in the root directory or <http://www.gnu.org/licenses/gpl-3.0.html>.
00010 ******************************************************************************/
00011 
00012 #include "poly_line.hpp"
00013 
00014 inline double square (double x) { return x*x; }
00015 
00016 /******************************************************************************
00017 * Extra routines for points
00018 ******************************************************************************/
00019 
00020 double
00021 min (point p) {
00022   ASSERT (N(p) > 0, "non zero length expected");
00023   double r= p[0];
00024   for (int i=1; i<N(p); i++)
00025     r= min (r, p[i]);
00026   return r;
00027 }
00028 
00029 double
00030 max (point p) {
00031   ASSERT (N(p) > 0, "non zero length expected");
00032   double r= p[0];
00033   for (int i=1; i<N(p); i++)
00034     r= max (r, p[i]);
00035   return r;
00036 }
00037 
00038 double
00039 l2_norm (point p) {
00040   double s= 0.0;
00041   for (int i=0; i<N(p); i++)
00042     s += square (p[i]);
00043   return sqrt (s);
00044 }
00045 
00046 double
00047 inner (point p, point q) {
00048   ASSERT (N(p) == N(q), "unequal lengths");
00049   double s= 0.0;
00050   for (int i=0; i<N(p); i++)
00051     s += p[i] * q[i];
00052   return s;
00053 }
00054 
00055 double
00056 distance (point p, point q) {
00057   ASSERT (N(p) == N(q), "unequal lengths");
00058   double s= 0.0;
00059   for (int i=0; i<N(p); i++)
00060     s += square (q[i] - p[i]);
00061   return sqrt (s);
00062 }
00063 
00064 point
00065 project (point p, point q1, point q2) {
00066   ASSERT (N(p) == N(q1) && N(p) == N(q2), "unequal lengths");
00067   int i, n= N(p);
00068   double s= 0.0, t= 0.0;
00069   for (i=0; i<n; i++) {
00070     s += (q2[i] - q1[i]) * (p[i] - q1[i]);
00071     t += square (q2[i] - q1[i]);
00072   }
00073   double a= s / t;
00074   if (a < 0.0) return q1;
00075   if (a > 1.0) return q2;
00076   return q1 + a * (q2 - q1);
00077 }
00078 
00079 double
00080 distance (point p, point q1, point q2) {
00081   if (q1 == q2) return distance (p, q1);
00082   else return distance (p, project (p, q1, q2));
00083 }
00084 
00085 point
00086 inf (point p, point q) {
00087   ASSERT (N(p) == N(q), "unequal lengths");
00088   point r (N(p));
00089   for (int i=0; i<N(p); i++)
00090     r[i]= min (p[i], q[i]);
00091   return r;
00092 }
00093 
00094 point
00095 sup (point p, point q) {
00096   ASSERT (N(p) == N(q), "unequal lengths");
00097   point r (N(p));
00098   for (int i=0; i<N(p); i++)
00099     r[i]= max (p[i], q[i]);
00100   return r;
00101 }
00102 
00103 point
00104 abs (point p) {
00105   point r (N(p));
00106   for (int i=0; i<N(p); i++)
00107     r[i]= (p[i] > 0? p[i]: -p[i]);
00108   return r;
00109 }
00110 
00111 /******************************************************************************
00112 * Poly lines
00113 ******************************************************************************/
00114 
00115 double
00116 distance (point p, poly_line pl) {
00117   double m= 1.0e10;
00118   if (N(pl) == 1) return distance (p, pl[0]);
00119   for (int i=0; i+1<N(pl); i++)
00120     m= min (m, distance (p, pl[i], pl[i+1]));
00121   return m;
00122 }
00123 
00124 bool
00125 nearby (point p, poly_line pl) {
00126   return distance (p, pl) <= 5.0;
00127 }
00128 
00129 point
00130 inf (poly_line pl) {
00131   ASSERT (N(pl) > 0, "non zero length expected");
00132   point p= pl[0];
00133   for (int i=1; i<N(pl); i++)
00134     p= inf (p, pl[i]);
00135   return p;
00136 }
00137 
00138 point
00139 sup (poly_line pl) {
00140   ASSERT (N(pl) > 0, "non zero length expected");
00141   point p= pl[0];
00142   for (int i=1; i<N(pl); i++)
00143     p= sup (p, pl[i]);
00144   return p;
00145 }
00146 
00147 poly_line
00148 operator + (poly_line pl, point p) {
00149   int i, n= N(pl);
00150   poly_line r (n);
00151   for (i=0; i<n; i++)
00152     r[i]= pl[i] + p;
00153   return r;
00154 }
00155 
00156 poly_line
00157 operator - (poly_line pl, point p) {
00158   int i, n= N(pl);
00159   poly_line r (n);
00160   for (i=0; i<n; i++)
00161     r[i]= pl[i] - p;
00162   return r;
00163 }
00164 
00165 poly_line
00166 operator * (double x, poly_line pl) {
00167   int i, n= N(pl);
00168   poly_line r (n);
00169   for (i=0; i<n; i++)
00170     r[i]= x * pl[i];
00171   return r;
00172 }
00173 
00174 poly_line
00175 normalize (poly_line pl) {
00176   if (N(pl) == 0) return pl;
00177   pl= pl - inf (pl);
00178   if (N(pl) == 1) return pl;
00179   double sc= max (sup (pl));
00180   if (sc == 0.0) return pl;
00181   return (1.0 / sc) * pl;
00182 }
00183 
00184 double
00185 length (poly_line pl) {
00186   double len= 0.0;
00187   for (int i=1; i<N(pl); i++)
00188     len += distance (pl[i], pl[i-1]);
00189   return len;
00190 }
00191 
00192 point
00193 access (poly_line pl, double t) {
00194   if (t < 0) return pl[0];
00195   for (int i=1; i<N(pl); i++) {
00196     point dp= pl[i] - pl[i-1];
00197     double len= l2_norm (dp);
00198     if (t < len) return pl[i-1] + (t / len) * dp;
00199     t -= len;
00200   }
00201   return pl[N(pl)-1];
00202 }
00203 
00204 /******************************************************************************
00205 * Contours
00206 ******************************************************************************/
00207 
00208 double
00209 distance (point p, contours gl) {
00210   double m= 1.0e10;
00211   for (int i=0; i<N(gl); i++)
00212     m= min (m, distance (p, gl[i]));
00213   return m;
00214 }
00215 
00216 bool
00217 nearby (point p, contours gl) {
00218   return distance (p, gl) <= 5.0;
00219 }
00220 
00221 point
00222 inf (contours gl) {
00223   ASSERT (N(gl) > 0, "non zero length expected");
00224   point p= inf (gl[0]);
00225   for (int i=1; i<N(gl); i++)
00226     p= inf (p, inf (gl[i]));
00227   return p;
00228 }
00229 
00230 point
00231 sup (contours gl) {
00232   ASSERT (N(gl) > 0, "non zero length expected");
00233   point p= sup (gl[0]);
00234   for (int i=1; i<N(gl); i++)
00235     p= sup (p, sup (gl[i]));
00236   return p;
00237 }
00238 
00239 contours
00240 operator + (contours gl, point p) {
00241   int i, n= N(gl);
00242   contours r (n);
00243   for (i=0; i<n; i++)
00244     r[i]= gl[i] + p;
00245   return r;
00246 }
00247 
00248 contours
00249 operator - (contours gl, point p) {
00250   int i, n= N(gl);
00251   contours r (n);
00252   for (i=0; i<n; i++)
00253     r[i]= gl[i] - p;
00254   return r;
00255 }
00256 
00257 contours
00258 operator * (double x, contours gl) {
00259   int i, n= N(gl);
00260   contours r (n);
00261   for (i=0; i<n; i++)
00262     r[i]= x * gl[i];
00263   return r;
00264 }
00265 
00266 contours
00267 normalize (contours gl) {
00268   if (N(gl) == 0) return gl;
00269   gl= gl - inf (gl);
00270   double sc= max (sup (gl));
00271   if (sc == 0.0) return gl;
00272   gl= (1.0 / sc) * gl;
00273   return gl;
00274 }
00275 
00276 /******************************************************************************
00277 * Vertex detection
00278 ******************************************************************************/
00279 
00280 array<double>
00281 vertices (poly_line pl) {
00282   pl= (1.0 / length (pl)) * pl;
00283   array<double> r;
00284   double t = 0.0;
00285   double dt= 0.025;
00286   r << 0.0;
00287   int    todo_i= -1;
00288   double todo_t= 0.0;
00289   double todo_p= 0.0;
00290   for (int i=0; i+1<N(pl); i++) {
00291     if (t >= r[N(r)-1] + dt && t <= 1.0 - dt) {
00292       double t1= max (t - dt, 0.000000001);
00293       double t2= min (t + dt, 0.999999999);
00294       point  p = pl[i];
00295       point  p1= access (pl, t1);
00296       point  p2= access (pl, t2);
00297       double pr= inner (p1 - p, p2 - p);
00298       if (pr >= 0 && (todo_i < 0 || pr > todo_p)) {
00299         todo_i= i;
00300         todo_t= t;
00301         todo_p= pr;
00302       }
00303     }
00304     t += distance (pl[i+1], pl[i]);
00305     if (todo_i >= 0 && t >= todo_t + dt) {
00306       r << todo_t;
00307       todo_i= -1;
00308     }
00309   }
00310   r << 1.0;
00311   return r;
00312 }
00313 
00314 /******************************************************************************
00315 * Associate invariants to glyphs
00316 ******************************************************************************/
00317 
00318 void
00319 invariants (poly_line pl, int level,
00320             array<tree>& disc, array<double>& cont)
00321 {
00322   double l= length (pl);
00323   int pieces= 20;
00324   for (int i=0; i<=pieces; i++) {
00325     double t= (0.999999999 * i) / pieces;
00326     point  p= access (pl, t * l);
00327     cont << p;
00328   }
00329 
00330   if (level <= 1) {
00331     array<double> ts= vertices (pl);
00332     disc << tree (as_string (N(ts)));
00333     cont << (2.5 * ts);
00334   }
00335 }
00336 
00337 void
00338 invariants (contours gl, int level,
00339             array<tree>& disc, array<double>& cont)
00340 {
00341   gl= normalize (gl);
00342   disc << tree (as_string (N(gl)));
00343   for (int i=0; i<N(gl); i++)
00344     invariants (gl[i], level, disc, cont);
00345 }