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

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

~ [ 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.