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
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.