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

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

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