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

Linux Cross Reference
Tina6/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 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 

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