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