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

Linux Cross Reference
Tina5/tina-libs/tina/geometry/geomCurve_conic_5pt.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_conic_5pt.c,v $
 23  * Date    :  $Date: 2008/12/07 04:33:58 $
 24  * Version :  $Revision: 1.2 $
 25  * CVS Id  :  $Id: geomCurve_conic_5pt.c,v 1.2 2008/12/07 04:33:58 paul Exp $
 26  *
 27  * Author  : Legacy TINA
 28  *
 29  * Notes :
 30  *
 31  * recover coeffs of conic through 5pts by
 32  * four line method
 33  *
 34  *********
 35 */
 36 
 37 #include "geomCurve_conic_5pt.h"
 38 
 39 #if HAVE_CONFIG_H
 40   #include <config.h>
 41 #endif
 42 
 43 #include <math.h>
 44 #include <tina/sys/sysDef.h>
 45 #include <tina/sys/sysPro.h>
 46 #include <tina/math/mathDef.h>
 47 #include <tina/math/mathPro.h>
 48 #include <tina/geometry/geomDef.h>
 49 #include <tina/geometry/geom_CurveDef.h>
 50 #include <tina/geometry/geomCurve_conic.h>
 51 #include <tina/geometry/geomCurve_con_util.h>
 52 
 53 
 54 /**
 55 coefficients of line  L = l0*x+l0*y+l2 = 0
 56 
 57 through  (x0 y0)  and  (x1 y1)
 58 **/
 59 static void line_coeffs(double *line, Vec2 p1, Vec2 p2)
 60 {
 61     line[0] = vec2_y(p2) - vec2_y(p1);
 62     line[1] = vec2_x(p1) - vec2_x(p2);
 63     line[2] = vec2_x(p2) * vec2_y(p1) - vec2_x(p1) * vec2_y(p2);
 64 }
 65 
 66 /**
 67 
 68 coefficients of degenerate conic (pair of lines)
 69 C = L0*L1
 70 **/
 71 static Conic *conic_from_2_lines(double *l0, double *l1)
 72 {
 73     Conic  *conic = conic_alloc(0);
 74 
 75     conic->a = l0[0] * l1[0];
 76     conic->b = 0.5 * (l0[0] * l1[1] + l0[1] * l1[0]);
 77     conic->c = l0[1] * l1[1];
 78     conic->d = 0.5 * (l0[0] * l1[2] + l0[2] * l1[0]);
 79     conic->e = 0.5 * (l0[1] * l1[2] + l0[2] * l1[1]);
 80     conic->f = l0[2] * l1[2];
 81     return (conic);
 82 }
 83 
 84 /**
 85 linear combination  C = t*C0+(1-t)*C1  of two conics
 86 **/
 87 static Conic *conic_combine(Conic * conic0, Conic * conic1, double t)
 88 {
 89     Conic  *conic = conic_alloc(0);
 90     double  t1 = 1.0 - t;
 91 
 92     conic->a = t * conic0->a + t1 * conic1->a;
 93     conic->b = t * conic0->b + t1 * conic1->b;
 94     conic->c = t * conic0->c + t1 * conic1->c;
 95     conic->d = t * conic0->d + t1 * conic1->d;
 96     conic->e = t * conic0->e + t1 * conic1->e;
 97     conic->f = t * conic0->f + t1 * conic1->f;
 98     return (conic);
 99 }
