~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina5/tina-libs/tina/geometry/geomCurve_con_util.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**********
  2  * 
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU Lesser General Public License as 
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA is distributed in the hope that it will be useful,
 11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
 12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13  * GNU Lesser General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU Lesser General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc., 
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  **********
 20  * 
 21  * Program :    TINA
 22  * File    :  $Source: /home/tina/cvs/tina-libs/tina/geometry/geomCurve_con_util.c,v $
 23  * Date    :  $Date: 2008/10/30 09:43:59 $
 24  * Version :  $Revision: 1.2 $
 25  * CVS Id  :  $Id: geomCurve_con_util.c,v 1.2 2008/10/30 09:43:59 neil Exp $
 26  *
 27  * Author  : Legacy TINA
 28  *
 29  * Notes : utility functions for conics
 30  *
 31  *********
 32 */
 33 
 34 #include "geomCurve_con_util.h"
 35 
 36 #if HAVE_CONFIG_H
 37   #include <config.h>
 38 #endif
 39 
 40 #include <math.h>
 41 #include <tina/sys/sysDef.h>
 42 #include <tina/sys/sysPro.h>
 43 #include <tina/math/mathDef.h>
 44 #include <tina/math/mathPro.h>
 45 #include <tina/geometry/geomDef.h>
 46 #include <tina/geometry/geom_CurveDef.h>
 47 
 48 
 49 /**
 50 scale coefficients so  a+c=1
 51 **/
 52 void    conic_normalise(Conic * conic)
 53 {
 54     double  k;
 55 
 56     if (conic == NULL)
 57         return;
 58 
 59     k = 1.0 / (conic->a + conic->c);
 60 
 61     conic->a *= k;
 62     conic->b *= k;
 63     conic->c *= k;
 64     conic->d *= k;
 65     conic->e *= k;
 66     conic->f *= k;
 67 }
 68 
 69 /**
 70 
 71 evaluate conic quadratic form at point (x y)
 72 **/
 73 double  conic_eval(Conic * conic, Vec2 p)
 74 {
 75     double  a = conic->a, b = conic->b, c = conic->c;
 76     double  d = conic->d, e = conic->e, f = conic->f;
 77     double  x = vec2_x(p), y = vec2_y(p);
 78 
 79     return (a * x * x + 2.0 * b * x * y + c * y * y + 2.0 * d * x + 2.0 * e * y + f);
 80 }
 81 
 82 /**
 83 
 84 evaluate conic gradient vector at point (x y)
 85 **/
 86 Vec2    conic_grad(Conic * conic, Vec2 p)
 87 {
 88     double  x = vec2_x(p), y = vec2_y(p);
 89 
 90     return (vec2(2.0 * (conic->a * x + conic->b * y + conic->d),
 91                  2.0 * (conic->b * x + conic->c * y + conic->e)));
 92 }
 93 
 94 /**
 95 discriminant of conic - positive for ellipse, negative for hyperbola
 96 **/
 97 double  conic_discrim(Conic * conic)
 98 {
 99     double  a = conic->a, b = conic->b, c = conic->c;
100 
101     return (a * c - b * b);
102 }
103 
104 /**
105 eigendirection and eigenvalues of symmetric 2by2 matrix
106 |a b|
107 |b c|
108 **/
109 void    conic_eigen(double a, double b, double c, double *theta, double *lambda1, double *lambda2)
110 {
111     double  p = (a + c) / 2.0, q = (a - c) / 2.0, r;
112 
113     r = sqrt(q * q + b * b);
114     if (b == 0.0 && q == 0.0)
115         *theta = 0.0;
116     else
117         *theta = 0.5 * atan2(b, q);
118     *lambda1 = p + r;
119     *lambda2 = p - r;
120 }
121 
122 /**
123 
124 value of  lambda  such that  lambda*C0+(1-lambda)*C1 = 0  at  (x y)
125 **/
126 double  conic_lambda(Conic * conic0, Conic * conic1, Vec2 p)
127 {
128     double  c0 = conic_eval(conic0, p);
129     double  c1 = conic_eval(conic1, p);
130 
131     return (c1 / (c1 - c0));
132 }
133 
134 /**
135 point on ellipse with angular parameter t
136 **/
137 Vec2    ellipse_point(Conic * conic, double t)
138 {
139     double  c, s, x, y;
140 
141     c = cos(conic->theta);
142     s = sin(conic->theta);
143     x = conic->alpha * cos(t);
144     y = conic->beta * sin(t);
145     return (vec2_sum(conic->center, vec2(x * c - y * s, x * s + y * c)));
146 }
147 
148 /**
149 angular parameter of point near ellipse
150 
151 return values in [0, twopi)
152 **/
153 double  ellipse_param(Conic * conic, Vec2 p)
154 {
155     double  c, s, x, y;
156     double  t;
157 
158     p = vec2_diff(p, conic->center);
159     c = cos(conic->theta);
160     s = sin(conic->theta);
161     x = vec2_x(p) * c + vec2_y(p) * s;
162     y = -vec2_x(p) * s + vec2_y(p) * c;
163 
164     t = atan2(conic->alpha * y, conic->beta * x);
165     if (t < 0.0)
166         t += TWOPI;
167     return (t);
168 }
169 
170 /**
171 point on hyperbola with angular parameter t
172 **/
173 Vec2    hyperbola_point(Conic * conic, double t)
174 {
175     double  c, s, x, y, xp, yp;
176     Vec2 p;
177 
178     c = cos(conic->theta);
179     s = sin(conic->theta);
180     x = conic->branch * conic->alpha * cosh(t);
181     y = conic->beta * sinh(t);
182     xp = x * c - y * s;
183     yp =  x * s + y * c;
184     p = vec2_sum(conic->center, vec2(xp, yp));
185     return (p);
186 }
187 
188 /**
189 angular parameter of point near hyperbola
190 **/
191 double  hyperbola_param(Conic * conic, Vec2 p)
192 {
193     double  c, s;
194 
195 
196     p = vec2_diff(p, conic->center);
197     c = cos(conic->theta);
198     s = sin(conic->theta);
199     return (asinh((-vec2_x(p) * s + vec2_y(p) * c) / conic->beta));
200 }
201 
202 int     hyperbola_branch(Conic * conic, Vec2 p)
203 {
204     double  c, s;
205 
206     p = vec2_diff(p, conic->center);
207     c = cos(conic->theta);
208     s = sin(conic->theta);
209     return (((vec2_x(p) * c + vec2_y(p) * s) > 0.0) ? 1 : -1);
210 }
211 
212 /**
213 approximate squared distance of point from conic
214 **/
215 double  conic_approx_sqrdist(Conic * conic, Vec2 p)
216 {
217     double  a = conic->a, b = conic->b, c = conic->c;
218     double  d = conic->d, e = conic->e, f = conic->f;
219     double  q0, qx, qy;
220     double  x = vec2_x(p), y = vec2_y(p);
221 
222     q0 = a * x * x + 2.0 * b * x * y + c * y * y + 2.0 * d * x + 2.0 * e * y + f;
223     qx = 2.0 * (a * x + b * y + d);
224     qy = 2.0 * (b * x + c * y + e);
225     return (q0 * q0 / (qx * qx + qy * qy));
226 }
227 
228 /**
229 set type, orientation and semiaxes of conic and return center
230 **/
231 void    conic_setup(Conic * conic)
232 {
233     double  a, b, c, d, e;
234     double  lambda1, lambda2;
235     double  q0, disc;
236 
237     if (conic == NULL)
238         return;
239 
240     a = conic->a;
241     b = conic->b;
242     c = conic->c;
243     d = conic->d;
244     e = conic->e;
245 
246     disc = conic_discrim(conic);
247     if (disc > 0.0)
248     {
249         conic->type = ELLIPSE;
250         conic->t1 = 0.0;
251         conic->t2 = TWOPI;
252     } else if (disc < 0.0)
253     {
254         conic->type = HYPERBOLA;
255         conic->t1 = -2.0;
256         conic->t2 = 2.0;
257         conic->branch = 1;
258     } else
259     {
260         conic->type = DEGENERATE;
261         return;
262     }
263     conic->center = vec2(-(c * d - b * e) / disc, -(a * e - b * d) / disc);
264 
265     q0 = conic_eval(conic, conic->center);
266     conic_eigen(a, b, c, &conic->theta, &lambda1, &lambda2);
267     if (lambda1 == 0.0 || lambda2 == 0.0 || q0 == 0.0)
268     {
269         conic->type = DEGENERATE;
270         return;
271     }
272     lambda1 /= -q0;
273     lambda2 /= -q0;
274     conic->alpha = 1.0 / sqrt(fabs(lambda1));
275     conic->beta = 1.0 / sqrt(fabs(lambda2));
276 }
277 
278 Vec2    conic_point(Conic * conic, double t)
279 {
280     Vec2    vec_2 = {Vec2_id};
281 
282     switch (conic->type)
283     {
284     case ELLIPSE:
285         vec_2 = (ellipse_point(conic, t));
286         break;
287     case HYPERBOLA:
288         vec_2 = (hyperbola_point(conic, t));
289         break;
290     case DEGENERATE:
291     default:
292         vec_2 = (vec2_zero());
293     }
294     return vec_2;
295 }
296 
297 double  conic_param(Conic * conic, Vec2 p)
298 {
299     switch (conic->type)
300     {
301         case ELLIPSE:
302         return (ellipse_param(conic, p));
303     case HYPERBOLA:
304         return (hyperbola_param(conic, p));
305     case DEGENERATE:
306     default:
307         return (0.0);
308     }
309 }
310 
311 static int ordered(double t1, double t2, double t3)
312 {
313     if( ((t1 < t2) && (t2 < t3))) return(1);
314     if( ((t2 < t3) && (t3 < t1))) return(1);
315     if( ((t3 < t1) && (t1 < t2))) return(1);
316 
317     return (0);
318 }
319 
320 void    conic_set_ends(Conic * conic, Vec2 p1, Vec2 p2, Vec2 pmid)
321 
322 /* set end points to be from p1 to p2 */
323 /* another point defining the sence of the curve */
324 {
325     double  t1, t2, tmid;
326 
327     if (conic == NULL || conic->type == DEGENERATE)
328         return;
329 
330     t1 = conic_param(conic, p1);
331     t2 = conic_param(conic, p2);
332 
333     if (conic->type == HYPERBOLA)
334     {
335         conic->branch = hyperbola_branch(conic, p1);
336         conic->t1 = t1;
337         conic->t2 = t2;
338         return;
339     }
340     tmid = conic_param(conic, pmid);
341 
342     if (ordered(t1, tmid, t2))
343     {
344         if (t1 > t2)
345             t2 += TWOPI;
346         conic->t1 = t1;
347         conic->t2 = t2;
348     } else
349     {
350         if (t2 >= t1)
351             t1 += TWOPI;
352         conic->t1 = t2;
353         conic->t2 = t1;
354     }
355 }
356 
357 double  conic_approx_length(Conic * conic, int n)
358 {
359     int     i;
360     Vec2    p = {Vec2_id};
361     Vec2    q = {Vec2_id};
362     double  dt = (conic->t2 - conic->t1) / n, len = 0.0;
363 
364     p = conic_point(conic, conic->t1);
365     for (i = 1; i <= n; i++)
366     {
367         q = conic_point(conic, conic->t1 + i * dt);
368         len += vec2_dist(p, q);
369         p = q;
370     }
371     return (len);
372 }
373 
374 double  conic_param_length(Conic * conic, Vec2 p1, Vec2 p2, Vec2 pmid)
375 {
376     double  t1, t2, tmid;
377 
378     if (conic == NULL || conic->type == DEGENERATE)
379         return (0.0);
380 
381     t1 = conic_param(conic, p1);
382     t2 = conic_param(conic, p2);
383 
384     if (conic->type == HYPERBOLA)
385         return (fabs(t2 - t1));
386 
387     tmid = conic_param(conic, pmid);
388     if (ordered(t1, tmid, t2))
389     {
390         if (t1 > t2)
391             t2 += TWOPI;
392     } else
393     {
394         if (t2 > t1)
395             t1 += TWOPI;
396     }
397     return (fabs(t2 - t1));
398 }
399 
400 
401 void    ellipse_format(Conic * conic)
402 {
403     Vec2    center = {Vec2_id};
404 
405     if (conic == NULL)
406         return;
407 
408     center = conic->center;
409     format("ellipse  : %15d\n", conic->label);
410     format("center   : %15.6f %15.6f\n", vec2_x(center), vec2_y(center));
411     format("theta    : %15.6f\n", conic->theta);
412     format("alpha    : %15.6f\n", conic->alpha);
413     format("beta     : %15.6f\n", conic->beta);
414 }
415 
416 void    hyperbola_format(Conic * conic)
417 {
418     Vec2    center = {Vec2_id};
419 
420     if (conic == NULL)
421         return;
422 
423     center = conic->center;
424     format("hyperbola: %15d\n", conic->label);
425     format("branch   : %3d\n", conic->branch);
426     format("center   : %15.6f %15.6f\n", vec2_x(center), vec2_y(center));
427     format("theta    : %15.6f\n", conic->theta);
428     format("alpha    : %15.6f\n", conic->alpha);
429     format("beta     : %15.6f\n", conic->beta);
430 }
431 
432 void    conic_format(Conic * conic)
433 {
434     Vec2    p = {Vec2_id};
435 
436     if (conic == NULL)
437         return;
438     switch (conic->type)
439     {
440     case ELLIPSE:
441         ellipse_format(conic);
442         break;
443     case HYPERBOLA:
444         hyperbola_format(conic);
445         break;
446     case DEGENERATE:
447     default:
448         return;
449     }
450     p = conic_point(conic, conic->t1);
451     format("p1       : %15.6f %15.6f\n", vec2_x(p), vec2_y(p));
452     p = conic_point(conic, conic->t2);
453     format("p2       : %15.6f %15.6f\n", vec2_x(p), vec2_y(p));
454     format("\n");
455 }
456 
457 /** ratio of minor to major axes (< 1.0) **/
458 double  conic_aratio(Conic * conic)
459 {
460     double  aratio;
461 
462     aratio = conic->alpha / conic->beta;
463     aratio = MIN(aratio, 1.0 / aratio);
464     return (aratio);
465 }
466 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.