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