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_test.c,v $
23 * Date : $Date: 2008/12/07 04:33:58 $
24 * Version : $Revision: 1.4 $
25 * CVS Id : $Id: geomCurve_con_test.c,v 1.4 2008/12/07 04:33:58 paul Exp $
26 *
27 * Author : Legacy TINA
28 *
29 * Notes :
30 *
31 *********
32 */
33
34 #include "geomCurve_con_test.h"
35
36 #if HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #include <math.h>
41 #include <string.h>
42 #include <tina/sys/sysDef.h>
43 #include <tina/sys/sysPro.h>
44 #include <tina/math/mathDef.h>
45 #include <tina/math/mathPro.h>
46 #include <tina/image/imgDef.h>
47 #include <tina/image/imgPro.h>
48 #include <tina/geometry/geomDef.h>
49 #include <tina/geometry/geom_GenPro.h>
50 #include <tina/geometry/geom_EdgeDef.h>
51 #include <tina/geometry/geom_EdgePro.h>
52 #include <tina/geometry/geom_CurveDef.h>
53 #include <tina/geometry/geom_IndxDef.h>
54 #include <tina/geometry/geom_IndxPro.h>
55 #include <tina/geometry/geomCurve_conic.h>
56 #include <tina/geometry/geomCurve_conicprox.h>
57 #include <tina/geometry/geomCurve_con_util.h>
58 #include <tina/geometry/geomCurve_con_fltr.h>
59 #include <tina/geometry/geomCurve_conic_5pt.h>
60 #ifdef _PCC
61 #include <memory.h>
62 #endif
63
64 double conic_chi2(Vec2 p, Conic * conic, Conic_stat * stats)
65 {
66 double px = vec2_x(p);
67 double py = vec2_y(p);
68 double z, h[5], uh[5];
69 double sum;
70 int i, j;
71
72 if (conic == NULL || stats == NULL)
73 return (0.0);
74
75 z = conic_eval(conic, p);
76 h[0] = px * px - py * py;
77 h[1] = 2.0 * px * py;
78 h[2] = 2.0 * px;
79 h[3] = 2.0 * py;
80 h[4] = 1.0;
81
82 for (i = 0; i < 5; i++)
83 {
84 sum = 0.0;
85 for (j = 0; j < i; j++)
86 sum += stats->u[j][i] * h[j];
87 sum += h[i];
88 uh[i] = sum;
89 }
90
91 sum = 0.0;
92 for (i = 0; i < 5; i++)
93 sum += uh[i] * stats->d[i] * uh[i];
94
95 return (chisq(z * z / sum, 1));
96 }
97
98 static double neglength(void *geom, int type)
99 {
100 switch (type)
101 {
102 case CONIC2:
103 return (-conic_approx_length((Conic *) geom, 2));
104 default:
105 return (0.0);
106 }
107 }
108
109 #if 0
110 static double closeness(void *geom, int type, Conic * conic1)
111 {
112 Vec2 p1 = {Vec2_id};
113 Vec2 p2 = {Vec2_id};
114 Vec2 q1 = {Vec2_id};
115 Vec2 q2 = {Vec2_id};
116 double d11, d12, d21, d22;
117
118 p1 = conic_point(conic1, conic1->t1);
119 p2 = conic_point(conic1, conic1->t2);
120
121 switch (type)
122 {
123 case CONIC2:
124 {
125 Conic *conic = (Conic *) geom;
126
127 q1 = conic_point(conic, conic->t1);
128 q2 = conic_point(conic, conic->t2);
129 break;
130 }
131 case LINE2:
132 {
133 Line2 *line = (Line2 *) geom;
134
135 q1 = line->p1;
136 q2 = line->p2;
137 break;
138 }
139 default:
140 return (1e10);
141 }
142
143 d11 = vec2_dist(p1, q1);
144 d12 = vec2_dist(p1, q2);
145 d21 = vec2_dist(p2, q1);
146 d22 = vec2_dist(p2, q2);
147 return (MIN4(d11, d12, d21, d22));
148 }
149
150 #endif
151
152 static Conic *conic_of_line(Line2 * line)
153 {
154 Tstring *str = (Tstring *) prop_get(line->props, STRING);
155 Vec2 p1 = {Vec2_id};
156 Vec2 p2 = {Vec2_id};
157 Vec2 p3 = {Vec2_id};
158
159 if (str->count < 3)
160 return (NULL);
161
162 DD_GET_POS2(str->start, p1);
163 DD_GET_POS2(str_mid_point(str), p2);
164 DD_GET_POS2(str->end, p3);
165
166 return (conic_circ_3pt(p1, p2, p3));
167 }
168
169 Bool conic_param_between(double t, double t1, double t2)
170 /* always in range 0 to 2PI */
171 /* t2 always > t1 and t1 is less than 2 PI */
172 {
173 if (t < t1)
174 t += TWOPI;
175
176 return ((t < t2) ? true : false);
177 }
178
179 #define MAXRATIO(a, b) MAX((a)/(b), (b)/(a))
180
181 /* check if a pair of conics are compatible for combination */
182 static Bool conics_compatible(Conic * conic1, Conic * conic2, double low_conf, double hi_conf)
183 {
184 Vec2 p1 = {Vec2_id};
185 Vec2 g1 = {Vec2_id};
186 Vec2 e1 = {Vec2_id};
187 Vec2 p2 = {Vec2_id};
188 Vec2 g2 = {Vec2_id};
189 Vec2 e2 = {Vec2_id};
190 Vec2 e11 = {Vec2_id};
191 Vec2 e12 = {Vec2_id};
192 Vec2 e21 = {Vec2_id};
193 Vec2 e22 = {Vec2_id};
194 Conic_stat *stats1;
195 Conic_stat *stats2;
196 double chi1, chi2;
197 double d11, d12, d21, d22;
198 double len1, len2;
199 double rad1, rad2;
200 double t, min_dist;
201
202 if (conic1 == NULL || conic2 == NULL || conic1->type != conic2->type)
203 {
204 /* this is a sign of corrupted data so watch out!
205 NAT 31/12/99 */
206 if (conic1->type != conic2->type)
207 error("incompatable structures in conics_compatible()",warning);
208 return (false);
209 }
210
211 /** conic mid points **/
212 p1 = conic_point(conic1, 0.5 * (conic1->t1 + conic1->t2));
213 p2 = conic_point(conic2, 0.5 * (conic2->t1 + conic2->t2));
214
215 stats1 = (Conic_stat *) prop_get(conic1->props, STATS);
216 stats2 = (Conic_stat *) prop_get(conic2->props, STATS);
217
218 chi1 = conic_chi2(p1, conic2, stats2);
219 chi2 = conic_chi2(p2, conic1, stats1);
220
221 /* test for evidence of continuation from chi test */
222 if (stats2 != NULL && stats1 != NULL)
223 {
224 if (chi1 < low_conf && chi2 < low_conf)
225 return (false);
226 } else if (stats2 != NULL)
227 {
228 if (chi1 < low_conf)
229 return (false);
230 } else if (stats1 != NULL)
231 {
232 if (chi2 < low_conf)
233 return (false);
234 }
235 /* get end points */
236 e11 = conic_point(conic1, conic1->t1);
237 e12 = conic_point(conic1, conic1->t2);
238 e21 = conic_point(conic2, conic2->t1);
239 e22 = conic_point(conic2, conic2->t2);
240
241 /* check the conics do not overlap */
242
243 t = conic_param(conic2, e11);
244 if (conic_param_between(t, conic2->t1, conic2->t2))
245 return (false);
246
247 t = conic_param(conic2, e12);
248 if (conic_param_between(t, conic2->t1, conic2->t2))
249 return (false);
250
251 t = conic_param(conic1, e21);
252 if (conic_param_between(t, conic1->t1, conic1->t2))
253 return (false);
254
255 t = conic_param(conic1, e21);
256 if (conic_param_between(t, conic1->t1, conic1->t2))
257 return (false);
258
259 /* distances between end points */
260 d11 = vec2_dist(e11, e21);
261 d12 = vec2_dist(e11, e22);
262 d21 = vec2_dist(e12, e21);
263 d22 = vec2_dist(e12, e22);
264 min_dist = MIN4(d11, d12, d21, d22);
265
266 /* if very close end points then TRUE */
267 if (min_dist < 10)
268 return (true);
269
270 if (min_dist == d11)
271 {
272 e1 = e11;
273 e2 = e21;
274 } else if (min_dist == d12)
275 {
276 e1 = e11;
277 e2 = e22;
278 } else if (min_dist == d21)
279 {
280 e1 = e12;
281 e2 = e21;
282 } else
283 {
284 e1 = e12;
285 e2 = e22;
286 }
287
288 len1 = conic_approx_length(conic1, 5);
289 len2 = conic_approx_length(conic2, 5);
290
291 /* if either of the lines is short then FALSE */
292 if (len1 < 20 || len2 < 20)
293 return (false);
294
295 /* if the gap is shorter than both of the lines then TRUE */
296 if (min_dist > len1 || min_dist > len2)
297 return (false);
298
299 g1 = vec2_unit(conic_grad(conic1, e1));
300 g2 = vec2_unit(conic_grad(conic2, e2));
301
302 /* if ends are co-tangental then TRUE */
303 if (fabs(vec2_dot(g1, g2)) > 0.9)
304 {
305 Vec2 t = {Vec2_id};
306
307 t = vec2_unit(vec2_diff(e2, e1));
308 if (fabs(vec2_dot(g1, t)) < 0.15 && fabs(vec2_dot(g2, t)) < 0.15)
309 return (true);
310 }
311 rad1 = MAX(conic1->alpha, conic1->beta);
312 rad2 = MAX(conic2->alpha, conic2->beta);
313
314 /** gross test for closeness **/
315 if (vec2_dist(p2, conic1->center) > 3.0 * rad1 &&
316 vec2_dist(p1, conic2->center) > 3.0 * rad2)
317 return (false);
318
319 if (chi1 > hi_conf || chi2 > hi_conf)
320 return (true);
321
322 /** gradients and concavity sign at test point match approximately **/
323 g1 = vec2_unit(conic_grad(conic1, p2));
324 g2 = vec2_unit(conic_grad(conic2, p2));
325 if (vec2_dot(g1, g2) > 0.95)
326 return (true);
327 g1 = vec2_unit(conic_grad(conic1, p1));
328 g2 = vec2_unit(conic_grad(conic2, p1));
329 if (vec2_dot(g1, g2) > 0.95)
330 return (true);
331
332 return (false);
333 }
334
335 static Bool conic_and_geom_compatible(Conic * conic1, void *geom, int type, double low_conf, double hi_conf)
336 {
337 switch (type)
338 {
339 case CONIC2:
340 return (conics_compatible(conic1, (Conic *) geom, low_conf, hi_conf));
341 case LINE2:
342 {
343 Conic *conic2 = conic_of_line((Line2 *) geom);
344 Bool result;
345
346 if (conic2 == NULL)
347 return (false);
348 result = conics_compatible(conic1, conic2, low_conf, hi_conf);
349 conic_free(conic2);
350 return (result);
351 }
352 default:
353 return (false);
354 }
355 }
356
357 static Conic *conic_join_strings(Tstring * str1, Tstring * str2, double lo_th, double hi_th, int max_div)
358 {
359 Tstring *str;
360 Conic *conic;
361 List *strings;
362
363 str = es_combine(str1, str2);
364 conic = (Conic *) conic_ellipse_fit(str->start, str->end, 0.2);
365 if (conic == NULL)
366 {
367 str_rm(str, (void (*) ()) NULL);
368 return (NULL);
369 }
370 conic->props = proplist_add(conic->props, (void *) str, STRING, str_rm_only_str);
371
372 strings = list_of(STRING, str1, STRING, str2, NULL);
373 if (!conic_check_list(conic, strings, lo_th, hi_th, max_div))
374 {
375 conic_free(conic);
376 list_rm_links(strings);
377 return (NULL);
378 }
379 list_rm_links(strings);
380 return (conic);
381 }
382
383 static Conic *conic_join_list(List * list, double lo_th, double hi_th, int max_div)
384 {
385 Tstring *str;
386 Conic *conic;
387
388 list = list_copy(list, (void *(*) ()) NULL, NULL);
389 list = es_list_order(list);
390 str = es_list_cat(list); /* total string */
391
392 conic = (Conic *) conic_ellipse_fit(str->start, str->end, 0.2);
393 if (conic == NULL)
394 {
395 str_rm(str, (void (*) ()) NULL);
396 list_rm_links(list);
397 return (NULL);
398 }
399 conic->props = proplist_add(conic->props, (void *) str, STRING, str_rm_only_str);
400
401 if (!conic_check_list(conic, list, lo_th, hi_th, max_div))
402 {
403 conic_free(conic);
404 list_rm_links(list);
405 return (NULL);
406 }
407 list_rm_links(list);
408 return (conic);
409 }
410
411 static Conic *conic_join_array(Tstring * str, int *t_array, Tstring ** s_array, int n, double lo_th, double hi_th, int max_div)
412 {
413 List *list = NULL;
414 Conic *conic;
415 int i;
416
417 for (i = 0; i < n; ++i)
418 if (t_array[i])
419 list = ref_addtostart(list, (void *) s_array[i], STRING);
420
421 if (list == NULL)
422 return (NULL);
423
424 list = ref_addtostart(list, (void *) str, STRING);
425 list = es_list_order(list);
426 str = es_list_cat(list); /* total string */
427
428 conic = (Conic *) conic_ellipse_fit(str->start, str->end, 0.2);
429 if (conic == NULL)
430 {
431 list_rm_links(list);
432 str_rm(str, (void (*) ()) NULL);
433 return (NULL);
434 }
435 conic->props = proplist_add(conic->props, (void *) str, STRING, str_rm_only_str);
436
437 if (!conic_check_list(conic, list, lo_th, hi_th, max_div))
438 {
439 conic_free(conic);
440 list_rm_links(list);
441 return (NULL);
442 }
443 list_rm_links(list);
444 return (conic);
445 }
446
447 static int rec_max_div; /* static data! */
448 static double rec_lo_thres, rec_hi_thres, max; /* static data! */
449 static Conic *best_conic=NULL; /* static data! */
450
451 static int *b_array=NULL; /* static data! */
452 static int *t_array=NULL; /* static data! */
453 static Tstring **s_array=NULL; /* static data! */
454 static Tstring *rec_str=NULL; /* static data! */
455 static int rec_n; /* static data! */
456
457 static void conic_rec_choose(int i, double val)
458 {
459 if (i == rec_n)
460 {
461 Conic *conic;
462
463 if (val < max)
464 return;
465 conic = conic_join_array(rec_str, t_array, s_array, rec_n,
466 rec_lo_thres, rec_hi_thres, rec_max_div);
467 if (conic == NULL)
468 return;
469 conic_free(best_conic);
470 best_conic = conic;
471 (void) memcpy((char *) b_array, (char *) t_array, rec_n * sizeof(int));
472 max = val;
473 return;
474 }
475 t_array[i] = 1;
476 conic_rec_choose(i + 1, val + s_array[i]->count);
477 t_array[i] = 0;
478 conic_rec_choose(i + 1, val);
479 }
480
481 static Conic *best_from_list(Conic * conic, List * list, List ** con_list, double lo_thres, double hi_thres, int max_div)
482 {
483 Conic *new_con;
484 List *ptr;
485 List *slist;
486 int i;
487
488 if (conic == NULL || list == NULL)
489 return (NULL);
490
491 slist = reclist_list_update(list, geom_prop_update_get, 0, (void *) STRING);
492 slist = ref_addtostart(slist, (void *) prop_get(conic->props, STRING), STRING);
493 new_con = conic_join_list(slist, lo_thres, hi_thres, max_div);
494 list_rm_links(slist);
495
496 if (new_con != NULL)
497 {
498 /* the whole list is consistent */
499 *con_list = list_copy(list, (void *(*) ()) NULL, NULL);
500 return (new_con);
501 }
502 best_conic = NULL;
503 rec_lo_thres = lo_thres;
504 rec_hi_thres = hi_thres;
505 rec_max_div = max_div;
506 rec_str = (Tstring *) prop_get(conic->props, STRING);
507 rec_n = list_length(list);
508 t_array = ivector_alloc(0, rec_n);
509 b_array = ivector_alloc(0, rec_n);
510 s_array = (Tstring **) pvector_alloc(0, rec_n);
511
512 for (i = 0, ptr = list; ptr != NULL; i++, ptr = ptr->next)
513 {
514 t_array[i] = b_array[i] = 0;
515 s_array[i] = (Tstring *) geom_prop_get(ptr->to, ptr->type, STRING);
516 }
517 #ifdef _ICC
518 {
519 double zero = 0.0;
520
521 conic_rec_choose(0, zero);
522 }
523 #else
524 conic_rec_choose(0, 0.0);
525 #endif
526
527 *con_list = NULL;
528 for (i = 0, ptr = list; ptr != NULL; i++, ptr = ptr->next)
529 if (b_array[i])
530 *con_list = ref_addtostart(*con_list, (void *) ptr->to, ptr->type);
531
532 ivector_free( t_array, 0);
533 ivector_free( b_array, 0);
534 pvector_free( s_array, 0);
535
536 return (best_conic);
537 }
538
539 static double str_tres; /* static data! */
540
541 void *conic_pos_thres(void *pos, int type, Conic * conic)
542 {
543 Vec2 p = {Vec2_id};
544 double t;
545
546 GET_POS2(pos, type, p);
547 t = conic_param(conic, p);
548 t = vec2_sqrdist(p, conic_point(conic, t));
549 if (t > str_tres)
550 return (NULL);
551 return (pos);
552 }
553
554 Tstring *conic_threshold_string(Conic * conic, Tstring * string, double thres)
555 {
556 str_tres = thres;
557 return (reclist_string_copy(string, conic_pos_thres, 0, (void *) conic));
558 }
559
560 Vec2 geom2_mid_point(void *geom, int type)
561 {
562 switch (type)
563 {
564 case POINT2:
565 return (((Point2 *) geom)->p);
566 case LINE2:
567 {
568 Vec2 p1 = {Vec2_id};
569 Vec2 p2 = {Vec2_id};
570
571 p1 = ((Line2 *) geom)->p1;
572 p2 = ((Line2 *) geom)->p2;
573
574 return (vec2_times(0.5, vec2_sum(p1, p2)));
575 }
576 case CONIC2:
577 {
578 double t1 = ((Conic *) geom)->t1;
579 double t2 = ((Conic *) geom)->t2;
580
581 return (conic_point((Conic *) geom, (t1 + t2) * 0.5));
582 }
583 default:
584 return (vec2_zero()); /* could be MAXFLOAT */
585 }
586 }
587
588 Vec2 geom2_p1(void *geom, int type)
589 {
590 switch (type)
591 {
592 case POINT2:
593 return (((Point2 *) geom)->p);
594 case LINE2:
595 return (((Line2 *) geom)->p1);
596 case CONIC2:
597 return (conic_point((Conic *) geom, ((Conic *) geom)->t1));
598 default:
599 return (vec2_zero()); /* could be MAXFLOAT */
600 }
601 }
602
603 Vec2 geom2_p2(void *geom, int type)
604 {
605 switch (type)
606 {
607 case POINT2:
608 return (((Point2 *) geom)->p);
609 case LINE2:
610 return (((Line2 *) geom)->p2);
611 case CONIC2:
612 return (conic_point((Conic *) geom, ((Conic *) geom)->t2));
613 default:
614 return (vec2_zero()); /* could be MAXFLOAT */
615 }
616 }
617
618 static void geom2_mid_add_to_index(void *geom, int type, Windex * index)
619 {
620 Ipos pos = {Ipos_id};
621
622 pos = wx_get_index(index, geom2_mid_point(geom, type));
623 wx_add_entry(index, geom, type, ipos_y(pos), ipos_x(pos));
624 }
625
626 static void geom2_mid_rm_from_index(void *geom, int type, Windex * index)
627 {
628 Ipos pos = {Ipos_id};
629
630 pos = wx_get_index(index, geom2_mid_point(geom, type));
631 wx_rm_entry(index, geom, ipos_y(pos), ipos_x(pos));
632 }
633
634 /* ARGSUSED quieten lint */
635 static void vec2_inc_region(Vec2 p, int type, Imregion * roi)
636 {
637 int x, y;
638
639 x = (int)floor(vec2_x(p));
640 y = (int)floor(vec2_y(p));
641
642 if (roi->lx == LARGEST_INT)
643 {
644 roi->lx = roi->ux = x;
645 roi->ly = roi->uy = y;
646 return;
647 }
648 if (x < roi->lx)
649 roi->lx = x;
650 if (y < roi->ly)
651 roi->ly = y;
652 if (x > roi->ux)
653 roi->ux = x;
654 if (y > roi->uy)
655 roi->uy = y;
656 }
657
658
659 static void geom_inc_region(void *geom, int type, Imregion * roi)
660 {
661 vec2_inc_region(geom2_mid_point(geom, type), VEC2, roi);
662 vec2_inc_region(geom2_p1(geom, type), VEC2, roi);
663 vec2_inc_region(geom2_p2(geom, type), VEC2, roi);
664 }
665
666 Windex *geom_mid_point_index_build(List * geom, Imregion * region)
667 {
668 Windex *index;
669 int m, n;
670
671 if (region == NULL)
672 {
673 region = roi_alloc(LARGEST_INT, LARGEST_INT, LARGEST_INT, LARGEST_INT);
674 reclist_list_apply(geom, geom_inc_region, 0, (void *) region);
675 region->lx -= 20;
676 region->ly -= 20;
677 region->ux += 20;
678 region->uy += 20;
679 }
680 m = (region->uy - region->ly) / 20 + 1;
681 n = (region->ux - region->lx) / 20 + 1;
682 index = wx_alloc(region, m, n, GEOMETRY_2);
683 reclist_list_apply(geom, geom2_mid_add_to_index, 0, (void *) index);
684 return (index);
685 }
686
687 List *geom_from_index(Conic * conic, Conic_stat * stats, Windex * index, char **mask, Ipos p, double conf)
688 {
689 int i, ii, j, jj;
690 Vec2 v = {Vec2_id};
691 Vec2 v1 = {Vec2_id};
692 Vec2 v2 = {Vec2_id};
693 List *list = NULL;
694 double t;
695
696 i = ipos_y(p);
697 j = ipos_x(p);
698
699 if (wx_in_index(index, i, j) == false || mask[i][j])
700 return (NULL);
701
702 mask[i][j] = 1;
703
704 v = wx_get_mid_pos2(index, p);
705 t = conic_param(conic, v);
706 if (t < conic->t1)
707 {
708 if (t + TWOPI < conic->t2)
709 return (NULL);
710 } else if (t < conic->t2)
711 return (NULL);
712
713 v1 = wx_get_pos2(index, p);
714 v2 = wx_get_pos2(index, ipos(j + 1, i + 1));
715
716 if (conic_chi2(v, conic, stats) < conf &&
717 conic_chi2(v1, conic, stats) < conf &&
718 conic_chi2(v2, conic, stats) < conf &&
719 conic_chi2(wx_get_pos2(index, ipos(j, i + 1)), conic, stats) < conf &&
720 conic_chi2(wx_get_pos2(index, ipos(j + 1, i)), conic, stats) < conf)
721 {
722 v = conic_point(conic, t); /* closest point on conic */
723
724 /* does the conic NOT pass through this grid point */
725 if (!BETWEEN(vec2_x(v), vec2_x(v1), vec2_x(v2)) ||
726 !BETWEEN(vec2_y(v), vec2_y(v1), vec2_y(v2)))
727 return (NULL);
728 }
729 list = list_copy((List *) wx_get(index, i, j), (void *(*) ()) NULL, NULL);
730
731 for (ii = -1; ii <= 1; ++ii)
732 for (jj = -1; jj <= 1; ++jj)
733 if (ii != 0 || jj != 0)
734 {
735 List *sub;
736
737 p = ipos(j + jj, i + ii);
738 sub = geom_from_index(conic, stats, index, mask, p, conf);
739 list = list_append(sub, list);
740 }
741 return (list);
742 }
743
744 static Bool conic_below_thres(Conic * c, double ang_th, double len_th)
745
746 /* angle threshold in radians */
747 /* angle threshold in radians */
748 {
749 if (c == NULL)
750 return (true);
751
752 if (c->t2 - c->t1 > ang_th)
753 return (false);
754
755 if (conic_approx_length(c, 10) > len_th)
756 return (false);
757
758 return (true);
759 }
760
761 List *conic_join(List * geom, Imregion * roi, double conf, double lo_thres, double hi_thres, int max_div)
762 {
763 List *newgeom;
764 List *ptr1;
765 Windex *index;
766 char **mask;
767 int m, n;
768 int i, j;
769
770 if (geom == NULL)
771 return (NULL);
772
773 /** flatten input list, then sort by approx length **/
774 newgeom = reclist_list_flat(geom, (void *(*) ()) NULL, 0, NULL);
775 newgeom = sort_list(newgeom, neglength, NULL);
776
777 index = geom_mid_point_index_build(newgeom, roi);
778 m = index->m;
779 n = index->n;
780 mask = carray_alloc(0, 0, m, n);
781
782 ptr1 = newgeom;
783 while (ptr1 != NULL)
784 {
785 Conic *conic;
786 Conic *new_conic;
787 Conic_stat *stats;
788 List *ptr2;
789 List *test_geom;
790 Tstring *str1, *str2;
791 Bool changed = false;
792 Vec2 v = {Vec2_id};
793 Ipos p = {Ipos_id};
794
795 if (ptr1->type != CONIC2 || conic_below_thres((Conic *) ptr1->to, 0.5, 20.0))
796 {
797 ptr1 = ptr1->next;
798 continue;
799 }
800 conic = (Conic *) ptr1->to;
801
802 stats = (Conic_stat *) prop_get(conic->props, STATS);
803 if (stats == NULL)
804 {
805 ptr1 = ptr1->next;
806 continue;
807 }
808 for (i = 0; i < m; ++i)
809 for (j = 0; j < n; ++j)
810 mask[i][j] = 0;
811
812 v = conic_point(conic, PI + (conic->t1 + conic->t2) * 0.5);
813 p = wx_get_index(index, v);
814 test_geom = geom_from_index(conic, stats, index, mask, p, conf);
815 v = conic_point(conic, conic->t1 - 0.1);
816 p = wx_get_index(index, v);
817 ptr2 = geom_from_index(conic, stats, index, mask, p, conf);
818 test_geom = list_append(ptr2, test_geom);
819 v = conic_point(conic, conic->t2 + 0.1);
820 p = wx_get_index(index, v);
821 ptr2 = geom_from_index(conic, stats, index, mask, p, conf);
822 test_geom = list_append(ptr2, test_geom);
823
824 do
825 {
826 List *compat = NULL;
827 List *chosen = NULL;
828
829 str1 = (Tstring *) prop_get(conic->props, STRING);
830
831 /** find candidates for joining **/
832 for (ptr2 = test_geom; ptr2 != NULL; ptr2 = ptr2->next)
833 {
834 if (ptr1->to == ptr2->to) /* the same conic */
835 continue;
836
837 if (conic_and_geom_compatible(conic, ptr2->to, ptr2->type, conf, 0.1))
838 {
839 str2 = (Tstring *) geom_prop_get(ptr2->to, ptr2->type, STRING);
840 new_conic = conic_join_strings(str1, str2, lo_thres, hi_thres, max_div);
841 if (new_conic == NULL)
842 continue;
843
844 conic_free(new_conic);
845 compat = ref_addtostart((List *) compat, (void *) ptr2->to, ptr2->type);
846 }
847 }
848
849 new_conic = best_from_list(conic, compat, &chosen, lo_thres, hi_thres, max_div);
850 list_rm_links(compat);
851
852 if (new_conic != NULL) /** new fit suceeded **/
853 {
854 changed = true;
855 for (ptr2 = chosen; ptr2 != NULL; ptr2 = ptr2->next)
856 {
857 test_geom = list_rm_ref(test_geom, ptr2->to, (void (*) ()) NULL);
858 geom2_mid_rm_from_index(ptr2->to, ptr2->type, index);
859 newgeom = list_rm_ref(newgeom, ptr2->to, (void (*) ()) geom_free);
860 }
861 list_rm_links(chosen);
862
863 str2 = (Tstring *) prop_get(new_conic->props, STRING);
864 str2 = conic_threshold_string(new_conic, str2, lo_thres);
865 (void) prop_set(new_conic->props, (void *) str2, STRING, true);
866
867 conic_free(conic);
868
869 geom2_mid_rm_from_index(ptr1->to, ptr1->type, index);
870 ptr1->to = conic = new_conic;
871 geom2_mid_add_to_index(ptr1->to, ptr1->type, index);
872 }
873 } while (new_conic != NULL);
874 list_rm_links(test_geom);
875 if (changed == false)
876 ptr1 = ptr1->next;
877 }
878 (void) reclist_list_free(geom, (void (*) ()) NULL, 0, NULL);
879 wx_free(index, list_rm_links);
880
881 return (newgeom);
882 }
883
884 List *conic_compatible(Conic * conic, List * geom, double lo_thres, double hi_thres, int max_div)
885 {
886 Conic *new_conic;
887 List *ptr;
888 List *candidates = NULL;
889 Tstring *str1, *str2;
890
891 geom = reclist_list_flat(geom, (void *(*) ()) NULL, 0, NULL);
892 str1 = (Tstring *) prop_get(conic->props, STRING);
893
894 /** find candidates for joining **/
895 for (ptr = geom; ptr != NULL; ptr = ptr->next)
896 {
897 if (conic == ptr->to) /* the same conic */
898 continue;
899
900 if (conic_and_geom_compatible(conic, ptr->to, ptr->type, 0.01, 0.1))
901 {
902 str2 = (Tstring *) geom_prop_get(ptr->to, ptr->type, STRING);
903 new_conic = conic_join_strings(str1, str2, lo_thres, hi_thres, max_div);
904 if (new_conic == NULL)
905 continue;
906 conic_free(new_conic);
907 candidates = ref_addtostart((List *) candidates, (void *) ptr->to, ptr->type);
908 }
909 }
910 list_rm_links(geom);
911 return (candidates);
912 }
913
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.