100 
101 /**
102 returns circular conic through three points
103 **/
104 Conic  *conic_circ_3pt(Vec2 p1, Vec2 p2, Vec2 p3)
105 {
106     double  r;
107     Vec2    c = {Vec2_id};
108     Conic  *conic;
109 
110     vec2_circ_3_points(p1, p2, p3, &c, &r);
111     conic = conic_make(ELLIPSE, c, 0.0, r, r, 0.0, TWOPI, 0);
112     conic_set_ends(conic, p1, p3, p2);
113     return (conic);
114 }
115 
116 /**
117 conic passing through  5  given points
118 calculated by the 4-line method
119 **/
120 Conic  *conic_5pt(Vec2 p1, Vec2 p2, Vec2 p3, Vec2 p4, Vec2 p5)
121 {
122     Conic  *conic;
123     Conic  *conic0;
124     Conic  *conic1;
125     double  l0[3], l1[3], l2[3], l3[3];
126 
127     line_coeffs(l0, p1, p2);
128     line_coeffs(l1, p2, p3);
129     line_coeffs(l2, p3, p4);
130     line_coeffs(l3, p4, p1);
131     conic0 = conic_from_2_lines(l0, l2);
132     conic1 = conic_from_2_lines(l1, l3);
133     conic = conic_combine(conic0, conic1, conic_lambda(conic0, conic1, p5));
134     conic_free(conic0);
135     conic_free(conic1);
136     conic_setup(conic);
137     return (conic);
138 }
139 
140 /* Ellipse passing through  5  given points if it exists aspect ratio
141  * is too small or 5pt fit is ellipse returns 3pt circle */
142 Conic  *conic_ellipse_5pt(Vec2 p1, Vec2 p2, Vec2 p3, Vec2 p4, Vec2 p5, double min_aratio)
143 {
144     Conic  *conic = conic_5pt(p1, p2, p3, p4, p5);
145 
146     if (conic->type == HYPERBOLA || conic_aratio(conic) < min_aratio)
147     {
148         conic_free(conic);
149         return (conic_circ_3pt(p1, p3, p5));
150     } else
151         return (conic);
152 }
153 
154 /**
155 conic given 2  tangencies and one other point
156 **/
157 Conic  *conic_3pt(Vec2 p1, Vec2 v1, Vec2 p2, Vec2 v2, Vec2 p3)
158 {
159     Conic  *conic;
160     Conic  *conic0;
161     Conic  *conic1;
162     double  l0[3], l1[3], l2[3];
163 
164     line_coeffs(l0, p1, vec2_sum(p1, v1));
165     line_coeffs(l1, p2, vec2_sum(p2, v2));
166     line_coeffs(l2, p1, p2);
167     conic0 = conic_from_2_lines(l0, l1);
168     conic1 = conic_from_2_lines(l2, l2);
169     conic = conic_combine(conic0, conic1, conic_lambda(conic0, conic1, p3));
170     conic_free(conic0);
171     conic_free(conic1);
172     conic_setup(conic);
173     return (conic);
174 }
175 
176 /**
177 ellipse given 2  tangencies and one other point if it exists
178 aspect ratio is too small or 5pt fit is ellipse returns
179 3pt circle
180 **/
181 Conic  *conic_ellipse_3pt(Vec2 p1, Vec2 v1, Vec2 p2, Vec2 v2, Vec2 p3, double min_aratio)
182 {
183     Conic  *conic = conic_3pt(p1, v1, p2, v2, p3);
184 
185     /** temporary change **/
186     if (conic->type == HYPERBOLA || conic_aratio(conic) < min_aratio)
187     {
188         conic_free(conic);
189         return (conic_circ_3pt(p1, p2, p3));
190     } else
191         return (conic);
192 }
193 
194 /**
195 5 point  conic fit to dd string -
196 looks for well separated points
197 **/
198 Conic  *ddstr_conic_5pt(List * start, List * end)
199 {
200     List *m1;
201     List *m2;
202     List *m3;
203     Vec2    p0 = {Vec2_id};
204     Vec2    p1 = {Vec2_id};
205     Vec2    p2 = {Vec2_id};
206     Vec2    p3 = {Vec2_id};
207     Vec2    p4 = {Vec2_id};
208     Conic  *conic;
209 
210     if (start == end)
211         return (NULL);
212 
213     m2 = ddstr_mid_point(start, end);
214     m1 = ddstr_mid_point(start, m2);
215     m3 = ddstr_mid_point(m2, end);
216 
217     if (start == m1 || m1 == m2 || m2 == m3 || m3 == end)       /* not distinct points */
218         return (NULL);
219 
220     DD_GET_POS2(start, p0);
221     DD_GET_POS2(m1, p1);
222     DD_GET_POS2(m2, p2);
223     DD_GET_POS2(m3, p3);
224     DD_GET_POS2(end, p4);
225 
226     if (vec2_mod(vec2_diff(p0, p4)) < vec2_mod(vec2_diff(p0, p1)))
227         conic = ddstr_conic_5pt(m1, end);       /* end points too close */
228     else
229         conic = conic_5pt(p0, p1, p2, p3, p4);
230     conic_set_ends(conic, p0, p4, p2);
231     return (conic);
232 }
233 
234 /* 5 point ellipse fit to point string - looks for well separated
235  * points */
236 Conic  *ddstr_conic_ellipse_5pt(List * start, List * end, double min_aratio)
237 {
238     List *m1;
239     List *m2;
240     List *m3;
241     Vec2    p0 = {Vec2_id};
242     Vec2    p1 = {Vec2_id};
243     Vec2    p2 = {Vec2_id};
244     Vec2    p3 = {Vec2_id};
245     Vec2    p4 = {Vec2_id};
246     Conic  *conic;
247 
248     if (start == end)
249         return (NULL);
250 
251     m2 = ddstr_mid_point(start, end);
252     m1 = ddstr_mid_point(start, m2);
253     m3 = ddstr_mid_point(m2, end);
254 
255     if (start == m1 || m1 == m2 || m2 == m3 || m3 == end)       /* not distinct points */
256         return (NULL);
257 
258     DD_GET_POS2(start, p0);
259     DD_GET_POS2(m1, p1);
260     DD_GET_POS2(m2, p2);
261     DD_GET_POS2(m3, p3);
262     DD_GET_POS2(end, p4);
263 
264     if (vec2_mod(vec2_diff(p0, p4)) < vec2_mod(vec2_diff(p0, p1)))
265         conic = ddstr_conic_ellipse_5pt(m1, end, min_aratio);
266     else
267         conic = conic_ellipse_5pt(p0, p1, p2, p3, p4, min_aratio);
268     conic_set_ends(conic, p0, p4, p2);
269     return (conic);
270 }
271 
272 /**
273 3 point circle fit to point string -
274 looks for well separated points
275 **/
276 Conic  *ddstr_conic_circ_3pt(List * start, List * end)
277 {
278     List *mid;
279     Vec2    p0 = {Vec2_id};
280     Vec2    p1 = {Vec2_id};
281     Vec2    p2 = {Vec2_id};
282     Conic  *conic;
283 
284     if (start == end)
285         return (NULL);
286 
287     mid = ddstr_mid_point(start, end);
288 
289     if (start == mid || mid == end)     /* not distinct points */
290         return (NULL);
291 
292     DD_GET_POS2(start, p0);
293     DD_GET_POS2(mid, p1);
294     DD_GET_POS2(end, p2);
295 
296     if (vec2_mod(vec2_diff(p0, p2)) < vec2_mod(vec2_diff(p0, p1)) / 2)
297         conic = ddstr_conic_circ_3pt(start, ddstr_mid_point(mid, end));
298     else
299         conic = conic_circ_3pt(p0, p1, p2);
300     conic_set_ends(conic, p0, p2, p1);
301     return (conic);
302 }
303 
304 /** chooses same points on string as conic fitting **/
305 
306 Bool    ddstr_5pts_choose(List * start, List * end, Vec2 * p)
307 {
308     List *m1;
309     List *m2;
310     List *m3;
311 
312     if (start == end)
313         return (false);
314 
315     m2 = ddstr_mid_point(start, end);
316     m1 = ddstr_mid_point(start, m2);
317     m3 = ddstr_mid_point(m2, end);
318 
319     if (start == m1 || m1 == m2 || m2 == m3 || m3 == end)
320         /* not distinct points */
321         return (false);
322 
323     DD_GET_POS2(start, p[0]);
324     DD_GET_POS2(m1, p[1]);
325     DD_GET_POS2(m2, p[2]);
326     DD_GET_POS2(m3, p[3]);
327     DD_GET_POS2(end, p[4]);
328 
329     if (vec2_mod(vec2_diff(p[0], p[4])) < vec2_mod(vec2_diff(p[0], p[1])))
330         /* end points too close */
331         return (ddstr_5pts_choose(m1, end, p));
332     return (true);
333 }
334 
335 double  line_error(Line2 * line, Vec2 p)
336 {
337     return (vec2_cross(line->v, vec2_diff(p, line->p)));
338 }
339 
340 void    line_errors_check(Line2 * line, double *sum1, double *sum2)
341 {
342     Tstring *str;
343     List *ptr;
344 
345     if (line == NULL)
346         return;
347     str = prop_get(line->props, STRING);
348     if (str == NULL)
349         return;
350 
351     for (ptr = str->start;; ptr = ptr->next)
352     {
353         double  e;
354         Vec2    p = {Vec2_id};
355 
356         DD_GET_POS2(ptr, p);
357         e = line_error(line, p);
358         *sum1 += e;
359         *sum2 += e * e;
360         if (ptr == str->end)
361             break;
362     }
363     *sum1 /= str->count;
364     *sum2 /= str->count;
365 }
366 

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