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

Linux Cross Reference
Tina5/tina-libs/tina/geometry/geomCurve_con_test.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 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 

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